My understanding is that it's fine to read the text & mathematical description, but you shouldn't look at the code.
Best, --Tim On Saturday, January 30, 2016 10:02:25 AM Isaiah Norton wrote: > > Numerical Recipes in C: The Art of Scientific Computing, 2nd edition. > > Please note that code derived from this book cannot be included in BSD, > LGPL, or GPL licensed libraries (which is to say, most Julia packages). > Distribution is restricted to compiled binaries only, with commercial and > non-commercial license variants. Source redistribution is prohibited. > > Numerical Recipes should not be relied on for any implementations that are > submitted to Julia core or the open source Julia package ecosystem. > > On Sat, Jan 30, 2016 at 8:59 AM, <pokerhontas...@googlemail.com> wrote: > > @Tomas: maybe check out Numerical Recipes in C: The Art of Scientific > > Computing, 2nd edition. There is also an edition for Fortran. The code > > that > > I use in C is basically from there. > > > > @Andrew: The xprime needs to be in the loop. I just made it ones to > > simplify but normally it changes every iteration. (In the DP problem, the > > loop is calculating an expecation and xprime is the possible future value > > of the state variable for each state of the world). Concerning the Dierckx > > package. I don't know about the general behaviour but for my particular > > problem (irregular grid + cubic spline) it is very slow. Run the following > > code: > > > > using Dierckx > > > > spacing=1.5 > > Nxx = 300 > > Naa = 350 > > Nalal = 200 > > sigma = 10 > > NoIter = 10000 > > > > xx=Array(Float64,Nxx) > > xmin = 0.01 > > xmax = 400 > > xx[1] = xmin > > for i=2:Nxx > > > > xx[i] = xx[i-1] + (xmax-xx[i-1])/((Nxx-i+1)^spacing) > > > > end > > > > f_util(c) = c.^(1-sigma)/(1-sigma) > > W=Array(Float64,Nxx,1) > > W[:,1] = f_util(xx) > > > > > > spl = Spline1D(xx,W[:,1]) > > > > function performance2(NoIter::Int64) > > > > W_temp = Array(Float64,Nalal*Naa) > > W_temp2 = Array(Float64,Nalal,Naa) > > xprime=Array(Float64,Nalal,Naa) > > for banana=1:NoIter > > > > xprime=ones(Nalal,Naa) > > W_temp = spl(xprime[:]) > > > > end > > W_temp2 = reshape(W_temp,Nalal,Naa) > > > > end > > > > @time performance2(10000) > > > > 30.878093 seconds (100.01 k allocations: 15.651 GB, 2.19% gc time) > > > > > > > > That's why I went on and asked my friend to help me out in the first > > place. I think the mnspline is really (not saying it's THE fastest) fast > > in doing the interpolation itself (magnitudes faster than MATLAB). But > > then > > I just don't understand how MATLAB can catch up by just looping through > > the > > same operation over and over. Intuitively (maybe I'm wrong) it should be > > somewhat proportional. If my code in Julia is 10 times faster whitin a > > loop, and then I just repeat the operation in that particular loop very > > often, how can it turn out to be only equally fast as MATLAB. Again, the > > mnspline uses all my threads maybe it has something to do with overhead, > > whatever. I don't know, hints appreciated.
