Could be there a bug in _jitcode.py running jacobian many times?

[erik-stormy]

With not heigh knowledge of python class programming I tried to edit file _jitcode.py to be able to solve large system of ODE. Before that I tried simple example of 3 dimensional Roessler system. And change differentiation of function to actual analytically calculated values. I show the code here:

def _jac_from_f_with_helpers(f, helpers, simplify, n):
    dependent_helpers = [
            find_dependent_helpers(helpers,y(i))
            for i in range(n)
        ]
    a = 0.2; c = 5.7
    Jac = [[ 0, -1, -1], [1, a ,0], [y(2), 0, y(0) - c] ]
#   def line(f_entry):
    def line(p):
        for j in range(n):

            #print('pred diff')         
            entry = Jac[p][j]
            # entry = f_entry.diff(y(j))
            #for helper in dependent_helpers[j]:
                #entry += f_entry.diff(helper[0]) * helper[1]
            if simplify:
                entry = entry.simplify(ratio=1.0)
            print('before entry j+1', j+1)
            print(entry)
            yield entry
    #for f_entry in f():
        #print()
        #print('line(f_entry)')
        #yield line(f_entry)
    for p in range(n):
        print()
        print('line(p) of Jacobian', p)
        yield line(p)

I believe these changes are very simple. When I run this new _jitcode.py with something like:

ODR = [ -y(1)-y(2), y(0)+a*y(1), b+y(2)*(y(0)-c)]
I = jitcode_lyap(ODR, n_lyap = 3 )
I.generate_f_C(simplify=False, do_cse=False)

There are 3 copies of Jacobian. And this changes with length of equation (n to n). So basically I think there is some unnecessary run of _jac_from_f_with_helpers() ? I dont understand the code more, and I know that if it is mistake it would improve run of the compilation for jitcode_lyap.

(I don't know if there is link, but when I was compilling large ODE system, in temp folder I found very large file (hundred of MB) of C code - when compilation crushed, and very high index value of y function)

[Wrzlprmft]

First of all, for every Lyapunov exponent computed with the Benettin method (which JiTCODE uses), you need to solve another set of differential equations as large as the original one. So, to calculate all Lyapunov exponents of your Rössler system, you need to solve a twelve-dimensional ODE.

Now, there are two ways to obtain the Jacobian for this:

• You generate it once, store it and re-use it. This can be problematic for large differential equations, since the Jacobian can be very large.
• You do not store the Jacobian at all, but use generator functions. This is what JiTCODE does. However this means, you need to compute the Jacobian from stretch, every time you use it.

I chose the latter since it does not make a difference for small system and for large systems you rarely ever want to compute more than one Lyapunov exponent. There is room for optimisation once you want two or more Lyapunov exponents for a large system, but it requires quite a lot of restructuring of the code as the Lyapunov functionality builds upon normal JiTCODE. Do you try to compute more than one Lyapunov exponent for your large system and if yes, do you really need that?

Finally, as for the Megabytes of C code: Yes, this is how large things get and not really avoidable. My record is 60 MB for a system with 40000 dynamical variables, IIRC.