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
> >>>>
> >>>
> >>
> >
>
