Have you looked at `Base.Cartesian`? This contains macros to e.g. unroll loops, which should give you an additional speed benefit over using explicit loops if the loop size is small.
-erik On Sun, Dec 27, 2015 at 2:49 PM, Cedric St-Jean <[email protected]> wrote: > Hey Stuart, > > I used to do research on whether the Maunder Minimum was caused by the sun > being pushed out of its chaotic basin of attraction, so this was interesting > to me! I don't think you need code generation. Here is the fastest code that > I could come up with, I'm sure that others might make it faster still: > > function explicit_f(X::Vector, N) > mu = X[N+1] > z = zeros(X) > t = 1.0 > @inbounds for i in 1:N > z[i] = mu*X[i]*(1-X[i]) - X[mod(i, N)+1] > t *= (1-2*X[i])*mu > end > z[N+1] = t*t-1 > return z > end > > For N=4, this is twice as slow as your unrolled loop. However, for larger N, > it becomes faster (eg. for N=15 my version is 50% faster). If anyone wants > to try his hand at optimizing my code, here is the N=4 unrolled code: > > function f(X::Vector) > z = zeros(X); > x1 = X[1]; > x2 = X[2]; > x3 = X[3]; > x4 = X[4]; > mu = X[5]; > z[1] = mu*x1*(1-x1) - x2; > z[2] = mu*x2*(1-x2) - x3; > z[3] = mu*x3*(1-x3) - x4; > z[4] = mu*x4*(1-x4) - x1; > mu4 = mu*mu*mu*mu; > xx1 = 1-2*x1; > xx2 = 1-2*x2; > xx3 = 1-2*x3; > xx4 = 1-2*x4; > t = mu4*xx1*xx2*xx3*xx4; > z[5] = t*t-1; > return z; > end; > > > You actually do not have a problem with globals in either version, since > you're passing the array to the function. Note, however, that if you want to > use my version, you should use FastAnonymous like in the Roots.jl > documentation: `@anon x->explicit_f(x, N)`. Julia's native closures are > still slow at the moment. > > Best, > > Cédric > > > On Sunday, December 27, 2015 at 8:45:17 AM UTC-5, Stuart Brorson wrote: >> >> Thanks for everybody's help. I got something working which satisfies >> me. >> >> Regarding Cedric's question, I am trying to compute the Feigenbaum >> number delta to arbitrary high precision. >> >> https://en.wikipedia.org/wiki/Feigenbaum_constants >> >> I am using a brute-force method to compute delta. The basic idea is >> that each bifurcation point may be found by solving a polynomial >> system of order N+1, where N is the number of periods. Delta is then >> found in the limit via a ratio of bifurcation points. I'm using >> Newton's method to solve the system giving the bifurcation points. My >> method is not particularly clever, but it's direct and easily >> understandable for a math-dud like me. >> >> Julia has two features which make this computational method possible: >> 1. I can generate code on the fly using metaprogramming. That's how I >> generate a function implementing my N+1st order system (and the >> associated Jacobean) which I can then call from Newton's method >> at each bifurcation point. 2. Ability to easily use BigFloat as my >> variable type, giving aribtrary precision. I do a lot of work in >> Matlab, and you can't do either of these things in Matlab. >> >> Regarding generating the code, I fiddled around with directly >> generating the AST briefly. Then I decided that it was a little too >> tricky for my knowledge of Julia at this point. I decided that I >> liked the approach where I generate a string, and then parse and eval >> it. The reason is that I can easily inspect and verify a textual >> representation of my program since it is simply human-readable Julia >> code. The AST is one level of indirection away, and I'm not an expert >> enough to be able to read it as easily as Julia code itself. >> >> Here's what I eventually did. This builds the textual representation >> of my polynomial system: >> >> function build_f(N) >> >> mymod(x, y) = mod(x-1, y)+1 >> >> P = string("function f(X::Vector)\n") >> >> P = string(P, "z = zeros(X);\n") >> >> for i in 1:N >> P = string(P, "x$i = X[$i];\n"); >> end >> P = string(P, "mu = X[$(N+1)];\n") >> >> for i in 1:N >> ip1 = mymod(i+1, N) >> P = string(P, "z[$i] = mu*x$i*(1-x$i) - x$ip1;\n") >> end >> >> Q = "mu" >> for i = 2:N >> Q = string(Q, "*mu") >> end >> P = string(P, "mu$N = $Q;\n") >> >> for i in 1:N >> P = string(P, "xx$i = 1-2*x$i;\n"); >> end >> >> Q = "mu$N" >> for i in 1:N >> Q = string(Q, "*xx$i") >> end >> P = string(P, "t = $Q;\n") >> P = string(P, "z[$(N+1)] = t*t-1;\n") >> >> P = string(P, "return z;\n"); >> P = string(P, "end;\n") >> >> return P >> >> end >> >> And this is the result of calling it for N=4: >> >> julia> y = build_f(4) >> "function f(X::Vector)\nz = zeros(X);\nx1 = X[1];\nx2 = X[2];\nx3 = >> X[3];\nx4 = X[4];\nmu = X[5];\nz[1] = mu*x1*(1-x1) - x2;\nz[2] = >> mu*x2*(1-x2) - x3;\nz[3] = mu*x3*(1-x3) - x4;\nz[4] = mu*x4*(1-x4) - >> x1;\nmu4 = mu*mu*mu*mu;\nxx1 = 1-2*x1;\nxx2 = 1-2*x2;\nxx3 = >> 1-2*x3;\nxx4 = 1-2*x4;\nt = mu4*xx1*xx2*xx3*xx4;\nz[5] = >> t*t-1;\nreturn z;\nend;\n" >> >> julia> println(y) >> function f(X::Vector) >> z = zeros(X); >> x1 = X[1]; >> x2 = X[2]; >> x3 = X[3]; >> x4 = X[4]; >> mu = X[5]; >> z[1] = mu*x1*(1-x1) - x2; >> z[2] = mu*x2*(1-x2) - x3; >> z[3] = mu*x3*(1-x3) - x4; >> z[4] = mu*x4*(1-x4) - x1; >> mu4 = mu*mu*mu*mu; >> xx1 = 1-2*x1; >> xx2 = 1-2*x2; >> xx3 = 1-2*x3; >> xx4 = 1-2*x4; >> t = mu4*xx1*xx2*xx3*xx4; >> z[5] = t*t-1; >> return z; >> end; >> >> I'm very pleased I can inspect the code in this format and make sure I >> haven't made any syntax errors. >> >> To solve the system, I pass the fcns to Newton's method like this: >> >> # Create fcn whose roots we want to find. >> y = build_f(N); >> Pf = eval(parse(y)); >> # Build Jacobean of fcn. >> y = build_J(N); >> PJ = eval(parse(y)); >> >> # Do N dimensional Newton's method. >> Xstar = newtonND(f, J, X0, tol); >> >> I probably don't need to assign the result of eval(parse()) to a >> variable Pf since eval(parse()) drops the name of the function into my >> current namespace, right? >> >> The other question is: does this approach create performance-degrading >> global objects? If so, how do I mitigate that problem? >> >> At this point I am struggling with getting Newton's method to converge >> properly for high N. However, that has to do with selecting the >> correct initial condition -- i.e. a typical numerical analysis problem >> -- and is unrelated to the specifics of Julia. What's cool is that I >> have spent only 1 1/2 days on this project, and am quite far along on >> completing it. If I had to stitch everything together using, say, C++ >> or Mathematica, I'd probably still be struggling with the question of >> how to create the system I need to solve. >> >> Stuart >> >> >> >> On Sat, 26 Dec 2015, Cedric St-Jean wrote: >> >> > >> >> >> >> That is, I'm solving a set of N polynomial equations, where N is a >> >> parameter. >> >> >> > >> > That sounds interesting, could you post some more details? At this >> > stage, >> > it's not clear if you really need to generate a new program at all. But >> > if >> > you do want to do that, then for sure, building the AST will be easier >> > and >> > more reliable than generating text. >> > >> > Best, >> > >> > C?dric >> > >> > On Saturday, December 26, 2015 at 3:24:51 PM UTC-5, Stuart Brorson >> > wrote: >> >> >> >> Cool. Thanks for the code example. I will play with it to see what I >> >> can do with it. >> >> >> >> Thanks, >> >> >> >> Stuart >> >> >> >> >> >> On Sat, 26 Dec 2015, Josh Langsfeld wrote: >> >> >> >>> As Yichao is hinting, you may find a macro to be a cleaner way of >> >>> making >> >>> your functions instead of constructing and parsing a string. >> >>> >> >>> macro return_fcn(N) >> >>> xexprs = Expr[:($(symbol(:x,i)) = X[$i]) for i=1:N] >> >>> return esc(:( >> >>> function $(symbol(:f,N))(X::Vector) >> >>> $(xexprs...) >> >>> mu = X[$(N+1)] >> >>> end >> >>> )) >> >>> end >> >>> >> >>> julia> macroexpand(:(@return_fcn(2))) >> >>> :(function f2(X::Vector) >> >>> x1 = X[1] >> >>> x2 = X[2] >> >>> mu = X[3] >> >>> end) >> >>> >> >>> julia> macroexpand(:(@return_fcn(4))) >> >>> :(function f4(X::Vector) >> >>> x1 = X[1] >> >>> x2 = X[2] >> >>> x3 = X[3] >> >>> x4 = X[4] >> >>> mu = X[5] >> >>> end) >> >>> >> >>> >> >>> >> >>> On Saturday, December 26, 2015 at 2:03:56 PM UTC-5, Stuart Brorson >> >> wrote: >> >>>> >> >>>> Julia users, >> >>>> >> >>>> I'm fiddling around with Julia's strings & metaprogramming. I am >> >>>> constructing a function by concatenating a bunch of strings together >> >>>> to create my function, like this: >> >>>> >> >>>> function return_fcn(N) >> >>>> P = string("function f$N(X::Vector)\n") >> >>>> for i in 1:N >> >>>> P = string(P, "x$i = X[$i];\n"); >> >>>> end >> >>>> P = string(P, "mu = X[$(N+1)];\n") >> >>>> etc.... >> >>>> >> >>>> When I execute this code, I get: >> >>>> >> >>>> julia> y = return_fcn(2) >> >>>> "function f2(X::Vector)\nx1 = X[1];\nx2 = X[2];\nmu = X[3];\n" >> >>>> >> >>>> However, what I really want to see is >> >>>> >> >>>> function f2(X::Vector) >> >>>> x1 = X[1]; >> >>>> x2 = X[2]; >> >>>> mu = X[3]; >> >>>> >> >>>> "show(y)" doesn't seem to do what I want. Later, when I do >> >>>> >> >>>> eval(parse(y)) >> >>>> >> >>>> then I get a function which executes correctly. My problem is simply >> >>>> that I can't get Julia to give me a string I can read easily. This >> >>>> will be a very big issue for me when N -> 1024, 2048, etc.... >> >>>> >> >>>> Questions: >> >>>> >> >>>> 1. How can I escape the \n to get a real <CR><LF> in my displayed >> >>>> string? >> >>>> >> >>>> 2. Is this the optimal way to construct a program for later >> >>>> execution >> >>>> (i.e. metaprogramming)? >> >>>> >> >>>> Thanks for all wisdom you have to offer. >> >>>> >> >>>> Stuart >> >>>> >> >>> >> >> >> > -- Erik Schnetter <[email protected]> http://www.perimeterinstitute.ca/personal/eschnetter/
