[julia-users] Re: Julia computational efficiency vs C vs Java vs Python vs Cython

[Simon Kornblith]

With a 64-bit build, Julia integers are 64-bit unless otherwise specified. 
In C, you use ints, which are 32-bit. Changing them to long long makes the 
C code perform similarly to the Julia code on my system. Unfortunately, 
it's hard to operate on 32-bit integers in Julia, since + promotes to 
64-bit by default (am I missing something)?

Simon

On Tuesday, January 14, 2014 4:32:16 PM UTC-5, Przemyslaw Szufel wrote:
>
> Dear Julia users,
>
> I am considering using Julia for computational projects. 
> As a first to get a feeling of the new language a I tried to benchmark 
> Julia speed against other popular languages.
> I used an example code from the Cython tutorial: 
> http://docs.cython.org/src/tutorial/cython_tutorial.html [ the code for 
> finding n first prime numbers]. 
>
> Rewriting the code in different languages and measuring the times on my 
> Windows laptop gave me the following results:
>
> Language | Time in seconds (less=better)
>
> Python: 65.5
> Cython (with MinGW): 0.82
> Java : 0.64
> Java (with -server option) : 0.64
> C (with MinGW): 0.64
> Julia (0.2): 2.1
> Julia (0.3 nightly build): 2.1
>
> All the codes for my experiments are attached to this post (Cython i 
> Python are both being run starting from the prim.py file)
>
> The thing that worries me is that Julia takes much much longer than Cython 
> ,,,
> I am a beginner to Julia and would like to kindly ask what am I doing 
> wrong with my code. 
> I start Julia console and use the command  include ("prime.jl") to execute 
> it.
>
> This code looks very simple and I think the compiler should be able to 
> optimise it to at least the speed of Cython?
> Maybe I my code has been written in non-Julia style way and the compiler 
> has problems with it?
>
> I will be grateful for any answers or comments.
>
> Best regards,
> Przemyslaw Szufel
>

[julia-users] Re: Julia computational efficiency vs C vs Java vs Python vs Cython

[Simon Kornblith]

In C long is only guaranteed to be at least 32 bits (IIRC it's 64 bits on 
64-bit *nix but 32-bit on 64-bit Windows). long long is guaranteed to be at 
least 64 bits (and is 64 bits on all systems I know of).

Simon

On Tuesday, January 14, 2014 5:46:04 PM UTC-5, Przemyslaw Szufel wrote:
>
> Simon,
> Thanks for the explanation!
> In Java int is 32 bit as well. 
> I have just replaced ints with longs in Java and found out that now I get 
> the Java speed also very similar to Julia. 
>
> However I tried in Cython:
> def primes_list(int kmax):
> cdef int k, i
> cdef long n
> cdef long p[2]
> ...
>
> and surprisingly the speed did not change...at first I thought that maybe 
> something did not compile or is in cache - but I made sure - it's not the 
> cache. 
>  Cython speed remains unchanged regardles using int or long? 
> I know that now it becomes other language question...but maybe someone can 
> explain?
>
> All best,
> Przemyslaw Szufel
>
>
> On Tuesday, 14 January 2014 23:29:40 UTC+1, Simon Kornblith wrote:
>>
>> With a 64-bit build, Julia integers are 64-bit unless otherwise 
>> specified. In C, you use ints, which are 32-bit. Changing them to long long 
>> makes the C code perform similarly to the Julia code on my system. 
>> Unfortunately, it's hard to operate on 32-bit integers in Julia, since + 
>> promotes to 64-bit by default (am I missing something)?
>>
>> Simon
>>
>> On Tuesday, January 14, 2014 4:32:16 PM UTC-5, Przemyslaw Szufel wrote:
>>>
>>> Dear Julia users,
>>>
>>> I am considering using Julia for computational projects. 
>>> As a first to get a feeling of the new language a I tried to benchmark 
>>> Julia speed against other popular languages.
>>> I used an example code from the Cython tutorial: 
>>> http://docs.cython.org/src/tutorial/cython_tutorial.html [ the code for 
>>> finding n first prime numbers]. 
>>>
>>> Rewriting the code in different languages and measuring the times on my 
>>> Windows laptop gave me the following results:
>>>
>>> Language | Time in seconds (less=better)
>>>
>>> Python: 65.5
>>> Cython (with MinGW): 0.82
>>> Java : 0.64
>>> Java (with -server option) : 0.64
>>> C (with MinGW): 0.64
>>> Julia (0.2): 2.1
>>> Julia (0.3 nightly build): 2.1
>>>
>>> All the codes for my experiments are attached to this post (Cython i 
>>> Python are both being run starting from the prim.py file)
>>>
>>> The thing that worries me is that Julia takes much much longer than 
>>> Cython ,,,
>>> I am a beginner to Julia and would like to kindly ask what am I doing 
>>> wrong with my code. 
>>> I start Julia console and use the command  include ("prime.jl") to 
>>> execute it.
>>>
>>> This code looks very simple and I think the compiler should be able to 
>>> optimise it to at least the speed of Cython?
>>> Maybe I my code has been written in non-Julia style way and the compiler 
>>> has problems with it?
>>>
>>> I will be grateful for any answers or comments.
>>>
>>> Best regards,
>>> Przemyslaw Szufel
>>>
>>

[julia-users] Re: Julia computational efficiency vs C vs Java vs Python vs Cython

[Simon Kornblith]

Unfortunately I think this will just slow things down, since +(Int32,Int32) 
converts its arguments to Int64 and returns an Int64 (presumably to avoid 
overflow):

julia> code_typed(+, (Int32,Int32))
1-element Array{Any,1}:
 :($(Expr(:lambda, {:x,:y}, {{},{{:x,Int32,0},{:y,Int32,0}},{}}, quote  # 
int.jl, line 16:
return 
box(Int64,add_int(box($(Int64),sext_int($(Int64),x::Int32))::Int64,box($(Int64),sext_int($(Int64),y::Int32))::Int64))::Int64
end)))

which means the types of variables that start out as Int32 will be 
unstable, since you are adding to them over the course of the loop. I am 
not actually sure how to get Julia to do 32-bit addition.

Simon

On Tuesday, January 14, 2014 6:25:20 PM UTC-5, Földes László wrote:
>
> You can force the literals by enclosing them in int32():
>
> p = [int32(0) for i=1:2]
> result = [int32(0) for i=1:2]   
> k = int32(0)
> n = int32(2)
> while k < int32(2)
> i = int32(0)
>
>
>
> On Wednesday, January 15, 2014 12:04:23 AM UTC+1, Przemyslaw Szufel wrote:
>>
>> Simon,
>> Thanks!
>> I changed in Cython to 
>> def primes_list(int kmax):
>> cdef int k, i
>> cdef long long n
>> cdef long long p[2]
>> and now I am getting 2.1 seconds - exactly the same time as Julia and 
>> Java with longs...
>>
>> Since the computational difference between 64bit longs and 32bit ints is 
>> soo high - is there any way to rewrite my toy example to force Julia to do 
>> 32 bit int calculations?
>>
>> All best,
>> Przemyslaw Szufel
>>
>>
>> On Tuesday, 14 January 2014 23:55:12 UTC+1, Simon Kornblith wrote:
>>>
>>> In C long is only guaranteed to be at least 32 bits (IIRC it's 64 bits 
>>> on 64-bit *nix but 32-bit on 64-bit Windows). long long is guaranteed 
>>> to be at least 64 bits (and is 64 bits on all systems I know of).
>>>
>>> Simon
>>>
>>> On Tuesday, January 14, 2014 5:46:04 PM UTC-5, Przemyslaw Szufel wrote:
>>>>
>>>> Simon,
>>>> Thanks for the explanation!
>>>> In Java int is 32 bit as well. 
>>>> I have just replaced ints with longs in Java and found out that now I 
>>>> get the Java speed also very similar to Julia. 
>>>>
>>>> However I tried in Cython:
>>>> def primes_list(int kmax):
>>>> cdef int k, i
>>>> cdef long n
>>>> cdef long p[2]
>>>> ...
>>>>
>>>> and surprisingly the speed did not change...at first I thought that 
>>>> maybe something did not compile or is in cache - but I made sure - it's 
>>>> not 
>>>> the cache. 
>>>>  Cython speed remains unchanged regardles using int or long? 
>>>> I know that now it becomes other language question...but maybe someone 
>>>> can explain?
>>>>
>>>> All best,
>>>> Przemyslaw Szufel
>>>>
>>>>
>>>> On Tuesday, 14 January 2014 23:29:40 UTC+1, Simon Kornblith wrote:
>>>>>
>>>>> With a 64-bit build, Julia integers are 64-bit unless otherwise 
>>>>> specified. In C, you use ints, which are 32-bit. Changing them to long 
>>>>> long 
>>>>> makes the C code perform similarly to the Julia code on my system. 
>>>>> Unfortunately, it's hard to operate on 32-bit integers in Julia, since + 
>>>>> promotes to 64-bit by default (am I missing something)?
>>>>>
>>>>> Simon
>>>>>
>>>>> On Tuesday, January 14, 2014 4:32:16 PM UTC-5, Przemyslaw Szufel wrote:
>>>>>>
>>>>>> Dear Julia users,
>>>>>>
>>>>>> I am considering using Julia for computational projects. 
>>>>>> As a first to get a feeling of the new language a I tried to 
>>>>>> benchmark Julia speed against other popular languages.
>>>>>> I used an example code from the Cython tutorial: 
>>>>>> http://docs.cython.org/src/tutorial/cython_tutorial.html [ the code 
>>>>>> for finding n first prime numbers]. 
>>>>>>
>>>>>> Rewriting the code in different languages and measuring the times on 
>>>>>> my Windows laptop gave me the following results:
>>>>>>
>>>>>> Language | Time in seconds (less=better)
>>>>>>
>>>>>> Python: 65.5
>>>>>> Cython (with MinGW): 0.82
>>>>>> Java : 0.64
>>>>>> Java (with -server option) : 0.64
>>>>>> C (with MinGW): 0.64
>>>>>> Julia (0.2): 2.1
>>>>>> Julia (0.3 nightly build): 2.1
>>>>>>
>>>>>> All the codes for my experiments are attached to this post (Cython i 
>>>>>> Python are both being run starting from the prim.py file)
>>>>>>
>>>>>> The thing that worries me is that Julia takes much much longer than 
>>>>>> Cython ,,,
>>>>>> I am a beginner to Julia and would like to kindly ask what am I doing 
>>>>>> wrong with my code. 
>>>>>> I start Julia console and use the command  include ("prime.jl") to 
>>>>>> execute it.
>>>>>>
>>>>>> This code looks very simple and I think the compiler should be able 
>>>>>> to optimise it to at least the speed of Cython?
>>>>>> Maybe I my code has been written in non-Julia style way and the 
>>>>>> compiler has problems with it?
>>>>>>
>>>>>> I will be grateful for any answers or comments.
>>>>>>
>>>>>> Best regards,
>>>>>> Przemyslaw Szufel
>>>>>>
>>>>>

[julia-users] Re: DataFrames changes

[Simon Kornblith]

I believe two identical symbols are the same object, which implies that 
Dict lookup shouldn't require hashing. I haven't benchmarked this, though.

Since this is a huge change (although one that I am in favor of) that 
presumably affects a lot of existing code, any objection if I add some 
deprecation warnings?

Simon

On Wednesday, January 29, 2014 8:46:19 PM UTC-5, Cristóvão Duarte Sousa 
wrote:
>
> BTW, is there some documentation about the choice of symbols vs strings 
> for this kind of stuff (dictionary keys, optional function args, etc.)? Are 
> symbols more efficient for this?
>
>
>
> On Wednesday, January 29, 2014 5:11:20 PM UTC, John Myles White wrote:
>>
>> As we continue trying to prune DataFrames down to the essentials that we 
>> can reasonably commit to maintaining for the long-term future, we've 
>> decided to start using only symbols for the names of columns and remove all 
>> uses of strings. 
>>
>> This change will go live on master today, so please don't pull from 
>> master until you're ready to update your code. 
>>
>>  -- John 
>>
>>

Re: [julia-users] Re: DataFrames changes

[Simon Kornblith]

Done.

On Wednesday, January 29, 2014 10:03:39 PM UTC-5, John Myles White wrote:
>
> Please go ahead and add deprecation warnings.
>
>  — John
>
> On Jan 29, 2014, at 6:51 PM, Simon Kornblith 
> > 
> wrote:
>
> I believe two identical symbols are the same object, which implies that 
> Dict lookup shouldn't require hashing. I haven't benchmarked this, though.
>
> Since this is a huge change (although one that I am in favor of) that 
> presumably affects a lot of existing code, any objection if I add some 
> deprecation warnings?
>
> Simon
>
> On Wednesday, January 29, 2014 8:46:19 PM UTC-5, Cristóvão Duarte Sousa 
> wrote:
>>
>> BTW, is there some documentation about the choice of symbols vs strings 
>> for this kind of stuff (dictionary keys, optional function args, etc.)? Are 
>> symbols more efficient for this?
>>
>>
>>
>> On Wednesday, January 29, 2014 5:11:20 PM UTC, John Myles White wrote:
>>>
>>> As we continue trying to prune DataFrames down to the essentials that we 
>>> can reasonably commit to maintaining for the long-term future, we've 
>>> decided to start using only symbols for the names of columns and remove all 
>>> uses of strings. 
>>>
>>> This change will go live on master today, so please don't pull from 
>>> master until you're ready to update your code. 
>>>
>>>  -- John 
>>>
>>>
>

[julia-users] Re: problem with open(f::SubString)

[Simon Kornblith]

ccall automatically calls convert(Ptr{Uint8}, fname), which fails if the 
string is a non-terminal SubString not because the conversion is 
impossible, but because it would need to return a pointer to a copy but not 
the copy itself, there's no guarantee that copy wouldn't be 
garbage-collected before the pointer could be used.

I filed a PR  to pass 
convert(Vector{Uint8}, 
fname) to ccall, which should work for any string. Just need someone to 
confirm that an argument to ccall can't be collected until the ccall is 
complete...

Simon

On Tuesday, February 4, 2014 10:09:56 AM UTC-5, David van Leeuwen wrote:
>
> Hello, 
>
> On Saturday, February 1, 2014 8:42:07 PM UTC+1, Eric Davies wrote:
>>
>> Please forgive me if I'm giving too much obvious info; I'm not aware of 
>> your background and I want to answer this as completely as possible.
>>
>> Thanks, it is fine, if only for future reference for others who happen to 
> stumble upon something similar...
>  
>
>> This has to do with how Julia stirngs are converted to C strings. C 
>> strings are just pointers to an array of bytes where the null byte '\0' 
>> denotes the termination of a string. Julia strings are a wrapper around 
>> this structure. For a regular String, converting a Julia string to a C 
>> string simply returns the underlying pointer.
>>
>> The following is my educated guess:
>>
>> SubStrings, however, refer to the same area in memory as the String they 
>> are a section of--they just contain indicators that denote where the start 
>> and end of the substring should be.
>>
>> This is my guess, too.  So I understand where the problems comes from, 
> but I think that the general interface to open(s::String) is either too 
> general, or open() is not fully implemented. 
>
> since SubString <: String, open(s::String) claims it can handle all 
> Strings but the ccall chokes on the SubString.  It took a while before I 
> figured out methods like ascii() and utf8() convert a SubString to a 
> concrete String type.  However, I could not find a way to copy a substring 
> into a string for, e.g., the purpose of an open() call.  For instance, both 
> string(), convert(String,) and copy() did not work, the last two for 
> obvious reasons.  Now that I'm busy with it again, I think the function
>
> unsubstring{S}(s::SubString{S}) = convert(S,s)
>
>
> (if such functionality doesn't already exist) could be useful in methods 
> like open(s::String). 
>
>> Returning a pointer to the location of the SubString is problematic, as 
>> as far as C routines are concerned, it ends when the original string ends: 
>> at the null byte. Not a problem if the SubString is supposed to end at the 
>> same place, but will not work if it is supposed to end earlier than that. 
>> So, any call to a C function must create a copy of any SubString that ends 
>> earlier than the end of the original String (bytestring, ascii, utf8, copy, 
>> etc. will do this for you). 
>>
>> copy() results in a copy of a substring, which still beaks open()
>  
>
>> On Saturday, 1 February 2014 07:01:50 UTC-6, David van Leeuwen wrote:
>>>
>>> Hello, 
>>>
>>> I have problems opening a file of which the name is specified as a 
>>> SubString.  The substring is a result of a readdlm, which is fine by me, 
>>> but I don't know how to get rid of its substringness.  Well, perhaps 
>>> ascii() will do the job---but that is not the point of this message. 
>>>
>>> The problem is isolated in the following code:
>>>
>>> julia> open(match(r"ab", "abc").match, "r") ## should try to open a 
>>> file called "ab"
>>> ERROR: a SubString must coincide with the end of the original string to 
>>> be convertible to pointer 
>>>  in convert at string.jl:659 
>>>  in open at io.jl:351 
>>>  in open at io.jl:363
>>>
>>>
>>> I believe the offending line is in a ccall
>>> systemerror("opening file $fname",
>>> ccall(:ios_file, Ptr{Void},
>>>   (Ptr{Uint8}, Ptr{Uint8}, Int32, Int32, Int32, 
>>> Int32),
>>>   s.ios, fname, rd, wr, cr, tr) == C_NULL)
>>>
>>>
>>> Cheers, 
>>>
>>> ---david
>>>
>>>
>>>

Re: [julia-users] SharedArray arithmetic

[Simon Kornblith]

There are also a couple of more concise options in Base.

For multiplication by a constant, you can also use scale!(b, 2), which will 
basically do the same thing as Tim's loop for SharedArrays, but may be 
faster for arrays of BLAS types since it calls BLAS scal!.

For elementwise operations on multiple arrays, you can use broadcast!(fn, 
dest, A, B, ...), where fn is the operator, dest is the array where the 
result is to be stored (which may be A or B, and may be a SharedArray), and 
A and B are input arrays. (This works for any number of input arrays.)

I'm not sure if there's a solution for adding a scalar that doesn't involve 
a loop.

On Thursday, February 20, 2014 12:09:32 PM UTC-5, Jonathan Malmaud wrote:
>
> The devectorize package https://github.com/lindahua/Devectorize.jl can 
> deal with this:
> f(x) = @devec x[:]=2.*x
>
> 'f' won't allocate any new memory.
>
> On Wednesday, February 19, 2014 12:11:02 PM UTC-8, Tim Holy wrote:
>>
>> Unfortunately, I'm pretty sure you can't do what you're hoping for. 
>> Consider: 
>>
>> julia> A = rand(2,2) 
>> 2x2 Array{Float64,2}: 
>>  0.63259   0.109017 
>>  0.667425  0.112597 
>>
>> julia> pointer(A) 
>> Ptr{Float64} @0x0446f6a0 
>>
>> julia> A *= 2 
>> 2x2 Array{Float64,2}: 
>>  1.26518  0.218033 
>>  1.33485  0.225195 
>>
>> julia> pointer(A) 
>> Ptr{Float64} @0x05696650 
>>
>> I think A *= 2 gets translated by the parser into A = A*2. There are some 
>> old 
>> posts that can do a better job than I explaining why this is important. 
>>
>> So I think your only bet is to do something like this: 
>>
>> for i = 1:length(A) 
>> A[i] *= 2 
>> end 
>>
>> which does work in-place. 
>>
>> --Tim 
>>
>> On Wednesday, February 19, 2014 11:26:48 AM Madeleine Udell wrote: 
>> > Currently in-place arithmetic operations turn shared arrays into normal 
>> > ones: eg 
>> > 
>> > b = Base.shmem_randn(10) 
>> > b *= 2 
>> > 
>> > results in b becoming a (not-shared) array. I'd like to overload these 
>> > in-place arithmetic operations so that instead of writing 
>> > 
>> > typeof(b)==SharedArray? b.s *= 2 : b *= 2 
>> > 
>> > I can just write 
>> > 
>> > b *= 2 
>> > 
>> > and know that the shared status of my array will be preserved. How 
>> might 
>> > one overload these operations? 
>> > 
>> > (Note: I don't want c = b + b to automatically create a shared array c, 
>> > since that allocates new shared memory.) 
>>
>

[julia-users] Re: kw arguments and expression question

[Simon Kornblith]

As you have surmised, the b=5 kw argument is not actually the same as 
:(b=5)(you can see this using 
xdump(ex.args)). Instead of :(b=7), use  Expr(:kw, :b, 7). 

On Friday, February 21, 2014 6:17:48 PM UTC-5, Bassem Youssef wrote:
>
> Suppose i have a simple function as follows:
>
> julia> func(a;b=3) = a + b
> func (generic function with 1 method)
>
> as expected:
> julia> func(3,b=5)
> 8
>
> julia> ex = :(func(3,b=5))
> :(func(3,b=5))
>
> julia> eval(ex)
> 8
>
> julia> ex.args
> 3-element Array{Any,1}:
>   :func
>  3
>   :(b=5)
>
> julia> ex.args[3]
> :(b=5)
>
> However, this doesn't work:
> julia> ex.args[3] = :(b=7)
> :(b = 7)
>
> Notice, the extra spacing between b,= and 7
> julia> ex.args
> 3-element Array{Any,1}:
>   :func
>  3
>   :(b = 7)
>
> now this doesn't work:
> julia> eval(ex)
> ERROR: no method func(Int64, Int64)
>
> regards,
> Bassem
>
>  
>  

Re: [julia-users] Re: Control system library for Julia?

[Simon Kornblith]

This is great! ss2tf and potentially other functionality is also relevant 
to DSP more generally. We currently have conversions between most standard 
filter representations in DSP.jl (see 
https://github.com/JuliaDSP/DSP.jl/blob/master/src/filter_design.jl) but no 
state space representation. There is also a Polynomial package 
(https://github.com/vtjnash/Polynomial.jl) that may be a good place for 
polynomial-related functions to live.

Simon

On Friday, February 21, 2014 1:11:19 AM UTC-5, James Crist wrote:
>
> Hello all,
>
> Glad to see there's some interest in this. I'm the one behind 
> https://github.com/jcrist/Control.jl.
>
> I've been working through the getting the base types defined, before 
> starting work on the analysis functions. The first goal is to get something 
> simliar to matlab's basic control toolbox, with all the commands that an 
> undergrad would use in an intro course. After that, higher level stuff will 
> be tackled. Both Python control and Octave's control toolbox have been 
> serving as inspiration. It's surprising (not really actually) how easily 
> most of this transposes into julia.
>
> I need to get ready for a seminar I'm giving tomorrow, but over the 
> weekend I plan to commit a major refactoring of the base types 
> (TransferFunction and StateSpace) to make them more julia-friendly (python 
> has pythonic, what's the julia version?). After that, it should be fairly 
> trivial for others to write functions that work on these types.
>
> Slicot will be used to do all the heavy lifting, because it's free* and 
> why bother reinventing the wheel. I have a set of wrappers that I generated 
> for the raw interface that still need a human to look over them. I was 
> planning on doing it as I got to using individual functions, but that'd be 
> an easy thing to look through for others.
>
> - - - - -
>
> Major quesion of the moment: what plotting library is best plotting 
> library? I'm coming from heavy python usage, so winston's syntax is more 
> friendly to me.  But gadfly looks great as well. I'd rather not use pyplot 
> - I'd like to keep it as much in julia as possible. Thoughts?
>
> -Jim
>
> *Slicot (per they're website) is no longer GPL after version 4.5. However, 
> the debian repo has 5.0, and the tar ball I got contains a GPL2 license. 
> Not sure what to make of this. The most recent free version should 
> definitely be the one used.
>  
> On Thursday, February 20, 2014 9:25:10 PM UTC-6, Jeremy West wrote:
>>
>> I guess somebody got impatient with my disappearance :) I'll probably 
>> contribute to that instead, it looks like a similar roadmap I had in mind 
>> before things got messy.
>> On Feb 20, 2014 8:00 PM, "Tony Kelman"  wrote:
>>
>>> Well and at least a start on some useful functionality, to get everybody 
>>> in the same place and not duplicating the initial effort.
>>>
>>> Those kite power projects are so incredibly cool! I imagine you're using 
>>> some combination of Casadi, Acado, and/or Optimica?
>>>
>>> I do model predictive control at Berkeley, we have our own custom 
>>> Matlab/Simulink tools that work pretty well for our uses but longer-term 
>>> I'd rather have something more elegant (and in an open environment) that 
>>> doesn't have to work around Matlab's limitations and Simulink's 10+ 
>>> subtly-incompatible but still-in-common-use versions.
>>>
>>>
>>> On Thursday, February 20, 2014 4:08:54 PM UTC-8, Uwe Fechner wrote:

  Hi,

 this looks already promising. The important thing is to get started and 
 to have an issue tracker, and with this
 git repo this is already in place.

 I am currently working on automated control of kite-power systems. A 
 little video about our
 project: http://www.youtube.com/watch?v=FJmlt3_dOuA

 Best regards:

 Uwe

 Am 21.02.2014 00:24, schrieb Tony Kelman:
  
 Have a look here, https://github.com/jcrist/Control.jl is making 
 better progress than anything else I've found in the topic. He has 
 wrappers 
 to Slicot as well.

 On Thursday, February 20, 2014 1:56:20 PM UTC-8, Uwe Fechner wrote: 
>
> Hello,
>
> I could not find any control system library for Julia yet. Would that 
> make sense?
> There is a control system library available for Python:
> http://www.cds.caltech.edu/~murray/wiki/index.php/Python-control
>
> Perhaps this could be used as starting point? I think that 
> implementing this in Julia
> should be easier and faster than in Python.
>
> Any comments?
> Should I open a feature request?
>
> Uwe Fechner, TU Delft, The Netherlands
>  
  
  

Re: [julia-users] Re: Control system library for Julia?

[Simon Kornblith]

On Saturday, February 22, 2014 1:14:56 PM UTC-5, James Crist wrote:
>
> @Tony,
>
> I agree on pulling from the debian repo. Doesn't make sense to rehost 
> code. The only things that need to be changed are the build process, and 
> the alternate xerbla.f file. Perhaps the makefiles and the xerbla 
> subroutine could live in the git repo, and julia could just switch them out?
>
> @Simon,
>
> I initially started off using the polynomial package, but stopped due to 
> its limitations and some design decisions. I may try to integrate my added 
> functions and send them a pull request.
>

That would probably be useful to many people. The current lack of 
polynomial division in Polynomial.jl is pretty conspicuous.
 

> As for the DSP library, nice work! The python control package integrates 
> nicely with scipy's dsp functions, perhaps there's a way we could do the 
> same? Two thoughts on this:
>
> 1.) Conversion function between DSP types and control types
> - or -
> 2.) Same base type for both DSP TF and control TF (they're basically the 
> same)
>
> For ease of use, I'm leaning towards the first one. Does DSP.jl support 
> mimo systems?
>

DSP.jl doesn't presently support MIMO systems. I'm not sure if MIMO is 
common enough outside of control systems to be worth supporting in DSP.jl, 
so it may make sense to have separate types with conversions. But it would 
be great if we could "borrow" your conversions between state space and 
transfer function representations in DSP.jl, since there are DSP filter 
transforms that are best implemented in state space. Additionally I suspect 
that we will both eventually want to implement many of the same plots.

 
> For those interested in helping, perhaps we should move this discussion 
> somewhere else? Not sure if github issues are the best way to communicate, 
> or if a google group should be started up.
>
> -Jim
>
>
> On Friday, February 21, 2014 5:40:52 PM UTC-6, Simon Kornblith wrote:
>>
>> This is great! ss2tf and potentially other functionality is also relevant 
>> to DSP more generally. We currently have conversions between most standard 
>> filter representations in DSP.jl (see 
>> https://github.com/JuliaDSP/DSP.jl/blob/master/src/filter_design.jl) but 
>> no state space representation. There is also a Polynomial package (
>> https://github.com/vtjnash/Polynomial.jl) that may be a good place for 
>> polynomial-related functions to live.
>>
>> Simon
>>
>> On Friday, February 21, 2014 1:11:19 AM UTC-5, James Crist wrote:
>>>
>>> Hello all,
>>>
>>> Glad to see there's some interest in this. I'm the one behind 
>>> https://github.com/jcrist/Control.jl.
>>>
>>> I've been working through the getting the base types defined, before 
>>> starting work on the analysis functions. The first goal is to get something 
>>> simliar to matlab's basic control toolbox, with all the commands that an 
>>> undergrad would use in an intro course. After that, higher level stuff will 
>>> be tackled. Both Python control and Octave's control toolbox have been 
>>> serving as inspiration. It's surprising (not really actually) how easily 
>>> most of this transposes into julia.
>>>
>>> I need to get ready for a seminar I'm giving tomorrow, but over the 
>>> weekend I plan to commit a major refactoring of the base types 
>>> (TransferFunction and StateSpace) to make them more julia-friendly (python 
>>> has pythonic, what's the julia version?). After that, it should be fairly 
>>> trivial for others to write functions that work on these types.
>>>
>>> Slicot will be used to do all the heavy lifting, because it's free* and 
>>> why bother reinventing the wheel. I have a set of wrappers that I generated 
>>> for the raw interface that still need a human to look over them. I was 
>>> planning on doing it as I got to using individual functions, but that'd be 
>>> an easy thing to look through for others.
>>>
>>> - - - - -
>>>
>>> Major quesion of the moment: what plotting library is best plotting 
>>> library? I'm coming from heavy python usage, so winston's syntax is more 
>>> friendly to me.  But gadfly looks great as well. I'd rather not use pyplot 
>>> - I'd like to keep it as much in julia as possible. Thoughts?
>>>
>>> -Jim
>>>
>>> *Slicot (per they're website) is no longer GPL after version 4.5. 
>>> However, the debian repo has 5.0, and the tar ball I got contains a

Re: [julia-users] Immutable types and caching

[Simon Kornblith]

To make type inference for memoized functions suck less, all we'd need is a 
way to get a function's inferred return type for a given set of inputs in a 
way that can be used by type inference, and then we could use that to put a 
typeassert on the result. This doesn't actually seem that hard given that 
we use a similar trick for typing comprehensions. There are other potential 
uses for such functionality as well, e.g. we could call similar with the 
correct type instead of using comprehensions for vectorized operators in 
Base and then they'd work for arbitrary DenseArrays.

Simon

On Tuesday, February 25, 2014 2:17:45 PM UTC-5, Stefan Karpinski wrote:
>
> I was talking about fully automatic memoization, which does not play 
> nicely with the type system. Thus, it's a speedup in a language where every 
> function is loosely typed anyway, but a distinct slow down in a language 
> with nice type behavior (and an implementation smart enough to leverage 
> it). That's why I called automatic memoization "worse than useless" – you 
> take a possibly slow, but nicely typed function and turn it into poorly 
> typed, memoized function. There might be a way to make it work better, but 
> I suspect it would require some very deep language support and would 
> generally be inferior to just having people use an explicit cache 
> structure. In the case of npartitions, this is precisely what's happening: 
> there's an explicit, correctly typed cache.
>
>
> On Tue, Feb 25, 2014 at 1:36 PM, Kevin Squire 
> 
> > wrote:
>
>> What you're suggesting is memoization, which has come up a number of 
>>> times. The short version is that it is an easy way to speed up a slow 
>>> language but is worse than useless in fast languages. You don't see people 
>>> doing memoization in C or Fortran, do you?
>>>
>>  
>> That's a little, um, strong, don't you think?  Certainly it's not common, 
>> but probably not "worse than useless".  ;-)  (The lack of a standard 
>> map/dictionary/hash table in these languages probably contributes as well.)
>>
>> A couple of examples in Julia where it's useful (and would be just as 
>> useful in C or Fortran):
>>
>>
>> https://github.com/JuliaLang/julia/blob/master/base/combinatorics.jl#L304-L321
>>
>> https://github.com/JuliaLang/julia/blob/master/base/combinatorics.jl#L368-L379
>>
>> (Those should probably be wrapped in let statements.  I'll do that.)
>>
>> These are not common functions, of course, but memoization here is 
>> definitely useful.
>>
>> Cheers!
>>Kevin
>>
>
>

[julia-users] Re: trouble running julia

[Simon Kornblith]

Looks like https://github.com/JuliaLang/julia/issues/5750. Deleting 
sys.dylib as described there should fix this, but it will make Julia 
startup slower. If you compile yourself, it should fix this and startup 
should be fast. I don't think there's a major reason to run the prerelease 
binaries instead of compiling master (with or without Homebrew) except that 
compiling takes time and requires dependencies.

On Tuesday, February 25, 2014 5:48:32 PM UTC-5, Jon Norberg wrote:
>
> I have been running julia for a while lately with ijulia and everything 
> worked nicely. then I started building julia with brew --HEAD to keep up to 
> date with the latest changes, worked well enough. But now something 
> happened and nothing works, test run won't work even though julia installs 
> without errorseven when I download the prerelease the process just 
> shuts down
>
> exec 
> '/Applications/Julia-0.3.0-prerelease-a673e4c4de.app/Contents/Resources/julia/bin/julia'
>
>
> [Process completed]
>
> So I have a couple of questions:
>
> 1) whats likely wrong here? is it some dependencies? how do I make a 
> REALLY clean install
> 2) is it worth brewing --HEAD julia? How often does/should it fail?
> 3) when sticking to a stable release I notice that the number of warnings 
> increases as time goes towards the next release. How to avoid? 
> 4) Whats the best strategy to have a healthy and reasonably updated julia
>
> Many thanks (again)
>

[julia-users] ANN: VML.jl

[Simon Kornblith]

I created a package that makes Julia use Intel's Vector Math Library for 
operations on arrays of Float32/Float64. VML provides vectorized versions 
of many of the functions in openlibm with equivalent (<1 ulp) accuracy (in 
VML_HA mode). The speedup is often quite stunning; see benchmarks at 
https://github.com/simonster/VML.jl.

Simon

[julia-users] Re: ANN: VML.jl

[Simon Kornblith]

Some of the poorest performers here are trunc, ceil, floor, and round (as 
of LLVM 3.4). We currently call out to libm for these, but there are LLVM 
intrinsics that are optimized into a single instruction. It looks like the 
loop vectorizer may even vectorize these intrinsics automatically.

Simon

On Thursday, February 27, 2014 9:04:18 PM UTC-5, Simon Kornblith wrote:
>
> I created a package that makes Julia use Intel's Vector Math Library for 
> operations on arrays of Float32/Float64. VML provides vectorized versions 
> of many of the functions in openlibm with equivalent (<1 ulp) accuracy (in 
> VML_HA mode). The speedup is often quite stunning; see benchmarks at 
> https://github.com/simonster/VML.jl.
>
> Simon
>

Re: [julia-users] Re: ANN: VML.jl

[Simon Kornblith]

For the arithmetic functions, the mutating versions are probably not worth 
it. For add!, subtract!, and multiply! the VML implementations on Float64 
are not any faster, and with the loop vectorizer I am hoping we can erase 
the remainder of the gap for Float32. In any case these operations are not 
expensive compared to load/store and having them might actually encourage 
people to write code that performs worse. Similarly I think that ceil, 
floor, round, and trunc should be quite fast once we switch to using LLVM 
intrinsics instead of openlibm. I am not sure how VML beats the naive 
implementation for division of Float64.

For the transcendental functions, the mutating versions may make more 
sense, since the performance difference will be much larger and I do not 
think the loop vectorizer will transform a loop into a call to the 
vectorized function.

Simon

On Friday, February 28, 2014 10:14:01 AM UTC-5, John Myles White wrote:
>
> If you get them all to export the same API, you could, in principle, just 
> switch `using VML` to `using Yeppp`.
>
> My question: are we finally conceding that add! and co. is probably worth 
> having?
>
>  — John
>
> On Feb 28, 2014, at 7:10 AM, Dahua Lin > 
> wrote:
>
> This is very nice.
>
> Now that we have several "back-ends" for vectorized computation, VML, 
> Yeppp, Julia's builtin functions, as well as the @simd-ized versions, I am 
> considering whether there is a way to switch back-end without affecting the 
> client codes.
>
> - Dahua
>
>
> On Thursday, February 27, 2014 11:58:20 PM UTC-6, Stefan Karpinski wrote:
>>
>> We really need to stop using libm for those.
>>
>>
>> On Fri, Feb 28, 2014 at 12:40 AM, Simon Kornblith wrote:
>>
>>> Some of the poorest performers here are trunc, ceil, floor, and round 
>>> (as of LLVM 3.4). We currently call out to libm for these, but there are 
>>> LLVM intrinsics that are optimized into a single instruction. It looks like 
>>> the loop vectorizer may even vectorize these intrinsics automatically.
>>>
>>> Simon
>>>
>>>
>>> On Thursday, February 27, 2014 9:04:18 PM UTC-5, Simon Kornblith wrote:
>>>>
>>>> I created a package that makes Julia use Intel's Vector Math Library 
>>>> for operations on arrays of Float32/Float64. VML provides vectorized 
>>>> versions of many of the functions in openlibm with equivalent (<1 ulp) 
>>>> accuracy (in VML_HA mode). The speedup is often quite stunning; see 
>>>> benchmarks at https://github.com/simonster/VML.jl.
>>>>
>>>> Simon
>>>>
>>>
>>
>

Re: [julia-users] Re: ANN: VML.jl

[Simon Kornblith]

I have benchmarked with short vectors, and for allocating functions, the 
overhead of allocation dominates everything else, so it doesn't really 
matter what implementation is used. I haven't tested mutating functions, 
but we don't presently have those in Base. I think it would be fine to tell 
people that, if their vectors are very short and they want maximal 
performance for mutating transcendentals, then they should use a loop.

I'm not sure why we shouldn't use the VML functions everywhere if they are 
available. In VML_HA mode they are equally accurate to openlibm. The only 
weird part is that sometimes it is possible that e.g. sin([x])[1] != 
sin(x), but I don't think there's much code out there that relies on that. 
With Yeppp there may be additional concerns since the functions are not as 
accurate as openlibm.

Simon

On Friday, February 28, 2014 12:30:21 PM UTC-5, Dahua Lin wrote:
>
> While VML is generally much faster for big arrays, the overhead is 
> considerable or even dominant for small ones. 
>
> I think the ideal way is to let people to explicit articulate their 
> intention to use these instead of changing the behavior silently.  For 
> example, I may want to be able to use these over several large arrays, but 
> not at each place with ``a + b``.
>
> Think about it more, I think the best way to do this is for each package 
> to define their own functions (not extending the functions in the base) 
> with consistent naming. Then in the client code, we will be able to do
> VML.exp(a1)
> IPP.exp(a2)
> Yeppp.exp(a3)
>
> I think what we only have to do is to coordinate the naming across several 
> packages. So I can easily replace ``Yeppp.add`` with ``VML.add`` or the 
> other way round.
>
> - Dahua
>
>
>
> On Friday, February 28, 2014 9:14:01 AM UTC-6, John Myles White wrote:
>>
>> If you get them all to export the same API, you could, in principle, just 
>> switch `using VML` to `using Yeppp`.
>>
>> My question: are we finally conceding that add! and co. is probably worth 
>> having?
>>
>>  — John
>>
>> On Feb 28, 2014, at 7:10 AM, Dahua Lin  wrote:
>>
>> This is very nice.
>>
>> Now that we have several "back-ends" for vectorized computation, VML, 
>> Yeppp, Julia's builtin functions, as well as the @simd-ized versions, I am 
>> considering whether there is a way to switch back-end without affecting the 
>> client codes.
>>
>> - Dahua
>>
>>
>> On Thursday, February 27, 2014 11:58:20 PM UTC-6, Stefan Karpinski wrote:
>>>
>>> We really need to stop using libm for those.
>>>
>>>
>>> On Fri, Feb 28, 2014 at 12:40 AM, Simon Kornblith 
>>> wrote:
>>>
>>>> Some of the poorest performers here are trunc, ceil, floor, and round 
>>>> (as of LLVM 3.4). We currently call out to libm for these, but there are 
>>>> LLVM intrinsics that are optimized into a single instruction. It looks 
>>>> like 
>>>> the loop vectorizer may even vectorize these intrinsics automatically.
>>>>
>>>> Simon
>>>>
>>>>
>>>> On Thursday, February 27, 2014 9:04:18 PM UTC-5, Simon Kornblith wrote:
>>>>>
>>>>> I created a package that makes Julia use Intel's Vector Math Library 
>>>>> for operations on arrays of Float32/Float64. VML provides vectorized 
>>>>> versions of many of the functions in openlibm with equivalent (<1 ulp) 
>>>>> accuracy (in VML_HA mode). The speedup is often quite stunning; see 
>>>>> benchmarks at https://github.com/simonster/VML.jl.
>>>>>
>>>>> Simon
>>>>>
>>>>
>>>
>>

[julia-users] Re: Ask for help in GLMNet Package

[Simon Kornblith]

It seems like you don't have a GLMNet binary. This was supposed to be built 
when you ran Pkg.add("GLMNet"); maybe you got an error message? Make sure 
you have gfortran installed and then run Pkg.build("GLMNet").

Simon

On Sunday, March 2, 2014 10:23:40 PM UTC-5, Bina Caroline wrote:
>
>
>
> I run the test file in GLMNet package, but it gives me errors like this. 
> Could someone help me figure out the problem? 
> Also, the GLMNet package was listed on the available package website 
> before but it disappear today. So confused.
>
> *julia> **path = glmnet(X, y)*
>
> *ERROR: error compiling glmnet: error compiling __glmnet#131__: error 
> compiling glmnet: error compiling __glmnet#132__: error compiling glmnet: 
> error compiling __glmnet#130__: error compiling glmnet!: error compiling 
> __glmnet!#126__: could not load module 
> /Users/yuchen/.julia/GLMNet/deps/libglmnet.so: 
> dlopen(/Users/yuchen/.julia/GLMNet/deps/libglmnet.so.dylib, 1): image not 
> found*
>
>
>
>

[julia-users] Re: String to Float64 performance question

[Simon Kornblith]

I've only taken a brief look, but one possibility is garbage collection 
overhead. If you have a lot of strings living in memory, each collection 
may take much longer. There are plans to make this perform better, but you 
are probably best off avoiding large arrays of strings to begin with. If 
you are reading from a text file, you can convert strings to floats as you 
read instead of reading a bunch of strings and then converting them all to 
floats.

Simon

On Monday, March 3, 2014 4:43:43 PM UTC-5, Keith Campbell wrote:
>
> Hi all, 
>
> While working to bring down the timing of a text IO oriented job, I seem 
> to be getting string-to-float conversion times that grow more than linearly 
> with respect to the length of the string array I'm converting.
>
> I'd appreciate any help you could provide in speeding up the Julia 
> version.  
>
> The context here is a port from Python.  Julia's times are quite 
> competitive for an array of 1000,000 strings, but for 15,000,000 strings 
> the Python job is substantially faster.  The Julia version is much faster 
> in several other parts of the code, so it would be great to nail down this 
> last detail.
>
> Essentially I'm just doing [float(v) for v in fstrings]
>
> More coding details, along with a table and plot of run time vs array 
> length are posted as an IJulia notebook that can be viewed at:
>
>
> http://nbviewer.ipython.org/gist/catawbasam
> /9335068
>
> I'm currently running: 
>  Version 0.3.0-prerelease (2014-02-28 04:44 UTC)
>   master/4e48d5b* (fork: -1 commits, 3 days) 
>  x86_64-redhat-linux
>
> thanks,
> Keith
>
>

[julia-users] Re: Immutable objects and identity

[Simon Kornblith]

If you have e.g. a function that always returns Int, then you can specify 
the associative type as Dict{Any,Int}, which allows type inference to 
determine that anything pulled out of it is an Int, so the return type of 
the memoized function can be inferred:

julia> using Memoize

julia> @memoize Dict{Any,Int} f(y) = y + 1

julia> Base.return_types(f, (Int,))
1-element Array{Any,1}:
 Int64

Of course this doesn't work if you don't specify the Dict type:

julia> @memoize f(y) = y + 1

julia> Base.return_types(f, (Int,))
1-element Array{Any,1}:
 Any

Simon

On Wednesday, March 5, 2014 2:53:09 PM UTC-5, andrew cooke wrote:
>
> we're talking about  
> https://github.com/simonster/Memoize.jl/blob/master/src/Memoize.jl ?
>
> i must be missing something.  where is the type?
>
> thanks,
> andrew
>
> On Wednesday, 5 March 2014 16:24:34 UTC-3, Pierre-Yves Gérardy wrote:
>>
>> It already accepts a type as an optional parameter.
>>
>> On Wednesday, March 5, 2014 7:41:42 PM UTC+1, andrew cooke wrote:
>>>
>>>
>>> it seemed to me (when i looked at the code) that would be trivial to add 
>>> to Memoize.jl by adding a type parameter to the macro and reading from the 
>>> cache into a typed variable.
>>>
>>> but since i spilt coke all over my netbook yesterday my julia 
>>> programming is on hiatus :o(
>>>
>>> andrew
>>>
>>>
>>> On Wednesday, 5 March 2014 10:26:37 UTC-3, Toivo Henningsson wrote:

 The trouble with memoization and type inference is that when type 
 inference sees a memoized function, it will see that the return value is 
 either computed, or looked up in the cache. The type of the cached value 
 might be anything, unless the user supplies an explicit type restriction. 
 There's no way to tell type inference that the cached must have been 
 computed by the function previously, so that it can be disregarded when it 
 comes to inferring the return type.

 On Wednesday, 5 March 2014 01:34:10 UTC+1, Pierre-Yves Gérardy wrote:
>
>
>
> On Tuesday, March 4, 2014 3:01:13 AM UTC+1, andrew cooke wrote:
>>
>>
>> i assumed that passing by value was an optimization, but i don't know 
>> what triggers it.  it may be that it's only for immutable types without 
>> pointers to mutable types.  or it may be based on size.  it might be 
>> something to worry about if your instances are large.
>>
>> there's a memoize.jl module somewhere that defines a macro that 
>> automatically memoizes the function it's applied to.  so you could use 
>> that 
>> on your constructor function, but:
>>
>>  1 - i'm not sure how nesting affects things
>>
>  2 - when i looked at memoize.jl, by default it used a simple map (no 
>> weak keys or lru eviction)
>>
>
> Yes, but you can pass a collection of your choice as a cache. Julia 
> provides a WeakKeyDict, but if your function depends on multiple inputs, 
> using tuples as keys makes it moot.
>
> I'll go ahead and skip this whole identity thing for now :-)
>
>  3 - someone said that using memoize.jl broke type inference of return 
>> arguments (i don't understand why) which had a performance impact.
>>
>
> That was Stephan Karpinski. I wonder if it holds when using a 
> `memoize` macro that creates a named function and a typed cache.
>  
> —Pierre-Yves
>

[julia-users] Re: Passing generic functions is slow? (Kernel Density Estimation)

[Simon Kornblith]

See "method lookup hosting" and "function-valued argument inlining" in 
https://github.com/JuliaLang/julia/issues/3440. One good workaround is to 
pass a type instead of a function. If you define:

type X; end
evaluate(::X) = ...

then you can pass X() to a function and evaluate(t) and will be as fast as 
if you defined a function in the global scope.

In this particular case, you might consider just using the kde function 
from StatsBase, which uses an FFT-based algorithm that should be much 
faster.

Simon

On Wednesday, March 5, 2014 8:57:09 PM UTC-5, Antonio Linero wrote:
>
> Below is a (highly non-optimal) implementation of kernel density 
> estimation. In principle, I'd like to pass in different kernel functions 
> (e.g. use the same code with a Gaussian kernel or triweight kernel). 
> However, I get a huge performance hit if I try to pass the kernel function 
> in. The following code calculates, for the Gaussian kernel, the KDE on a 
> grid of size 1000 given 1000 training data points. The first implementation 
> gets the kernel function from the global environment while the second takes 
> the kernel function passed in as an argument.
>
> Why does this give different performance and is there any way to repair 
> it? The version that does not pass the kernel function is (1) at least 
> three times faster and (2) allocates far more memory.
>
> using Rmath
>
>
> function kern(x::Float64, y::Float64, h::Float64)
> return Rmath.dnorm(x, y, h)
> end
>
>
> function eval_kern(x::Array{Float64, 1}, X::Array{Float64, 1}, h::Float64)
> N   = length(x)
> K   = length(X)
> out = zeros(N)
> for i in 1:N
> for k in 1:K
> out[i] += kern(x[i], X[k], h)
> end
> end
> return out / N
> end
>
>
> function eval_kern(x::Array{Float64, 1}, X::Array{Float64, 1}, h::Float64, 
> Kern::Function)
> N   = length(x)
> K   = length(X)
> out = zeros(N)
> for i in 1:N
> for k in 1:K
> out[i] += Kern(x[i], X[k], h)
> end
> end
> return out / N
> end
>
> N = 1000
> X = randn(N)
> h = std(X) * N ^ (-.2) * 1.06
> x_gridd = linspace(-4.0, 4.0, 1000);
>
> @time y_gridd = eval_kern(x_gridd, X, h);
> @time y_gridd = eval_kern(x_gridd, X, h, kern);
>
> This gives the following results for me:
>
> julia> @time y_gridd = eval_kern(x_gridd, X, h);
> elapsed time: 0.062120766 seconds (422240 bytes allocated)
>
> julia> @time y_gridd = eval_kern(x_gridd, X, h, kern);
> elapsed time: 0.249457403 seconds (104151016 bytes allocated)
>
>
>
> I've experimented with a few things. Using Rmath.dnorm directly rather 
> than kern cuts down on the memory allocation for the first implementation 
> but doesn't seem to change the run time.. If I use
> out[i] += Kern(x[i]::Float64, X[k]::Float64, h::Float64)::Float64
> instead I get within a factor of 2 in terms of speed and a 50% decreasing 
> in memory allocation. 
>

[julia-users] Re: Julia optimization for mapping data classifications

[Simon Kornblith]

Your algorithm looks fine. The problems are entirely in your testing 
script. The first issue is that JSON.parse returns a Vector{Any}, which 
deoptimizes everything. Try:

f = open("test.json")
data = convert(Vector{Float64}, JSON.parse(f))

The second issue is that you're including compilation time in your test. 
Try:

jenks(data, 5)
@time jenks(data, 5)

With these changes, I get:

➜  jenks.jl git:(master) ✗ julia jenks-run.jl 
elapsed time: 0.013948242 seconds (199084 bytes allocated)

On Tuesday, March 25, 2014 8:45:43 PM UTC-4, Shaun Walbridge wrote:
>
> Hello,
>
> I've been playing around with Julia for some data classifiers commonly 
> used in mapping, such as the Jenks "natural breaks" 
> algorithm. 
> (For background and a Javascript implementation, I highly recommend Tom 
> MacWright's literate programming 
> implementation.) 
>  Recently, a friend created a Cython Python 
> version of 
> Jenks which is very performant. I've created my Julia version largely based 
> on his version, and based on posts like this one on isotonic regression 
> in Julia , I had hoped the 
> Julia one would be in a similar ballpark in terms of performance to Cython. 
> Here's my very basic timing results:
>
> Matt Perry's Cython version:
>   In [15]: %timeit jenks(data, 5)
>   100 loops, best of 3: 13.9 ms per loop
>
> My Julia version, running in 0.2.1 (running against master produced slower 
> results):
> julia> @time jenks(data, 5)
> elapsed time: 1.046356641 seconds (646397168 bytes allocated)
>
> In comparison, an implementation in Python which only uses 
> lists runs 
> in about 2.8 seconds on my machine. So I imagine that I must be doing 
> something wrong, because I imagine the performance different should not be 
> in favor of the Cython version by a factor of 75, and should handily 
> dispatch the implementation which uses ill-fitting data structures. My 
> Julia code is a rather crude port and is by no means idiomatic, I did do 
> some basic profiling, and most of the runtime seems to come from the basic 
> math performed in each loop, I've only done some minor performance 
> optimization. Any and all advice appreciated on what would make this code 
> more performant for this particular task. 
>
> Code and data at:
>   https://github.com/scw/jenks.jl
>
> Thanks for your help,
> Shaun
>

Re: [julia-users] Cubic spline interpolation

[Simon Kornblith]

You could write a C interface to the C++ spline library and then ccall that 
from Julia. That's probably not too much work, but I can't promise that the 
interface plus the port to Julia would be less work than fixing the 
segfaults.

Simon

On Wednesday, March 26, 2014 10:30:52 AM UTC-4, Tomas Lycken wrote:
>
> I'd happily use quadratic interpolation if it weren't for that exact 
> limitation - I *do* need continuous second derivatives in the 2D 
> interpolation...
>
> I've never actually implemented spline interpolation myself, so I'd 
> probably have to spend quite a lot of time delving into the current code to 
> understand how it works and how to add to it. Might be good exercise, but 
> unfortunately not something I have time for at the moment. It's an 
> interesting learning project, though, so I'll definitely consider coming 
> back to it later =)
>
> The reason I started this thread was that I'm currently working on a quite 
> large project (my masters thesis project) in C++, but after realizing today 
> that my current codebase has some deeply nested memory management issues 
> that I'm having difficulties troubleshooting, I wanted to see if it would 
> be possible to port the project to Julia in a relatively short amount of 
> time. For that purpose, I wanted to see if all the stuff I'm doing through 
> external C++ libraries was implemented in Julia packages as well, so that I 
> wouldn't have to write any Julia code for stuff that I haven't had to write 
> C++ code for. Cubic splines was the only missing piece of the puzzle :P 
> However, implementing a cubic spline interpolation routine is, 
> unfortunately, well out of scope for my thesis, so it will have to wait 
> until I have time to spend on it. Until then, I'd better get back to those 
> segfaults...
>
> // T
>
> On Wednesday, March 26, 2014 2:54:55 PM UTC+1, Tim Holy wrote:
>>
>> I should also add that unless you need a continuous second derivative, 
>> quadratic (while not as popular) is IMO nicer than cubic in several 
>> ways---for 
>> example, by requiring fewer points (27 rather than 64 in 3 dimensions), 
>> it's 
>> faster to evaluate. I suspect the main reason quadratic isn't popular is 
>> that 
>> it's "non-interpolating" (in the terminology of P. Thevenaz, T. Blu, and 
>> M. 
>> Unser, 2000), but the infrastructure in Grid (based on interp_invert!) 
>> compensates for that. 
>>
>> --Tim 
>>
>> On Wednesday, March 26, 2014 06:26:12 AM Tomas Lycken wrote: 
>> > Hi, 
>> > 
>> > Is there a (maintained) package somewhere with cubic spline 
>> capabilities? I 
>> > need something that fulfills the following requirements: 
>> > 
>> > * Scalar-valued functions of one variable, f(x), specified on uniform 
>> or 
>> > non-uniform x-grids 
>> > * Scalar-valued functions of two variables, f(x,y), at least specified 
>> on 
>> > uniform grids that don't need to have the same spacing in x and y (i.e. 
>> > rectangular, but not necessarily quadratic, grid cells) 
>> > * Evaluation of function value and evaluate up to at least second order 
>> > derivatives, i.e. both f'x, f'y, f'xx, f'xy and f'yy in the 2D case 
>> > 
>> > The only packages I've found that seem to approach this functionality 
>> are 
>> > 
>> > * https://github.com/timholy/Grid.jl - only up to quadratic splines, 
>> as far 
>> > as I can tell from the readme; also unsure on if it can evaluate second 
>> > order derivatives 
>> > * https://github.com/gusl/BSplines.jl - only 1D interpolation, and 
>> base 
>> > splines rather than exact cubic splines 
>> > * https://github.com/EconForge/splines - Specifies Julia 0.2- in its 
>> > require and hasn't been touched in four months => I doubt it works with 
>> my 
>> > Julia installation which is at master. It would also probably take 
>> quite a 
>> > lot of work to learn to use it, since it has no documentation at all. 
>> > 
>> > Are there any others that I've missed? Is there any non-official effort 
>> > toward creating this functionality in Julia? 
>> > 
>> > Thanks in advance, 
>> > 
>> > // Tomas 
>>
>

[julia-users] least squares residuals from QR ldiv?

[Simon Kornblith]

Is there a way to get the least squares residuals for qrfact(X)\y without 
having to compute X*beta, as is given by gels? I believe this should be 
||Q2'y||; Q2'y appears to be computed but explicitly zeroed in A_ldiv_B!.

Simon

[julia-users] Re: least squares residuals from QR ldiv?

[Simon Kornblith]

(To clarify, I mean the sum of the squared residuals, and there's a ^2 
missing below.)

On Monday, April 7, 2014 5:48:49 PM UTC-4, Simon Kornblith wrote:
>
> Is there a way to get the least squares residuals for qrfact(X)\y without 
> having to compute X*beta, as is given by gels? I believe this should be 
> ||Q2'y||; Q2'y appears to be computed but explicitly zeroed in A_ldiv_B!.
>
> Simon
>

Re: [julia-users] Re: multi-taper F-test and Slepian sequences

[Simon Kornblith]

The dpss function now lives in DSP.jl (https://github.com/JuliaDSP/DSP.jl)

Simon

On Wednesday, April 9, 2014 7:47:38 PM UTC-4, Kevin Squire wrote:
>
> Hi Simon,
>
> I noticed that you removed windows.jl from FrequencyDomainAnalysis.jl (now 
> Synchrony.jl).  Are they by chance available elsewhere?
>
> Cheers,
>Kevin
>
>
> On Fri, Sep 20, 2013 at 8:05 AM, Simon Kornblith 
> 
> > wrote:
>
>> I have some code for this here:
>>
>>
>> https://github.com/simonster/FrequencyDomainAnalysis.jl/blob/master/src/windows.jl
>>
>> There's also code for multitaper PSD/coherence estimates and a few other 
>> useful things in that repository. I will get around to 
>> finishing/documenting/releasing this soon...
>>
>>
>> On Friday, September 20, 2013 8:59:21 AM UTC-4, Ben Arthur wrote:
>>
>>> folks,
>>>
>>> i'm new to julia and very impressed with it's speed.
>>>
>>> i have a signal processing project involving thomson's multi-taper 
>>> F-test that i'm about to port from matlab.  before doing so i just wanted 
>>> to confirm that no julia code existed yet to compute slepian's discrete 
>>> prolate spheroidal sequences.  i assume not since i don't even see hamming 
>>> and hann window functions in the standard library, but please let me know 
>>> if i'm wrong!
>>>
>>> thanks,
>>>
>>> ben arthur
>>> HHMI janelia farm
>>>
>>
>

Re: [julia-users] Re: Non-GPL Julia?

[Simon Kornblith]

I'm not a lawyer either, but I don't think this is any more problematic 
than the current situation with inclusion of Rmath etc. The combined work 
is GPL (or LGPL, once we get rid of the GPL parts), but the vast majority 
of the source files are (also) MIT. This seems fine since MIT gives all of 
the permissions of (L)GPL and then some, and just because an author 
licenses something as (L)GPL doesn't mean that author can't also make that 
file available under another license.

On Thursday, April 10, 2014 1:31:24 PM UTC-4, Mike Nolta wrote:
>
> On Thu, Apr 10, 2014 at 11:07 AM, Tony Kelman > 
> wrote: 
> > MPFR and GMP are LGPL, which is not quite as problematic or viral. 
> > 
> > Some of SuiteSparse is GPL, parts of it are LGPL, and at least one file 
> of 
> > the Julia code in base for sparse matrices that is based on parts of 
> > SuiteSparse is also LGPL. 
> > 
>
> IANAL, but doesn't this mean we've crossed the line from "work that 
> uses the library" to "derivative work", and thus are in violation of 
> the LGPL? Specifically, this part of section (2): 
>
> "These requirements apply to the modified work as a whole. If 
> identifiable sections of that work are not derived from the Library, 
> and can be reasonably considered independent and separate works in 
> themselves, then this License, and its terms, do not apply to those 
> sections when you distribute them as separate works. But when you 
> distribute the same sections as part of a whole which is a work based 
> on the Library, the distribution of the whole must be on the terms of 
> this License, whose permissions for other licensees extend to the 
> entire whole, and thus to each and every part regardless of who wrote 
> it." 
>
> -Mike 
>
> > 
> > On Thursday, April 10, 2014 7:18:45 AM UTC-7, Jake Bolewski wrote: 
> >> 
> >> Are there any more liberally licensed libraries that get close to the 
> >> functionality in MPFR, I know that there are some non-GPL BigInt 
> >> implementations out there although I don't think any match GMP's 
> >> performance. 
> >> 
> >> On Thursday, April 10, 2014 10:13:20 AM UTC-4, Stefan Karpinski wrote: 
> >>> 
> >>> There's been a lot of talk about turning various bits of functionality 
> >>> like GMP, MPFR, FFTW and such into packages that simply happen to be 
> >>> pre-loaded by default, making it easy to get a much more spare basic 
> Julia 
> >>> version. This will definitely happen over the summer. Note that 
> OpenBLAS is 
> >>> BSD-license, so using MKL is not necessary for non-GPL Julia. 
> >>> 
> >>> 
> >>> On Thu, Apr 10, 2014 at 10:07 AM, Isaiah Norton  
> >>> wrote: 
>  
>  Rmath too, but I think that is on the way out. All of the licenses 
> are 
>  linked here: 
>  https://github.com/JuliaLang/julia/blob/master/LICENSE.md 
>  
>  
>  On Thu, Apr 10, 2014 at 10:05 AM, Jake Bolewski  
>
>  wrote: 
> > 
> > As readline is now removed I think GMP (BigInt's) and MPFR 
> (arbitrary 
> > precision floating points) are the only GNU ibraries left if you are 
> able to 
> > use MKL and don't require FFTW.  Steven also working on a branch 
> where he 
> > provides FFT support in pure julia. 
> > 
> > Best, 
> > Jake 
> > 
> > 
> > On Thursday, April 10, 2014 9:26:55 AM UTC-4, Jay Kickliter wrote: 
> >> 
> >> There are bits and pieces in Github issues and posts, but can post 
> a 
> >> definitive list of what needs to be replaced/removed to make Julia 
> non GPL? 
> >> Will any functionality be missing? From what I understand I can use 
> MKL for 
> >> some stuff. I've read that MKL has the ability to mimic FFTW, but 
> will Julia 
> >> use that interface? 
> >> 
> >> For the record I'm not anti-GPL. I'd like to pitch Julia to my 
> company 
> >> as alternative to Matlab and C++. But our customers can't accept a 
> project 
> >> built with GPL. It's not a problem now, but I'm looking down the 
> road when 
> >> Julia can be compiled in to executables. 
>  
>  
> >>> 
> > 
>

Re: [julia-users] Re: Export/view movie of plots

[Simon Kornblith]

Assuming avconv or ffmpeg is available on your system, you can open a pipe 
to it:

pipe, process = writesto(`avconv -y -f rawvideo -pix_fmt gray -s 100x100
 -r 30 -i - -an -c:v libx264 -pix_fmt yuv420p movie.mp4`)

The options are detailed in the docs; -s is the movie size and -r is the 
frame rate.

(If you have ffmpeg, s/avconv/ffmpeg/ and I think this should still work.) 
Now every 100x100 block of Uint8s that you write to pipe creates a frame in 
the movie (you can obviously change the resolution). Note that the first 
dimension of the array is the x dimension and the second is the y 
dimension, so you may want to transpose the array before writing it. For 
example, to create a movie with 1000 frames of static on the top and black 
on the bottom:

for i = 1:1000
write(pipe, [rand(Uint8, 100, 50) zeros(Uint8, 100, 50)])
end
close(pipe)

If you have a 3D array of Uint8s where the dimensions are ordered as x, y, 
and time, then you can also write that directly. If you write something 
besides Uint8s, you will get garbage. You can tweak the avconv/ffmpeg 
encoding settings as necessary and also change the pixel format to add 
color if you want.

Hope this helps,
Simon

On Friday, April 11, 2014 8:45:43 AM UTC-4, Sheehan Olver wrote:
>
> Was hoping for something one line on Julia
>
> Sent from my iPhone
>
> On Apr 11, 2014, at 10:18 PM, Arnim Bleier > 
> wrote:
>
> Hi,
>
> create pictures for each frame [1:N] ... (lets say pic_n.png)
>
> then
>
> mencoder mf://*.png -mf w=200:h=150:fps=25:type=png -ovc copy -oac copy -o 
> out.avi
>
> now you should have a "out.avi" 
>
>
> Best
>
> Arnim
>
>
>
> On Friday, April 11, 2014 1:20:18 PM UTC+2, Sheehan Olver wrote:
>>
>> I'm trying to do movies of an evolving solution to a PDE.  Let's say that 
>> the solution at time steps is stored as columns of an array.Is this 
>> possible?   Preferably in IJulia, but if that's not possible export from 
>> Julia.
>>
>>

[julia-users] Re: Function naming idioms in Julia

[Simon Kornblith]

It's rarely used, but |> may be what you're looking for:

julia> X = zeros(50, 50);

julia> X|>size
(50,50)

julia> X|>length
2500

On Friday, April 11, 2014 6:49:16 PM UTC-4, Ben Racine wrote:
>
> Hi all,
>
> I understand most of the mental maps between the common object-oriented 
> systems (Python, JavaScript) to Julia's multiple dispatch system.
>
> But, when a function really does logically belong to its first argument, I 
> (sometimes) find myself missing the function namespacing inherent to those 
> systems. I find myself wanting to do what one can do in R and inject a '.' 
> into the function name just for the illusion of namespacing. 
>
> I suspect there must be a Julia idiom right under my nose that I'm 
> forgetting.
>
> Any help greatly appreciated.
>
> Thanks!
> Ben Racine
>
>

[julia-users] Re: HypothesisTests question

[Simon Kornblith]

There was a missing method for ApproximateMannWhitneyUTest. If you run 
Pkg.update(), this should be fixed. (I also added some tests to make sure 
that show() isn't completely broken for any of the types.)

Simon

On Saturday, April 12, 2014 12:29:22 PM UTC-4, Iain Gallagher wrote:
>
> Hi 
>
> I'm new to julia - test driving as it were. I'm trying to carry out a 
> Mann-Wihtney U test on two columns of a DataFrame.
>
> # data available at 
> http://www.openintro.org/download.php?file=os2_data&referrer=/stat/extras.php
> using DataFrames
> data_in = readtable("marathon.csv", header=true, separator = ',')
>
> male_data = subset(data_in, :(Gender .== "m"))
> female_data = subset(data_in, :(Gender .=="f"))
>
> # type DataArray has to be converted to type Array (use typeof function to 
> see this)
> # for the MannWhitneyUTest
>
> using HypothesisTests
> mv = vector(male_data["Time"]) 
> fv = vector(female_data["Time"])
> MannWhitneyUTest(mv, fv)
>
> Using the above I get an error after running the test.
>
> no method testname(ApproximateMannWhitneyUTest,) 
> in show at /home/iain/.julia/v0.2/HypothesisTests/src/HypothesisTests.jl:66 
> in anonymous at no file:945 in with_output_limit at show.jl:922 
> in showlimited at show.jl:944 in writemime at repl.jl:2 
> in writemime at multimedia.jl:41 in sprint at io.jl:434 
> in display_dict at /home/iain/.julia/v0.2/IJulia/src/execute_request.jl:24 
>
> HypothesisTests is installed and loaded.
>
> If I try the same without the vector conversion (I can't remember where I 
> read about the need to convert to Array from DataArray for this!) I get the 
> same error. I'm using julia 0.2.1 on Linux Mint (Petra).
>
> Can anyone help with this error?
>
> Thanks.
>
> Iain 
>
>
>

Re: [julia-users] Huge performance boost with USE_SYSTEM_LIBM=1

[Simon Kornblith]

They still look like function calls to me, but given the performance 
difference, I would be surprised if they are as accurate. sin is >4x faster 
for Float64, which is on par with VML with 1 ulp accuracy, but VML has the 
benefit of vectorization. Indeed, if 
thisis
 the right file:

If the rounding mode is round-to-nearest, return sine(x) within a
few ULP.  The maximum error of this routine is not precisely
known.  The maximum error of the reduction might be around 3 ULP,
although this is partly a guess.  The polynomials have small
errors.  The polynomial evaluation might have an error under 1
ULP.  So the worst error for this routine might be under 4 ULP.


On Thursday, April 17, 2014 4:26:07 PM UTC-4, Stefan Karpinski wrote:
>
> Apple's libm is very good, which may explain some of this, but my guess is 
> that LLVM is clever enough to notice calls into the system libm and replace 
> them with x86 instructions for the same functions. The hardware 
> implementations are fast but don't always have the same accuracy as 
> openlibm. Other operating systems don't have nearly as good default libm's 
> – they are often inaccurate, slow, or both.
>
>
> On Thu, Apr 17, 2014 at 4:07 PM, John Travers 
> > wrote:
>
>> For a function which contains many calls to exp, sin, cos I got about a 
>> 4x performance boost by compiling with USE_SYSTEM_LIBM=1. Is this expected? 
>> Why exactly does Julia bundle and default to using its own libm? For 
>> consistency and accuracy? I'm on the latest MacBook Pro with a Haswell 
>> CPU - maybe this is to do with AVX?
>>
>>
>

Re: [julia-users] array with different column types

[Simon Kornblith]

The most performant approach would be to store the columns as vectors in a 
tuple or immutable. DataFrames can be nearly as performant if you:

- Extract columns (df[:mycol]) and index into them whenever possible 
instead of indexing individual elements (df[1, :mycol])
- Add typeasserts when you perform indexing operations 
(df[:mycol]::Vector{Int}), or pass the columns to another function

Otherwise you will incur a slowdown because the compiler doesn't know the 
types.

Simon

On Thursday, April 17, 2014 5:34:24 PM UTC-4, John Myles White wrote:
>
> It's actually possible to place pure Julia vectors in a DataFrame, which 
> might be convenient in this case. But you could always just store the 
> columns in a Vector{Any}, which is what the DataFrame does behind the 
> scenes anyway.
>
>  -- John
>
> On Apr 17, 2014, at 2:27 PM, Stefan Karpinski 
> > 
> wrote:
>
> A DataFrame does seem like a good option, but those have NA support that 
> you may not need. Can you elaborate a little more on the use case? Is it a 
> fixed set of column names and types? Or will you need to support different 
> schemas?
>
>
> On Thu, Apr 17, 2014 at 5:16 PM, Stéphane Laurent 
> 
> > wrote:
>
>> Hello,
>>
>>  I need to deal with some objects represented as arrays whose some 
>> columns are BigFloat, some columns are Int, some columns are logical. Is it 
>> a good idea to use a DataFrame ? Is there a better solution ?This is for a 
>> computationally intensive program.
>>
>
>
>

[julia-users] Re: Bindling a command line utility with a package

[Simon Kornblith]

I think you should be able to put:

#!/usr/bin/env julia

at the top of the file and chmod a+rx. Then you can symlink it to 
/usr/local/bin, add the directory to PATH, etc.

On Tuesday, April 22, 2014 8:07:27 PM UTC-4, andrew cooke wrote:
>
>
> I have a package that also, in a separate file, has a command line 
> utility.  The utility uses ArgParse (which works nicely) and is run by
>
> julia path/file.jl -opt1 -opt2 arg1 arg2
>
> That works fine, but is a bit clunky when the path is something like 
> /home/andrew/.julia/v0.3/CRC/src/crc.jl
>
> Is there a better way of arranging this?  I guess in the far future there 
> will be standalone executables and the like.  But for now, how can I make 
> this easy for a user?
>
> More context - 
> https://github.com/andrewcooke/CRC.jl#from-the-command-line
>
> Thanks,
> Andrew
>
>

Re: [julia-users] Surprising range behavior

[Simon Kornblith]

pi*(0:0.01:1) or similar should work.

On Wednesday, April 23, 2014 7:12:58 PM UTC-4, Peter Simon wrote:
>
> Thanks for the explanation--it makes sense now.  This question arose for 
> me because of the example presented in 
> https://groups.google.com/d/msg/julia-users/CNYaDUYog8w/QH9L_Q9Su9YJ :
>
> x = [0:0.01:pi]
>
> used as the set of x-coordinates for sampling a function to be integrated 
> (ideally over the interval (0,pi)).  But the range defined in x has a last 
> entry of 3.14, which will contribute to the inaccuracy of the integral 
> being sought in that example.  As an exercise, I was trying to implement 
> the interpolation solution described later in that thread by Cameron 
> McBride:  "BTW, another possibility is to use a spline interpolation on the 
> original data and integrate the spline evaluation  with quadgk()".  It 
> seems that one cannot use e.g. linspace(0,pi,200) for the x values, because 
> CoordInterpGrid will not accept an array as its first argument, so you have 
> to use a range object.  But the range object has a built-in error for the 
> last point because of the present issue.  Any suggestions?
>
> Thanks,
>
> --Peter
>
> On Wednesday, April 23, 2014 3:24:10 PM UTC-7, Stefan Karpinski wrote:
>>
>> The issue is that float(pi) < 100*(pi/100). The fact that pi is not 
>> rational – or rather, that float64(pi) cannot be expressed as the division 
>> of two 24-bit integers as a 64-bit float – prevents rational lifting of the 
>> range from kicking in. I worried about this kind of issue when I was 
>> working on FloatRanges, but I'm not sure what you can really do, aside from 
>> hacks where you just decide that things are "close enough" based on some ad 
>> hoc notion of close enough (Matlab uses 3 ulps). For example, you can't 
>> notice that pi/(pi/100) is an integer – because it isn't:
>>
>> julia> pi/(pi/100)
>> 99.99
>>
>>
>> One approach is to try to find a real value x such that float64(x/100) == 
>> float64(pi)/100 and float64(x) == float64(pi). If any such value exists, it 
>> makes sense to do a lifted FloatRange instead of the default naive stepping 
>> seen here. In this case there obviously exists such a real number – π 
>> itself is one such value. However, that doesn't quite solve the problem 
>> since many such values exist and they don't necessarily all produce the 
>> same range values – which one should be used? In this case, π is a good 
>> guess, but only because we know that's a special and important number. 
>> Adding in ad hoc special values isn't really satisfying or acceptable. It 
>> would be nice to give the right behavior in cases where there is only one 
>> possible range that could have been intended (despite there being many 
>> values of x), but I haven't figured out how determine if that is the case 
>> or not. The current code handles the relatively straightforward case where 
>> the start, step and stop values are all rational.
>>
>>
>> On Wed, Apr 23, 2014 at 5:59 PM, Peter Simon  wrote:
>>
>>> The first three results below are what I expected.  The fourth result 
>>> surprised me:
>>>
>>> julia> (0:pi:pi)[end] 
>>> 3.141592653589793 
>>>   
>>> julia> (0:pi/2:pi)[end]   
>>> 3.141592653589793 
>>>   
>>> julia> (0:pi/3:pi)[end]   
>>> 3.141592653589793 
>>>   
>>> julia> (0:pi/100:pi)[end] 
>>> 3.1101767270538954 
>>>
>>> Is this behavior correct? 
>>>
>>> Version info:
>>> julia> versioninfo() 
>>> Julia Version 0.3.0-prerelease+2703  
>>> Commit 942ae42* (2014-04-22 18:57 UTC)   
>>> Platform Info:   
>>>   System: Windows (x86_64-w64-mingw32)   
>>>   CPU: Intel(R) Core(TM) i7 CPU 860  @ 2.80GHz   
>>>   WORD_SIZE: 64  
>>>   BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)   
>>>   LAPACK: libopenblas
>>>   LIBM: libopenlibm  
>>>
>>>
>>> --Peter
>>>
>>>
>>

Re: [julia-users] Surprising range behavior

[Simon Kornblith]

I believe that bug only applies to multiplication/division of integer 
ranges, but it will certainly be good to have it fixed.

On Wednesday, April 23, 2014 10:50:57 PM UTC-4, Steven G. Johnson wrote:
>
>
>
> On Wednesday, April 23, 2014 10:17:23 PM UTC-4, Simon Kornblith wrote:
>>
>> pi*(0:0.01:1) or similar should work.
>>
>
> Actually, that may not work because of 
> https://github.com/JuliaLang/julia/issues/6364
>
> However, this should be fixable (and in fact I just submitted a PR for 
> it). 
>

Re: [julia-users] Surprising range behavior

[Simon Kornblith]

Julia's computation is more accurate because we use pairwise summation 
instead of naive summation. This gives the same results as MATLAB:

x = 0.0; for i = 1:100; x += pi/100; end; x - pi

We don't use summation in Range construction, though. The reason the range 
ends up different is that pi/(pi/100) is very slightly less than 100, and 
MATLAB fudges it because the difference is less than 3 ulp.

Simon

On Thursday, April 24, 2014 2:05:51 PM UTC-4, Hans W Borchers wrote:
>
> There is a strange difference in results between Matlab and Julia that 
> could be the reason why  [0:pi/100:pi]  ends up differently:
>
> matlab> x = ones(100, 1) * (pi/100);
> matlab> sum(x) - pi
> ans =
>   -4.4409e-15
>
>
> julia> x = ones(100)*(pi/100);
> julia> sum(x) - pi
> 1.3322676295501878e-15
>
>
> So it appears that summing pi/100 100-times is smaller than pi in Matlab 
> and larger than pi in Julia.
>
>
> On Thursday, April 24, 2014 5:47:34 PM UTC+2, Stefan Karpinski wrote:
>>
>> On Thu, Apr 24, 2014 at 11:43 AM, andrew cooke  wrote:
>>
>>>
>>> sorry, rational is a stupid suggestion for irrational numbers...
>>>
>>> what i am trying to say is that there is no perfect solution to this 
>>> problem (that i know of).  if matlab works better then it is either luck or 
>>> because they are actually fudging things.  fudging things is appealing but 
>>> usually bites you somewhere down the line.
>>>
>>
>> This is exactly how I feel about it. I still think it may be possible to 
>> improve our FloatRange behavior even further, but I'd like to avoid 
>> "fudging it".
>>
>

[julia-users] Re: idiomatic way to generate variable length list or tuple?

[Simon Kornblith]

You can do:

[1, 2, (flag ? 3 : [])]

or:

tuple(1, 2, (flag ? (3,) : ())...)

On Friday, April 25, 2014 7:35:49 PM UTC-4, andrew cooke wrote:
>
>
> really i'm asking if there's an idiomatic way to do the kind of thing you 
> do with linked lists (usually, in functional languages) in julia...
>
> On Friday, 25 April 2014 20:34:43 UTC-3, andrew cooke wrote:
>>
>> oh, cute.
>>
>> no, it's not what i was looking for, unfortunately.  i just used 1, 2 and 
>> 3 as placeholders.  they could be anything.
>>
>> thanks,
>> andrew
>>
>

Re: [julia-users] Octave diag and Julia diag

[Simon Kornblith]

If diag is passed a vector rather than a matrix, we already give a good 
error message:

julia> diag([1, 2, 3, 4])
ERROR: use diagm instead of diag to construct a diagonal matrix
 in diag at linalg/generic.jl:49

It wouldn't hurt to have this in the docs, though.

On Sunday, April 27, 2014 4:07:52 PM UTC-4, Andreas Noack Jensen wrote:
>
> I agree. It would probably avoid some confusion if the documentation was a 
> little longer and pointed to diagm and Diagonal.
>
>
> 2014-04-27 22:02 GMT+02:00 Ivar Nesje >:
>
>> This difference should be explained in the documentation for diag
>>
>> The current documentation is kind of short:
>>
>> Base.diag(M[, k]) 
>> The "k"-th diagonal of a matrix, as a vector.
>>
>> Ivar
>>
>> kl. 21:54:43 UTC+2 søndag 27. april 2014 skrev John Code følgende:
>>
>>> Thank you.
>>>
>>> On Sunday, April 27, 2014 11:49:12 PM UTC+4, Andreas Noack Jensen wrote:

 Hi John

 In julia, the function diag extract the diagonal of a matrix and if the 
 matrix is rectangular, it extracts the diagonal of the largest square sub 
 matrix. Note that in julia, [1 2 3 4] is not vector but a matrix. To 
 construct a matrix from a vector you can either use the function diagm, 
 which does what you expected diag did,

 julia> diagm([1,2,3,4])
 4x4 Array{Int64,2}:
  1  0  0  0
  0  2  0  0
  0  0  3  0
  0  0  0  4

 but it is often better to use Diagonal, which creates a special 
 Diagonal matrix,

 julia> Diagonal([1,2,3,4])

 4x4 Diagonal{Int64}:
  1  0  0  0
  0  2  0  0
  0  0  3  0
  0  0  0  4


 2014-04-27 21:40 GMT+02:00 John Code :
 >
 > Hi all,
 > I would like to ask why there is a difference between Octave diag 
 function
 > and the function that julia provide. For example, in the following 
 Octave session I get:
 >
 > 
 > octave:1> v = [1 2 3 4]
 > v =
 >
 >1   2   3   4
 >
 > octave:2> a = diag(v)
 > a =
 >
 > Diagonal Matrix
 >
 >1   0   0   0
 >0   2   0   0
 >0   0   3   0
 >0   0   0   4
 > =
 >
 > But in Julia I get:
 >
 > =
 > julia> v = [1 2 3 4]
 > 1x4 Array{Int64,2}:
 >  1  2  3  4
 >
 > julia> a = diag(v)
 > 1-element Array{Int64,1}:
 >  1
 >
 >
 > =
 >
 >
 > Why is this the case and how to get a similar effect of the octave 
 code.
 > Thank you.




 --
 Med venlig hilsen

 Andreas Noack Jensen

>>>
>
>
> -- 
> Med venlig hilsen
>
> Andreas Noack Jensen
>  

Re: [julia-users] Re: Array/Cell - a useful distinction, or not?

[Simon Kornblith]

It seems to me that the massive difference in performance between 
homogeneous and heterogeneous arrays is at least in part a characteristic 
of the implementation and not the language. We currently store 
heterogeneous arrays as arrays of boxed pointers and perform function calls 
on values taken out of them using jl_apply_generic, but I think this could 
be made more efficient. For arrays of types that are sufficiently small or 
sufficiently close in size, e.g. Array{Union(Float32,Float64)}, we could 
store type information and value in the array instead of resorting to 
boxing, and we could create code that branches based on type unions instead 
of resorting to jl_apply_generic. Then (some?) heterogeneous arrays could 
be nearly as fast as homogeneous arrays, modulo branch prediction for the 
type. This is basically how we store DataArrays, which could just be 
Array{Union(T,NA)} if the compiler did this. This is easier for arrays of 
unions of concrete types than it is for arrays of abstract types, since for 
abstract types someone could subtype the element type of the array after 
the array has been constructed and both storage and code need to be able to 
accommodate that. But I'm not sure it's right to structure the language 
based on the performance characteristics of the current implementation 
unless we think they cannot be improved.

Simon

On Wednesday, April 30, 2014 10:08:15 AM UTC-4, Oliver Woodford wrote:
>
>
>
> On Wednesday, April 30, 2014 2:31:43 PM UTC+1, Patrick O'Leary wrote:
>>
>>
>> It's a flexible type system, but it doesn't provide the power of ML or 
>> Haskell. If you really, really want this, do a runtime check:
>>
>> reduce((==), [typeof(el) for el in a])
>>
>
> I feel that the difference between homogeneous and heterogeneous arrays is 
> a very important distinction, and not some odd thing that you might only 
> rarely care about. The distinction has a massive impact on speed. The point 
> of Julia is to be a fast dynamic language. Hiding the distinction under the 
> carpet seems contrary to one of the aims of Julia.
>
> Ivar's suggestion of:
> @assert isleaftype(T)
> is nice, but it doesn't feel quite right to me.
>

Re: [julia-users] what's the maximum number of dimensions for julia arrays?

[Simon Kornblith]

See https://github.com/lindahua/ArrayViews.jl/issues/3

On Thursday, May 1, 2014 4:05:02 PM UTC-4, Tim Holy wrote:
>
> No, there's no limit. You seem to be using NumericExtensions (check the 
> error 
> message carefully), which appears to overwrite the implementation of 
> maximum 
> in base. When I use base Julia, your example works for me. 
>
> --Tim 
>
> On Thursday, May 01, 2014 11:38:49 AM Florian Oswald wrote: 
> > is there anything like a maximum number of dimensions that are 
> supported? 
> > 
> > I'm asking because I get this: 
> > 
> > A = rand(4,4,4,4,4) 
> > maximum(A,5) 
> > ERROR: no method contrank(Colon, Colon, Colon, Colon, Int64) 
> >  in _reducedimimpl! at 
> > /Users/florianoswald/.julia/v0.3/NumericExtensions/src/reducedim.jl:174 
>

[julia-users] Re: Parallelizing file processing using a sharedarray of Uint8's

[Simon Kornblith]

Using different blocks should be fine. Unlike with DArrays, all processes 
have access to all of the memory. I believe localindexes etc. are just for 
convenience. Ideally, you would put the results of parsing the csv file 
directly into a SharedArray, but I'm not sure how easy that would be.

Using mmap_array ought to be faster, though. If you mmap the same file on 
all the workers, they're all supposed to share the same cache. To do this I 
believe you have to call mmap on each worker; I'm pretty sure that anything 
that serializes an mmapped array will copy it, which is not what you want.

If you plan to read the same file many times, you might also consider 
reading it once and then saving it as a JLD file. This should be more 
efficient, because binary data is smaller and doesn't need to be parsed. 
You can use the mmaparrays flag or the hyperslab interface if you only need 
to access chunks of the data at a time. Again, if you want to access chunks 
of the file from several workers on the same system, the best strategy is 
probably read the file using mmaparrays on each worker.

Simon

On Friday, May 2, 2014 5:53:09 PM UTC-4, Douglas Bates wrote:
>
> In another thread I mentioned looking at some largish data sets (10's of 
> gigabytes in the form of a .csv file).  I have made some progress treating 
> the file as a memory-mapped Uint8 array but that hasn't been as effective 
> as I would have hoped.  Using a shared array and multiple processes seems 
> an effective way to parallelize the initial reduction of the .csv file.
>
> The best way I have come up with of getting a file's contents as a shared 
> array is
>
> sm = convert(SharedArray, open(readbytes,"./kaggle/trainHistory.csv"))
>
> It would be convenient to process the contents on line boundaries.  I can 
> determine suitable ranges with something like
>
> function blocks(v::SharedVector{Uint8})
> np = length(v.pids)
> len = length(v)
> bsz = div(len,np)
> blks = Array(UnitRange{Int},np)
> low = 1
> for i in 1:np-1
> eolpos = findnext(v, '\n', i*bsz)
> blks[i] = UnitRange(low, eolpos)
> low = eolpos + 1
> end
> blks[np] = UnitRange(low,len)
> blks
> end
>
> which in this case produces
>
> julia> blocks(sm)
> 8-element Array{UnitRange{Int64},1}:
>  1:794390   
>  794391:1588775 
>  1588776:2383151
>  2383152:3177538
>  3177539:3971942
>  3971943:4766322
>  4766323:5560686
>  5560687:6355060
>
>
> (This is a smaller file that I am using for testing.  The real files are 
> much larger.)
>
> These blocks will be different from what I would get with 
> sm.loc_subarray_1d.  It seems to me that I should be able to use these 
> blocks rather than the .loc_subarray_1d blocks if I do enough juggling with 
> @spawnat, fetch, etc.  Is there anything that would stand in the way of 
> doing so?
>

[julia-users] Re: readcsv returns Array{T, 1}

[Simon Kornblith]

vec should have minimal overhead. (Unlike [:], it doesn't make a copy of 
the underlying data.)

On Sunday, May 11, 2014 6:00:10 PM UTC-4, Ethan Anderes wrote:
>
> Thanks for the response. Still not sure I understand what you mean. 
> readcsv returns Array{T, 2} where T is determined form the file ... so the 
> type of the output does change based on the file. Since column vectors are 
> thought of as 1-d arrays in Julia I would have assumed the Julian way to 
> load a column is to return Array{T,1}. I was just suggesting that readcsv 
> would just return Array{T,d} where T and d are appropriately determined 
> form the file.
>
> Anyhoo, I guess there must be a reason it works the way it does but the 
> whole readcsv(...)[:] thing seemed an unnecessary overhead so I thought I 
> would check to see if I'm using the wrong command.
>

Re: [julia-users] Re: tanh() speed / multi-threading

[Simon Kornblith]

According to my VML.jl 
benchmarks,VML tanh is ~5x 
faster for large arrays even with only a single core. 
(VML.jl is currently only single threaded because I haven't figured out how 
to get multithreading without ccalling into MKL, which could lead to 
problems if Julia is linked against OpenBLAS.) It would be great to have 
similar performance in base Julia, but I don't know of any open source 
vector math libraries that provide <1 ulp accuracy.

Simon

On Sunday, May 18, 2014 10:48:16 AM UTC-4, Andreas Noack Jensen wrote:
>
> The computation of `tanh` is done in openlibm, not openblas, and it is not 
> multithreaded.  Probably, MATLAB uses Intel's Vectorized Mathematical 
> Functions (VML) in MKL. If you have MKL you can do that yourself. It makes 
> a big difference as you also saw in MATLAB. With openlibm I get
>
> julia> @time y = tanh(x);
> elapsed time: 1.229392453 seconds (16096 bytes allocated)
>
> and with VML I get
>
> julia> @time (ymkl=similar(x);ccall((:vdtanh_,Base.libblas_name), Void, 
> (Ptr{Int}, Ptr{Float64}, Ptr{Float64}), &length(x), x, ymkl))
> elapsed time: 0.086282489 seconds (16112 bytes allocated)
>
> It appears that we can get something similar with Tobias' work, which is 
> cool.
>
>
> 2014-05-18 16:35 GMT+02:00 Carlos Becker 
> >:
>
>> Sounds great!
>> I just gave it a try, and with 16 threads I get 0.07sec which is 
>> impressive.
>>
>> That is when I tried it in isolated code. When put together with other 
>> julia code I have, it segfaults. Have you experienced this as well?
>>  Le 18 mai 2014 16:05, "Tobias Knopp" 
>> > 
>> a écrit :
>>
>>  sure, the function is Base.parapply though. I had explicitly imported 
>>> it.
>>>
>>> In the case of vectorize_1arg it would be great to automatically 
>>> parallelize comprehensions. If someone could tell me where the actual 
>>> looping happens, this would be great. I have not found that yet. Seems to 
>>> be somewhere in the parser.
>>>
>>> Am Sonntag, 18. Mai 2014 14:30:49 UTC+2 schrieb Carlos Becker:

 btw, the code you just sent works as is with your pull request branch?


 --
 Carlos
  

 On Sun, May 18, 2014 at 1:04 PM, Carlos Becker wrote:

> HI Tobias, I saw your pull request and have been following it closely, 
> nice work ;)
>
> Though, in the case of element-wise matrix operations, like tanh, 
> there is no need for extra allocations, since the buffer should be 
> allocated only once.
>
> From your first code snippet, is julia smart enough to pre-compute 
> i*N/2 ?
> In such cases, creating a kind of array view on the original data 
> would probably be faster, right? (though I don't know how allocations 
> work 
> here).
>
> For vectorize_1arg_openmp, I was thinking of "hard-coding" it for 
> known operations such as trigonometric ones, that benefit a lot from 
> multi-threading.
> I know this is a hack, but it is quick to implement and brings an 
> amazing speed up (8x in the case of the code I posted above).
>
>
>
>
> --
> Carlos
>  
>
> On Sun, May 18, 2014 at 12:30 PM, Tobias Knopp <
> tobias...@googlemail.com> wrote:
>
>> Hi Carlos,
>>
>> I am working on something that will allow to do multithreading on 
>> Julia functions (https://github.com/JuliaLang/julia/pull/6741). 
>> Implementing vectorize_1arg_openmp is actually a lot less trivial as the 
>> Julia runtime is not thread safe (yet)
>>
>> Your example is great. I first got a slowdown of 10 because the 
>> example revealed a locking issue. With a little trick I now get a 
>> speedup 
>> of 1.75 on a 2 core machine. Not to bad taking into account that memory 
>> allocation cannot be parallelized.
>>
>> The tweaked code looks like
>>
>> function tanh_core(x,y,i)
>>
>> N=length(x)
>>
>> for l=1:N/2
>>
>>   y[l+i*N/2] = tanh(x[l+i*N/2])
>>
>> end
>>
>> end
>>
>>
>> function ptanh(x;numthreads=2)
>>
>> y = similar(x)
>>
>> N = length(x)
>>
>> parapply(tanh_core,(x,y), 0:1, numthreads=numthreads)
>>
>> y
>>
>> end
>>
>>
>> I actually want this to be also fast for
>>
>>
>> function tanh_core(x,y,i)
>>
>> y[i] = tanh(x[i])
>>
>> end
>>
>>
>> function ptanh(x;numthreads=2)
>>
>> y = similar(x)
>>
>> N = length(x)
>>
>> parapply(tanh_core,(x,y), 1:N, numthreads=numthreads)
>>
>> y
>>
>> end
>>
>> Am Sonntag, 18. Mai 2014 11:40:13 UTC+2 schrieb Carlos Becker:
>>
>>> now that I think about it, maybe openblas has nothing to do here, 
>>> since @which tanh(y) leads to a ca

Re: [julia-users] Re: Addition of (unused & uninstatiated) type, slows down program by 25%

[Simon Kornblith]

It may not be desirable, but I think it's expected. The relevant comment in 
inference.jl says:

# don't consider more than N methods. this trades off between
# compiler performance and generated code performance.
# typically, considering many methods means spending lots of time
# obtaining poor type information.
# It is important for N to be >= the number of methods in the error()
# function, so we can still know that error() is always None.
# here I picked 4.


In this case, performance is 20% worse because the compiler can't infer the 
return type of intersect if there are >4 possible matching methods. You can 
restore it to the same level if you add a typeassert to the call to 
intersect.

OTOH, the performance you're losing to failed type inference is nothing 
compared to the performance you're losing to dynamic dispatch. You don't 
know the type of object at compile time (it could be Sphere, Cube, or 
Plane), so the compiler can't optimize the call to intersect into a direct 
call. The following code does not have this problem, and is a little over 
150x faster on my system.

f(object, objects...) = intersect(object) + f(objects...)
f() = 0.0
 
function main()
objects = (Sphere(1.0), Cube(2.0), Plane(3.0))
 
ss = 0.0
 
for i in 1:4000
ss += f(objects...)
end
println("Sum: ",ss)
end

On Saturday, May 31, 2014 12:04:02 PM UTC+2, Mike Innes wrote:
>
> You should definitely open an issue about this – if your timings are right 
> it's definitely not desirable behaviour.
>
> https://github.com/JuliaLang/julia/issues?state=open
>
>
> On 31 May 2014 11:00, mike c > wrote:
>
>> I've narrowed down the problem.  It's not a profiling problem.  Julia 
>> seems to have a step-change in speed when there are too many functions of a 
>> similar signature.
>>
>> I've made a short example that reproduces this slowdown: 
>> http://pastebin.com/iHAa2Cws
>>
>> Run the code once as-is, and then uncomment the intersect() function 
>> which is currently disabled and run it again.   I see a 20% drop in speed. 
>> Note: This intersect function is NEVER actually being called. And the type 
>> it is related to is NEVER INSTANTIATED.
>>
>> I think this probably qualifies as a bug, but it may just be the price to 
>> pay for multiple dispatch when there are too many functions (in this case 5 
>> functions) to choose from.
>>
>
>

Re: [julia-users] Re: Addition of (unused & uninstatiated) type, slows down program by 25%

[Simon Kornblith]

On Saturday, May 31, 2014 3:11:40 PM UTC+2, mike c wrote:

>
>
> On Saturday, 31 May 2014 21:07:58 UTC+10, Simon Kornblith wrote:
>>
>> It may not be desirable, but I think it's expected. The relevant comment 
>> in inference.jl says:
>> In this case, performance is 20% worse because the compiler can't infer 
>> the return type of intersect if there are >4 possible matching methods. 
>> You can restore it to the same level if you add a typeassert to the call to 
>> intersect.
>>
> Thanks.  So return types are inferred when there's up to 4 functions with 
> similar signatures. After that there's the performance drop off. 
>
>>
>> OTOH, the performance you're losing to failed type inference is nothing 
>> compared to the performance you're losing to dynamic dispatch. You don't 
>> know the type of object at compile time (it could be Sphere, Cube, or 
>> Plane), so the compiler can't optimize the call to intersect into a 
>> direct call. The following code does not have this problem, and is a little 
>> over 150x faster on my system.
>>
>
> Unfortunately, my example to demonstrate the performance drop off, doesn't 
> really encapsulate what the intersect methods are doing. The intersect() 
> functions are returning values dependant upon the inputs.  It's more like.
>
> function intersect(self::Sphere, i)
>   return self.a + i
> end
>
> where the same "i" is being sent to each intersect().  And instead of 
> summing, I'm really finding the minimum i.e.   min( intersect(object, i) ) 
> 
>
> In reality, it's even more complex than that, with the input "i" really 
> being a random-ish vector representing a light ray.
>

Yes, this is kind of a toy solution to a toy example. It sounds like what 
you describe is implementable as:

f(i, object, objects..) = min(intersect(object, i), f(i, objects...))
f(i) = 0.0

However, there is a caveat here. The compiler is generating specialized 
code for each potential set of object types, which is fine as long as you 
don't expect too many. Performance will fall for more than a few (7?) 
objects, since the compiler will cease to generate specialized code.

(It's also possible that dynamic dispatch isn't actually a bottleneck in 
your code, depending on how much is happening in the intersect function.)


>> f(object, objects...) = intersect(object) + f(objects...)
>> f() = 0.0
>> ...
>> ss += f(objects...)
>>
>  
> In this simple case, is this just pre-calculating the sum of intersect()s 
> once and then using it repeatedly?
>

Yes, it looks like the optimizer is smart enough to do this, but not quite 
smart enough to elide the loop entirely. Still, even if there were 3 adds 
in the inner loop instead of 1, this code would still be much faster than 
the original.
 

> Thanks for the info & code. Changing objects to a tuple speeds up that 
> original code I posted by a factor of 3!
>

Yes, I noticed that too. Type inference doesn't appear to know that objects 
= [Sphere(1.0), Cube(2.0), Plane(3.0)] will be an array at compile time, so 
iteration becomes slow. You will get the same speedup if you pass that 
array into the function, or if you define it as objects = 
Object[Sphere(1.0), Cube(2.0), Plane(3.0)].

Simon

[julia-users] Re: temporary arrays in loops

[Simon Kornblith]

There may be a more elegant solution, but the sledgehammer approach would 
be to replace the reshape operation with:

tmp = pointer_to_array(pointer(buffer), currentdims)

Obviously don't access tmp after resize! but before setting its value or 
else Julia may segfault.

Simon

On Monday, June 2, 2014 11:34:13 PM UTC+2, Jutho wrote:
>
> Dear Julia users,
>
> I often need to use large temporary multidimensional arrays within loops. 
> These temporary arrays don't necessarily have a fixed size, so allocating 
> them just once before the start of the loop is not really an option. 
> Allocating them within the loop triggers a lot of gc activity, which can 
> have a non-negible impact on the performance.
>
> I was now trying to use the following strategy:
> buffer=zeros(T,prod(firstdims))
> tmp=reshape(buffer,firstdims)
> # initialize tmp
> for i in iter
>  # do something with tmp of previous iteration
>
>  # from here on, previous tmp is no longer needed
>  resize!(buffer,prod(currentdims))
>  tmp=reshape(buffer,currentdims)
>  # initialise new tmp
> end
>
> The idea would be that by using resize!, memory would only be allocated if 
> the current length of tmp exceeds any length encountered before. In 
> addition, I think resize! is directly using memory allocation calls in C 
> and is not passing along the GC. Unfortunately, I get the following error:
>
> *ERROR: cannot resize array with shared data*
>
> * in resize! at array.jl:496*
>
> which makes sense. the old buffer is still used in the previous value of 
> tmp, and I have not expressed that I will no longer need tmp. Is there a 
> possibility to express this and to make this work, or is there a better 
> suggestion of how to tackle this kind of problem?
>
>

[julia-users] ANN: StructsOfArrays.jl

[Simon Kornblith]

Yichao, Oscar, and I were unhappy with the current state of vectorization 
of operations involving complex numbers and other immutables so I decided 
to do something about it. I'm pleased to announce StructsOfArrays.jl 
, which performs the Array 
of Structures -> Structure of Arrays memory layout optimization without 
requiring code changes. This alternative memory layout permits SIMD 
optimizations for immutables for which such optimizations would not 
otherwise be possible or profitable, either because of limitations of the 
Julia codegen and LLVM optimizer or because of the type of the operations 
performed. The benchmark in the README shows that StructsOfArrays can give 
non-negligible speedups for simple operations involving arrays of complex 
numbers.

Simon

[julia-users] Re: ArrayViews setindex! error

[Simon Kornblith]

ArrayViews doesn't support indexing with ranges, vectors, etc. on Julia 
0.3, although this should work on Julia 0.4. (Also on 0.4, SubArrays 
created with sub/slice should be equally fast as ArrayViews, and both 
should be faster than on 0.3.) On 0.3 you need to write an explicit loop to 
set multiple indices with ArrayViews.

Simon

On Thursday, August 13, 2015 at 6:52:03 PM UTC-4, Yee Whye Teh wrote:
>
> Apologies, newbie to Julia here.  
>
> I'm trying to use ArrayViews and ran into simple problem, where only the 
> simplest setindex! usage (1 entry) works with a view, and anything 
> involving more than one entry does not:
>
> *julia> **x = ones(4)*
>
> *4-element Array{Float64,1}:*
>
> * 1.0*
>
> * 1.0*
>
> * 1.0*
>
> * 1.0*
>
>
> *julia> **y=view(x,:)*
>
> *4-element ContiguousView{Float64,1,Array{Float64,1}}:*
>
> * 1.0*
>
> * 1.0*
>
> * 1.0*
>
> * 1.0*
>
>
> *julia> **x[1:2]=5 # ok*
>
> *5*
>
>
> *julia> **y[1]=5 # ok*
>
> *5*
>
>
> *julia> **y[1:2]=6 # not ok*
>
> *ERROR: `setindex!` has no method matching 
> setindex!(::ContiguousView{Float64,1,Array{Float64,1}}, ::Int64, 
> ::UnitRange{Int64})*
>
>
>
> (I realize that there was a previous post related, but for complex 
> setindex! use.  This one is very simple and I'd think should work?)
>
> cheers,
> -yw
>

[julia-users] Re: automatic conversion Pair to Tuple -- failure of type inference?

[Simon Kornblith]

There is no conversion from Pair to Tuple. The construction:

(a, b, c) = d

works for any iterable collection d. The same holds for the for loop 
construction.

There used to be a type inference issue, but I fixed it in 
https://github.com/JuliaLang/julia/pull/12493. The output of code_warntype 
is a little wonky but you can see it knows the types of k and v here:

  s1 = (Base.box)(Base.Int,(Base.add_int)(s1::Int64,k::Int64)::Any)::
Int64 # none, line 7:
  s2 = (Base.box)(Base.Float64,(Base.add_float)(s2::Float64,v::Float64
)::Any)::Float64

The LLVM IR also indicates it's being optimized correctly.

Simon

On Tuesday, September 8, 2015 at 10:10:21 AM UTC-4, vav...@uwaterloo.ca 
wrote:
>
> The following code loops over a Dict:
>
> function test_paircvt1()
> d = Dict(1=>5.5, 3=>-2.2)
> s1 = 0
> s2 = 0.0
> for (k,v) in d
> s1 += k
> s2 += v
> end
> s1, s2
> end
>   
> Two issues to ask about:
>
> (1) Since eltype(d)==Pair{Int,Float64}, Julia is automatically converting 
> Pair{Int,Float64} to Tuple{Int,Float64}.  How is this done?  I looked for 
> the appropriate convert() method but couldn't find it. 
>
> (2)  I ran @code_warntype on the above segment and found something 
> worrisome: the compiler apparently thinks that both k and v are of type 
> Union{Int,Float64}.  Why can't it infer the correct types?  Is there a 
> performance loss from the compiler's inability to know the correct types?
>
> The reason I ask is that I have developed SortedDict for the 
> DataStructures.jl package, and currently eltype{SortedDict{K,V}) is defined 
> to be Tuple{K,V} rather than Pair{K,V} so it is incompatible with Dict. 
>  I'm trying to understand the relevant issues in order to fix this.  I am 
> running Julia 0.4, 15-day-old master.
>
> Thanks,
> Steve Vavasis
>
>

[julia-users] Re: How to capture the output of a system command?

[Simon Kornblith]

readall(`cat test`) or similar

On Friday, September 11, 2015 at 7:56:43 PM UTC-4, J Luis wrote:
>
> Ok, I've spend about an hour around "run" "open", "run(pipeline(..." but 
> no way.
> In Matlab I would do
>
> [status, cmdout] = system(cmd);
>
> but in Julia the most a can reach is to run the command
>
> com = "C:/programs/GraphicsMagick/gm.exe compare -density 200 ...
>
>run(`com')
>
> but I need the result of that execution. 
> How to?
>
> Thanks
>
>
>
>
>

Re: [julia-users] AES (Rijndael) block cipher implementation, performance

[Simon Kornblith]

With some tweaks I got a 12x speedup. Original (using Iain's bench with 
10 iterations):

 0.639475 seconds (4.10 M allocations: 279.236 MB, 1.96% gc time) 
 0.634781 seconds (4.10 M allocations: 279.236 MB, 1.90% gc time)

With 9ab84caa046d687928642a27c30c85336efc876c 

 
from my fork (which avoids allocations, adds inbounds, inlines gf_mult, and 
avoids some branches):

 0.091223 seconds 
 0.090931 seconds

With 3694517e7737fe35f59172666da9971f701189ab 
,
 
which uses a lookup table for gf_mult:

 0.062077 seconds
 0.062132 seconds

With 6e05894856e2bec372b75cd52ae91f36731d2096 
,
 
which uglifies shift_rows! for performance:

 0.052652 seconds 
 0.052450 seconds

There is probably a way to make gf_mult faster without using a lookup 
table, since in many cases it's probably doing the same work several times, 
but I didn't put much thought into it.

Simon

On Saturday, September 12, 2015 at 2:29:49 PM UTC-4, Kristoffer Carlsson 
wrote:
>
> I played around with devectorization and made it allocation free but only 
> got the time down by a factor of 2.
>
> Most of the time is spent in gf_mult anyway and I don't know how to 
> optimize that one. If the C library is using a similar function, maybe 
> looking at the generated code to see what is different.
>

Re: [julia-users] AES (Rijndael) block cipher implementation, performance

[Simon Kornblith]

Values in tuples are either stored in registers or on the stack as decided 
by the compiler, whereas Arrays are currently always allocated on the heap. 
Allocating on the heap is more expensive to begin with, and heap allocated 
objects have to be garbage collected when they're no longer in use.

I did this with Julia 0.4, which as Iain noted above appears to be 
substantially faster at this than 0.3 even without changes.

Simon

On Saturday, September 12, 2015 at 8:33:02 PM UTC-4, Corey Moncure wrote:
>
> Frankly honored and flattered.  That is clever to use the result of the == 
> comparator, which yields a 1 or a 0, to avoid an if/else branch.  Also, 
> using assignment to a tuple of matrix indices to avoid allocating a 
> temporary vector.  The difficulty of mix_columns! is that each calculation 
> depends on all four values in that column of the state, so the state cannot 
> be updated until all the results have been obtained.  I wonder, does this 
> tuple assignment method work substantially differently at the code level?  
> How is the allocation avoided?  By storing in a register?
>
> It looks like you might have access to a different version of Julia than 
> the one I'm using, so I can't verify it myself with your code, but in 0.3 I 
> found using the @inbounds tag decreased performance each time.  
>
> Thanks for your attention. 
>
> On Saturday, September 12, 2015 at 4:38:50 PM UTC-4, Simon Kornblith wrote:
>>
>> With some tweaks I got a 12x speedup. Original (using Iain's bench with 
>> 10 iterations):
>>
>>  0.639475 seconds (4.10 M allocations: 279.236 MB, 1.96% gc time) 
>>  0.634781 seconds (4.10 M allocations: 279.236 MB, 1.90% gc time)
>>
>> With 9ab84caa046d687928642a27c30c85336efc876c 
>> <https://github.com/simonster/crypto/commit/9ab84caa046d687928642a27c30c85336efc876c>
>>  
>> from my fork (which avoids allocations, adds inbounds, inlines gf_mult, and 
>> avoids some branches):
>>
>>  0.091223 seconds 
>>  0.090931 seconds
>>
>> With 3694517e7737fe35f59172666da9971f701189ab 
>> <https://github.com/simonster/crypto/commit/3694517e7737fe35f59172666da9971f701189ab>,
>>  
>> which uses a lookup table for gf_mult:
>>
>>  0.062077 seconds
>>  0.062132 seconds
>>
>> With 6e05894856e2bec372b75cd52ae91f36731d2096 
>> <https://github.com/simonster/crypto/commit/6e05894856e2bec372b75cd52ae91f36731d2096>,
>>  
>> which uglifies shift_rows! for performance:
>>
>>  0.052652 seconds 
>>  0.052450 seconds
>>
>> There is probably a way to make gf_mult faster without using a lookup 
>> table, since in many cases it's probably doing the same work several times, 
>> but I didn't put much thought into it.
>>
>> Simon
>>
>> On Saturday, September 12, 2015 at 2:29:49 PM UTC-4, Kristoffer Carlsson 
>> wrote:
>>>
>>> I played around with devectorization and made it allocation free but 
>>> only got the time down by a factor of 2.
>>>
>>> Most of the time is spent in gf_mult anyway and I don't know how to 
>>> optimize that one. If the C library is using a similar function, maybe 
>>> looking at the generated code to see what is different.
>>>
>>

[julia-users] Re: Regression on a singular matrix

[Simon Kornblith]

\ is supposed to use pivoted QR by default, which a standard way to solve 
least squares when X is not full rank. It gives me a LAPACKException on 0.3 
too, but it works for me on 0.4:

julia> X2 = [1 2 3;
 1 2 5
 1 2 -8
 2 4 23];

julia> y = [2.3, 5.6, 9.2, 10.5];

julia> X2\y
3-element Array{Float64,1}:
  1.28112 
  2.56223 
 -0.146502

You can explicitly tell Julia to solve using pivoted QR, which works on 0.3:

julia> qrfact(X2, pivot=true)\y # on Julia 0.4: qrfact(X2, Val{true})\y
3-element Array{Float64,1}:
  1.28112 
  2.56223 
 -0.146502

Or you can solve using SVD:

julia> svdfact(X2)\y
3-element Array{Float64,1}:
  1.28112 
  2.56223 
 -0.146502

Pivoted QR should be slightly faster than SVD, but both should be much 
faster than rref.

Simon

On Wednesday, September 16, 2015 at 10:02:49 AM UTC-4, Cedric St-Jean wrote:
>
> Sorry, "singular" was the wrong word. I meant that the matrix does not 
> have full rank: the columns are not linearly independent.
>
> # Linearly independent columns
> X1 = [1 2 3;
>   1 2 5
>   1 2 -8
>   1 4 23] 
> # Linearly dependent columns (4th row changed
> X2 = [1 2 3;
>   1 2 5
>   1 2 -8
>   2 4 23]
> y = [2.3, 5.6, 9.2, 10.5]
>
> X1 \ y
> >3-element Array{Float64,1}:
>
> > -8.18265 
> >  6.94133 
> > -0.394898
>
> X2 \ y
> >LAPACKException(2)
>
>
> This is because there is no unique solution to the linear regression 
> problem when the matrix does not have full rank. The usual solution is 
> regularization (eg. ridge regression), but selecting a linearly independent 
> subset of the columns seems like a valid strategy too, if only prediction 
> accuracy matters.
>
> Hence: isn't that a valid non-pedagogical use of rref?
>
> Cédric
>
> On Wednesday, September 16, 2015 at 12:13:55 AM UTC-4, Steven G. Johnson 
> wrote:
>>
>>
>>
>> On Tuesday, September 15, 2015 at 10:22:47 AM UTC-4, Cedric St-Jean wrote:
>>>
>>> In the discussion of issue 9804 
>>> , 
>>> the general consensus seemed to be:
>>>
>>> > Right. "rref" has no use in real numerical coding, but every use in 
>>> pedagogy. 
>>>
>>>
>>>  I had a linear regression yield a PosDefError (the X matrix contains 
>>> only integers and was singular). I solved it by selecting a linearly 
>>> independent subset of columns 
>>> .
>>>  
>>> This uses the Row-Echelon-Form, rref. Is there a simpler way of doing 
>>> this that I'm missing?
>>>
>>
>> A \ x, where A is a non-square matrix, finds least-square solution via 
>> QR.  This is normally the "right" way to do regression.
>>
>> (You can also compute the QR factorization of A explicitly to get a rank 
>> estimate; google "rank revealing QR" for more info on this sort of thing.)
>>
>

[julia-users] Re: Repmat speed

[Simon Kornblith]

Keyword arguments can add additional overhead when calling a function and 
prevent its return type from being inferred properly, but the code that 
runs is specialized for the types of the keyword arguments, so I don't 
think keyword arguments alone explain this. But it looks like there is some 
type instability inside repeat because it creates arrays where the number 
of dimensions can't be inferred.

Simon

On Monday, March 2, 2015 at 5:35:15 AM UTC-5, Ivar Nesje wrote:
>
> The speed difference is because Julia doesn't specialize for specific 
> types for keyword arguments, and type unstable code creates really slow 
> loops in Julia.
>
> Thanks for reporting, I'll work on a fix.
>
> mandag 2. mars 2015 09.09.04 UTC+1 skrev antony schutz følgende:
>>
>> Hi Steven,
>>
>> Thanks for your answer but my question is a bit different.
>> I'm not asking about creating a mesh grid and how or why doing it. 
>> My question is more "naive": 
>> Why the first method is faster than the 2 with only (well written) 1 line 
>> command.
>>
>> and a more general question is: why "repeat" is slower than "repmat"   , 
>> for example: 
>>
>> *tic(); repmat([1:10],1000,1000); toc()*
>>
>> elapsed time: 0.035878247 seconds
>>
>>
>> *tic(); **repeat([1:10],outer=[1000,1000]);** toc()*
>>
>> elapsed time: 1.176858309 seconds
>>
>> Thanks
>>
>>
>> Le vendredi 27 février 2015 15:11:10 UTC+1, antony schutz a écrit :
>>>
>>> Hello, 
>>>
>>> I have a question about the best way to implement a grid similar to a 
>>> mesh grid: 
>>> My first intuition was to do the following: 
>>>
>>> nx = 256
>>> nb = 195
>>>
>>> kx = [-nx/2:-1+nx/2]
>>>
>>> tic()
>>> k1   = repmat(kx,1,nx)
>>> k1v = vec(k1)'#/nx
>>> k1m = repmat(k1v,nb,1)
>>> toc()
>>> # 0.0256 sec 
>>>
>>> Then I tried in one operation the following: 
>>> tic()
>>> ka = repeat(kx,outer=[nx,nb])'# reshape( 
>>> repeat(repeat(kx,inner=[nb]),outer=[nx]) ,nb,nx*nx )
>>> toc()
>>> # 0.477
>>>
>>> Does somebody knows why the repeat is ~20 times slower than the 
>>> repeat/vec/repmat 
>>> Is it the best way to do this ?
>>>
>>> Thanks in advance
>>>
>>> Bests.
>>>
>>> Antony
>>>
>>>

[julia-users] Re: Passing in Dict for Variable Number of Keyword Arguments like Python

[Simon Kornblith]

julia> foo(;d...)
2-element Array{Any,1}:
 (:a,97)
 (:b,95)

Note the semicolon. Otherwise the Dict is splatted as positional arguments.

Simon

On Tuesday, March 31, 2015 at 3:29:47 PM UTC-4, Michael Turok wrote:
>
> Is there a way to pass in a dictionary for the keyword arguments that have 
> a varargs specification, much like python?
>
> Note that this works in julia for passing in arrays/tuples for the 
> non-keyword arguments.
>
> *julia> foo(;kwargs...) = return kwargs*
> *foo (generic function with 1 method)*
>
> *julia> foo(x=1,y=2)*
> *2-element Array{Any,1}:*
> * (:x,1)*
> * (:y,2)*
>
> *julia> d = { :a=> 97, :b => 95 }*
> *Dict{Any,Any} with 2 entries:*
> *  :b => 95*
> *  :a => 97*
>
> *julia> foo(d...)*
> *ERROR: `foo` has no method matching foo(::(Symbol,Int64), 
> ::(Symbol,Int64))*
>
> Python example: 
>  
>
> >>> def foo(**kwargs):  return kwargs
> ... 
> >>> foo(x=1,y=2)
> {'y': 2, 'x': 1}
> >>> dict = { "a" : 97, "b" : 95 }
> >>> foo(**dict)
> {'a': 97, 'b': 95}
> >>> 
>
>

[julia-users] Re: warning with @load in JLD

[Simon Kornblith]

The reason this particular form of @load now shows a warning is that it 
can't be made to work properly in a function. If you can't know either the 
file name or the variable names ahead of time, there are two things you 
could do:

- Use the load function (e.g. x = load("data_run$(run).jld")) rather than 
the macro. This returns a Dict of all variables in the file, so you'd 
access them as x["a"] etc.
- Use @eval @load $("data_run$(run).jld")

Simon 

On Tuesday, August 30, 2016 at 3:05:58 PM UTC-4, Ethan Anderes wrote:
>
> I often use @load as a way to conveniently load a bunch of variables 
> obtained from a simulation run.
> For example, the simulation script might look like this
>
> using JLD
> run   = 1 # <--- simulation parameter
> a,b,c = 1,2,3
> @save "data_run$(run).jld" a b c
>
> Then, to view these simulations later, I would use the load macro. In v0.4 
> my REPL session would look something like this:
>
> julia> using JLD
>
> julia> run = 1  # <--- specify which simulation run to load
> 1
>
> julia> @load "data_run$(run).jld"
> 3-element Array{Symbol,1}:
>  :a
>  :b
>  :c
>
> julia> a
> 1
>
> julia> b
> 2
>
> julia> c
> 3
>
> However, in v0.5.0-rc3+0 , I get the following warning.
>
> julia> using JLD
>
> julia> run = 1  # <--- specify which simulation run to load
> 1
>
> julia> @load "data_run$(run).jld"
> WARNING: @load-ing a file without specifying the variables to be loaded may 
> produce
> unexpected behavior unless the file is specified as a string literal. Future
> versions of JLD will require that the file is specified as a string literal
> in this case.
> 3-element Array{Symbol,1}:
>  :a
>  :b
>  :c
>
> julia> a
> 1
>
> julia> b
> 2
>
> julia> c
> 3
>
> julia> Pkg.status("JLD")
>  - JLD   0.6.3+ master
>
> I was going to file an issue with JLD but I figured I would post here 
> first to see if I’m using an ill-advised workflow or missing some simple 
> work around.
>
> Note: there are a few features of my situation which makes @load 
> particularly convenient.
>
> - Each simulation produces a bunch of variables and some older runs have a 
> smaller subset of variables than newer ones. This makes it clumsy to use the 
> load function (rather than the macro) where one needs to specify the variable 
> names ahead of time.
> - I need to load in multiple `.jld` files, sometimes in a large plotting 
> script, so printing out each file string in the REPL and pasting it into the 
> REPL after `@load` is also unwieldy.
>
> ​
>

Re: [julia-users] Re: Performance of accessing object data

[Simon Kornblith]

There shouldn't be any allocation in the last 3 calls. There is only 
allocation because, at the top of buftest(), you define:

_ = with_alloc(cont)
_ = without_alloc(cont)
_ = without_alloc_with_temp(cont)
_ = using_parameter(cont, buf)

This doesn't actually introduce uncertainty about the types, but it does 
deoptimize things a bit: it causes Julia to box _ when you use it as an 
iteration variable later in buftest(), which is the source of the 
allocation. Simply deleting the LHS here removes the allocation in your 
tests of those three functions, although it doesn't make much of a 
difference in performance.

Also note that you put gc() within the @time block, which may not be what 
you want. With the change above, there isn't really a point to calling gc() 
at all, since only with_alloc allocates.

Simon

On Wednesday, June 25, 2014 3:24:38 PM UTC-4, Spencer Russell wrote:
>
> Mystery solved. In #3 I was missing the indexing `[i]` so I was adding a 
> constant to the whole array instead of just incrementing that item each 
> time.
>
> My final results were (note I'm doing way more trials now to offset the GC 
> time):
>
> ```
> Doing 10 Iterations of each...
> ---
> Timing with allocation each call
> elapsed time: 0.474479476 seconds (417591824 bytes allocated, 42.90% gc 
> time)
> Timing without allocation each call
> elapsed time: 0.35568662 seconds (1591824 bytes allocated, 2.88% gc time)
> Timing without allocation using a temp buffer each call
> elapsed time: 0.18984845 seconds (1591824 bytes allocated, 5.64% gc time)
> Timing passing array as a parameter
> elapsed time: 0.173606393 seconds (1591824 bytes allocated, 6.03% gc time)
> ```
>
> So the really good news is that assigning the object field to a temp 
> variable before my tight loop gives about a 2x speed increase, without 
> modifying my method's API.
>
> Thanks all for the pointers!
>
>
> peace,
> s
>
>
> On Wed, Jun 25, 2014 at 2:57 PM, Spencer Russell  > wrote:
>
>> Hi Thomas,
>>
>> Thanks for the pointer towards @time and the GC info.
>>
>> I also just realized I broke a golden performance rule in my test, which 
>> was referencing variables in global scope.
>>
>> Putting the whole test inside a function gives more reasonable results in 
>> the sense that #2 and #4 do the exact same amount of allocation, and #2 is 
>> a bit faster than #1, but not as fast as #4.
>>
>> ```
>> Timing with allocation each call
>> elapsed time: 0.043889071 seconds (41751824 bytes allocated)
>> Timing without allocation each call
>> elapsed time: 0.026565517 seconds (151824 bytes allocated)
>> Timing without allocation using a temp buffer each call
>> elapsed time: 29.461950105 seconds (42762391824 bytes allocated, 59.40% 
>> gc time)
>> Timing passing array as a parameter
>> elapsed time: 0.01580412 seconds (151824 bytes allocated)
>> ```
>>
>> I'm still a bit surprised that #2 is that much slower than #4, as it 
>> seems like it's just another pointer dereference, and that #3 isn't a fix 
>> for that.
>>  
>> peace,
>> s
>>
>>
>> On Wed, Jun 25, 2014 at 2:39 PM, Tomas Lycken > > wrote:
>>
>>> If you measure time using the @time macro instead of with tic()/toc(), 
>>> you also get information about memory allocation and garbage collection. 
>>> Doing that, I find
>>>
>>> Timing with allocation each call
>>> elapsed time: 0.004325641 seconds (4167824 bytes allocated)
>>> Timing without allocation each call
>>> elapsed time: 0.53675596 seconds (98399824 bytes allocated, 7.60% gc 
>>> time)
>>> Timing without allocation using a temp buffer each call
>>> elapsed time: 2.165323004 seconds (4309087824 bytes allocated, 54.22% 
>>> gc time)
>>> Timing passing array as a parameter
>>> elapsed time: 0.001356721 seconds (7824 bytes allocated)
>>>
>>> so you see that the third method is terribly memory-inefficient, both 
>>> allocating and garbage collecting way more than any other method. The last 
>>> method is much faster since it barely allocates any new memory.
>>>
>>> // T
>>>
>>> On Wednesday, June 25, 2014 7:57:18 PM UTC+2, Spencer Russell wrote:

 I'm having some trouble understanding some performance issues. I put 
 together a minimal example here:

 https://gist.github.com/ssfrr/8934c14d8d2703a3d203

 I had some methods that were allocating arrays on each call, which I 
 figured wasn't very efficient.

 My first attempt to improve things was to allocate an array in my main 
 container type on initialization, and then share that between function 
 calls.

 Suprisingly (to me) this slowed things down by about 60x.

 I wondered if maybe this was because of the extra dereference to get 
 the array (though the slowdown seemed too dramatic for that) so I saved 
 the 
 reference to the array in a temp variable before my tight loop.

 This slowed down by an additional 7x (more surprises!).

 Passing the array as a parameter directly to each function invocation 
>>

[julia-users] Re: include("otherfile.jl") inside a function - otherfile.jl doesn't seem to know about variables defined before include

[Simon Kornblith]

include evaluates at top-level, so this would only work if foo were a 
global variable. It not possible to include in a function context for the 
same reason it is not possible to eval in a function context.

Simon

On Thursday, June 26, 2014 1:03:00 PM UTC-4, Tomas Lycken wrote:
>
> I have the following two files:
>
> *includetest1.jl*:
>
> module IncludeTest
>
> function testinclude()
> foo = "foo"
> println(foo)
> include("includetest2.jl")
> end
>
> end
>
> *includetest2.jl*
>
> println(foo)
>
> If I now try to execute this the function from the REPL, I get errors 
> stating that foo is not defined:
>
> julia> include("includetest1.jl")
>
> julia> IncludeTest.testinclude()
> foo
> ERROR: foo not defined
>  in include at boot.jl:244
> while loading [...]/includetest2.jl, in expression starting on line 1
>
> I thought include was supposed to just insert the contents of the file in 
> whatever context you’re in? If include is not the way to do this, is 
> there another?
>
> For completeness:
>
>
> julia> versioninfo()
> Julia Version 0.3.0-prerelease+3884
> Commit 3e6a6c7* (2014-06-25 10:41 UTC)
> Platform Info:
>   System: Linux (x86_64-linux-gnu)
>   CPU: Intel(R) Core(TM) i7-2630QM CPU @ 2.00GHz
>   WORD_SIZE: 64
>   BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)
>   LAPACK: libopenblas
>   LIBM: libopenlibm
>
> // T
> ​
>

[julia-users] Re: Succinct syntax for removing a column from a DataFrame using the column name

[Simon Kornblith]

delete!(A, [:C,:D])

On Tuesday, July 1, 2014 7:38:00 PM UTC-4, Andre P. wrote:
>
> I figured out one way to do this...
>
> B = A[:, setdiff(names(A),[symbol("D"), symbol("E")])] # removes columns C 
> & D using column names
>
> A less verbose way?
>
> On Wednesday, July 2, 2014 7:01:56 AM UTC+9, Andre P. wrote:
>>
>> I have a data frame with a large number of columns and I would like to 
>> remove certain unnecessary columns as part of the initial data wrangling.
>>
>> I can do this succinctly using the numeric index...
>>
>> using DataFrames
>> A = DataFrame(A=1:3, B=4:6, C=7:9, D=10:12, E=13:15)
>> B = A[:, setdiff(1:end,[3,4])] # returns DFs without columns C & D
>>
>> Is there a way to do this using the column names rather than their 
>> numeric indices?
>>
>> Andre
>>
>>
>>
>>

Re: [julia-users] type matching Array{Tuple}

[Simon Kornblith]

See also https://github.com/JuliaLang/julia/issues/6561

On Monday, July 14, 2014 10:48:42 AM UTC-4, Mauro wrote:
>
> Using type parameters works for me: 
>
> julia> f{T<:Float64}(a::Array{(T,T)}) = eltype(a) 
> f (generic function with 1 method) 
>
> julia> f([(3.,4.), (4.,5.)]) 
> (Float64,Float64) 
>
> julia> f{T<:Tuple}(a::Array{T}) = eltype(a) 
> f (generic function with 2 methods) 
>
> julia> f([(3,4), (4,5)]) 
> (Int64,Int64) 
>
> The reason is invariance of parametric types.  Tuple is an abstact type, 
> this means that for parametric types (which arrays also belong to): 
> "Even though (Float64,Float64) <: Tuple we DO NOT have 
> Array{(Float64,Float64)} <: Array{Tuple}" 
>
> (citing the manual: 
>
> http://docs.julialang.org/en/latest/manual/types/#parametric-composite-types) 
>
>
> (A side note: tuples themselves are covariant.  E.g. (Int,) <: (Real,) ) 
>
> Last, why this worked for tuples before, I don't know.  Presumably it 
> was a bug and now got fixed. 
>
> On Mon, 2014-07-14 at 12:01, David van Leeuwen  > wrote: 
> > Hello, 
> > 
> > Somewhere in the last 150 days the value of 
> > 
> > Array{(Float64,Float64)} <: Array{Tuple} 
> > 
> > changed from "true" to "false".  There may be good reasons for that, I 
> > can't tell. 
> > 
> > I would like to specify as a function argument something of the type 
> > Array{Tuple}, but I can't get that matched any more.  I don't know how 
> to 
> > make the function type more specific so that it will match again, i.e., 
> I 
> > don't know how to specify the types withing the tuple.  A naive 
> > Array{Tuple{Float64,Float64}} does not work. 
> > 
> > Any ideas? 
> > 
> > Thanks 
> > 
> > ---david 
>

Re: [julia-users] Append! resulting in #undef elements. Intended behavior?

[Simon Kornblith]

The fact that append! grows the array on failure seems like a bug 
nonetheless. If convert throws it seems preferable to leave the array as 
is. I'll file an issue.

Simon

On Thursday, July 17, 2014 9:34:21 AM UTC-4, Jacob Quinn wrote:
>
> Hi Jan,
>
> You have your syntax a little mixed up. The usage of:
>
> Type[]
>
> actually declares an empty array with element type of `Type`. So you're 
> first line:
>
> x = Array[]
>
> is actually creating an array of arrays.
>
> Likewise, you're seeing the error in
>
> Array[1]
>
> Because you're trying to put an Int[1] array into an array of arrays (and 
> julia doesn't know how to make that conversion).
>
> The last error is unfortunate, because it seems that the call to `append!` 
> is allocating space for the array you're appending but then fails when 
> actually trying to put the 2nd array into the newly allocated space. 
> Because `append!` is mutating, you're left with your original array mutated 
> with the extra space, with the #undef.
>
> I think the semantics you're looking for are as follows:
>
> x = Int[]  # declares an empty 1-d array with element type of `Int`
>
> y = [1, 2, 3]  # literal syntax for creating an array with elements, (1, 
> 2, 3). Type inference figures out that the elements are of type `Int`
>
> append!(x,y)  # mutates the `x` array by appending all the elements of y 
> to it; this works because they're both of the same type
>
> push!(x, 5)  # adds a single element, 5,  to the `x` array 
>
> Feel free to read through the section in teh manual on arrays to get a 
> better feel for what's going on.
>
> http://docs.julialang.org/en/latest/manual/arrays/
>
> Cheers,
>
> -Jacob
>
>
> On Thu, Jul 17, 2014 at 8:21 AM, Jan Strube  > wrote:
>
>>
>>
>> Hi List!
>>
>> I'm a particle physicist just getting started with some julia. I've been 
>> using some python in the past but now I'm trying to use analyzing some lab 
>> data as an excuse to learn Julia.
>> So far it looks like I'm going to stick with it for a while.
>>
>> I've been trying to play with basic image analysis, and I've come across 
>> the following behavior: append! complains that it doesn't find a suitable 
>> method for my call signature, but it does append an #undef. Is this 
>> intended?
>>
>> Please see below for a session.
>>
>>
>>
>>_   _ _(_)_ |  A fresh approach to technical computing
>>   (_) | (_) (_)|  Documentation: http://docs.julialang.org
>>_ _   _| |_  __ _   |  Type "help()" to list help topics
>>   | | | | | | |/ _` |  |
>>   | | |_| | | | (_| |  |  Version 0.3.0-rc1 (2014-07-14 02:04 UTC)
>>  _/ |\__'_|_|_|\__'_|  |  
>> |__/   |  x86_64-apple-darwin13.3.0
>>
>> julia> x = Array[]
>> 0-element Array{Array{T,N},1}
>>
>> julia> append!(x, Array[1])
>> ERROR: `convert` has no method matching convert(::Type{Array{T,N}}, 
>> ::Int64)
>>  in getindex at array.jl:121
>>
>> julia> x
>> 0-element Array{Array{T,N},1}
>>
>> julia> append!(x, Array[[1]])
>> 1-element Array{Array{T,N},1}:
>>  [1]
>>
>> julia> append!(x, [1])
>> ERROR: `convert` has no method matching convert(::Type{Array{T,N}}, 
>> ::Int64)
>>  in copy! at abstractarray.jl:197
>>  in append! at array.jl:475
>>
>> julia> x
>> 2-element Array{Array{T,N},1}:
>> [1]
>>  #undef
>>
>
>

Re: [julia-users] comprehension and closures in 0.3.0-rc1

[Simon Kornblith]

Bug filed: https://github.com/JuliaLang/julia/issues/7679

On Sunday, July 20, 2014 6:08:46 PM UTC-4, Keno Fischer wrote:
>
> I don't recall that being changed again intentionally. 
>
> On Sun, Jul 20, 2014 at 2:18 PM, Cyril Slobin  > wrote: 
> > Sorry, I've did it wrong the first time, the second attempt. Seems like 
> the 
> > behavior of comprehensions in 0.3.0-rc1 was changed (compared to stable 
> > 0.2.1), and changed in a wrong way. The expected behavior is documented 
> > here: 
> > 
> > 
> https://github.com/JuliaLang/julia/blob/master/NEWS.md#new-language-features-1
>  
> > 
> > What I see on my installation: 
> > 
> > julia> VERSION 
> > v"0.2.1" 
> > 
> > julia> map(f->f(), { ()->i for i=1:3 }) 
> > 3-element Array{Any,1}: 
> >  1 
> >  2 
> >  3 
> > 
> > But: 
> > 
> > julia> VERSION 
> > v"0.3.0-rc1+61" 
> > 
> > julia> map(f->f(), { ()->i for i=1:3 }) 
> > 3-element Array{Any,1}: 
> >  3 
> >  3 
> >  3 
> > 
> > Is it a temporary regression in this specific rc, or an intentional 
> language 
> > change? Or I am just misunderstanding something? Thanks in advance! 
>

[julia-users] Re: Checking for Undefined Reference.

[Simon Kornblith]

isdefined(myobject, :myfield)

On Tuesday, July 22, 2014 6:54:00 PM UTC-4, Ben Ward wrote:
>
> I have a type that contains references to other instances of the same type 
> - like a doubly linked list or - as it's intended use, like a tree like 
> structure. They contain references to a parent of the same types and an 
> array of references to children of the same type. However a root or singly 
> constructed instance of the type will not start out with a parent or 
> children and so these fields will be #undef, how can I check for #undef 
> fields? It didn't say in the types section in the manual, although it did 
> say accessing it is an automatic error so I guess a try block might work - 
> but I don't know if try blocks are preferable in Julia to simply 
> anticipating the #undef. I've normally heard the mantra of anticipating an 
> expected exception rather than make excessive try catches.
>
> Thanks. 
>

Re: [julia-users] simultaneous sin and cos: sincos in julia from openlibm?

[Simon Kornblith]

Presumably we could use the same global or let-scoped array rather than 
allocating a new array on each call.

On Monday, July 28, 2014 11:39:44 AM UTC-4, Simon Byrne wrote:
>
> Yes, I would agree: as Elliot mentioned, you might get some gain by only 
> doing the range-reduction once.
>
> Looking at the source of openlibm/src/s_sincos.c, it seems that's what it 
> does, as well as calculating z = x*x and w= z*z once...
>
> Is there a more efficient way to get return reference values than using 
> arrays?
>
>
> On Monday, 28 July 2014 16:29:05 UTC+1, Stuart Brorson wrote:
>>
>> Yup, you certainly don't want to use CORDIC.  My e-mail is just a 
>> historical note. 
>>
>> The standard (fdlibm) implementation of sin & cos involves 
>> summing a polynomial with 6 or 8 coefficients (after folding down to 
>> the first quadrant).  Here's the kernel: 
>>
>> http://www.netlib.org/fdlibm/k_sin.c 
>>
>> This impl seems pretty quick, IMO.  I'd wager that 
>> there's probably not much room for improvement over the existing, 
>> separate implementations. 
>>
>> Stuart 
>>
>>
>> On Mon, 28 Jul 2014, Stefan Karpinski wrote: 
>>
>> > This CORDIC explanation seems quite plausible. However, it seems like 
>> we 
>> > probably don't want to use this algorithm, which means that we're not 
>> at 
>> > this point able to get any kind of speedup for sincos(x) relative to 
>> > sin(x), cos(x). 
>> > 
>> > 
>> > On Mon, Jul 28, 2014 at 10:13 AM, Stuart Brorson  
>> wrote: 
>> > 
>> >> A bell rang in the back of my head as I was on my way to work this 
>> >> morning.  I was thinking about sincos again, and remembered something 
>> >> about CORDIC algorithms from the distant past.  These are add and 
>> >> shift algorithms used to compute certain trig and other elementary 
>> >> functions.  They were very popular for scientific calculators back 
>> when 
>> >> hand-held calculators were new since they are easily implementable in 
>> >> hardware, and don't require floating point multiply.  A couple of 
>> >> references: 
>> >> 
>> >> http://en.wikipedia.org/wiki/CORDIC 
>> >> http://ww1.microchip.com/downloads/en/AppNotes/01061A.pdf 
>> >> 
>> >> It appears that a common CORDIC computation will product sin & cos 
>> >> together.  I'll bet dollars to doughnuts (without actually knowing) 
>> >> that the x87 assembly instruction mentioned below was doing a CORDIC 
>> >> computation, and it made sense to return both sin & cos since 
>> >> they were computed together. 
>> >> 
>> >> The paper by Jeannerod & JourdanLu refer to CORDIC methods, but is 
>> >> apparently an extension, as far as I can tell. 
>> >> 
>> >> Stuart 
>> >> 
>> >> 
>> >> 
>> >> 
>> >> On Mon, 28 Jul 2014, Simon Byrne wrote: 
>> >> 
>> >>  I've often wondered this myself. As I understand it, the purpose of 
>> the 
>> >>> sincos function was to call the FSINCOS assembly instruction for x87 
>> FPU. 
>> >>> On modern processors however, it is generally acknowledged that 
>> calling a 
>> >>> well-written math library compiled to use SSE instructions is 
>> typically 
>> >>> faster (and can be more accurate) than using x87 trig instructions. 
>> See 
>> >>> this discussion: 
>> >>> http://stackoverflow.com/questions/12485190/calling- 
>> >>> 
>> fsincos-instruction-in-llvm-slower-than-calling-libc-sin-cos-functions 
>> >>> 
>> >>> Calling sincos using Isaiah's method seems to be about 9 times slower 
>> than 
>> >>> calling the sin and cos separately: 
>> >>> https://gist.github.com/734dcacde1f107397b3b.git 
>> >>> though a lot of this seems to be due to the overhead in creating and 
>> >>> destroying the arrays for return values. 
>> >>> 
>> >>> 
>> >>> On Monday, 28 July 2014 13:34:53 UTC+1, Kevin Squire wrote: 
>> >>> 
>>  
>>  This paper seems relevant, though possibly only for 32-bit: 
>>  
>>  
>>  http://hal.archives-ouvertes.fr/docs/00/67/23/27/PDF/ 
>>  Jeannerod-JourdanLu.pdf 
>>  
>>  Cheers, 
>> Kevin 
>>  
>>  On Monday, July 28, 2014, Stuart Brorson > > 
>>  
>>  wrote: 
>>  
>>   On Sun, 27 Jul 2014, Viral Shah wrote: 
>> > 
>> >  Is sincos a standard libm function? 
>> > 
>> >> 
>> >> 
>> > Out of curiosity I looked into sincos since I had never heard of 
>> it. 
>> > A quick check shows there's no sincos in fdlibm 
>> > (on netlib).  However, a little Googling reveals an old Sun math 
>> > library libsunmath seems to implement it. 
>> > 
>> > I did find a couple of libm variants which implemented sincos. 
>> > However, they simply called sin & cos separately.  As Stephan says 
>> > upthread, no performance improvement. 
>> > 
>> > As far as I know, sin & cos are usually computed using mod to fold 
>> > the input x down to the first quadrant, and then using a power 
>> series 
>> > (needs only 6 or 8 terms IIRC) to compute the function.  Perhaps 
>> > libsunmath computed e.g. cos first, and then did sin = sqrt(1 - 
>> 

[julia-users] Re: random machine freeze

[Simon Kornblith]

I suspect you are running out of RAM and your system is thrashing 
.

On Thursday, July 31, 2014 3:06:31 AM UTC-4, K leo wrote:
>
> Sorry, but really don't know what is going on.  Today I had two sessions 
> of julia each running a program (single processing), then the machine 
> froze with no mouse and keyboard.  This type of thing happened a few 
> times in the past though not every time I run 2 julia sessions.  I don't 
> know how to determine the real cause, but suspect it has to do with 
> julia. Hasn't anyone had similar experiences? 
>
> I am running Xubuntu 14.04. 
>
> $ julia 
> _ 
> _   _ _(_)_ |  A fresh approach to technical computing 
>(_) | (_) (_)|  Documentation: http://docs.julialang.org 
> _ _   _| |_  __ _   |  Type "help()" to list help topics 
>| | | | | | |/ _` |  | 
>| | |_| | | | (_| |  |  Version 0.3.0-rc1+54 (2014-07-17 05:40 UTC) 
>   _/ |\__'_|_|_|\__'_|  |  Commit 4e92487 (14 days old master) 
> |__/   |  x86_64-linux-gnu 
>
>
>

Re: [julia-users] OptionTypes.jl

[Simon Kornblith]

 types 
> currently are being boxed because their type could be easily inferred by 
> the compiler.
>
> If you're not familiar with the requirement for breaking our 
> current abstraction barriers, note that indexing into a DataArray at the 
> moment poisons the performance of every program you write, because the 
> result has an uncertain type that the compiler doesn't know how to 
> optimize. To work around this, you have to write code that effectively 
> accesses the raw .na and .data fields of a DataArray. Simon's done some 
> great work to make this easier to do, but I'm not sure that's the right 
> direction for us to head in the long-run.
>
> In general, I think an approach to missing data built on the combination 
> of Option{T} and DataArray{T} provides an interface that's simple and 
> consistent (everything happens in terms of isna/get) even if it's an 
> interface that's somewhat unfamilar to R/Python folks. Most important to me 
> is that using Option types efficiently doesn't require a deep understanding 
> of Julia's type system, whereas our current abstractions require you to 
> understand how to work around problems raised by Union types when 
> programming for the current compiler. An Option type is basically a forcing 
> function that says: "Julia is aggressively typed. If you want to work with 
> a missing value of type T, you need to explicitly say how you're going to 
> handle any missingness so that the system only interacts with values of 
> type T."
>
> I should note that I'm not very sure the use of Options is the right 
> approach: Simon Kornblith has argued very persuasively for waiting for the 
> compiler to improve its ability to handle tagged unions like those 
> generated by indexing into DataArrays. My personal feeling is that it's 
> easier to expose a simple abstraction that doesn't assume the compiler will 
> change dramatically. This is based on my aesthetic sense that Julia's power 
> comes from making optimizations transparent and explicit, rather than 
> making the compiler smarter.
>
> It's also worth noting that there are problems for which the use of 
> Option types isn't very helpful: the computation of medians, for example, 
> isn't defined in terms of scalars, so having a better abstraction for 
> expressing missing scalars won't get us anywhere.
>
> One other caveat is that I'm not providing an `unsafe_get` method, which 
> means that the `isna` followed by `get` idiom I showed above does two 
> checks when you could get away with only one. I still haven't figured out 
> how I'd like to handle that issue.
>
> So there's still a lot of design work needed, but I wanted to let people 
> see how this interface could work were we to choose it.
>
>  -- John
>
> On Jul 31, 2014, at 8:15 AM, Harlan Harris  > wrote:
>
> John, how might this interact with DataArrays? This design, unlike 
> DataArrays, requires that you use an entire byte to store the missingness, 
> so it's not likely to be as fast if you're manipulating a lot of them. Is 
> the intended use case here something sum(DataArray) yields Option{Float64}? 
>
>
> On Thu, Jul 31, 2014 at 11:01 AM, Stefan Karpinski  > wrote:
>
>> This looks quite promising. The `get` interface looks like a very nice 
>> way to generically deal with missing values – provide a default or get an 
>> error. Automatic conversion of the default to the type of the Option value 
>> is particularly nice. This seems like it will be a pleasant and efficient 
>> API.
>>
>> Minor question: how come the non-NA constructor for Option takes both the 
>> `na` and `value` arguments? Doesn't supplying a value imply that it's not 
>> NA while not supplying a value implies that it is NA?
>>
>>
>> On Thursday, July 31, 2014, John Myles White > > wrote:
>>
>>> At JuliaCon, I described the idea of using Option types instead of 
>>> NAtype to make it easier to write type-stable code that interacts with 
>>> NA’s. To help facilitate a debate about the utility of this approach, I 
>>> just wrote up a minimal package that implements Option types: 
>>> https://github.com/johnmyleswhite/OptionTypes.jl
>>>
>>> Essentially, an Option type is a scalar version of a DataArray: it 
>>> contains a value and a Boolean mask indicating whether that value is NA or 
>>> not. Because Option types are parametric, they allow us to express the 
>>> variants of NA that R has, such as a Float64 specific NA and an Int64 
>>> specific NA.
>>>
>>>  — John
>>>
>>>  
>
>

Re: [julia-users] OptionTypes.jl

[Simon Kornblith]

On Friday, August 1, 2014 6:23:59 AM UTC-4, Milan Bouchet-Valat wrote:
>
>  Le jeudi 31 juillet 2014 à 21:19 -0700, John Myles White a écrit :
>
> To address Simon’s general points, which are really good reasons to avoid 
> jumping on the Option{T} bandwagon too soon: 
>
>  
>
>  * I agree that most languages use tagged union types for Option{T} 
> rather than a wrapper type that contains a Boolean value. It’s also totally 
> true that many compilers are able to make those constructs more efficient 
> than Julia currently does. But what we should expect from Julia in the 
> coming years isn’t so clear to me. (And I personally think we need to 
> settle on a solution for representing missing data that’s viable in a year 
> rather than viable in five years.) This is an issue that I’d really like to 
> have input on from Jeff, Keno, Jameson or someone else involved with the 
> internals of the compiler. Getting input from the broader community is the 
> main reason I wanted to put a demo of OptionTypes.jl out in front of other 
> folks. 
>
>  
>
>  * I’m not clear how we could come to know that a datum is not missing 
> without a resolution step that’s effectively equivalent to the get() 
> function for Option{T}. I agree that the enforced use of get() means that 
> you can’t hope to use generic functions like sum on collections of 
> Option{T}. But I’m also not sure that’s such a bad thing: I think the 
> easiest way to express to the compiler that you know that all of the 
> entries of a DataArray are not NA is to convert the DataArray to a straight 
> Array. But maybe you have other mechanisms for expressing this knowledge. 
> Certainly my proposal to do conversions to Arrays isn’t the most elegant 
> strategy. It’s just all that I’ve got so far. 
>
>  
>
>  * I kind of like the idea of Option{T} standing outside of the main type 
> system in a kind of mirror type system. I’m less happy about Union(NA, T) 
> being a super type of T, even though there are some good reasons that you’d 
> like to view T as a specialization of Union(NA, T). But I agree that I 
> don’t have a good feel about where missing data belongs in the type 
> hierarchy. This is another question for which I’d love to get input from 
> others. 
>
>  
>
>  In regard to Simon’s performance points: 
>
>  
>
>  * Yes, memory usage alone argues strongly for working with DataArray{T} 
> rather than Array{Option{T}}. 
>
>  
>
>  * Exploting tricks that make operations like anyna() faster is another 
> good argument for keeping DataArray{T} around. 
>
>  
>
>  * I’m not sure how to deal with inlining concerns or the undefined 
> reference checks. Do you have ideas for improving this within DataArrays or 
> do we need supporting changes in the compiler? 
>
> Actually it seems it would be possible to make Array{Union(NAtype, T)} 
> more similar to and as efficient as DataArray{T}, by handling a few 
> things in the compiler. This would create a generalization of DataArray 
> to any kind of union type, which could be useful in other contexts. But 
> more importantly, it would make missing values integrate seamlessly into 
> Julia, getting rid of any hacks.
>
> More specifically, the following features would need to be supported:
> - a way of telling the compiler to store the data as two arrays of 
> concrete types (here T and NAtype), instead of as an array of boxed 
> values, so that:
> * efficient operations can be performed on the T values (by skipping 
> the missing ones manually)
> * T values are stored as a dense array and can be converted to 
> Array{T} without any copy or passed to BLAS when no missing values are 
> present
>* NA values can be packed in a BitArray to save memory and make NA 
> detection faster (see below)
> - a fonction to check whether a given element of the array is of type T 
> rather than of NAtype (generalization of isna())
> - a fonction to check whether all elements of the array are of type T 
> rather than of NAtype (generalization of anyna(), more efficient than 
> calling the previous function on all elements thanks to the packing of NAs 
> in a BitArray) 
> In this scheme, what is missing is how to allow the compiler to pack NAs 
> in a BitArray. Somehow, NAtype would have to be defined as a 1-bit 
> object. Maybe by making it an enum-like immutable with a 1-bit field inside 
> it?
>
> How does it sound? 
>

I've thought a bit about this, but it seems like it would be too much 
complexity in the compiler. Storing arrays as something besides contiguous 
elements and interaction between the codegen in C/C++ and the BitArray code 
in Julia both seem likely to be painful, although Jeff, Keno, and Jameson 
would know better than I. Additionally, this optimization (of storage of 
arrays of unions of a singleton type and a bits type) seems pretty specific 
to DataArrays, but the actual advantages in terms of performance and 
expressibility would be small or non-existent. (This is in contrast to 
optimizing sto

[julia-users] Re: Using a callback from a macro

[Simon Kornblith]

Assuming you're generating code that calls fn, as opposed to trying to call 
it when generating code in the macro (usually a bad idea), you should use 
esc(fn) in place of fn

On Friday, August 1, 2014 10:50:15 AM UTC-4, Abe Schneider wrote:
>
> I've come across a problem where I have macro to define  a grammar 
> (similar to how PEGParser currently works):
>
> @grammar foo begin
>   rule[fn] = "some" + ("rule" | "test")
> end
>
> where the `[fn]` next to the rule defines a function to call on the 
> results (in this case it's an expansion). The issue is that parsing the 
> Expr tree, `fn` is given as a Symbol (which makes sense).
>
> However, if I try to turn `fn` into a function I run into the namespace 
> issue I've had previously. If `fn` is defined in my module, it works 
> without problem. If it's defined in the code that imports the module, it 
> will not work because that function doesn't exist in the namespace of the 
> module.
>
> I'm guessing there isn't an easy solution to fix this problem, but I 
> thought I'd check to see if someone had an idea.
>
> Thanks!
>
>

[julia-users] Re: Using a callback from a macro

[Simon Kornblith]

That looks like the output of :(esc(fn)), but you don't want to quote that, 
you want to evaluate it when generating the code, i.e. use $(esc(fn)) 
wherever you were previously using $fn

On Friday, August 1, 2014 11:38:44 AM UTC-4, Abe Schneider wrote:
>
> That's correct, I'm generating code that keeps a pointer to that function 
> to call later on.
>
> If I do that, I have the type `esc(fn)` in the macro (this is what I get 
> when printing it out), and taking a dump of that argument gives:
>
> Expr 
>   head: Symbol call
>   args: Array(Any,(2,))
> 1: Symbol esc
> 2: Symbol fn
>   typ: Any
>
>
> Does that mean instead of carrying around a function I will need to carry 
> around the Expr?
>
> Thanks!
>
> On Friday, August 1, 2014 11:17:50 AM UTC-4, Simon Kornblith wrote:
>>
>> Assuming you're generating code that calls fn, as opposed to trying to 
>> call it when generating code in the macro (usually a bad idea), you should 
>> use esc(fn) in place of fn
>>
>> On Friday, August 1, 2014 10:50:15 AM UTC-4, Abe Schneider wrote:
>>>
>>> I've come across a problem where I have macro to define  a grammar 
>>> (similar to how PEGParser currently works):
>>>
>>> @grammar foo begin
>>>   rule[fn] = "some" + ("rule" | "test")
>>> end
>>>
>>> where the `[fn]` next to the rule defines a function to call on the 
>>> results (in this case it's an expansion). The issue is that parsing the 
>>> Expr tree, `fn` is given as a Symbol (which makes sense).
>>>
>>> However, if I try to turn `fn` into a function I run into the namespace 
>>> issue I've had previously. If `fn` is defined in my module, it works 
>>> without problem. If it's defined in the code that imports the module, it 
>>> will not work because that function doesn't exist in the namespace of the 
>>> module.
>>>
>>> I'm guessing there isn't an easy solution to fix this problem, but I 
>>> thought I'd check to see if someone had an idea.
>>>
>>> Thanks!
>>>
>>>

Re: [julia-users] OptionTypes.jl

[Simon Kornblith]

Is there a reason we can't change the way unions of small bits types are 
represented, so that if we know something is a Union(Float64,NA) it can 
live in registers or on the stack instead of having to be heap allocated?

On Friday, August 1, 2014 2:54:49 PM UTC-4, Jameson wrote:
>
> I could (and probably will, someday) revive that commit. At the time, 
> though, I seemed to find that it provided little performance benefit -- the 
> gc cost of allocating boxes was far greater (for type uncertainty involving 
> bitstypes) and the type dispatch wasn't as much of a performance impact as 
> I had previously assumed. 
>
>
> On Friday, August 1, 2014, Keno Fischer  > wrote:
>
>> It is possible to do generic compiler improvements for Union types
>> (Jameson had a branch at some point that did callsite splitting if we
>> inferred a Union type). However, I think the best way to go here is to
>> maintain the current separation of two arrays (one of the values one
>> for the NAs), but give an option type on access. The option type would
>> then most likely be in memory and wouldn't have much overhead. Please
>> let me know if there's anything specific I should explain how the
>> compiler will handle it, I admit I have only skimmed this thread.
>>
>> On Fri, Aug 1, 2014 at 9:18 AM, Simon Kornblith  
>> wrote:
>> > On Friday, August 1, 2014 6:23:59 AM UTC-4, Milan Bouchet-Valat wrote:
>> >>
>> >> Le jeudi 31 juillet 2014 à 21:19 -0700, John Myles White a écrit :
>> >>
>> >> To address Simon’s general points, which are really good reasons to 
>> avoid
>> >> jumping on the Option{T} bandwagon too soon:
>> >>
>> >>
>> >>
>> >> * I agree that most languages use tagged union types for Option{T} 
>> rather
>> >> than a wrapper type that contains a Boolean value. It’s also totally 
>> true
>> >> that many compilers are able to make those constructs more efficient 
>> than
>> >> Julia currently does. But what we should expect from Julia in the 
>> coming
>> >> years isn’t so clear to me. (And I personally think we need to settle 
>> on a
>> >> solution for representing missing data that’s viable in a year rather 
>> than
>> >> viable in five years.) This is an issue that I’d really like to have 
>> input
>> >> on from Jeff, Keno, Jameson or someone else involved with the 
>> internals of
>> >> the compiler. Getting input from the broader community is the main 
>> reason I
>> >> wanted to put a demo of OptionTypes.jl out in front of other folks.
>> >>
>> >>
>> >>
>> >> * I’m not clear how we could come to know that a datum is not missing
>> >> without a resolution step that’s effectively equivalent to the get()
>> >> function for Option{T}. I agree that the enforced use of get() means 
>> that
>> >> you can’t hope to use generic functions like sum on collections of
>> >> Option{T}. But I’m also not sure that’s such a bad thing: I think the
>> >> easiest way to express to the compiler that you know that all of the 
>> entries
>> >> of a DataArray are not NA is to convert the DataArray to a straight 
>> Array.
>> >> But maybe you have other mechanisms for expressing this knowledge. 
>> Certainly
>> >> my proposal to do conversions to Arrays isn’t the most elegant 
>> strategy.
>> >> It’s just all that I’ve got so far.
>> >>
>> >>
>> >>
>> >> * I kind of like the idea of Option{T} standing outside of the main 
>> type
>> >> system in a kind of mirror type system. I’m less happy about Union(NA, 
>> T)
>> >> being a super type of T, even though there are some good reasons that 
>> you’d
>> >> like to view T as a specialization of Union(NA, T). But I agree that I 
>> don’t
>> >> have a good feel about where missing data belongs in the type 
>> hierarchy.
>> >> This is another question for which I’d love to get input from others.
>> >>
>> >>
>> >>
>> >> In regard to Simon’s performance points:
>> >>
>> >>
>> >>
>> >> * Yes, memory usage alone argues strongly for working with DataArray{T}
>> >> rather than Array{Option{T}}.
>> >>
>> >>
>> >>
>> >> * Exploting tricks that make operations like anyna() faster is another
>> >> good argument for keeping DataA

Re: [julia-users] bit wrangling advice

[Simon Kornblith]

Assuming you have enough memory to write a BitArray to the JLD file 
initially, if you later open the JLD file with mmaparrays=true and read it, 
JLD will mmap the underlying Vector{Uint64} so that pieces are read from 
the disk as they are accessed. (The actual specifics of how this works is 
up to the OS, but generally it works well.) In principle you can also 
modify the BitArray the changes will be saved to the disk, although I'm not 
sure how well that works since I don't do it in my own code. There is no 
easy way to resize the BitArray if you do this, though.

Simon

On Tuesday, August 5, 2014 5:06:16 PM UTC-4, Tim Holy wrote:
>
> To me it sounds like you've come up with the main options: BitArray or 
> Array{Bool}. Since a BitArray is, underneath, a Vector{Uint64} with 
> different 
> indexing semantics, it seems you could probably come up with a reasonable 
> way 
> to update just part of it. But even if you use Array{Bool}, you're "only" 
> talking a few hundred megabytes, which is not a catastrophically large. 
> Also 
> consider keeping everything in memory; with 100GB of RAM you could store 
> an 
> awful lot of selections. 
>
> --Tim 
>
> On Tuesday, August 05, 2014 12:01:58 PM g wrote: 
> > Hello, 
> > 
> > I have an application where I have a few hundred million events, and I'd 
> > like to make and work with different selections of sets of those events. 
> > The events each have various values associated with them, say for 
> > simplicity color, timestamp, and loudness. Say one selection includes 
> all 
> > the events within 5 minutes after a blue event.  Or I want to select all 
> > events that aren't above some loudness threshold. I'd like to be able to 
> > save these selections in a JLD file for later use on some or all events. 
> I 
> > also need to be able update the selections as I observe more events. 
> > 
> > My baseline plane it to have an integer associated with each event and 
> each 
> > bit in the integer i corresponds to a selection.  So bit 1 is true for 
> > events within 5 minutes and bit 2 is true for events above the loudness 
> > threshold.  Then for each event's integer I can do bits(i)[1] or 
> bits(i)[2] 
> > to figure out if it is included in each selection. That seems like it 
> would 
> > be inefficient since bits() returns a string, so I would have to call 
> > bool(bits(i)[1]).  I could use a bitwise mask of some sort like 1&i==0 
> for 
> > the first bit and 2&i==0 for the second bit. 
> > 
> > A BitArray seems like a decent choice, except that you can only 
> read/write 
> > the entire array from a JLD file, rather than just a part of it.  That 
> will 
> > be inefficient since I'll often want to look at only a small subset of 
> the 
> > total events. And every time I want to update for new events, I would 
> need 
> > to read the entire BitArray, extend it in memory, then delete the old 
> JLD 
> > object and replace it with a new JLD object.  It seems plausible I could 
> > figure out how to read/write part of a BitArray from a JLD as I've 
> already 
> > done some hacking on HDF5.jl, but that could be a large amount of work. 
> > 
> > An Array{Bool} works well with JLD, and seems just as well suited as a 
> > BitArray.  I think it's 8 times bigger than BitArray, and has a similar 
> > space ratio to an integer (depending on how many selections I actually 
> use) 
> > because bools are stored as 1 byte? I can probably live with that, 
> although 
> > again it seems sort of inefficient. 
> > 
> > Any advice on how I should go about deciding, or options I hadn't 
> > considered?  Also why does bits() return a string, instead of say 
> > Vector{Bool} or BitArray? 
>
>

[julia-users] Re: Ordinary least squares regression

[Simon Kornblith]

This doesn't seem quite right: assuming the model has an intercept, SStot 
is traditionally the sum of squares of the intercept only model (i.e., 
sumabs2(y 
- mean(y)). You can see this if you add a constant to the first column of dd, 
which should not change R^2 but instead results in an implausibly high 
value for random data. Two one-liners that should give correct results are:

r2 = 1-sumabs2(residuals(mod))/sumabs2(y - mean(y)) # or 
1-var(residuals(mod))/var(y)
adjr2 = 1-scale(mod, true)/var(y)

where y is the dependent variable and mod is the model. If the model is fit 
using a DataFrame, at least for now, adjr2 needs to be:

adjr2 = 1-scale(mod.model, true)/var(y)

Simon

On Wednesday, August 6, 2014 11:34:54 PM UTC-4, Taylor Maxwell wrote:
>
> I haven't seen code to do it yet but it is very simple to calculate with a 
> LinearModel from GLM  Below is some code to calculate r-squared or adjusted 
> r-squared from a linear model in GLM calculated from a dataframe or from a 
> model calculated without a data frame.  At the bottom is some simple code 
> to apply these functions.  It uses the Abs2Fun() from the NumericFuns 
> package to calculate the sum of squares quickly.
>
> using DataFrames
> using GLM
>
> import DataFrames.DataFrameRegressionModel
> import GLM.LinearModel
> import GLM.residuals
> using NumericFuns
>
>
> function 
> rsquared(mod::DataFrameRegressionModel{LinearModel{DensePredQR{Float64}},Float64})
> SStot=sum(Abs2Fun(),mod.model.rr.y)
> SSres=sum(Abs2Fun(),residuals(mod))
> return (1-(SSres/SStot))
> end
>
> function rsquared(mod::LinearModel{DensePredQR{Float64}})
> SStot=sum(Abs2Fun(),mod.rr.y)
> SSres=sum(Abs2Fun(),residuals(mod))
> return (1-(SSres/SStot))
> end
>
> function 
> adjrsquared(mod::DataFrameRegressionModel{LinearModel{DensePredQR{Float64}},Float64})
> SStot=sum(Abs2Fun(),mod.model.rr.y)
> SSres=sum(Abs2Fun(),residuals(mod))
> n=size(mod.model.rr.y,1)  #number of samples
> p=size(mod.mm.m,2)-1  #number of variables besides constant (assumes 
> intercept)
> return 1- ( (SSres/(n-p-1)) / (SStot/(n-1)) )
> end
>
> function adjrsquared(mod::LinearModel{DensePredQR{Float64}})
> SStot=sum(Abs2Fun(),mod.rr.y)
> SSres=sum(Abs2Fun(),residuals(mod))
> n=size(mod.rr.y,1)#number of samples
> p=size(mod.pp.X,2)-1  #number of variables besides constant (assumes 
> intercept)
> return 1- ( (SSres/(n-p-1)) / (SStot/(n-1)) )
> end
>
>
> dd=randn(1000,7);
> df=convert(DataFrame,dd);
> aa=fit(LinearModel,x1~x2+x3+x4+x5+x6,df)
> bb=fit(LinearModel,aa.mm.m,aa.model.rr.y)
>
> rsquared(aa)
> adjrsquared(aa)
>
> rsquared(bb)
> adjrsquared(bb)
>
>
> On Wednesday, August 6, 2014 5:33:15 PM UTC-6, Jason Solack wrote:
>>
>> Hello everyone.  I am looking for a way to run this regression and 
>> retrieve the corresponding r squared value.  Is there a package out there 
>> that does this?
>>
>> Thank you for the help!
>>
>> Jason
>>
>

Re: [julia-users] How to find the indices of the unique rows in a matrix

[Simon Kornblith]

unique with a dim argument actually computes this as byproduct but does not 
return it. All we need is an API.

On Monday, August 11, 2014 1:27:57 PM UTC-4, Jacob Quinn wrote:
>
> There's an open issue about it here: 
> https://github.com/JuliaLang/julia/issues/1845
>
> I've also played around with a simple version based on a Set:
>
> function uniqueind(C)
> out = Int[]
> seen = Set{eltype(C)}()
> for i = 1:length(C)
> @inbounds x = C[i]
> if !in(x, seen)
> push!(seen, x)
> push!(out, i)
> end
> end
> out
> end
>
>
> On Mon, Aug 11, 2014 at 1:22 PM, Peter Simon  > wrote:
>
>> Greetings.
>>
>> I didn't find any built-in function to return the indices for the unique 
>> rows in a matrix, so I tried rolling my own.  I looked at, but didn't 
>> understand, the highly macro-ized, arbitrary dimensional code in the 
>> built-in unique() function.  I did discover hashing, though, from looking 
>> at that code.  Here is what I have so far:
>>
>>  function uniquerowindices2{T}(a::Matrix{T})
>>  rowhash = Uint64[hash(a[k,:]) for k = 1:size(a,1)]
>>  hu = unique(rowhash)
>>  ind = zeros(Int,length(hu))
>>  k1 = 1
>>  for m = 1:length(ind)
>>  for k = k1:length(rowhash)
>>  if rowhash[k] == hu[m]
>>  ind[m] = k
>>  k1 = k
>>  break
>>  end
>>  end
>>  end
>>  return ind
>>  end
>>
>> This is quite a bit faster than my first effort, where I was comparing 
>> rows, rather than row hashes, but I'm betting someone out there has a much 
>> better way.  Thanks in advance,
>>
>> --Peter
>>
>
>

Re: [julia-users] How to find the indices of the unique rows in a matrix

[Simon Kornblith]

Here it is copy/pasted from Base, in case it's useful for you for now:

using Base.Cartesian, Base.Prehashed
@ngenerate N typeof(A) function indunique{T,N}(A::AbstractArray{T,N}, dim::
Int)
1 <= dim <= N || return copy(A)
hashes = zeros(Uint, size(A, dim))

# Compute hash for each row
k = 0
@nloops N i A d->(if d == dim; k = i_d; end) begin
   @inbounds hashes[k] = hash(hashes[k], hash((@nref N A i)))
end

# Collect index of first row for each hash
uniquerow = Array(Int, size(A, dim))
firstrow = Dict{Prehashed,Int}()
for k = 1:size(A, dim)
uniquerow[k] = get!(firstrow, Prehashed(hashes[k]), k)
end
uniquerows = collect(values(firstrow))

# Check for collisions
collided = falses(size(A, dim))
@inbounds begin
@nloops N i A d->(if d == dim
  k = i_d
  j_d = uniquerow[k]
  else
  j_d = i_d
  end) begin
if (@nref N A j) != (@nref N A i)
collided[k] = true
end
end
end

if any(collided)
nowcollided = BitArray(size(A, dim))
while any(collided)
# Collect index of first row for each collided hash
empty!(firstrow)
for j = 1:size(A, dim)
collided[j] || continue
uniquerow[j] = get!(firstrow, Prehashed(hashes[j]), j)
end
for v in values(firstrow)
push!(uniquerows, v)
end

# Check for collisions
fill!(nowcollided, false)
@nloops N i A d->begin
 if d == dim
 k = i_d
 j_d = uniquerow[k]
 (!collided[k] || j_d == k) && continue
 else
 j_d = i_d
 end
 end begin
if (@nref N A j) != (@nref N A i)
nowcollided[k] = true
end
end
(collided, nowcollided) = (nowcollided, collided)
end
end

sort!(uniquerows)
end

On Monday, August 11, 2014 2:43:24 PM UTC-4, Simon Kornblith wrote:
>
> unique with a dim argument actually computes this as byproduct but does 
> not return it. All we need is an API.
>
> On Monday, August 11, 2014 1:27:57 PM UTC-4, Jacob Quinn wrote:
>>
>> There's an open issue about it here: 
>> https://github.com/JuliaLang/julia/issues/1845
>>
>> I've also played around with a simple version based on a Set:
>>
>> function uniqueind(C)
>> out = Int[]
>> seen = Set{eltype(C)}()
>> for i = 1:length(C)
>> @inbounds x = C[i]
>> if !in(x, seen)
>> push!(seen, x)
>> push!(out, i)
>> end
>> end
>> out
>> end
>>
>>
>> On Mon, Aug 11, 2014 at 1:22 PM, Peter Simon  wrote:
>>
>>> Greetings.
>>>
>>> I didn't find any built-in function to return the indices for the unique 
>>> rows in a matrix, so I tried rolling my own.  I looked at, but didn't 
>>> understand, the highly macro-ized, arbitrary dimensional code in the 
>>> built-in unique() function.  I did discover hashing, though, from looking 
>>> at that code.  Here is what I have so far:
>>>
>>>  function uniquerowindices2{T}(a::Matrix{T})
>>>  rowhash = Uint64[hash(a[k,:]) for k = 1:size(a,1)]
>>>  hu = unique(rowhash)
>>>  ind = zeros(Int,length(hu))
>>>  k1 = 1
>>>  for m = 1:length(ind)
>>>  for k = k1:length(rowhash)
>>>  if rowhash[k] == hu[m]
>>>  ind[m] = k
>>>  k1 = k
>>>  break
>>>  end
>>>  end
>>>  end
>>>  return ind
>>>  end
>>>
>>> This is quite a bit faster than my first effort, where I was comparing 
>>> rows, rather than row hashes, but I'm betting someone out there has a much 
>>> better way.  Thanks in advance,
>>>
>>> --Peter
>>>
>>
>>

[julia-users] Re: Specifying Julia commit in REQUIRE file

[Simon Kornblith]

Does it still not work to use 0.4.0-dev+n as the version in the REQUIRE 
file? This used to almost work, but some of the nightlies were missing the 
commit number. It certainly seems easier than searching through all the 
hashes, although I don't know the git command to get the commit number for 
a given hash off the top of my head.

Simon

On Wednesday, August 20, 2014 1:36:55 PM UTC-4, Ivar Nesje wrote:
>
> That would only require storing ~2 sha1 hashes and search through, so 
> it is definitely computationally possible. If we limit the selection to all 
> commits since last tag, we will only have a few thousand.
>
> The big question is whether anyone will actually use this functionality 
> correctly. Package devs are usually not very interested in such tedious 
> work.
>
> You can use the commit number to warn users that the new version of your 
> package require a new version of Julia.
>
> if BASE.GIT_VERSION_INFO.build_number < 60
>error("please update Julia ") 
> end
>
>

Re: [julia-users] Re: Specifying Julia commit in REQUIRE file

[Simon Kornblith]

How did you test it? On my system it seems to work: I can get Pkg to 
upgrade a package if the REQUIRE file specifies julia 0.4.0-dev+274 but not 
if it specifies julia 0.4.0-dev+275. (I'm on 0.4.0-dev+274.)

Simon

On Wednesday, August 20, 2014 4:29:36 PM UTC-4, Elliot Saba wrote:
>
> @Simon I just tried it out and on my machine at least, it didn't work. 
>  The REQUIRE file doesn't honor the + symbol for some reason, so things 
> like 0.3.0-rc4 work, but things like 0.4.0-dev+158 don't work.  (They are 
> treated the same as 0.4.0-dev I believe)
> -E
>
>
> On Wed, Aug 20, 2014 at 3:32 PM, Júlio Hoffimann  > wrote:
>
>> What do you recommend for my use case?
>>
>> Could you please clarify the meaning of all the following version numbers?
>>
>> 0.3-
>> 0.3
>> 0.3+
>>
>> Thanks,
>> Júlio.
>>
>
>

[julia-users] Re: Specifying Julia commit in REQUIRE file

[Simon Kornblith]

>From that script it looks like it's

git rev-list  ^v0.3.0-rc3 | wc -l

On Wednesday, August 20, 2014 9:22:36 PM UTC-4, Gray Calhoun wrote:
>
> On Wednesday, August 20, 2014 12:56:42 PM UTC-5, Simon Kornblith wrote:
>>
>> Does it still not work to use 0.4.0-dev+n as the version in the REQUIRE 
>> file? This used to almost work, but some of the nightlies were missing the 
>> commit number. It certainly seems easier than searching through all the 
>> hashes, although I don't know the git command to get the commit number for 
>> a given hash off the top of my head.
>>
>
> It's not a simple git command. It looks like base/version_git.sh has a 
> script to get the commit number, which is then stored in 
> base/version_git.jl.
>

Re: [julia-users] Re: Can this magic square function be further optimized?

[Simon Kornblith]

Here's an interesting comparison:

https://gist.github.com/simonster/6195af68c6df33ca965d

idiv for 64-bit integers is one of the most expensive extant x86-64 
instructions. For 32-bit integers, it is much cheaper, and this function 
runs nearly twice as fast. When LLVM knows the divisor in advance, it can 
avoid idiv and perform the division as a multiplication and bit shift, 
which is faster still. But the biggest gain comes from avoiding rem in the 
inner loop entirely.

Of course, this is only tackling the odd case, and additional performance 
can be gained for all of these benchmarks by avoiding allocation.

Simon

On Tuesday, August 26, 2014 9:11:47 AM UTC-4, Phillip Berndt wrote:
>
>
> (1) allocate the output M outside of the core algorithm, and pass it as an 
>> input, i.e., 
>>
>
> I did that, though it can be argued that this is cheating given that the 
> competitors also have to allocate an array for each loop. With that version 
> (and some more slight optimization: Storing intermediate values in the for 
> loops, using column-major indexing and @simd) [
> https://gist.github.com/phillipberndt/7dc0aed7eb855f900f0d/21cce76664bdc59f6203ff6f3496e80e256f54cb],
>  
> the overall time for the N=3..1000 test case is down to 3.67s.
>
> (2) @time (for i = 1:100; magic!(M); end). Did it allocate any memory? 
>> Then 
>> you have a problem. Use the profiler, or run julia with --track- 
>> allocation=user, to find out where it occurs.  
>
>
> It does, about 3 Mb on line 2 (if n % 2 == 1). Doesn't make much sense so 
> I guess the profiler interfered with the optimizer here?! I doubt that 
> trying to get rid of the 3Mb will gain another second though.
>  
>
>> (3) Even if it's not allocating, you may have a bottleneck. Use the 
>> profiler to 
>> find it. 
>>
>
> The line where the most time is spent is line 11, filling the array in the 
> odd case. I don't see how it could be optimized any further, so that's 
> probably as far as one gets?!
>
> - Phillip
>

Re: [julia-users] Re: Can this magic square function be further optimized?

[Simon Kornblith]

See my answer. If you don't mind a copy/pasting a bunch of bit twiddly code 
from Base, this can be efficiently devectorized, and that turns out to be a 
decent perf win. But in this case the devectorized code is quite difficult 
to decipher and Stefan Schwarz's point definitely applies.

Simon

On Tuesday, August 26, 2014 4:42:59 PM UTC-4, Iain Dunning wrote:
>
> To muddy the water in this now somewhat confusing email thread: 
> devectorizing/avoiding memory allocation isn't always a win:
>
>
> http://stackoverflow.com/questions/25412323/how-can-i-do-a-bitwise-or-reduction-along-an-axis-of-a-boolean-array-in-julia/25414595#25414595
>
>
>
> On Tuesday, August 26, 2014 4:12:49 PM UTC-4, gael@gmail.com wrote:
>>
>> >I'm a little surprised that you have found the performance implications 
>> unconvincing in discussions where the OGs advocate devectorization.
>>
>> A 10x speed-up on a 10 ms calculation is generally considered 
>> unconvincing. An unknown gain on an unprofiled code for an undefined 
>> context is definitely unconvincing, especially if the vectorized form is 
>> widely clearer. 
>>
>> I'm not talking about this specific case for which I find explicit loops 
>> at least as clear. But I, too, sometimes find the devectorization advice 
>> not convincing. But fortunately, this is usually for small scripts, not for 
>> actual Julia libs.
>>
>>

[julia-users] Re: Performance issue with @eval and outer constructor?

[Simon Kornblith]

It seems to have to do with the way the unparametrized type gets 
interpolated into the AST. Changing

for (fp, fpc) in [(DoublelengthFloat{Float64}, DoublelengthFloat)]


to

for (fp, fpc) in [(DoublelengthFloat{Float64}, :DoublelengthFloat)]


on line 108 makes g and h perform equivalently for me.

Simon

On Tuesday, August 26, 2014 7:22:25 PM UTC-4, Jiahao Chen wrote:
>
> Sorry for the vague title, but I'm not really sure what's going on in this 
> code. 
>
> This gist 
>
> https://gist.github.com/jiahao/f2e20f2b42118ae06b79 
>
> contains an implementation of emulated double double precision 
> arithmetic (Dekker, 1971; 
> http://link.springer.com/article/10.1007%2FBF01397083) 
>
> The last part of the gist starting at line 94 runs a simple benchmark 
> on a typical arithmetic problem. Three versions of the same kernel are 
> defined: 
>
> 1. f contains a manually written kernel calling the outer constructor 
> (on line 8) directly 
> 2. g is generated at compile time using @eval and calls the same outer 
> constructor 
> 3. h is generated at compile time using @eval and calls `convert` (on 
> line 9) which calls the same outer constructor. 
>
> The resulting runtimes are: 
>
> elapsed time: 0.80593951 seconds (13904 bytes allocated) 
> elapsed time: 2.402564403 seconds (120024 bytes allocated, 24.93% gc 
> time) 
> elapsed time: 0.74216029 seconds (24 bytes allocated) 
>
> It's not clear to me why there is a big performance gap between g and 
> the other versions. At first I thought the constructor might not be 
> inlined in g, resulting in type instability during the summation. 
> However, the annotations from code_typed on the variable thesum are 
> the same, i.e. Any. 
>

[julia-users] Re: A common idiom

[Simon Kornblith]

It doesn't look like the first snippet is inlining hash(p, zero(Uint)) to 
me. Rather, it looks like that code isn't specialized on the input 
argument. There's still only one call outside of the call to jl_box_uint64 
and no other work being done, so it can't really be calling hash twice, and 
it's making the call via jl_apply_generic rather than directly. I'm not 
sure why it wouldn't be specialized, but it seems likely to be an issue 
with @code_llvm rather than the compiler.

Simon

On Thursday, August 28, 2014 6:21:17 AM UTC-4, Steven G. Johnson wrote:
>
> Interesting.   Defining foo(x, y=0) is supposed to be equivalent to 
> defining foo(x, y) and foo(x)=foo(x,0), but it looks like the latter form 
> triggers inlining and the former does not. The latter form is more 
> common in older Julia code because the former syntax was only added in 
> Julia 0.2.
>
> However, I can't reproduce it on simpler examples; in most cases that I 
> try, the two definition forms produce the same code.
>

[julia-users] Re: "self" dot product of all columns in a matrix: Julia vs. Octave

[Simon Kornblith]

Actually, if you want this to be fast I don't think you can avoid the if y 
== 2 branch, although ideally it should be outside the loop and not inside 
it. LLVM will optimize x^2, but if it doesn't know what the y in x^y is at 
compile time, it will just call libm. It's not going to emit a branch to 
handle the case where y == 2, since it has no idea what arguments the 
function will actually be called with and it doesn't want to slow things 
down if it's not called with y == 2. It just so happens that X.^2 is pretty 
common in real life.

In the Octave source I see that Octave special-cases X.^2, X.^3, and 
X.^(-1). (Notice how X.^4 is ~6x slower.) We could do the same, or inline 
the loop so that LLVM's optimizers can get a crack at it in the constant 
power case, or dynamically generate a specialized loop that uses multiply 
operations for arbitrary integer powers (but see here 
 for 
a caveat: this will be faster than calling libm, but it will be less 
accurate). But if you want a fast version of X.^2, right now you can just 
call abs2(X).

Tony is right that sumabs2(X, 1) will be faster, both versus sum(abs2(X), 1) 
and versus anything you could achieve with Octave, since it avoids 
allocating a temporary array and exploits BLAS. (On my machine it's about 
9x faster than sum(X.^2, 1) in Octave.) But on my machine sum(abs2(X), 1) 
in Julia is only a few percent slower than sum(X.^2, 1) in Octave if that's 
all you want.

Simon

On Monday, September 8, 2014 12:37:16 PM UTC-4, Tony Kelman wrote:
>
> Calm down. It's very easy to match Octave here, by writing the standard 
> library in C. Julia chose not to go that route. Of course everyone would 
> like the concise vectorized syntax to be as fast as the for-loop syntax. If 
> it were easy to do that *without* writing the standard library in C, it 
> would've been done already. In the meantime it isn't easy and hasn't been 
> done yet, so when someone asks how to write code that's faster than Octave, 
> this is the answer.
>
> There have been various LLVM bugs regarding automatic optimizations of 
> converting small integer powers to multiplications. "if y == 2" code really 
> shouldn't be necessary here, LLVM is smarter than that but sometimes fickle 
> and buggy under the demands Julia puts on it.
>
> And in this case the answer is clearly a reduction for which Julia has a 
> very good set of tools to express efficiently. If your output is much less 
> data than your input, who cares how fast you can calculate temporaries that 
> you don't need to save?
>
>
> On Monday, September 8, 2014 9:18:02 AM UTC-7, Jeff Waller wrote:
>>
>>
>>
>> On Monday, September 8, 2014 7:37:52 AM UTC-4, Mohammed El-Beltagy wrote:
>>>
>>> Octave seems to be doing a pretty good job for this type of calculation. 
>>> If you want to squeeze a bit more performance out of Julia, you can try to 
>>> use explicit loops (devectorize as in 
>>> http://julialang.org/blog/2013/09/fast-numeric/). Then you might also 
>>> remove bounds checking in a loop for faster performance. 
>>>
>>> Try this function:
>>> function doCalc(x::Array{Float64,2})
>>> XX=Array(Float64, 7000,7000)
>>> for j=1:7000,i=1:7000
>>> @inbounds XX[i,j]=x[i,j]*x[i,j]
>>> end
>>> XX
>>> end
>>>
>>> Followed by 
>>> @time XX=doCalc(X);
>>>
>>
>>
>> I am totally not cool with this. Just simply making things fast at 
>> whatever the cost is not good enough. You can't tell someone that "Oh sure, 
>> Julia does support that extremely convenient syntax, but because it's 6 
>> times slower (and it is), you need to essentially write C".  There's no 
>> reason that Julia for this should be any less fast than Octave except for 
>> bugs which will eventually be fixed.  I think it's perfectly ok to say 
>> devectorize if you want to do even better, but it must be at least equal.
>>
>> Here's the implementation of .^
>>
>>
>> .^(x::StridedArray, y::Number) = reshape([ x[i] ^ y for i=1:length(x) ], 
>> size(x))
>> how is that made less vectorized?
>>
>> Yes R, Matlab, Octave do conflate vectorization for speed and style and 
>> that's a mistake.  Julia tries not to do that that's better.  But do you 
>> expect people to throw away the Matlab syntax?  No way.
>>
>> for me:
>>
>> *julia> **x = rand(7000,7000);*
>>
>> *julia> **@time x.^2;*
>>
>> elapsed time: 1.340131226 seconds (392000256 bytes allocated)
>>
>> And in Octave
>>
>> >> x=rand(7000);
>>
>> >> tic;y=x.^2;toc
>>
>> Elapsed time is 0.201613 seconds.
>>
>> Comprehensions might be the culprit, or simply use of the generic x^y 
>>  here's an alternate implementation I threw together
>>
>> *julia> **function .^{T}(x::Array{T},y::Number)*
>>
>>*   z = *(size(x)...)*
>>
>>*   new_x = Array(T,z);*
>>
>>*   tmp_x = reshape(x,z);*
>>
>>*   for i in 1:z*
>>
>>*  @inbounds if y == 2*
>>
>>*   

[julia-users] Re: cld not defined in Julia 0.3.0

[Simon Kornblith]

This was just added to Julia 0.4 yesterday, which is why it's not defined 
in Julia 0.3 :). In general if you're using Julia 0.3 you should be looking 
at the release-0.3 docs and not latest (click at the lower right on 
julia.readthedocs.org to change).

Simon

On Monday, September 8, 2014 4:59:28 PM UTC-4, curiou...@gmail.com wrote:
>
>
> Is it a bug that when I call *cld* in Julia 0.3.0 on Mac OSX it says:
>
> ERROR: cld not defined
>
> Thanks.
>

Re: [julia-users] Re: "self" dot product of all columns in a matrix: Julia vs. Octave

[Simon Kornblith]

I don't think there's anything wrong with specializing small powers. We 
already specialize small powers in openlibm's pow and even in inference.jl. 
If you want to compute x^2 at maximum speed, something, somewhere needs 
that when the power is 2 it should just call mulsd instead of doing 
whatever magic it does for arbitrary powers. Translating integer exponents 
to multiplies only kind of works for powers higher than 3 anyway; it will 
be faster but not as accurate as calling libm.

Another option would be to use a vector math library, which would make 
vectorized versions of everything faster. As far as I'm aware VML is still 
the only option that claims 1 ulp accuracy, which does indeed provide a 
substantial 
performance boost  (and 
also seems to special case x^2). I don't think we could actually use it in 
Base, though; even if the VML license allows it, I don't think we could 
link against Rmath, UMFPACK, and whatever other GPL libraries if we did.

Simon

On Monday, September 8, 2014 4:02:34 PM UTC-4, Stefan Karpinski wrote:
>
> On Mon, Sep 8, 2014 at 7:48 PM, Tony Kelman  > wrote:
>
>>
>> Expectations based on the way Octave/Matlab work are really not 
>> applicable to a system that works completely differently.
>>
>
> That said, we do want to make this faster – but we don't want to cheat to 
> do it, and special casing small powers is kind of cheating. In this case it 
> may be worth it, but what we'd really like is a more general kind of 
> optimization here that has the same effect in this particular case.
>

Re: [julia-users] type generation

[Simon Kornblith]

Yup. The reason this works when oVector is a Vector{Float64} is that Julia 
makes a copy in the process of converting it to a Vector{Int} when you 
construct the Cell.

Simon

On Thursday, September 11, 2014 12:53:52 AM UTC-4, John Myles White wrote:
>
> This sure looks like you're not making any copies when you seem to want 
> copies. 
>
> In particular, this line: 
>
> >   cellList[i] = Cell(i, oVector) 
>
> probably needs to be 
>
> >   cellList[i] = Cell(i, copy(oVector)) 
>
>  -- John 
>
> On Sep 10, 2014, at 6:41 PM, Andre Bieler  > wrote: 
>
> > can anyone tell me why the following code does not work as (I) expected? 
> > I have a simple type Cell which only has an index and a origin array. 
> > When creating multiple instances of this Cell type in a loop and 
> assigning 
> > different origin arrays to them, in the end all instances have the same 
> > origin array. I am using julia 0.3 
> > 
> > (It does work if the commented line is un-commented though..) 
> > 
> > 
> > type Cell 
> >   index::Int64 
> >   origin::Array{Int64,1} 
> > end 
> > 
> > nDim = int(3) 
> > nCells = int(2) 
> > 
> > cellList = Array(Cell, nCells) 
> > oVector = Array(Int64, nDim)   # does not work as expected 
> > #oVector = Array(Float64, nDim) # does work 
> > 
> > for i=1:nCells 
> >   for j=1:nDim 
> > oVector[j] = i 
> >   end 
> >   cellList[i] = Cell(i, oVector) 
> > end 
> > 
> > # after the loop, both cells have same origin... 
> > println(cellList[1].index, ": ", cellList[1].origin) 
> > println(cellList[2].index, ": ", cellList[2].origin) 
> > 
> > 
> > thanks in advance 
> > 
> > andre 
>
>

Re: [julia-users] Why do bools have functions like abs and sign defined?

[Simon Kornblith]

But then the question is why we define specialized versions for Bool rather 
than using the methods for Real in number.jl. For all but abs LLVM is smart 
enough to optimize the methods for Real to no-ops.

Simon

On Thursday, September 11, 2014 5:27:10 AM UTC-4, Tim Holy wrote:
>
> It's because 
>
> julia> super(Bool) 
> Integer 
>
> julia> super(Integer) 
> Real 
>
> and you'd expect these to be defined for Integers. If this weren't the 
> case 
> then 
>
> function myalgorithm{T<:Real}(x::AbstractArray{T}) 
> # do something involving abs 
> end 
>
> would fail if the user supplied an array of Bools. I could create a 
> special 
> version of myalgorithm for Array{Bool} that omitted calls to abs, but that 
> would be a pain in the neck, and I'd have to worry about every single 
> possible 
> element type ("does Uint16 have an abs method? let's check"). Since Julia 
> compiles a new version of myalgorithm specificially for Array{Bool}, 
> there's 
> also no performance hit (such trivial functions just get inlined out into 
> no- 
> ops). 
>
> --Tim 
>
> On Wednesday, September 10, 2014 10:20:42 PM Dan Luu wrote: 
> > In bool.jl, there are things like 
> > 
> > signbit(x::Bool) = false 
> > sign(x::Bool) = x 
> > abs(x::Bool) = x 
> > abs2(x::Bool) = x 
> > 
> > Are these because Bool is a subtype of something where you'd expect 
> > these to be defined, or are these useful in their own right, and if 
> > so, when would you use these? 
> > 
> > 
> > Thanks, 
> > Dan 
>
>

[julia-users] Re: Julia strange behaviour with GC

[Simon Kornblith]

You're right that microbenchmarks often do not reflect real-world 
performance, but in defense of using the sample mean for benchmarking, it's 
a good estimator of the population mean (provided the distribution has 
finite variance), and if you perform an operation n times, it will take 
about nμ seconds. If you want your code to run as fast as possible that's 
probably what matters.

The reason Julia's behavior is periodic and Octave's is not is almost 
certainly Octave uses reference counting, which incurs overhead whenever an 
array is allocated or destroyed, while Julia uses garbage collection, which 
incurs overhead only once a certain amount of memory has been allocated. 
People are certainly aware that Julia's garbage collector is suboptimal 
, and there is work being 
done  to improve it. As long 
as you run an operation enough times, the mean will accurately reflect 
garbage collection overhead. (That isn't just by chance.) But the frequency 
of garbage collection will depend not only on the amount of memory 
allocated, but also on the total amount of memory the system has, so that 
is not necessarily useful unless you benchmark on the same machine.

Simon

On Friday, September 12, 2014 11:11:13 AM UTC-4, gael@gmail.com wrote:
>
> A few days ago, Ján Dolinský posted some microbenchmark 
> s 
> about Julia being slower than Octave for squaring quite large matrices. At 
> that time, I wanted to answer with the traditional *microbenchmarks are 
> evil*™ type of answer.
>
> I wanted to emphases that:
>
>1. The mean is hardly a good distribution location estimator in the 
>general (non-normal) case.
>2. Neglecting that, what he calculated was a typical execution time 
>for the squaring of *one and only one* random matrix filled with 
>7000x7000 random Float64 numbers with his computer.
>3. Assuming the 7000x7000 matrix size represents the typical matrix 
>size for his work, that means that Julia is likely *fast enough*: it's 
>bellow 1 s, the difference between Julia and Octave is not noticeable by 
>anything except a computer.
>
>
> But from his latter concerns, I guess this was not really what he had in 
> mind (but I may have misinterpreted his words). I supposed that what he 
> really wanted was to calculate *many *of such squared matrices, so many 
> of them that the observed difference in execution times would matter a lot.
>
>
> What would be a not-too-bad benchmark adapted to this use case? I guess 
> repeating the calculation is required. But instead of measuring time 
> outside of the loop, it might be very interesting to get it inside of the 
> loop: for each one of the squarings. These individual measurements, if not 
> too short, would allow us:
>
>- to plot the complete time distribution instead of assuming normality 
>and providing only one distribution parameter;
>- to get the cumulative distribution to represent the case of multiple 
>consecutive matrix squarings;
>- to follow the execution time evolution through time or iterations.
>
> So here is the snippet of code that I used:
>
> t = Float64[]
> for i = 1:2020
> X = rand(1000,1000)
> tic()
> XX = X.^2
> b = toq()
> push!(t, b)
> end
>
> The matrix generation is in-loop to avoid any smart caching.
>
>
> *1- Distributions*
>
> So here are the distributions (mind the Y axis is log):
>
>
> 
>
> Octave execution times are small and very close from each other. The 
> behavior of those two Julia repetitions of the same calculation are much 
> weirder : there is not just one peak, but 3 or 4. That's a strong hint 
> *against* mean calculations. Furthermore, complete distributions contain 
> much more information. I really think that calculating such estimations of 
> the distributions is the bare minimum when benchmarking.
>
> Ooops, I just said it: *estimations*. The methods used to estimate 
> distributions usually need arbitrary parameters to be set (here, the 
> Gaussian kernel and its bandwidth). So the peaks above are of almost 
> arbitrary width. Who do I know whether these peaks represent something or 
> are just produced by poor bandwidth choice? Repetitions of the whole 
> process is one good way to test that. But here, I'll take another approach: 
> plotting execution time vs. the iteration numbers.
>
> *2- Execution times vs. iteration number*
>
> So let's look at these execution times:
>
>
> 
>
> Ohoh! (Not so) Surprising! While Octave seems to behave as expected 
> (random distribution around a single value), Julia exhibits a completely 

[julia-users] Re: Julia strange behaviour with GC

[Simon Kornblith]

In the world of technical computing, we are mostly interested in running 
the same computationally intensive functions over and over, and we usually 
only care about how long that will take. The best way to approximate how 
long that will take is to run the function over and over and see how long 
it takes. Moreover, the central limit theorem guarantees that, given a 
function that takes a mean time μ to run, the amount of time it takes to 
run that function n times is approximately normally distributed with mean 
nμ for sufficiently large n regardless of the distribution of the 
individual function evaluation times as long as it has finite variance. Of 
course if something takes longer to run the more times you run it the 
central limit theorem will not apply, but then a histogram won't be much 
use then either, although plotting function eval time vs. clock time would 
be.

Your effect from running the function many times is curious but it's not 
reproducible for me. After running for 20 minutes I still just see two 
bands. The difference between the means of the 2020 samples at the 
beginning, middle, and end is <2%. Can you reproduce that?

Simon

On Friday, September 12, 2014 3:03:49 PM UTC-4, gael@gmail.com wrote:
>
> You're right that microbenchmarks often do not reflect real-world 
>> performance, but in defense of using the sample mean for benchmarking, it's 
>> a good estimator of the population mean (provided the distribution has 
>> finite variance), and if you perform an operation n times, it will take 
>> about nμ seconds. If you want your code to run as fast as possible that's 
>> probably what matters.
>>
> I don't agree. In general, whole programs with I/O has approximate 
> log-normal distribution (from my experience). In which case the geometric 
> mean is a far better estimator of the location and the mode of the 
> distribution.
>
> If you want your code to run as fast as possible, you should profile it. 
>
> If you want to profile code snippets, to get the best local performances, 
> then, and it's my whole point, you should at least plot the distributions.
>
> Means can be deceptive. A quick example:
>
>
> 
> Imagine the experimental data are obtained from function `f1` with a mean 
> of 550 s (asymmetric distribution). 
> Imagine the normal data are obtained from function `f2` with a mean of 490 
> s.
>
> You would choose `f2` without any doubt. Now, Stefan, who has the same 
> computer as you but invested in an SSD found that you made a poor choice in 
> choosing function `f2` because on his computer, function `f1` is normal and 
> gives a mean of 460 s vs. 490 s for `f1`. Who's correct?
>
> (We could prove that the tail was produced by inefficient IO)
>
>  
>
>> As long as you run an operation enough times, the mean will accurately 
>> reflect garbage collection overhead. (That isn't just by chance.) 
>>
> What was by chance is the almost linear trend which I proved wrong by 
> sampling values a while longer.
>
>

Re: [julia-users] slow julia version of c code

[Simon Kornblith]

There is most certainly a type problem. You're not getting type information 
for sparse_grid.lvl_l which deoptimizes a lot of things. In your code, you 
have:

type sparse_gridd::Int64q::Int64n::Int64grid
::Array{Float64}ind::Array{Int64}lvl::Array{Int64}
lvl_l::Array{Int64}lvl_s::Array{Float64}bounds::Array{
Float64}gtypeBBfubend

Try specifying the number of dimensions for the arrays, e.g. 
grid::Array{Float64,2}. Also using any of the untyped fields will slow code 
down, but I don't see them in that function.

A hackish way to find type issues is:

filter!(x->!isleaftype(x[2]), @code_typed(f(args))[1].args[2][2])

which returns a list of variables where type inference couldn't infer a 
concrete/leaf type. In some cases these will be "hidden" variables so you 
may still need to inspect the output of code_typed. TypeCheck.jl and 
Lint.jl may also provide more sophisticated tools for this purpose.

Simon

On Monday, September 15, 2014 5:16:36 AM UTC-4, Zac wrote:
>
> https://gist.github.com/Zac12345/519bd7a503a1fd1b8d98 has the updated 
> function and code_typed output
>
> Pkg.clone("https://github.com/Zac12345/Sparse";) should be fine for the 
> test as the gist does'nt require functions using the shared lib.
>
> Are there any specific things i should be looking for in the code_typed 
> results?
>

Re: [julia-users] Re: Control system library for Julia?

[Simon Kornblith]

Manually managing memory is far easier than implementing a hard real-time 
garbage collector, but hard real-time garbage collectors do exist: the Staccato 
paper 
 
claims worst-case pause times of under 1 ms.

Simon

On Monday, September 29, 2014 4:55:41 AM UTC-4, Steven Sagaert wrote:
>
> Hi Uwe, 
> In the Java world, Azul Zing is something along the lines of the first ref 
> but I think it is still only soft real time. If you want that kind of 
> garantees it's probably better not to keep garbage lying around to be 
> collected later but free the memory immediately when no longer needed like 
> via smart pointers as in C++. In fact that gives you almost everything a GC 
> does except for collecting cycles but circular datastructures are rare and 
> you can always break the cycles yourself so this isn't a big deal in 
> practice.
>
> On Sunday, September 28, 2014 9:10:39 PM UTC+2, Uwe Fechner wrote:
>>
>> Event though pre-allocation is preferred for real-time
>> application, often there are some parts of the code
>> that are hard to implement without any garbage
>> collection.
>>
>> Hard real-time is possible even if you are using a garbage 
>> collector.
>>
>> Some references for this point of view:
>>
>> http://michaelrbernste.in/2013/06/03/real-time-garbage-collection-is-real.html
>>
>> https://lwn.net/images/conf/rtlws-2011/proc/Klotzbuecher.pdf
>>
>> Best regards:
>>
>> Uwe Fechner
>>
>> On Sunday, September 28, 2014 4:53:48 PM UTC+2, Steven Sagaert wrote:
>>>
>>> GC will always be non-deterministic. For "hard" real time you just need 
>>> to manage memory yourself. That's the approach used by real time Java 
>>> http://www.rtsj.org/
>>>
>>> On Monday, September 15, 2014 10:25:07 AM UTC+2, Uwe Fechner wrote:

 Hi,
 I am working an airborne wind energy as well. I wrote a kite-power 
 system simulator in Python, where also one of the controllers (the winch 
 controller) is 
 implemented in Python. ( https://bitbucket.org/ufechner/freekitesim )

 With Python you can reach a jitter of less than 2 ms in the 20Hz 
 control loop quite easily (on low-latency Linux). In my case this is 
 sufficient for prototyping,
 but the real flight control system should run at a higher update 
 frequency (perhaps 200 Hz).

 In contrast to Julia Python is using reference counting, and in my 
 Python applications I just turn off the garbage collection.

 For Julia (and I would love to rewrite the simulator in Julia) this is 
 probably not an option. A better garbage collector 
 (which is in the pipeline, see: 
 https://github.com/JuliaLang/julia/pull/5227 ) would definitely help. 

 Generating embedded controllers in LLVM IR would be great!

 Best regards:

 Uwe

 On Monday, September 15, 2014 8:32:02 AM UTC+2, Andrew Wagner wrote:
>
> Hi Spencer!
>
> My job in airborne wind energy is ending soon so I don't have a 
> specific application (aside from control), but I would want to stay 
> sub-ms 
> for anything in-process.  I have been using Orocos extensively for the 
> last 
> few years.  It's the best control middleware in the open source world, 
> but 
> I think a lot of things could be improved if it was re-implemented in a 
> language with a better typesystem and introspection... one example would 
> be 
> that adding a new type to the system requires quite a bit of boilerplate 
> code, creating an incentive for users to just pass data in flat arrays, 
> subverting type safety.  
>
> Cheers,
> Andrew
>
> On Mon, Sep 15, 2014 at 7:03 AM, Spencer Russell <
> spencer@gmail.com> wrote:
>
>> Hi Andrew,
>>
>> What are your realtime deadlines? I'm working on live audio 
>> processing stuff with Julia, where I'd like to get the audio latency 
>> down 
>> into a few ms. Julia definitely isn't there yet (and might never get 
>> true 
>> hard-realtime), but there's some promising work being done on the GC to 
>> reduce pause time for lower-latency applications. It's also helpful to 
>> profile the code to reduce allocations (and the need for GC) down to a 
>> minimum. I haven't yet gotten down to zero-allocation code in my render 
>> loop, but once I got it down below 100 bytes I moved on to other more 
>> pressing features. At some point I'll dig deeper to see if I can get rid 
>> of 
>> the last few allocations.
>>
>> I'd definitely be happy if there are some more folks out there 
>> driving demand for lower-latency Julia. :)
>>
>> peace,
>> s
>>  
>> On Sun, Sep 14, 2014 at 3:58 PM, Andrew Wagner  
>> wrote:
>>
>>> Hello again Uwe!  
>>>
>>> It's fun running into someone I know on a language geek forum :) 
>>>  I'm helping one of our bac

Re: [julia-users] Re: Control system library for Julia?

[Simon Kornblith]

A collector with a hard constraint on the worst-case GC time is a hard 
real-time collector. I don't think many collectors can guarantee 2-4 ms 
worst case GC times (although evidently Staccato can). There is also a 
tradeoff between having a real-time collector and having the fastest 
possible collector for non-real-time applications that has been discussed 
elsewhere.

Simon

On Monday, September 29, 2014 1:37:13 PM UTC-4, Uwe Fechner wrote:
>
> For me it would be important to reach the real-time performance of Python/ 
> Numba with Julia:
>
> I would like to port my flight simulator (
> https://bitbucket.org/ufechner/freekitesim) to Julia,
> to improve the performance, make maintenance easier and to start working 
> with automated
> differentiation.
>
> I think this performance (max. 2-4 ms time for a call of an incremental 
> garbage collector) should
> be possible to achieve by implementing an improved garbage collector for 
> Julia.
>
> This would not be hard real-time, but already very useful for many 
> practical applications.
>
> Uwe
>
>
> On Monday, September 29, 2014 4:03:43 PM UTC+2, Simon Kornblith wrote:
>>
>> Manually managing memory is far easier than implementing a hard real-time 
>> garbage collector, but hard real-time garbage collectors do exist: the 
>> Staccato 
>> paper 
>> <http://researcher.watson.ibm.com/researcher/files/us-groved/rc24504.pdf> 
>> claims worst-case pause times of under 1 ms.
>>
>> Simon
>>
>> On Monday, September 29, 2014 4:55:41 AM UTC-4, Steven Sagaert wrote:
>>>
>>> Hi Uwe, 
>>> In the Java world, Azul Zing is something along the lines of the first 
>>> ref but I think it is still only soft real time. If you want that kind of 
>>> garantees it's probably better not to keep garbage lying around to be 
>>> collected later but free the memory immediately when no longer needed like 
>>> via smart pointers as in C++. In fact that gives you almost everything a GC 
>>> does except for collecting cycles but circular datastructures are rare and 
>>> you can always break the cycles yourself so this isn't a big deal in 
>>> practice.
>>>
>>> On Sunday, September 28, 2014 9:10:39 PM UTC+2, Uwe Fechner wrote:
>>>>
>>>> Event though pre-allocation is preferred for real-time
>>>> application, often there are some parts of the code
>>>> that are hard to implement without any garbage
>>>> collection.
>>>>
>>>> Hard real-time is possible even if you are using a garbage 
>>>> collector.
>>>>
>>>> Some references for this point of view:
>>>>
>>>> http://michaelrbernste.in/2013/06/03/real-time-garbage-collection-is-real.html
>>>>
>>>> https://lwn.net/images/conf/rtlws-2011/proc/Klotzbuecher.pdf
>>>>
>>>> Best regards:
>>>>
>>>> Uwe Fechner
>>>>
>>>> On Sunday, September 28, 2014 4:53:48 PM UTC+2, Steven Sagaert wrote:
>>>>>
>>>>> GC will always be non-deterministic. For "hard" real time you just 
>>>>> need to manage memory yourself. That's the approach used by real time 
>>>>> Java 
>>>>> http://www.rtsj.org/
>>>>>
>>>>> On Monday, September 15, 2014 10:25:07 AM UTC+2, Uwe Fechner wrote:
>>>>>>
>>>>>> Hi,
>>>>>> I am working an airborne wind energy as well. I wrote a kite-power 
>>>>>> system simulator in Python, where also one of the controllers (the winch 
>>>>>> controller) is 
>>>>>> implemented in Python. ( https://bitbucket.org/ufechner/freekitesim )
>>>>>>
>>>>>> With Python you can reach a jitter of less than 2 ms in the 20Hz 
>>>>>> control loop quite easily (on low-latency Linux). In my case this is 
>>>>>> sufficient for prototyping,
>>>>>> but the real flight control system should run at a higher update 
>>>>>> frequency (perhaps 200 Hz).
>>>>>>
>>>>>> In contrast to Julia Python is using reference counting, and in my 
>>>>>> Python applications I just turn off the garbage collection.
>>>>>>
>>>>>> For Julia (and I would love to rewrite the simulator in Julia) this 
>>>>>> is probably not an option. A better garbage collector 
>>>>>> (which is in the pipeline, see: 
>>>>>> https://github.com/JuliaLang/juli

Re: [julia-users] Re: error of std function

[Simon Kornblith]

Or alternatively:

std(reshape(A, 10, div(length(A), 10)), 1)

Simon

On Thursday, October 9, 2014 7:10:11 PM UTC-4, Patrick O'Leary wrote:
>
> "optionally *along dimensions in region*" (emphasis mine). You are 
> attempting to read along the tenth dimension of the array.
>
> You're trying to split the array into groups of ten elements, it sounds 
> like.
>
> [std(A[10(n-1)+1:10n]) for n in 1:length(A)./10]
>
> On Thursday, October 9, 2014 5:56:01 PM UTC-5, K leo wrote:
>>
>> I am hoping to get the std's of every 10 consecutive elements in A. 
>>
>> std(v[, region]) 
>> Compute the sample standard deviation of a vector or array v, optionally 
>> along dimensions in region. The algorithm returns an estimator of the 
>> generative distribution’s standard deviation under the assumption that 
>> each entry of v is an IID drawn from that generative distribution. This 
>> computation is equivalent to calculating sqrt(sum((v - mean(v)).^2) / 
>> (length(v) - 1)). Note: Julia does not ignore NaN values in the 
>> computation. For applications requiring the handling of missing data, 
>> the DataArray package is recommended. 
>>
>> On 2014年10月10日 06:49, Patrick O'Leary wrote: 
>> > On Thursday, October 9, 2014 5:42:40 PM UTC-5, K leo wrote: 
>> > 
>> > julia> std(A, 10) 
>> > 
>> > 
>> > A only has elements along the first dimension. What behavior do you 
>> > expect here? 
>>
>>

[julia-users] Re: Is there an idiomatic way to discard values returned from a function

[Simon Kornblith]

It's not really very inefficient. Compare:

julia> f(x, y, z) = return x^2, y^3, z^4;

julia> g(a, b, c) = f(a, b, c)[3];

julia> @code_llvm g(1.0, 1.0, 1.0)

define double @"julia_g;767954"(double, double, double) {
top:
  %3 = call double @pow(double %2, double 4.00e+00), !dbg !525
  ret double %3, !dbg !525
}

julia> function g(a, b, c)
   _, _, v = f(a, b, c)
   v
   end;

julia> @code_llvm g(1.0, 1.0, 1.0)

define double @"julia_g;767955"(double, double, double) {
top:
  %3 = call double @pow(double %2, double 4.00e+00), !dbg !528
  ret double %3, !dbg !530
}

If the return values are bits types, the optimizer is smart enough that 
there's no overhead to pulling unused variables out of tuples. (In this 
case, f was inlined, so g doesn't even compute x^2 or y^3.) At present,  _, 
a, b = f(x, y, z) will actually be more efficient than a, b = f(x, y, 
z)[2:3], since the latter will allocate a new tuple on the heap, although 
that could probably be improved.

If the return values are not bits types, or if the return type of f is not 
inferrable, there will be slight overhead to pulling an unused value out of 
a tuple, but in the former case the time to allocate the unused return 
values will be vastly greater, and in the latter case the type instability 
would be a performance bottleneck no matter what.

Simon

On Wednesday, October 15, 2014 11:23:20 AM UTC-4, Gray Calhoun wrote:
>
> Hi everyone, if a function returns multiple values, is there an idiomatic 
> way to discard all but (say) the 3rd value returned? For example:
>
> f(x, y, z) = return x^2, y^3, z^4
> _, _, v = f(1, 2, 3)
>
> expresses the intent to only use the third value returned, but surely is 
> inefficient.
>
> Thanks!
>

[julia-users] Re: circshift

[Simon Kornblith]

As has been discussed here in the past, "dimension" may be an ambiguous 
term. A matrix has two dimensions, so you're circularly shifting the first 
dimension by 2 (which has no effect since the rows are identical) and the 
second by 1 (which shifts [1 2 3 4 5] to [5 1 2 3 4]). There are no third, 
fourth, or fifth dimensions, so those dimensions are not altered.

If you want to circularly shift each row, you can do that with:

shifts = [2, 1, 1, 1, 1]
d = similar(c)
for i = 1:size(c, 2)
d[i, :] = circshift(c[i, :], (0, shifts[i]))
end

julia> d
5x5 Array{Int64,2}:
 4  5  1  2  3
 5  1  2  3  4
 5  1  2  3  4
 5  1  2  3  4
 5  1  2  3  4

(This could also be done more efficiently with an explicit loop.)

Simon

On Thursday, October 16, 2014 7:33:32 PM UTC-4, Arch Call wrote:
>
> I am trying to use this function little luck.
>
>  manual definition 
> circshift(A, shifts)
> Circularly shift the data in an array. The second argument is a vector 
> giving the amount to shift in each dimension
> ==
>
> Sample script with test data
> ==
> module MyModule
> #-- attempting to use circshift to rearrange columns
> numrows = 5
> a1 = fill(1,numrows,1)
> a2 = fill(2,numrows,1)
> a3 = fill(3,numrows,1)
> a4 = fill(4,numrows,1)
> a5 = fill(5,numrows,1)
> c = hcat(a1,a2,a3,a4,a5)
> println(c)
> d = circshift(c,[2,1,1,1,1])
> println(d)
> end
>
> I just cannot get my head around on how the shifts work for each dimension.
>
> Can anyone elaborate?...Thanks Archie
>
> Sample output
> 
>
> 5 1 2 3 4
>  5 1 2 3 4
>  5 1 2 3 4
>  5 1 2 3 4
>  5 1 2 3 4
>
>
>

[julia-users] Re: Bug in deserialize?

[Simon Kornblith]

Probably https://github.com/JuliaLang/julia/issues/8631 (although I think 
this example alone runs fine; maybe you're loading Color.jl somewhere?)

Simon

On Friday, October 17, 2014 5:21:30 PM UTC-4, Spencer Lyon wrote:
>
> Consider this very simple example
>
> if nprocs() == 1
> addprocs(2)
> end
>
> ll2(F::Matrix, G::Matrix, H::Matrix, obs::Matrix) = sum(F*G*obs*obs'*H)
>
> function create_ll2_args()
> obs = rand(3, 50)
> [(randn(3, 3) * .5, randn(3, 3)*.2, eye(3),  obs) for i=1:2]
> end
>
> function works_with_map()
> args = create_ll2_args()
> map(x->ll2(x...), args)
> end
>
> function craps_out_with_pmap()
> args = create_ll2_args()
> pmap(x->ll2(x...), args)
> end
>
> Now see what happens when I call the last two functions
>
> julia> works_with_map()
> 2-element Array{Float64,1}:
>  6.60313
>  8.82563
>
> julia> workers()  # make sure I have workers
> 2-element Array{Int64,1}:
>  2
>  3
>
> julia> craps_out_with_pmap()
> fatal error on fatal error on 23: : ERROR: `convert` has no method matching 
> convert(::Type{Int64...}, ::Int64)
>  in convert at base.jl:13
>  in convert at ./base.jl:21
>  in deserialize at serialize.jl:447
>  in handle_deserialize at serialize.jl:351
>  in deserialize at serialize.jl:334
>  in anonymous at serialize.jl:354
>  in ntuple at tuple.jl:30
>  in deserialize_tuple at serialize.jl:354
>  in handle_deserialize at serialize.jl:346
>  in deserialize at serialize.jl:334
>  in anonymous at serialize.jl:354
>  in ntuple at tuple.jl:30
>  in deserialize_tuple at serialize.jl:354
>  in handle_deserialize at serialize.jl:346
>  in anonymous at task.jl:853
> ERROR: `convert` has no method matching convert(::Type{Int64...}, ::Int64)
>  in convert at base.jl:13
>  in convert at ./base.jl:21
>  in deserialize at serialize.jl:447
>  in handle_deserialize at serialize.jl:351
>  in deserialize at serialize.jl:334
>  in anonymous at serialize.jl:354
>  in ntuple at tuple.jl:30
>  in deserialize_tuple at serialize.jl:354
>  in handle_deserialize at serialize.jl:346
>  in deserialize at serialize.jl:334
>  in anonymous at serialize.jl:354
>  in ntuple at tuple.jl:30
>  in deserialize_tuple at serialize.jl:354
>  in handle_deserialize at serialize.jl:346
>  in anonymous at task.jl:853
> Worker 3 terminated.
> Worker 2 terminated.
> 2-element Array{Any,1}:
>  ProcessExitedException()
>  ProcessExitedException()
>
> I’m not sure what the problem is. I tried to follow the traceback and 
> ended up at this line 
> , 
> but I’m not really sure what is going wrong or how to fix it. Any ideas?
>
> If this really is a bug we can create an issue for it…
>
> Thanks
> ​
>

[julia-users] Re: General questions/advice - user types, constructors, pmap

[Simon Kornblith]

On Sunday, October 19, 2014 5:25:34 PM UTC-4, Greg Plowman wrote:
>
> Hi,
>
> I have several general questions that came up in my first foray into Julia.
>
> Julia seems such a delight to work with, things seems to work magically 
> and lots of details are not required or implicitly assumed. 
> Whilst this is great for programming, it does mean I'm a little unsure 
> about some things, especially about types, efficiency and what I get for 
> free.   
> In any case, here are some questions:
>
>
> I want to do some parallel simulations and combine/reduce the set of 
> results into a single result. 
>
> As a minimal starting point, I defined a composite type to hold sim 
> results, a single no-argument constructor, and a + function for reducing.
>
>  
>
> type Counters
>
> freqBase::Array{Int64,1}
>
> freqFeature::Array{Int64,1}
>
> freqWin::Array{Int64,1}
>
> freqPrize::Array{Int64,1}
>
> freqCombination::Array{Int64,2}
>
>  
>
> # no-argument constructor
>
> function Counters()
>
> this = new()
>
> this.freqBase = zeros(Int64, 10)
>
> this.freqFeature = zeros(Int64, 10)
>
> this.freqWin = zeros(Int64, 10)
>
> this.freqPrize = zeros(Int64, maxPrize)
>
> this.freqCombination = zeros(Int64, numSymbols, maxState)
>
> return this
>
> end
>
> end
>
>  
>
> function +(c1::Counters, c2::Counters)
>
> c = Counters()
>
> c.freqBase= c1.freqBase + c2.freqBase
>
> c.freqFeature = c1.freqFeature  + c2.freqFeature
>
> c.freqWin = c1.freqWin  + c2.freqWin
>
> c.freqPrize   = c1.freqPrize+ c2.freqPrize
>
> c.freqCombination = c1.freqCombination  + c2.freqCombination
>
> return c
>
> end
>
>  
>
> To my surprise this was sufficient to provide the functionality I needed 
> to return sim result from pmap and reduce to a single result:
>
>  
>
> const numProcessors = 4
>
>  
>
> if nprocs() < numProcessors
>
> addprocs(numProcessors - nprocs())
>
> end
>
>  
>
> const numTrials = 100
>
> const numPlays = 100
>
>  
>
> trialCounts = pmap(Simulation, fill(numPlays, numTrials))
>
> totalCounts = sum(trialCounts)
>
> 
>
> PrintCountersSummary(trialCounts) # print summary for each trial
>
> PrintCounters(totalCounts)# print total combined results
>
>  
>
> Surprisingly (to me) all this worked.
>
>  
>
> Q1. Why does pmap return Vector{Any} rather than Vector{Counters}, when 
> the return type from Simulation() is my user-defined type Counters? 
>
>  
>
> I inserted an extra line:
>
> trialCounts = convert(Vector{Counters}, trialCounts)
>
>  
>
> I was surprised this even worked, because I didn’t define convert.
>

pmap returns a Vector{Any} because the remote data is inherently untyped 
and it doesn't try to tighten the type of the array after retrieving the 
items. convert{T,n,S}(::Type{Array{T,n}}, x::Array{S,n}) is defined here 
.
 
(You can find this using which).
 

> Q2. Is there any advantage to using convert? Is it more efficient? E.g. 
> PrintCounters could be defined to accept argument with Vector{Counters}
>

In some cases it will be more efficient since function lookup can be static 
instead of dynamic, although if the cost of the function call is small 
relative to the cost of running the function it doesn't matter.
 

> Q3. trialCounts is Vector{Any} but sum() works? Presumably sum uses run 
> time type of actual elements of vector?
>

sum is basically just doing something like:

a = x[1]
for i = 2:length(x)
a += x[i]
end

which doesn't actually need type information to work, although type 
information can make it faster.
 

> Q4. Presumably sum() uses my definition of operator +. I also noted that 
> += works. Where are these defined? What else do I get for "free"?
>

x += y is rewritten to x = x + y in the frontend
 

> It occurred to me that my implementation of + could be improved by 
> defining a copy constructor.
>
>  
>
> function Counters(c::Counters)
>
> this = new()
>
> this.freqBase = copy(c.freqBase)
>
> this.freqFeature = copy(c.freqFeature)
>
> this.freqWin = copy(c.freqWin)
>
> this.freqPrize = copy(c.freqPrize)
>
> this.freqCombination = copy(c.freqCombination)
>
> return this
>
> end
>
>  
>
> Then define + as:
>
>  
>
> function +(c1, c2)
>
> c = Counters(c1)
>
>c1.freqBase += c2.freqBase
>
>...
>
>return c
>
> end
>
>  
>
> This would seem to eliminate first initialising with zeros.
>
> However, there was no improvement in practice. Maybe allocation is 
> insignificant compared to addition.
>

copying still requires allocation. This just avoids initializing with 
zeros, which is not usually very expensive.
 

> Then it occurred to me that summing by my definition creates

Re: [julia-users] BLAS a' * b

[Simon Kornblith]

On Tuesday, October 21, 2014 6:41:24 PM UTC-4, David van Leeuwen wrote:
 

> I replaced the `γ = broadcast()` with the lines below that.  No globals, 
> but perhaps the field type gmm.μ is spoiling things.  I am not sure if this 
> is a case of an abstractly typed field
>
> type GMM{T<:FloatingPoint}
> ...
>μ::Array{T}
> ...
> end
>
> Should I have written GMM{T2} in the declaration of stats()?
>

The declaration of stats() is fine, but μ is abstract because it's missing 
the number of dimensions. It should be Array{T,2} or something like that, 
or else GMM should be parametrized on the number of dimensions.

Simon

[julia-users] Displaying a polygon mesh

[Simon Kornblith]

Is there an easy way to display a polygon mesh in Julia, i.e., vertices and 
faces loaded from an STL file or created by marching tetrahedra using 
Meshes.jl? So far, I see:

   - PyPlot/matplotlib, which seems to be surprisingly difficult to 
   convince to do this.
   - GLPlot, which doesn't currently work for me on 0.4. (I haven't tried 
   very hard yet.)
   - ihnorton's VTK bindings, which aren't registered in METADATA.jl. 

Is there another option I'm missing? If not, can I convince one of these 
packages to show my mesh with minimal time investment, or should I use a 
separate volume viewer (or maybe a Python package via PyPlot)?

Thanks,
Simon

Re: [julia-users] Security problem with unitialized memory

[Simon Kornblith]

In general, arrays cannot be assumed to be 16-byte aligned because it's 
always possible to create one that isn't using pointer_to_array. However, 
from Intel's AVX introduction 

:

Intel® AVX has relaxed some memory alignment requirements, so now Intel AVX 
by default allows unaligned access; however, this access may come at a 
performance slowdown, so the old rule of designing your data to be memory 
aligned is still good practice (16-byte aligned for 128-bit access and 
32-byte aligned for 256-bit access).

On Monday, November 24, 2014 10:01:45 PM UTC-5, Erik Schnetter wrote:
>
> On Mon, Nov 24, 2014 at 9:30 PM, Steven G. Johnson 
> > wrote: 
> > Unfortunately, Julia allocates 16-byte aligned data by default (to help 
> SIMD 
> > code), and there is no calloc version of posix_memalign as far as I 
> know. 
>
> The generated machine code I've seen does not make use of this. All 
> the load/store instructions in vectorized or unrolled loops assume 
> unaligned pointers. (Plus, with AVX one should align to 32 bytes 
> instead.) 
>
> -erik 
>
> -- 
> Erik Schnetter > 
> http://www.perimeterinstitute.ca/personal/eschnetter/ 
>

[julia-users] Re: Package name for embedding R within Julia

[Simon Kornblith]

On Friday, January 2, 2015 2:59:04 PM UTC-5, Douglas Bates wrote:
>
> For many statistics-oriented Julia users there is a great advantage in 
> being able to piggy-back on R development and to use at least the data sets 
> from R packages.  Hence the RDatasets package and the read_rda function in 
> the DataFrames package for reading saved R data.
>
> Over the last couple of days I have been experimenting with running an 
> embedded R within Julia and calling R functions from Julia. This is similar 
> in scope to the Rif package except that this code is written in Julia and 
> not as a set of wrapper functions written in C. The R API is a C API and, 
> in some ways, very simple. Everything in R is represented as a "symbolic 
> expression" or SEXPREC and passed around as pointers to such expressions 
> (called an SEXP type).  Most functions take one or more SEXP values as 
> arguments and return an SEXP.
>
> I have avoided reading the code for Rif for two reasons:
>  1. It is GPL3 licensed
>  2. I already know a fair bit of the R API and where to find API function 
> signatures.
>

AFAICT, Rif.jl is GPLv2+. I'm not sure how much a less restrictive license 
helps here. My understanding is that, because R is GPLv2+, code that links 
against it must be redistributed under GPLv2+ or a less restrictive 
license, i.e., while it would be legal to redistribute code that uses 
either Rif.jl or RCall.jl under GPLv2+ or MIT, neither could be used in 
closed source software.

Simon

[julia-users] Re: Calculate 95% confidential intervals of p in binomial distribution

[Simon Kornblith]

You can use:

using HypothesisTests
ci(BinomialTest(x, n))

Several methods of constructing binomial confidence intervals are 
implemented; see the docs 

.

Simon

On Wednesday, December 31, 2014 12:21:24 PM UTC-5, Jerry Xiong wrote:
>
> I want to calculate the 95% confidential intervals of the parameter p of a 
> binomial distribution, for a given x and n. I known that in MATLAT, it 
> could be got in the 2nd output of binofit(x,n,0.95)
> Is there any way to do it in Julia, using Distributions.jl or any python 
> package?
>