JiTCODE: How to implement if condition?

Question

I'm trying to solve a dynamic food web with JiTCODE. One aspect of the network is that populations which undergo a threshold are set to zero. So I'm getting a not differentiable equation. Is there a way to implement that in JiTCODE? Another similar problem is a Heaviside dependency of the network.

Example code:

import numpy as np
from jitcode import jitcode, y, t
def f():
    for i in range(N):
        if i <5:
            #if y(N-1) > y(N-2): #Heavyside, how to make the if-statement
                #yield (y(i)*y(N-2))**(0.02)
            #else:
                yield (y(i)*y(N-1))**(0.02)
        else:
            #if y(i) > thr:
                #yield y(i)**(0.2) #?? how to set the population to 0 ??
            #else:
                yield y(i)**(0.3)

N = 10
thr = 0.0001 
initial_value = np.zeros(N)+1

ODE = jitcode(f)
ODE.set_integrator("vode",interpolate=True)
ODE.set_initial_value(initial_value,0.0)        

Answer 1

Python conditionals will be evaluated during the code-generation and not during the simulation (using the generated code). Therefore you cannot use them here. Instead you need to use special conditional objects that provide a reasonably smooth approximation of a step function (or build such a thing yourself):

def f():
    for i in range(N):
        if i<5:
            yield ( y(i)*conditional(y(N-1),y(N-2),y(N-2),y(N-1)) )**0.2
        else:
            yield y(i)**conditional(y(i),thr,0.2,0.3)

For example, you can treat conditional(y(i),thr,0.2,0.3) to be evaluated as 0.2 if y(i)>thr and 0.3 otherwise (at simulation time).

how to set the population to 0 ??

You cannot do such a discontinuous jump within JiTCODE or the framework of differential equations in general. Usually, you would use a sharp population decline to simulate this, possibly introducing a delay (and thus JiTCDDE). If you really need this, you can either:

• Detect threshold crossings after each integration step and reinitialise the integrator with respective initial conditions. If you just want to fully kill populations that went below a reproductive threshold, this seems to be a valid solution.

• Implement a binary-switch dynamical variable.

Also see this GitHub issue.