Re: [julia-users] Re: How to built an array with some keys of a dictionary ?

[Dan]

a = Vector{String}



is also a good option.

On Monday, November 7, 2016 at 5:27:49 PM UTC+2, Steven G. Johnson wrote:
>
>
>
> On Monday, November 7, 2016 at 9:45:37 AM UTC-5, Mauro wrote:
>>
>> On Mon, 2016-11-07 at 15:27, Fred  wrote: 
>> > Hi Steven, 
>> > 
>> > 
>> > I also tried a = Array{String}() unfortunately it produces errors as 
>> well. 
>>
>> Array{String}(0) 
>>
>
> Whoops, sorry about the typo. 
>

Re: [julia-users] Re: How to built an array with some keys of a dictionary ?

[Dan]

I meant,

a = Vector{String}()

Dan

On Tuesday, November 8, 2016 at 9:39:53 AM UTC+2, Dan wrote:
>
> a = Vector{String}
>
>
>
> is also a good option.
>
> On Monday, November 7, 2016 at 5:27:49 PM UTC+2, Steven G. Johnson wrote:
>>
>>
>>
>> On Monday, November 7, 2016 at 9:45:37 AM UTC-5, Mauro wrote:
>>>
>>> On Mon, 2016-11-07 at 15:27, Fred  wrote: 
>>> > Hi Steven, 
>>> > 
>>> > 
>>> > I also tried a = Array{String}() unfortunately it produces errors as 
>>> well. 
>>>
>>> Array{String}(0) 
>>>
>>
>> Whoops, sorry about the typo. 
>>
>

[julia-users] Re: jump ahead in Mersenne Twister

[Dan]

Hi Joaquim,

The *jumppoly* parameter to randjump is a string which depends on the 
specific parameters of the MersenneTwister. By default Julia uses the 19937 
parameters in SFMT. To calculate these strings a utility was (or is still) 
supplied with SFMT called calc-jump. Its homepage is 
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/JUMP/ . Compiling it 
is possible (I've managed to do so now), and then the command:
./calc-jump 1 poly.19937.txt
outputs the jumppoly used by Julia to setup multiple "independent" random 
generators (10^20 steps apart in the generator).
The smallest jump, by 1 step has the simple jumppoly equal to the string 
"2".
Finally, note that the Julia random generator can get two Float64s from 
each step of the Twister, and there is also a cache of about 382 Float64s, 
so jumping the Twister will not influence the generated random unless the 
cache is cleared and getting to the right position might require a *2 
somewhere.
To setup a fixed *jumppoly* parameter of your choice, you should be able to 
use calc-jump. To do it dynamically in Julia would require ccall-ing the 
SFMT routine in SFMT Jump which is not too hard as well.

G'day, Dan Getz
 
On Friday, October 21, 2016 at 3:43:05 AM UTC+3, Joaquim Masset Lacombe 
Dias Garcia wrote:
>
> Hi,
>
> Have anybody implemented an arbitrary jump ahead for Mersenne Twister in 
> julia?
>
> I saw that the randjump documentation that goes as following:
>
> randjump(*r::MersenneTwister*, *jumps*[, *jumppoly*]) → 
> Vector{MersenneTwister}
>
> Create an array of the size jumps of initialized MersenneTwister RNG 
> objects where the first RNG object given as a parameter and following 
> MersenneTwister RNGs in the array initialized such that a state of the 
> RNG object in the array would be moved forward (without generating numbers) 
> from a previous RNG object array element on a particular number of steps 
> encoded by the jump polynomial jumppoly
>
> Given that, my question would reduce to: how do I encode an arbitrary 
> number of step (N) in a *jumppoly* ?
>
> thanks!
>

Re: [julia-users] Please help me to convert a BitArray to a Integer array in Julia

[Dan]

If saving characters is a thing, then:
julia> a = rand(Bool,3,2)
3×2 Array{Bool,2}:
 false  false
  true  false
 false   true


julia> 1a
3×2 Array{Int64,2}:
 0  0
 1  0
 0  1



On Monday, October 17, 2016 at 2:08:44 PM UTC+3, Scott Jones wrote:
>
> Tim, do you know if there is any difference in performance between the two 
> methods?
>
> Note, Sujoy, the first method that Tim showed is only available on v0.5 
> and later (some of the nice stuff added to entice people to move off of 
> v0.4.x to v0.5.x ;-) )
>
> On Monday, October 17, 2016 at 5:48:20 AM UTC-4, Tim Holy wrote:
>>
>> julia> a = bitrand(3,5)
>> 3×5 BitArray{2}:
>>   true  false  false  true   true 
>>
>  false   true   true  true  false
>>   true   true   true  true   true
>>
>> julia> Int.(a)
>> 3×5 Array{Int64,2}:
>>  1  0  0  1  1
>>  0  1  1  1  0
>>  1  1  1  1  1
>>
>> julia> convert(Array{Int}, a)
>> 3×5 Array{Int64,2}:
>>  1  0  0  1  1
>>  0  1  1  1  0
>>  1  1  1  1  1
>>
>>
>> On Mon, Oct 17, 2016 at 2:55 AM, Sujoy Datta  wrote:
>>
>>>
>>>
>>> I am a new user of Julia. Please help me to convert a nxm BitArray to an 
>>> nxm IntegerArray.
>>> What I want is to print 1 for 'true' and 0 for 'false'.
>>> Thank you in advance.
>>>
>>
>>

[julia-users] Re: julz: a command line utility for creating scalable Julia applications

[Dan Segal]

This is what @dpsanders brought up in issue #22 
. 

Summarizing my response there (repeated below), I wanted to get a prototype 
working to gauge the likelihood of adoption.



"
it probably could be, the truth is:
+ I wanted to bang it out in a weekend and get a working prototype
+ I found a good tutorial for building a python cli 

"


On Thursday, October 6, 2016 at 1:49:37 PM UTC-4, datn...@gmail.com wrote:
>
> Looks great. Out of curiosity, why didn't you want to write it in Julia?

[julia-users] Re: julz: a command line utility for creating scalable Julia applications

[Dan Segal]

Sorry this wasn't clear!

   - Julz  is a framework for building 
   Julia projects. 
   - It is *not* designed for making webpages, I was just trying to 
   highlight the analogy  
   
The reason for the weird structure of the project is that: 

   - this is a CLI tool written in Python
   - that creates Julia projects (an example project being the dummy folder 
   on the github page )

Finally,

   - Gemfile and Gemfile.lock were lifted from Ruby because I think they 
   provide a extensive package solution 
   - (I should probably just change them to REQUIRE and REQUIRE.lock though)

[julia-users] Re: julz: a command line utility for creating scalable Julia applications

[Dan Segal]

One more thing I wanted to say is, 

   - I tried to do with julz what ember-cli  did 
   for ember.js 
   - I really think ember-cli shaped their community into what it is today

--
 

> *Code*
>
>- github repo: https://github.com/djsegal/julz
>- (see readme.rst as well as included dummy Julia project)
>
> *Summary*
>
> Julz is a command line utility for creating ambitious Julia applications. 
> Its goal is to remove many of the petty decisions involved in starting a 
> project and establish a standard that can be adopted by the entire Julia 
> ecosystem. In this way it channels Ruby on Rails’ mantra of “convention 
> over configuration”: tell people where files should go, but allow them to 
> tweak it if they so desire.
>
> *Description*
>
> The problem Julz attempts to resolve is architecting a Julia project that 
> is resilient to both the frequent on-boarding/turnover of engineers as well 
> as the mass accumulation of files. In this way, a person reading a codebase 
> for the first time does not get bogged down by the sheer number of files 
> that have accumulated over time inside a shallow depth src directory.
>
> The way to alleviate this problem is to treat the src folder and the test 
> folder as boxes that hold smaller boxes of like content. In a Julia 
> context, this means that the src and test folder each have corresponding 
> subdirectories for: methods, types, functions, macros, etc. Each 
> subdirectory then has an associated generator command that can be used to 
> make files of that format. For example, generating a Rational type with two 
> integer fields from the command line would use boiler-plate templates to 
> compile:
>
>- a working rational.jl file in the ‘src/types’ directory
>- an associated rational_test.jl in the ‘test/types’ directory
>
> This may not seem like a big idea at first, but getting every project to 
> adopt this level-deeper folder structure is the cornerstone to a successful 
> framework. With it, packages can be created that expand on the native 
> collection of data structures, i.e. to produce singletons, compilers, etc. 
> Here, the addition of generators just make everything less user-error prone 
> and promotes adequate test coverage of code.
>
>

[julia-users] Re: julz: a command line utility for creating scalable Julia applications

[Dan Segal]

   
   - One final note is I tried to do with julz what ember-cli 
    did for Ember.js 
   - I really think the addition of ember-cli into the ecosystem changed 
   the community in a way they never thought possible

[julia-users] Re: julz: a command line utility for creating scalable Julia applications

[Dan Segal]

   
   - One final note is I tried to do with julz what ember-cli 
   <https://ember-cli.com/> did for Ember.js <http://emberjs.com/>
   - I really think the addition of ember-cli into the ecosystem changed 
   the community in a way they never thought possible
   

On Thursday, October 6, 2016 at 10:12:21 AM UTC-4, Dan Segal wrote:
>
> *Code*
>
>- github repo: https://github.com/djsegal/julz
>- (see readme.rst as well as included dummy Julia project)
>
> *Summary*
>
> Julz is a command line utility for creating ambitious Julia applications. 
> Its goal is to remove many of the petty decisions involved in starting a 
> project and establish a standard that can be adopted by the entire Julia 
> ecosystem. In this way it channels Ruby on Rails’ mantra of “convention 
> over configuration”: tell people where files should go, but allow them to 
> tweak it if they so desire.
>
> *Description*
>
> The problem Julz attempts to resolve is architecting a Julia project that 
> is resilient to both the frequent on-boarding/turnover of engineers as well 
> as the mass accumulation of files. In this way, a person reading a codebase 
> for the first time does not get bogged down by the sheer number of files 
> that have accumulated over time inside a shallow depth src directory.
>
> The way to alleviate this problem is to treat the src folder and the test 
> folder as boxes that hold smaller boxes of like content. In a Julia 
> context, this means that the src and test folder each have corresponding 
> subdirectories for: methods, types, functions, macros, etc. Each 
> subdirectory then has an associated generator command that can be used to 
> make files of that format. For example, generating a Rational type with two 
> integer fields from the command line would use boiler-plate templates to 
> compile:
>
>- a working rational.jl file in the ‘src/types’ directory
>- an associated rational_test.jl in the ‘test/types’ directory
>
> This may not seem like a big idea at first, but getting every project to 
> adopt this level-deeper folder structure is the cornerstone to a successful 
> framework. With it, packages can be created that expand on the native 
> collection of data structures, i.e. to produce singletons, compilers, etc. 
> Here, the addition of generators just make everything less user-error prone 
> and promotes adequate test coverage of code.
>
>

[julia-users] julz: a command line utility for creating scalable Julia applications

[Dan Segal]

*Code*

   - github repo: https://github.com/djsegal/julz
   - (see readme.rst as well as included dummy Julia project)

*Summary*

Julz is a command line utility for creating ambitious Julia applications. 
Its goal is to remove many of the petty decisions involved in starting a 
project and establish a standard that can be adopted by the entire Julia 
ecosystem. In this way it channels Ruby on Rails’ mantra of “convention 
over configuration”: tell people where files should go, but allow them to 
tweak it if they so desire.

*Description*

The problem Julz attempts to resolve is architecting a Julia project that 
is resilient to both the frequent on-boarding/turnover of engineers as well 
as the mass accumulation of files. In this way, a person reading a codebase 
for the first time does not get bogged down by the sheer number of files 
that have accumulated over time inside a shallow depth src directory.

The way to alleviate this problem is to treat the src folder and the test 
folder as boxes that hold smaller boxes of like content. In a Julia 
context, this means that the src and test folder each have corresponding 
subdirectories for: methods, types, functions, macros, etc. Each 
subdirectory then has an associated generator command that can be used to 
make files of that format. For example, generating a Rational type with two 
integer fields from the command line would use boiler-plate templates to 
compile:

   - a working rational.jl file in the ‘src/types’ directory
   - an associated rational_test.jl in the ‘test/types’ directory

This may not seem like a big idea at first, but getting every project to 
adopt this level-deeper folder structure is the cornerstone to a successful 
framework. With it, packages can be created that expand on the native 
collection of data structures, i.e. to produce singletons, compilers, etc. 
Here, the addition of generators just make everything less user-error prone 
and promotes adequate test coverage of code.

[julia-users] julz: a command line utility for creating scalable Julia applications

[Dan Segal]

Summary

Julz is a command line utility for creating ambitious Julia applications. 
Its goal is to remove many of the petty decisions involved in starting a 
project and establish a standard that can be adopted by the entire Julia 
ecosystem. In this way it channels Ruby on Rails’ mantra of “convention 
over configuration”: tell people where files should go, but allow them to 
tweak it if they so desire.


   - 
   
   Code available at: https://github.com/djsegal/julz 
   - 
   
   (see readme.rst as well as included dummy Julia project)
   

Description

>> julz new   // start new julz project

The problem Julz attempts to resolve is architecting a Julia project that 
is resilient to both the frequent on-boarding/turnover of engineers as well 
as the mass accumulation of files. In this way, a person reading a codebase 
for the first time does not get bogged down by the sheer number of files 
that have accumulated over time inside a shallow depth src directory.

>> julz generate// generate julia files

The way to alleviate this problem is to treat the src folder and the test 
folder as boxes that hold smaller boxes of like content. In a Julia 
context, this means that the src and test folder each have corresponding 
subdirectories for: methods, types, functions, macros, etc. Each 
subdirectory then has an associated generator command that can be used to 
make files of that format. For example, generating a SymArrow type with two 
array fields from the command line would use boiler-plate templates to 
compile:

   - 
   
   a working sym_arrow.jl file in the ‘src/types’ directory 
   - 
   
   an associated sym_arrow_test.jl in the ‘test/types’ directory
   

>> julz test  // test julz project

This may not seem like a big idea at first, but getting every project to 
adopt this level-deeper folder structure is the cornerstone to a successful 
framework. With it, packages can be created that expand on the native 
collection of data structures, i.e. to produce singletons, compilers, etc. 
Here, the addition of generators just make everything less user-error prone 
and promotes adequate test coverage of code.

>> julz install  // install Julia packages

One aspect of Julia that is currently lacking is a comprehensive solution 
to package management. The REQUIRE file and the currently unreleased 
“.lock” file are obviously the right ways to go, but I personally think 
reinventing the wheel for this case is a little extreme. Therefore julz 
uses a “Gemfile” and “Gemfile.lock” approach taken directly from the Ruby 
community. Imitation is the finest form of flattery and their backend 
support of package installation is second to none. Out of the box it 
supports local packages, multiple package hosts, version compatibility 
resolution, and environmental controls.

>> julz run  // run julz project

The final piece of the puzzle julz adds is a config folder for user’s to 
both alter the default conventions of the language/framework *and* separate 
environmental logic. Altering the default conventions could be something:

   - 
   
   as language specific as allowing -0 to equal +0
   - 
   
   or something more architectural, like adding a “services” subdirectory 
   to both the /src and /test folders. 
   

Environmental logic, on the other hand, is used in its web development 
context in which labels are applied to the different ways you run your 
code: 


   - 
   
   development (local machine) 
   - 
   
   test (testing suite)
   - 
   
   production (network clusters)
   

Then, decisions can be made built off this knowledge. For example, you 
could default to single precision for your test suite, double precision on 
your local box, or quadruple precision in production — all set without 
explicit reference inside your src directory.

>> julz send  // send julz project elsewhere (unimplemented)

This is where julz stands. Admittedly, it is still very much in its 
infancy. The goal therefore is to accomplish what ember-cli did for the 
ember.js community and get everyone on the same page when it comes to 
project layout and development. In this process, hopefully we can foster a 
sense of community for our users and lessen the fatigue of navigating 
through a project.

Re: [julia-users] Return structure-like elements from a function

[Dan]

That last line was probably meant to be

paths = Vector{Matrix{Float64}}(4)

On Tuesday, September 27, 2016 at 7:01:14 PM UTC-4, Yichao Yu wrote:
>
> On Tue, Sep 27, 2016 at 6:28 PM, Zhilong Liu  > wrote: 
> > Hello All, 
> > 
> > I am trying to return values to elements inside a structure-like matrix. 
> > Here is the code snippet: 
> > 
>
> What you actually meant to write is 
>
> paths  = [Array{Float64,2}() for i in 1:4] 
>
> And you might as well do 
>
> paths = Matrix{Float64}(4) 
>
> > for i = 1 : 4 
> >   paths[i] = compute_path(i) 
> > end 
> > 
> > I would like each path[i] element to be a Nx3 matrix. But I got an error 
> > saying 
> > 
> > MethodError: `convert` has no method matching 
> > convert(::Type{Array{DataType,1}}, ::Array{Float64,2}) 
> > 
> > 
> > Do anyone know how to do this properly? 
> > 
> > 
> > Thanks! 
> > 
> > Zhilong Liu 
> > 
> > 
>

[julia-users] Re: Broadcast slices

[Dan]

Nice.
It is easier when the payoffs are in vector form. My last iteration:

is_nash_equilibrium(po) = 
!reduce(|,falses(po[1]),(broadcast(<,po[i],mapslices(maximum,po[i],i)) for 
i=1:length(po)))

A one-liner :)

On Sunday, September 25, 2016 at 2:25:45 PM UTC-4, Brandon Taylor wrote:
>
> Cool! The implementation I have is:
>
> equals_max(x) = x .== maximum(x)
>
> best_response_dimension(payoff_matrix, dimension) =
> mapslices(equals_max, payoff_matrix, dimension)
>
> is_nash_equilibrium(payoffs) = @chain begin
> payoffs
> broadcast(best_response_dimension, _, 1:length(_) )
> zip(_...)
> map(all, _)
> end
>
> On Sunday, September 25, 2016 at 11:57:47 AM UTC-4, Dan wrote:
>>
>> Oops, that `cat` code was supposed to be:
>>
>> cat(1,map(x->reshape(x,1,size(x)...),array_of_array)...)
>>
>> Mew!
>>
>> On Sunday, September 25, 2016 at 11:54:43 AM UTC-4, Dan wrote:
>>>
>>> OK. So, to get the array to have the first dim as the player selector, 
>>> you can go:
>>>
>>> cat(1,map(x->reshape(1,size(x)),array_of_arrays)
>>>
>>>
>>> Anyway, keeping with the same payoff_matrix as before, I realized you 
>>> might just want a boolean array which is true if entry is a best response 
>>> (for the appropriate player according to last dim). It is the same flavor 
>>> of my previous one-liner, with `maximum` replacing `indmax` and a `.==`:
>>>
>>> isbr = payoff_matrix .== cat(nplayers+1,(mapslices(x->fill(maximum(x),
>>> size(payoff_matrix,i)), payoff_matrix[fill(:,nplayers)...,i],i) for i=1:
>>> nplayers)...)
>>>
>>> Anyway, gotta go now. Have a good one.
>>>
>>> On Sunday, September 25, 2016 at 11:46:26 AM UTC-4, Brandon Taylor wrote:
>>>>
>>>> For now, I have an array of arrays. 1 payoff array for each player. The 
>>>> arrays can be zipped to get the strategy profiles. It seems to work, but 
>>>> having everything in 1 array just seems so much more neat. Which is why I 
>>>> was looking for a neat implementation of broadcast_slices to match.
>>>>
>>>> On Sunday, September 25, 2016 at 10:53:57 AM UTC-4, Dan wrote:
>>>>>
>>>>> Have you found the right implementation?
>>>>>
>>>>> Fiddling a bit, I tend to agree with Steven G. Johnson `for` loops 
>>>>> would be the most efficient and probably the most understandable 
>>>>> implementation.
>>>>>
>>>>> Also, would it not be easier to have the first index in the 
>>>>> `payoff_matrix` determine which player's payoff are we using?
>>>>>
>>>>> Finally, following is an implementation using `mapslices` which seems 
>>>>> to work:
>>>>>
>>>>> nplayers = last(size(payoff_matrix));
>>>>>
>>>>> bestresponse = cat(nplayers+1,(mapslices(x->fill(indmax(x),size(
>>>>> payoff_matrix,i)), payoff_matrix[fill(:,nplayers)...,i],i) for i=1:
>>>>> nplayers)...)
>>>>>
>>>>> The `bestresponse` array is the same shape as `payoff_matrix`, with 
>>>>> each entry in `bestresponse[..,..,..,..,i]` denoting the strategy number 
>>>>> which is a best response to the others choices for player `i` (chosen in 
>>>>> the last index). The other player's strategies are determined by all the 
>>>>> `..,...,..` indices before, with the choice of player `i` immaterial 
>>>>> (since 
>>>>> a single best response is chosen by the `indmax` function.
>>>>>
>>>>> This is a good exercise, perhaps another question on Stackoverflow 
>>>>> would yield interesting variations.   
>>>>>
>>>>> On Saturday, September 24, 2016 at 9:40:54 PM UTC-4, Brandon Taylor 
>>>>> wrote:
>>>>>>
>>>>>> Or I guess that should be
>>>>>>
>>>>>> broadcast_slices(best_response_dimension, player_dimension, 
>>>>>> payoff_matrix, players)
>>>>>>
>>>>>> On Saturday, September 24, 2016 at 9:38:55 PM UTC-4, Brandon Taylor 
>>>>>> wrote:
>>>>>>>
>>>>>>> I guess, but I'm trying to write a generic program where I don't 
>>>>>>> know the size of the array? I'm trying to find Nash Equilibrium for an 
>>>>>>> n 
>>>>>>> dimensional array, where the player strategies are along dimensions 
>>>>>>> 1:n-1, 
>>>>>>> and the players are along dimension n. So:
>>>>>>>
>>>>>>> equals_max(x) = x .== maximum(x)
>>>>>>>
>>>>>>> best_response_dimension(payoff_matrix, dimension) =
>>>>>>> mapslices(equals_max, payoff_matrix, dimension)
>>>>>>>
>>>>>>> I'd want to do something like this:
>>>>>>>
>>>>>>> player_dimension = ndims(payoff_matrix)
>>>>>>> other_dimensions = repeat([1], inner = player_dimension - 1)
>>>>>>> number_of_players = size(payoff_matrix)[player_dimension]
>>>>>>>
>>>>>>>
>>>>>>> players = reshape(1:number_of_players, other_dimensions..., 
>>>>>>> number_of_players)
>>>>>>>
>>>>>>> broadcast_slices(best_response_dimension, payoff_matrix, players)
>>>>>>>
>>>>>>> On Thursday, September 22, 2016 at 9:00:51 PM UTC-4, Steven G. 
>>>>>>> Johnson wrote:
>>>>>>>>
>>>>>>>> At some point, it is simpler to just write loops than to try and 
>>>>>>>> express a complicated operation in terms of higher-order functions 
>>>>>>>> like 
>>>>>>>> broadcast.
>>>>>>>>
>>>>>>>

[julia-users] Re: Broadcast slices

[Dan]

Oops, that `cat` code was supposed to be:

cat(1,map(x->reshape(x,1,size(x)...),array_of_array)...)

Mew!

On Sunday, September 25, 2016 at 11:54:43 AM UTC-4, Dan wrote:
>
> OK. So, to get the array to have the first dim as the player selector, you 
> can go:
>
> cat(1,map(x->reshape(1,size(x)),array_of_arrays)
>
>
> Anyway, keeping with the same payoff_matrix as before, I realized you 
> might just want a boolean array which is true if entry is a best response 
> (for the appropriate player according to last dim). It is the same flavor 
> of my previous one-liner, with `maximum` replacing `indmax` and a `.==`:
>
> isbr = payoff_matrix .== cat(nplayers+1,(mapslices(x->fill(maximum(x),size
> (payoff_matrix,i)), payoff_matrix[fill(:,nplayers)...,i],i) for i=1:
> nplayers)...)
>
> Anyway, gotta go now. Have a good one.
>
> On Sunday, September 25, 2016 at 11:46:26 AM UTC-4, Brandon Taylor wrote:
>>
>> For now, I have an array of arrays. 1 payoff array for each player. The 
>> arrays can be zipped to get the strategy profiles. It seems to work, but 
>> having everything in 1 array just seems so much more neat. Which is why I 
>> was looking for a neat implementation of broadcast_slices to match.
>>
>> On Sunday, September 25, 2016 at 10:53:57 AM UTC-4, Dan wrote:
>>>
>>> Have you found the right implementation?
>>>
>>> Fiddling a bit, I tend to agree with Steven G. Johnson `for` loops would 
>>> be the most efficient and probably the most understandable implementation.
>>>
>>> Also, would it not be easier to have the first index in the 
>>> `payoff_matrix` determine which player's payoff are we using?
>>>
>>> Finally, following is an implementation using `mapslices` which seems to 
>>> work:
>>>
>>> nplayers = last(size(payoff_matrix));
>>>
>>> bestresponse = cat(nplayers+1,(mapslices(x->fill(indmax(x),size(
>>> payoff_matrix,i)), payoff_matrix[fill(:,nplayers)...,i],i) for i=1:
>>> nplayers)...)
>>>
>>> The `bestresponse` array is the same shape as `payoff_matrix`, with each 
>>> entry in `bestresponse[..,..,..,..,i]` denoting the strategy number which 
>>> is a best response to the others choices for player `i` (chosen in the last 
>>> index). The other player's strategies are determined by all the `..,...,..` 
>>> indices before, with the choice of player `i` immaterial (since a single 
>>> best response is chosen by the `indmax` function.
>>>
>>> This is a good exercise, perhaps another question on Stackoverflow would 
>>> yield interesting variations.   
>>>
>>> On Saturday, September 24, 2016 at 9:40:54 PM UTC-4, Brandon Taylor 
>>> wrote:
>>>>
>>>> Or I guess that should be
>>>>
>>>> broadcast_slices(best_response_dimension, player_dimension, 
>>>> payoff_matrix, players)
>>>>
>>>> On Saturday, September 24, 2016 at 9:38:55 PM UTC-4, Brandon Taylor 
>>>> wrote:
>>>>>
>>>>> I guess, but I'm trying to write a generic program where I don't know 
>>>>> the size of the array? I'm trying to find Nash Equilibrium for an n 
>>>>> dimensional array, where the player strategies are along dimensions 
>>>>> 1:n-1, 
>>>>> and the players are along dimension n. So:
>>>>>
>>>>> equals_max(x) = x .== maximum(x)
>>>>>
>>>>> best_response_dimension(payoff_matrix, dimension) =
>>>>> mapslices(equals_max, payoff_matrix, dimension)
>>>>>
>>>>> I'd want to do something like this:
>>>>>
>>>>> player_dimension = ndims(payoff_matrix)
>>>>> other_dimensions = repeat([1], inner = player_dimension - 1)
>>>>> number_of_players = size(payoff_matrix)[player_dimension]
>>>>>
>>>>>
>>>>> players = reshape(1:number_of_players, other_dimensions..., 
>>>>> number_of_players)
>>>>>
>>>>> broadcast_slices(best_response_dimension, payoff_matrix, players)
>>>>>
>>>>> On Thursday, September 22, 2016 at 9:00:51 PM UTC-4, Steven G. Johnson 
>>>>> wrote:
>>>>>>
>>>>>> At some point, it is simpler to just write loops than to try and 
>>>>>> express a complicated operation in terms of higher-order functions like 
>>>>>> broadcast.
>>>>>>
>>>>>

[julia-users] Re: Broadcast slices

[Dan]

OK. So, to get the array to have the first dim as the player selector, you 
can go:

cat(1,map(x->reshape(1,size(x)),array_of_arrays)


Anyway, keeping with the same payoff_matrix as before, I realized you might 
just want a boolean array which is true if entry is a best response (for 
the appropriate player according to last dim). It is the same flavor of my 
previous one-liner, with `maximum` replacing `indmax` and a `.==`:

isbr = payoff_matrix .== cat(nplayers+1,(mapslices(x->fill(maximum(x),size(
payoff_matrix,i)), payoff_matrix[fill(:,nplayers)...,i],i) for i=1:nplayers
)...)

Anyway, gotta go now. Have a good one.

On Sunday, September 25, 2016 at 11:46:26 AM UTC-4, Brandon Taylor wrote:
>
> For now, I have an array of arrays. 1 payoff array for each player. The 
> arrays can be zipped to get the strategy profiles. It seems to work, but 
> having everything in 1 array just seems so much more neat. Which is why I 
> was looking for a neat implementation of broadcast_slices to match.
>
> On Sunday, September 25, 2016 at 10:53:57 AM UTC-4, Dan wrote:
>>
>> Have you found the right implementation?
>>
>> Fiddling a bit, I tend to agree with Steven G. Johnson `for` loops would 
>> be the most efficient and probably the most understandable implementation.
>>
>> Also, would it not be easier to have the first index in the 
>> `payoff_matrix` determine which player's payoff are we using?
>>
>> Finally, following is an implementation using `mapslices` which seems to 
>> work:
>>
>> nplayers = last(size(payoff_matrix));
>>
>> bestresponse = cat(nplayers+1,(mapslices(x->fill(indmax(x),size(
>> payoff_matrix,i)), payoff_matrix[fill(:,nplayers)...,i],i) for i=1:
>> nplayers)...)
>>
>> The `bestresponse` array is the same shape as `payoff_matrix`, with each 
>> entry in `bestresponse[..,..,..,..,i]` denoting the strategy number which 
>> is a best response to the others choices for player `i` (chosen in the last 
>> index). The other player's strategies are determined by all the `..,...,..` 
>> indices before, with the choice of player `i` immaterial (since a single 
>> best response is chosen by the `indmax` function.
>>
>> This is a good exercise, perhaps another question on Stackoverflow would 
>> yield interesting variations.   
>>
>> On Saturday, September 24, 2016 at 9:40:54 PM UTC-4, Brandon Taylor wrote:
>>>
>>> Or I guess that should be
>>>
>>> broadcast_slices(best_response_dimension, player_dimension, 
>>> payoff_matrix, players)
>>>
>>> On Saturday, September 24, 2016 at 9:38:55 PM UTC-4, Brandon Taylor 
>>> wrote:
>>>>
>>>> I guess, but I'm trying to write a generic program where I don't know 
>>>> the size of the array? I'm trying to find Nash Equilibrium for an n 
>>>> dimensional array, where the player strategies are along dimensions 1:n-1, 
>>>> and the players are along dimension n. So:
>>>>
>>>> equals_max(x) = x .== maximum(x)
>>>>
>>>> best_response_dimension(payoff_matrix, dimension) =
>>>> mapslices(equals_max, payoff_matrix, dimension)
>>>>
>>>> I'd want to do something like this:
>>>>
>>>> player_dimension = ndims(payoff_matrix)
>>>> other_dimensions = repeat([1], inner = player_dimension - 1)
>>>> number_of_players = size(payoff_matrix)[player_dimension]
>>>>
>>>>
>>>> players = reshape(1:number_of_players, other_dimensions..., 
>>>> number_of_players)
>>>>
>>>> broadcast_slices(best_response_dimension, payoff_matrix, players)
>>>>
>>>> On Thursday, September 22, 2016 at 9:00:51 PM UTC-4, Steven G. Johnson 
>>>> wrote:
>>>>>
>>>>> At some point, it is simpler to just write loops than to try and 
>>>>> express a complicated operation in terms of higher-order functions like 
>>>>> broadcast.
>>>>>
>>>>

[julia-users] Re: Broadcast slices

[Dan]

Have you found the right implementation?

Fiddling a bit, I tend to agree with Steven G. Johnson `for` loops would be 
the most efficient and probably the most understandable implementation.

Also, would it not be easier to have the first index in the `payoff_matrix` 
determine which player's payoff are we using?

Finally, following is an implementation using `mapslices` which seems to 
work:

nplayers = last(size(payoff_matrix));

bestresponse = cat(nplayers+1,(mapslices(x->fill(indmax(x),size(
payoff_matrix,i)), payoff_matrix[fill(:,nplayers)...,i],i) for i=1:nplayers
)...)

The `bestresponse` array is the same shape as `payoff_matrix`, with each 
entry in `bestresponse[..,..,..,..,i]` denoting the strategy number which 
is a best response to the others choices for player `i` (chosen in the last 
index). The other player's strategies are determined by all the `..,...,..` 
indices before, with the choice of player `i` immaterial (since a single 
best response is chosen by the `indmax` function.

This is a good exercise, perhaps another question on Stackoverflow would 
yield interesting variations.   

On Saturday, September 24, 2016 at 9:40:54 PM UTC-4, Brandon Taylor wrote:
>
> Or I guess that should be
>
> broadcast_slices(best_response_dimension, player_dimension, payoff_matrix, 
> players)
>
> On Saturday, September 24, 2016 at 9:38:55 PM UTC-4, Brandon Taylor wrote:
>>
>> I guess, but I'm trying to write a generic program where I don't know the 
>> size of the array? I'm trying to find Nash Equilibrium for an n dimensional 
>> array, where the player strategies are along dimensions 1:n-1, and the 
>> players are along dimension n. So:
>>
>> equals_max(x) = x .== maximum(x)
>>
>> best_response_dimension(payoff_matrix, dimension) =
>> mapslices(equals_max, payoff_matrix, dimension)
>>
>> I'd want to do something like this:
>>
>> player_dimension = ndims(payoff_matrix)
>> other_dimensions = repeat([1], inner = player_dimension - 1)
>> number_of_players = size(payoff_matrix)[player_dimension]
>>
>>
>> players = reshape(1:number_of_players, other_dimensions..., 
>> number_of_players)
>>
>> broadcast_slices(best_response_dimension, payoff_matrix, players)
>>
>> On Thursday, September 22, 2016 at 9:00:51 PM UTC-4, Steven G. Johnson 
>> wrote:
>>>
>>> At some point, it is simpler to just write loops than to try and express 
>>> a complicated operation in terms of higher-order functions like broadcast.
>>>
>>

[julia-users] Re: dependent types in julia?

[Dan]

Maybe the following is the form you are looking for:

julia> decomplexify{T}(::Type{Complex{T}}) = T
decomplexify (generic function with 1 method)


julia> type bar{S,T} 
   sum::S 
   sumsqr::T
   function bar(s,ss)
   if typeof(ss) != decomplexify(typeof(s))
   error("Yaiks")
   end
   new(s,ss)
   end
   end


julia> bar{Complex{Float64},Float64}(1.5+2.0im,1.0)
bar{Complex{Float64},Float64}(1.5 + 2.0im,1.0)


julia> bar{S,T}(x::S,y::T) = bar{S,T}(x,y)
bar{S,T}


julia> bar(1.5+2.0im,1.0)
bar{Complex{Float64},Float64}(1.5 + 2.0im,1.0)


The outer constructor is necessary to get the last line working. The inner 
constructor basically maintains the constraint between S and T of: T == 
Complex{S}.

On Wednesday, September 14, 2016 at 3:38:53 PM UTC-4, Neal Becker wrote:
>
> Evan Fields wrote: 
>
> > How about something like the following? 
> > 
> > type CT{T} 
> > ctsum::Complex{T} 
> > ctsumsq::T 
> > end 
> > 
>
> I'm aware that it's easier to make the type parameter the scalar type, 
> allowing writing as you show, but as a learning exercise I'd like to know 
> how Julia would go the other way. 
>
> In c++ this would be done with template metaprogramming, but as Julia 
> claims 
> to have types as 1st class objects, I thought there should be some elegant 
> way to do this. 
>
> That is, given T is Complex{U}, I need the type U so I can write 
> ctsumsq::U 
>
>

Re: [julia-users] Julia Zmq Server creates Memory Leaks

[Dan]

This did the trick. Thanks! 

[julia-users] Julia Zmq Server creates Memory Leaks

[Dan]

 

Hi, 


I am trying to create a Julia server instance by using ZMQ with the Req/Rep 
pattern. My Julia code for the Rep-server looks like this:

module JServer

# import zmq
using ZMQ
import Base.gc

# create context
ctx=Context();
# create socket
socket=Socket(ctx, REP);
# msg
msg = 0;
# repeatedly collect garbage every second
timer = Timer(gc, 1, 1);

function connect(host, port)

# connect
ZMQ.connect(socket, string("tcp://", host, ":", port));
println(string("Julia server connected to ", host, " on port ", port));
end

 

function startServer()
  # notify
  println("Waiting for request...")
  i = 1
  while i <= 2
 # invoke garbage collection
 # gc();
 # Poll For message
 msg = receive();
 println("Message received");
 # Reply
 send();
 println("Message acknowledged");
end
end

 

function receive()
   msg = ZMQ.recv(socket);
end

 

function send()
ZMQ.send(socket, Message("test request"));
end

 

function gc(event)
gc();
end

end

I can successfully start the server by runnung this script (indicating host 
and port via ARGS):

# include server module
include(string(dirname(Base.source_path()), "\JServer.jl"));
# connect server
JServer.connect(ARGS[1], ARGS[2]);
# start server
JServer.startServer();

My tests seem to be odd, however: Whenever I send large messages from 
client, i.e. ~1MB per message, Julia seems to properly keep tracking about 
memory usage: The overall process memory remains constant. But when I am 
sending small messaged (i.e. 10k per message) Julia memory usage is growing 
permanenlty. In addition response time slows down over my requests. This 
does not happen if I explicitly invoke the garbage collector inside the 
while loop, but this would slow down the replies significantly. 


I am using Julia 0.4.6.


Is my code ok? Is there a callback solution in order to be able to avoid 
the while loop?

[julia-users] Re: Unexpected slowdown

[Dan]

Try keeping the loop over all primes and breaking when `a > maxp` where 
maxp is calculated once before the loop as `maxp = floor(Int,sqrt(n))`.

On Monday, August 22, 2016 at 4:58:30 PM UTC-4, Achu wrote:
>
> I have a simple piece of code that finds me 10001 prime numbers. 
>
> function a()
> pl=[2]
> n=3
> ct=1
> while(ct<10001)
> isnprime=true
> for a in pl  
> if n%a==0
> isnprime=false
> break
> end
> end
> if isnprime
> push!(pl,n)
> ct+=1
> end
> n+=2
> end
> return pl
> end
>
> When I tweaked the code to check only prime factors less than the sqrt of 
> the number, it slowed it down by a factor of 3.
>
> function a()
> pl=[2]
> n=3
> ct=1
> while(ct<10001)
> isnprime=true
> for a in pl[pl. if n%a==0
> isnprime=false
> break
> end
> end
> if isnprime
> push!(pl,n)
> ct+=1
> end
> n+=2
> end
> return pl
> end
>
> Why is that?
>
>

[julia-users] Re: Unexpected Performance Behaviour

[Dan]

The following standalone version recovers the speed of `test2` using the 
@fastmath macro. Integer exponents have trade-offs in implementation using 
repeated squaring and multiplication and intrinsic power instructions. Not 
really sure how to control the implementation used in every instance.

function test3(N)
r = 0.234; s = 0.0
for n = 1:N
s += @fastmath r^3 + r^5
end
return s
end




On Tuesday, August 2, 2016 at 12:33:24 AM UTC-4, Christoph Ortner wrote:
>
> Below are two tests, in the first a simple polynomial is "hard-coded", in 
> the second it is passed as a function. I would expect the two to be 
> equivalent, but the second case is significantly faster. Can anybody 
> explain what is going on?  @code_warntype doesn't show anything that would 
> explain it? 
>
> function test1(N)
>
> r = 0.234; s = 0.0
> for n = 1:N
> s += r^3 + r^5
> end
> return s
> end
>
>
> function test2(N, f1)
> r = 0.234; s = 0.0
> for n = 1:N
> s += f1(r)
> end
> return s
> end
>
>
> g1(r) = r^3 + r^5
>
>
> test1(10)
> test2(10, g1)
>
>
> println("Test1: hard-coded functions")
> @time test1(1_000_000)
> @time test1(1_000_000)
> @time test1(1_000_000)
>
>
> println("Test2: pass functions")
> @time test2(1_000_000, g1)
> @time test2(1_000_000, g1)
> @time test2(1_000_000, g1)
>
>
>
>
> # $ julia5 weird_test2.jl
> # Test1: hard-coded functions
> #   0.086683 seconds (4.00 M allocations: 61.043 MB, 50.75% gc time)
> #   0.142487 seconds (4.00 M allocations: 61.035 MB, 76.91% gc time)
> #   0.025388 seconds (4.00 M allocations: 61.035 MB, 4.28% gc time)
> # Test2: pass functions
> #   0.000912 seconds (5 allocations: 176 bytes)
> #   0.000860 seconds (5 allocations: 176 bytes)
> #   0.000846 seconds (5 allocations: 176 bytes)
>
>
>
>
>
>

[julia-users] Re: Linear-Equation-System Solver("\") does not work reliably for the "Rational" Type

[Dan]

You can also try to use the pseudo-inverse $A^{+}$ in LaTeX. It is also 
called the Moore-Penrose inverse.
In this case (with variable names in thread and `AP` used for the inverse) 
it is:
 
AP = A'*inv(A*A')

And then we have ('x' is the solution to the equations):

x = AP*v

and as it should be

A*x == v

is true.
Hope this helps.
-- Dan

On Friday, July 15, 2016 at 2:50:01 PM UTC-4, Kurolong wrote:
>
> Hey Guys and Gals  ^^
> I'm having trouble with Julia. Maybe one of you can help me out
> The program i am writing requires the 
> linear-equation-system-solver(command "\") in the folder 
> "Julia-0.4.5\share\julia\base\linalg". 
> Now my problem is that i need the equation solved strictly on integers, an 
> approximate solution is of no use to me, which is why i tried working with 
> the "Rational" Type. 
> Curiously, the solver seems to work for the Rational Type if and only if 
> no pivoting is required. Otherwise, the following Error is thrown:
>
> WARNING: pivoting only implemented for Float32, Float64, Complex64 and 
> Complex128
> ERROR: LoadError: OverflowError()
>  in * at rational.jl:188
>  [inlined code] from linalg/generic.jl:471
>  in reflector! at no file:0
>  in A_ldiv_B! at linalg/qr.jl:346
>  in \ at linalg/qr.jl:398
>  in \ at linalg/dense.jl:450
>  in include at boot.jl:261
>  in include_from_node1 at loading.jl:320
> while loading C:\users\ omitted>\documents\julia-0.4.5\sg_project\v0.3\Test02.jl, in expression 
> starting on line 3
>
> Code to reproduce the Error:
>
> A=[[2//1,0//1] [3//1,3//2] [1//1,3//2]]
> v=[2//1,0//1]
> A\v
>
> Now i am under the impression that i could build a relatively simple 
> workaround if i knew exactly where in the developer's code the actual 
> problem is, but i am confused by the whole organization of the thing. 
> Building my own solver could certainly be done, but this will most likely 
> be ineffecient. Unfortunately, runtime is very important here, as the 
> method i am currently devising is supposed to be embedded in some very 
> expensive loops.
>
> I am a very unexperienced programmer and will likely not understand much 
> programmertalk, but i am a mathematician and will probably have little 
> problems in this respect.
>
> Maybe you have some ideas how to handle this problem.
>
> Thank you in advance for taking the time to look at my problem :)
>

[julia-users] Re: Want: "Deprecation" badge for Github

[Dan]

+1

There was also a suggestion at JuliaCon the badge could count the number of 
warnings displayed - perhaps, this would give some incentive for 
progressive cleanup/levelingup of code.

On Tuesday, July 12, 2016 at 1:15:24 PM UTC-4, Erik Schnetter wrote:
>
> Our Julia projects live on Github, and many projects proudly present 
> badges for Travis, Appveyor, Codecov, and probably a few more. TL;DR: There 
> should be a badge for "no compatibility warnings" as well.
>
> In the past weeks, Julia 0.5 deprecated a certain number of features. 
> These continue to work fine, but lead to long trails of "WARNING: XYZ is 
> deprecated". Usually, some package depends on some other packages, and one 
> of these isn't using Compat in the right way yet, and the result is a long 
> string of "Base.OS_NAME is deprecated in deps/build.jl:21", repeated 27 
> times, in some package for which you have only a tangential interest.
>
> I've occasionally began to fix this for packages where I contribute (HDF5, 
> MPI), and for packages they use (BinDeps, Compat), but truth be told, 
> there's a lot of these warnings, and it's not always easy to figure out 
> where they come from, or whether they are already fixed on the package's 
> master branch, which just hasn't been tagged yet.
>
> So -- here is my idea: Create a badge for "my project doesn't show any 
> compatibility warnings for the Julia master branch". I'd happily add it to 
> all my projects.
>
> -erik
>
> -- 
> Erik Schnetter > 
> http://www.perimeterinstitute.ca/personal/eschnetter/
>

[julia-users] Re: caching the pivot ordering for a sparse cholesky factorization?

[Dan]

OK, details.

Since SuiteSparse is part of the Julia dependencies, if you compile Julia 
from source (recommended) then you already have KLU and SuiteSparse 
usually. To use KLU we need to make it into a shared library (like cholmod 
is). I've patched some makefiles to do this (just 3 makefiles).

After installing the klu shared library it can be ccalled. I've put 
together a test which basically replicates klu_simple.c from the tests in 
KLU directory of SuiteSparse.

Finally, added a little benchmark to see if there is actually a speedup 
when pre-analyzing a matrix. There is - on my computer it was a 5x speedup, 
but this speedup really depends on the use patterns.

The attachments include a patch file detailing the makefile changes and the 
klu_simple.jl file.

Would love to hear if you put the KLU library to use and what kind of 
speedup you get. Also, if you found a different solution it would be 
interesting to hear about it.

On Sunday, July 10, 2016 at 9:38:44 PM UTC-4, Gabriel Goh wrote:
>
> I'd love more details Dan!
>
> On Saturday, July 9, 2016 at 8:53:16 PM UTC-7, Dan wrote:
>>
>> Hi Gabe,
>>
>> SuiteSparse which comes together with Julia includes a library KLU which 
>> does sparse LU factorization. It has the option to return the fill-in 
>> reducing order it manages to find for a given sparsity pattern. With this 
>> returned structure it can subsequently factorize many matrices (and solve 
>> equations). After some tinkering, I have got it to work, so it is possible.
>>
>> If this is still an issue and using KLU fits the bill, I'll provide more 
>> details of my tinkering.
>>
>> DG 
>>
>> On Friday, July 8, 2016 at 4:40:45 AM UTC-4, Gabriel Goh wrote:
>>>
>>> Hey,
>>>
>>> I have a sequence of sparse matrix factorizations I need to do, each one 
>>> a different matrix but with the same sparsity structure. Is there a way I 
>>> can save the AMD (or any other) ordering that sparsesuite returns, it does 
>>> not need to be recomputed each time?
>>>
>>> Thanks a lot for the help!
>>>
>>> Gabe
>>>
>>

klu_simple.jl
Description: Binary data
diff --git a/Makefile b/Makefile
index f9a02c6..44227c7 100644
--- a/Makefile
+++ b/Makefile
@@ -279,7 +279,7 @@ JL_PRIVATE_LIBS += arpack
 endif
 ifeq ($(USE_SYSTEM_SUITESPARSE),0)
 ifeq ($(USE_GPL_LIBS), 1)
-JL_PRIVATE_LIBS += amd camd ccolamd cholmod colamd umfpack spqr suitesparseconfig
+JL_PRIVATE_LIBS += amd camd ccolamd cholmod klu colamd umfpack spqr suitesparseconfig
 endif
 endif
 ifeq ($(USE_SYSTEM_LLVM),0)
diff --git a/contrib/repackage_system_suitesparse4.make b/contrib/repackage_system_suitesparse4.make
index cb39df7..535a3fc 100755
--- a/contrib/repackage_system_suitesparse4.make
+++ b/contrib/repackage_system_suitesparse4.make
@@ -24,13 +24,15 @@ default:
 	mkdir -p $(build_libdir)
 	mkdir -p $(JULIAHOME)/deps/build/SuiteSparse-SYSTEM/lib
 	cd $(JULIAHOME)/deps/build/SuiteSparse-SYSTEM/lib && \
-	rm -f $(build_libdir)/lib{amd,cholmod,colamd,spqr,umfpack}.$(SHLIB_EXT)
+	rm -f $(build_libdir)/lib{amd,cholmod,btf,klu,colamd,spqr,umfpack}.$(SHLIB_EXT)
 	$(CC) -shared $(WHOLE_ARCHIVE) $(SS_LIB)/libsuitesparseconfig.a $(NO_WHOLE_ARCHIVE) -o $(build_libdir)/libsuitesparseconfig.$(SHLIB_EXT)
 	$(INSTALL_NAME_CMD)libsuitesparseconfig.$(SHLIB_EXT) $(build_libdir)/libsuitesparseconfig.$(SHLIB_EXT)
 	$(CC) -shared $(WHOLE_ARCHIVE) $(SS_LIB)/libamd.a $(NO_WHOLE_ARCHIVE) -o $(build_libdir)/libamd.$(SHLIB_EXT) $(LDFLAGS) -L$(build_libdir) -lsuitesparseconfig $(RPATH_ORIGIN)
 	$(INSTALL_NAME_CMD)libamd.$(SHLIB_EXT) $(build_libdir)/libamd.$(SHLIB_EXT)
 	$(CC) -shared $(WHOLE_ARCHIVE) $(SS_LIB)/libcolamd.a $(NO_WHOLE_ARCHIVE) -o $(build_libdir)/libcolamd.$(SHLIB_EXT) $(LDFLAGS) $(LIBBLAS) -L$(build_libdir) -lsuitesparseconfig $(RPATH_ORIGIN)
 	$(INSTALL_NAME_CMD)libcolamd.$(SHLIB_EXT) $(build_libdir)/libcolamd.$(SHLIB_EXT)
+	$(CXX) -shared $(WHOLE_ARCHIVE) $(SS_LIB)/libbtf.a $(NO_WHOLE_ARCHIVE) -o $(build_libdir)/libbtf.$(SHLIB_EXT) $(LDFLAGS) $(LIBBLAS) -L$(build_libdir) -lsuitesparseconfig -lcolamd -lccolamd -lcamd -lamd $(RPATH_ORIGIN)
+	$(CXX) -shared $(WHOLE_ARCHIVE) $(SS_LIB)/libklu.a $(NO_WHOLE_ARCHIVE) -o $(build_libdir)/libklu.$(SHLIB_EXT) $(LDFLAGS) $(LIBBLAS) -L$(build_libdir) -lsuitesparseconfig -lcolamd -lccolamd -lcamd -lamd $(RPATH_ORIGIN)
 	$(CXX) -shared $(WHOLE_ARCHIVE) $(SS_LIB)/libcholmod.a $(NO_WHOLE_ARCHIVE) -o $(build_libdir)/libcholmod.$(SHLIB_EXT) $(LDFLAGS) $(LIBBLAS) -L$(build_libdir) -lsuitesparseconfig -lcolamd -lccolamd -lcamd -lamd $(RPATH_ORIGIN)
 	$(INSTALL_NAME_CMD)libcholmod.$(SHLIB_EXT) $(build_libdir)/libcholmod.$(SHLIB_EXT)
 	$(CXX) -shared $(WHOLE_ARCHIVE) $(SS_LIB)/libumfpack.a $(NO_WHOLE_ARCHIVE) -o $(build_libdir)/libumfpack.$(SHLIB_EXT) $(LDFLAGS) $(LIBBLAS) -L$

[julia-users] Re: caching the pivot ordering for a sparse cholesky factorization?

[Dan]

Hi Gabe,

SuiteSparse which comes together with Julia includes a library KLU which 
does sparse LU factorization. It has the option to return the fill-in 
reducing order it manages to find for a given sparsity pattern. With this 
returned structure it can subsequently factorize many matrices (and solve 
equations). After some tinkering, I have got it to work, so it is possible.

If this is still an issue and using KLU fits the bill, I'll provide more 
details of my tinkering.

DG 

On Friday, July 8, 2016 at 4:40:45 AM UTC-4, Gabriel Goh wrote:
>
> Hey,
>
> I have a sequence of sparse matrix factorizations I need to do, each one a 
> different matrix but with the same sparsity structure. Is there a way I can 
> save the AMD (or any other) ordering that sparsesuite returns, it does not 
> need to be recomputed each time?
>
> Thanks a lot for the help!
>
> Gabe
>

[julia-users] Re: changing record separator to read text files

[Dan]

`readuntil` is your friend.

`readline` uses `readuntil`. In fact, it is defined as:

readuntil(s,'\n')

 
On Saturday, July 9, 2016 at 7:40:37 AM UTC-4, Fred wrote:
>
> Hi !
>
> It is often very useful to read a text files by blocks of lines, using 
> another line separator than '\n'.
>
> Especially in bio-informatics, for example DNA or Protein FASTA sequences 
> are separated by '\n>' (see uniprot.txt attached).
>
> In Perl,  it is possible to change the line separator using :
>  local $/ = '\n>'
> for example.
>
> In Julia I did not found how to do that :
>
> line by line :
>
> julia> f = open("uniprot.txt")
> IOStream()
>
> julia> readline(f, '\n>' )
> ERROR: syntax: invalid character literal
>
>
>
>
> all lines in an array (I prefer line by line because some files do not fit 
> into RAM) :
>
> julia> readdlm("uniprot.txt", '\n>' )
> ERROR: syntax: invalid character literal
>
> readdlm("uniprot.txt", '>' ) # works but does not give the expected result
>
>
> So I suppose that this feature is currently not implemented in Julia ? 
>
> Thanks in advance for your comments !
>
>
>
>

[julia-users] Re: Another nested function scope issue

[Dan]

Can you explain why it modified the one from the second block?

On Friday, July 1, 2016 at 9:58:26 AM UTC-4, FANG Colin wrote:
>
> Never mind, I get it.
>
> On the third run, the `state` variable that `counter` modified is the one 
> from the second let block.
>

[julia-users] Re: readdlm() pound sign at beginning of string is not read

[Dan]

Yes, this is the intended functionality because # is used a character which 
starts a comment in the file which is not parsed. The simplest way around 
this is to disble comments using the named parameter `comments` like this:
readdlm(file,'\t',comments=false)

You can also change the character used for comments to something else, 
using named parameter `comment_char`.

The Julia documentation about `readdlm` should explain this 
(http://docs.julialang.org/en/release-0.4/stdlib/io-network/#Base.readdlm).

Perhaps if the fields are quoted, the comment character does not have 
effect, but haven't checked this.

On Saturday, June 25, 2016 at 4:41:55 PM UTC-4, Anonymous wrote:
>
> Let's say I have a file of the form
>
> 1#hello
> 2hello
>
> Where the space between the number and the ascii string is intended to be 
> a tab, then I run
>
> my_file = readdlm(file, '\t')
>
> this reads in the first line as "".  Is this the desired functionality? 
>  It works fine with other characters like @ or &.
>

Re: [julia-users] How to make a Matrix of Matrix's?

[Dan]

[[M] [M]] works.
And is the same as Matrix{Float64}[[M] [M]]

But concur it is unintuitive.

On Thursday, June 23, 2016 at 10:28:39 PM UTC-4, Sheehan Olver wrote:
>
>
> [M,M]  will do a Vector of Matrices in 0.5.   But [M M] still does 
> concatenation.  So the question is how to do Matrices of Matrices without 
> concatinating. 
>
>
>
>
> > On 24 Jun 2016, at 12:05 PM, Lutfullah Tomak  > wrote: 
> > 
> > By default [M M] (without delimeter like , or ;) means concatenation so 
> it throws error. But I think in julia 0.5 [M, M] should do Matrix of 
> Matricies. 
>
>

[julia-users] Re: Value written to file different from value printed

[Dan]

The print to screen statement is rather unequivocal. By elimination, the -1 
in the file is the culprit. My guess is the -1 in the file is the value of 
a different variable. When writing to the file, add some text to identify 
which line is which var. Notice all vars are initialized to -1.

P.S. Initializing to -1 also gives an Int type to the vars, later replaced 
by Float. This is not recommended practice in Julia since it messes with 
type inference logic of the compiler (probably not recommended in any 
language).

On Friday, June 17, 2016 at 11:29:58 AM UTC-4, James Noeckel wrote:
>
> This is driving me crazy. In the below code, you can see the print 
> statement followed directly by a file write. The default value of the 
> variable mhullvol is -1. This is what gets written to the file, despite it 
> being assigned a value by that time, as demonstrated by the print, which 
> shows 0.827. Is there some behavior I should know about? Is the compiler 
> changing the order of execution? Even stranger, this is only happening to 
> mhullvol, and none of the other variables.
>
> include("LCIOUtils.jl")
> using MeshTools
> using LCIOUtils
> const globalscale = 0.005
>
> if length(ARGS) != 10
>   println("Usage: PROGRAM filename collection startevent endevent resX 
> resY resZ blurRadius bufferpercent thresholdPercent")
>   exit()
> end
>
> filename = ARGS[1]
> outputName = match(r".*/(.*).slcio", filename)[1]
> collectionname = ARGS[2]
> eventstart = parse(Int, ARGS[3])
> eventend = parse(Int, ARGS[4])
> resX, resY, resZ = parse(Int, ARGS[5]), parse(Int, ARGS[6]), parse(Int, 
> ARGS[7])
> stdev = parse(Int, ARGS[8])
> bufferpercent = parse(Int, ARGS[9])
> thresholdPercent = parse(Int, ARGS[10])
>
>
> type Collector
>   offset :: Int
>   densityfilename
>   ratiofilename
>   meshdensityfilename
>   energyfilename
>   mhullvolumefilename
>   Collector(num, offset) = new(offset, 
> "output/$(outputName)-$(collectionname)-events_$(eventstart)_to_$(eventend)-$(resX)x$(resY)x$(resZ)-blur_$(stdev)-buffer_$(bufferpercent)-threshold_$(thresholdPercent)-hulldensity.txt",
>   
>  
> "output/$(outputName)-$(collectionname)-events_$(eventstart)_to_$(eventend)-$(resX)x$(resY)x$(resZ)-blur_$(stdev)-buffer_$(bufferpercent)-threshold_$(thresholdPercent)-volumeratio.txt",
>   
>  
> "output/$(outputName)-$(collectionname)-events_$(eventstart)_to_$(eventend)-$(resX)x$(resY)x$(resZ)-blur_$(stdev)-buffer_$(bufferpercent)-threshold_$(thresholdPercent)-meshdensity.txt",
>   
>  
> "output/$(outputName)-$(collectionname)-events_$(eventstart)_to_$(eventend)-$(resX)x$(resY)x$(resZ)-blur_$(stdev)-buffer_$(bufferpercent)-threshold_$(thresholdPercent)-mhullvolume.txt",
>   
>  
> "output/$(outputName)-$(collectionname)-events_$(eventstart)_to_$(eventend)-MCEnergy.txt")
> end
>
> import Base.call
>
> function call(collector :: Collector, positions :: Matrix, i :: Int, 
> maxEnergy::Float64)
>   print("processing event $i")
>   mhull_energydensity = -1
>   ratio = -1
>   mesh_energydensity = -1
>   maxEnergy = -1
>   mhullvol = -1
>   if size(positions)[2] == 0
> println("no points to load!")
>   else
> grid = MeshGrid(resX, resY, resZ, positions, true, stdev, 
> bufferpercent / 100.0)
> maxEnergyDensity = maximum(grid.data)
> minEnergyDensity = minimum(grid.data)
> print(".")
>
> verts, indices = createMesh(grid.data, (grid.maxPt - grid.minPt)..., 
> grid.minPt..., (maxEnergyDensity - minEnergyDensity) / 100.0 * 
> thresholdPercent + minEnergyDensity, 1, globalscale)
> print(".")
> verts, indices = removeDoubles(verts, indices)
> print(".")
> verts_mhull, indices_mhull = convexhull(verts)
> print(".")
> vol = volume(verts, indices)
> mhullvol = volume(verts_mhull, indices_mhull)
>
> totalE = sum(positions[4, :])
> mhull_energydensity = totalE/mhullvol
> ratio = vol/mhullvol
> mesh_energydensity = totalE/vol
> print(".\n")
>
>   end
>
>   densityfile = open(collector.densityfilename, "a")
>   ratiofile = open(collector.ratiofilename, "a")
>   meshdensityfile = open(collector.meshdensityfilename, "a")
>   energyfile = open(collector.energyfilename, "a")
>   mhullvolumefile = open(collector.mhullvolumefilename, "a")
>   write(densityfile, "$(mhull_energydensity)\n")
>   write(ratiofile, "$(ratio)\n")
>   write(meshdensityfile, "$(mesh_energydensity)\n")
>   write(energyfile, "$(maxEnergy)\n")
>   println("Hull volume: $(mhullvol)")
>   write(mhullvolumefile, "$(mhullvol)\n")
>
>   close(densityfile)
>   close(ratiofile)
>   close(meshdensityfile)
>   close(energyfile)
>   close(mhullvolumefile)
> end
>
> mappositions(Collector(eventend - eventstart + 1, eventstart-1), filename, 
> eventstart, eventend)
>
> As you can see, the print(mhullvol) and write("$(mhullvol)\n") are right 
> next to each other, so the discrepancy makes no 

Re: [julia-users] Abstract version of rng=1:end ?

[Dan]

A reason such 'extended' ranges might be good, is the ability to dispatch 
on them for efficiency. For example, `deleteat!(vec,5:end)`  could be 
implemented faster than just any `deleteat!` of any range. This would be 
applicable to other such structures with edge effects. Indeed, we can learn 
from mathematics which often embraced infinity for ease of use and 
expression.
Again, this complication might not be worth the benefits.

So how about a half-bounded types of the form: 
typeof(5:end) == ExtendedUnitRange{BU}(5) # the BU stands for 
Bounded-Unbounded

a UnboundedUnbounded Unit range could be like `:` meaning unbounded in both 
directions.

To summarize: there are algorithmic optimizations which are enabled by 
knowing the operated range spans to the end (or even just up close to the 
end by 
> On Mon, Jun 13, 2016 at 10:47 AM, Matthew Pearce  > wrote: 
> > Hello 
> > 
> > I find myself frequently wishing to define functions that accept a range 
> > argument to operate on like: 
> > 
> > function foo{T}(M::Array{T,1}, rng::UnitRange) 
> > return sum(M[rng].^2) 
> > end 
> > 
> > It would be nice to be able to pass in `rng = 1:end`, but this isn't 
> legal. 
> > So I need a signature like, 
>
> rng=1:endof(M) 
>
> > 
> > function foo{T}(M::Array{T,1}) 
> > foo(M, 1:length(M)) 
> > end 
> > 
> > Is there a way to set `rng` to achieve my intended effect without having 
> to 
> > resort to declaring additional functions? 
> > 
> > Cheers 
> > 
> > Matthew 
>

[julia-users] Re: local package, howto?

[Dan]

Not sure, but I may have had a similar problem and doing (in the shell):

git config --global url."https://".insteadOf git://

solved it.

On Sunday, June 12, 2016 at 1:42:13 PM UTC-4, Andreas Lobinger wrote:
>
> Hello colleagues,
>
> like with similar problems, i don't know how i ended up there, but my 
> question is how to avoid
> Being asked for githup authentification for a package that exists only 
> locally (i.e. has no remote on github or similar)?
>
> julia v0.4.5 was not updated for long, but yesterday i ran a Pkg.update() 
> and today with a Pkg.update i'm asked for githup username and passwd?
>
> Any idea?
>
> Wishing a happy day,
>   Andreas
>

[julia-users] Re: Keeping heads and tails

[Dan]

Great.
This is exactly what I needed.

On Thursday, June 9, 2016 at 5:11:05 PM UTC-4, Kristoffer Carlsson wrote:
>
> Use resize! maybe.

[julia-users] Keeping heads and tails

[Dan]

Hello All,

Given a vector (v) of length n, if I want to truncate it to its first k 
elements. How does one do this?

My code currently uses `deleteat!(v,(k+1):length(v))`. But this is somewhat 
cumbersome and might be inefficient since `deleteat!` has to take care of 
all other ranges of elements to delete.

Would it be appropriate to add specialized function for manipulating the 
ends of a vector.
Suggested name could be `head` and `head!` or perhaps overload `first`.

Anyone, read up to here and wants to encourage this addition?

Regarding efficiency - going through the code of `deleteat!` shows that a 
few checks and an addition and a subtraction could be avoided until a call 
to `jl_array_del_end` is dispatched.

Re: [julia-users] confusion with Julia -p 2

[Dan Y]

The most funny part - I managed to reduce the running time to 3 seconds, 
rofl.

In my code I use a lot of computations of the dot(A*B*f,f)  for matrices 
A,B and vectors f filled with complex{float64}.  
But this is the same as dot(B*f,A'*f) !!

So, multiplication of matrices was the most time-consuming part. Even 
though matrices were sparse with sizes around 50x50.

Re: [julia-users] confusion with Julia -p 2

[Dan Y]

I must add that making code parallel was worth it! 
In my particular case I've got the following numbers:

21s  -  parallel code, blas_set_num_threads doesn't affect this much
32s  -  single-threaded code with blas_set_num_threads(2)
94s  -  single-threaded code with blas_set_num_threads(1)

Also parallel code used significantly lesser amount of allocations  - 17MB 
vs 8.6GB (which looks strange to me, I guess allocations within workers 
weren't counted).

So, while BLAS can utilize multiple cores to achieve decent results, "tame" 
parallel code is still considerably better. 
And it is not so hard to write it in Julia)

Re: [julia-users] confusion with Julia -p 2

[Dan Y]

Thanks! With blas_set_num_threads(CPU_CORES) the running time is equal.

On Saturday, April 23, 2016 at 9:07:06 PM UTC+3, Jason Eckstein wrote:
>
> Yes, that is correct.  If you want to force each worker to use every core 
> for linear algebra operations you need to use the following 
> line: blas_set_num_threads(CPU_CORES)
>
> I've done this myself and had parallel julia processes each using more 
> than one core.  Make use that line of code is used in each parallel 
> instance being called.
>
> On Saturday, April 23, 2016 at 4:29:40 AM UTC-6, Dan Y wrote:
>>
>> Maybe it is because I use matrix multiplication in my code. 
>> I guess that base library which is responsible for linear algebra uses 
>> multiple cores automatically.
>> But with "-p 2" it is forced to work only on 1 core for some reason.
>> It's only a guess, though.
>>
>

Re: [julia-users] confusion with Julia -p 2

[Dan Y]

Maybe it is because I use matrix multiplication in my code. 
I guess that base library which is responsible for linear algebra uses 
multiple cores automatically.
But with "-p 2" it is forced to work only on 1 core for some reason.
It's only a guess, though.

[julia-users] Re: Should Iterators be nonmutating?

[Dan]

The OS resource cost of creating a new stream for each Eachline does not 
seem like a problem, since the stream is closed at the end of the iteration 
and massive parallel Eachlines on the same file are the only blocking 
usecase. Additionally, since the `done` function for Eachline closes the 
stream, it seems weird not to create a new one on `start`.

So, IOStream should be clonable in the sense you suggest (not totally sure 
I got what you meant).

Furthermore, the semantics of an IOStream is to allow stream access to an 
underlying resource (a file for example). This resource should have a 
"name", and IOStream should allow you to get the name. With this name, you 
could try to open a new IOStream to the same resource (possibly with 
different permissions).

Currently, a file IOStream has a `name` field with the filename in some 
format. But perhaps a `getURI` (or `getURL` or `name`) method should return 
a representation of the underlying resource of IOStream which could be used 
to create another IOStream for that resource.

On Saturday, April 23, 2016 at 6:36:36 AM UTC+3, Lyndon White wrote:
>
> I'm not sure changing this behavior is possible,
> at least not without making IOStream clonable.
> And I suspect there are filesystem constraints against this.
>
> Well there is a scary way, when the state just holds the current file 
> position,
> then seeks to it when next is called,
> reads a line,
> then seeks back to where it was before.
> So it looks like we are not mutating the underlying object.
> (cos we undo our changes)
> This would not be threadsafe though.
>
>
> On Saturday, 23 April 2016 05:22:11 UTC+8, Dan wrote:
>>
>> Looks like it might be appropriate to make the `stream` in `Eachline` be 
>> the state of the iterator (instead of the `nada` returned by `next`).
>>
>> Opening an issue (perhaps with a PR of a fix) could be the next step. 
>> Unless someone finds some other logic for keeping `Eachline` the way it is.
>>
>> On Friday, April 22, 2016 at 7:09:13 AM UTC+3, Lyndon White wrote:
>>>
>>> I am writing some iterators that work with `IO`.
>>>
>>> I thought, so that my iterators would not mutate,
>>> I would `deepcopy` the IO item, and store it as part of my state.
>>> (then only my state would mutate, and no-one cares probably (unless they 
>>> were trying to `tee` for free))
>>>
>>> But `deepcopying` IO items apparently does not actually cause them to 
>>> become separate.
>>> At least not for IOStreams.
>>>
>>>
>>>
>>> io=open("./data/text8/text8","r")
>>> ia = deepcopy(io)
>>> >IOStream()
>>>
>>>
>>> position(io)
>>> >0
>>>
>>> position(seekend(io))
>>> >1
>>>
>>> position(ia)
>>> >1
>>>
>>>
>>>
>>> I was under the assumption that iterators never mutated, as there is no 
>>> `!` in `next`.
>>> (then again there is also no `!` in `seekend`)
>>>
>>> So I was a bit stuck for how to implement.
>>>
>>> So I thought I would see how `Eachline` is implemented. 
>>> <https://github.com/JuliaLang/julia/blob/master/base/io.jl#L356>
>>> It is mutating.
>>> So maybe my assumptions are wrong, and iterators can be mutating?
>>> They would be a lot easier to write that way, I guess.
>>>
>>> But julia does not have a `tee` function 
>>> <https://docs.python.org/2/library/itertools.html#itertools.tee>
>>>
>>>
>>>
>>>

[julia-users] Re: Should Iterators be nonmutating?

[Dan]

Looks like it might be appropriate to make the `stream` in `Eachline` be 
the state of the iterator (instead of the `nada` returned by `next`).

Opening an issue (perhaps with a PR of a fix) could be the next step. 
Unless someone finds some other logic for keeping `Eachline` the way it is.

On Friday, April 22, 2016 at 7:09:13 AM UTC+3, Lyndon White wrote:
>
> I am writing some iterators that work with `IO`.
>
> I thought, so that my iterators would not mutate,
> I would `deepcopy` the IO item, and store it as part of my state.
> (then only my state would mutate, and no-one cares probably (unless they 
> were trying to `tee` for free))
>
> But `deepcopying` IO items apparently does not actually cause them to 
> become separate.
> At least not for IOStreams.
>
>
>
> io=open("./data/text8/text8","r")
> ia = deepcopy(io)
> >IOStream()
>
>
> position(io)
> >0
>
> position(seekend(io))
> >1
>
> position(ia)
> >1
>
>
>
> I was under the assumption that iterators never mutated, as there is no 
> `!` in `next`.
> (then again there is also no `!` in `seekend`)
>
> So I was a bit stuck for how to implement.
>
> So I thought I would see how `Eachline` is implemented. 
> 
> It is mutating.
> So maybe my assumptions are wrong, and iterators can be mutating?
> They would be a lot easier to write that way, I guess.
>
> But julia does not have a `tee` function 
> 
>
>
>
>

[julia-users] confusion with Julia -p 2

[Dan Y]

I use julia-0.4.5 on windows10 with pentium B970 (2 cores, 2 threads)

I have a code with no parallel instructions.
When I run it in a standard session with no options my CPU utilization is 
about 80-90%.
But when I run it in a session with "Julia -p 2" the CPU utilization is 
only about 20-25% .
And the code actually works about 3x times slower.

I wonder why. 
I didn't try to make code parallel yet, but I doubt it'll worth it.

[julia-users] Re: Modifying large vectors of repeated data

[Dan]

Similar `rle` and `inverse_rle` functions exists in StatsBase.jl 
https://github.com/JuliaStats/StatsBase.jl/blob/master/src/misc.jl#L15 .
These function may be enough for the application involved.

On Thursday, April 14, 2016 at 6:48:40 PM UTC+3, Kenta Sato wrote:
>
> I guess what you want is RLEVector of RLEVectors.jl: 
> https://github.com/phaverty/RLEVectors.jl 
> 
> .
> This does run-length encoding for repeated elements and hence applying a 
> function to all elements would be fast.
>
>
> julia> v = RLEVector(vcat(zeros(1), ones(1)))
> RLEVectors.RLEVector{Float64,Int64}
> run values: [0.0,1.0]
> run ends:   [1,2]
>
> julia> map!(x -> x + 1, v)
> RLEVectors.RLEVector{Float64,Int64}
> run values: [1.0,2.0]
> run ends:   [1,2]
>
>
>
> On Thursday, April 14, 2016 at 11:48:48 PM UTC+9, Terry Seaward wrote:
>>
>> Hi,
>>
>> I have a large vector of repeated data (dates in my case, but I'd like to 
>> keep this generic). In R I can apply a function to the elements efficiently 
>> using something like this:
>>
>> myfun <- function(x, fun, ...) {
>> ux <- unique(x); 
>> fun(ux, ...)[match(x, ux)]
>> }
>>
>> myfun(x, as.character)
>>
>> So basically apply the function to the unique elements then use match to 
>> recreate x. 
>>
>> What's the best way to do this in Julia?
>>
>> -- Terry
>>
>

[julia-users] Re: Fast coroutines

[Dan]

Hello Carl and welcome to the Julia community (and fan club).

Regarding your question:
1. `normal_sum` version is the one of the most optimized construct you can 
imagine. So it is good it's working so fast, even relative with a less 
common Task construct.
2. `values` seems to be a global variable (and not const) which is a source 
of trouble for Julia type-inference based optimization (it is mentioned in 
the Performance Tips). So try to wrap up things in a function or make them 
const (but make sure they are not global which cannot be relied to keep 
their types).

Dan

On Saturday, April 2, 2016 at 2:55:51 PM UTC+3, Carl wrote:
>
> Hello,
>
> Julia is great.  I'm new fan and trying to figure out how to write simple 
> coroutines that are fast.  So far I have this:
>
>
> function task_iter(vec)
> @task begin
> i = 1
> for v in vec
> produce(v)
> end
> end
> end
>
> function task_sum(vec)
> s = 0.0
> for val in task_iter(vec)
> s += val
> end
> return s
> end
>
> function normal_sum(vec)
> s = 0.0
> for val in vec
> s += val
> end
> return s
> end
>
> values = rand(10^6)
> task_sum(values)
> normal_sum(values)
>
> @time task_sum(values)
> @time normal_sum(values)
>
>
>
>   1.067081 seconds (2.00 M allocations: 30.535 MB, 1.95% gc time)
>   0.006656 seconds (5 allocations: 176 bytes)
>
> I was hoping to be able to get the speeds to match (as close as possible). 
>  I've read the performance tips and I can't find anything I'm doing wrong. 
>  I also tried out 0.5 thinking that maybe it would be faster with 
> supporting fast anonymous functions but it was slower (1.5 seconds).
>
>
> Carl
>
>

[julia-users] Re: How do I, as efficiently as possible, iterate over the first and last elements along a given dimension of an array?

[Dan]

Sorry for the bad formatting (getting tangled with the editor).

Additional useful code to supplement the previous note:

function dropdim(A,i)
  v = collect(size(A))
  splice!(A,i)
  return tuple(v...)
end

# now something like this works:
# which gives the sum of the two edge planes of array A3 as a matrix 
reshape(map(x->A3[x[1]]+A3[x[2]],mapedges(A3,3)),dropdim(A3,3))






On Thursday, March 24, 2016 at 3:21:02 AM UTC+2, Dan wrote:
>
> # something along these lines:
>
> topleft(A,d) = tuple(ones(Int,ndims(A))...)
> bottomleft(A,d) = tuple([i==d ? 1 : e for (i,e) in enumerate(size(A))]...)
>
> topright(A,d) = tuple([i==d ? e : 1 for (i,e) in enumerate(size(A))]...)
> bottomright(A,d) = size(A)
> mapedges(A,d) = zip(
>   CartesianRange{CartesianIndex{ndims(A)}}(
> CartesianIndex{ndims(A)}(topleft(A,d)...),
> CartesianIndex{ndims(A)}(bottomleft(A,d)...)),
>   CartesianRange{CartesianIndex{ndims(A)}}(
> CartesianIndex{ndims(A)}(topright(A,d)...),
> CartesianIndex{ndims(A)}(bottomright(A,d)...)
>   ))
>
> A3 = rand(10,15,8)
>
> julia> mapedges(A3,1) |> length
> 120
>
> julia> mapedges(A3,2) |> length
> 80
>
> julia> mapedges(A3,3) |> length
> 150
>
>
> On Thursday, March 24, 2016 at 2:53:48 AM UTC+2, Tomas Lycken wrote:
>>
>> …but not really. Reading the docstring more carefully:
>>
>> Transform the given dimensions of array A using function f. f is called 
>> on each slice of A of the form A[…,:,…,:,…]. dims is an integer vector 
>> specifying where the colons go in this expression. The results are 
>> concatenated along the remaining dimensions. For example, if dims is [1,2] 
>> and A is 4-dimensional, f is called on A[:,:,i,j] for all i and j.
>>
>> What I want to do, is rather call f on A[:,:,1,:] and A[:,:,end,:], but 
>> nothing in between 1 and end for that dimension. mapslices still 
>> eventually visit the entire array (either by slicing, or by iteration), but 
>> I only want to visit the “edges”. I might be missing something, though.
>>
>> // T
>>
>> On Thursday, March 24, 2016 at 1:48:36 AM UTC+1, Tomas Lycken wrote:
>>
>> Yes, probably - thanks for the tip! I'll see if I can cook something up...
>>>
>>> On Thursday, March 24, 2016 at 1:45:32 AM UTC+1, Benjamin Deonovic wrote:
>>>>
>>>> Can mapslices help here?
>>>>
>>>>
>>>> On Wednesday, March 23, 2016 at 6:59:59 PM UTC-5, Tomas Lycken wrote:
>>>>>
>>>>> Is there an effective pattern to iterate over the “endpoints” of an 
>>>>> array along a given dimension?
>>>>>
>>>>> What I eventually want to accomplish is to apply a function (in this 
>>>>> case an equality test) to the two end points along a particular dimension 
>>>>> of an array. I think the pattern is easiest explained by considering 1D, 
>>>>> 2D 
>>>>> and 3D:
>>>>>
>>>>> # assume the existence of some scalar-valued function f(x,y)
>>>>>
>>>>> A1 = rand(10)
>>>>> f(A1[1], A1[end]) # d == 1 (the only possible value) -> one evaluation
>>>>>
>>>>> A2 = rand(10, 15)
>>>>> map(f, A2[1,:], A2[end,:]) # d == 1 -> 15 evaluations
>>>>> map(f, A2[:,1], A2[:,end]) # d == 2 -> 10 evaluations
>>>>>
>>>>> A3 = rand(10, 15, 8)
>>>>> map(f, A3[1,:,:], A3[end,:,:]) # d == 1 -> 15x8 evaluations
>>>>> map(f, A3[:,1,:], A3[:,end,:]) # d == 2 -> 10x8 evaluations
>>>>> map(f, A3[:,:,1], A3[:,:,end]) # d == 3 -> 10x15 evaluations
>>>>>
>>>>> I just want to consider one dimension at a time, so given A and d, 
>>>>> and in this specific use case I don’t need to collect the results, so a 
>>>>> for-loop without an allocated place for the answer instead of a map 
>>>>> is just fine (probably preferrable, but it’s easier to go in that 
>>>>> direction 
>>>>> than in the other). What I’m struggling with, is how to generally 
>>>>> formulate 
>>>>> the indexing expressions (like [, 1, >>>> instances of :>], but not in pseudo-code…). I assume this can be done 
>>>>> somehow using CartesianIndexes and/or CartesianRanges, but I can’t 
>>>>> get my mind around to how. Any help is much appreciated.
>>>>>
>>>>> // T
>>>>> ​
>>>>>
>>>> ​
>>
>

[julia-users] Re: How do I, as efficiently as possible, iterate over the first and last elements along a given dimension of an array?

[Dan]

# something along these lines:

topleft(A,d) = tuple(ones(Int,ndims(A))...)
bottomleft(A,d) = tuple([i==d ? 1 : e for (i,e) in enumerate(size(A))]...)

topright(A,d) = tuple([i==d ? e : 1 for (i,e) in enumerate(size(A))]...)
bottomright(A,d) = size(A)
mapedges(A,d) = zip(
  CartesianRange{CartesianIndex{ndims(A)}}(
CartesianIndex{ndims(A)}(topleft(A,d)...),
CartesianIndex{ndims(A)}(bottomleft(A,d)...)),
  CartesianRange{CartesianIndex{ndims(A)}}(
CartesianIndex{ndims(A)}(topright(A,d)...),
CartesianIndex{ndims(A)}(bottomright(A,d)...)
  ))

A3 = rand(10,15,8)

julia> mapedges(A3,1) |> length
120

julia> mapedges(A3,2) |> length
80

julia> mapedges(A3,3) |> length
150


On Thursday, March 24, 2016 at 2:53:48 AM UTC+2, Tomas Lycken wrote:
>
> …but not really. Reading the docstring more carefully:
>
> Transform the given dimensions of array A using function f. f is called on 
> each slice of A of the form A[…,:,…,:,…]. dims is an integer vector 
> specifying where the colons go in this expression. The results are 
> concatenated along the remaining dimensions. For example, if dims is [1,2] 
> and A is 4-dimensional, f is called on A[:,:,i,j] for all i and j.
>
> What I want to do, is rather call f on A[:,:,1,:] and A[:,:,end,:], but 
> nothing in between 1 and end for that dimension. mapslices still 
> eventually visit the entire array (either by slicing, or by iteration), but 
> I only want to visit the “edges”. I might be missing something, though.
>
> // T
>
> On Thursday, March 24, 2016 at 1:48:36 AM UTC+1, Tomas Lycken wrote:
>
> Yes, probably - thanks for the tip! I'll see if I can cook something up...
>>
>> On Thursday, March 24, 2016 at 1:45:32 AM UTC+1, Benjamin Deonovic wrote:
>>>
>>> Can mapslices help here?
>>>
>>>
>>> On Wednesday, March 23, 2016 at 6:59:59 PM UTC-5, Tomas Lycken wrote:

 Is there an effective pattern to iterate over the “endpoints” of an 
 array along a given dimension?

 What I eventually want to accomplish is to apply a function (in this 
 case an equality test) to the two end points along a particular dimension 
 of an array. I think the pattern is easiest explained by considering 1D, 
 2D 
 and 3D:

 # assume the existence of some scalar-valued function f(x,y)

 A1 = rand(10)
 f(A1[1], A1[end]) # d == 1 (the only possible value) -> one evaluation

 A2 = rand(10, 15)
 map(f, A2[1,:], A2[end,:]) # d == 1 -> 15 evaluations
 map(f, A2[:,1], A2[:,end]) # d == 2 -> 10 evaluations

 A3 = rand(10, 15, 8)
 map(f, A3[1,:,:], A3[end,:,:]) # d == 1 -> 15x8 evaluations
 map(f, A3[:,1,:], A3[:,end,:]) # d == 2 -> 10x8 evaluations
 map(f, A3[:,:,1], A3[:,:,end]) # d == 3 -> 10x15 evaluations

 I just want to consider one dimension at a time, so given A and d, and 
 in this specific use case I don’t need to collect the results, so a 
 for-loop without an allocated place for the answer instead of a map is 
 just fine (probably preferrable, but it’s easier to go in that direction 
 than in the other). What I’m struggling with, is how to generally 
 formulate 
 the indexing expressions (like [, 1, >>> instances of :>], but not in pseudo-code…). I assume this can be done 
 somehow using CartesianIndexes and/or CartesianRanges, but I can’t get 
 my mind around to how. Any help is much appreciated.

 // T
 ​

>>> ​
>

Re: [julia-users] Re: Mmap allocation question

[Dan]

Yep, `peCounters`, `paCounters` and `dims` are not type-stable. They are 
one type by their default values and then assigned another. Perhaps rename 
the default parameters, and copy them to `peCounters`, `paCounters` and 
`dims` only if they are set to something other than `0`.

Also, `mdhcounters` might not return a definite type (need to check that 
function).
Fixing these should make the loop efficient.

On Saturday, March 12, 2016 at 4:05:39 PM UTC+2, Tim Loderhose wrote:
>
> Here's the actual code: 
> https://gist.github.com/timlod/0f607e311d0464fd6c63
>
> I am running the code from the REPL, may that be a problem? (As I read in 
> the REPL everything is global). In the file nothing is global.
> Also, the counters are UInt16s, but that shouldnt matter I guess.
>
> Thanks for the help so far!
>
> On Saturday, 12 March 2016 14:22:38 UTC+1, Dan wrote:
>>
>> It's better to have code which actually runs in the post. In any case, 
>> the allocations at the `for` lines is suspicious - the for should basically 
>> only allocate a counter. Are there any global variables? Is `counter1` or 
>> `counter2` or `dims` global? Globals are always a good source of confusion 
>> to the type-inference engine.
>>
>> On Saturday, March 12, 2016 at 2:28:51 PM UTC+2, Tim Loderhose wrote:
>>>
>>> The code is in a function. I changed the names a bit to make it more 
>>> understandable. The actual function is longer and has different variable 
>>> names.
>>>
>>> On Saturday, 12 March 2016 13:01:28 UTC+1, tshort wrote:
>>>>
>>>> Is that code in a function? (It should be.) Also, one of your variable 
>>>> names changed to `counter1s`. Suspect a type instability.
>>>> On Mar 12, 2016 4:12 AM, "Tim Loderhose"  wrote:
>>>>
>>>>> I tried around with that a bit, but then it gets much worse: From ~1s 
>>>>> to ~6s, allocation as shown:
>>>>>
>>>>> 153710487 mat = Array{Complex64}(dims...)
>>>>>   4722450   file = Mmap.mmap(filename, Array{Complex64,2}, 
>>>>> (dims[2],length(counter1)))
>>>>>  9568  for i = 1:dims[2]
>>>>>  4000 for j = 1:length(counter1)
>>>>> 1690462534  mat[counter1s[j],i,counter2[j]] = file[i,j]
>>>>> - end
>>>>>
>>>>> I swapped the for loops around here, but that didn't matter. I can 
>>>>> gain a little bit by indexing i into the first dimension of mat, but it 
>>>>> still lags far behind.
>>>>> Any other ideas?
>>>>>
>>>>> On Saturday, 12 March 2016 03:15:33 UTC+1, Greg Plowman wrote:
>>>>>>
>>>>>> I think array slices (on right hand side of assignment) create new 
>>>>>> arrays, hence the allocation.
>>>>>> Try writing an explicit loop instead, something like:
>>>>>>
>>>>>> for j = 1:length(counter1)
>>>>>>for i = 1:size(file,1)
>>>>>>mat[counter1[j],i,counter2[j]] = file[i,j]
>>>>>>end
>>>>>> end
>>>>>>
>>>>>>
>>>>>> On Saturday, March 12, 2016 at 12:25:00 PM UTC+11, Tim Loderhose 
>>>>>> wrote:
>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> I have a question regarding some allocation in my code I would like 
>>>>>>> to get rid of.
>>>>>>> I am memory mapping a file (which could be very large) which is part 
>>>>>>> of a complex 3D matrix, and then put its contents into the preallocated 
>>>>>>> matrix along the second dimension. I need the counters because the 
>>>>>>> contents 
>>>>>>> of file are only a subset of the full matrix.
>>>>>>>
>>>>>>> Here's a profiled snippet, where the file which is loaded has 
>>>>>>> 120619520 bytes.
>>>>>>>
>>>>>>> 153705063 mat = Array{Complex64}(dims...)
>>>>>>>  4721282file = Mmap.mmap(filename, Array{Complex64,2}, 
>>>>>>> (dims[2],length(counter1)))
>>>>>>> 16   for i = 1:length(counter1)
>>>>>>> 148179531   mat[counter1[i],:,counter2[i]] = file[:,i]
>>>>>>> -  end
>>>>>>>
>>>>>>> Why does the code allocate so much memory inside the for-loop (even 
>>>>>>> more bytes than the contents of file)?
>>>>>>> It seems like this is a trivial matter, right now I just can't get 
>>>>>>> my head around it, any help is appreciated :)
>>>>>>>
>>>>>>> Thanks,
>>>>>>> Tim
>>>>>>>
>>>>>>

Re: [julia-users] Re: Mmap allocation question

[Dan]

It's better to have code which actually runs in the post. In any case, the 
allocations at the `for` lines is suspicious - the for should basically 
only allocate a counter. Are there any global variables? Is `counter1` or 
`counter2` or `dims` global? Globals are always a good source of confusion 
to the type-inference engine.

On Saturday, March 12, 2016 at 2:28:51 PM UTC+2, Tim Loderhose wrote:
>
> The code is in a function. I changed the names a bit to make it more 
> understandable. The actual function is longer and has different variable 
> names.
>
> On Saturday, 12 March 2016 13:01:28 UTC+1, tshort wrote:
>>
>> Is that code in a function? (It should be.) Also, one of your variable 
>> names changed to `counter1s`. Suspect a type instability.
>> On Mar 12, 2016 4:12 AM, "Tim Loderhose"  wrote:
>>
>>> I tried around with that a bit, but then it gets much worse: From ~1s to 
>>> ~6s, allocation as shown:
>>>
>>> 153710487 mat = Array{Complex64}(dims...)
>>>   4722450   file = Mmap.mmap(filename, Array{Complex64,2}, 
>>> (dims[2],length(counter1)))
>>>  9568  for i = 1:dims[2]
>>>  4000 for j = 1:length(counter1)
>>> 1690462534  mat[counter1s[j],i,counter2[j]] = file[i,j]
>>> - end
>>>
>>> I swapped the for loops around here, but that didn't matter. I can gain 
>>> a little bit by indexing i into the first dimension of mat, but it still 
>>> lags far behind.
>>> Any other ideas?
>>>
>>> On Saturday, 12 March 2016 03:15:33 UTC+1, Greg Plowman wrote:

 I think array slices (on right hand side of assignment) create new 
 arrays, hence the allocation.
 Try writing an explicit loop instead, something like:

 for j = 1:length(counter1)
for i = 1:size(file,1)
mat[counter1[j],i,counter2[j]] = file[i,j]
end
 end


 On Saturday, March 12, 2016 at 12:25:00 PM UTC+11, Tim Loderhose wrote:

> Hi,
>
> I have a question regarding some allocation in my code I would like to 
> get rid of.
> I am memory mapping a file (which could be very large) which is part 
> of a complex 3D matrix, and then put its contents into the preallocated 
> matrix along the second dimension. I need the counters because the 
> contents 
> of file are only a subset of the full matrix.
>
> Here's a profiled snippet, where the file which is loaded has 
> 120619520 bytes.
>
> 153705063 mat = Array{Complex64}(dims...)
>  4721282file = Mmap.mmap(filename, Array{Complex64,2}, 
> (dims[2],length(counter1)))
> 16   for i = 1:length(counter1)
> 148179531   mat[counter1[i],:,counter2[i]] = file[:,i]
> -  end
>
> Why does the code allocate so much memory inside the for-loop (even 
> more bytes than the contents of file)?
> It seems like this is a trivial matter, right now I just can't get my 
> head around it, any help is appreciated :)
>
> Thanks,
> Tim
>

[julia-users] Re: what name would you give to the functional 'complement' of trunc()?

[Dan]

trunc_to_inf and add an alias trunc_to_zero for regular trunc ??

On Saturday, March 12, 2016 at 1:49:41 PM UTC+2, Jeffrey Sarnoff wrote:
>
> x = trunc(fp) yields fp rounded towards zero to the closest integer valued 
> float.
> What name would you give to the complementary function, rounding away from 
> zero to the closest integer valued float?
>
> stretch(fp)  = signbit(fp) ? floor(fp) : ceil(fp)   # abs(trunc(fp) - 
> stretch(fp)) is 0.0 when fp is integer valued and is 1.0 otherwise
>
>
>
>

[julia-users] Re: What wrong in (D'*D)/k.-(mean(D,1)'*mean(D,1)) ?

[Dan]

The matrix `D` is based on type `Int8` which overflows and wraps (without 
error) very quickly. The `cov` function automatically converts to a `Float64` 
but the manual calculation does not. To make the calculations match, try:
D = convert(Matrix{Int},D)

On Friday, March 4, 2016 at 10:42:02 AM UTC+2, paul analyst wrote:
>
> Why (D'*D)/k.-(mean(D,1)'*mean(D,1)) compute var = -0.018 if var must 
> be 0.124326 ?**
>
>
> julia> using HDF5, JLD
>
> julia> D=load("D_test.jld","D_test");
>
> julia> k,l=size(D)
> (100,10)
>
> julia> cov(D[:,1],D[:,1],corrected=false)
> 0.12432634086422582
>
> julia> cov(D,corrected=false)
> 10x10 Array{Float64,2}:
> 0.124326 0.0239205 -0.0119771 -0.0201163 0.013
> 0.0239205 0.154538 -0.0156671 -0.0263335 -0.023
> -0.0119771 -0.0156671 0.0853087 -0.0158088 0.002
> -
>
> julia> (D'*D)/k.-(mean(D,1)'*mean(D,1))
> 10x10 Array{Float64,2}:
> -0.018 -0.0244635 -0.0119771 -0.0247243 -0.01791
> -0.0244635 -0.0320863 -0.0156671 -0.0324775 -0.02332
> -0.0119771 -0.0156671 -0.00761925 -0.0158088 -0.01140
>
> julia> (D[:,1]'*D[:,1])/k.-(mean(D[:,1],1)'*mean(D[:,1],1))
> 1-element Array{Float64,1}:
> 0.124326
>
> julia> mean(D,1)'
> 10x1 Array{Float64,2}:
> 0.136944
> 0.179135
> 0.08746
> 0.180755
> 0.130206
> 0.015541
> 0.036076
> 0.101094
> 0.155723
> 0.045505
>
> julia> mean(D,1)
> 1x10 Array{Float64,2}:
> 0.136944 0.179135 0.08746 0.180755 0.130206 0.015541 0
>
> julia> ee=mean(D[:,1],1)'*mean(D[:,1],1)
> 1-element Array{Float64,1}:
> 0.0187537
>
> julia> dd=(D[:,1]'*D[:,1])/k
> 1-element Array{Float64,1}:
> 0.14308
>
> julia> dd-ee
> 1-element Array{Float64,1}:
> 0.124326
>
>
> file: 
>
> https://drive.google.com/open?id=0B9xW5VtANWhDOHlwVlZxZk1xanc
>
>
>

Re: [julia-users] Adding element to array and filtering multi-dimensional arrays.

[Dan]

On Monday, February 22, 2016 at 12:18:02 PM UTC+2, Tamas Papp wrote:
>
>
> On Mon, Feb 22 2016, Aleksandr Mikheev wrote: 
>
> > Hi everyone. I have two simple questions, answer to which I cannot find 
> by 
> > myself. 
> > 
> > 1. Is there any way to add an element to array? For example, I have 
> > Array{...}(n,m) or Vector{...}(n). And I'd like to expand it, i.e. I 
> want 
> > to have now Array{n+1,m} or Vector{...}(n+1). How can I do this? 
>
> See push! and resize! for vectors. For general arrays, I would just use 
> hcat to make a new one if I am only doing this a few times, but you can 
> find messages on the lists with tricks to avoid reallocation if 
> performance is critical (grow them by more than a single row, eg 
> doubling). Alternatively, you can use a vector of vectors and 
> concatenate them when you have the result. 
>
> > 2. Imagine I have an array [1 2 3 4 5 6 7 8 9 10; 12 4 5 0 2 0 45 8 0 
> 7]. 
> > Is there any easy way to remove all columns with 0 in second line? 
> > Something like filter!(...)? 
>
> Try something like 
>
> A[:,mapslices(col->all(col.!=0),A,1)] 
>
> This could be slightly faster (not benchmarked): 
>
> function columnswithzeros{T <: Number}(A::AbstractArray{T,2}) 
> n,m = size(A) 
> z = zeros(Bool,m) 
> for col = 1:m 
> for row = 1:n 
> if A[row,col] == 0 
> z[col] = true 
> break 
> end 
> end 
> end 
> z 
> end 
>
> A[:,~columnswithzeros(A)] 
>
> Best, 
>
> Tamas 
>

For question #1: hcat(...) and cat(...) as Tamas suggested.
For question #2:
a = [1 2 3 4 5 6 7 8 9 10; 12 4 5 0 2 0 45 8 0 7]
a = a[:,a[2,:].!=0]
a

2x7 Array{Int64,2}:
  1  2  3  5   7  8  10
 12  4  5  2  45  8   7

Re: [julia-users] Adding a mtrix to the third dimension of a multi-dimensional array

[Dan]

>From the help for `cat`:
cat(dims, A...)

  Concatenate the input arrays along the specified dimensions in the 
iterable dims

And indeed, if 
size(M1)=(3,3,3) 
and 
size(M2)=(3,3)
Then,
size(cat(3,M1,M2)) = (3,3,4)

This method may not be efficient (though in terms of memory layout it could 
be).

On Saturday, February 13, 2016 at 1:06:02 PM UTC+2, Mauro wrote:
>
> I think this is not possible.  Instead use a Vector, append! to that and 
> then reshape in the end.  The reshape should result in a view and not a 
> copy, thus will be fast. 
>
> Mauro 
>
> On Sat, 2016-02-13 at 10:51, Vishnu Raj  > wrote: 
> > Hi, 
> > 
> > I have a three dimensional array (say 'z' ) like this : 
> > 
> > julia> z 
> > 4x3x3 Array{Int64,3}: 
> > [:, :, 1] = 
> >  1  5   9 
> >  2  6  10 
> >  3  7  11 
> >  4  8  12 
> > 
> > [:, :, 2] = 
> >  13  17  21 
> >  14  18  22 
> >  15  19  23 
> >  16  20  24 
> > 
> > [:, :, 3] = 
> >  25  29  33 
> >  26  30  34 
> >  27  31  35 
> >  28  32  36 
> > 
> > Now I have a matrix 'z1' as 
> > julia> z1 
> > 4x3 Array{Int64,2}: 
> >  37  41  45 
> >  38  42  46 
> >  39  43  47 
> >  40  44  48 
> > 
> > I want to put z1 at the end of third dimension of 'z'. So that the 
> output 
> > will be like 
> > julia> z[*:,:,4*] 
> > 4x3 Array{Int64,2}: 
> >  37  41  45 
> >  38  42  46 
> >  39  43  47 
> >  40  44  48 
> > 
> > I tried push!() and append!(), both gives me errors. 
> > 
> > Kindly suggest a way to do this. I want 'z' to grow in third dimension 
> as 
> > simulation progresses. 
>

Re: [julia-users] randperm run time is slow

[Dan]

A good description of random permutation generation is given in Wikipedia. 
It's the Fisher-Yates algorithm and its variants.

On Saturday, January 23, 2016 at 4:43:36 PM UTC+2, Kevin Squire wrote:
>
> Hi Dan, 
>
> It would also be good if you deleted that post (or edited it and removed 
> the code), for the same reason Viral mentioned: if someone reads the post 
> and then creates a pull request for changing the Julia implementation, it 
> could be argued that that implementation was derived from GPL code, which 
> makes it GPL.  A clean-room implementation (such as creating the patch from 
> a higher level description of the code) is okay.
>
> (Viral, it would probably be good for you to remove the code from your 
> post as well.  I did so in this post below.)
>
> Cheers,
>Kevin
>
> On Saturday, January 23, 2016 at 6:30:07 AM UTC-8, Viral Shah wrote:
>>
>> Unfortunately, this makes the implementation GPL. If you can describe the 
>> algorithm in an issue on github, someone can do a cleanroom implementation.
>>
>> -viral
>>
>> On Saturday, January 23, 2016 at 3:17:19 PM UTC+5:30, Dan wrote:
>>>
>>> Another way to go about it, is to look at R's implementation of a random 
>>> permutation and recreate it in Julia, hoping for similar performance. 
>>> Having done so, the resulting Julia code is:
>>>
>>> 
>>> A size 1,000,000 permutation was generated x2 faster with this method.
>>> But, this method uses uniform floating random numbers, which might not 
>>> be uniform on integers for large enough numbers. In general, it should be 
>>> possible to optimize a more correct implementation. But if R envy is a 
>>> driver, R performance is recovered with a pure Julia implementation (the R 
>>> implementation is in C).
>>>  
>>> On Saturday, January 23, 2016 at 7:26:49 AM UTC+2, Rafael Fourquet wrote:
>>>>
>>>> > Let's capture this as a Julia performance issue on github, 
>>>> > if we can't figure out an easy way to speed this up right away. 
>>>>
>>>> I think I remember having identified a potentially sub-optimal 
>>>> implementation of this function few weeks back (perhaps no more than 
>>>> what Tim suggested) and had planned to investigate further (when time 
>>>> permits...) 
>>>>
>>>

Re: [julia-users] randperm run time is slow

[Dan]

As requested I deleted the post. But algorithms for creating permutations 
are standard and very much in the public domain (what does 
TheArtOfProgramming say?). If someone really invests a little effort into 
it, a more formally correct algorithm can be implemented.

On Saturday, January 23, 2016 at 4:43:36 PM UTC+2, Kevin Squire wrote:
>
> Hi Dan, 
>
> It would also be good if you deleted that post (or edited it and removed 
> the code), for the same reason Viral mentioned: if someone reads the post 
> and then creates a pull request for changing the Julia implementation, it 
> could be argued that that implementation was derived from GPL code, which 
> makes it GPL.  A clean-room implementation (such as creating the patch from 
> a higher level description of the code) is okay.
>
> (Viral, it would probably be good for you to remove the code from your 
> post as well.  I did so in this post below.)
>
> Cheers,
>Kevin
>
> On Saturday, January 23, 2016 at 6:30:07 AM UTC-8, Viral Shah wrote:
>>
>> Unfortunately, this makes the implementation GPL. If you can describe the 
>> algorithm in an issue on github, someone can do a cleanroom implementation.
>>
>> -viral
>>
>> On Saturday, January 23, 2016 at 3:17:19 PM UTC+5:30, Dan wrote:
>>>
>>> Another way to go about it, is to look at R's implementation of a random 
>>> permutation and recreate it in Julia, hoping for similar performance. 
>>> Having done so, the resulting Julia code is:
>>>
>>> 
>>> A size 1,000,000 permutation was generated x2 faster with this method.
>>> But, this method uses uniform floating random numbers, which might not 
>>> be uniform on integers for large enough numbers. In general, it should be 
>>> possible to optimize a more correct implementation. But if R envy is a 
>>> driver, R performance is recovered with a pure Julia implementation (the R 
>>> implementation is in C).
>>>  
>>> On Saturday, January 23, 2016 at 7:26:49 AM UTC+2, Rafael Fourquet wrote:
>>>>
>>>> > Let's capture this as a Julia performance issue on github, 
>>>> > if we can't figure out an easy way to speed this up right away. 
>>>>
>>>> I think I remember having identified a potentially sub-optimal 
>>>> implementation of this function few weeks back (perhaps no more than 
>>>> what Tim suggested) and had planned to investigate further (when time 
>>>> permits...) 
>>>>
>>>

Re: [julia-users] randperm run time is slow

[Dan]

Another way to go about it, is to look at R's implementation of a random 
permutation and recreate it in Julia, hoping for similar performance. 
Having done so, the resulting Julia code is:

function nrandperm(r::AbstractRNG, n::Integer)
   res = Array(typeof(n),n)
   if n == 0
  return res
   end
   a = typeof(n)[i for i=1:n]
   nn = n
   @inbounds for i = 1:n
   j = floor(Int,nn*rand(r))+1
   res[i] = a[j]
   a[j] = a[nn]
   nn -= 1
   end
   return res
end
nrandperm(n::Integer) = nrandperm(Base.Random.GLOBAL_RNG, n)

A size 1,000,000 permutation was generated x2 faster with this method.
But, this method uses uniform floating random numbers, which might not be 
uniform on integers for large enough numbers. In general, it should be 
possible to optimize a more correct implementation. But if R envy is a 
driver, R performance is recovered with a pure Julia implementation (the R 
implementation is in C).
 
On Saturday, January 23, 2016 at 7:26:49 AM UTC+2, Rafael Fourquet wrote:
>
> > Let's capture this as a Julia performance issue on github, 
> > if we can't figure out an easy way to speed this up right away. 
>
> I think I remember having identified a potentially sub-optimal 
> implementation of this function few weeks back (perhaps no more than 
> what Tim suggested) and had planned to investigate further (when time 
> permits...) 
>

[julia-users] Re: What's the difference between unique(a) and union(a)

[Dan]

Essentially, not much. `union` accepts multiple arguments with multiple 
types, but `unique` accepts a single collection.
For example:
union(1,2,3,3) == unique([1,2,3,3])
This would probably make `union` slightly slower.

On Wednesday, January 13, 2016 at 11:58:18 PM UTC+2, Sheehan Olver wrote:
>
> Is there a difference between these two commands:
>
> *julia> **unique([1,2,3,3])*
>
> *3-element Array{Int64,1}:*
>
> * 1*
>
> * 2*
>
> * 3*
>
>
> *julia> **union([1,2,3,3])*
>
> *3-element Array{Int64,1}:*
>
> * 1*
>
> * 2*
>
> * 3*
>
>

Re: [julia-users] Iterating over leafs of a tree: use Cartesian?

[Dan]

For the special problem of getting all the betas in the post, they can be 
calculated with combinatorics: 

betas(levels,parts) = 
  map(x->([x;levels+parts].-[0;x]-1)/parts, combinations(1:(levels+parts-1
),(levels-1)))

The idea is to choose which parts of (0.0,0.2), (0.2,0.4)... (0.8,1.0) go 
into which level. The function does it for any number of levels and parts. 
For example for 4 parts of (0,1) and a 2-level tree:

julia> betas(2,4)

 [0.0,1.0]  
 [0.25,0.75]
 [0.5,0.5]  
 [0.75,0.25]
 [1.0,0.0]  


Needless to say this would be faster and more memory efficient than the 
other solutions.

On Friday, January 8, 2016 at 1:23:13 AM UTC+2, Asbjørn Nilsen Riseth wrote:
>
> Ah, that makes sense. I've mostly been working with Floats, so I didn't 
> realise there was an Int type.
>
> On Thu, 7 Jan 2016 at 23:02 Kristoffer Carlsson  > wrote:
>
>> Just an unsolicited tip here.
>>
>> Try not to "over type" your functions. Your function here would be 
>> cumbersome to use for someone on a 32-bit system for example. Unless you 
>> have a good reason, prefer Int instead of explicit Int32 / Int64.
>>
>>
>> On Thursday, January 7, 2016 at 7:18:31 PM UTC+1, Asbjørn Nilsen Riseth 
>> wrote:
>>>
>>> For future reference, I used the following:
>>>
>>> function traverse(numpoints::Int64, β::Array{Int64}, level::Int64 = 1)
>>> N = length(β)
>>> @assert 0 < level < N
>>>
>>> for i = 0:numpoints-sum(β[1:level-1])
>>> β[level] = i
>>> if level < N-2
>>> traverse(numpoints, β, level+1)
>>> else
>>> for j = 0:numpoints-sum(β[1:N-2])
>>> β[N-1] = j
>>> β[N] = numpoints - sum(β[1:N-1])
>>> @show β/numpoints
>>> end
>>> end
>>> end
>>> end
>>> traverse(5, Array{Int64}(4))
>>>
>>> Asbjørn
>>>
>>>
>>> On Thursday, 7 January 2016 17:14:03 UTC, Asbjørn Nilsen Riseth wrote:

 Thank you both. I'll have a look at a recursion approach.


 On Thursday, 7 January 2016 13:47:02 UTC, Tim Holy wrote:
>
> This is a problem for recursion. See, for example, the LightGraphs or 
> Graphs 
> packages. 
>
> Best, 
> --Tim 
>
> On Thursday, January 07, 2016 02:10:18 AM Asbjørn Nilsen Riseth wrote: 
> > I would like to loop over all the \beta_i values of the following 
> tree. 
> > 
> > <
> https://lh3.googleusercontent.com/-EBjanFs-BC4/Vo409IFAiTI/MEw/U7iq
>  
> > Ev9nPmU/s1600/tree.png> 
> > 
> > What is the best Julia way of doing this? I'm also interested in 
> > generalising this to N levels. 
> > 
> > 
> > For the given tree, the following works. *Is there a way to use 
> @nloops to 
> > shorten this?* 
> > n = Int(1/0.2) 
> > for i_3 = 0:n 
> > for i_2 = 0:n-i_3 
> > for i_1 = 0:n-i_2-i_3 
> > i_0 = n-sum(@ntuple 3 i) 
> > β = [i_3, i_2, i_1, i_0]/n 
> > @show β 
> > end 
> > end 
> > end 
>
>

Re: [julia-users] Re: Optimize sum of elements in an Array according to some conditions

[Dan]

Nice. In this case, you can go faster with binary searching for the start 
and end indices of the sum. This is easily implemented with 
`searchsortedfirst` Julia built-in (haven't coded and benchmarked it, but 
will be faster, especially if vectors are longer).

On Wednesday, December 30, 2015 at 11:54:27 PM UTC+2, Charles Santana wrote:
>
> Thanks, Dan! Indeed my "x" vector is sorted and your suggestion is really 
> fast!
>
> Best,
>
> Charles
>
> On 30 December 2015 at 12:08, Dan > wrote:
>
>> Thanks Kristoffer, turns out there is always interesting stuff in the bag 
>> of optimization tricks.
>> Regarding the original function, a cheat could make it faster: The `x` 
>> vector is sorted, which means:
>>
>> function calcSum4(x::Array{Float64,1}, y::Array{Float64,1}, Ei::Float64, 
>> Ef::Float64, N::Int64)
>>mysum=0.0::Float64;
>>i=1
>>@inbounds while i<=N && x[i]<=Ei i+=1 ; end
>>j=i
>>@inbounds while j<=N && x[j]<=Ef j+=1 ; end
>>@inbounds @simd for k=i:(j-1) mysum += y[k] ; end
>>return(mysum);
>>end 
>>
>> returns the same answer.
>> there are always more options ;)
>>
>> On Wednesday, December 30, 2015 at 12:19:41 PM UTC+2, Charles Santana 
>> wrote:
>>>
>>> The magic of @inbounds and @simd :)
>>>
>>> Thanks, Kristoffer!
>>>
>>> Charles
>>>
>>>
>>> On Wednesday, December 30, 2015, Kristoffer Carlsson  
>>> wrote:
>>>
>>>> If you want to get an even faster version you could do something like:
>>>>
>>>> function calcSum_simd{T}(x::Vector{T}, y::Vector{T}, Ei::T, Ef::T)
>>>> mysum = zero(T)
>>>> @inbounds @simd for i in eachindex(x, y)
>>>>  mysum += ifelse(Ei < x[i] <= Ef, y[i], zero(T)) 
>>>>  
>>>> end
>>>> return mysum
>>>> end
>>>>
>>>> which would use SIMD instructions.
>>>>
>>>> Timing difference:
>>>>
>>>> N = 1000
>>>> y = rand(N);
>>>> x = rand(N)
>>>> Ei = 0.2;
>>>> Ef = 0.7;
>>>>
>>>> julia> @time calcSum_simd(x,y,Ei, Ef);
>>>>   0.021155 seconds (5 allocations: 176 bytes)
>>>>
>>>>
>>>> julia> @time calcSum(x,y,Ei, Ef)
>>>>   0.069911 seconds (5 allocations: 176 bytes)
>>>>
>>>>
>>>> Regarding map being slow. That is worked on here 
>>>> https://github.com/JuliaLang/julia/pull/13412
>>>>
>>>>
>>>> On Wednesday, December 30, 2015 at 3:05:47 AM UTC+1, Charles Santana 
>>>> wrote:
>>>>>
>>>>> Sorry, there was a typo in the function calcSum2. Please consider the 
>>>>> following code:
>>>>>
>>>>> function calcSum2(x::Array{Float64,1}, y::Array{Float64,1}, 
>>>>> Ei::Float64, Ef::Float64, N::Int64)
>>>>>
>>>>> return sum(y[map(v -> Ei < v <= Ef, x)]);
>>>>> end
>>>>>   
>>>>>
>>>>> And so the results of the calls for this function change a bit (but 
>>>>> not the performance):
>>>>>
>>>>> @time calcSum2(x,y,Ei,Ef,N)
>>>>>   0.000110 seconds (1.01 k allocations: 20.969 KB)
>>>>> 246.1975746121703
>>>>>
>>>>> @time calcSum2(x,y,Ei,Ef,N)
>>>>>   0.79 seconds (1.01 k allocations: 20.969 KB)
>>>>> 246.1975746121703
>>>>>
>>>>> @time calcSum2(x,y,Ei,Ef,N)
>>>>>   0.51 seconds (1.01 k allocations: 20.969 KB)
>>>>> 246.1975746121703
>>>>>
>>>>>
>>>>> Thanks again, sorry for this inconvenience!
>>>>>
>>>>> Charles
>>>>>
>>>>> On 30 December 2015 at 03:00, Charles Novaes de Santana <
>>>>> charles...@gmail.com> wrote:
>>>>>
>>>>>> Dear all,
>>>>>>
>>>>>> In a project I am developing a @profile shows me that the slowest 
>>>>>> part of the code is the sum of elements of an Array that follow some 
>>>>>> conditions. 
>>>>>>
>>>>>>

[julia-users] Re: Optimize sum of elements in an Array according to some conditions

[Dan]

Thanks Kristoffer, turns out there is always interesting stuff in the bag 
of optimization tricks.
Regarding the original function, a cheat could make it faster: The `x` 
vector is sorted, which means:

function calcSum4(x::Array{Float64,1}, y::Array{Float64,1}, Ei::Float64, 
Ef::Float64, N::Int64)
   mysum=0.0::Float64;
   i=1
   @inbounds while i<=N && x[i]<=Ei i+=1 ; end
   j=i
   @inbounds while j<=N && x[j]<=Ef j+=1 ; end
   @inbounds @simd for k=i:(j-1) mysum += y[k] ; end
   return(mysum);
   end 

returns the same answer.
there are always more options ;)

On Wednesday, December 30, 2015 at 12:19:41 PM UTC+2, Charles Santana wrote:
>
> The magic of @inbounds and @simd :)
>
> Thanks, Kristoffer!
>
> Charles
>
> On Wednesday, December 30, 2015, Kristoffer Carlsson  > wrote:
>
>> If you want to get an even faster version you could do something like:
>>
>> function calcSum_simd{T}(x::Vector{T}, y::Vector{T}, Ei::T, Ef::T)
>> mysum = zero(T)
>> @inbounds @simd for i in eachindex(x, y)
>>  mysum += ifelse(Ei < x[i] <= Ef, y[i], zero(T)) 
>>  
>> end
>> return mysum
>> end
>>
>> which would use SIMD instructions.
>>
>> Timing difference:
>>
>> N = 1000
>> y = rand(N);
>> x = rand(N)
>> Ei = 0.2;
>> Ef = 0.7;
>>
>> julia> @time calcSum_simd(x,y,Ei, Ef);
>>   0.021155 seconds (5 allocations: 176 bytes)
>>
>>
>> julia> @time calcSum(x,y,Ei, Ef)
>>   0.069911 seconds (5 allocations: 176 bytes)
>>
>>
>> Regarding map being slow. That is worked on here 
>> https://github.com/JuliaLang/julia/pull/13412
>>
>>
>> On Wednesday, December 30, 2015 at 3:05:47 AM UTC+1, Charles Santana 
>> wrote:
>>>
>>> Sorry, there was a typo in the function calcSum2. Please consider the 
>>> following code:
>>>
>>> function calcSum2(x::Array{Float64,1}, y::Array{Float64,1}, Ei::Float64, 
>>> Ef::Float64, N::Int64)
>>>
>>> return sum(y[map(v -> Ei < v <= Ef, x)]);
>>> end
>>>   
>>>
>>> And so the results of the calls for this function change a bit (but not 
>>> the performance):
>>>
>>> @time calcSum2(x,y,Ei,Ef,N)
>>>   0.000110 seconds (1.01 k allocations: 20.969 KB)
>>> 246.1975746121703
>>>
>>> @time calcSum2(x,y,Ei,Ef,N)
>>>   0.79 seconds (1.01 k allocations: 20.969 KB)
>>> 246.1975746121703
>>>
>>> @time calcSum2(x,y,Ei,Ef,N)
>>>   0.51 seconds (1.01 k allocations: 20.969 KB)
>>> 246.1975746121703
>>>
>>>
>>> Thanks again, sorry for this inconvenience!
>>>
>>> Charles
>>>
>>> On 30 December 2015 at 03:00, Charles Novaes de Santana <
>>> charles...@gmail.com> wrote:
>>>
 Dear all,

 In a project I am developing a @profile shows me that the slowest part 
 of the code is the sum of elements of an Array that follow some 
 conditions. 

 Please consider the following code:

 y = rand(1000);
 x = collect(0.0:0.001:0.999);
 Ei = 0.2;
 Ef = 0.7;
 N = length(x)

 I want to calculate the sum of elements in "y" for which elements the 
 respective values in "x" are between "Ei" and "Ef". If I was using R, for 
 example, I would use something like:

 mysum = sum(y[which((x < Ef)&&(x > Ei))]); #(not tested in R, but I 
 suppose that is the way to do it)

 In Julia, I can think in at least two ways to calculate it:

 function calcSum(x::Array{Float64,1}, y::Array{Float64,1}, Ei::Float64, 
 Ef::Float64, N::Int64)
 mysum=0.0::Float64;
 for(i in 1:N)
  if( Ei < x[i] <= Ef)
  mysum += y[i];
  end
 end
 return(mysum);
 end

 function calcSum2(x::Array{Float64,1}, y::Array{Float64,1}, 
 Ei::Float64, Ef::Float64, N::Int64)
 return sum(y[map(v -> Ei < v < Ef, x)]);
 end

 As you can see below, for the first function (calcSum) I got a much 
 better performance than for the second one (minimum 10x faster). 


  @time calcSum(x,y,Ei,Ef,N)
   0.003986 seconds (2.56 k allocations: 125.168 KB)
 246.19757461217014

 @time calcSum(x,y,Ei,Ef,N)
   0.03 seconds (5 allocations: 176 bytes)
 246.19757461217014

 @time calcSum(x,y,Ei,Ef,N)
   0.02 seconds (5 allocations: 176 bytes)
 246.19757461217014

 @time calcSum2(x,y,Ei,Ef,N)
   0.003762 seconds (1.61 k allocations: 53.743 KB)
 245.48156534879303

 @time calcSum2(x,y,Ei,Ef,N)
   0.50 seconds (1.01 k allocations: 20.969 KB)
 245.48156534879303

 @time calcSum2(x,y,Ei,Ef,N)
   0.000183 seconds (1.01 k allocations: 20.969 KB)
 245.48156534879303

 Does any one have an idea about how to impr

[julia-users] Re: python a[range(x), y] equivalent in Julia

[Dan]

Is the following red expression a good Julia equivalent:
 

julia> a = [ -1 2; 3 -4 ]

 ⋮ 

julia> y = [ 1, 2 ]

 ⋮

julia> [a[i,y[i]] for i=1:2]
2-element Array{Any,1}:
 -1
 -4

(it's also fast)

On Sunday, December 20, 2015 at 2:44:20 AM UTC+2, Lex wrote:
>
> I am expecting [-1 -4] in the Julia example. 
>
> On Saturday, December 19, 2015 at 4:42:07 PM UTC-8, Lex wrote:
>>
>> Hi
>>
>> In Python, I am able to select values like this:
>>
>> >>> a = np.matrix([[-1,2],[3,-4]])
>> >>> y = np.matrix([1,0])
>> >>> a[range(2),y]
>> matrix([[2, 3]])
>> >>>
>>
>> Is there any equivalent in Julia or using a loop is the only way?
>>
>> julia> a
>> 2x2 Array{Int64,2}:
>>  -1   2
>>   3  -4
>>
>> julia> y
>> 1x2 Array{Int64,2}:
>>  1  2
>>
>> julia> a[:, y]
>> ERROR: MethodError: `index_shape_dim` has no method matching 
>> index_shape_dim(::Array{Int64,2}, ::Int64, ::Array{Int64,2})
>>
>> You might have used a 2d row vector where a 1d column vector was required.
>> Note the difference between 1d column vector [1,2,3] and 2d row vector [1 
>> 2 3].
>> You can convert to a column vector with the vec() function.
>> Closest candidates are:
>>   index_shape_dim(::Any, ::Any, ::Real...)
>>   index_shape_dim(::Any, ::Any, ::Colon)
>>   index_shape_dim(::Any, ::Any, ::Colon, ::Any, ::Any...)
>>   ...
>>  in getindex at abstractarray.jl:488
>>
>> julia> a[:, y[]]
>> 2-element Array{Int64,1}:
>>  -1
>>   3
>>
>> julia> a[:, y[:]]
>> 2x2 Array{Int64,2}:
>>  -1   2
>>   3  -4
>>
>> julia>
>>
>>
>>

[julia-users] Re: Using callable types or FastAnonymous with Sundials

[Dan]

The parameter for the `cvode` function is `f` and so is the function you 
want to use. These get confused, and it tries to use the "function" `J` 
instead. Changing the parameter name to something other than `{f}` should 
work

On Wednesday, December 16, 2015 at 4:26:41 PM UTC+2, Simon Frost wrote:
>
> Dear Julia Users,
>
> I'm trying to speed up some code that employs passing functions as 
> arguments. One part of the code solves an ODE; if I use CVODE from 
> Sundials, and rewrite the function to accept callable types, I get a 
> ReadOnlyMemoryError - as I don't know the Sundials API, can someone help me 
> with where I'm going wrong? Code below.
>
> Best
> Simon
>
> **
>
> using Sundials
>
> function cvode{f}(::Type{f}, y0::Vector{Float64}, t::Vector{Float64}; 
> reltol::Float64=1e-4, abstol::Float64=1e-6)
> neq = length(y0)
> mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
> flag = Sundials.CVodeInit(mem, cfunction(Sundials.cvodefun, Int32, 
> (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})), 
> t[1], Sundials.nvector(y0))
> flag = Sundials.CVodeSetUserData(mem, f)
> flag = Sundials.CVodeSStolerances(mem, reltol, abstol)
> flag = Sundials.CVDense(mem, neq)
> yres = zeros(length(t), length(y0))
> yres[1,:] = y0
> y = copy(y0)
> tout = [0.0]
> for k in 2:length(t)
> flag = Sundials.CVode(mem, t[k], y, tout, Sundials.CV_NORMAL)
> yres[k,:] = y
> end
> Sundials.CVodeFree([mem])
> return yres
> end
>
> function f(t, y, ydot)
> ydot[1] = 0.1*(-72-y[1])+0.1*1.4*exp((y[1]+48)/1.4)+10
> ydot[3] = 0.
> ydot[2] = 0.
> end
>
> immutable J; end
> call(::Type{J},t, y, ydot) = f(t, y, ydot)
>
> t = [0.1:0.0001:1]
> res = cvode(J, [-60.0, 0.0, 0.0], t);
>

[julia-users] Re: Compute selected eigenvectors with eigvecs

[Dan]

Linear algebra tells us Symmetric matrices have real eigenvalues and 
general matrices have complex eigenvalues in general. `eigvecs` can be 
called with an upper and lower bound on the eigenvalues of the eigenvectors 
to be found, but upper and lower bounds are applicable only to Real 
numbers, which means we need to use the Symmetric matrix version of 
`eigvals` and `eigvecs` and then things work out better.

In code:
julia> A = rand(3,3)
3x3 Array{Float64,2}:
 0.636680.252549  0.467667 
 0.0372453  0.288091  0.0918441
 0.383552   0.478298  0.148384 

julia> B = Symmetric(A+A')
3x3 Symmetric{Float64,Array{Float64,2}}:
 1.27336   0.289794  0.851219
 0.289794  0.576181  0.570142
 0.851219  0.570142  0.296767

julia> eigvals(B)
3-element Array{Float64,1}:
 -0.349693
  0.516073
  1.97993 

julia> eigvecs(B,1.9,2.0)
3x1 Array{Float64,2}:
 -0.772992
 -0.369141
 -0.515963


Because of numerical errors, when looking for a specific eigenvalue it 
might be prudent to use enlarge the range by some small factor. This info 
is the outcome of using `help` and `methods` in the REPL, so more 
information might be available with a little more investigation.

On Wednesday, December 16, 2015 at 11:43:42 AM UTC+2, Robert DJ wrote:
>
> I would like to compute the eigenvector of a particular eigenvalue and 
> according to the docs, eigvecs is up for the job.
> However, I get errors. The following call should give the same output as 
> eigvecs(A):
>
> julia> A = rand(3,3)
> 3x3 Array{Float64,2}:
>  0.165111  0.995671  0.393572
>  0.496099  0.640361  0.595876
>  0.632881  0.39090.723987
>
> julia> v = eigvals(A)
> 3-element Array{Float64,1}:
>   1.68982
>  -0.230301
>   0.0699443
>
> julia> eigvecs(A, v)
> ERROR: MethodError: `eigfact` has no method matching 
> eigfact(::Array{Float64,2}, ::Array{Float64,1})
> Closest candidates are:
>   
> eigfact{T}(::Union{DenseArray{T,2},SubArray{T,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64,LD}})
>   eigfact{TA,TB}(::AbstractArray{TA,2}, ::AbstractArray{TB,2})
>  in eigvecs at linalg/eigen.jl:70
>
> julia> eigvecs(A, v)
>
> I have
>
> julia> versioninfo()
> Julia Version 0.4.0
> Commit 0ff703b (2015-10-08 06:20 UTC)
> Platform Info:
>   System: Darwin (x86_64-apple-darwin14.5.0)
>   CPU: Intel(R) Core(TM)2 Duo CPU P7350  @ 2.00GHz
>   WORD_SIZE: 64
>   BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Penryn)
>   LAPACK: libopenblas64_
>   LIBM: libopenlibm
>   LLVM: libLLVM-3.3
>
>
>

[julia-users] Re: DataFrame melt question

[Dan]

the following re-reads the header and generates a dictionary which assigns 
the original column name to the converted one, in a one-liner-ish:

df = readtable("/the/file.csv")

h = Dict(zip(keys(df.colindex.lookup),split(open("/tmp/file.csv") do f 
chomp(readline(f)) ; end,",")[collect(values(df.colindex.lookup))]))


now aside from using `h` in other ways, you can do:

melteddf[:Region] = [h[r] for r in melteddf[:Region]]


to fix the `melteddf`.

On Wednesday, December 16, 2015 at 2:39:57 AM UTC+2, David Anthoff wrote:
>
> Hi,
>
>  
>
> I have a csv file that roughly looks like this:
>
>  
>
>  
>
> Year,Name of country 1, Name of country 2
>
> 1950, 5., 6.
>
> 1951, 6., 8.
>
>  
>
> The real file has more columns and rows.
>
>  
>
> I want to bring this into tidy format, so that I have a DataFrame that 
> looks like this:
>
>  
>
> Year, Region, Value
>
> 1950, Name of country 1, 5.
>
> 1950, Name of country 2, 6.
>
> 1951, Name of country 1, 6.
>
> 1951, Name of country 2, 8.
>
>  
>
> Right now I read the file with readtable into a DataFrame and then use
>
>  
>
> melt(df, :Year)
>
>  
>
> This gives me the right structure, but now all the country names are 
> messed up, e.g. they look like “Name_of_country_1” instead of “Name of 
> country 1”.
>
>  
>
> I understand why that is the case, i.e. readtable converts strings into 
> symbols and has to insert these underscores, but I’m wondering whether the 
> original string value is preserved somewhere, and could be used in the melt 
> operation in some way?
>
>  
>
> Thanks,
>
> David
>
>  
>
> --
>
> David Anthoff
>
> University of California, Berkeley
>
>  
>
> http://www.david-anthoff.com
>
>  
>

Re: [julia-users] Re: A non-breaking Julia control-flow change to separate the One language from mortal languages, doomed to die

[Dan]

Good idea regarding macro.

Not really trying to convince people, just to entice them and gauge the 
response. The traditional method of achieving the same results is OK, but 
I'm just thinking this could be better. Having said all of this, stability 
of a language has much merit as well.

On Wednesday, December 9, 2015 at 6:27:00 PM UTC+2, Erik Schnetter wrote:
>
> One way to convince people that a new feature is worthwhile is to 
> implement an approximation of the feature, e.g. via a macro. If it then 
> sees widespread use, you have a strong argument in your favour.
>
> Here's a possibility, using the customary "macro-in-front, 
> begin-end-around-everything" approach:
>
> @loop begin for i in 1:10
> ...
> if i==8 break end
> end
> finally
> info("did not break")
> end
>
> I'm sure the syntax needs cleaning up.
>
> -erik
>
>
>
> On Wed, Dec 9, 2015 at 10:47 AM, Seth  > wrote:
>
>> Seems to me that this is just a more complex repeat...until or do...while 
>> bottom-conditional loop, which has already been discussed and rejected (to 
>> my disappointment).
>>
>>
>> On Wednesday, December 9, 2015 at 7:05:43 AM UTC-8, Dan wrote:
>>>
>>> *A non-breaking Julia control-flow change to separate the One language from 
>>> *
>>>
>>> *mortal languages, doomed to die
>>> *
>>> From a design patterns and DRY perspective it seems `for` loops together 
>>> with the `break` statement 
>>>
>>> are serving as a useful `try`/`catch` mechanism for non-error situations. A 
>>> common pattern defines 
>>>
>>> a variable before a `for` loop and tries to discern, using said variable, 
>>> the circumstances of exiting 
>>>
>>> the `for` loop. Leveraging the embedded mental construct of `try`/`catch` 
>>> blocks, we can simplify 
>>>
>>> this repeating pattern and make it easier and more readable.
>>>
>>> The syntax can look as follows:
>>>
>>> for ...(loop vars here)... 
>>>
>>> ...(local vars here)...
>>>
>>> ...
>>>
>>> if ...
>>>
>>> break 
>>>
>>> end
>>>
>>> ...
>>>
>>> catch 
>>>
>>> ...(access local/loop vars and  here)...
>>>
>>> 
>>>
>>> finally
>>>
>>> ...(access local/loop vars here)...
>>>
>>> 
>>>
>>> end 
>>>
>>> Syntax changes worthy of consideration should be self-explanatory to 
>>> many, as I hope this suggestion is.
>>> But an example should help:
>>>
>>> After change:
>>>
>>> for line in eachline(f) 
>>>
>>> fields = split(line) 
>>>
>>> status = fileds[2] 
>>>
>>> if status=="FAIL" 
>>>
>>> break 
>>>
>>> end 
>>>
>>> catch 
>>>
>>> println("Test $(fields[1]) failed with reason $(fields[3])") 
>>>
>>> finally 
>>>
>>> println("Test successful!") 
>>>
>>> end 
>>>
>>> Before change:
>>>
>>>  
>>>
>>> fields = [] 
>>>
>>> didfail = false 
>>>
>>> for line in eachline(f) 
>>>
>>> fields = split(line) 
>>>
>>> status = fields[2] 
>>>
>>> if status=="FAIL" 
>>>
>>> didfail = true 
>>>
>>> break 
>>>
>>> end 
>>>
>>> end 
>>>
>>> if didfail 
>>>
>>> println("Test $(fields[1]) failed with reason $(fields[3])") 
>>>
>>> else 
>>>
>>> println("Test successful!") 
>>>
>>> end 
>>>
>>> The above example does not use the `break ` suggestion. 
>>> With multiple `break` circumstances it should
>>>
>>> be useful to have separate code for each circumstance. One option is to 
>>> replicate the `try`/`catch` form 
>>>
>>> with a `` passed as ``, another would be to use `catch 
>>> ` to run catch block 
>>>
>>> if "`break  == catch `" (the `` can be a human 
>>> readable symbol). A direct jump in the 
>>>
>>> LLVM code to the appropriate `catch` block would be an easy optimization to 
>>> implement.
>>>
>>> Another issue to consider is the return value of a `for` block. The 
>>> suggestion above adds the option to change 
>>>
>>> the current default `nothing` return with `` in the new optional 
>>> blocks.
>>>
>>> Notice that the `catch` and `finally` blocks are completely optional and 
>>> therefore backward compatibility 
>>>
>>> is easy. Additionally there are no additional keywords and block nesting 
>>> would not conflict with `try`/`catch`.
>>>
>>> Hey, Julia is still < 0.9 so there is still time for some change.
>>>
>>> Obligatory Disclaimerism:
>>>
>>> 1. This has already been considered before by **name**/**issue**.
>>> 2. Before suggesting I should have studied past discussions more.
>>> 3. Other langauges actually already have this... bla bla.
>>> 4. This post is what it is.
>>>
>>>
>>> This suggestion is also readable in github markdown in the gist: 
>>> https://gist.github.com/getzdan/75136f130c6f8ffeb60e
>>>
>>>
>>> Feedback/Pointers/Gotchas welcome, of course.
>>>
>>>
>
>
> -- 
> Erik Schnetter > 
> http://www.perimeterinstitute.ca/personal/eschnetter/
>

[julia-users] A non-breaking Julia control-flow change to separate the One language from mortal languages, doomed to die

[Dan]

*A non-breaking Julia control-flow change to separate the One language from *

*mortal languages, doomed to die
*
>From a design patterns and DRY perspective it seems `for` loops together with 
>the `break` statement 

are serving as a useful `try`/`catch` mechanism for non-error situations. A 
common pattern defines 

a variable before a `for` loop and tries to discern, using said variable, the 
circumstances of exiting 

the `for` loop. Leveraging the embedded mental construct of `try`/`catch` 
blocks, we can simplify 

this repeating pattern and make it easier and more readable.

The syntax can look as follows:

for ...(loop vars here)... 

...(local vars here)...

...

if ...

break 

end

...

catch 

...(access local/loop vars and  here)...



finally

...(access local/loop vars here)...



end 

Syntax changes worthy of consideration should be self-explanatory to 
many, as I hope this suggestion is.
But an example should help:

After change:

for line in eachline(f) 

fields = split(line) 

status = fileds[2] 

if status=="FAIL" 

break 

end 

catch 

println("Test $(fields[1]) failed with reason $(fields[3])") 

finally 

println("Test successful!") 

end 

Before change:

 

fields = [] 

didfail = false 

for line in eachline(f) 

fields = split(line) 

status = fields[2] 

if status=="FAIL" 

didfail = true 

break 

end 

end 

if didfail 

println("Test $(fields[1]) failed with reason $(fields[3])") 

else 

println("Test successful!") 

end 

The above example does not use the `break ` suggestion. 
With multiple `break` circumstances it should

be useful to have separate code for each circumstance. One option is to 
replicate the `try`/`catch` form 

with a `` passed as ``, another would be to use `catch ` 
to run catch block 

if "`break  == catch `" (the `` can be a human readable 
symbol). A direct jump in the 

LLVM code to the appropriate `catch` block would be an easy optimization to 
implement.

Another issue to consider is the return value of a `for` block. The suggestion 
above adds the option to change 

the current default `nothing` return with `` in the new optional 
blocks.

Notice that the `catch` and `finally` blocks are completely optional and 
therefore backward compatibility 

is easy. Additionally there are no additional keywords and block nesting would 
not conflict with `try`/`catch`.

Hey, Julia is still < 0.9 so there is still time for some change.

Obligatory Disclaimerism:

1. This has already been considered before by **name**/**issue**.
2. Before suggesting I should have studied past discussions more.
3. Other langauges actually already have this... bla bla.
4. This post is what it is.


This suggestion is also readable in github markdown in the gist: 
https://gist.github.com/getzdan/75136f130c6f8ffeb60e


Feedback/Pointers/Gotchas welcome, of course.

[julia-users] A non-breaking Julia control-flow change to separate the One language from mortal languages, doomed to die

[Dan]

*A non-breaking Julia control-flow change to separate the One language from *

*mortal languages, doomed to die
*
>From a design patterns and DRY perspective it seems `for` loops together with 
>the `break` statement 

are serving as a useful `try`/`catch` mechanism for non-error situations. A 
common pattern defines 

a variable before a `for` loop and tries to discern, using said variable, the 
circumstances of exiting 

the `for` loop. Leveraging the embedded mental construct of `try`/`catch` 
blocks, we can simplify 

this repeating pattern and make it easier and more readable.

The syntax can look as follows:


for ...(loop vars here)... 

...(local vars here)...

...

if ...

break 

end

...

catch 

...(access local/loop vars and  here)...



finally

...(access local/loop vars here)...



end 


Syntax changes worthy of consideration should be self-explanatory to many, as I 
hope this suggestion is.
But an example should help:

After change:


for line in eachline(f) 

fields = split(line) 

status = fileds[2] 

if status=="FAIL" 

break 

end 

catch 

println("Test $(fields[1]) failed with reason $(fields[3])") 

finally 

println("Test successful!") 

end 

Before change:
 

fields = [] 

didfail = false 

for line in eachline(f) 

fields = split(line) 

status = fields[2] 

if status=="FAIL" 

didfail = true 

break 

end 

end 

if didfail 

println("Test $(fields[1]) failed with reason $(fields[3])") 

else 

println("Test successful!") 

end 


The above example does not use the `break ` suggestion. With multiple 
`break` circumstances it should

be useful to have separate code for each circumstance. One option is to 
replicate the `try`/`catch` form 

with a `` passed as ``, another would be to use `catch ` 
to run catch block 

if "`break  == catch `" (the `` can be a human readable 
symbol). A direct jump in the 

LLVM code to the appropriate `catch` block would be an easy optimization to 
implement.

Another issue to consider is the return value of a `for` block. The suggestion 
above adds the option to change 

the current default `nothing` return with `` in the new optional 
blocks.

Notice that the `catch` and `finally` blocks are completely optional and 
therefore backward compatibility 

is easy. Additionally there are no additional keywords and block nesting would 
not conflict with `try`/`catch`.

Hey, Julia is still < 0.9 so there is still time for some change.

Obligatory Disclaimerism:

1. This has already been considered before by **name**/**issue**.
2. Before suggesting I should have studied past discussions more.
3. Other langauges actually already have this... bla bla.
4. This post is what it is.


This suggestion is also readable in github markdown in the gist: 
https://gist.github.com/getzdan/75136f130c6f8ffeb60e


Feedback/Pointers/Gotchas welcome, of course.

[julia-users] A non-breaking Julia control-flow change to separate the One language from mortal languages, doomed to die

[Dan]

*A non-breaking Julia control-flow change to separate the One language from 
mortal languages, doomed to die
*
>From a design patterns and DRY perspective it seems `for` loops together with 
>the `break` statement are serving
as a useful `try`/`catch` mechanism for non-error situations. A common pattern 
defines a variable before a `for` loop and tries to discern, using said 
variable, the circumstances of exiting the `for` loop. Leveraging the embedded 
mental construct of `try`/`catch` blocks, we can simplify this repeating 
pattern and make it easier and more readable.

The syntax can look as follows:


for ...(loop vars here)... 

...(local vars here)...

...

if ...

break 

end

...

catch 

...(access local/loop vars and  here)...



finally

...(access local/loop vars here)...



end 


Syntax changes worthy of consideration should be self-explanatory to many, as I 
hope this suggestion is.
But an example should help:

After change:


for line in eachline(f) 

fields = split(line) 

status = fileds[2] 

if status=="FAIL" 

break 

end 

catch 

println("Test $(fields[1]) failed with reason $(fields[3])") 

finally 

println("Test successful!") 

end 

Before change:
 

fields = [] 

didfail = false 

for line in eachline(f) 

fields = split(line) 

status = fields[2] 

if status=="FAIL" 

didfail = true 

break 

end 

end 

if didfail 

println("Test $(fields[1]) failed with reason $(fields[3])") 

else 

println("Test successful!") 

end 


The above example does not use the `break ` suggestion. With multiple 
`break` circumstances it should be useful to have separate code for each 
circumstance. One option is to replicate the `try`/`catch` form with a 
`` passed as ``, another would be to use `catch ` to run 
catch block if "`break  == catch `" (the `` can be a human 
readable symbol). A direct jump in the LLVM code to the appropriate `catch` 
block would be an easy optimization to implement.

Another issue to consider is the return value of a `for` block. The suggestion 
above adds the option to change the current default `nothing` return with 
`` in the new optional blocks.

Notice that the `catch` and `finally` blocks are completely optional and 
therefore backward compatibility is easy. Additionally there are no additional 
keywords and block nesting would not conflict with `try`/`catch`.

Hey, Julia is still < 0.9 so there is still time for some change.

Obligatory Disclaimerism:

1. This has already been considered before by **name**/**issue**.
2. Before suggesting I should have studied past discussions more.
3. Other langauges actually already have this... bla bla.
4. This post is what it is.


This suggestion is also readable in github markdown in the gist: 
https://gist.github.com/getzdan/75136f130c6f8ffeb60e


Feedback/Pointers/Gotchas welcome, of course.

[julia-users] Re: hex to bit notation

[Dan]

The shift and AND operators are usually used in these situations (>> and 
&). With the following definition:

getbit(x::UInt64,i) = Int((x >> (64-i))&1)

You can get any bit from a 64-bit unsigned integer. For example 
getbit(0xa7f1d92a82c8d8fe,29).

On Wednesday, December 2, 2015 at 11:41:51 AM UTC+2, Martin Somers wrote:
>
> Just wondering 
>
>  
>
> I have a 32 bit hex number as following -  0xa7f6d92a82c8d8f2
>
>  
>
> How do I acccertain the 29th bit 
>
>  
>
> Functions bin and bits all return an ascii string 
>
>  
>
>  
>
> bin(0xa7f1d92a82c8d8fe)
>
> "1010011100011101100100101010101011001000110110001110"
>
>  
>
> what is the most efficient way of finding a bit of a hex value 
>
>  
>
> bin(0xa7f1d92a82c8d8fe)[29] and then convert this string to a bool?
>
>  
>
> Interested to see the alternatives
>
>
> cheers
>
> M
>
>  
>
>  
>

Re: [julia-users] Re: `findmin` with arbitrary comparison operator?

[Dan]

suggestion: call the named argument `lt` to match the argument of `sort`.

On Tuesday, December 1, 2015 at 6:23:07 PM UTC+2, Júlio Hoffimann wrote:
>
> Hi Erik,
>
> What I meant was a derived type from Julia's Tuple type. Though I agree 
> with you that this may be too much of an overhead for such a simple task.
>
> I'll open an issue on GitHub to ask if others are ok with an additional 
> predicate option a la C++:
>
> findmin(array, pred=<)
>
> If they are all positive, I'll add it for you.
>
> Best,
> -Júlio
>

Re: [julia-users] Re: slow wordcount

[Dan]

and replacing
println("$(t.first)\t$(t.second)")

with
@printf("%s\t%d\n",t.first,t.second)

also halves the print time (which might or might not make a big difference 
but definitely > 1 second)

On Monday, November 30, 2015 at 6:20:54 PM UTC+2, Dan wrote:
>
> Using the `customlt` function for the sort order cut the time in half on a 
> test file. So try:
>
> customlt(a,b) = (b.second < a.second) ? true : b.second == a.second ? 
> a.first < b.first : false
>
> function main()
> wc = Dict{UTF8String,Int64}()
> for l in eachline(STDIN)
> for w in split(l)
> wc[w]=get(wc, w, 0) + 1
> end
> end
>
> v = collect(wc) 
> sort!(v,lt=customlt) # in-place sort saves a memory copy
> for t in v
> println("$(t.first)\t$(t.second)")
> end
> end
>
> main()
>
>
>
> On Monday, November 30, 2015 at 5:31:20 PM UTC+2, Attila Zséder wrote:
>>
>> Hi,
>>
>>
>> The data I'm using is part of (Hungarian) Wikipedia dump with 5M lines of 
>> text. On this data, python runs for 65 seconds, cpp for 35 seconds, julia 
>> baseline for 340 seconds, julia with FastAnonymous.jl for 280 seconds. (See 
>> https://github.com/juditacs/wordcount#leaderboard for details)
>>
>> Dan:
>> I can use external packages, it's not a big issue. However, FastAnonymous 
>> didn't give results comparable to python.
>> The baseline python code I compare to is here:
>> https://github.com/juditacs/wordcount/blob/master/python/wordcount_py2.py
>>
>> 2) The community is part of the language, so it should be regarded when 
>> making considerations.
>> What do you mean by this?
>> My (our) purpose is not to evaluate this language or that one is 
>> better/faster/etc because it is faster in unoptimized word counting. So I 
>> don't want to make any judgements, considerations and anything like this. 
>> This is just for fun. And even though it looks like _my_ julia 
>> implementation of wc is not fast right now, I didn't lose interest in 
>> following what's going on with this language.
>>
>>
>> Your other points:
>> 1) I do this with all the other languages as well. The test runs for 
>> about 30-300 seconds. If Julia load time or any other thing takes serious 
>> amount of time, then it does. This test is not precise, I didn't include 
>> c++ compile time for example, but it took less than a second. But I felt 
>> like my implementation is dummy, and other things take my time, not Julia 
>> load.
>> 2) What if my test is about IO + dictionary storage? Then I have to 
>> include printouts into my test.
>> 3) I think 5m lines of text file is enough to avoid this noises.
>>
>>
>>
>> Tim:
>> Yes, I did this code split, and with larger files it looked like after 
>> sorting, dictionary manipulation (including hashes) took most of the time, 
>> and printing was less of an issue. But I do have to analyze this more 
>> precisely, seeing your numbers.
>>
>>
>> Thank you all for your help!
>>
>> Attila
>>
>> On Mon, Nov 30, 2015 at 4:20 PM, Tim Holy  wrote:
>>
>>> If you don't want to figure out how to use the profiler, your next best 
>>> bet is
>>> to split out the pieces so you can understand where the bottleneck is. 
>>> For
>>> example:
>>>
>>> function docount(io)
>>> wc = Dict{AbstractString,Int64}()
>>> for l in eachline(io)
>>> for w in split(l)
>>> wc[w]=get(wc, w, 0) + 1
>>> end
>>> end
>>> wc
>>> end
>>>
>>> @time open("somefile.tex") do io
>>>docount(io)
>>>end;
>>>   0.010617 seconds (27.70 k allocations: 1.459 MB)
>>>
>>> vs
>>>
>>> @time open("somefile.tex") do io
>>>main(io)
>>>end;
>>> # < lots of printed output >
>>>   1.233154 seconds (330.59 k allocations: 10.829 MB, 1.53% gc time)
>>>
>>> (I modified your `main` to take an io input.)
>>>
>>> So it's the sorting and printing which is taking 99% of the time. Most 
>>> of that
>>> turns out to be the printing.
>>>
>>> --Tim
>>>
>>>
>>

Re: [julia-users] Re: slow wordcount

[Dan]

Using the `customlt` function for the sort order cut the time in half on a 
test file. So try:

customlt(a,b) = (b.second < a.second) ? true : b.second == a.second ? 
a.first < b.first : false

function main()
wc = Dict{UTF8String,Int64}()
for l in eachline(STDIN)
for w in split(l)
wc[w]=get(wc, w, 0) + 1
end
end

v = collect(wc) 
sort!(v,lt=customlt) # in-place sort saves a memory copy
for t in v
println("$(t.first)\t$(t.second)")
end
end

main()



On Monday, November 30, 2015 at 5:31:20 PM UTC+2, Attila Zséder wrote:
>
> Hi,
>
>
> The data I'm using is part of (Hungarian) Wikipedia dump with 5M lines of 
> text. On this data, python runs for 65 seconds, cpp for 35 seconds, julia 
> baseline for 340 seconds, julia with FastAnonymous.jl for 280 seconds. (See 
> https://github.com/juditacs/wordcount#leaderboard for details)
>
> Dan:
> I can use external packages, it's not a big issue. However, FastAnonymous 
> didn't give results comparable to python.
> The baseline python code I compare to is here:
> https://github.com/juditacs/wordcount/blob/master/python/wordcount_py2.py
>
> 2) The community is part of the language, so it should be regarded when 
> making considerations.
> What do you mean by this?
> My (our) purpose is not to evaluate this language or that one is 
> better/faster/etc because it is faster in unoptimized word counting. So I 
> don't want to make any judgements, considerations and anything like this. 
> This is just for fun. And even though it looks like _my_ julia 
> implementation of wc is not fast right now, I didn't lose interest in 
> following what's going on with this language.
>
>
> Your other points:
> 1) I do this with all the other languages as well. The test runs for about 
> 30-300 seconds. If Julia load time or any other thing takes serious amount 
> of time, then it does. This test is not precise, I didn't include c++ 
> compile time for example, but it took less than a second. But I felt like 
> my implementation is dummy, and other things take my time, not Julia load.
> 2) What if my test is about IO + dictionary storage? Then I have to 
> include printouts into my test.
> 3) I think 5m lines of text file is enough to avoid this noises.
>
>
>
> Tim:
> Yes, I did this code split, and with larger files it looked like after 
> sorting, dictionary manipulation (including hashes) took most of the time, 
> and printing was less of an issue. But I do have to analyze this more 
> precisely, seeing your numbers.
>
>
> Thank you all for your help!
>
> Attila
>
> On Mon, Nov 30, 2015 at 4:20 PM, Tim Holy 
> > wrote:
>
>> If you don't want to figure out how to use the profiler, your next best 
>> bet is
>> to split out the pieces so you can understand where the bottleneck is. For
>> example:
>>
>> function docount(io)
>> wc = Dict{AbstractString,Int64}()
>> for l in eachline(io)
>> for w in split(l)
>> wc[w]=get(wc, w, 0) + 1
>> end
>> end
>> wc
>> end
>>
>> @time open("somefile.tex") do io
>>docount(io)
>>end;
>>   0.010617 seconds (27.70 k allocations: 1.459 MB)
>>
>> vs
>>
>> @time open("somefile.tex") do io
>>main(io)
>>end;
>> # < lots of printed output >
>>   1.233154 seconds (330.59 k allocations: 10.829 MB, 1.53% gc time)
>>
>> (I modified your `main` to take an io input.)
>>
>> So it's the sorting and printing which is taking 99% of the time. Most of 
>> that
>> turns out to be the printing.
>>
>> --Tim
>>
>>
>

Re: [julia-users] Re: slow wordcount

[Dan]

My suggestions would be to replace
for t in sort(collect(wc), by=x -> (-x.second, x.first)) 
println(t.first, "\t", t.second)
end

with
customlt(a,b) = (b.second < a.second) ? true : b.second == a.second ? a.first 
< b.first : false

function main()
:
:
for t in sort(collect(wc), lt=customlt)
println(t.first, "\t", t.second)
end
end



On Monday, November 30, 2015 at 5:08:56 PM UTC+2, Dan wrote:
>
> Can you provide the comparable python code? Perhaps even the data used for 
> testing?
>
> Since you are evaluating Julia, there are two important points to remember:
> 1) In Julia because the language is fast enough to implement basic 
> functionality in Julia, then the distinction between Base Julia and 
> additional packages is small. Opting to use 'just' the core makes less 
> sense - the core is just a pre-compiled package.
> 2) The community is part of the language, so it should be regarded when 
> making considerations.
>
> On Monday, November 30, 2015 at 4:21:51 PM UTC+2, Attila Zséder wrote:
>>
>> Hi,
>>
>> Thank you all for the responses.
>>
>> 1. I tried simple profiling, but its output was difficult me to 
>> interpret, maybe if i put more time in it. I will try ProfileView later.
>> 2. FastAnonymous gave me a serious speedup (20-30%). (But since it is an 
>> external dependency, it's kind of cheating, seeing the purpose of this 
>> small word count test)
>> 3. Using ASCIIString is not a good option right now, since there are 
>> unicode characters there. I am trying with both UTF8String and 
>> AbstractString, I don't see any difference in performance right now.
>> 4. Using ht_keyindex() is out of scope for me right now, because this is 
>> a pet project, I just wanted to see how fast current implementation is, 
>> without these kind of solutions.
>>
>> I think I will keep trying with later versions of julia, but with 
>> sticking to the standard library only, without using any external packages.
>>
>> Attila
>>
>> 2015. november 29., vasárnap 17:59:42 UTC+1 időpontban Yichao Yu a 
>> következőt írta:
>>>
>>> On Sun, Nov 29, 2015 at 11:42 AM, Milan Bouchet-Valat  
>>> wrote: 
>>> > Le dimanche 29 novembre 2015 à 08:28 -0800, Cedric St-Jean a écrit : 
>>> >> What I would try: 
>>> >> 
>>> >> 1. ProfileView to pinpoint the bottleneck further 
>>> >> 2. FastAnonymous to fix the lambda 
>>> >> 3. http://julia-demo.readthedocs.org/en/latest/manual/performance-tip 
>>> >> s.html In particular, you may check `code_typed`. I don't have 
>>> >> experience with `split` and `eachline`. It's possible that they are 
>>> >> not type stable (the compiler can't predict their output's type). I 
>>> >> would try `for w::ASCIIString in ...` 
>>> >> 4. Dict{ASCIIString, Int}() 
>>> >> 5. Your loop will hash each string twice. I don't know how to fix 
>>> >> that, anyone? 
>>> > You can use the unexported Base.ht_keyindex() function like this: 
>>> > 
>>> https://github.com/nalimilan/FreqTables.jl/blob/7884c000e6797d7ec621e07 
>>> > b8da58e7939e39867/src/freqtable.jl#L36 
>>> > 
>>> > But this is at your own risk, as it may change without warning in a 
>>> > future Julia release. 
>>> > 
>>> > We really need a public API for it. 
>>>
>>> IIUC, https://github.com/JuliaLang/julia/issues/12157 
>>>
>>> > 
>>> > 
>>> > Regards 
>>> > 
>>> >> 
>>> >> Good luck, 
>>> >> 
>>> >> Cédric 
>>> >> 
>>> >> On Saturday, November 28, 2015 at 8:08:49 PM UTC-5, Lampkld wrote: 
>>> >> > Maybe it's the lambda? These are slow in julia right now. 
>>>
>>

Re: [julia-users] Re: slow wordcount

[Dan]

Some points:
1) If you are timing the program by running from shell (i.e. including the 
whole Julia load time) - it's not a good idea for a bunch of reasons.
2) The printouts are also a little dubious to be used in a test.
3) Too short an input makes the test noisy.


On Monday, November 30, 2015 at 4:21:51 PM UTC+2, Attila Zséder wrote:
>
> Hi,
>
> Thank you all for the responses.
>
> 1. I tried simple profiling, but its output was difficult me to interpret, 
> maybe if i put more time in it. I will try ProfileView later.
> 2. FastAnonymous gave me a serious speedup (20-30%). (But since it is an 
> external dependency, it's kind of cheating, seeing the purpose of this 
> small word count test)
> 3. Using ASCIIString is not a good option right now, since there are 
> unicode characters there. I am trying with both UTF8String and 
> AbstractString, I don't see any difference in performance right now.
> 4. Using ht_keyindex() is out of scope for me right now, because this is a 
> pet project, I just wanted to see how fast current implementation is, 
> without these kind of solutions.
>
> I think I will keep trying with later versions of julia, but with sticking 
> to the standard library only, without using any external packages.
>
> Attila
>
> 2015. november 29., vasárnap 17:59:42 UTC+1 időpontban Yichao Yu a 
> következőt írta:
>>
>> On Sun, Nov 29, 2015 at 11:42 AM, Milan Bouchet-Valat  
>> wrote: 
>> > Le dimanche 29 novembre 2015 à 08:28 -0800, Cedric St-Jean a écrit : 
>> >> What I would try: 
>> >> 
>> >> 1. ProfileView to pinpoint the bottleneck further 
>> >> 2. FastAnonymous to fix the lambda 
>> >> 3. http://julia-demo.readthedocs.org/en/latest/manual/performance-tip 
>> >> s.html In particular, you may check `code_typed`. I don't have 
>> >> experience with `split` and `eachline`. It's possible that they are 
>> >> not type stable (the compiler can't predict their output's type). I 
>> >> would try `for w::ASCIIString in ...` 
>> >> 4. Dict{ASCIIString, Int}() 
>> >> 5. Your loop will hash each string twice. I don't know how to fix 
>> >> that, anyone? 
>> > You can use the unexported Base.ht_keyindex() function like this: 
>> > https://github.com/nalimilan/FreqTables.jl/blob/7884c000e6797d7ec621e07 
>> > b8da58e7939e39867/src/freqtable.jl#L36 
>> > 
>> > But this is at your own risk, as it may change without warning in a 
>> > future Julia release. 
>> > 
>> > We really need a public API for it. 
>>
>> IIUC, https://github.com/JuliaLang/julia/issues/12157 
>>
>> > 
>> > 
>> > Regards 
>> > 
>> >> 
>> >> Good luck, 
>> >> 
>> >> Cédric 
>> >> 
>> >> On Saturday, November 28, 2015 at 8:08:49 PM UTC-5, Lampkld wrote: 
>> >> > Maybe it's the lambda? These are slow in julia right now. 
>>
>

Re: [julia-users] Re: slow wordcount

[Dan]

Can you provide the comparable python code? Perhaps even the data used for 
testing?

Since you are evaluating Julia, there are two important points to remember:
1) In Julia because the language is fast enough to implement basic 
functionality in Julia, then the distinction between Base Julia and 
additional packages is small. Opting to use 'just' the core makes less 
sense - the core is just a pre-compiled package.
2) The community is part of the language, so it should be regarded when 
making considerations.

On Monday, November 30, 2015 at 4:21:51 PM UTC+2, Attila Zséder wrote:
>
> Hi,
>
> Thank you all for the responses.
>
> 1. I tried simple profiling, but its output was difficult me to interpret, 
> maybe if i put more time in it. I will try ProfileView later.
> 2. FastAnonymous gave me a serious speedup (20-30%). (But since it is an 
> external dependency, it's kind of cheating, seeing the purpose of this 
> small word count test)
> 3. Using ASCIIString is not a good option right now, since there are 
> unicode characters there. I am trying with both UTF8String and 
> AbstractString, I don't see any difference in performance right now.
> 4. Using ht_keyindex() is out of scope for me right now, because this is a 
> pet project, I just wanted to see how fast current implementation is, 
> without these kind of solutions.
>
> I think I will keep trying with later versions of julia, but with sticking 
> to the standard library only, without using any external packages.
>
> Attila
>
> 2015. november 29., vasárnap 17:59:42 UTC+1 időpontban Yichao Yu a 
> következőt írta:
>>
>> On Sun, Nov 29, 2015 at 11:42 AM, Milan Bouchet-Valat  
>> wrote: 
>> > Le dimanche 29 novembre 2015 à 08:28 -0800, Cedric St-Jean a écrit : 
>> >> What I would try: 
>> >> 
>> >> 1. ProfileView to pinpoint the bottleneck further 
>> >> 2. FastAnonymous to fix the lambda 
>> >> 3. http://julia-demo.readthedocs.org/en/latest/manual/performance-tip 
>> >> s.html In particular, you may check `code_typed`. I don't have 
>> >> experience with `split` and `eachline`. It's possible that they are 
>> >> not type stable (the compiler can't predict their output's type). I 
>> >> would try `for w::ASCIIString in ...` 
>> >> 4. Dict{ASCIIString, Int}() 
>> >> 5. Your loop will hash each string twice. I don't know how to fix 
>> >> that, anyone? 
>> > You can use the unexported Base.ht_keyindex() function like this: 
>> > https://github.com/nalimilan/FreqTables.jl/blob/7884c000e6797d7ec621e07 
>> > b8da58e7939e39867/src/freqtable.jl#L36 
>> > 
>> > But this is at your own risk, as it may change without warning in a 
>> > future Julia release. 
>> > 
>> > We really need a public API for it. 
>>
>> IIUC, https://github.com/JuliaLang/julia/issues/12157 
>>
>> > 
>> > 
>> > Regards 
>> > 
>> >> 
>> >> Good luck, 
>> >> 
>> >> Cédric 
>> >> 
>> >> On Saturday, November 28, 2015 at 8:08:49 PM UTC-5, Lampkld wrote: 
>> >> > Maybe it's the lambda? These are slow in julia right now. 
>>
>

Re: [julia-users] Re: Arrays of arrays.

[Dan]

comprehension is indeed fine (and fast). should remember to read/search the 
docs and issues more (as usual).
fill(function,dims) is also good.

wondering, if it might be productive to make fill warn or even err when 
there is ambiguity and make versions:
`shallowfill`,`copyfill`,`deepcopyfill`. all this is unnecessary, but might 
help for readability purposes.

On Sunday, November 29, 2015 at 3:49:19 PM UTC+2, Milan Bouchet-Valat wrote:
>
> Le dimanche 29 novembre 2015 à 04:21 -0800, Dan a écrit : 
> > Since this situation comes up once in a while, it might be nice to 
> > have a copyfill which does a copy on each cell. In code: 
> > 
> > function copyfill!{T}(a::Array{T}, x) 
> > xT = convert(T, x) 
> > for i in eachindex(a) 
> > @inbounds a[i] = copy(xT) 
> > end 
> > return a 
> > end 
> > 
> > copyfill(v, dims::Dims)   = copyfill!(Array(typeof(v), dims), v) 
> > copyfill(v, dims::Integer...) = copyfill!(Array(typeof(v), dims...), 
> > v) 
> > 
> > And then we have: 
> > 
> > julia> m = copyfill(Int[1],3,3) 
> > 3x3 Array{Array{Int64,1},2}: 
> >  [1]  [1]  [1] 
> >  [1]  [1]  [1] 
> >  [1]  [1]  [1] 
> > 
> > julia> push!(m[1,1],2) 
> > 2-element Array{Int64,1}: 
> >  1 
> >  2 
> > 
> > julia> m 
> > 3x3 Array{Array{Int64,1},2}: 
> >  [1,2]  [1]  [1] 
> >  [1][1]  [1] 
> >  [1][1]  [1] 
> > 
> > As some would expect. Is there a +1 on this? 
> This has been discussed in depth here: 
> https://github.com/JuliaLang/julia/pull/8759 
>
> One issue is that sometimes copy() is enough, but sometimes you want 
> deepcopy() instead (when the object itself contains references). Thus, 
> a more general solution would be to take as first argument a function 
> returning the value: 
> https://groups.google.com/d/msg/julia-users/L5fXkHPduBo/UNoq6NTQOnEJ 
>
> But as Stefan noted the comprehension form is quite good too. Maybe the 
> docs should point to it after the warning about shared references. 
>
>
> Regards 
>
> > On Sunday, November 29, 2015 at 1:36:03 PM UTC+2, Kristoffer Carlsson 
> > wrote: 
> > > I guess the simplest would be: 
> > > 
> > > [Int[] for i = 1:3, j=1:3, k=1:3] 
> > > 
> > > 
> > > And to repeat what Milan already said, you don't want fill! because 
> > > then all your arrays point to the same memory location. 
> > > 
> > > On Sunday, November 29, 2015 at 11:51:28 AM UTC+1, Aleksandr 
> > > Mikheev wrote: 
> > > > Hi all. Once again I have some questions about Julia. 
> > > > 
> > > > I know that there is a possibility to create a list of arrays. 
> > > > For exmple: 
> > > > 
> > > > s = fill(Array(Int64,1),4) 
> > > > 
> > > > And then I can do something like this: 
> > > > 
> > > > s[1] = [1; 2] 
> > > > s[1] = [s[1]; 5] 
> > > > 
> > > > By parity of reasoning I did this: 
> > > > 
> > > > s = fill(Array(Int64,1),4,4,4) 
> > > > 
> > > > And it worked fine. But in both cases I had initial elements in s 
> > > > (like when I construct arrays with Array{Int64}(m,n)): 
> > > > 
> > > > julia> s = fill(Array(Int64,1),3,3,3) 
> > > > 3x3x3 Array{Array{Int64,1},3}: 
> > > > [:, :, 1] = 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > > 
> > > > 
> > > > [:, :, 2] = 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > > 
> > > > 
> > > > [:, :, 3] = 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > >  [2221816832]  [2221816832]  [2221816832] 
> > > > 
> > > > 
> > > > Is there something I could do to prevent this? I know I could 
> > > > easily fix it by: 
> > > > 
> > > > 
> > > > for i = 1:3 
> > > > for j = 1:3 
> > > > for k = 1:3 
> > > > s[i,j,k] = [] 
> > > > end 
> > > > end 
> > > > end 
> > > > 
> > > > But I guess this is a weird solution. 
> > > > 
> > > > Thank you in advance! 
> > > > 
> > > > 
>

[julia-users] Re: Arrays of arrays.

[Dan]

Since this situation comes up once in a while, it might be nice to have a 
copyfill which does a copy on each cell. In code:

function copyfill!{T}(a::Array{T}, x)
xT = convert(T, x)
for i in eachindex(a)
@inbounds a[i] = copy(xT)
end
return a
end

copyfill(v, dims::Dims)   = copyfill!(Array(typeof(v), dims), v)
copyfill(v, dims::Integer...) = copyfill!(Array(typeof(v), dims...), v)

And then we have:

julia> m = copyfill(Int[1],3,3)
3x3 Array{Array{Int64,1},2}:
 [1]  [1]  [1]
 [1]  [1]  [1]
 [1]  [1]  [1]

julia> push!(m[1,1],2)
2-element Array{Int64,1}:
 1
 2

julia> m
3x3 Array{Array{Int64,1},2}:
 [1,2]  [1]  [1]
 [1][1]  [1]
 [1][1]  [1]

As some would expect. Is there a +1 on this?

On Sunday, November 29, 2015 at 1:36:03 PM UTC+2, Kristoffer Carlsson wrote:
>
> I guess the simplest would be:
>
> [Int[] for i = 1:3, j=1:3, k=1:3]
>
>
> And to repeat what Milan already said, you don't want fill! because then 
> all your arrays point to the same memory location.
>
> On Sunday, November 29, 2015 at 11:51:28 AM UTC+1, Aleksandr Mikheev wrote:
>>
>> Hi all. Once again I have some questions about Julia.
>>
>> I know that there is a possibility to create a list of arrays. For exmple:
>>
>> s = fill(Array(Int64,1),4)
>>
>> And then I can do something like this:
>>
>> s[1] = [1; 2]
>> s[1] = [s[1]; 5]
>>
>> By parity of reasoning I did this:
>>
>> s = fill(Array(Int64,1),4,4,4)
>>
>> And it worked fine. But in both cases I had initial elements in s (like 
>> when I construct arrays with Array{Int64}(m,n)):
>>
>> julia> s = fill(Array(Int64,1),3,3,3)
>> 3x3x3 Array{Array{Int64,1},3}:
>> [:, :, 1] =
>>  [2221816832]  [2221816832]  [2221816832]
>>  [2221816832]  [2221816832]  [2221816832]
>>  [2221816832]  [2221816832]  [2221816832]
>>
>>
>> [:, :, 2] =
>>  [2221816832]  [2221816832]  [2221816832]
>>  [2221816832]  [2221816832]  [2221816832]
>>  [2221816832]  [2221816832]  [2221816832]
>>
>>
>> [:, :, 3] =
>>  [2221816832]  [2221816832]  [2221816832]
>>  [2221816832]  [2221816832]  [2221816832]
>>  [2221816832]  [2221816832]  [2221816832]
>>
>>
>> Is there something I could do to prevent this? I know I could easily fix 
>> it by:
>>
>>
>> for i = 1:3
>> for j = 1:3
>> for k = 1:3
>> s[i,j,k] = []
>> end
>> end
>> end
>>
>> But I guess this is a weird solution.
>>
>> Thank you in advance!
>>
>>

[julia-users] Re: is it possible for collect() to return a collection of a specific type?

[Dan]

Looking at the collect code, it seems you really should define Base.eltype for 
your iterator and things will work out.

On Saturday, November 28, 2015 at 3:15:18 AM UTC+2, Seth wrote:
>
> I guess that makes sense, though I struggle to see why one would create an 
> iterator that produces multiple types.
>
> On Friday, November 27, 2015 at 4:45:26 PM UTC-8, ele...@gmail.com wrote:
>>
>> Collect only requires that the collection is iterable, for which 
>> providing an eltype() function is optional.  I don't know if it is possible 
>> to check at runtime if eltype() exists for the collection then it could use 
>> that instead of Any, otherwise it would have to iterate the collection to 
>> find all the types and either accumulate the results in an Any collection 
>> as it goes and copy them to the right type collection later, or iterate 
>> twice.
>>
>> On Saturday, November 28, 2015 at 10:17:09 AM UTC+10, Seth wrote:
>>>
>>> Well, I just found collect(Edge, edges(g)) works, but it would be nice 
>>> if collect() returned a vector of Edge by default. Any ideas?
>>>
>>> On Friday, November 27, 2015 at 3:46:58 PM UTC-8, Seth wrote:

 I have implemented my own iterator that works against edges in a graph, 
 and edges(g) returns the iterator. However, collect(edges(g)) returns 
 an array of Any,1. I'd like it to return an array of Edge, 1. What am 
 I missing?

[julia-users] Re: is it possible for collect() to return a collection of a specific type?

[Dan]

have you tried running collect within a function? sometimes the Global 
scope is bad for type inference.

On Saturday, November 28, 2015 at 2:45:26 AM UTC+2, ele...@gmail.com wrote:
>
> Collect only requires that the collection is iterable, for which providing 
> an eltype() function is optional.  I don't know if it is possible to check 
> at runtime if eltype() exists for the collection then it could use that 
> instead of Any, otherwise it would have to iterate the collection to find 
> all the types and either accumulate the results in an Any collection as it 
> goes and copy them to the right type collection later, or iterate twice.
>
> On Saturday, November 28, 2015 at 10:17:09 AM UTC+10, Seth wrote:
>>
>> Well, I just found collect(Edge, edges(g)) works, but it would be nice 
>> if collect() returned a vector of Edge by default. Any ideas?
>>
>> On Friday, November 27, 2015 at 3:46:58 PM UTC-8, Seth wrote:
>>>
>>> I have implemented my own iterator that works against edges in a graph, 
>>> and edges(g) returns the iterator. However, collect(edges(g)) returns 
>>> an array of Any,1. I'd like it to return an array of Edge, 1. What am I 
>>> missing?
>>>
>>>

[julia-users] Re: Precompilation and functions with keyword arguments

[Dan]

Further investigation suggests recompilation is happening. This raises 
questions:
1) why is there recompilation?
2) why are the memory allocations associated not reported?
3) is there some memory leak?

The following run suggests the recompilation:

  | | |_| | | | (_| |  |  Version 0.4.2-pre+15 (2015-11-20 12:12 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit c4bbb89* (5 days old release-0.4)
|__/   |  x86_64-linux-gnu


julia> testfunction2(n) = (x = 0.0; for i=1:n x += sin(i) ; end ; x)
testfunction2 (generic function with 1 method)


julia> @time precompile(testfunction2,(Int64,))
  0.004197 seconds (2.84 k allocations: 146.004 KB)


julia> a = Base._dump_function(testfunction2,(Int64,),false,false,true,false
)
"\ndefine double @julia_testfunction2_21229(i64) {\ntop:\n  %1 = icmp sgt 
i64 %0, 0\n  %2 = select i1 %1, i64 %0, i64 0\n  %3 = icmp eq i64 %2, 0\n 
 br i1 %3, label %L3, label %L\n\nL:   
 ; preds = %pass, %top\n  %x.0 = phi double [ %11, %pass ], [ 
0.00e+00, %top ]\n  %\"#s1.0\" = phi i64 [ %10, %pass ], [ 1, %top ]\n 
 %4 = sitofp i64 %\"#s1.0\" to double\n  %5 = call double inttoptr (i64 
139659687619792 to double (double)*)(double inreg %4)\n  %6 = fcmp uno 
double %4, 0.00e+00\n  %7 = fcmp ord double %5, 0.00e+00\n  %8 = or 
i1 %7, %6\n  br i1 %8, label %pass, label %fail\n\nfail:   
  ; preds = %L\n  %9 = load %jl_value_t** 
@jl_domain_exception, align 8\n  call void 
@jl_throw_with_superfluous_argument(%jl_value_t* %9, i32 1)\n 
 unreachable\n\npass: ; preds = 
%L\n  %10 = add i64 %\"#s1.0\", 1\n  %11 = fadd double %x.0, %5\n  %12 = 
icmp eq i64 %\"#s1.0\", %2\n  br i1 %12, label %L3, label %L\n\nL3: 
  ; preds = %pass, %top\n  %x.1 = phi 
double [ 0.00e+00, %top ], [ %11, %pass ]\n  ret double %x.1\n}\n"


julia> @time testfunction2(100_000)
  0.005485 seconds (5 allocations: 176 bytes)
1.84103630354


julia> a = Base._dump_function(testfunction2,(Int64,),false,false,true,false
)
"\ndefine double @julia_testfunction2_21293(i64) {\ntop:\n  %1 = icmp sgt 
i64 %0, 0\n  %2 = select i1 %1, i64 %0, i64 0\n  %3 = icmp eq i64 %2, 0\n 
 br i1 %3, label %L3, label %L\n\nL:   
 ; preds = %pass, %top\n  %x.0 = phi double [ %11, %pass ], [ 
0.00e+00, %top ]\n  %\"#s1.0\" = phi i64 [ %10, %pass ], [ 1, %top ]\n 
 %4 = sitofp i64 %\"#s1.0\" to double\n  %5 = call double inttoptr (i64 
139659687619792 to double (double)*)(double inreg %4)\n  %6 = fcmp uno 
double %4, 0.00e+00\n  %7 = fcmp ord double %5, 0.00e+00\n  %8 = or 
i1 %7, %6\n  br i1 %8, label %pass, label %fail\n\nfail:   
  ; preds = %L\n  %9 = load %jl_value_t** 
@jl_domain_exception, align 8\n  call void 
@jl_throw_with_superfluous_argument(%jl_value_t* %9, i32 1)\n 
 unreachable\n\npass: ; preds = 
%L\n  %10 = add i64 %\"#s1.0\", 1\n  %11 = fadd double %x.0, %5\n  %12 = 
icmp eq i64 %\"#s1.0\", %2\n  br i1 %12, label %L3, label %L\n\nL3: 
  ; preds = %pass, %top\n  %x.1 = phi 
double [ 0.00e+00, %top ], [ %11, %pass ]\n  ret double %x.1\n}\n"


As can be seen, the first version is julia_testfunction2_21229 and the 
second is julia_testfunction2_21293. The signatures are the same.

Disclaimer - I'm very unfamiliar with the Julia codegen code (and much of 
the rest).
Intuition: It might be related to the following snippet from src/codegen.cpp
:

extern "C" DLLEXPORT
void *jl_get_llvmf(jl_function_t *f, jl_tupletype_t *tt, bool getwrapper)
{
:
:
:
if (sf->linfo->functionObject != NULL) {
// found in the system image: force a recompile
Function *llvmf = (Function*)sf->linfo->functionObject;
if (llvmf->isDeclaration()) {
sf->linfo->specFunctionObject = NULL;
sf->linfo->functionObject = NULL;
}
}


Note the force a recompile bit.


On Wednesday, November 25, 2015 at 10:23:04 PM UTC+2, Dan wrote:
>
> I've managed to reproduce the timing peculiarities you described, with a 
> different function. In fact, I get the same first run slower behavior even 
> with regular functions with -no- keywords. It is almost as if `precompile` 
> is not doing the entire job.
> Have you tried, just testing the basic precompile-1st run-2nd run timing 
> of a no keyword generic function? 
>
> On Tuesday, November 24, 2015 at 9:43:41 PM UTC+2, Tim Holy wrote:
>>
>> I've been experimenting further with SnoopCompile and Immerse/Gadfly, 
>> trying to 
>> shave off more of the time-to-first-plot. If I 

[julia-users] Re: Precompilation and functions with keyword arguments

[Dan]

Aha! tested the same sequence in a non-REPL run and the 1st and 2nd 
runtimes after a `precompile` are the same.

On Wednesday, November 25, 2015 at 10:23:04 PM UTC+2, Dan wrote:
>
> I've managed to reproduce the timing peculiarities you described, with a 
> different function. In fact, I get the same first run slower behavior even 
> with regular functions with -no- keywords. It is almost as if `precompile` 
> is not doing the entire job.
> Have you tried, just testing the basic precompile-1st run-2nd run timing 
> of a no keyword generic function? 
>
> On Tuesday, November 24, 2015 at 9:43:41 PM UTC+2, Tim Holy wrote:
>>
>> I've been experimenting further with SnoopCompile and Immerse/Gadfly, 
>> trying to 
>> shave off more of the time-to-first-plot. If I pull out all the stops 
>> (using the 
>> "userimg.jl" strategy), I can get it down to about 2.5 seconds. However, 
>> doing 
>> the same plot a second time is about 0.02s. This indicates that despite 
>> my 
>> efforts, there's still a lot that's not being cached. 
>>
>> About 0.3s of that total (not much, but it's all I have data on) can be 
>> observed via snooping as re-compilation of functions that you might 
>> imagine 
>> should have been precompiled. The biggest offenders are all functions 
>> with 
>> keyword arguments. In miniature, I think you can see this here: 
>>
>> julia> function foo(X; thin=true) 
>>svdfact(X) 
>>end 
>> foo (generic function with 1 method) 
>>
>> # Before compiling this, let's make sure the work compiling foo isn't 
>> hidden 
>> # by other compilation needs: 
>> julia> A = rand(3,3) 
>> 3x3 Array{Float64,2}: 
>>  0.570780.33557   0.56497   
>>  0.0679035  0.944406  0.816098 
>>  0.0922775  0.404697  0.0900726 
>>
>> julia> svdfact(A) 
>> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>>  -0.507226   0.861331   0.0288001 
>>  -0.825227  -0.475789  -0.304344 
>>  -0.248438  -0.178138   0.952127 , 
>> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
>> Array{Float64,2}: 
>>  -0.248222  -0.707397  -0.661797 
>>   0.873751  -0.458478   0.162348 
>>   0.418264   0.537948  -0.731893) 
>>
>> # OK, let's precompile foo 
>> julia> @time precompile(foo, (Matrix{Float64},)) 
>>   0.000469 seconds (541 allocations: 35.650 KB) 
>>
>> julia> @time foo(A) 
>>   0.001174 seconds (18 allocations: 3.063 KB) 
>> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>>  -0.507226   0.861331   0.0288001 
>>  -0.825227  -0.475789  -0.304344 
>>  -0.248438  -0.178138   0.952127 , 
>> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
>> Array{Float64,2}: 
>>  -0.248222  -0.707397  -0.661797 
>>   0.873751  -0.458478   0.162348 
>>   0.418264   0.537948  -0.731893) 
>>
>> # Note the 2nd call is 10x faster, despite precompilation 
>> julia> @time foo(A) 
>>   0.000164 seconds (18 allocations: 3.063 KB) 
>> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>>  -0.507226   0.861331   0.0288001 
>>  -0.825227  -0.475789  -0.304344 
>>  -0.248438  -0.178138   0.952127 , 
>> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
>> Array{Float64,2}: 
>>  -0.248222  -0.707397  -0.661797 
>>   0.873751  -0.458478   0.162348 
>>   0.418264   0.537948  -0.731893) 
>>
>> # Note adding a keyword argument to the call causes a further 10x 
>> slowdown... 
>> julia> @time foo(A; thin=true) 
>>   0.014787 seconds (3.36 k allocations: 166.622 KB) 
>> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>>  -0.507226   0.861331   0.0288001 
>>  -0.825227  -0.475789  -0.304344 
>>  -0.248438  -0.178138   0.952127 , 
>> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
>> Array{Float64,2}: 
>>  -0.248222  -0.707397  -0.661797 
>>   0.873751  -0.458478   0.162348 
>>   0.418264   0.537948  -0.731893) 
>>
>> # ...but only for the first call 
>> julia> @time foo(A; thin=true) 
>>   0.000209 seconds (19 allocations: 3.141 KB) 
>> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>>  -0.507226   0.861331   0.0288001 
>>  -0.825227  -0.475789  -0.304344 
>>  -0.248438  -0.178138   0.952127 , 
>> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
>> Array{Float64,2}: 
>>  -0.248222  -0.707397  -0.661797 
>>   0.873751  -0.458478   0.162348 
>>   0.418264   0.537948  -0.731893) 
>>
>>
>> Obviously the times here don't add up to much, but for a project the size 
>> of 
>> Gadfly it might matter. 
>>
>> I should add that 
>> precompile(foo, (Vector{Any}, Matrix{Float64})) 
>> doesn't seem to do anything useful. 
>>
>> Any ideas? 
>>
>> Best, 
>> --Tim 
>>
>

[julia-users] Re: Precompilation and functions with keyword arguments

[Dan]

I've managed to reproduce the timing peculiarities you described, with a 
different function. In fact, I get the same first run slower behavior even 
with regular functions with -no- keywords. It is almost as if `precompile` 
is not doing the entire job.
Have you tried, just testing the basic precompile-1st run-2nd run timing of 
a no keyword generic function? 

On Tuesday, November 24, 2015 at 9:43:41 PM UTC+2, Tim Holy wrote:
>
> I've been experimenting further with SnoopCompile and Immerse/Gadfly, 
> trying to 
> shave off more of the time-to-first-plot. If I pull out all the stops 
> (using the 
> "userimg.jl" strategy), I can get it down to about 2.5 seconds. However, 
> doing 
> the same plot a second time is about 0.02s. This indicates that despite my 
> efforts, there's still a lot that's not being cached. 
>
> About 0.3s of that total (not much, but it's all I have data on) can be 
> observed via snooping as re-compilation of functions that you might 
> imagine 
> should have been precompiled. The biggest offenders are all functions with 
> keyword arguments. In miniature, I think you can see this here: 
>
> julia> function foo(X; thin=true) 
>svdfact(X) 
>end 
> foo (generic function with 1 method) 
>
> # Before compiling this, let's make sure the work compiling foo isn't 
> hidden 
> # by other compilation needs: 
> julia> A = rand(3,3) 
> 3x3 Array{Float64,2}: 
>  0.570780.33557   0.56497   
>  0.0679035  0.944406  0.816098 
>  0.0922775  0.404697  0.0900726 
>
> julia> svdfact(A) 
> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>  -0.507226   0.861331   0.0288001 
>  -0.825227  -0.475789  -0.304344 
>  -0.248438  -0.178138   0.952127 , 
> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
> Array{Float64,2}: 
>  -0.248222  -0.707397  -0.661797 
>   0.873751  -0.458478   0.162348 
>   0.418264   0.537948  -0.731893) 
>
> # OK, let's precompile foo 
> julia> @time precompile(foo, (Matrix{Float64},)) 
>   0.000469 seconds (541 allocations: 35.650 KB) 
>
> julia> @time foo(A) 
>   0.001174 seconds (18 allocations: 3.063 KB) 
> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>  -0.507226   0.861331   0.0288001 
>  -0.825227  -0.475789  -0.304344 
>  -0.248438  -0.178138   0.952127 , 
> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
> Array{Float64,2}: 
>  -0.248222  -0.707397  -0.661797 
>   0.873751  -0.458478   0.162348 
>   0.418264   0.537948  -0.731893) 
>
> # Note the 2nd call is 10x faster, despite precompilation 
> julia> @time foo(A) 
>   0.000164 seconds (18 allocations: 3.063 KB) 
> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>  -0.507226   0.861331   0.0288001 
>  -0.825227  -0.475789  -0.304344 
>  -0.248438  -0.178138   0.952127 , 
> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
> Array{Float64,2}: 
>  -0.248222  -0.707397  -0.661797 
>   0.873751  -0.458478   0.162348 
>   0.418264   0.537948  -0.731893) 
>
> # Note adding a keyword argument to the call causes a further 10x 
> slowdown... 
> julia> @time foo(A; thin=true) 
>   0.014787 seconds (3.36 k allocations: 166.622 KB) 
> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>  -0.507226   0.861331   0.0288001 
>  -0.825227  -0.475789  -0.304344 
>  -0.248438  -0.178138   0.952127 , 
> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
> Array{Float64,2}: 
>  -0.248222  -0.707397  -0.661797 
>   0.873751  -0.458478   0.162348 
>   0.418264   0.537948  -0.731893) 
>
> # ...but only for the first call 
> julia> @time foo(A; thin=true) 
>   0.000209 seconds (19 allocations: 3.141 KB) 
> Base.LinAlg.SVD{Float64,Float64,Array{Float64,2}}(3x3 Array{Float64,2}: 
>  -0.507226   0.861331   0.0288001 
>  -0.825227  -0.475789  -0.304344 
>  -0.248438  -0.178138   0.952127 , 
> [1.4844598265207638,0.5068781079827415,0.19995120630810712],3x3 
> Array{Float64,2}: 
>  -0.248222  -0.707397  -0.661797 
>   0.873751  -0.458478   0.162348 
>   0.418264   0.537948  -0.731893) 
>
>
> Obviously the times here don't add up to much, but for a project the size 
> of 
> Gadfly it might matter. 
>
> I should add that 
> precompile(foo, (Vector{Any}, Matrix{Float64})) 
> doesn't seem to do anything useful. 
>
> Any ideas? 
>
> Best, 
> --Tim 
>

Re: [julia-users] Concatenation without splatting

[Dan]

Using `append!` instead of `push!` and letting efficient dynamic 
reallocation of vector do the resizing:

function vcat_nosplat2a(y)
result = Array(eltype(y[1]), 0)
for a in y
append!(result, a)
end
result
end

(@benchmark shows way less allocations and 2-3x time)
BTW the splat version is just as quick, but perhaps allocation on stack is 
problematic (anybody check the limit?)

On Monday, November 23, 2015 at 12:07:18 PM UTC+2, Tomas Lycken wrote:
>
> That doesn’t quite seem to do what you want, Mauro:
>
> julia> arr_of_arr = Vector{Int}[[1],[2,3],[4,5]]
> 3-element Array{Array{Int64,1},1}:
>  [1]
>  [2,3]
>  [4,5]
>
> julia> vcat_nosplat(y) = eltype(y[1])[el[1] for el in y]
> vcat_nosplat (generic function with 1 method)
>
> julia> vcat_nosplat(arr_of_arr)
> 3-element Array{Int64,1}:
>  1
>  2
>  4
>
> It’s quite trivial to achieve the desired result with only a few lines of 
> code, though:
>
> julia> function vcat_nosplat2(y)
>result = Array(eltype(y[1]), 0)
>sizehint!(result, sum(map(length, y))) #skip if iterating is more 
> expensive than reallcoation
>
>for a in y
>for x in a
>push!(result, x)
>end
>end
>
>result
>end
> vcat_nosplat2 (generic function with 1 method)
>
> julia> vcat_nosplat2(arr_of_arr)
> 5-element Array{Int64,1}:
>  1
>  2
>  3
>  4
>  5
>
> // T
>
> On Sunday, November 22, 2015 at 9:12:55 PM UTC+1, Mauro wrote:
>
> In ODE.jl, I've used 
>>
>> vcat_nosplat(y) = eltype(y[1])[el[1] for el in y] # Does vcat(y...) 
>> without the splatting 
>>
>> I think the eltype might not be needed.  There may be better ways though. 
>>
>> On Sun, 2015-11-22 at 14:04, Cedric St-Jean  wrote: 
>> > I have a big vector of vectors. Is there any way to vcat/hcat them 
>> without 
>> > splatting? 
>> > 
>> > arr_of_arr = Vector[[1],[2,3],[4,5]] 
>> > vcat(arr_of_arr...) 
>> > 
>> > I'm asking because splatting big arrays is a performance issue (and 
>> IIRC it 
>> > blows the stack at some point). 
>>
> ​
>

Re: [julia-users] Re: Better alternative to find all permutations?

[Dan]

If we are already at this, it would be cool to have versions with Gray ordering.
Something like:

   for (i1,i2) in permutations(vector,indices=true,Gray=true)
   theperm[[i1,i2]] = theperm[[i2,i1]]
   @show theperm
   end

would go over all permutations.
The syntax would be different and this is not a Gray order use-case.
But Gray ordering is useful and compact.

Re: [julia-users] Re: split a utf-8 string

[Dan]

`split` already has hooks for using other splitters. To achieve the UTF8 
space splitting functionality, one can leverage the `isspace` function and 
add some decorations. For example:

type FuncSplitter
   pred::Function
end
function Base.search(s::AbstractString, splt::FuncSplitter, start::Int = 1)
   i = start
   n = endof(s)
   while i<=endof(s) 
   if splt.pred(s[i])
   return i
   end
   i = nextind(s,i)
   end
   return 0
end
spacesplitter = FuncSplitter(x->isspace(x))

s = "Ｔｉｍｅ ｆｌｉｅｓ ｌｉｋｅ ａｎ ａｒｒｏｗ"
split(s,spacesplitter)

Gives:
5-element Array{SubString{UTF8String},1}:
 "Ｔｉｍｅ" 
 "ｆｌｉｅｓ"
 "ｌｉｋｅ" 
 "ａｎ"   
 "ａｒｒｏｗ"

This isn't fully optimized, but probably suffices for many uses.

On Sunday, November 22, 2015 at 1:29:46 PM UTC+2, Pontus Stenetorp wrote:
>
> On 22 November 2015 at 01:46,  > wrote: 
> > 
> > On Sunday, November 22, 2015 at 10:02:03 AM UTC+10, James Gilbert wrote: 
> >> 
> >> The spaces in your string are '\u3000' the ideographic space. 
> >> isspace('\u3000') returns true, and split(s) is supposed to split on 
> all 
> >> space characters, so I think this might be a julia bug. 
> > 
> > Or a documentation bug, the actual default is only the ASCII spaces 
> > https://github.com/JuliaLang/julia/blob/master/base/strings/util.jl#L62 
>
> It should probably be pointed out that at least Python3 (but not 
> Python2) gets it "right". 
>
> > python3 
> Python 3.4.3+ (default, Oct 14 2015, 16:03:50) 
> [GCC 5.2.1 20151010] on linux 
> Type "help", "copyright", "credits" or "license" for more information. 
> >>> "Ｔｉｍｅ ｆｌｉｅｓ ｌｉｋｅ ａｎ ａｒｒｏｗ".split() 
> ['Ｔｉｍｅ', 'ｆｌｉｅｓ', 'ｌｉｋｅ', 'ａｎ', 'ａｒｒｏｗ'] 
>
> I would argue that Unicode is a first class citizen and that Julia 
> should also get this "right".  This would require some fairly 
> straightforward, yet not trivial, tinkering and would be an excellent 
> first contribution if someone wants to take a stab at it. 
>
> Pontus 
>

[julia-users] Re: What is the best way to get two by two tables in Julia?

[Dan]

Hi Arin,
It would be helpful to have more details about the input (a dataframe?) and 
output (a two-by-two table or a table indexed by categories?). Some code to 
give context to the question would be even more help (possibly in another 
language, such as R).

Having said this, here is a starting point for some code:

If these packages are missing Pkg.add works:

using NamedArrays
using DataFrames
using RDatasets

Gets the dataset and makes some categorical variables in DataFrames style:

iris = dataset("datasets","iris")
iris[:PetalWidth] = PooledDataArray(iris[:PetalWidth])
iris[:Species] = PooledDataArray(iris[:Species])

Define function for a `twobytwo` and a general categorical table 
`crosstable`:

function twobytwo(data::DataFrame,cond1,cond2)
   nres= NamedArray(zeros(Int,2,2),Any[[false,true],[false,true]],[
"cond1","cond2"])
   for i=1:nrow(data)
   nres[Int(cond1(data[i,:]))+1,Int(cond2(data[i,:]))+1] += 1
   end
   nres
end

function crosstable(data::DataFrame,col1,col2)
   @assert isa(data[col1],PooledDataArray)
   @assert isa(data[col2],PooledDataArray)
   nres= NamedArray(zeros(Int,length(data[col1].pool),length(data[col2].
pool)),Any[data[col1].pool,data[col2].pool],[col1,col2])
   for i=1:nrow(data)
   nres[data[col1].refs[i],data[col2].refs[i]] += 1
   end
   nres
end

Finally, using the functions, make some tables:

tbt = twobytwo(iris,r->r[1,:Species]=="setosa",r->r[1,:PetalWidth]>=1.5)
ct = crosstable(iris,:PetalWidth,:Species)

My summary and conclusions:
1) Julia is general purpose and with a little familiarity any data handling 
is possible.
2) This is a basic data exploration operation and there must be some easy 
way to do this.

Waiting for more opinions/solutions on this question, as it is also basic 
for my needs.

Thanks for the question.

On Sunday, November 22, 2015 at 3:34:56 AM UTC+2, Arin Basu wrote:
>
> Hi All,
>
> Can you kindly advise how to get a simple way to do two by two tables in 
> Julia with two categorical variables. I have tried split-apply-combine (by 
> function) and it works with single variables, but with two or more 
> variables, I cannot get the table I want. 
>
> This is really an issue if we need to do statistical data analysis in 
> Epidemiology. 
>
> Any help or advice will be greatly appreciated.
>
> Arin Basu
>

[julia-users] Re: Coveralls/Codecov report certain lines are not covered, but test for the lines do appear to exist in code

[Dan]

Regarding the last printing point only
print(string(eltype(8)) * ", " * string(eltype(8)))

is probably what you meant. AbstractString(Int64), tries to create an 
AbstractString from an Int64, while string(Int64) asks for a string 
representation. The latter semantics is more accurate, and works.

On Sunday, November 22, 2015 at 7:46:17 AM UTC+2, Zygmunt Szpak wrote:
>
> Hi All,
>
> I've been perusing the Coveralls/Codecov reports looking for opportunities 
> to write some additional tests.
> I was looking at the promote_rule in int.jl and noticed that some 
> promote_rules were apparently missing tests. 
>
> See for example:
>  promote_rule(::Type{Int64}, ::Type{Int8} ) = Int64 
>
>  
>
> https://codecov.io/github/JuliaLang/julia/base/int.jl?ref=e6cd2bc8c6ab2444d7750cf42d43b4c45e8f3545
>
> Apparently, this line is not covered by any test as indicated by the red 
> highlight.
>
> However, the test file for int.jl
>
> https://github.com/JuliaLang/julia/blob/master/test/int.jl
>
> contains the following tests:
>
> UItypes = (UInt8, UInt16, UInt32, UInt64, UInt128)
> SItypes = (Int8, Int16, Int32, Int64, Int128)
>
>
> for T in UItypes, S in UItypes
> @test promote(S(3), T(3)) === (sizeof(T) < sizeof(S) ? (S(3), S(3)) : 
> (T(3), T(3)))
> end
>
>
> for T in SItypes, S in SItypes
> @test promote(S(3), T(3)) === (sizeof(T) < sizeof(S) ? (S(3), S(3)) : 
> (T(3), T(3)))
> end
>
>
> for T in SItypes, S in UItypes
> R = sizeof(S) < sizeof(Int) ? Int : S
> @test promote(R(3), T(3)) === (sizeof(R) < sizeof(T) ? (T(3), T(3)) : 
> (R(3), R(3)))
> end
>
> Looking deeper at one of these tests using the following code:
>
> for  T in SItypes, S in SItypes
> print(eltype(S(3)) )
> print(", ")
> print(eltype(T(3)) )
> print(" -> ")
> print(eltype(promote(S(3), T(3
> print(" | ")
> print(@which promote(S(3), T(3)))
> println(" ")
> end
>
> which yields:
>
> Int8, Int8 -> Int8 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int16, Int8 -> Int16 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int32, Int8 -> Int32 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> *Int64, Int8 -> Int64 | promote{T,S}(x::T, y::S) at promotion.jl:132 *
>> Int128, Int8 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int8, Int16 -> Int16 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int16, Int16 -> Int16 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int32, Int16 -> Int32 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int64, Int16 -> Int64 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int128, Int16 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int8, Int32 -> Int32 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int16, Int32 -> Int32 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int32, Int32 -> Int32 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int64, Int32 -> Int64 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int128, Int32 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int8, Int64 -> Int64 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int16, Int64 -> Int64 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int32, Int64 -> Int64 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int64, Int64 -> Int64 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int128, Int64 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int8, Int128 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int16, Int128 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int32, Int128 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int64, Int128 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>> Int128, Int128 -> Int128 | promote{T,S}(x::T, y::S) at promotion.jl:132 
>
>
> it would appear that the line *highlighted in blue* above is indeed 
> testing:
>  promote_rule(::Type{Int64}, ::Type{Int8} ) = Int64 
>
> contrary to what the Coveralls/Codecov is reporting. I'm not sure how 
> to reconcile these results, and would appreciate any insight.
>
> On a sightly different note, could anyone suggest a neat way of printing 
> element types on a single
> line such as:
>
>> Int64, Int64 -> Int64 
>
>
> For example, is there a way to convert Tuples to Strings? 
>
> I've tried things such as:
> print(AbstractString(eltype(8)) * ", " * AbstractString(eltype(8)))
>
> LoadError: MethodError: `convert` has no method matching 
>> convert(::Type{AbstractString}, ::Type{Int64})
>> This may have arisen from a call to the constructor AbstractString(...),
>> since type constructors fall back to convert methods.
>> Closest candidates are:
>> call{T}(::Type{T}, ::Any)
>> convert{T<:AbstractString,S<:Union{Char,Int32,UInt32}}(::Type{T<:AbstractString},
>>  
>> !Matched::AbstractArray{S<:Union{Char,Int32,UInt32},1})
>> convert{S<:AbstractString}(::Type{S<:AbstractString}, 
>> !Matched::Base.UTF8proc.GraphemeIterator{S<:AbstractString})
>> ...
>> while loading In[128], in expression starting on line 1 in call at 

[julia-users] Re: Better alternative to find all permutations?

[Dan]

Knuth's lexicographic permutations version is best (in terms of speed, memory. 
perhaps less readable).
The function name in the code snippet even has a little permutation by itself ;)

[julia-users] Re: Better alternative to find all permutations?

[Dan]

Well, the previous listing worked only for Int vectors. Following is a more 
properly typed version.
Various optimizations are possible: @inbounds, reducing allocations, etc.


function uniquepermutations(base)
zpos = Vector{Vector{Vector{Int}}}()
zval = Vector{eltype(base)}()
left = length(base)
for (v,c) in countmap(base)
push!(zpos,collect(subsets(collect(1:left),c)))
push!(zval,v)
left -= c
end
res = Vector{Vector{eltype(base)}}()
for a in product(zpos...)
slots = collect(1:length(base))
perm = similar(base)
for (val,b) in zip(zval,a)
perm[slots[b]] = val
slots[b] = 0
slots = filter(x->x>0,slots)
end
push!(res,perm)
end
res
end

[julia-users] Re: Better alternative to find all permutations?

[Dan]

How about this:

using StatsBase
using Iterators

function uniquepermutations(base)
zpos = Vector{Vector{Vector{Int}}}()
zval = Vector{Int}()
left = length(base)
for (v,c) in countmap(base)
push!(zpos,collect(subsets(collect(1:left),c)))
push!(zval,v)
left -= c
end
res = Vector{Vector{Int}}()
for a in product(zpos...)
slots = collect(1:length(base))
perm = zeros(length(base))
for (val,b) in zip(zval,a)
perm[slots[b]] = val
slots[b] = 0
slots = filter(x->x>0,slots)
end
push!(res,perm)
end
res
end

Familiarity with `countmap`,`subsets`,`product`,`zip` from StatsBase and 
Iterators helps. The logic is straight-forward.

[julia-users] Re: Looking for a few good algorithms

[Dan]

OK, forget the earlier stuff. Calculating A'A was a disappointment. 
Calculating AA' is the reward.
First, lets look at the function:

function AAtr(sp::SparseMatrixCSC)
res = zeros(eltype(sp),sp.m,sp.m)
for i=1:sp.n
for j=sp.colptr[i]:(sp.colptr[i+1]-1)
for k=sp.colptr[i]:(sp.colptr[i+1]-1)
res[sp.rowval[j],sp.rowval[k]] += sp.nzval[j]*sp.nzval[k]
end
end
end
res
end


Nice and compact.
The benchmarks give a substantial improvement over `full(A*A')`.
But, A'A = A'*(A')', so we also get improvement for the original problem.
The trick is doing the calculation in the right direction for the sparse 
storage, and it also uses the fact A multiplied by its own transpose. The 
exact speedup depends on the density and size. But I got 4x on 
`sprand(200,400,0.1)`

On Wednesday, November 18, 2015 at 4:41:36 AM UTC+2, Dan wrote:
>
> Embarrassingly, there is yet another bug.
>  while sp.colptr[col+1] should be
>  while sp.colptr[col+1]<=i
> (Arrg... where is the unsend/edit button)
> In case someone ever, copy-pastes this. The (currently) correct version is:
>
> function AtrA(sp::SparseMatrixCSC)
> res = zeros(eltype(sp),sp.n,sp.n)
> col = 1
> @inbounds for i=1:length(sp.nzval)
> while sp.colptr[col+1]<=i
> col+= 1
> end
> col2 = 1
> @inbounds for j in 1:length(sp.nzval)
> if sp.rowval[j] != sp.rowval[i]
> continue
> end
> while j>=sp.colptr[col2+1]
> col2+= 1
> end
> res[col,col2] += sp.nzval[i]*sp.nzval[j]
> end
> end
> res
> end
>
> On Wednesday, November 18, 2015 at 4:18:37 AM UTC+2, Dan wrote:
>>
>> It works, but it is slow (compared to `full(A'*A)`)
>>
>> On Wednesday, November 18, 2015 at 3:56:31 AM UTC+2, Dan wrote:
>>>
>>> And finally, getting rid of debug print.
>>>
>>> function AtrA(sp::SparseMatrixCSC)
>>> res = zeros(eltype(sp),sp.n,sp.n)
>>> col = 1
>>> @inbounds for i=1:length(sp.nzval)
>>> while sp.colptr[col+1]>> col+= 1
>>> end
>>> col2 = 1
>>> @inbounds for j in 1:length(sp.nzval)
>>> if sp.rowval[j] != sp.rowval[i]
>>> continue
>>> end
>>> while j>=sp.colptr[col2+1]
>>> col2+= 1
>>> end
>>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>>> end
>>> end
>>> res
>>> end
>>>
>>>
>>>
>>>
>>> On Wednesday, November 18, 2015 at 3:54:56 AM UTC+2, Dan wrote:
>>>>
>>>> Oops... a bug (which showed itself when matrix was sparser than my 
>>>> previous tests):
>>>> [when columns where completely empty, the first `if` should have been a 
>>>> `while`]
>>>>
>>>> function AtrA(sp::SparseMatrixCSC)
>>>> res = zeros(eltype(sp),sp.n,sp.n)
>>>> col = 1
>>>> @inbounds for i=1:length(sp.nzval)
>>>> while sp.colptr[col+1]>>> col+= 1
>>>> end
>>>> col2 = 1
>>>> @inbounds for j in 1:length(sp.nzval)
>>>> if sp.rowval[j] != sp.rowval[i]
>>>> continue
>>>> end
>>>> while j>=sp.colptr[col2+1]
>>>> col2+= 1
>>>> end
>>>> @show col,col2,sp.nzval[i],sp.nzval[j]
>>>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>>>> end
>>>> end
>>>> res
>>>> end
>>>>
>>>>
>>>>
>>>>
>>>> On Wednesday, November 18, 2015 at 3:35:30 AM UTC+2, Dan wrote:
>>>>>
>>>>> This is an O(nnz(A)^2) implementation. But, benchmarks will determine 
>>>>> the best solution.
>>>>>
>>>>> function AtrA(sp::SparseMatrixCSC)
>>>>> res = zeros(eltype(sp),sp.n,sp.n)
>>>>> col = 1
>>>>> @inbounds for i=1:length(sp.nzval)
>>>>> if sp.colptr[col+1]==i
>>>>> col+= 1
>>>>> end
>>>>> col2 = 1
>>>>> @inbounds for j in 1:length(sp.nzval)
>>>>> if sp.rowval[j] != sp.r

[julia-users] Re: Looking for a few good algorithms

[Dan]

Embarrassingly, there is yet another bug.
 while sp.colptr[col+1]=sp.colptr[col2+1]
col2+= 1
end
res[col,col2] += sp.nzval[i]*sp.nzval[j]
end
end
res
end

On Wednesday, November 18, 2015 at 4:18:37 AM UTC+2, Dan wrote:
>
> It works, but it is slow (compared to `full(A'*A)`)
>
> On Wednesday, November 18, 2015 at 3:56:31 AM UTC+2, Dan wrote:
>>
>> And finally, getting rid of debug print.
>>
>> function AtrA(sp::SparseMatrixCSC)
>> res = zeros(eltype(sp),sp.n,sp.n)
>> col = 1
>> @inbounds for i=1:length(sp.nzval)
>> while sp.colptr[col+1]> col+= 1
>> end
>> col2 = 1
>> @inbounds for j in 1:length(sp.nzval)
>> if sp.rowval[j] != sp.rowval[i]
>> continue
>> end
>> while j>=sp.colptr[col2+1]
>> col2+= 1
>> end
>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>> end
>> end
>> res
>> end
>>
>>
>>
>>
>> On Wednesday, November 18, 2015 at 3:54:56 AM UTC+2, Dan wrote:
>>>
>>> Oops... a bug (which showed itself when matrix was sparser than my 
>>> previous tests):
>>> [when columns where completely empty, the first `if` should have been a 
>>> `while`]
>>>
>>> function AtrA(sp::SparseMatrixCSC)
>>> res = zeros(eltype(sp),sp.n,sp.n)
>>> col = 1
>>> @inbounds for i=1:length(sp.nzval)
>>> while sp.colptr[col+1]>> col+= 1
>>> end
>>> col2 = 1
>>> @inbounds for j in 1:length(sp.nzval)
>>> if sp.rowval[j] != sp.rowval[i]
>>> continue
>>> end
>>> while j>=sp.colptr[col2+1]
>>> col2+= 1
>>> end
>>> @show col,col2,sp.nzval[i],sp.nzval[j]
>>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>>> end
>>> end
>>> res
>>> end
>>>
>>>
>>>
>>>
>>> On Wednesday, November 18, 2015 at 3:35:30 AM UTC+2, Dan wrote:
>>>>
>>>> This is an O(nnz(A)^2) implementation. But, benchmarks will determine 
>>>> the best solution.
>>>>
>>>> function AtrA(sp::SparseMatrixCSC)
>>>> res = zeros(eltype(sp),sp.n,sp.n)
>>>> col = 1
>>>> @inbounds for i=1:length(sp.nzval)
>>>> if sp.colptr[col+1]==i
>>>> col+= 1
>>>> end
>>>> col2 = 1
>>>> @inbounds for j in 1:length(sp.nzval)
>>>> if sp.rowval[j] != sp.rowval[i]
>>>> continue
>>>> end
>>>> while j>=sp.colptr[col2+1]
>>>> col2+= 1
>>>> end
>>>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>>>> end
>>>> end
>>>> res
>>>> end
>>>>
>>>>
>>>>
>>>>
>>>> On Wednesday, November 18, 2015 at 2:30:14 AM UTC+2, Dan wrote:
>>>>>
>>>>> Well, its not as simple as that because of everything being sparse 
>>>>> (adding columns is not just a loop over sum). Actually, it would be 
>>>>> interesting to see a clever solution. 
>>>>>
>>>>> On Wednesday, November 18, 2015 at 2:14:39 AM UTC+2, Dan wrote:
>>>>>>
>>>>>> the columns of A'A when A is sparse are just sparse linear 
>>>>>> combinations of columns of A. the coefficients being stored in the 
>>>>>> column-based sparse matrix representation (internally a vector of 
>>>>>> coefficients with pointers to new-column indices). so simply generating 
>>>>>> column by column with a for loop should be as-quick-as-it-gets without 
>>>>>> GPUing or some Fourier trickery.
>>>>>> did you already consider this? (i can try and make this description 
>>>>>> into code)
>>>>>>
>>>>>> On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates 
>>>>>> wrote:
>>>>>>>
>>>>>>> More careful reading of the spmatmul function in 
>>>>>>> base/sparse/linalg.jl provided a method directly downdates the dense 
>>>>>>> square 
>>>>>>> matrix C (i.e. does not form a sparse A'A then downdate) but I still 
>>>>>>> need 
>>>>>>> to form a transpose.  If anyone knows of a fast way of forming A'A 
>>>>>>> without 
>>>>>>> creating the transpose I would be happy to hear of it.
>>>>>>>
>>>>>>

[julia-users] Re: Looking for a few good algorithms

[Dan]

It works, but it is slow (compared to `full(A'*A)`)

On Wednesday, November 18, 2015 at 3:56:31 AM UTC+2, Dan wrote:
>
> And finally, getting rid of debug print.
>
> function AtrA(sp::SparseMatrixCSC)
> res = zeros(eltype(sp),sp.n,sp.n)
> col = 1
> @inbounds for i=1:length(sp.nzval)
> while sp.colptr[col+1] col+= 1
> end
> col2 = 1
> @inbounds for j in 1:length(sp.nzval)
> if sp.rowval[j] != sp.rowval[i]
> continue
> end
> while j>=sp.colptr[col2+1]
> col2+= 1
> end
> res[col,col2] += sp.nzval[i]*sp.nzval[j]
> end
> end
> res
> end
>
>
>
>
> On Wednesday, November 18, 2015 at 3:54:56 AM UTC+2, Dan wrote:
>>
>> Oops... a bug (which showed itself when matrix was sparser than my 
>> previous tests):
>> [when columns where completely empty, the first `if` should have been a 
>> `while`]
>>
>> function AtrA(sp::SparseMatrixCSC)
>> res = zeros(eltype(sp),sp.n,sp.n)
>> col = 1
>> @inbounds for i=1:length(sp.nzval)
>> while sp.colptr[col+1]> col+= 1
>> end
>> col2 = 1
>> @inbounds for j in 1:length(sp.nzval)
>> if sp.rowval[j] != sp.rowval[i]
>> continue
>> end
>> while j>=sp.colptr[col2+1]
>> col2+= 1
>> end
>>     @show col,col2,sp.nzval[i],sp.nzval[j]
>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>> end
>> end
>> res
>> end
>>
>>
>>
>>
>> On Wednesday, November 18, 2015 at 3:35:30 AM UTC+2, Dan wrote:
>>>
>>> This is an O(nnz(A)^2) implementation. But, benchmarks will determine 
>>> the best solution.
>>>
>>> function AtrA(sp::SparseMatrixCSC)
>>> res = zeros(eltype(sp),sp.n,sp.n)
>>> col = 1
>>> @inbounds for i=1:length(sp.nzval)
>>> if sp.colptr[col+1]==i
>>> col+= 1
>>> end
>>> col2 = 1
>>> @inbounds for j in 1:length(sp.nzval)
>>> if sp.rowval[j] != sp.rowval[i]
>>> continue
>>> end
>>> while j>=sp.colptr[col2+1]
>>> col2+= 1
>>>     end
>>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>>> end
>>> end
>>> res
>>> end
>>>
>>>
>>>
>>>
>>> On Wednesday, November 18, 2015 at 2:30:14 AM UTC+2, Dan wrote:
>>>>
>>>> Well, its not as simple as that because of everything being sparse 
>>>> (adding columns is not just a loop over sum). Actually, it would be 
>>>> interesting to see a clever solution. 
>>>>
>>>> On Wednesday, November 18, 2015 at 2:14:39 AM UTC+2, Dan wrote:
>>>>>
>>>>> the columns of A'A when A is sparse are just sparse linear 
>>>>> combinations of columns of A. the coefficients being stored in the 
>>>>> column-based sparse matrix representation (internally a vector of 
>>>>> coefficients with pointers to new-column indices). so simply generating 
>>>>> column by column with a for loop should be as-quick-as-it-gets without 
>>>>> GPUing or some Fourier trickery.
>>>>> did you already consider this? (i can try and make this description 
>>>>> into code)
>>>>>
>>>>> On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates 
>>>>> wrote:
>>>>>>
>>>>>> More careful reading of the spmatmul function in 
>>>>>> base/sparse/linalg.jl provided a method directly downdates the dense 
>>>>>> square 
>>>>>> matrix C (i.e. does not form a sparse A'A then downdate) but I still 
>>>>>> need 
>>>>>> to form a transpose.  If anyone knows of a fast way of forming A'A 
>>>>>> without 
>>>>>> creating the transpose I would be happy to hear of it.
>>>>>>
>>>>>

[julia-users] Re: Looking for a few good algorithms

[Dan]

And finally, getting rid of debug print.

function AtrA(sp::SparseMatrixCSC)
res = zeros(eltype(sp),sp.n,sp.n)
col = 1
@inbounds for i=1:length(sp.nzval)
while sp.colptr[col+1]=sp.colptr[col2+1]
col2+= 1
end
res[col,col2] += sp.nzval[i]*sp.nzval[j]
end
end
res
end




On Wednesday, November 18, 2015 at 3:54:56 AM UTC+2, Dan wrote:
>
> Oops... a bug (which showed itself when matrix was sparser than my 
> previous tests):
> [when columns where completely empty, the first `if` should have been a 
> `while`]
>
> function AtrA(sp::SparseMatrixCSC)
> res = zeros(eltype(sp),sp.n,sp.n)
> col = 1
> @inbounds for i=1:length(sp.nzval)
> while sp.colptr[col+1] col+= 1
> end
> col2 = 1
> @inbounds for j in 1:length(sp.nzval)
> if sp.rowval[j] != sp.rowval[i]
> continue
> end
> while j>=sp.colptr[col2+1]
> col2+= 1
> end
> @show col,col2,sp.nzval[i],sp.nzval[j]
> res[col,col2] += sp.nzval[i]*sp.nzval[j]
> end
> end
>     res
> end
>
>
>
>
> On Wednesday, November 18, 2015 at 3:35:30 AM UTC+2, Dan wrote:
>>
>> This is an O(nnz(A)^2) implementation. But, benchmarks will determine the 
>> best solution.
>>
>> function AtrA(sp::SparseMatrixCSC)
>> res = zeros(eltype(sp),sp.n,sp.n)
>> col = 1
>> @inbounds for i=1:length(sp.nzval)
>> if sp.colptr[col+1]==i
>> col+= 1
>> end
>> col2 = 1
>> @inbounds for j in 1:length(sp.nzval)
>> if sp.rowval[j] != sp.rowval[i]
>> continue
>> end
>> while j>=sp.colptr[col2+1]
>>     col2+= 1
>> end
>> res[col,col2] += sp.nzval[i]*sp.nzval[j]
>> end
>> end
>> res
>> end
>>
>>
>>
>>
>> On Wednesday, November 18, 2015 at 2:30:14 AM UTC+2, Dan wrote:
>>>
>>> Well, its not as simple as that because of everything being sparse 
>>> (adding columns is not just a loop over sum). Actually, it would be 
>>> interesting to see a clever solution. 
>>>
>>> On Wednesday, November 18, 2015 at 2:14:39 AM UTC+2, Dan wrote:
>>>>
>>>> the columns of A'A when A is sparse are just sparse linear combinations 
>>>> of columns of A. the coefficients being stored in the column-based sparse 
>>>> matrix representation (internally a vector of coefficients with pointers 
>>>> to 
>>>> new-column indices). so simply generating column by column with a for loop 
>>>> should be as-quick-as-it-gets without GPUing or some Fourier trickery.
>>>> did you already consider this? (i can try and make this description 
>>>> into code)
>>>>
>>>> On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates wrote:
>>>>>
>>>>> More careful reading of the spmatmul function in base/sparse/linalg.jl 
>>>>> provided a method directly downdates the dense square matrix C (i.e. does 
>>>>> not form a sparse A'A then downdate) but I still need to form a 
>>>>> transpose. 
>>>>>  If anyone knows of a fast way of forming A'A without creating the 
>>>>> transpose I would be happy to hear of it.
>>>>>
>>>>

[julia-users] Re: Looking for a few good algorithms

[Dan]

Oops... a bug (which showed itself when matrix was sparser than my previous 
tests):
[when columns where completely empty, the first `if` should have been a 
`while`]

function AtrA(sp::SparseMatrixCSC)
res = zeros(eltype(sp),sp.n,sp.n)
col = 1
@inbounds for i=1:length(sp.nzval)
while sp.colptr[col+1]=sp.colptr[col2+1]
col2+= 1
end
@show col,col2,sp.nzval[i],sp.nzval[j]
res[col,col2] += sp.nzval[i]*sp.nzval[j]
end
end
res
end




On Wednesday, November 18, 2015 at 3:35:30 AM UTC+2, Dan wrote:
>
> This is an O(nnz(A)^2) implementation. But, benchmarks will determine the 
> best solution.
>
> function AtrA(sp::SparseMatrixCSC)
> res = zeros(eltype(sp),sp.n,sp.n)
> col = 1
> @inbounds for i=1:length(sp.nzval)
> if sp.colptr[col+1]==i
> col+= 1
> end
> col2 = 1
> @inbounds for j in 1:length(sp.nzval)
> if sp.rowval[j] != sp.rowval[i]
> continue
> end
> while j>=sp.colptr[col2+1]
> col2+= 1
> end
> res[col,col2] += sp.nzval[i]*sp.nzval[j]
> end
> end
> res
> end
>
>
>
>
> On Wednesday, November 18, 2015 at 2:30:14 AM UTC+2, Dan wrote:
>>
>> Well, its not as simple as that because of everything being sparse 
>> (adding columns is not just a loop over sum). Actually, it would be 
>> interesting to see a clever solution. 
>>
>> On Wednesday, November 18, 2015 at 2:14:39 AM UTC+2, Dan wrote:
>>>
>>> the columns of A'A when A is sparse are just sparse linear combinations 
>>> of columns of A. the coefficients being stored in the column-based sparse 
>>> matrix representation (internally a vector of coefficients with pointers to 
>>> new-column indices). so simply generating column by column with a for loop 
>>> should be as-quick-as-it-gets without GPUing or some Fourier trickery.
>>> did you already consider this? (i can try and make this description into 
>>> code)
>>>
>>> On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates wrote:
>>>>
>>>> More careful reading of the spmatmul function in base/sparse/linalg.jl 
>>>> provided a method directly downdates the dense square matrix C (i.e. does 
>>>> not form a sparse A'A then downdate) but I still need to form a transpose. 
>>>>  If anyone knows of a fast way of forming A'A without creating the 
>>>> transpose I would be happy to hear of it.
>>>>
>>>

[julia-users] Re: Looking for a few good algorithms

[Dan]

This is an O(nnz(A)^2) implementation. But, benchmarks will determine the 
best solution.

function AtrA(sp::SparseMatrixCSC)
res = zeros(eltype(sp),sp.n,sp.n)
col = 1
@inbounds for i=1:length(sp.nzval)
if sp.colptr[col+1]==i
col+= 1
end
col2 = 1
@inbounds for j in 1:length(sp.nzval)
if sp.rowval[j] != sp.rowval[i]
continue
end
while j>=sp.colptr[col2+1]
col2+= 1
end
res[col,col2] += sp.nzval[i]*sp.nzval[j]
end
end
res
end




On Wednesday, November 18, 2015 at 2:30:14 AM UTC+2, Dan wrote:
>
> Well, its not as simple as that because of everything being sparse (adding 
> columns is not just a loop over sum). Actually, it would be interesting to 
> see a clever solution. 
>
> On Wednesday, November 18, 2015 at 2:14:39 AM UTC+2, Dan wrote:
>>
>> the columns of A'A when A is sparse are just sparse linear combinations 
>> of columns of A. the coefficients being stored in the column-based sparse 
>> matrix representation (internally a vector of coefficients with pointers to 
>> new-column indices). so simply generating column by column with a for loop 
>> should be as-quick-as-it-gets without GPUing or some Fourier trickery.
>> did you already consider this? (i can try and make this description into 
>> code)
>>
>> On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates wrote:
>>>
>>> More careful reading of the spmatmul function in base/sparse/linalg.jl 
>>> provided a method directly downdates the dense square matrix C (i.e. does 
>>> not form a sparse A'A then downdate) but I still need to form a transpose. 
>>>  If anyone knows of a fast way of forming A'A without creating the 
>>> transpose I would be happy to hear of it.
>>>
>>

[julia-users] Re: Looking for a few good algorithms

[Dan]

Well, its not as simple as that because of everything being sparse (adding 
columns is not just a loop over sum). Actually, it would be interesting to 
see a clever solution. 

On Wednesday, November 18, 2015 at 2:14:39 AM UTC+2, Dan wrote:
>
> the columns of A'A when A is sparse are just sparse linear combinations of 
> columns of A. the coefficients being stored in the column-based sparse 
> matrix representation (internally a vector of coefficients with pointers to 
> new-column indices). so simply generating column by column with a for loop 
> should be as-quick-as-it-gets without GPUing or some Fourier trickery.
> did you already consider this? (i can try and make this description into 
> code)
>
> On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates wrote:
>>
>> More careful reading of the spmatmul function in base/sparse/linalg.jl 
>> provided a method directly downdates the dense square matrix C (i.e. does 
>> not form a sparse A'A then downdate) but I still need to form a transpose. 
>>  If anyone knows of a fast way of forming A'A without creating the 
>> transpose I would be happy to hear of it.
>>
>

[julia-users] Re: Converting ASCIIString to and from Array{UInt8,1}

[Dan]

`convert` is the go-to equivalent in Julia of casting. In this case:

uint_array = convert(Vector{UInt8},str)

but there are many ways of getting similar results!
BTW `convert(ASCIIString,uint_array)` gets back.

On Wednesday, November 18, 2015 at 1:49:36 AM UTC+2, James Gilbert wrote:
>
> Solved!  After finding base/ascii.jl I find I can:
>
> str = "hello"
> uint_array = str.data
> str = ASCIIString(uint_array)
>
> though I probably shouldn't be calling str.data outside of the Julia 
> language files...
>
>

[julia-users] Re: Looking for a few good algorithms

[Dan]

the columns of A'A when A is sparse are just sparse linear combinations of 
columns of A. the coefficients being stored in the column-based sparse 
matrix representation (internally a vector of coefficients with pointers to 
new-column indices). so simply generating column by column with a for loop 
should be as-quick-as-it-gets without GPUing or some Fourier trickery.
did you already consider this? (i can try and make this description into 
code)

On Tuesday, November 17, 2015 at 10:26:21 PM UTC+2, Douglas Bates wrote:
>
> More careful reading of the spmatmul function in base/sparse/linalg.jl 
> provided a method directly downdates the dense square matrix C (i.e. does 
> not form a sparse A'A then downdate) but I still need to form a transpose. 
>  If anyone knows of a fast way of forming A'A without creating the 
> transpose I would be happy to hear of it.
>

[julia-users] Re: push! performance

[Dan]

AFAIK the implementation of `push!` in Julia only occasionally allocates a 
bigger buffer for the new vector. This results in a very low average time 
for adding multiple elements (O(1) amortized time in complexity jargon).

Notice, `a = [a new_element ]` would in Julia do something else than 
`push!` - it would create a new array (with allocation), put the old `a` 
and `new_element` in it, and then rebind the symbol `a` to the new array, 
and let the garbage collector free the old array. This would be costly, so 
it's better to use push!.

Even better, is to do the experiment and possibly report the answers on the 
list with working code, timings and versions. Hope this helps.

On Sunday, November 15, 2015 at 6:51:30 AM UTC+2, Achu wrote:
>
> So in Matlab there's a pretty big difference in performance between 
> preallocating an array, and growing it in a loop with a=[a new_element];
> How big is the performance difference in Julia between preallocating and 
> using push! to grow the array in a loop?
>
> Thanks!
> Achu
>

[julia-users] Find

[Dan]

`x .> 3` is a bit-array. The function form, which is probably what fits, is 
`x->x.>3` (the back-ticks are for quotation, not part of the code).   

[julia-users] Re: Median filter

[Dan]

As a quick stopgap, you could use the following:

medianfilter1(v,ws) = [median(v[i:(i+ws-1)]) for i=1:(length(v)-ws+1)]

medianfilter(v) = medianfilter1(vcat(0,v,0),3)


`medianfilter` does a 3 window median filter. It is not the fastest 
implementation, but it works (hopefully). 

On Tuesday, November 10, 2015 at 4:31:42 PM UTC+2, Rajn wrote:
>
> I was wondering if a 'median filter' algorithm is available in Julia? 
> There was some discussion on this topic in 2013. Wonder if someone built it 
> eventually. I am looking for something similar to medflt1 of Matlab.
> Thanks
>