condition on state variable, reset if reach to a threshold

[Ziaeemehr]

I am wondering if it's possible to set a condition on the state variable if it reached a threshold reset it (at the definition of derivative). This is quite common in neural dynamics.

if I try to do this get the following error:

def derivative():
    if # condition on y(0) :
        y(0) = 0.0
    yield  # some code

SyntaxError: cannot assign to function call
``

[Wrzlprmft]

Quick answer: Please have a look at this. If this doesn’t solve your problem, please describe in more detail what you want to achieve.

[Ziaeemehr]

yes, I have seen and used it #35.

For example, if membrane potential reaches a threshold, it needs to be reset to zero at LIF neuron, if I use a condition like the link above, and return 0 for that component,

def derivative():
    if  y(0) > threshold:
        yield 0 # for dv
    else:
        yield 1/tau * ...
    yield ... # for other components

this does not reset the membrane potential (right?) because y_i = y_i + dv * dt y_i still has nonzero value.

is that make sense or is this too obvious and I am missing something?

[Ziaeemehr]

I need this for the specific example of the Montbrio model.

This is the implementation of Heun stochastic in C++:

        derivative(y, k1, t);     // k1 is the slope and being updated in place by reference
        for (size_t i = 0; i < nn; ++i)
                tmp[i] = y[i] + dt * k1[i] + rNoise * rand();
        derivative(tmp, k2, t+dt); // k2 is the slope (derivative) and is updated in place by reference
        for (size_t i = 0; i < nn; ++i)
        {
                y[i] += 0.5 * dt * (k1[i] + k2[i]) + rNoise * rand();
                if (y[i] < 0)    // reset the value of r if negative  (I just simplified the code, it has some more code to seperate the r component)
                    y[i] = 0.0;
        }
                else
                    y[i] += 0.5 * dt * (k1[i] + k2[i]) + vNoise * dist(mt_rand);
        }

I have set the value of one of the components to zero if it becomes negative. this is the main part of the derivative just for completing the code but it does not matter,

    for (size_t i = 0; i < N; ++i)
        {
            double coupling_term = 0;
            for (size_t j = 0; j < adjlist[i].size(); ++j)
            {
                int k = adjlist[i][j];
                coupling_term += adj[i][k] * x[k];
            }
            dxdt[i] = 1.0 / tau * (delta / (tau * M_PI) + 2 * x[i] * x[i + N]);
            dxdt[i + N] = 1.0 / tau * (x[i + N] * x[i + N] + eta + iapp + J * tau * x[i] - (M_PI * M_PI * tau * tau * x[i] * x[i]) + G * coupling_term);
        }

I need to do the same thing in jitcsde.

because I get this behavior in C++ code: and this one from jitcsde:

however it seems r never touches negative values, but this is the only difference I could found that may be the source of the difference.

[Wrzlprmft]

Such discontinuous changes of dynamical variables are not possible in JiTC*DE. I don’t see a way to change this because the integrators, error estimators, etc. inherently rely on continuity for many features.

I see two ways for you to work around this:

• Detect whether the threshold was crossed in Python and reinitialise the integrator accordingly. Mind though that there is no feature to detect an exact threshold crossing.

• Use a binary-switch dynamical variable to approximate your logic. The downside of this is that it may considerably affect your performance. Here is an example that resets a stochastic growth to zero, whenever it reaches the threshold 10. y(1) is the switch. Whenever y(1) exceeds the threshold, it switches from −1 to 1 and only resets to −1, once y(0) is very close to zero:

  import numpy as np
  from symengine import sqrt
  from jitcsde import y, jitcsde
  from jitcxde_common import conditional
  
  big = 1e3
  small = 1e-3
  width = small
  
  threshold = 10
  
  f = [
        1 + 0.3*y(0) + conditional( y(1), 0, 0, -big*y(0) , width ),
        big*(
            - y(1)*(y(1)+1)*(y(1)-1)
            + conditional( y(0), threshold, 0,  1-y(1), width )
            + conditional( y(0),     small, -1-y(1), 0, width )
        ),
  ]
  
  g = [
        sqrt(0.2*y(0)),
        0
    ]
  
  SDE = jitcsde(f,g,verbose=False)
  SDE.set_initial_value([0,0])
  
  for time in np.linspace(0,10,1000):
    print(time,SDE.integrate(time)[0])