[julia-users] Re: Project organization and variable scope question

[Ralph Smith]

Unlike Matlab, Julia doesn't give special treatment to script files.
By "script" we usually mean a file containing a series of expressions
(which can include type and function definitions) intended for
evaluation in the Main module.  The Main module is where expressions
typed interactively at the REPL or IDE, or read from the file
specified on the Julia command line, are evaluated. Source files
implementing parts of other modules differ only in intention.

1) Each module has its own so-called "global" scope, which is where the
contents of "included" files are evaluated. Thus "include" does not
just paste text into your function, as it would in Fortran or C. In
fact, the inclusion is done when the surrounding code is *run*, so the
module being modified may not be the one where the function is
defined; "include" should not normally appear in functions. This
behavior is not yet clearly documented.

2) Global variables have some performance drawbacks, so are not
recommended for your situation.  What I do is to define a composite
type (e.g. "Context") containing most of the run parameters and
widely needed arrays, construct an instance of it in my "main()"
function, and pass the "context" variable as an argument to everything
else.  The "Parameters" package provides some nice macros and methods
for usage like this.

3) Most of us find it wise to organize the pieces of even mid-size
projects into modules. You don't need to make them full packages:
just "include" the top-level module files in your main script. Then
you can access their components with "MyModule.thing", or via "using"
and "import" if you prefer.

Finally a couple of syntax points: you need the word "function" in the
definition of "main()", and Julia doesn't have the "++" operator.



On Friday, November 18, 2016 at 1:25:05 PM UTC-5, Nicholas Mueschke wrote:
>
> Question #1:  If main entry point to run a calculation "main()" is a 
> function, it gets its own variable workspace, right?  Now, if I write a 
> script (not a function) and use include("some_script.jl") with main(), does 
> Julia just inline that code within main()?  In terms of scope, should the 
> script file be able to see all of the variables in the scope of main()?  In 
> Matlab that would be true.  In Fortran/C that wouldn't.  I guess, I'm not 
> sure what scope implications there are for Julia script files.
>
> Question #2:  If I've defined a bunch of functions as shown in the 
> pseudocode above, what is the most performant way to have the large 1D 
> arrays accessible within the scope of each function.  As you can tell, I'm 
> trying to avoid writing functions that accept a long list of input 
> parameters.  The old Fortran solution is to simply make the arrays global, 
> so that each function can access them as needed.  How terrible is that idea 
> within the Julia framework?  Also, how can I even do that?  I've tried 
> writing a script (not a function) to declare a long list of global 
> variables and then used include("DeclareGlobalVariables,jl) within my main. 
>  But, when I return to main(), those variables do not show up in the 
> workspace for main???  What am I missing?
>
> Question #3: I come from a VisualStudio IDE background, so I'm having 
> trouble figuring out how to organize a Juila project.  I'm trying out Atom 
> for my first Julia tests.  For a project that's bigger than just a script 
> or a few functions, should I be defining a defining main entry point 
> function within a module?  Why Does Julia force modules to be added as 
> packages so they can be loaded with the "using" command?  That seems 
> strange.  Or, should I just write everything as a collection of files with 
> functions in them and not worry about modules?  Simple REPL and one file 
> Julia examples are everywhere.  There are also large coding 
> projects/libraries/utilities on github as examples, but I'm having trouble 
> figuring out the structure of these larger projects.  I guess, I'm 
> somewhere in between these two cases, where I'm just want to crunch some 
> numbers, but I'm a little more complicated/sophisticated than the single 
> file examples.  What's the best way to proceed with such a project/file 
> structure?
>
> Thanks in advance for any help.
>
> Nick
>
>
>
>
>
>
>

[julia-users] Re: Why construction of Dict via generator is not type-stable?

[Ralph Smith]

If you can accept an extra call of your function t,


f2{T}(t::Function,a::AbstractArray{T})::Dict{T,typeof(t(a[1])} = Dict((x,t(x
)) for x in a)


credit to mauro3 in discussion of Julia issue 1090.


On Monday, November 14, 2016 at 6:17:38 AM UTC-5, bogumil@gmail.com 
wrote:
>
> Thank you.
> However, my problem is that t(x) does not have to have type T.
> And this is exactly the question - how to determine type of t(x) given 
> that we know that x has type T.
>
> On Monday, November 14, 2016 at 12:29:56 AM UTC+1, Ralph Smith wrote:
>>
>> Until the issue with generators is resolved, this may suffice:
>>
>> f2x{T}(t::Function,a::AbstractArray{T})::Dict{T,T} = Dict((x,t(x)) for x 
>> in a)
>>
>>
>> I think that diverts all the ambiguities to checks by conversion methods.
>>
>> On Sunday, November 13, 2016 at 11:16:20 AM UTC-5, bogumil@gmail.com 
>> wrote:
>>>
>>> I have the folowing questions:
>>>
>>>1. Why when f2 is used as written above the return type is not 
>>>Dict{Float64,Float64} (like in the case of f1 where it is 
>>>Array{Float64,1})?
>>>2. How to fix f2 so that Julia can infer the return type?
>>>3. Is there a way to change f2 so that first empty Dict{>>type of a>,}() variable is created, where >>type of a> and  are determined based on the 
>>>passed arguments, and only next this dictionary is populated with data?
>>>
>>>
>>>

[julia-users] Re: Why construction of Dict via generator is not type-stable?

[Ralph Smith]

Until the issue with generators is resolved, this may suffice:

f2x{T}(t::Function,a::AbstractArray{T})::Dict{T,T} = Dict((x,t(x)) for x in 
a)


I think that diverts all the ambiguities to checks by conversion methods.

On Sunday, November 13, 2016 at 11:16:20 AM UTC-5, bogumil@gmail.com 
wrote:
>
> I have the folowing questions:
>
>1. Why when f2 is used as written above the return type is not 
>Dict{Float64,Float64} (like in the case of f1 where it is 
>Array{Float64,1})?
>2. How to fix f2 so that Julia can infer the return type?
>3. Is there a way to change f2 so that first empty Dict{of a>,}() variable is created, where type of a> and  are determined based on the 
>passed arguments, and only next this dictionary is populated with data?
>
>
>

Re: [julia-users] Getting parameter names from Method

[Ralph Smith]

julia> mt=methods(randjump)
# 2 methods for generic function "randjump":
randjump(r::MersenneTwister, jumps::Integer) at random.jl:145
randjump(mt::MersenneTwister, jumps::Integer, jumppoly::AbstractString) at 
random.jl:137

julia> mt.ms[2].lambda_template.slotnames
8-element Array{Any,1}:
 Symbol("#self#")
 :mt 
 :jumps  
 :jumppoly   
 :mts
 Symbol("#temp#")
 :i  
 :cmt

julia> mt.ms[2].lambda_template.nargs
4




On Friday, November 4, 2016 at 7:00:56 PM UTC-4, David Anthoff wrote:
>
> How do I get at this nargs field? On julia 0.5 Method doesn't seem to have 
> a field with that name? 
>
> > -Original Message- 
> > From: julia...@googlegroups.com  [mailto:julia- 
>  
> > us...@googlegroups.com ] On Behalf Of Yichao Yu 
> > Sent: Friday, November 4, 2016 3:54 PM 
> > To: Julia Users  
> > Subject: Re: [julia-users] Getting parameter names from Method 
> > 
> > On Fri, Nov 4, 2016 at 6:31 PM, David Anthoff  > 
> > wrote: 
> > > Is there a way to get the names of the parameters of a method from a 
> > > Method type instance on julia 0.5? 
> > 
> > Roughly: 
> > 
> > nargs tell you how many arguments the method accepts, the first one 
> being 
> > the object (function) being called. 
> > The local variables names are available in the CodeInfo.slotnames and 
> the 
> > first nargs ones are the parameter names. 
> > The CodeInfo can be found in the source field of the method for non- 
> > generated functions and unspecialized.inferred for generated function. 
> > 
> > > 
> > > 
> > > 
> > > Thanks, 
> > > 
> > > David 
>
>

Re: [julia-users] Question: Forcing readtable to create string type on import

[Ralph Smith]

Unless I misunderstand,

df1 = readtable(file1,eltypes=[String,String,String])


seems to be what you want.

If you're new to Julia, the fact that a "vector of types" really means 
exactly that may be surprising. 

Let us hope that the new versions of DataFrames include a parser that 
doesn't treat most 10-digit numbers as Int32 on systems like yours.

On Wednesday, November 2, 2016 at 4:15:20 PM UTC-4, LeAnthony Mathews wrote:
>
> Spoke too soon.  
> Again I simple want the CSV column that is read in to not be an int32, but 
> a string.
>
> Still having issues casting the CSV file back into a Dataframe.
> Its hard to understand why the Julia system is attempting to determine the 
> type of the columns when I use readtable and I have no control over this.
>
> Why can I not say:
> df1 = readtable(file1; types=Dict(1=>String)) # assuming your account 
> number is column # 1
>
> *Reading the Julia spec-Advanced Options for Reading CSV Files*
> *readtable accepts the following optional keyword arguments:*
>
> *eltypes::Vector{DataType} – Specify the types of all columns. Defaults to 
> [].*
>
>
> *df1 = readtable(file1, Int32::Vector(String))*
>
> I get 
> *ERROR: TypeError: typeassert: expected Array{String,1}, got Type{Int32}*
>
> Is this even an option?  Or how about convert the df1_CSV to 
> df1_dataframe?  
> *df1_dataframe = convert(dataframe, df1_CSV)*
> Since the CSV .read seems to give more granular control.
>
>
> On Tuesday, November 1, 2016 at 7:28:36 PM UTC-4, LeAnthony Mathews wrote:
>>
>> Great, that worked for forcing the column into a string type.
>> Thanks
>>
>> On Monday, October 31, 2016 at 3:26:14 PM UTC-4, Jacob Quinn wrote:
>>>
>>> You could use CSV.jl: http://juliadata.github.io/CSV.jl/stable/
>>>
>>> In this case, you'd do:
>>>
>>> df1 = CSV.read(file1; types=Dict(1=>String)) # assuming your account 
>>> number is column # 1
>>> df2 = CSV.read(file2; types=Dict(1=>String))
>>>
>>> -Jacob
>>>
>>>
>>> On Mon, Oct 31, 2016 at 12:50 PM, LeAnthony Mathews  
>>> wrote:
>>>
 Using v0.5.0
 I have two different 10,000 line CSV files that I am reading into two 
 different dataframe variables using the readtable function.
 Each table has in common a ten digit account_number that I would like 
 to use as an index and join into one master file.

 Here is the account number example in the original CSV from file1:
 8018884596
 8018893530
 8018909633

 When I do a readtable of this CSV into file1 then do a* 
 typeof(file1[:account_number])* I get:
 *DataArrays.DataArray(Int32,1)*
  -571049996
  -571041062
  -571024959

 when I do a 
 *typeof(file2[:account_number])*
 *DataArrays.DataArray(String,1)*


 *Question:  *
 My CSV files give no guidance that account_number should be Int32 or 
 string type.  How do I force it to make both account_number elements type 
 String?

 I would like this join command to work:
 *new_account_join = join(file1, file2, on =:account_number,kind = 
 :left)*

 But I am getting this error:
 *ERROR: TypeError: typeassert: expected Union{Array{Symbol,1},Symbol}, 
 got Array{*
 *Array{Symbol,1},1}*
 * in (::Base.#kw##join)(::Array{Any,1}, ::Base.#join, 
 ::DataFrames.DataFrame, ::D*
 *ataFrames.DataFrame) at .\:0*


 Any help would be appreciated.  



>>>

[julia-users] Re: Error calculating eigs of pentadiagonal matrix.

[Ralph Smith]

Eigs uses shift-and-invert for the :sm variant, so it tries to solve M x = b. 
Try adding a small sigma parameter.

[julia-users] Re: ANN: ParallelAccelerator v0.2 for Julia 0.5 released.

[Ralph Smith]

I looked a bit deeper (i.e. found a machine where I have access to an Intel 
compiler, albeit not up to date - my shop is cursed by budget cuts).  ICC 
breaks up a loop like
for (i=0; i<n; i++) {
  a[i] = exp(cos(b[i]));
  s += a[i];
}
into calls to vector math library functions and a separate loop for the 
sum. The library is bundled with ICC; it's not MKL, but its domain overlaps 
with MKL -
hence my misapprehension - so your point stands. Something like 
blackscholes benefits from these vector library calls, and GCC doesn't do 
that.

It would be nice if Julia's LLVM system included an optimization pass which 
invoked a vector math library when appropriate. I guess that's a challenge 
outside
the scope of ParallelAccelerator, but maybe good ground for some other 
project.

On Thursday, October 27, 2016 at 1:04:33 PM UTC-4, Todd Anderson wrote:
>
> That's interesting.  I generally don't test with gcc  and my experiments 
> with ICC/C have shown something like 20% slower for LLVM/native threads for 
> some class of benchmarks (like blackscholes) but 2-4x slower for some other 
> benchmarks (like laplace-3d).  The 20% may be attributable to ICC being 
> better (including at vectorization like you mention) but certainly not the 
> 2-4x.  These larger differences are still under investigation.
>
> I guess something we have said in the docs or our postings have created 
> this impression that our performance gains are somehow related to MKL or 
> blas in general.  If you have MKL then you can compile Julia to use it 
> through its LLVM path.  ParallelAccelerator does not insert calls to MKL 
> where they didn't exist in the incoming IR and I don't think ICC does 
> either.  If MKL calls exist in the incoming IR then we don't modify them 
> either.  
>
> On Wednesday, October 26, 2016 at 7:51:33 PM UTC-7, Ralph Smith wrote:
>>
>> This is great stuff.  Initial observations (under Linux/GCC) are that 
>> native threads are about 20% faster than OpenMP, so I surmise you are 
>> feeding LLVM some very tasty
>> code. (I tested long loops with straightforward memory access.)
>>
>> On the other hand, some of the earlier posts make me think that you were 
>> leveraging the strong vector optimization of the Intel C compiler and its 
>> tight coupling to
>> MKL libraries.   If so, is there any prospect of getting LLVM to take 
>> advantage of MKL?
>>
>>
>> On Wednesday, October 26, 2016 at 8:13:38 PM UTC-4, Todd Anderson wrote:
>>>
>>> Okay, METADATA with ParallelAccelerator verison 0.2 has been merged so 
>>> if you do a standard Pkg.add() or update() you should get the latest 
>>> version.
>>>
>>> For native threads, please note that we've identified some issues with 
>>> reductions and stencils that have been fixed and we will shortly be 
>>> released in version 0.2.1.  I will post here again when that release takes 
>>> place.
>>>
>>> Again, please give it a try and report back with experiences or file 
>>> bugs.
>>>
>>> thanks!
>>>
>>> Todd
>>>
>>

Re: [julia-users] Recursive data structures with Julia

[Ralph Smith]

Conversion is done by methods listed in base/nullable.jl

I would like to know if there is any drawback to an alternative like

abstract Bst

immutable NullNode <: Bst end

type BstNode <: Bst
val::Int
left::Bst
right::Bst
end

isnull(t::Bst) = isa(t,NullNode)

BstNode(key::Int) = BstNode(key, NullNode(), NullNode())

which appears to be good for type-safety, and is (sometimes) slightly 
faster and less cumbersome than Nullables.

On Sunday, October 30, 2016 at 6:24:42 PM UTC-4, Ángel de Vicente wrote:
>
> Hi, 
>
> by searching in the web I found 
> (http://stackoverflow.com/questions/36383517/how-to-implement-bst-in-julia) 
>
> a way to make my BST code much cleaner (as posted below). Nevertheless, 
> I don't find this very ellegant, since the head node is of type Bst, 
> while the children are of type Nullable{Bst} (is this the 'canonical' way 
> of building recursive data structures with Julia?). 
>
> But when I first read the code in SO, I thought that it was probably 
> wrong, since it does: 
>
> node.left = BST(key) 
> where node.left is of type Nullable{BST}. 
>
> Then I realized that automatic conversion from BST to Nullable{BST} is 
> done when assigning to node.left, so all is good. Coming from Fortran, 
> this is a bit odd for me... what are the rules for automatic conversion?   
>
>
> Thanks a lot, 
> Ángel de Vicente 
>
>
>
>
>
> , 
> | module b 
> | 
> | type Bst 
> | val::Int 
> | left::Nullable{Bst} 
> | right::Nullable{Bst} 
> | end 
> | Bst(key::Int) = Bst(key, Nullable{Bst}(), Nullable{Bst}())   
> | 
> | "Given an array of Ints, it will create a BST tree, type: Bst" 
> | function build_bst(list::Array{Int,1}) 
> | head = list[1] 
> | tree = Bst(head) 
> | for e in list[2:end] 
> | place_bst(tree,e) 
> | end 
> | return tree 
> | end 
> | 
> | function place_bst(tree::Bst,e::Int) 
> | if e == tree.val 
> | println("Dropping $(e). No repeated values allowed") 
> | elseif e < tree.val 
> | if (isnull(tree.left)) 
> | tree.left = Bst(e) 
> | else 
> | place_bst(tree.left.value,e) 
> | end 
> | else 
> | if (isnull(tree.right)) 
> | tree.right = Bst(e) 
> | else 
> | place_bst(tree.right.value,e) 
> | end 
> | end 
> | end 
> | 
> | function print_bst(tree::Bst) 
> | if !isnull(tree.left) print_bst(tree.left.value) end 
> | println(tree.val) 
> | if !isnull(tree.right) print_bst(tree.right.value) end 
> | end 
> | 
> | end 
> ` 
>
> , 
> | julia> include("bst.jl") 
> | 
> | julia> b.print_bst( b.build_bst([4,5,10,3,20,-1,10])) 
> | Dropping 10. No repeated values allowed 
> | -1 
> | 3 
> | 4 
> | 5 
> | 10 
> | 20 
> | 
> | julia> 
> ` 
>
>
> -- 
> Ángel de Vicente 
> http://www.iac.es/galeria/angelv/   
>

[julia-users] Re: ANN: ParallelAccelerator v0.2 for Julia 0.5 released.

[Ralph Smith]

This is great stuff.  Initial observations (under Linux/GCC) are that 
native threads are about 20% faster than OpenMP, so I surmise you are 
feeding LLVM some very tasty
code. (I tested long loops with straightforward memory access.)

On the other hand, some of the earlier posts make me think that you were 
leveraging the strong vector optimization of the Intel C compiler and its 
tight coupling to
MKL libraries.   If so, is there any prospect of getting LLVM to take 
advantage of MKL?


On Wednesday, October 26, 2016 at 8:13:38 PM UTC-4, Todd Anderson wrote:
>
> Okay, METADATA with ParallelAccelerator verison 0.2 has been merged so if 
> you do a standard Pkg.add() or update() you should get the latest version.
>
> For native threads, please note that we've identified some issues with 
> reductions and stencils that have been fixed and we will shortly be 
> released in version 0.2.1.  I will post here again when that release takes 
> place.
>
> Again, please give it a try and report back with experiences or file bugs.
>
> thanks!
>
> Todd
>

Re: [julia-users] Re: BLAS.set_num_threads() and peakflops scaling

[Ralph Smith]

On looking more carefully, I believe I was mistaken about thread assignment 
to cores - that seems to be done well in OpenBLAS (and maybe Linux in 
general nowadays).  Perhaps the erratic benchmarks under hyperthreading - 
even after heap management is tamed - arise when the operating system 
detects idle virtual cores and schedules disruptive processes there.

On Friday, October 21, 2016 at 12:09:07 AM UTC-4, Ralph Smith wrote:
>
> That's interesting, I see the code in OpenBLAS. However, on the Linux 
> systems I use, when I had hyperthreading enabled the allocations looked 
> random, and I generally got less consistent benchmarks.  I'll have to check 
> that again.
>
> You can also avoid the memory allocation effects by something like
> using BenchmarkTools
> a = rand(n,n); b=rand(n,n); c = similar(a);
> @benchmark A_mul_B!($c,$a,$b)
>
> Of course this is only directly relevant to your real workload if that is 
> dominated by sections where you can optimize away allocations and memory 
> latency.
>
>
> On Thursday, October 20, 2016 at 11:00:41 PM UTC-4, Thomas Covert wrote:
>>
>> Thanks - I will try to figure out how to do that.  I will note, however, 
>> that the OpenBLAS FAQ suggests that OpenBLAS tries to avoid allocating 
>> threads to the same physical core on machines with hyper threading, so 
>> perhaps this is not the cause:
>>
>> https://github.com/xianyi/OpenBLAS/blob/master/GotoBLAS_03FAQ.txt
>>
>>
>>
>> On Thursday, October 20, 2016 at 4:45:51 PM UTC-5, Stefan Karpinski wrote:
>>>
>>> I think Ralph is suggesting that you disable the CPU's hyperthreading if 
>>> you run this kind of code often. We've done that on our benchmarking 
>>> machines, for example.
>>>
>>> On Wed, Oct 19, 2016 at 11:47 PM, Thomas Covert <thom@gmail.com> 
>>> wrote:
>>>
>>>> So are you suggesting that real numerical workloads under 
>>>> BLAS.set_num_threads(4) will indeed be faster than 
>>>> under BLAS.set_num_threads(2)?  That hasn't been my experience and I 
>>>> figured the peakflops() example would constitute an MWE.  Is there another 
>>>> workload you would suggest I try to figure out if this is just a peak 
>>>> flops() aberration or a real issue?
>>>>
>>>>
>>>> On Wednesday, October 19, 2016 at 8:28:16 PM UTC-5, Ralph Smith wrote:
>>>>>
>>>>> At least 2 things contribute to erratic results from peakflops(). With 
>>>>> hyperthreading, the threads are not always put on separate cores. 
>>>>> Secondly, 
>>>>> the measured time includes
>>>>> the allocation of the result matrix, so garbage collection affects 
>>>>> some of the results.  Most available advice says to disable 
>>>>> hyperthreading 
>>>>> on dedicated number crunchers
>>>>> (most full loads work slightly more efficiently without the extra 
>>>>> context switching).  The GC issue seems to be a mistake, if "peak" is to 
>>>>> be 
>>>>> taken seriously.
>>>>>
>>>>> On Wednesday, October 19, 2016 at 12:04:00 PM UTC-4, Thomas Covert 
>>>>> wrote:
>>>>>>
>>>>>> I have a recent iMac with 4 logical cores (and 8 hyper threads).  I 
>>>>>> would have thought that peakflops(N) for a large enough N should be 
>>>>>> increasing in the number of threads I allow BLAS to use.  I do find that 
>>>>>> peakflops(N) with 1 thread is about half as high as peakflops(N) with 2 
>>>>>> threads, but there is no gain to 4 threads.  Are my expectations wrong 
>>>>>> here, or is it possible that BLAS is somehow configured incorrectly on 
>>>>>> my 
>>>>>> machine?  In the example below, N = 6755, a number relevant for my work, 
>>>>>> but the results are similar with 5000 or 1.
>>>>>>
>>>>>> here is my versioninfo()
>>>>>> julia> versioninfo()
>>>>>> Julia Version 0.5.0
>>>>>> Commit 3c9d753* (2016-09-19 18:14 UTC)
>>>>>> Platform Info:
>>>>>>   System: Darwin (x86_64-apple-darwin15.6.0)
>>>>>>   CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
>>>>>>   WORD_SIZE: 64
>>>>>>   BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Haswell)
>>>>>>   LAPACK: libopenblas
>>>>>>   LIBM: libopenlibm
>>>>>>   LLVM: libLLVM-3.7.1 (ORCJIT, haswell)
>>>>>>
>>>>>> here is an example peakflops() exercise:
>>>>>> julia> BLAS.set_num_threads(1)
>>>>>>
>>>>>> julia> mean(peakflops(6755) for i=1:10)
>>>>>> 5.225580459387056e10
>>>>>>
>>>>>> julia> BLAS.set_num_threads(2)
>>>>>>
>>>>>> julia> mean(peakflops(6755) for i=1:10)
>>>>>> 1.004317640281997e11
>>>>>>
>>>>>> julia> BLAS.set_num_threads(4)
>>>>>>
>>>>>> julia> mean(peakflops(6755) for i=1:10)
>>>>>> 9.838116463900085e10
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>

Re: [julia-users] Re: BLAS.set_num_threads() and peakflops scaling

[Ralph Smith]

That's interesting, I see the code in OpenBLAS. However, on the Linux 
systems I use, when I had hyperthreading enabled the allocations looked 
random, and I generally got less consistent benchmarks.  I'll have to check 
that again.

You can also avoid the memory allocation effects by something like
using BenchmarkTools
a = rand(n,n); b=rand(n,n); c = similar(a);
@benchmark A_mul_B!($c,$a,$b)

Of course this is only directly relevant to your real workload if that is 
dominated by sections where you can optimize away allocations and memory 
latency.


On Thursday, October 20, 2016 at 11:00:41 PM UTC-4, Thomas Covert wrote:
>
> Thanks - I will try to figure out how to do that.  I will note, however, 
> that the OpenBLAS FAQ suggests that OpenBLAS tries to avoid allocating 
> threads to the same physical core on machines with hyper threading, so 
> perhaps this is not the cause:
>
> https://github.com/xianyi/OpenBLAS/blob/master/GotoBLAS_03FAQ.txt
>
>
>
> On Thursday, October 20, 2016 at 4:45:51 PM UTC-5, Stefan Karpinski wrote:
>>
>> I think Ralph is suggesting that you disable the CPU's hyperthreading if 
>> you run this kind of code often. We've done that on our benchmarking 
>> machines, for example.
>>
>> On Wed, Oct 19, 2016 at 11:47 PM, Thomas Covert <thom@gmail.com> 
>> wrote:
>>
>>> So are you suggesting that real numerical workloads under 
>>> BLAS.set_num_threads(4) will indeed be faster than 
>>> under BLAS.set_num_threads(2)?  That hasn't been my experience and I 
>>> figured the peakflops() example would constitute an MWE.  Is there another 
>>> workload you would suggest I try to figure out if this is just a peak 
>>> flops() aberration or a real issue?
>>>
>>>
>>> On Wednesday, October 19, 2016 at 8:28:16 PM UTC-5, Ralph Smith wrote:
>>>>
>>>> At least 2 things contribute to erratic results from peakflops(). With 
>>>> hyperthreading, the threads are not always put on separate cores. 
>>>> Secondly, 
>>>> the measured time includes
>>>> the allocation of the result matrix, so garbage collection affects some 
>>>> of the results.  Most available advice says to disable hyperthreading on 
>>>> dedicated number crunchers
>>>> (most full loads work slightly more efficiently without the extra 
>>>> context switching).  The GC issue seems to be a mistake, if "peak" is to 
>>>> be 
>>>> taken seriously.
>>>>
>>>> On Wednesday, October 19, 2016 at 12:04:00 PM UTC-4, Thomas Covert 
>>>> wrote:
>>>>>
>>>>> I have a recent iMac with 4 logical cores (and 8 hyper threads).  I 
>>>>> would have thought that peakflops(N) for a large enough N should be 
>>>>> increasing in the number of threads I allow BLAS to use.  I do find that 
>>>>> peakflops(N) with 1 thread is about half as high as peakflops(N) with 2 
>>>>> threads, but there is no gain to 4 threads.  Are my expectations wrong 
>>>>> here, or is it possible that BLAS is somehow configured incorrectly on my 
>>>>> machine?  In the example below, N = 6755, a number relevant for my work, 
>>>>> but the results are similar with 5000 or 1.
>>>>>
>>>>> here is my versioninfo()
>>>>> julia> versioninfo()
>>>>> Julia Version 0.5.0
>>>>> Commit 3c9d753* (2016-09-19 18:14 UTC)
>>>>> Platform Info:
>>>>>   System: Darwin (x86_64-apple-darwin15.6.0)
>>>>>   CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
>>>>>   WORD_SIZE: 64
>>>>>   BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Haswell)
>>>>>   LAPACK: libopenblas
>>>>>   LIBM: libopenlibm
>>>>>   LLVM: libLLVM-3.7.1 (ORCJIT, haswell)
>>>>>
>>>>> here is an example peakflops() exercise:
>>>>> julia> BLAS.set_num_threads(1)
>>>>>
>>>>> julia> mean(peakflops(6755) for i=1:10)
>>>>> 5.225580459387056e10
>>>>>
>>>>> julia> BLAS.set_num_threads(2)
>>>>>
>>>>> julia> mean(peakflops(6755) for i=1:10)
>>>>> 1.004317640281997e11
>>>>>
>>>>> julia> BLAS.set_num_threads(4)
>>>>>
>>>>> julia> mean(peakflops(6755) for i=1:10)
>>>>> 9.838116463900085e10
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>

[julia-users] Re: BLAS.set_num_threads() and peakflops scaling

[Ralph Smith]

At least 2 things contribute to erratic results from peakflops(). With 
hyperthreading, the threads are not always put on separate cores. Secondly, 
the measured time includes
the allocation of the result matrix, so garbage collection affects some of 
the results.  Most available advice says to disable hyperthreading on 
dedicated number crunchers
(most full loads work slightly more efficiently without the extra context 
switching).  The GC issue seems to be a mistake, if "peak" is to be taken 
seriously.

On Wednesday, October 19, 2016 at 12:04:00 PM UTC-4, Thomas Covert wrote:
>
> I have a recent iMac with 4 logical cores (and 8 hyper threads).  I would 
> have thought that peakflops(N) for a large enough N should be increasing in 
> the number of threads I allow BLAS to use.  I do find that peakflops(N) 
> with 1 thread is about half as high as peakflops(N) with 2 threads, but 
> there is no gain to 4 threads.  Are my expectations wrong here, or is it 
> possible that BLAS is somehow configured incorrectly on my machine?  In the 
> example below, N = 6755, a number relevant for my work, but the results are 
> similar with 5000 or 1.
>
> here is my versioninfo()
> julia> versioninfo()
> Julia Version 0.5.0
> Commit 3c9d753* (2016-09-19 18:14 UTC)
> Platform Info:
>   System: Darwin (x86_64-apple-darwin15.6.0)
>   CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
>   WORD_SIZE: 64
>   BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY Haswell)
>   LAPACK: libopenblas
>   LIBM: libopenlibm
>   LLVM: libLLVM-3.7.1 (ORCJIT, haswell)
>
> here is an example peakflops() exercise:
> julia> BLAS.set_num_threads(1)
>
> julia> mean(peakflops(6755) for i=1:10)
> 5.225580459387056e10
>
> julia> BLAS.set_num_threads(2)
>
> julia> mean(peakflops(6755) for i=1:10)
> 1.004317640281997e11
>
> julia> BLAS.set_num_threads(4)
>
> julia> mean(peakflops(6755) for i=1:10)
> 9.838116463900085e10
>
>
>
>
>
>

[julia-users] Re: mktemp doesn't accept anonymous function?

[Ralph Smith]

specifically: 

>   mktemp(parent) Returns (path, io), where path is the path of a new 
>> temporary file in parent
>
>   and io is an open file object for this path.
>
>
>> so you would normally use the second entry (x in my example) in the 
function block. 

On Wednesday, October 19, 2016 at 5:45:56 PM UTC-4, Ralph Smith wrote:
>
> mktemp returns a 2-tuple, so (at least in v0.5)
>
> julia> mktemp() do y,x
>println(y)
>end
> /tmp/tmpS63NPs
>
> julia> mktemp((y,x) -> println(y))
> /tmp/tmpA2p0Y1
>
> julia> 
>
>
>
>
> On Wednesday, October 19, 2016 at 5:06:55 PM UTC-4, John Brock wrote:
>>
>>
>>
>> Am I doing something wrong or is this a bug? It seems like if the 1st 
>> version works, the 2nd and 3rd should, too.
>>
>> julia> mktemp(println)
>> /var/folders/jd/1skd5rh11hnc_s19lmx93zywgp/T/tmpf7HaUHIOStream(> >)
>>
>>
>> julia> mktemp(x->println(x))
>> ERROR: wrong number of arguments
>>  in anonymous at none:1
>>  in mktemp at file.jl:218
>>  in mktemp at file.jl:216
>>
>>
>> julia> mktemp() do x
>>  println(x)
>>end
>> ERROR: wrong number of arguments
>>  in anonymous at none:2
>>  in mktemp at file.jl:218
>>  in mktemp at file.jl:216
>>  
>>
>>
>> julia> versioninfo()
>> Julia Version 0.4.7
>> Commit ae26b25* (2016-09-18 16:17 UTC)
>> Platform Info:
>>   System: Darwin (x86_64-apple-darwin15.6.0)
>>   CPU: Intel(R) Core(TM) i7-4980HQ CPU @ 2.80GHz
>>   WORD_SIZE: 64
>>   BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
>>   LAPACK: libopenblas64_
>>   LIBM: libopenlibm
>>   LLVM: libLLVM-3.3
>>
>>

[julia-users] Re: mktemp doesn't accept anonymous function?

[Ralph Smith]

mktemp returns a 2-tuple, so (at least in v0.5)

julia> mktemp() do y,x
   println(y)
   end
/tmp/tmpS63NPs

julia> mktemp((y,x) -> println(y))
/tmp/tmpA2p0Y1

julia> 




On Wednesday, October 19, 2016 at 5:06:55 PM UTC-4, John Brock wrote:
>
>
>
> Am I doing something wrong or is this a bug? It seems like if the 1st 
> version works, the 2nd and 3rd should, too.
>
> julia> mktemp(println)
> /var/folders/jd/1skd5rh11hnc_s19lmx93zywgp/T/tmpf7HaUHIOStream( >)
>
>
> julia> mktemp(x->println(x))
> ERROR: wrong number of arguments
>  in anonymous at none:1
>  in mktemp at file.jl:218
>  in mktemp at file.jl:216
>
>
> julia> mktemp() do x
>  println(x)
>end
> ERROR: wrong number of arguments
>  in anonymous at none:2
>  in mktemp at file.jl:218
>  in mktemp at file.jl:216
>  
>
>
> julia> versioninfo()
> Julia Version 0.4.7
> Commit ae26b25* (2016-09-18 16:17 UTC)
> Platform Info:
>   System: Darwin (x86_64-apple-darwin15.6.0)
>   CPU: Intel(R) Core(TM) i7-4980HQ CPU @ 2.80GHz
>   WORD_SIZE: 64
>   BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
>   LAPACK: libopenblas64_
>   LIBM: libopenlibm
>   LLVM: libLLVM-3.3
>
>

[julia-users] Re: Parallel file access

[Ralph Smith]

How are the processes supposed to interact with the database?  Without 
extra synchronization logic, SQLite.jl gives (occasionally)
ERROR: LoadError: On worker 2:
SQLite.SQLiteException("database is locked")
which on the face of it suggests that all workers are using the same 
connection, although I opened the DB separately in each process.
(I think we should get "busy" instead of "locked", but then still have no 
good way to test for this and wait for a wake-up signal.)
So we seem to be at least as badly off as the original post, except with DB 
calls instead of simple writes.

We shouldn't have to stand up a separate multithreaded DB server just for 
this. Would you be kind enough to give us an example of simple (i.e. not 
client-server) multiprocess DB access in Julia?

On Saturday, October 15, 2016 at 9:40:17 AM UTC-4, Steven Sagaert wrote:
>
> It still surprises me how in the scientific computing field people still 
> refuse to learn about databases and then replicate database functionality 
> in files in a complicated and probably buggy way. HDF5  is one example, 
> there are many others. If you want to to fancy search (i.e. speedup search 
> via indices) or do things like parallel writes/concurrency you REALLY 
> should use databases. That's what they were invented for decades ago. 
> Nowadays there a bigger choice than ever: Relational or non-relational 
> (NOSQL), single host or distributed, web interface or not,  disk-based or 
> in-memory,... There really is no excuse anymore not to use a database if 
> you want to go beyond just reading in a bunch of data in one go in memory.
>
> On Monday, October 10, 2016 at 5:09:39 PM UTC+2, Zachary Roth wrote:
>>
>> Hi, everyone,
>>
>> I'm trying to save to a single file from multiple worker processes, but 
>> don't know of a nice way to coordinate this.  When I don't coordinate, 
>> saving works fine much of the time.  But I sometimes get errors with 
>> reading/writing of files, which I'm assuming is happening because multiple 
>> processes are trying to use the same file simultaneously.
>>
>> I tried to coordinate this with a queue/channel of `Condition`s managed 
>> by a task running in process 1, but this isn't working for me.  I've tried 
>> to simiplify this to track down the problem.  At least part of the issue 
>> seems to be writing to the channel from process 2.  Specifically, when I 
>> `put!` something onto a channel (or `push!` onto an array) from process 2, 
>> the channel/array is still empty back on process 1.  I feel like I'm 
>> missing something simple.  Is there an easier way to go about coordinating 
>> multiple processes that are trying to access the same file?  If not, does 
>> anyone have any tips?
>>
>> Thanks for any help you can offer.
>>
>> Cheers,
>> ---Zachary
>>
>

[julia-users] Re: Parallel file access

[Ralph Smith]

There are good synchronization primitives for Tasks, and a bit for threads, 
but not much for parallel processes. (One could use named system semaphores 
on Linux and Windows, but there's no Julia wrapper yet AFAIK.)
I also found management of parallel processes confusing, and good 
nontrivial examples are not obvious.  So here's my humble offering:
https://gist.github.com/RalphAS/2a37b1c631923efa30ac4a7f02a2ee9d

(I just happen to be working on some synchronization problems this month; 
if a real expert has a better solution, let's hear it.)

On Friday, October 14, 2016 at 3:19:25 PM UTC-4, Zachary Roth wrote:
>
> Thanks for the reply and suggestion, Ralph.  I tried to get this working 
> with semaphores/mutexes/locks/etc.  But I've not been having any luck.
>
> Here's a simplified, incomplete version of what I'm trying to do.  I'm 
> hoping that someone can offer a suggestion if they see some sample code.  
>
> function localfunction()
> files = listfiles()
> locks = [Threads.SpinLock() for _ in files]
> ranges = getindexranges(length(files))
>
> pmap(pairs(ranges)) do rows_and_cols
> rows, cols = rows_and_cols
> workerfunction(files, locks, rows, cols)
> end
> end
>
> function workerfunction(files, locks, rows, cols)
> data = kindofexpensive(...)
> pairs = pairs(rows, cols)
> 
> @sync for idx in unique([rows; cols])
> @async begin
> lock(locks[idx])
> updatefile(files[idx], data[idx])
> unlock(locks[idx])
> end
> end
> end
>
> This (obviously) does not work.  I think that the problem is that the 
> locks are being copied when the function is spawned on each process.  I've 
> tried wrapping the locks/semaphores in Futures/RemoteChannels, but that 
> also hasn't worked for me.
>
> I found that I could do the sort of coordination that I need by starting 
> Tasks on the local process.  More specifically, each file would have an 
> associated Task to handle the coordination between processes.  But this 
> only worked for me in a simplified situation with the Tasks being declared 
> globally.  When I tried to implement this coordination within localfunction 
> above, I got an error (really a bunch of errors) that said that a running 
> Task cannot be serialized.
>
> Sorry for the long post, but I'm really hoping that someone can help me 
> out.  I have a feeling that I'm missing something pretty simple.
>
> ---Zachary
>
>
>
>
> On Tuesday, October 11, 2016 at 10:15:06 AM UTC-4, Ralph Smith wrote:
>>
>> You can do it with 2 (e.g. integer) channels per worker (requests and 
>> replies) and a task for each pair in the main process. That's so ugly I'd 
>> be tempted to write an
>>  interface to named system semaphores. Or just use a separate file for 
>> each worker.
>>
>> On Monday, October 10, 2016 at 11:09:39 AM UTC-4, Zachary Roth wrote:
>>>
>>> Hi, everyone,
>>>
>>> I'm trying to save to a single file from multiple worker processes, but 
>>> don't know of a nice way to coordinate this.  When I don't coordinate, 
>>> saving works fine much of the time.  But I sometimes get errors with 
>>> reading/writing of files, which I'm assuming is happening because multiple 
>>> processes are trying to use the same file simultaneously.
>>>
>>> I tried to coordinate this with a queue/channel of `Condition`s managed 
>>> by a task running in process 1, but this isn't working for me.  I've tried 
>>> to simiplify this to track down the problem.  At least part of the issue 
>>> seems to be writing to the channel from process 2.  Specifically, when I 
>>> `put!` something onto a channel (or `push!` onto an array) from process 2, 
>>> the channel/array is still empty back on process 1.  I feel like I'm 
>>> missing something simple.  Is there an easier way to go about coordinating 
>>> multiple processes that are trying to access the same file?  If not, does 
>>> anyone have any tips?
>>>
>>> Thanks for any help you can offer.
>>>
>>> Cheers,
>>> ---Zachary
>>>
>>

Re: [julia-users] python-like generators?

[Ralph Smith]

According to the manual, Tasks are coroutines. The Wikipedia article 
 relates them to generators.  For 
any recovering C programmers in the audience, this related article 
 might remind 
you why you prefer to be here in the land of honest macros.

On Thursday, October 13, 2016 at 6:02:18 PM UTC-4, David Anthoff wrote:
>
> I don't think so. I think the python generator functions are modeled after 
> the C# yield keyword, which is a way to implement an iterator. I believe 
> tasks are quite a different beast. At least in C# tasks interact with the 
> async/await story really well. I think the julia tasks might be similar to 
> that, but I'm not sure. 
>
> In any case, having a yield keyword to implement standard iterators would 
> be 
> fantastic. Simply a way to get the start, next and done methods 
> implemented 
> automatically. 
>
> > -Original Message- 
> > From: julia...@googlegroups.com  [mailto:julia- 
>  
> > us...@googlegroups.com ] On Behalf Of Mauro 
> > Sent: Thursday, October 13, 2016 6:00 AM 
> > To: julia...@googlegroups.com  
> > Subject: Re: [julia-users] python-like generators? 
> > 
> > Isn't this the same as tasks in Julia? 
> > http://docs.julialang.org/en/release-0.5/stdlib/parallel/ 
> > 
> > Although, note that their performance is not on par with start-next-done 
> > iteration. 
> > 
> > On Thu, 2016-10-13 at 14:46, Neal Becker  > wrote: 
> > > julia-0.5 supports generator expressions, ala python, which is very 
> nice. 
> > > 
> > > Any thoughts on supporting the more general python generator 
> > > functions, as described e.g.: https://wiki.python.org/moin/Generators? 
>
> > > I haven't used them much myself (well, once), but they seem a really 
> cool 
> > idea. 
>
>

[julia-users] Re: Parallel file access

[Ralph Smith]

You can do it with 2 (e.g. integer) channels per worker (requests and 
replies) and a task for each pair in the main process. That's so ugly I'd 
be tempted to write an
 interface to named system semaphores. Or just use a separate file for each 
worker.

On Monday, October 10, 2016 at 11:09:39 AM UTC-4, Zachary Roth wrote:
>
> Hi, everyone,
>
> I'm trying to save to a single file from multiple worker processes, but 
> don't know of a nice way to coordinate this.  When I don't coordinate, 
> saving works fine much of the time.  But I sometimes get errors with 
> reading/writing of files, which I'm assuming is happening because multiple 
> processes are trying to use the same file simultaneously.
>
> I tried to coordinate this with a queue/channel of `Condition`s managed by 
> a task running in process 1, but this isn't working for me.  I've tried to 
> simiplify this to track down the problem.  At least part of the issue seems 
> to be writing to the channel from process 2.  Specifically, when I `put!` 
> something onto a channel (or `push!` onto an array) from process 2, the 
> channel/array is still empty back on process 1.  I feel like I'm missing 
> something simple.  Is there an easier way to go about coordinating multiple 
> processes that are trying to access the same file?  If not, does anyone 
> have any tips?
>
> Thanks for any help you can offer.
>
> Cheers,
> ---Zachary
>

[julia-users] Re: Warning since 0.5

[Ralph Smith]

This is indeed fairly confusing, partly because the terminological 
conventions of functional programming are somewhat arbitrary. Some of my 
discussion upthread
uses these conventions, but some (e.g. my use of "attach" and 
"associate/dissociate") may not be conventional.

In Julia, and as Steven J & I used it above, "name" is a field of a 
function object. Thus
f(x)=x^2
sets the name and other fields all at once, so it's never anonymous. In 
particular, it attaches a method for argument type "Any" (since none was 
specified).

Methods are effectively rules for applying the function to tuples of 
argument types.
Methods in turn have "specializations" which seem to be (essentially) the 
compiled code for the specific leaf types of arguments, so
a = f(2.0)
gets Julia to compile a specialization of the (f,Any) method for 
(f,Float64), and save it for later use by stashing it in a subfield of the 
function object.

f(x::Int)=x^2+1
would add another method to the function named "f", just to be used with 
integer arguments.  But
f(x)=x^3
replaces the first method (because it has the same nominal argument type, 
i.e. Any) That is, it replaces it within the object we set up with the 
original definition.
Operations like map, broadcast, and closure construction use and store 
references to parts of the function object (in particular, methods), and 
associate the name "f"
with those references. In Julia v0.5, these references are not updated 
correctly on redefinition. That's what I meant when I said the 
disassociation is incomplete.


On Monday, October 10, 2016 at 3:09:24 PM UTC-4, digxx wrote:
>
> The function object then points into
>> method tables. You can't assign the name "f" to a different function 
>> object, just attach different methods to it.  The troubles seem to arise 
>> from cached references to
>> orphaned method table entries, which are not completely dissociated from 
>> the object named "f". 
>> f(x)=x^2
>>
>
>  so this definition associates the anonymous function with a name (and 
> then it is not anoymous anymore?)
> However how do these method tables look like? In particular the name f is 
> then associated to these methods and for some reason one does not want to 
> be able to change f to refer to different methods?
> How would different mathod attaching look like?
> when I write
>
> f(x)=x^2
> and then
> f(x)=x^3
>
> the old method (is the method here the operation x^2 or x^3?) is still 
> there and x^3 is just attached to it?!?
> is it possible to completely dissociate the old method (in this case as I 
> understand it correctly it is x^2?) with the name f?
>

[julia-users] Re: Warning since 0.5

[Ralph Smith]

By "module file" I just meant a source code file where all definitions are 
enclosed in modules, so if you "include" it, it replaces the whole module. 
Thus

module M
function f(x)
   x^2
end
end


then references to M.f are more likely to be consistent than a bare 
function f defined at the top level (the Main module).

As you surmised, 

g = x -> x^2


binds the mutable variable with symbol :g  to an anonymous function (i.e. 
one with a hidden name and type).  If you bind :g to a different function, 
references to :g should
always use the new association.  Defining a regular function (like M.f 
above) embeds "f" as a name in the function object itself. The function 
object then points into
method tables. You can't assign the name "f" to a different function 
object, just attach different methods to it.  The troubles seem to arise 
from cached references to
orphaned method table entries, which are not completely dissociated from 
the object named "f". 

On Saturday, October 8, 2016 at 5:09:19 PM UTC-4, digxx wrote:
 | what do u mean by "module files"?

 | So what is the difference between
 | f(x)=x^2
 | and 
 | f=x->x^2

is f(x)=x^2 not an anonymous function?!?!
>

[julia-users] Re: Warning since 0.5

[Ralph Smith]

The problem is that some Julia processing stores references to definitions 
in hidden locations which are not updated consistently, so you get 
inconsistency like this:

julia> f(x)=x^2
f (generic function with 1 method)


julia> map(f,[1,2,3])
3-element Array{Int64,1}:
 1
 4
 9


julia> f(x)=x^2+1
WARNING: Method definition f(Any) in module Main at REPL[7]:1 overwritten 
at REPL[9]:1.
f (generic function with 1 method)


julia> [f(x) for x in [1,2,3]]
3-element Array{Int64,1}:
 1
 4
 9


julia> for x in [1,2,3]; println(f(x)); end
2
5
10


julia> map(f,[1,2,3])
3-element Array{Int64,1}:
 1
 4
 9


(Thanks to fcard for pointing this out.  See this issue 
 for more discussion.) 
 Fixing this and related problems is hard, but there is hope for v0.6.

Anyway, my point was that "function files" are more dangerous than "module 
files", and that's why they generate more warnings.

In the REPL one can also drop *all* old definitions & bindings by using 
workspace(), but I find that more painful than just using modules and 
restarting Julia when
necessary.

If you're convinced that all this doesn't apply to your case, you can 
suppress messages: 
https://github.com/cstjean/ClobberingReload.jl#silencing-warnings 

On Saturday, October 8, 2016 at 9:35:10 AM UTC-4, digxx wrote:
>
> Hey,
> Thx for ur answer. So The first time I call my program which includes a 
> file with function definitions there is no problem. I do this because with 
> 0.4 parallel loops didnt work with functions which are defined in the same 
> file even though an @everywhere is prefixed.
> I still dont understand why this should happen at all. Isnt it totally 
> natural that a function file where things like
> g(x)=x.^2
> are defined and not much more complex it might happen that I want to 
> change it to 
> g(x)=x.^2+1
> without restarting Julia in order to not get the error.
> This also happens in the REPL when I overwrite a simple function. So it is 
> probably not correlated to the @everywhere.
> How can I avoid the warning or is it possible to "free" all allocated 
> function definitions before I redefine them?
>

[julia-users] Re: Warning since 0.5

[Ralph Smith]

TL;DR: put everything in modules.

"WARNING: Method definition ... overwritten..."

This warning was extended to more cases in v0.5.  It has caused some 
confusion, and some available explanations are themselves confusing.  Here 
is another try.

Overwriting a method (or type) has always been a formally undefined action 
- that is, there is no guarantee of consistency when you refer to the 
method after
overwriting it (  issue 265 
).  In v0.5 there were 
changes to the implementation of functions which make it more likely that 
this will cause conflicts and problems, 
so it is good that the warning is more prevalent.

If you put your method definitions in modules, and reload the module files, 
you just get warnings for the module redefinitions. You would then also 
have to reload
any modules which use the redefined ones, and so on. It is possible to 
cause conflicts with consistently redefined modules, but it takes some 
effort - unless you
are foolish enough to type
using MyVolatileModule
in interactive sessions; that is asking for trouble.  If even your global 
variables are in modules, consistent reloading will redefine them in terms 
of
updated used modules. For production work (i.e. if you care about valid 
results) you should not overwrite definitions.  If you use packages which 
themselves have
overwrites, please file an issue.

In your case, if the `@everywhere` clause includes definitions which were 
really not changed (not just the file text, but also any global context it 
refers to), you might not
encounter conflicts. It is nevertheless a deprecated usage.  It would seem 
to be either pointless or dangerous.

Caveat: I am not a Julia developer, just a user with a tendency to stumble 
into edge cases.  I happily defer to corrections from the experts, but 
regret that they didn't
explain this issue in the regular documentation.

On Thursday, October 6, 2016 at 7:42:36 PM UTC-4, digxx wrote:
>
> Something has changed since the code which ran with 0.4 does return a 
> warning now:
>
> include("C:\\Users\\Diger\\Documents\\Julia\\Diffusion\\bessel0_expansion_nlsolve.jl")
> WARNING: Method definition WARNING: Method definition WARNING: Method 
> definition WARNING: Method definition WARNING: Method definition g_r(Any) 
> in module Main at 
> C:\Users\Diger\Documents\Julia\Rohrstroemung_function.jl:11 overwritten at 
> C:\Users\Diger\Documents\Julia\Rohrstroemung_function.jl:11.
> WARNING: Method definition g_i(Any, Any) in module Main at 
> C:\Users\Diger\Documents\Julia\Rohrstroemung_function.jl:36 overwritten at 
> C:\Users\Diger\Documents\Julia\Rohrstroemung_function.jl:36.
>
>
> I had a function file which i included with @everywhere
> Now when I read it after the first time it seems to be that he doesnt like 
> to "overwrite" it again?
> Why?
>

[julia-users] Re: Tracking down type instability

[Ralph Smith]

A correction: the native Haswell build does use packed conversions for 
32-bit integers, and valiantly mixes unpacked conversions with packed 
division for 64-bit.  I didn't see any of this from the pre-built binaries.

On Friday, September 16, 2016 at 10:22:16 PM UTC-4, Ralph Smith wrote:
>
> The integer vectors d and l have to be converted to floats before 
> division. Vectorization at this point means that the compiler writes code 
> to use wide registers to do more than one at a time (and lots of extra 
> logic to handle alignment and stragglers).
> Something in the chain (LLVM front end, I suspect) seems ignorant or 
> unhappy about the packed conversion instructions on recent Intel chips 
> (they are not emitted on my build of v0.5rc4 from source, targeted for 
> native Haswell architecture).  What you show should still compile to a 
> pretty fast loop.
>
> Surely this is not your bottleneck?
>
> On Friday, September 16, 2016 at 8:12:53 PM UTC-4, Ben Ward wrote:
>>
>> function 
>> distance{T<:MutationType,A<:NucleotideAlphabet}(::Type{Proportion{T}}, 
>> seqs::Vector{BioSequence{A}})
>> d, l = distance(Count{T}, seqs)
>> D = Vector{Float64}(length(d))
>> @inbounds @simd for i in 1:length(D)
>> D[i] = d[i] / l[i]
>> end
>> return D, l
>> end
>>
>> *julia> **@code_llvm distance(Proportion{AnyMutation}, dnas2)*
>>
>>
>> define %jl_value_t* @julia_distance_70450(%jl_value_t*, %jl_value_t*) #0 {
>>
>> top:
>>
>>   %2 = call %jl_value_t*** @jl_get_ptls_states() #1
>>
>>   %3 = alloca [17 x %jl_value_t*], align 8
>>
>>   %.sub = getelementptr inbounds [17 x %jl_value_t*], [17 x 
>> %jl_value_t*]* %3, i64 0, i64 0
>>
>>   %4 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 0, 
>> i64 2
>>
>>   %5 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 0, 
>> i64 3
>>
>>   %6 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 0, 
>> i64 4
>>
>>   %7 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 0, 
>> i64 5
>>
>>   %8 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 0, 
>> i64 6
>>
>>   %9 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 0, 
>> i64 7
>>
>>   %10 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 
>> 0, i64 8
>>
>>   %11 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 
>> 0, i64 9
>>
>>   %12 = bitcast [17 x %jl_value_t*]* %3 to i64*
>>
>>   %13 = bitcast %jl_value_t** %4 to i8*
>>
>>   call void @llvm.memset.p0i8.i64(i8* %13, i8 0, i64 120, i32 8, i1 false)
>>
>>   store i64 30, i64* %12, align 8
>>
>>   %14 = bitcast %jl_value_t*** %2 to i64*
>>
>>   %15 = load i64, i64* %14, align 8
>>
>>   %16 = getelementptr [17 x %jl_value_t*], [17 x %jl_value_t*]* %3, i64 
>> 0, i64 1
>>
>>   %17 = bitcast %jl_value_t** %16 to i64*
>>
>>   store i64 %15, i64* %17, align 8
>>
>>   store %jl_value_t** %.sub, %jl_value_t*** %2, align 8
>>
>>   %18 = call %jl_value_t* @julia_seqmatrix_70452(%jl_value_t* %1, 
>> %jl_value_t* inttoptr (i64 13078874632 to %jl_value_t*)) #0
>>
>>   store %jl_value_t* %18, %jl_value_t** %4, align 8
>>
>>   store %jl_value_t* %18, %jl_value_t** %5, align 8
>>
>>   %19 = call %jl_value_t* @julia_count_mutations_70454(%jl_value_t* 
>> inttoptr (i64 4687660080 to %jl_value_t*), %jl_value_t* %18) #0
>>
>>   store %jl_value_t* %19, %jl_value_t** %6, align 8
>>
>>   %20 = getelementptr inbounds %jl_value_t, %jl_value_t* %19, i64 0, i32 0
>>
>>   %21 = load %jl_value_t*, %jl_value_t** %20, align 8
>>
>>   store %jl_value_t* %21, %jl_value_t** %7, align 8
>>
>>   %22 = getelementptr inbounds %jl_value_t, %jl_value_t* %19, i64 1, i32 0
>>
>>   %23 = load %jl_value_t*, %jl_value_t** %22, align 8
>>
>>   store %jl_value_t* %23, %jl_value_t** %8, align 8
>>
>>   store %jl_value_t* %21, %jl_value_t** %9, align 8
>>
>>   %24 = getelementptr inbounds %jl_value_t, %jl_value_t* %21, i64 1
>>
>>   %25 = bitcast %jl_value_t* %24 to i64*
>>
>>   %26 = load i64, i64* %25, align 8
>>
>>   %27 = call %jl_value_t* inttoptr (i64 4388816272 to %jl_value_t* 
>> (%jl_value_t*, i64)*)(%jl_value_t* inttoptr (i64 4473018288 to 
>> %jl_value_t*), i64 %26)
>>
>>   store %jl_value_t* %27, %jl_value_t** %10, align 8
>>
>>   store %jl_value_t* %27, %jl_value_t** 

[julia-users] Re: Tracking down type instability

[Ralph Smith]

e_warntype which seems ok:
>>>
>>> *julia> **@code_warntype distance(Proportion{AnyMutation}, dnas2)*
>>>
>>> Variables:
>>>
>>>   #self#::Bio.Var.#distance
>>>
>>>   #unused#::Type{Bio.Var.Proportion{Bio.Var.AnyMutation}}
>>>
>>>   seqs@_3::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}
>>>
>>>   d::Array{Int64,1}
>>>
>>>   l::Array{Int64,1}
>>>
>>>   #temp#@_6::Int64
>>>
>>>   D::Array{Float64,1}
>>>
>>>   #temp#@_8::Int64
>>>
>>>   i::Int64
>>>
>>>   seqs@_10::Array{Bio.Seq.DNANucleotide,2}
>>>
>>>
>>> Body:
>>>
>>>   begin 
>>>
>>>   $(Expr(:inbounds, false))
>>>
>>>   # meta: location /Users/bward/.julia/v0.5/Bio/src/var/distances.jl 
>>> distance 133
>>>
>>>   # meta: location 
>>> /Users/bward/.julia/v0.5/Bio/src/var/mutation_counting.jl count_mutations 
>>> 263
>>>
>>>   seqs@_10::Array{Bio.Seq.DNANucleotide,2} = $(Expr(:invoke, 
>>> LambdaInfo for 
>>> seqmatrix(::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}, 
>>> ::Symbol), :(Bio.Var.seqmatrix), :(seqs@_3), :(:seq)))
>>>
>>>   # meta: pop location
>>>
>>>   # meta: pop location
>>>
>>>   $(Expr(:inbounds, :pop))
>>>
>>>   SSAValue(0) = $(Expr(:invoke, LambdaInfo for 
>>> count_mutations(::Type{Bio.Var.AnyMutation}, 
>>> ::Array{Bio.Seq.DNANucleotide,2}), :(Bio.Var.count_mutations), 
>>> Bio.Var.AnyMutation, :(seqs@_10)))
>>>
>>>   #temp#@_6::Int64 = $(QuoteNode(1))
>>>
>>>   SSAValue(9) = (Base.getfield)(SSAValue(0),1)::Array{Int64,1}
>>>
>>>   SSAValue(10) = (Base.box)(Int64,(Base.add_int)(1,1))
>>>
>>>   d::Array{Int64,1} = SSAValue(9)
>>>
>>>   #temp#@_6::Int64 = SSAValue(10)
>>>
>>>   SSAValue(11) = (Base.getfield)(SSAValue(0),2)::Array{Int64,1}
>>>
>>>   SSAValue(12) = (Base.box)(Int64,(Base.add_int)(2,1))
>>>
>>>   l::Array{Int64,1} = SSAValue(11)
>>>
>>>   #temp#@_6::Int64 = SSAValue(12) # line 256:
>>>
>>>   SSAValue(6) = (Base.arraylen)(d::Array{Int64,1})::Int64
>>>
>>>   D::Array{Float64,1} = 
>>> (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Float64,1)::Type{Array{Float64,1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,Array{Float64,1},0,SSAValue(6),0)::Array{Float64,1}
>>>  
>>> # line 257:
>>>
>>>   $(Expr(:inbounds, true))
>>>
>>>   SSAValue(8) = (Base.arraylen)(D::Array{Float64,1})::Int64
>>>
>>>   SSAValue(13) = 
>>> (Base.select_value)((Base.sle_int)(1,SSAValue(8))::Bool,SSAValue(8),(Base.box)(Int64,(Base.sub_int)(1,1)))::Int64
>>>
>>>   #temp#@_8::Int64 = 1
>>>
>>>   26: 
>>>
>>>   unless (Base.box)(Base.Bool,(Base.not_int)((#temp#@_8::Int64 === 
>>> (Base.box)(Int64,(Base.add_int)(SSAValue(13),1)))::Bool)) goto 37
>>>
>>>   SSAValue(14) = #temp#@_8::Int64
>>>
>>>   SSAValue(15) = (Base.box)(Int64,(Base.add_int)(#temp#@_8::Int64,1))
>>>
>>>   i::Int64 = SSAValue(14)
>>>
>>>   #temp#@_8::Int64 = SSAValue(15) # line 258:
>>>
>>>   SSAValue(5) = 
>>> (Base.box)(Base.Float64,(Base.div_float)((Base.box)(Float64,(Base.sitofp)(Float64,(Base.arrayref)(d::Array{Int64,1},i::Int64)::Int64)),(Base.box)(Float64,(Base.sitofp)(Float64,(Base.arrayref)(l::Array{Int64,1},i::Int64)::Int64
>>>
>>>   
>>> (Base.arrayset)(D::Array{Float64,1},SSAValue(5),i::Int64)::Array{Float64,1}
>>>
>>>   35: 
>>>
>>>   goto 26
>>>
>>>   37: 
>>>
>>>   $(Expr(:inbounds, :pop)) # line 260:
>>>
>>>   return 
>>> (Core.tuple)(D::Array{Float64,1},l::Array{Int64,1})::Tuple{Array{Float64,1},Array{Int64,1}}
>>>
>>>   end::Tuple{Array{Float64,1},Array{Int64,1}}
>>>
>>> But then again the loop in this function is also not vectorised which I 
>>> struggle to see why... unless division can't be.
>>>
>>> On Friday, September 16, 2016 at 3:27:47 AM UTC+1, Ralph Smith wrote:
>>>>
>>>> SSAValue(15) = (Base.getfield)(SSAValue(0),1)
>>>> *::Union{Array{Float64,1},Array

[julia-users] Re: Tracking down type instability

[Ralph Smith]

SSAValue(15) = (Base.getfield)(SSAValue(0),1)
*::Union{Array{Float64,1},Array{Int64,1}}*


indicates that the first element of SSAValue(0) is ambiguous. Earlier it 
shows that this means p from

p, l = distance(Proportion{AnyMutation}, seqs)

which we can't analyze from what you show here.

On Thursday, September 15, 2016 at 10:08:16 AM UTC-4, Ben Ward wrote:
>
> Hi I have two functions and a function which calls them:
>
> @inline function expected_distance(::Type{JukesCantor69}, p::Float64)
> return -0.75 * log(1 - 4 * p / 3)
> end
>
> @inline function variance(::Type{JukesCantor69}, p::Float64, l::Int64)
> return p * (1 - p) / (((1 - 4 * p / 3) ^ 2) * l)
> end
>
> function distance{A<:NucleotideAlphabet}(::Type{JukesCantor69}, 
> seqs::Vector{BioSequence{A}})
> p, l = distance(Proportion{AnyMutation}, seqs)
> D = Vector{Float64}(length(p))
> V = Vector{Float64}(length(p))
> @inbounds for i in 1:length(p)
> D[i] = expected_distance(JukesCantor69, p[i])
> V[i] = variance(JukesCantor69, p[i], l[i])
> end
> return D, V
> end
>
> But I'm seeing type uncertainty:
>
> *@code_warntype distance(JukesCantor69, dnas)*
>
> Variables:
>
>   #self#::Bio.Var.#distance
>
>   #unused#::Type{Bio.Var.JukesCantor69}
>
>   seqs::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}
>
>   p::Array{Float64,1}
>
>   l::Array{Int64,1}
>
>   #temp#@_6::Int64
>
>   D::Array{Float64,1}
>
>   V::Array{Float64,1}
>
>   #temp#@_9::Int64
>
>   i::Int64
>
>
> Body:
>
>   begin 
>
>   SSAValue(0) = $(Expr(:invoke, LambdaInfo for 
> distance(::Type{Bio.Var.Proportion{Bio.Var.AnyMutation}}, 
> ::Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}), 
> :(Bio.Var.distance), Bio.Var.Proportion{Bio.Var.AnyMutation}, :(seqs)))
>
>   #temp#@_6::Int64 = $(QuoteNode(1))
>
>   SSAValue(15) = (Base.getfield)(SSAValue(0),1)
> *::Union{Array{Float64,1},Array{Int64,1}}*
>
>   SSAValue(16) = (Base.box)(Int64,(Base.add_int)(1,1))
>
>   p::Array{Float64,1} = SSAValue(15)
>
>   #temp#@_6::Int64 = SSAValue(16)
>
>   SSAValue(17) = (Base.getfield)(SSAValue(0),2)
> *::Union{Array{Float64,1},Array{Int64,1}}*
>
>   SSAValue(18) = (Base.box)(Int64,(Base.add_int)(2,1))
>
>   l::Array{Int64,1} = SSAValue(17)
>
>   #temp#@_6::Int64 = SSAValue(18) # line 314:
>
>   SSAValue(7) = (Base.arraylen)(p::Array{Float64,1})::Int64
>
>   D::Array{Float64,1} = 
> (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Float64,1)::Type{Array{Float64,1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,Array{Float64,1},0,SSAValue(7),0)::Array{Float64,1}
>  
> # line 315:
>
>   SSAValue(9) = (Base.arraylen)(p::Array{Float64,1})::Int64
>
>   V::Array{Float64,1} = 
> (Core.ccall)(:jl_alloc_array_1d,(Core.apply_type)(Core.Array,Float64,1)::Type{Array{Float64,1}},(Core.svec)(Core.Any,Core.Int)::SimpleVector,Array{Float64,1},0,SSAValue(9),0)::Array{Float64,1}
>  
> # line 316:
>
>   $(Expr(:inbounds, true))
>
>   SSAValue(11) = (Base.arraylen)(p::Array{Float64,1})::Int64
>
>   SSAValue(19) = 
> (Base.select_value)((Base.sle_int)(1,SSAValue(11))::Bool,SSAValue(11),(Base.box)(Int64,(Base.sub_int)(1,1)))::Int64
>
>   #temp#@_9::Int64 = 1
>
>   22: 
>
>   unless (Base.box)(Base.Bool,(Base.not_int)((#temp#@_9::Int64 === 
> (Base.box)(Int64,(Base.add_int)(SSAValue(19),1)))::Bool)) goto 43
>
>   SSAValue(20) = #temp#@_9::Int64
>
>   SSAValue(21) = (Base.box)(Int64,(Base.add_int)(#temp#@_9::Int64,1))
>
>   i::Int64 = SSAValue(20)
>
>   #temp#@_9::Int64 = SSAValue(21) # line 317:
>
>   SSAValue(12) = (Base.arrayref)(p::Array{Float64,1},i::Int64)::Float64
>
>   $(Expr(:inbounds, false))
>
>   # meta: location /Users/bward/.julia/v0.5/Bio/src/var/distances.jl 
> expected_distance 69
>
>   SSAValue(13) = $(Expr(:invoke, LambdaInfo for log(::Float64), 
> :(Bio.Var.log), 
> :((Base.box)(Base.Float64,(Base.sub_float)((Base.box)(Float64,(Base.sitofp)(Float64,1)),(Base.box)(Base.Float64,(Base.div_float)((Base.box)(Base.Float64,(Base.mul_float)((Base.box)(Float64,(Base.sitofp)(Float64,4)),SSAValue(12))),(Base.box)(Float64,(Base.sitofp)(Float64,3)
>
>   # meta: pop location
>
>   $(Expr(:inbounds, :pop))
>
>   SSAValue(5) = 
> (Base.box)(Base.Float64,(Base.mul_float)(-0.75,SSAValue(13)))
>
>   
> (Base.arrayset)(D::Array{Float64,1},SSAValue(5),i::Int64)::Array{Float64,1} 
> # line 318:
>
>   SSAValue(14) = (Base.arrayref)(p::Array{Float64,1},i::Int64)::Float64
>
>   SSAValue(6) = 
> 

Re: [julia-users] Advice on (perhaps) chunking to HDF5

[Ralph Smith]

I have better luck with 

inds = fill(:,3)

By the way, if anyone appropriate is watching, can we have a sticky post 
about how to format Julia code here?
And is the comprehension form of a one-line "for" loop considered good 
style? I don't see it in the manual anywhere.

On Tuesday, September 13, 2016 at 9:36:58 PM UTC-4, sparrowhawker wrote:
>
> Cool! The colons approach makes sense to me, followed by splatting.
>
> I'm unfamiliar with the syntax here but when I try to create a tuple in 
> the REPL using
>
> inds = ((:) for i in 1:3)
>
> I get 
> ERROR: syntax: missing separator in tuple
>
>
>
> On 13 September 2016 at 17:27, Erik Schnetter  > wrote:
>
>> If you have a varying rank, then you should probably use something like 
>> `CartesianIndex` and `CartesianRange` to represent the indices, or possible 
>> tuples of integers. You would then use the splatting operator to create the 
>> indexing instructions:
>>
>> ```Julia
>> indrange = CartesianRange(xyz)
>> dset[indrange..., i] = slicedim
>> ```
>>
>> I don't know whether the expression `indrange...` works as-is, or whether 
>> you have to manually create a tuple of `UnitRange`s.
>>
>> If you want to use colons, then you'd write
>>
>> ```Julia
>> inds = ((:) for i in 1:rank)
>> dset[inds..., i] = xyz
>> ```
>>
>> -erik
>>
>>
>>
>>
>> On Tue, Sep 13, 2016 at 5:08 PM, Anandaroop Ray > > wrote:
>>
>>> Many thanks for your comprehensive recommendations. I think HDF5 views 
>>> are probably what I need to go with - will read up more and then ask.
>>>
>>> What I mean about dimension is rank, really. The shape is always the 
>>> same for all samples. One slice for storage, i.e., one sample, could be 
>>> chunked as dset[:,:,i] or dset[:,:,:,:,i] but always of the form, 
>>> dset[:,...,:i], depending on input to the code at runtime.
>>>
>>> Thanks
>>>
>>> On 13 September 2016 at 14:47, Erik Schnetter >> > wrote:
>>>
 On Tue, Sep 13, 2016 at 11:27 AM, sparrowhawker  wrote:

> Hi,
>
> I'm new to Julia, and have been able to accomplish a lot of what I 
> used to do in Matlab/Fortran, in very little time since I started using 
> Julia in the last three months. Here's my newest stumbling block.
>
> I have a process which creates nsamples within a loop. Each sample 
> takes a long time to compute as there are expensive finite difference 
> operations, which ultimately lead to a sample, say 1 to 10 seconds. I 
> have 
> to store each of the nsamples, and I know the size and dimensions of each 
> of the nsamples (all samples have the same size and dimensions). However, 
> depending on the run time parameters, each sample may be a 32x32 image or 
> perhaps a 64x64x64 voxset with 3 attributes, i.e., a 64x64x64x3 
> hyper-rectangle. To be clear, each sample can be an arbitrary dimension 
> hyper-rectangle, specified at run time.
>
> Obviously, since I don't want to lose computation and want to see 
> incremental progress, I'd like to do incremental saves of these samples 
> on 
> disk, instead of waiting to collect all nsamples at the end. For 
> instance, 
> if I had to store 1000 samples of size 64x64, I thought perhaps I could 
> chunk and save 64x64 slices to an HDF5 file 1000 times. Is this the right 
> approach? If so, here's a prototype program to do so, but it depends on 
> my 
> knowing the number of dimensions of the slice, which is not known until 
> runtime,
>
> using HDF5
>
> filename = "test.h5"
> # open file
> fmode ="w"
> # get a file object
> fid = h5open(filename, fmode)
> # matrix to write in chunks
> B = rand(64,64,1000)
> # figure out its dimensions
> sizeTuple = size(B)
> Ndims = length(sizeTuple)
> # set up to write in chunks of sizeArray
> sizeArray = ones(Int, Ndims)
> [sizeArray[i] = sizeTuple[i] for i in 1:(Ndims-1)] # last value of 
> size array is :...:,1
> # create a dataset models within root
> dset = d_create(fid, "models", datatype(Float64), dataspace(size(B)), 
> "chunk", sizeArray)
> [dset[:,:,i] = slicedim(B, Ndims, i) for i in 1:size(B, Ndims)]
> close(fid)
>
> This works, but the second last line, dset[:,:,i] requires syntax 
> specific to writing a slice of a dimension 3 array - but I don't know the 
> dimensions until run time. Of course I could just write to a flat binary 
> file incrementally, but HDF5.jl could make my life so much simpler!
>

 HDF5 supports "extensible datasets", which were created for use cases 
 such as this one. I don't recall the exact syntax, but if I recall 
 correctly, you can specify one dimension (the first one in C, the last one 
 in Julia) to be extensible, and then you can add more data as you go. You 
 will probably need to specify a chunk size, which 

[julia-users] Re: Workflow question - reloading modules

[Ralph Smith]

I have had good results with the following:  put test code in its own 
module, say TestModule.jl, which uses/imports ModuleA.  Then the REPL 
session looks like this:

include("ModuleA.jl"); include("ModuleB.jl"); include("TestModule.jl")
TestModule.run()

repeated over and over.  Done in this order, each module uses the latest 
version of its prerequisites.  Is there some edge case where this would not 
work?

On Sunday, September 11, 2016 at 4:43:46 PM UTC-4, Chris wrote:
>
> Hi all,
>
> If I'm debugging a single module within the REPL, it is sufficient to just 
> `import Module`, test the things I want to test, make changes to the code 
> in Module, then `reload("Module")`. 
>
> Now, I'm trying to debug code that uses two modules: ModuleA and ModuleB 
> (ModuleA uses ModuleB). Debugging involves changing code in both modules. 
> My usual strategy of `import ModuleA` and `reload("ModuleA")` doesn't work 
> - when I change code in ModuleB and do that, I get a MethodError for 
> ModuleA.function(). So, I have to restart my julia session, which gets very 
> tedious. What should I be doing here?
>
> Thanks.
>

[julia-users] Re: Error in call to eigs--0 index in arpack wrapper

[Ralph Smith]

In particular, the wrappers in v0.4 did not check the error returns from 
Arpack code, so if you passed in a zero starting vector (which is a 
mistake) you would get the unhelpful result you reported.  With v0.5 you 
will get a generic error exception.

The starting vector is used to seed the construction of an orthogonal 
subspace. For tricky problems convergence may vary a lot depending on where 
one starts, so a persistent user is given the chance to try anew if the 
default doesn't work, or to continue iterating from the last residual in 
case convergence was incomplete. In any case, zero is not a valid starting 
residual.

On Sunday, September 11, 2016 at 8:06:59 AM UTC-4, Andreas Noack wrote:
>
> It would be helpful if you could provide a self-contained example. Also, 
> would it be possible to try out the release candidate for 0.5. We have made 
> a few changes to the ARPACK wrappers so it would be useful to know if is 
> only happening on 0.4. Thanks.
>
> On Saturday, September 10, 2016 at 10:41:28 PM UTC-4, Steven White wrote:
>>
>> It only happens when I give an initial vector.  I am also specifying in 
>> LinearMap that it is both symmetric and hermitian, which are both true (the 
>> matrix is also real), but maybe one is not supposed to say hermitian if it 
>> is real?  Here is the call to LinearMap:
>>
>> HLinup = 
>> LinearMap(applyHamup,npnt*nj;ismutating=false,issym=true,ishermitian=true)
>>
>>
>> and the call to eigs:
>>
>> (valup,vecup,nconv,niter,nmult) = eigs(HLinup, nev=Nup, 
>> which=:SR,ritzvec=true,ncv=50,maxiter=2,tol=etol,v0=voup)
>>
>>
>> By the way, I was incorrectly assuming [ns]aupd wasn't modifying its 
>> arguments because there was no ! in its name!
>>
>> -Steve
>>
>> On Saturday, September 10, 2016 at 6:18:01 PM UTC-7, Ralph Smith wrote:
>>>
>>> `ipntr` is modified by the call to `[ns]aupd`, normally to make a sane 
>>> index set (the line you quote is there to initialize memory).  I find that 
>>> the example using `eigs()` at the bottom of the LinearMap README still 
>>> works, and have not encountered trouble in other use of Arpack under Julia. 
>>> Can you provide more specifics of the case where you had trouble?
>>>
>>> On Saturday, September 10, 2016 at 7:08:43 PM UTC-4, Steven White wrote:
>>>>
>>>> I am calling eigs and I am getting an out of bounds, which looks like 
>>>> an obvious mistake in the Arpack wrapper.  
>>>>
>>>> calling eigs
>>>>
>>>> ERROR: LoadError: LoadError: BoundsError: attempt to access 
>>>> 2169-element Array{Float64,1}:
>>>>
>>>>0.0 
>>>>
>>>>
>>>>   at index [0:722]
>>>>
>>>>  in throw_boundserror at abstractarray.jl:156
>>>>
>>>>  in getindex at array.jl:289
>>>>
>>>>  in aupd_wrapper at linalg/arpack.jl:55
>>>>
>>>>  in _eigs at linalg/arnoldi.jl:228
>>>>
>>>>
>>>> In the file arpack.jl I see line 55 is:
>>>>
>>>>
>>>> x = workd[load_idx]
>>>>
>>>>
>>>> Before that, load_idx is defined with:
>>>>
>>>> zernm1 = 0:(n-1)
>>>>
>>>> ...
>>>>
>>>> ipntr  = zeros(BlasInt, (sym && !cmplx) ? 11 : 14)
>>>>
>>>> ...
>>>>
>>>> load_idx = ipntr[1]+zernm1
>>>>
>>>>
>>>> and workd was previously initialized with:
>>>>
>>>> workd  = Array(T, 3*n)
>>>>
>>>>
>>>> So it looks like an attempt to get a section of an array with index
>>>>
>>>> starting at zero???
>>>>
>>>>
>>>> This is on a Mac, Julia version 0.4.6, in the directory:
>>>>
>>>>
>>>> /Applications/Julia-0.4.6.app/Contents/Resources/julia/share/julia/base/linalg
>>>>
>>>>
>>>> If it matters, I was calling eigs with a LinearMap object (Package 
>>>> LinearMaps)
>>>>
>>>

[julia-users] Re: Error in call to eigs--0 index in arpack wrapper

[Ralph Smith]

`ipntr` is modified by the call to `[ns]aupd`, normally to make a sane 
index set (the line you quote is there to initialize memory).  I find that 
the example using `eigs()` at the bottom of the LinearMap README still 
works, and have not encountered trouble in other use of Arpack under Julia. 
Can you provide more specifics of the case where you had trouble?

On Saturday, September 10, 2016 at 7:08:43 PM UTC-4, Steven White wrote:
>
> I am calling eigs and I am getting an out of bounds, which looks like an 
> obvious mistake in the Arpack wrapper.  
>
> calling eigs
>
> ERROR: LoadError: LoadError: BoundsError: attempt to access 2169-element 
> Array{Float64,1}:
>
>0.0 
>
>
>   at index [0:722]
>
>  in throw_boundserror at abstractarray.jl:156
>
>  in getindex at array.jl:289
>
>  in aupd_wrapper at linalg/arpack.jl:55
>
>  in _eigs at linalg/arnoldi.jl:228
>
>
> In the file arpack.jl I see line 55 is:
>
>
> x = workd[load_idx]
>
>
> Before that, load_idx is defined with:
>
> zernm1 = 0:(n-1)
>
> ...
>
> ipntr  = zeros(BlasInt, (sym && !cmplx) ? 11 : 14)
>
> ...
>
> load_idx = ipntr[1]+zernm1
>
>
> and workd was previously initialized with:
>
> workd  = Array(T, 3*n)
>
>
> So it looks like an attempt to get a section of an array with index
>
> starting at zero???
>
>
> This is on a Mac, Julia version 0.4.6, in the directory:
>
>
> /Applications/Julia-0.4.6.app/Contents/Resources/julia/share/julia/base/linalg
>
>
> If it matters, I was calling eigs with a LinearMap object (Package 
> LinearMaps)
>

[julia-users] PyPlot x(y)tick label fontsize

[Ralph Smith]

 ax[:tick_params]("both",labelsize=24)

See http://matplotlib.org/api/pyplot_api.html
for related functions and arguments.

[julia-users] @threads not providing as big speedup as expected

[Ralph Smith]

Dict can be slow. Try

  d_cl = Array{Array,1}
  for i=1:np
d_cl[i] = copy(inv_cl)
  end

I can't say if that leads to good threads  since nthreads() gives me only 1 
today.

[julia-users] Re: Avoiding duplicate code when defining concrete subtypes of an abstract type

[Ralph Smith]

If you don't need the flexible interrelations provided by J. Sarnoff's 
approach, you might use metaprogramming like this:

abstract Cash


function _validate(p::Cash)
p.b >= 100 && throw(ArgumentError("Too much..."))
(p.b < 0 || p.a < 0) && throw(ArgumentError("Cannot have negative 
amounts"))
end


macro make_currency(C)
esc(quote
  type $C <: Cash
a::Int64
b::Int64
function $C(a::Int64,b::Int64)
  this = new(a,b)
  _validate(this)
  this
end
  end
  # other boilerplate specific to C goes here
end)
end


for t = (:Usd,:Gbp,:Whole)
# show(macroexpand(:(@make_currency $t))) # for debugging
@eval @make_currency($t)
end

# specialize for the exceptional case:
function _validate(p::Whole)
p.a < 0 && throw(ArgumentError("Cannot have negative amounts"))
p.b != 0 && throw(ArgumentError("Must be whole amount"))
end






On Wednesday, August 24, 2016 at 1:35:19 PM UTC-4, Cliff wrote:
>
> As an example of what I mean, suppose that I'm trying to make a currency 
> type for various currencies:
>
> abstract cash
>
> type usd <: cash
>   a::Int64
>   b::Int64
>   function usd(a,b)
> b>=100 && error("Too much of smaller currency unit (over 99).")
> (b<0 || a<0) && error("Cannot have negative amounts.")
> new(a,b)
>   end
> end
> type gbp <: cash
>   a::Int64
>   b::Int64
>   function gbp(a,b)
> b>=100 && error("Too much of smaller currency unit (over 99).")
> (b<0 || a<0) && error("Cannot have negative amounts.")
> new(a,b)
>   end
> end
>
> function Base.(:+){T<:cash}(x::T, y::T)
>   total = 100x.a+100y.a+x.b+y.b
>   return T(div(total, 100), total%100)
> end
>
> I'm able to define addition once for any subtype of cash. On the other 
> hand, I have to write out a large amount of almost identical code each time 
> I want to define a new subtype of cash. Is there some way I can define this 
> once in an abstract manner (with any manually specified constructor 
> overriding the general one), so that we could have something like:
>
> abstract cash
> general type T <: cash
>   a::Int64
>   b::Int64
>   function T(a::Int64,b::Int64)
> b>=100 && error("Too much of smaller currency unit (over 99).")
> (b<0 || a<0) && error("Cannot have negative amounts.")
> new(a,b)
>   end
> end
>
> type usd <: cash end
> type gbp <: cash end
> type whole <: cash
>   function whole(a,b)
> a<0 && error("Cannot have negative amounts.")
> b != 0 && error("Must be whole amount.")
> new(a,b)
>   end
> end
>
>
>

[julia-users] Re: ANN: Documenter.jl 0.3

[Ralph Smith]

The make.jl in Documenter/docs has the appropriate incantation.  Try this:

cd(joinpath(Pkg.dir(),"Documenter","docs"))
include("make.jl")

Then point your browser at XXX/Documenter/docs/build/html/index.html (where 
XXX designates your Pkg directory).

The developers seem to have dumped their mkdocs.yml so more work is needed 
to build a local copy of the polished docs.


On Friday, August 19, 2016 at 7:18:37 PM UTC-4, Christoph Ortner wrote:
>
> I want to give this a try but I can't find the example of HTML output, 
> which is supposed to be in test/html? 
>
> Thank you. 
>
>
>
> On Friday, 19 August 2016 22:07:55 UTC+1, Morten Piibeleht wrote:
>>
>> We are happy to announce version 0.3 of Documenter.jl 
>>  – a package for building 
>> the documentation of other Julia packages. Documentation is available at 
>> https://juliadocs.github.io/Documenter.jl/stable/.
>>
>>
>> Compared to the 0.2-branch there should be no visible changes to the user 
>> other than some resolved bugs. However, there have been some changes to the 
>> code base so please don’t hesitate to report any issues on the issue 
>> tracker .
>>
>>
>> The largest bit of news is that experimental *native HTML output* is now 
>> available as an option [*]. This means that instead of relying on MkDocs, 
>> Documenter.jl is able to produce its own complete static HTML site. The 
>> generated site’s style borrows heavily from the Julia manual and is already 
>> being used to build Documenter’s own documentation.
>>
>>
>> We would like to encourage everyone to give the HTML output a try 
>> (including just browsing) and send feedback so that we could iron out any 
>> remaining issues and improve user experience. See the documentation 
>> 
>>  for 
>> more information on how to enable HTML output on your own repository. There 
>> is the meta-issue #212 
>>  where we are 
>> tracking future improvements of this feature.
>>
>>
>> Finally, a big thanks to everyone who have been testing the package and 
>> submitting issues so far. Keep ‘em coming!
>>
>>
>> Cheers, 
>> Morten & Mike
>>
>>
>> [*] Developed as part of Morten’s Google Summer of Code project 
>> 
>> .
>>
>

[julia-users] Re: where is libmpfr found?

[Ralph Smith]

For v0.5

joinpath(JULIA_HOME,Base.LIBDIR,"julia","libmpfr" * 
dll_extension_if_needed)

should work for standard installations. For v0.4, replace Base.LIBDIR with 
"..", "lib".
If using a system library, one would presumably need to query the system 
loader to get a reliable value.

On Monday, August 15, 2016 at 5:30:48 PM UTC-4, Jeffrey Sarnoff wrote:
>
> What is the generally applicable expression which, from within Julia, 
> gives the full path to libmpfr?
>

[julia-users] Re: Problems with A\b and BigFloat

[Ralph Smith]

The OP wants extremely high precision and indicated that he was willing to 
factor the matrix.  I recommended iterative refinement which converges very 
quickly, and exploits the state-of-the-art direct solvers.  The solvers in 
IterativeSolvers.jl are for a different domain, where the matrix is too 
large or expensive to factor.  To get high accuracy with them generally 
requires tailored preconditioners, which are not "out of the box". In fact 
one usually need a preconditioner to get any convergence with the 
non-symmetric ones for interesting ranks.  (I've been struggling for months 
to find a good preconditioner for an application of GMRES in my work, so 
this is a sore point.)

On Wednesday, August 10, 2016 at 10:56:22 PM UTC-4, Chris Rackauckas wrote:
>
> Yes textbook answer is, why do you want to use `\`? Iterative techniques 
> are likely better suited for the problem. There's no need to roll you own, 
> the package IterativeSolvers.jl has a good number of techniques implemented 
> which are well-suited for the problem since A is a large sparse matrix. 
> Their methods should work out of the box with Bigs, though you will likely 
> want to adjust the tolerances.
>
> On Wednesday, August 10, 2016 at 7:37:35 PM UTC-7, Ralph Smith wrote:
>>
>> Here is a textbook answer.  Appropriate choice of n depends on condition 
>> of A.
>>
>> """
>>
>> iterimprove(A,b,n=1,verbose=true)
>>
>>  
>>
>> Solve `A x = b` for `x` using iterative improvement 
>>
>> """ 
>>
>> function iterimprove{T<:AbstractFloat}(A::SparseMatrixCSC{T},
>>>b::Vector{T},n=1,verbose=true)
>>>
>>  eps(T) < eps(Float64) || throw(ArgumentError("wrong 
>>> implementation")) 
>>
>>  A0 = SparseMatrixCSC{Float64}(A)
>>> F = factorize(A0)
>>> x = zeros(b)
>>> r = copy(b)
>>> for iter = 1:n+1
>>> y = F \ Vector{Float64}(r)
>>> for i in eachindex(x)
>>> x[i] += y[i]
>>> end
>>> r = b - A * x
>>> if verbose
>>> @printf "at iter %d resnorm = %.3g\n" iter norm(r)
>>> end
>>> end
>>> x
>>> end
>>
>>
>>
>> On Wednesday, August 10, 2016 at 3:47:10 PM UTC-4, Nicklas Andersen wrote:
>>>
>>> Hello
>>>
>>> I'm trying to solve a large, sparse and unsymmetrical linear system Ax = 
>>> b.
>>> For this task I'm using Julias *SparseMatrixCSC *type for the 
>>> definition of my matrices and Julias built in backslash ' \ ' operator for 
>>> the solution of the system.
>>> I need *quadruple precision* and thus I've been trying to implement my 
>>> routine with the *BigFloat *type together with the SparseMatrixCSC type.
>>>
>>> To illustrate this, I give a simple example here:
>>> set_bigfloat_precision(128);
>>> A  = speye(BigFloat, 2, 2);
>>> b = ones(BigFloat, 2, 1);
>>> x = A\b;
>>>
>>> If I do this I either get a StackOverFlow error:
>>> ERROR: StackOverflowError:
>>>  in copy at array.jl:100
>>>  in float at sparse/sparsematrix.jl:234
>>>  in call at essentials.jl:57 (repeats 254 times)
>>>
>>> or the solver seems to run forever and never terminates. As the second 
>>> error indicates it seems like the sparse solver only accepts the normal 
>>> *float* types.
>>> My question is then, is there a way to get quadruple precision with the 
>>> standard solvers in Julia, in this case UMFpack I assume ? or should I look 
>>> for something else (in this case any suggestions :) ) ?
>>>
>>> Regards Nicklas A.
>>>
>>>

[julia-users] Re: Problems with A\b and BigFloat

[Ralph Smith]

Here is a textbook answer.  Appropriate choice of n depends on condition of 
A.

"""

iterimprove(A,b,n=1,verbose=true)

 

Solve `A x = b` for `x` using iterative improvement 

""" 

function iterimprove{T<:AbstractFloat}(A::SparseMatrixCSC{T},
>b::Vector{T},n=1,verbose=true)
>
 eps(T) < eps(Float64) || throw(ArgumentError("wrong implementation")) 

 A0 = SparseMatrixCSC{Float64}(A)
> F = factorize(A0)
> x = zeros(b)
> r = copy(b)
> for iter = 1:n+1
> y = F \ Vector{Float64}(r)
> for i in eachindex(x)
> x[i] += y[i]
> end
> r = b - A * x
> if verbose
> @printf "at iter %d resnorm = %.3g\n" iter norm(r)
> end
> end
> x
> end



On Wednesday, August 10, 2016 at 3:47:10 PM UTC-4, Nicklas Andersen wrote:
>
> Hello
>
> I'm trying to solve a large, sparse and unsymmetrical linear system Ax = b.
> For this task I'm using Julias *SparseMatrixCSC *type for the definition 
> of my matrices and Julias built in backslash ' \ ' operator for the 
> solution of the system.
> I need *quadruple precision* and thus I've been trying to implement my 
> routine with the *BigFloat *type together with the SparseMatrixCSC type.
>
> To illustrate this, I give a simple example here:
> set_bigfloat_precision(128);
> A  = speye(BigFloat, 2, 2);
> b = ones(BigFloat, 2, 1);
> x = A\b;
>
> If I do this I either get a StackOverFlow error:
> ERROR: StackOverflowError:
>  in copy at array.jl:100
>  in float at sparse/sparsematrix.jl:234
>  in call at essentials.jl:57 (repeats 254 times)
>
> or the solver seems to run forever and never terminates. As the second 
> error indicates it seems like the sparse solver only accepts the normal 
> *float* types.
> My question is then, is there a way to get quadruple precision with the 
> standard solvers in Julia, in this case UMFpack I assume ? or should I look 
> for something else (in this case any suggestions :) ) ?
>
> Regards Nicklas A.
>
>

[julia-users] Re: Save a REPL defined function to a file

[Ralph Smith]

Your definition is entered into your history file, probably when you end 
the session.  On Linux, the file is ${HOME}/.julia_history - perhaps 
someone else can report on OSX or Windows.  This does not appear to be 
mentioned in the main documentation.

On Monday, August 8, 2016 at 6:56:30 PM UTC-4, Naelson Douglas wrote:
>
> Is there a way for me to take a function I defined on REPL and then save 
> it to a file?
>
> Example
> julia > double(x) = 2*x
> double (generic function with 1 method)
> julia > savetofile(double, "myfile.jl")
>
> After that I would open myfile.jl and see it's writher "double(x) = 2*x" 
> there
>
>
>

Re: [julia-users] Re: eigs fails silently

[Ralph Smith]

Right, ARPACK only enforces the projection for the generalized problem. I 
think your example relies on roundoff error in the orthogonalization step 
to populate the null space.

With regard to semi-definite B (or C in the OP), the shifted version of the 
algorithm (mode 3 in ARPACK) will work for this case, but does require ncv 
<= rank(B). Here one needs to solve (A - sigma * B) x = y instead of B x = 
y, with only a nonsingularity requirement.  I got this to work by 
suppressing the isposdef() check in arnoldi.jl; passing sigma=0 selects the 
right mode.  The disadvantage is that when B is *indefinite* the algorithm 
produces nonsense without complaint. Do you know of a simple test for 
semi-definiteness?

Just for the record, the logic in arnoldi.jl is not tricky.  If one has a 
really large problem where the factorizations are impractical, one could 
readily adapt the arnoldi.jl code to invoke iterative linear solvers.

On Saturday, August 6, 2016 at 6:07:58 PM UTC-4, Andreas Noack wrote:
>
> Ralph, ARPACK seems to be okay with the singularity (of A) in the 
> non-generalized case. In the problem here B=I so it's only a generalized 
> problem by algorithm, not by math and 
>
> julia> eigs(A)[1]
> 6-element Array{Float64,1}:
>  -100.0
>   100.0
> 4.44392e-15
>-3.0281e-18
> 3.5842e-33
> 1.69212e-33
>
> so here ARPACK is doing something different. Maybe expanding with vectors 
> from the null space which I guess is fine as long the αs and βs are zeroed 
> and you reorthogonalize.
>
> On Saturday, August 6, 2016 at 5:02:42 PM UTC-4, Ralph Smith wrote:
>>
>> The 0.5 behavior is the same as Fortran, and is technically correct. 
>>  ARPACK cannot create more (orthogonal) Krylov vectors than the rank of the 
>> matrix A,
>> so your example needs the additional keyword ncv=2. This doesn't seem to 
>> be adequately documented in arpack-ng either. The 0.4 behavior looks 
>> unreliable.
>>
>> On Saturday, August 6, 2016 at 2:03:21 AM UTC-4, Madeleine Udell wrote:
>>>
>>> Actually: I'm not sure we should chalk the error up to ARPACK. Julia 0.4 
>>> is able to produce a (correct, I think) answer to
>>>
>>> eigs(A, C)
>>>
>>> but 0.5 gives an ARPACK error. I don't suppose ARPACK changed between 
>>> Julia versions...!?
>>>
>>> On Sat, Aug 6, 2016 at 1:54 AM, Madeleine Udell <madelei...@gmail.com> 
>>> wrote:
>>>
>>>> Andreas, thanks for the investigation. I'll use 0.5 for now, and hope 
>>>> the real problems I encounter are within the capabilities of ARPACK.
>>>>
>>>> It's embarrassing to be bested by Matlab...
>>>>
>>>> On Fri, Aug 5, 2016 at 9:23 PM, Andreas Noack <andreasno...@gmail.com> 
>>>> wrote:
>>>>
>>>>> I've looked a bit into this. I believe there is a bug in the Julia 
>>>>> wrappers on 0.4. The good news is that the bug appears to be fixed on 
>>>>> 0.5. 
>>>>> The bad news is the ARPACK seems to have a hard time with the problem. I 
>>>>> get
>>>>>
>>>>> julia> eigs(A,C,nev = 1, which = :LR)[1]
>>>>> ERROR: Base.LinAlg.ARPACKException("unspecified ARPACK error: -")
>>>>>  in aupd_wrapper(::Type{T}, 
>>>>> ::Base.LinAlg.#matvecA!#67{Array{Float64,2}}, ::Base.LinAlg.##64#71{Array
>>>>> {Float64,2}}, ::Base.LinAlg.##65#72, ::Int64, ::Bool, ::Bool, 
>>>>> ::String, ::Int64, ::Int64, ::String, :
>>>>> :Float64, ::Int64, ::Int64, ::Array{Float64,1}) at 
>>>>> ./linalg/arpack.jl:53
>>>>>  in #_eigs#60(::Int64, ::Int64, ::Symbol, ::Float64, ::Int64, ::Void, 
>>>>> ::Array{Float64,1}, ::Bool, ::B
>>>>> ase.LinAlg.#_eigs, ::Array{Float64,2}, ::Array{Float64,2}) at 
>>>>> ./linalg/arnoldi.jl:271
>>>>>  in (::Base.LinAlg.#kw##_eigs)(::Array{Any,1}, ::Base.LinAlg.#_eigs, 
>>>>> ::Array{Float64,2}, ::Array{Floa
>>>>> t64,2}) at ./:0
>>>>>  in #eigs#54(::Array{Any,1}, ::Function, ::Array{Float64,2}, 
>>>>> ::Array{Float64,2}) at ./linalg/arnoldi.
>>>>> jl:80
>>>>>  in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, 
>>>>> ::Array{Float64,2}, ::Array{Float6
>>>>> 4,2}) at ./:0
>>>>>
>>>>> and since SciPy ends up with the same conclusion I conclude that the 
>>>>> issue is ARPACK. Matlab is doing something else because they are able to 
>>>>> handle thi

Re: [julia-users] Re: eigs fails silently

[Ralph Smith]

The 0.5 behavior is the same as Fortran, and is technically correct. 
 ARPACK cannot create more (orthogonal) Krylov vectors than the rank of the 
matrix A,
so your example needs the additional keyword ncv=2. This doesn't seem to be 
adequately documented in arpack-ng either. The 0.4 behavior looks 
unreliable.

On Saturday, August 6, 2016 at 2:03:21 AM UTC-4, Madeleine Udell wrote:
>
> Actually: I'm not sure we should chalk the error up to ARPACK. Julia 0.4 
> is able to produce a (correct, I think) answer to
>
> eigs(A, C)
>
> but 0.5 gives an ARPACK error. I don't suppose ARPACK changed between 
> Julia versions...!?
>
> On Sat, Aug 6, 2016 at 1:54 AM, Madeleine Udell  > wrote:
>
>> Andreas, thanks for the investigation. I'll use 0.5 for now, and hope the 
>> real problems I encounter are within the capabilities of ARPACK.
>>
>> It's embarrassing to be bested by Matlab...
>>
>> On Fri, Aug 5, 2016 at 9:23 PM, Andreas Noack > > wrote:
>>
>>> I've looked a bit into this. I believe there is a bug in the Julia 
>>> wrappers on 0.4. The good news is that the bug appears to be fixed on 0.5. 
>>> The bad news is the ARPACK seems to have a hard time with the problem. I get
>>>
>>> julia> eigs(A,C,nev = 1, which = :LR)[1]
>>> ERROR: Base.LinAlg.ARPACKException("unspecified ARPACK error: -")
>>>  in aupd_wrapper(::Type{T}, 
>>> ::Base.LinAlg.#matvecA!#67{Array{Float64,2}}, ::Base.LinAlg.##64#71{Array
>>> {Float64,2}}, ::Base.LinAlg.##65#72, ::Int64, ::Bool, ::Bool, ::String, 
>>> ::Int64, ::Int64, ::String, :
>>> :Float64, ::Int64, ::Int64, ::Array{Float64,1}) at ./linalg/arpack.jl:53
>>>  in #_eigs#60(::Int64, ::Int64, ::Symbol, ::Float64, ::Int64, ::Void, 
>>> ::Array{Float64,1}, ::Bool, ::B
>>> ase.LinAlg.#_eigs, ::Array{Float64,2}, ::Array{Float64,2}) at 
>>> ./linalg/arnoldi.jl:271
>>>  in (::Base.LinAlg.#kw##_eigs)(::Array{Any,1}, ::Base.LinAlg.#_eigs, 
>>> ::Array{Float64,2}, ::Array{Floa
>>> t64,2}) at ./:0
>>>  in #eigs#54(::Array{Any,1}, ::Function, ::Array{Float64,2}, 
>>> ::Array{Float64,2}) at ./linalg/arnoldi.
>>> jl:80
>>>  in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs, 
>>> ::Array{Float64,2}, ::Array{Float6
>>> 4,2}) at ./:0
>>>
>>> and since SciPy ends up with the same conclusion I conclude that the 
>>> issue is ARPACK. Matlab is doing something else because they are able to 
>>> handle this problem.
>>>
>>> Given that 0.5 is almost release, I'll not spend more time on the issue 
>>> on 0.4. Thought, if anybody is able to figure out what is going on, please 
>>> let us know.
>>>
>>> On Friday, August 5, 2016 at 8:47:26 AM UTC-4, Madeleine Udell wrote:

 Setting `which=:LR, nev=1` does not return the generalized eigenvalue 
 with the largest real parts, and does not give a warning or error:

 n = 10
 C = eye(n)
 A = zeros(n,n)
 A[1] = 100
 A[end] = -100
 @assert eigs(A, C, nev=1, which=:LR)[1][1] == maximum(eigs(A, C)[1])

 Am I expected to set nev greater than the number of eigenvalues I truly 
 desire, based on my intuition as a numerical analyst? Or has eigs broken 
 its implicit guarantee?


>>
>>
>> -- 
>> Madeleine Udell
>> Assistant Professor, Operations Research and Information Engineering
>> Cornell University
>> https://people.orie.cornell.edu/mru8/
>> (415) 729-4115
>>
>
>
>
> -- 
> Madeleine Udell
> Assistant Professor, Operations Research and Information Engineering
> Cornell University
> https://people.orie.cornell.edu/mru8/
> (415) 729-4115
>