neurophysik / jitcdde Goto Github PK

I have installed jitcdde, chspy, and jitcxde_common - all from git source in Arch Linux.

Using jitcdde works as expected. However importing pylab from matplotlib leads to an error "tcmalloc.cc:333] Attempt to free invalid pointer"

minimal working example:

In [1]: from jitcdde import y, t; import matplotlib.pyplot as plt
src/tcmalloc.cc:333] Attempt to free invalid pointer 0x560a54c4cb10 
Aborted (core dumped)

Installed versions of all(?) relevant packages:
Kernel: 5.11.15-arch1-2
Python: Python 3.9.3
gcc 10.2.0-6
python-jinja 2.11.3-2
python-scipy 1.6.2-1
python-numpy 1.20.1-1
python-symengine 0.7.2-1
symengine 0.7.0-1
python-chspy-git 1.2.r0.gbd54f91-1
python-jitcdde-git 1.8.0.r4.ge732be5-2
python-jitcxde_common-git 1.5.2.r2.gced3d5d-1

Is there anything to do that would help me/us to locate the source of the error?

Having a dynamical system smoothly or sharply vary one of its parameters over time is not uncommon. The former is usually meant either as colored noise or an input, while the latter is usually a perturbation to the system. As far as I could gather from Open and Closed issues, the docu and examples, currently the intended way to do this jitcdde_input.

It is however is not obvious what tool to use to do a parameter perturbation, if you are not sure what you are looking for. I simply propose adding an example Mackey-Glass system where β has a sharp change in the middle of the integration. If this sounds good, I'd be willing to PR the following example if you think it's worth adding it to the repo. I can't think of a good filename that would point a newcomer to this example if they want to do parameter perturbations, though. I hope this complies with both Python and jitcdde good practices.

### Makey-Glass with time-dependent parameter

import numpy as np
from jitcdde import y, t, jitcdde, jitcdde_input, input
from chspy import CubicHermiteSpline
from matplotlib import pyplot as plt

# Defining parameters
# --------------

τ = 15
n = 10
β0 = 0.25
β1 = 0.4
γ = 0.1

# Using input to modify a parameter mid-run
# --------------

input_times = np.arange(0, 800)
t0 = 400 # time of perturbation

parameter_over_time = lambda t: β0 if t<t0 else β1
parameter_spline = CubicHermiteSpline(n=1)
parameter_spline.from_function(parameter_over_time, times_of_interest = input_times)

# Defining Dynamics
# -----------------

β = input(0)
f = [ β * y(0,t-τ) / (1 + y(0,t-τ)**n) - γ*y(0) ]

# Actual Integration
# ------------------

DDE = jitcdde_input(f, parameter_spline)
DDE.constant_past(1)
DDE.adjust_diff()
times = DDE.t + np.arange(0, input_times[-1])
data = np.vstack([DDE.integrate(time) for time in times])

# Plotting
# --------

fig, (ax1, ax2) = plt.subplots(2,1, sharex=True)

parameter_spline.plot(ax1)
ax1.set_title('Parameter over time')
ax1.set_ylabel('β')

ax2.plot(times, data, label='perturbed')
ax2.set_xlabel('time')

Output:

Suppose my model equations are

f = [
y(0)*(r0-r2*y(3)-r1*y(0))-alpha*y(0)*y(1)/(1+e1*y(1)),
y(1)*(s0-s2*y(2)-s1*y(1))+(g*alpha*y(0)*y(1))/(1+e1*y(1)),
-delta0*y(2)+sigma*y(1,t-tau)*y(3,t-tau)-beta*y(1)*y(2),
varphi-delta1*y(3)-sigma*y(1)*y(3)+nu*beta*y(1)*y(2),
]

We want to solve the equation with a history say (10,10,10,10) in [0,T].
When time say t1 in [0,T] is of the form n*k where k is fixed and n in non-negative integer, I want to update the solutions obtained by solver at the time just before n*k by the following:

  y(0)=y(0)+(1-p)y(0)
  y(1)=y(1)
  y(2)=y(2)
  y(3)=y(3)

i.e. the change is only in y(0) but all other remains as it is obtained from the solver. The solution to the DDE again will start with the newly changed variable at t1=n*k to next time t2=(n+1)*k where there again changes in the solution by the same above equations.

This solution process will go till final time T for each n=0,1,2,3, ......

The detailed problems are available in the attached file.
impulses.pdf

I was trying to do a small parameter search varying the delay, so I decided to use it as a control parameter in my system and change it using set_parameters instead of compiling a new model every time. The following pieces of code are not the actual code I use, but a minimal example that reproduces the same behavior.

The first naïve approach to the problem is just compiling a model without extra keywords and letting the library figure out what delays the system has:

import numpy as np
from jitcdde import jitcdde, t, y
from symengine import Symbol

tau = Symbol('tau')

b = 20
f = [- y(0) + b / (1 + y(0, t-tau)**2 )]

taus = [ 5, 10, 20]
DDE = jitcdde(f, control_pars=[tau])

for t in taus:
    DDE.constant_past([1])
    
    DDE.set_parameters(t)
    DDE.step_on_discontinuities()
    
    times = DDE.t + np.arange(500)
    X = np.vstack( [DDE.integrate(tt) for tt in times] )
    DDE.purge_past()

This raises the exception

ValueError: At least one delay depends on time or dynamics; cannot automatically determine steps.

on the line DDE.step_on_discontinuities(), which is fine. I understand that maybe set_parameters doesn't trigger the automatic delay calculations again, so you need to manually specify them. Since every run has a different delay, I though the smart thing was to manually set the attribute each cycle like so:

DDE = jitcdde(f, control_pars=[tau])

for t in taus:
    DDE.delays = [0, t]
    DDE.constant_past([1])
    
    DDE.set_parameters(t)
    DDE.step_on_discontinuities()
    
    times = DDE.t + np.arange(500)
    X = np.vstack( [DDE.integrate(tt) for tt in times] )
    DDE.purge_past()

(I know jitcdde appends a 0 to the delays, so I included that, but removing it doesn't alter the following behavior). Running this results in a crash of some sort: running on console it returns Segmentation fault (core dumped) and when trying in an interactive setting like Spyder I get

Kernel died, restarting
Restarting kernel... 

[SpyderKernelApp] WARNING | No such comm: 5538ec1a9be111ecb3f349eb681f0c91

which I don't think is relevant but may be of some help for you.

Even adding a list to the delays on initialization like so results in the same crash: DDE = jitcdde(f, control_pars=[tau], delays=taus).

Regarding the problem itself, it's easily solved if I abandon the idea of having the model use the delay it's actually using in that run and just give it the full list of delays I'm gonna use:

DDE = jitcdde(f, control_pars=[tau], delays=taus)

for t in taus:
    DDE.constant_past([1])
    
    DDE.set_parameters(t)
    DDE.step_on_discontinuities()
    
    times = DDE.t + np.arange(500)
    X = np.vstack( [DDE.integrate(tt) for tt in times] )
    DDE.purge_past()

which I'm not sure is the best call, as it may be worse or slower at integrating than the case where the model knows the actual delay used. I wouldn't expect the fix to allow the line DDE.delays = [0, t], but rather just prevent the user from setting the attribute.

On a separate note, since delays gets appended with an extra 0, and it actually uses the user-passed object, rather than a copy of it, the last bit of code I showed will run for delay values [0, 5, 10, 20]. I'm not sure this will still happen in the unreleased version that changes the very line that appends the 0, but if it does, it'd be reasonable to warn the user in the docs. or make a copy.

jitcdde version: 1.8.0 (I'm eagerly awaiting a new release!)

Dear developer,

this is rather an implementation question than an issue. I am using jitc*de for almost all my simulations (computational neuroscience) and I love the speed! Recently, I wanted to make a leap and rewrite all my models to jitc*de however this one turned out to be particularly challenging.
I do not want to bother you with all the equations, so I just state the necessary:

• each node in the network has some internal dynamics, which can be written as symbolic function acting on state vector y
• each node has also a firing rate variable, which is actually not part of the state vector, i.e. it can be computed in each timestep given other variables from the state vector, i.e. Q = f(y)
• now comes the problem: some variables inside the state vector y depend on the past of firing rate (in the very simple example with one node just the past of the firing rate, in a network actually on the past of firing rates from all other nodes), symbolically something like y=f(y, Q(t-tau))
• this, in turn, makes the firing rate itself dependent on its own past: Q = f(y, Q(t-tau))

is this doable to implement in jitcdde? Without the dependence of Q on Q(t-tau) I would just make Q helper and compute it in each timestep. But, I need to somehow save its state. Then I was thinking about using the "hack" you thaught me: adding one dummy state variable y(n) which would "store" firing rates, but I wasn't able to write anything to it during the integration itself. Is there a way how to add e.g. anchors (as in CHSPy) during integration to dummy state variable?

The second question is simple yes/no: is there a way how to implement stochastic differential equations with delays? Something like mixing jitcdde and jitcsde.

Thanks a lot! And as I said, love the jitc*de :)

Yes, versions:

jitcdde==1.5.0
jitcode==1.4.0
jitcsde==1.4.0
jitcxde-common==1.5.0

Hi:

I have a question relating to whether or not there is a way to speed up the function of setting parameters. I have followed the tutorial and the example code within #30 to setup the DDE for outer parameter optimization. The hopefully relevant part of the model code is as followed:

import numpy as np
from jitcdde import jitcdde,y,t, jitcdde_input, input
import sympy as sp
from ttictoc import tic,toc



IncubeD = 5
RecoverID = 10
RecoverHD = 15
DetectD = 7
maxT = 300
t_delay = 14

S,V,E,EV,I,IV,IQ,ID,IU,R,D,DC = [y(i) for i in range(12)]
IQd = y(6,t-t_delay)
IDd = y(7,t-t_delay)
control_pars = [alpha, r_gov, p_und, p_dth, r_dth, d_decay, p_d, k1, k2, k3] = symbols("alpha rgov pund pdth rdth ddecay pd k1 k2 k3")
gamma_t = (
    sp.Max(1 - r_gov * sp.Max(sp.log(IQd + IDd + 1)/ (np.log(N)), 0), 0)
)            
r_i = np.log(2) / IncubeD  # Rate of infection leaving incubation phase
r_d = np.log(2) / DetectD  # Rate of detection
r_ru = np.log(2) / RecoverID  # Rate of recovery not under infection
r_rq = np.log(2) / RecoverHD  # Rate of recovery under hospitalization
r_decay = 1 / d_decay  
t_predictions = [i for i in range(maxT)]
input_spline = CubicHermiteSpline.from_data(np.array(t_predictions),np.zeros(((len(t_predictions),1)))

model_long_covid = {
    S: -alpha * gamma_t * S * (I + p_und * IU)  / N - input(0) + r_decay * V + r_decay * R,
    V: input(0) - alpha * (1 - beta) * gamma_t * V * (I + p_und * IU) / N - r_decay * V,
    E:  alpha * gamma_t * S * (I + p_und * IU)  / N  - r_i * E,
    EV: alpha * (1 - beta) * gamma_t * V * (I + p_und * IU)  / N - r_i * EV,
    I: r_i * E - r_d * I,
    IV: r_i * EV - r_d * IV,
    IQ: p_d * (r_d * IV  + r_d * (1 - p_dth) * I) - r_rq * IQ,
    ID: r_d * p_dth * I - r_dth * ID,
    IU: (1 - p_d) * (r_d * IV  + r_d * (1 - p_dth) * I) - r_ru * IU,
    R: r_rq * IQ + r_ru * IU - r_decay * R,
    D: r_dth * ID,
    # Helper states (usually important for some kind of output)
    DC: p_d * (r_d * IV  + r_d * (1 - p_dth) * I) + r_d * p_dth * I
    }

DDE = jitcdde_input(model_long_covid, n=12, input = input_spline, control_pars=control_pars, max_delay=t_delay, verbose=False )
DDE.compile_C(simplify=True, do_cse=False, chunk_size=30, verbose=True) 

def residuals_totalcases(params) -> float:
    """
    Function that makes sure the parameters are in the right range during the fitting process and computes
    the loss function depending on the optimizer that has been chosen for this run as a global variable
    :param params: currently fitted values of the parameters during the fitting process
    :return: the value of the loss function as a float that is optimized against (in our case, minimized)
    """
    # Variables Initialization for the ODE system
    params =  tuple(params_validation(params))
    x_0_cases = get_initial_conditions(
        params_fitted=params, global_params_fixed=GLOBAL_PARAMS_FIXED
    )
    # Force params values to stay in a certain range during the optimization process with re-initializations
    tic()
    DDE.purge_past()
    print(f"purge past: {toc()}")
    tic()
    DDE.constant_past(x_0_cases)      
    print(f"setting constant_past: {toc()}")           
    tic()
    DDE.set_parameters(params)
    print(f"setting parameters: {toc()}")
    tic()
    DDE.adjust_diff()
    print(f"adjust diff: {toc()}")
    tic()
    x_sol = np.transpose(np.vstack([ DDE.integrate(time) for time in t_cases ]))
    print(f"actual integration time: {toc()}")
    weights = np.ones(len(cases_data_fit))
    residuals_value = get_residuals_value(
        optimizer=OPTIMIZER,
        balance=balance,
        x_sol=x_sol,
        cases_data_fit=cases_data_fit,
        deaths_data_fit=deaths_data_fit,
        weights=weights,
        balance_total_difference=balance_total_difference 
    )
    return residuals_value

Basically, I designed a wrapper function residuals_totalcases to wrap the dde model model_long_covid for downstream optimization purposes (in which I would not bore you here). My question/problem is that as you can see with my liberal use of the tictoc package, I have been trying to figure out what is taking most of the time as the DDE formulation of these equations scale roughly 10x slower than the equivalent set of equations in which there is no delay term (i.e. setting the constant delay to 0). And my timing experiments show the following (timing of each step in seconds):

As you can see, apparently the slowest step is the step that sets the parameters which take by far the majority of the time. That is somewhat perplexing to me as I would think setting parameters should not take up most of the time of the integration, so I am wondering if I am doing anything wrong and if there is any way to speed the integration up? Thanks a lot!

HI,
I saw many tutorials to learn how to use Jitcdde, but many of them use constant delays. I want to know if is possible to solve variable delay using Jitcdde. Something like this

$$ dy/ dt = y(t) \times ( G(t + \tau) - G(t)) $$

$$ d\tau/ dt = F(t) / F(t + d\tau) $$

Where G and F are functions. For exemplification both could be a $sin$

Please @Wrzlprmft could you give an example of how to implement this type of delay equations using JITCDDE.

Hello,

First of all, I would like to thank you for this package, it's really amazing (even a beginner in programming like me can appreciate the lot of features and the way they are implemented).

I would like also to ask you a really quick question: can a delay parameter become a "control" parameter in the sense specified in the docs? I would like to sweep over it and try different delays.

(Actually, i would like to minimise the difference between my model and some data, fitting the model to them with a minimisation of residuals, but if it's possible to sweep over the delay parameter, I should do it).

I am talking about a "constant" delay, not a state-dependent delay nor a time-dependent delay.

I continue my attempt to apply Jitcdde. Now I encounter a issue. While solving a DDE system with three equations and one constant delay. With a set of parameter values where time delay is 0.2 then the system is showing some math domain error. Can you please give me some suggestion regarding such problem. If you provide me your mail id i will send you the code. Here I am providing the error message.

File "", line 1, in
runfile('C:/Users/user/Downloads/jitcode-master/jitcode-master/examples/untitled1.py', wdir='C:/Users/user/Downloads/jitcode-master/jitcode-master/examples')

File "C:\ProgramData\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 710, in runfile
execfile(filename, namespace)

File "C:\ProgramData\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 101, in execfile
exec(compile(f.read(), filename, 'exec'), namespace)

File "C:/Users/user/Downloads/jitcode-master/jitcode-master/examples/untitled1.py", line 22, in
data.append( DDE.integrate(time) )

File "C:\ProgramData\Anaconda3\lib\site-packages\jitcdde_jitcdde.py", line 789, in integrate
self.DDE.get_next_step(self.dt)

File "C:\ProgramData\Anaconda3\lib\site-packages\jitcdde_python_core.py", line 225, in get_next_step
k_3 = self.eval_f(self.t + 0.75delta_t, self.y + 0.75delta_t*k_2)

File "C:\ProgramData\Anaconda3\lib\site-packages\jitcdde_python_core.py", line 218, in eval_f
return self.f(t, y, *self.parameters)

File "C:\ProgramData\Anaconda3\lib\site-packages\jitcdde_python_core.py", line 177, in
self.f = lambda *args: np.array(F(*args)).flatten()

File "", line 1, in

ValueError: math domain error

When I am running the example mackey_glass.py, I am getting the error

Traceback (most recent call last):
  File "mackey_glass.py", line 73, in <module>
    DDE.step_on_discontinuities()
  File "/Users/sarojkumar/Library/Python/3.7/lib/python/site-packages/jitcdde/_jitcdde.py", line 904, in step_on_discontinuities
    for delay in self.delays
  File "/Users/sarojkumar/Library/Python/3.7/lib/python/site-packages/jitcdde/_jitcdde.py", line 904, in <listcomp>
    for delay in self.delays
AttributeError: 'RealDouble' object has no attribute 'real_part'

I have gone through a similar problem where you have listed a solution to use

pip3 install symengine==0.3.1.dev1 --force-reinstall

which again shows the following error

Found existing installation: symengine 0.4.0
    Uninstalling symengine-0.4.0:
      Successfully uninstalled symengine-0.4.0
  Running setup.py install for symengine ... error
    ERROR: Complete output from command /Library/Frameworks/Python.framework/Versions/3.7/bin/python3.7 -u -c 'import setuptools, tokenize;__file__='"'"'/private/var/folders/py/fk3c24856kn0b81dvlf3mwl80000gn/T/pip-install-azaxqukr/symengine/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /private/var/folders/py/fk3c24856kn0b81dvlf3mwl80000gn/T/pip-record-g6dym_ro/install-record.txt --single-version-externally-managed --compile:
    ERROR: running install
    running build
    running build_ext
    error: [Errno 2] No such file or directory: 'cmake': 'cmake'
    ----------------------------------------
  Rolling back uninstall of symengine
  Moving to /Users/sarojkumar/Library/Python/3.7/lib/python/site-packages/symengine-0.4.0.dist-info/
   from /Users/sarojkumar/Library/Python/3.7/lib/python/site-packages/~ymengine-0.4.0.dist-info
  Moving to /Users/sarojkumar/Library/Python/3.7/lib/python/site-packages/symengine/
   from /Users/sarojkumar/Library/Python/3.7/lib/python/site-packages/~ymengine
ERROR: Command "/Library/Frameworks/Python.framework/Versions/3.7/bin/python3.7 -u -c 'import setuptools, tokenize;__file__='"'"'/private/var/folders/py/fk3c24856kn0b81dvlf3mwl80000gn/T/pip-install-azaxqukr/symengine/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /private/var/folders/py/fk3c24856kn0b81dvlf3mwl80000gn/T/pip-record-g6dym_ro/install-record.txt --single-version-externally-managed --compile" failed with error code 1 in /private/var/folders/py/fk3c24856kn0b81dvlf3mwl80000gn/T/pip-install-azaxqukr/symengine/

I am using Python 3.7.3 in Mac OSX 10.13.6 with Clang 6.0 (clang-600.0.57.

I'm reaching out about an issue I had when I tried to run the Mackey-Glass example code that you provided (I didn't modify your code at all). I downloaded Visual Studio Build Tools as suggested but I get the following error:

Generating, compiling, and loading C code.
jitced.c
C:\Program Files\Python38\include\pyconfig.h(59): fatal error C1083: Cannot open include file: 'io.h': No such file or directory
error: command 'C:\Program Files (x86)\Microsoft Visual Studio\2019\BuildTools\VC\Tools\MSVC\14.27.29110\bin\HostX86\x64\cl.exe' failed with exit status 2

I am using a Windows machine and jitcdde-1.8.0.

Is there something that I am doing incorrectly? Here is a picture of the packages that I downloaded:

I'd appreciate any help troubleshooting this.

Hello, if shortly the problem is the following - I am trying to implement a mathematical model of induction motor [https://www.mathworks.com/matlabcentral/fileexchange/63759-dynamic-model-of-3-ph-im-in-synchronously-rotating-reference-frame?s_tid=prof_contriblnk]. Based on previous discussions in the topic "problems about time-dependent delay and jitcdde parameter f_sym", I build a piecewise implementation (shown below), but unfortunately it was not functional. You can explain what the problem is, I think it would be useful to many.

from numpy import sin,cos,sqrt,dot,array,angle,linspace
from jitcdde import y,t,jitcdde
import matplotlib.pyplot as plt
import math
plt.rcParams['figure.figsize'] = (10.0, 4.0)
#===============================================================
#Voltage Source
#===============================================================
def volt_source(f,v_rms,t_i):
    omega=2*math.pi*f
    v_r=v_rms*sin(omega*t_i)
    v_y=v_rms*sin(omega*t_i-2*math.pi/3)
    v_b=v_rms*sin(omega*t_i+2*math.pi/3)
    A=array([[1,-0.5,-0.5],[0,math.sqrt(3)/2,-math.sqrt(3)/2]])
    [v_a,v_b]=dot(A,array([v_r,v_y,v_b]))*sqrt(2/3)
    theta=angle(complex(v_a,v_b))
    v_d=v_a*cos(theta)+v_b*sin(theta)
    v_q=-v_a*sin(theta)+v_b*cos(theta)
    return [v_d,v_q,theta]
#===============================================================
#Parameters
#===============================================================
f=50
v_rms=230
R_s=0.435 #Stator resistance
R_r=0.816 #Rotor resistance
L_ls=0.002 #Stator inductance
L_lr=0.002 #Rotor inductance
L_m=0.0693
M=L_m # Mutual inductance
P=2 #Number of Poles
J=0.089 #Inertia
B=0.01 #friction coefficient
T_l=0.0 #Load torque
W_e=2*math.pi*f
L_ss=L_ls+L_m
L_rr=L_lr+L_m
v_d=[]
v_q=[]
#===============================================================
#Initial conditions
#===============================================================
delta_t=1e-4
initial_state=array([.0,.0,.0,.0,.0])
n=1
num_sampl=10000
times=linspace(0,n,n*num_sampl)
for time in times:
    v_d.append(volt_source(f,v_rms,time)[0])
    v_q.append(volt_source(f,v_rms,time)[1])
#===============================================================
# Solution (piecewise)
#===============================================================
result=[]
for i,valu_time in enumerate(times):
   #==============================================================
    #Equations dynamic model induction motor with synchronously rotating reference frame
  #===============================================================
    v_q_t=v_q[i]
    v_d_t=v_d[i]
    equat_dyna_model=[
            (-((y(3)-y(3,t-delta_t))/delta_t)*L_m+v_q_t-W_e*(y(1)*
               (L_ls+L_m)+y(2)*L_m)-R_s*y(0))/(L_ls+L_m),
            (-((y(2)-y(2,t-delta_t))/delta_t)*L_m+v_d_t+W_e*(y(0)*
               (L_ls+L_m)+y(3)*L_m)-R_s*y(1))/(L_ls+L_m),
            (-((y(1)-y(1,t-delta_t))/delta_t)*L_m+(W_e-y(4)*2/P)*
             (y(3)*(L_lr+L_m)+y(0)*L_m)-R_r*y(2))/(L_lr+L_m),
            (-((y(0)-y(0,t-delta_t))/delta_t)*L_m-(W_e-y(4)*2/P)*
             (y(2)*(L_lr+L_m)+y(1)*L_m)-R_r*y(3))/(L_lr+L_m),
            ((((y(0)*(y(1)+y(2))*L_m+y(2)*L_lr-y(1)*(y(0)+y(3))*
                L_m+y(3)*L_lr)*(L_m/(L_lr+L_m))*(2/3*P/2))-T_l)*(1/J)-y(4)*B/J)
            ]
    solut=jitcdde(equat_dyna_model)
    solut.constant_past(initial_state,time=valu_time)
    result.append(solut.integrate(valu_time))
    initial_state=result[-1]
    time=valu_time

Guess this must be on the roadmap already, but just putting it out there since many of us will use Python 3.5.

Hi @Wrzlprmft ,

I continue my attempt to apply Jitcdde. Now I encounter a new issue that puzzles me. I'm trying to solve a DDE system. The system contains two equations and three discrete delays; and nine other parameters. With a set of parameter values, Jitcdde does solve the equation. However, with another of the parameter values, Jitcdde does not. My system is well-posted, that is, the right-hand side of the system is continuously differentiable and thus the solution with continuous initial conditions must exist and unique. Can you point me a direction to solve this issue?

Here is another issue I encountered. With the set of good parameter values that can solve the equation, increasing a delay parameter value does not produce a Hopf bifurcation. By theory and also our theoretical proof, when this delay increases, a Hopf bifurcation should occur for this system. Maybe I used the package incorrectly. Can you let me know to use Jitcdde, what else I should install or whatsoever?

Thank you very much!

Hang

I have a query as to why the jitcdde solution is different than obtained by Mathematica NDSolve when it applies to the same problem with the same set of parameters and the same initial history function. 1a.jpeg and 1b.jpeg appear to be the same except the initial solution abound time=0.


python.txt

Hi, sir. when I use the JiTcdde to get the Lyapunov exponent spectrum varying with the delay, the Lyapaunov exponent spectrum is not smooth, as shown in the diagram on the right. But when I input the time delays one by one, run them individually to save the final result, and then repeat this operation, I get a smooth Lyapunov exponent spectrum with changing time delays, as shown in the diagram on the left. How can I solve this problem? The following is my code

Code 1:

from jitcdde import jitcdde_lyap, y, t
import numpy as np
from symengine import sin, pi

tau = 8.0
a = 1.2
b = 0.15

def ff(t,y,tau,a,b):
   return [y(1, t - tau), y(2), -a * y(1) - a * y(2) + a * sin(2 * pi * b * y(0))]

n_lyap = 6
DDE = jitcdde_lyap(ff(t,y,tau,a,b), n_lyap = n_lyap)
DDE.constant_past( [1.0,0.0,0.0], time=0.0 )
DDE.step_on_discontinuities()
max_step = 0.05
DDE.set_integration_parameters(max_step=max_step, atol=1e-5)
pre_T = 10
DDE.integrate_blindly(pre_T, max_step)
data = []
lyaps = []
weights = []
LEs = []
for time in np.arange(0, 1000, 0.05):
# for time in np.arange(DDE.t, DDE.t+10000, 10):
	state, lyap, weight = DDE.integrate(time)
	data.append(state)
	lyaps.append(lyap)
	weights.append(weight)

lyaps = np.vstack(lyaps)
np.savetxt("LEs_t.dat", lyaps)
for i in range(n_lyap):
	Lyap = np.average(lyaps[400:, i], weights=weights[400:])
	LEs.append(Lyap)
	print(Lyap)

np.savetxt("LEs.dat", LEs)
with open('eye1.dat','ab') as f:
	np.savetxt(f, LEs, delimiter=" ")

Code 2:

from jitcdde import jitcdde_lyap, y, t
import numpy as np
from scipy.stats import sem
from symengine import sin, pi

a = 1.2
b = 0.15

def ff(t,y,tau,a,b):
   return [y(1, t - tau), y(2), -a * y(1) - a * y(2) + a * sin(2 * pi * b * y(0))]

Tau = np.linspace(0.1, 8, 159)

for tau in Tau:
	n_lyap = 6
	DDE = jitcdde_lyap(ff(t, y, tau, a, b), n_lyap=n_lyap)
	DDE.constant_past([1.0, 0.0, 0.0], time=0.0)
	DDE.step_on_discontinuities()
	max_step = 0.05
	DDE.set_integration_parameters(max_step=max_step, atol=1e-5)
	pre_T = 10
	DDE.integrate_blindly(pre_T, max_step)
	data = []
	lyaps = []
	weights = []
	LEs = []
	for time in np.arange(0, 1000, 0.05):
		state, lyap, weight = DDE.integrate(time)
		data.append(state)
		lyaps.append(lyap)
		weights.append(weight)

	lyaps = np.vstack(lyaps)
	np.savetxt("LEs_t.dat", lyaps)
	for i in range(n_lyap):
		Lyap = np.average(lyaps[400:, i], weights=weights[400:])
		LEs.append(Lyap)
		print(Lyap)

	np.savetxt("LEs.dat", LEs)
	with open('eye3.dat', 'ab') as f:
		np.savetxt(f, LEs, delimiter=" ")

Hi, sir. When I use the Jitcdde to solve the the famous Lang-Kobayashi semiconductor laser model with delayed optical feedback.

It doesn't work. Can you help me solve it? The script is following:

from jitcdde import t, y, jitcdde
import numpy as np
import pylab as plt
from symengine import cos, sin, pi, Symbol
#--------------------------the system parameters----------------------------
Gn = 8.4e-13
N0 = 1.4e24
taup = 1.927e-12
k = 6e9
τ = 1.5e-9
taus = 2.04e-9
w0 = 0
a = 3.0
J = 1.098e33

f = [
	0.5 * ( Gn * ( y(2) - N0 ) - 1/taup) * y(0) + k * y(0, t - τ) * cos(w0*τ + y(1) - y(1, t - τ)),
	a * 0.5 * ( Gn * ( y(2) - N0) - 1/taup) - k * y(0, t - τ) / y(0) * sin(w0*τ + y(1) - y(1, t - τ)),
	J - y(2)/taus - Gn * ( y(2) - N0)*(y(0)**2),
]

DDE = jitcdde(f)
DDE.constant_past( [1.3e10, 0.0, 1.9e24], time=0.0 )
DDE.step_on_discontinuities()

dt = 5e-12
times = np.arange(dt, 200e-9, dt)
Y = np.vstack([ DDE.integrate(T) for T in times ])

plt.figure()
plt.plot(Y[10000:50000,0],Y[10000:50000,1], 'r-', lw=0.5)
plt.xlabel('x', fontsize=15)
plt.ylabel('y', fontsize=15)
plt.tick_params(labelsize=15)
# plt.xlim(0, 7)
# plt.ylim(-1.5, 1.5)
plt.title('Phase portrait')
plt.show()

Dear Ansmann,
While I integrate time delayed equation with "tanh()" function it shows the following error: (imported numpy as np)

Traceback (most recent call last):
File "ddex.py", line 12, in
f = [ -ay(0)+b(-ny(0,t-tau)+m(np.tanh(l*y(0,t-tau))))]
AttributeError: 'Mul' object has no attribute 'tanh'

I used "sp.tanh" (sympy as sp) but the state diverges, which is not the case. Please help.

Debabrata Biswas ([email protected])

It is stated that state-dependent delays should work. Is there a minimal example provided?

I would like to use JiTCDDE to solve neutral DDE system of equations. This system represents a physical model of sound production in a recorder (musical instrument). I would like to reproduce results presented in an article by Auvray et al. https://asa.scitation.org/doi/10.1121/1.3672815
The system of equations looks like this:

And is a set of damped harmonic oscillator coupled by non-linear functions.

It would be great to be able to access past values of function derivative although it would lead to necessity of designing initial past values.

Dear developer,

I encountered a weird behaviour of jitcdde depending on the host operating system. These are actually two issues, I have a hunch they might be connected, but maybe they are not.

jitcdde + numba import order

I am codeveloping a framework for whole-brain brain simulations neurolib and we use two backends - Euler forward scheme written in numba-jitted functions and jitcdde. I noticed a weird behaviour on MacOS. When I import jitcdde first and then numba I get weird assertion failed and my ipython kernel instantly dies. This is not a problem on Linux, and also not a problem on MacOS when I do the imports in the opposite order - first numba, then jitcdde.
See screenshot:

On the left side I have a connection to the Ubuntu server - both orders work
On the right side I have my computer (MacOS) - the first order gives assert failed at weird address (/Users/ci/miniconda3) - I do not have this address on my computer, I don't use miniconda

This is not a huge issue, as I can just make sure I first import numba everytime and then jitcdde, it is just a bit annoying

Segfault on Ubuntu

Not sure if these two are connected, but this time it is the other way around. Within the above-mentioned neurolib I've written a bunch of tests, where I integrate the brain model using both backends (jitcdde and numba based). I was working reasonably well for 5-6 months. Yesterday I pushed some changes in the numba backend (not related to jitcdde in any way), but one of my tests fail - related to jitcdde. I tried then to downgrade packages, nothing helps. Even weirder is -- the segfault is not happening on MacOS, only on Ubuntu.

steps to reproduce:

1. have Ubuntu
2. git clone [email protected]:neurolib-dev/neurolib.git
3. cd neurolib
4. pip install -r requirements.txt
5. pip install pytest
6. pip install .
7. pytest tests/multimodel/test_aln.py

See screenshot:

left side again: Ubuntu - segfaults on _jitcdde.py jump function (called when setting past of the state vector with constant and then adjust_diff())
right side: MacOS, tests run without any problems..
week or two ago, everything worked on Linux as well, something changed and I don't know what, but I am getting this behaviour on all Linux systems, and not MacOS.

Thanks a lot for looking into it!
Best,
Nikola

I'm trying to solve a DDE multiple times with different initial conditions (all other parameters stay the same). The code generation part seems to be taking the bulk of the runtime for this and making the whole simulation unnecessarily slow.

I there some way to cache the generated code? What would be the best way to tackle this issue with JiTCDDE?

Python version: 3.11.5
JiTCDDE version: 1.8.1

Is there a way to add a time-dependent constant in the model?
Let's say for some reason I want to add a function of time t as a constant in the mackey_glass_jump example, notice that I added sin(t) at the end of the function, but that would not work. What is the correct syntax for it?

def time_dependent_constant(x):
	return 1.2 * x

τ = 15
n = 10
β = 0.25
γ = 0.1



f = [ β * y(0,t-τ) / (1 + y(0,t-τ)**n) - γ*y(0) + time_dependent_constant(t) ]
DDE = jitcdde(f,verbose=False)

DDE.constant_past([1.0])

DDE.step_on_discontinuities()

for time in numpy.arange(DDE.t, DDE.t+1000, 1):
	print( time, *DDE.integrate(time) )

DDE.jump(1,DDE.t)

for time in numpy.arange(DDE.t, DDE.t+1000, 1):
	print( time, *DDE.integrate(time) )

Dear developer,

is it possible to make a release containing the newest addition to jitcdde, that would be the callbacks? We are using jitcdde as the main integrator in modular part of the neurolib package (which is a package for computational neuroscience) and we would love to put jitcdde into our requirements, but ideally from PyPI... since we depend on callbacks, right now we must install the master version directly from GitHub and this doesn't play nicely if we wanna do the release and upload neurolib to PyPI.

Thanks for considering,
Best,

Nikola

It's Kuramoto model, I defined coupling g=Symbole("g") as parameter.
I wand to pass the coupling to kuramotos_f().
if else statement at main function just check for presence of compiled file.

one more question is that is it possible to pass the adjacency matrix also as a parameter?

import os
import numpy as np
import pylab as plt
from jitcode import jitcode, y
from numpy.random import choice
from numpy.random import uniform
from symengine import sin, Symbol


def kuramotos_f():  # kuramotos_f(g, A)??
    for i in range(n):
        coupling_sum = sum(sin(y(j)-y(i))
                           for j in range(n)
                           if A[i, j])
        yield omega[i] + g / (n-1)*coupling_sum


def main(coupling, time_simulation):
        cfilename = "data/jitced.so"
        if not os.path.exists(cfilename):

            I = jitcode(kuramotos_f, n=n, control_pars=[g])
            I.set_parameters(coupling)
            I.generate_f_C()
            I.compile_C()
            I.save_compiled(overwrite=True, destination=cfilename)
        else:
            I = jitcode(n=n, module_location=cfilename)
            I.set_parameters(coupling)

        I.set_integrator("dopri5", atol=1e-6, rtol=1e-5)
        I.set_initial_value(initial_state, time=0.0)

        times = range(0, time_simulation)
        phases = np.empty((len(times), n))

        for i in range(len(times)):
            phases[i, :] = I.integrate(times[i]) % (2*np.pi)
    
        return times, phases


n = 100
q = 0.2
g = Symbol("g")
time_simulation = 2001

if __name__ == "__main__":

    np.random.seed(1)

    A = choice([1, 0], size=(n, n), p=[q, 1-q])
    initial_state = uniform(0, 2*np.pi, n)
    omega = uniform(-0.5, 0.5, n)
    omega.sort()
    
    coupling = 3.0
    times, phases = main(coupling, time_simulation)

Thank you for any guide.

Might be this is just a simple usage error, a simple python test script throws an error:

Traceback (most recent call last):
  File "test.py", line 14, in <module>
    DDE.step_on_discontinuities()
  File "/home/ms/.local/lib/python3.6/site-packages/jitcdde/_jitcdde.py", line 906, in step_on_discontinuities
    for delay in self.delays
  File "/home/ms/.local/lib/python3.6/site-packages/jitcdde/_jitcdde.py", line 906, in <listcomp>
    for delay in self.delays
AttributeError: 'RealDouble' object has no attribute 'real_part'

my testfile is:

from jitcdde import jitcdde, y, t
import numpy

τ = 15
n = 10
β = 0.25
γ = 0.1

f = [ β * y(0,t-τ) / (1 + y(0,t-τ)**n) - γ*y(0) ]
DDE = jitcdde(f)

DDE.constant_past([1.0])

DDE.step_on_discontinuities()

data = []
for time in numpy.arange(DDE.t, DDE.t+10000, 10):
    data.append( DDE.integrate(time) )
numpy.savetxt("timeseries.dat", data)

And I am using newest versions by compiling from source with:

sudo apt install python3-numpy libgmp-dev  python3-pip cython3
sudo -H pip3 install --upgrade scipy sympy # mit oder ohne macht keinen Unterschied

git clone https://github.com/symengine/symengine.git && \
cd symengine/ && \
cmake . && \
make -j4 && \
sudo make install
# i'm going to install system wide
sudo -H pip3 install --upgrade git+git://github.com/symengine/symengine.py
# grab some coffee
sudo -H pip3 install --upgrade git+git://github.com/neurophysik/jitcxde_common git+git://github.com/neurophysik/jitcode git+git://github.com/neurophysik/jitcdde git+git://github.com/neurophysik/jitcsde

I copy and paste the error/warning:

Generating, compiling, and loading C code.
/Users/alejandro/.local/lib/python3.6/site-packages/jitcxde_common/_jitcxde.py:191: UserWarning: Traceback (most recent call last):
  File "/Users/alejandro/.local/lib/python3.6/site-packages/jitcxde_common/_jitcxde.py", line 189, in _attempt_compilation
    self.compile_C()
  File "/Users/alejandro/.local/lib/python3.6/site-packages/jitcdde/_jitcdde.py", line 515, in compile_C
    chunk_size = chunk_size, # only for OMP
  File "/Users/alejandro/.local/lib/python3.6/site-packages/jitcxde_common/_jitcxde.py", line 142, in _render_template
    codefile.write(template.render(kwargs))
UnicodeEncodeError: ascii codec can't encode character '\u201c' in position 2338: ordinal not in range(128)

  warn(format_exc())
/Users/alejandro/.local/lib/python3.6/site-packages/jitcxde_common/_jitcxde.py:193: UserWarning: 
============================================================
READ ME FIRST
============================================================
Generating compiled integrator failed; resorting to lambdified functions. If you can live with using the Python backend, you can call generate_lambdas to explicitly do this and bypass the compile attempt and error messages. Otherwise, you want to take care of fixing the above errors.
============================================================

============================================================

  warn(line + "READ ME FIRST" + line + "Generating compiled integrator failed; resorting to lambdified functions. If you can live with using the Python backend, you can call generate_lambdas to explicitly do this and bypass the compile attempt and error messages. Otherwise, you want to take care of fixing the above errors." + 2*line)

I suppose that I need to install an specific version/distribution of C.

A heads up: the following is not a well written issue report since I wasn't able to isolate the cause nor provide a minimal example that reproduces the issue.

Part of the use I give to this library is integrating multiple thousands of systems with slightly different parameters values to draw maps in the parameter space (think bifurcation diagrams, for example). This requires that I run the integrator in parallel using the built-in multiprocessing library.

The issue is, I've noticed that under some circumstances (I'm not sure which) jitcdde fails to delete the temp files created in /tmp, sometimes leaving thousands of them behind. I've noticed this issue in at least two different cases so far

1. when Ctrl+C-ing to interrupt the integration sometimes only some of the temp files are left behind
2. when running a large-ish set of parameters (say 30k) eventually execution grinds to a halt and I get the error message -bash: cannot create temp file for here-document: No space left on device, which is supposed to be related to either a full disk or a lack of inodes available, none of which are the case. I do however have around 2k temp files in /tmp created by jitcdde after having run around 6k integrations, so some were deleted.

I don't know if there's a chance that if the integration somehow fails because the system is ill-determined for a specific set of parameter values, the temp files won't be deleted.

I know this is very little info to go off of, but maybe there's a known issue relating to this that can help. If it's any sort of help, I just started noticing this after I started using helpers in my integrations. If there's any other info I can provide that may be of use, just let me know.

jitcdde version: 1.8.1
OS: openSUSE Leap (version 15)

Hi,
As I know there is only Runke kutta 23 method available for DDE.
Is it possible for you to add also Euler maruyama or if you don't have time
give some guide how to add it?

best,
Abolfazl

Hello,

I am getting this error on running the jitcdde

clang: error: the clang compiler does not support '-march=native'
error: command '/usr/bin/clang' failed with exit code 1

I am using the version jitcxde-common 1.5.4

The documentatio for this library encourages the user to give the model builder the delay and max_delay of a system when the system is large to cut down on overhead to the pre-processing times of the model. For this reason, I think the Kuramoto example should include this when defining the model:

I = jitcdde(kuramotos,n=n,verbose=False, max_delay=max(τ), delays=τ)

The line as written above doesn't work, and I'm not sure why; it throws

Traceback (most recent call last):

  File "/home/marcos/Documents/DeltaDimer/untitled5.py", line 30, in <module>
    I = jitcdde(kuramotos,n=n,verbose=False, max_delay=max(τ), delays=τ)

  File "/home/marcos/.local/lib/python3.9/site-packages/jitcdde/_jitcdde.py", line 200, in __init__
    self.delays = delays

  File "/home/marcos/.local/lib/python3.9/site-packages/jitcdde/_jitcdde.py", line 218, in delays
    self._delays.append(0)

AttributeError: 'numpy.ndarray' object has no attribute 'append'

The example should of course include a working version of the line above. Note that on my system adding only the max_delay parameter reduced pre-processing times by 20% (from 1.9 seconds average to around 1.55). I'd expect the addition of delay to make this even faster.

Running mackey_glass_lyap.py in examples folder causes Python crash.
mackey_glass.py is running okay.

Here is the output:

(python3) C:\Users\cyao\Documents\Python\jitcdde\examples>python mackey_glass_ly
ap.py
Generating, compiling, and loading C code.
jitced.c
C:\Users\cyao\AppData\Local\Temp\tmpvmly19ga\jitced.c(706): warning C4244: 'init
ializing': conversion from 'Py_ssize_t' to 'unsigned int', possible loss of data

   Creating library C:\Users\cyao\AppData\Local\Temp\tmpvmly19ga\Release\Users\c
yao\AppData\Local\Temp\tmpvmly19ga\jitced.cp36-win_amd64.lib and object C:\Users
\cyao\AppData\Local\Temp\tmpvmly19ga\Release\Users\cyao\AppData\Local\Temp\tmpvm
ly19ga\jitced.cp36-win_amd64.exp
Generating code
Finished generating code
Using default integration parameters.

Using Python 3.6, Windows 7, and VC++ Build Tool 14.0

Hi,

I want to use jitcdde to solve a dde in a physics problem and I am new to your package (and new to python scientific computing packages like sympy). Thank you for your excellent package first.
I come across following problems

1. My delay is a time-dependent, piecewise function, which means the symbolic expressions of parameter "f_sym" in jitcdde has to change with current time. I am not sure if this package is capable of my demand. I have tried to write a function and return different expression list, then use it as f_sym, like the following code

#y(0):Re[b],y(1):Im[b]
def fdelay(tnow):
    hnew = h
    if tnow-tstart<delta and tnow>tstart:
        if (tnow-tstart)/h > 10:    # adaptive step length
            m = n
            hnew = (tnow-tstart)/m
        else:
            m = np.ceil((tnow-tstart)*n**2/10/delta).astype('int')   
            hnew = 10*delta/n**2
        tvals = np.linspace(0,tnow,m+1)
        yvals = np.array([func.Tdis(tstart+sh,t)*y(0,tstart+sh) for sh in tvals])
    elif tnow-tstart>delta:
        tvals = np.linspace(0,delta,n+1)
        yvals = np.array([func.Tdis(t-delta+sh,t)*y(0,t-delta+sh) for sh in tvals])
    else:
        yvals = [1.0]
    rhr = func.we(t)*y(1,t)
    rhi = -func.we(t)*y(0,t) + 2*C*intg.simps(yvals,None,hnew)
    return [rhr, rhi]

data = []
tlist = []
for tnow in np.arange(tstart, tend, 100):
    DDE = jitcdde(fb)
    DDE.constant_past([re0,im0])
    DDE.step_on_discontinuities()
    re, im = DDE.integrate(tnow)
    data.append( np.array([re, im, re**2+im**2]) )
    tlist.append(tnow)

Considering that list "fb" will change with time "tnow", I try to create a DDE object with different "fb" in each "for" cycle. But it doesn't work. It reports errors (e.g. "SystemError: ..\Objects\moduleobject.c:449: bad argument to internal function"), or sometimes it just computes the first points.
So I want to know if there is a solution to deal with this piecewise function with varying symbolic expressions?

2. I need to determine current time in delay f(t,y), but I am not sure if f_sym can get current time (e.g. DDE.t) when it is called in DDE.integrate. That's why I create new DDE objects in each "for" cycle in above codes.

3. My delay term is actually an integration over a period of past time, like this, "\int_{t-t0}^{t}d\tau,T(\tau,,t)y(\tau)". Now I use scipy.integrate.simps to evaluate the delay, it does give a correct symbolic expression, but can jitcdde offer a better solution to this integral-like delay? or is it better to use sympy to evaluate the integral delay?

4. There are too few and too simple examples in documentation, so I sometimes can't find similar cases for complicated problems, or understand how some interface like f_sym works during the computation.

Thanks a lot in advance!

Wance

Hi,

I want to solve my delay differential equation system, I used ddeint, but the results are not good at all, I want to minimize my objective function based on the results of Model. I have one delay part in 'Id' :

def SEIRD(y, t, N, beta, kappa, gamma, mu, tau):
    S, E, I, R, D = y(t)
    Sd, Ed, Id, Rd, Dd = y(t-tau)

    dSdt = - beta(t) * S * I/N
    dEdt = beta(t) * S * I/N  - kappa * E
    dIdt = kappa * E - (1-mu)* gamma*I - mu*gamma*Id
    dRdt = (1-mu) * gamma*I 
    dDdt = mu*gamma*Id
    return np.array([dSdt, dEdt, dIdt, dRdt, dDdt])

days = 64
tau = 11.5
N0 = 80000000
kappa = 1/3
gamma = 1/10
mu = 0.02
E0 = 280
I0 = 114
R0 = 16
D0 = 0
S0 = N - E0 - I0 - R0 - D0

def Model(days, N,  beta0, beta1, beta2, beta3,  mu, E0, delta):

    def beta(t):
        xx=np.linspace(0,63,64)
        for x in xx[:16]:
            if t < x:
                return beta0
        
        for x in xx[16:22]:
            if t < x:
                return beta1
            
        for x in xx[22:52]:
            if t < x:
                return beta2
            
        for x in xx[52:65]:
            if t <= x or t > 63:
                return beta3
        print(t)

    # Historical initial conditions
    g = lambda t: np.array([S0, E0 , I0 *np.exp((-np.log(0.1)/tau )*(t-t0)), R0, D0])

    tt = np.linspace(0, days-1, days)
    yy = ddeint(SEIRD, g, tt, fargs=( N, beta, kappa, gamma, mu, tau))
    S, E, I, R, D = yy.T

    return t, S, E, I, R, D

I am a mathematician and working on my thesis, and I am forced to solve this part, but my programing language is not good enough. It will be great if you help me to solve it. thanks in advance

Best regards,
Nik

I'm getting a long and cryptic message when trying to compile a module to later save it to a file to make multiple short integrations in parallel. The error in question reads

CompileError: command 'gcc' failed with exit status 1


During handling of the above exception, another exception occurred:

SystemExit: error: command 'gcc' failed with exit status 1

/home/user/anaconda3/envs/dde/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3426: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.
  warn("To exit: use 'exit', 'quit', or Ctrl-D.", stacklevel=1)

with a really long traceback, which I can try to provide if needed. The system in question is 12-D composed of two 6-D coupled systems. Here's the code, which I fear I need to provide in full, since a more minimal example does work (see below). Assume get_maximum_delay is a function defined somewhere that retrieves the maximum value τ and τ*fc will have. You can set maxdelay=3 for testing purposes.

def make_model(parameters_lists):
    symbol_names = 'τ, βH, βC, βD, βN, η, α, σ, λm, λp, κC, κD, κE, fc'
    param_symbols = τ, βH, βC, βD, βN, η, α, σ, λm, λp, κC, κD, κE, fc = symengine.symbols(symbol_names)
    def model():
        h1, c1, d1, e1, n1, s1 = [y(i) for i in range(6)]
        h2, c2, d2, e2, n2, s2 = [y(i+6) for i in range(6)]
        h1d = y(0, t-τ)
        s1d = y(5, t-τ)
        h2d = y(6, t-τ)
        s2d = y(11, t-τ)
        h1c = y(0, t-τ*fc)
        h2c = y(6, t-τ*fc)    
        
        yield -h1 + βH * (1 + α * (σ * s1d)**η ) / ( (1 + (σ * s1d)**η) * (1 + h1d**η) )
        yield -c1 + βC / (1 + h1c**η) + λm * e1 - λp * c1 * d1 - κC * c1 * f2
        yield -d1 + βD + λm * e1 - λp * c1 * d1 - κD * d1 * f2
        yield -e1 - λm * e1 + λp * c1 * d1 - κE * e1 * f2
        yield -f1 + βN - f1 * (κC * c2 + κD * d2 + κE * e2)
        yield -s1 + f1 * (κC * c2 + κD * d2 + κE * e2)
        
        yield -h2 + βH * (1 + α * (σ * s2d)**η ) / ( (1 + (σ * s2d)**η) * (1 + h2d**η) )
        yield -c2 + βC / (1 + h2c**η) + λm * e2 - λp * c2 * d2 - κC * c2 * f1
        yield -d2 + βD + λm * e2 - λp * c2 * d2 - κD * d2 * f1
        yield -e2 - λm * e2 + λp * c2 * d2 - κE * e2 * f1
        yield -f2 + βN - f2 * (κC * c1 + κD * d1 + κE * e1)
        yield -s2 + f2 * (κC * c1 + κD * d1 + κE * e1)
    
    # initialize integrator
    maxdelay = get_maximum_delay(parameters_lists)
    DDE = jitcdde(model, verbose=True, max_delay=maxdelay, n=12, control_pars=param_symbols)
    model_location = DDE.save_compiled('model.so', overwrite=True)
    return model_location

After running this, which raises the error, the program is supposed to launch a bunch of processes to load the model and integrate it in parallel for a bunch of different parameter values.

Just as a sanity check, here's a very minimal example of a 2D system that does work when doing this. Ignore of course the fact that we are not multiprocessing here, since the error arises before calling the integration anyway.

import symengine
import numpy as np
from jitcdde import jitcdde, y, t

d = symengine.symbols('d')

def f():
    u = y(0)
    v = y(1)
    ud = y(0, t-d)
    vd = y(1, t-d)
    
    yield 0.5 * u * (1 - vd)
    yield -0.5 * v * (1 - ud)

DDE = jitcdde(f, control_pars=[d], max_delay=1, n=2)
file = DDE.save_compiled('testDDE', overwrite=True)

# this would be in a separate function in the actual case
for this_d in [0.1, 0.2, 0.3, 0.4]:
    I = jitcdde(module_location='testDDE.so', n=2, max_delay=1)
    I.constant_past([1,2])
    I.set_parameters(this_d)
    I.step_on_discontinuities()
    
    tt = I.t + np.arange(0, 30, 0.5)
    data = np.vstack( [I.integrate(this_t) for this_t in tt] )
    

Is there any way to use JITCDDE for DDEs that cannot be expressed symbolically? I am training a neural network using Tensorflow+Keras to learn a DDE, and I want to use a DDE solver to numerically integrate this learned DDE. However, Tensorflow+Keras can not take in symbolic inputs.

Hi! I'm trying to do a Particle Swarm Optimization (PSO) on DDE using jitcdde. Here is my code:

from pyswarm import pso
from jitcdde import jitcdde, y, t
import numpy as np
import matplotlib.pyplot as plt

def dde(τ, β, γ):
    f = [ γ * y(0) * (1 - y(0, t-τ)/β) ]
    DDE = jitcdde(f, verbose=False)

    DDE.add_past_point(-1.0, [.3], [0.0])
    DDE.add_past_point( 0.0, [.3], [0.0])
    DDE.step_on_discontinuities()

    DDE.set_integration_parameters(atol=1e-10,  pws_atol=1e-10)
    data = []
    for time in np.arange(0, 2000, 100):
        data.append( DDE.integrate(time) )
    cleandata = np.array(data).T[0]
    cleandata = cleandata[~np.isnan(cleandata)]
    time = np.arange(0, 2000, 100)
    print(DDE.t)
    plt.plot(time[:len(cleandata)], cleandata, 'o-')

    return cleandata

def object(x):
    x1 = x[0]
    x2 = x[1]
    x3 = x[2]
    x_res = dde(x1, x2, x3)
    y_trunc = y_data[0:len(x_res)]
    return np.sum(x_res - y_trunc) ** 2

x = dde(0.3, 500., 6.)
y_data = list(x)

lb = [0.2, 495., 4.]
ub = [0.4, 505., 8.]

xopt, fopt = pso(object, lb, ub, swarmsize=20, debug=True)

However, after few iterations, it is very easy to throw an error message:

<ipython-input-13-3fae17b8e6a0> in dde(τ, β, γ)
     24     data = []
     25     for time in np.arange(0, 2000, 100):
---> 26         data.append( DDE.integrate(time) )
     27     cleandata = np.array(data).T[0]
     28     cleandata = cleandata[~np.isnan(cleandata)]

/Users/HangYao/anaconda/envs/python3/lib/python3.5/site-packages/jitcdde/_jitcdde.py in integrate(self, target_time)
    704                         self.successful = False
    705                         if self.do_raise_exception:
--> 706                                 raise error
    707                         else:
    708                                 warn(str(error))

/Users/HangYao/anaconda/envs/python3/lib/python3.5/site-packages/jitcdde/_jitcdde.py in integrate(self, target_time)
    699                                                         continue
    700 
--> 701                                         self._adjust_step_size()
    702 
    703                 except UnsuccessfulIntegration as error:

/Users/HangYao/anaconda/envs/python3/lib/python3.5/site-packages/jitcdde/_jitcdde.py in _adjust_step_size(self)
    629                 if p > self.decrease_threshold:
    630                         self.dt *= max(self.safety_factor*p**(-1/self.q), self.min_factor)
--> 631                         self._control_for_min_step()
    632                 else:
    633                         self.successful = True

/Users/HangYao/anaconda/envs/python3/lib/python3.5/site-packages/jitcdde/_jitcdde.py in _control_for_min_step(self)
    615                                 "• The DDE is ill-posed or stiff.\n"
    616                                 "• You did not allow for an absolute error tolerance (atol) though your DDE calls for it. Even a very small absolute tolerance (1e-16) may sometimes help."
--> 617 				% (self.atol, self.rtol, self.min_step))
    618 
    619         def _increase_chance(self, new_dt):

UnsuccessfulIntegration: 
Could not integrate with the given tolerance parameters:

atol: 0.000000e+00
rtol: 1.000000e-05
min_step: 1.000000e-10

The most likely reasons for this are:
• You did not sufficiently address initial discontinuities.
• The DDE is ill-posed or stiff.
• You did not allow for an absolute error tolerance (atol) though your DDE calls for it. Even a very small absolute tolerance (1e-16) may sometimes help.

Is that something I can avoid by adjusting parameters, such as set_integration_parameters? Or is it something I can't avoid at all? Please advise. Thank you!

Great project, I am excited to see the final release!

By the way, do you plan to implement also (continuously) distributed delays?

My problem is the following:

1. I want to solve a system of DDEs say for example 2D prey-predator model with some given initial function in interval [0, Tf].

2. Next, I want to collect all solutions in a text file having three columns(in this case) where the 1st column is time steps, 2nd for prey population and 3rd for the predator population. Given examples only gives the dependent states eg. prey and predator. Although we know the time steps as we have given in for loop.

3. my second problem is for bifurcation diagram, where I want to collect all the maximum points in consecutive three points taken for each iterate and for each dependent variable in whole time sapce and at a partcular value of parameter 'p' i.e.

say y1(t_i) represent the solution of the system at ith time step.
y1'=f1(y1(t),y2(t),y1(t-tau1),y2(t-tau2),p)
y2'=f2(y1(t),y2(t),y1(t-tau1),y2(t-tau2),p)
at p, then program should do the following

For p=0.1:0.1: 10
solve y1 and y2
For each i from 600 to 1000
    a1=y1(i), a2=y1(i+1), a3=y1(i+2)
    if a1<a2 & a2>a3, then
          c=a2
       save p with c in file with 1st column as p and second column as c
     end if
  end for
end for

Finally, my second file should contain the 1st column as p and the second column as values c i.e. y1 which is maximum in the selected three points. If there is no maximum for given i, it should go to the next value of i till the last point without printing p and c.

4. The solver given by Thompson and Shampine has event locator facility inbuilt. Is there anything like this in jitcxde where it can be used for impulsive differential equation solution.

Hi, Sir. Here I am studying a laser system model, I want to calculate the Kolmogorov–Sinai entropy, which is the sum of all positive Lyapunov exponents of the chaotic system. I have no idea how to calculate the 30 Lyapunov exponents since the model with delay time. Can you help me get the first 30 Lyapunov exponents? The initial value is [0.01 0.01 0.01]

Hi,

I ran a simple test using an example code from (https://computationalmindset.com/en/mathematics/solving-delay-differential-equations-in-python-using-numerical-methods.html). I got the results, however, I also got a set of warnings related to the type of pointers. I did not find any posts or issues about these warnings. I ran the example on MacOS, jitcdde version = 1.8.1 on a virtual environment with python38.

I know they are only warnings, but maybe there is something wrong with my setup.

Best,
Christian

In file included from /var/folders/nz/2ft4b1c53sd0cwgd_4rfpqqr0000gn/T/jitcxde_7sb1fvpt/jitced_4.c:2:
In file included from /Users/me/.pyenv/versions/3.8.9/include/python3.8/Python.h:25:
In file included from /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/stdio.h:64:
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/_stdio.h:93:16: warning: pointer is missing a nullability type specifier (_Nonnull, _Nullable, or _Null_unspecified) [-Wnullability-completeness]
        unsigned char   *_base;
                        ^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/_stdio.h:93:16: note: insert '_Nullable' if the pointer may be null
        unsigned char   *_base;
                        ^
                          _Nullable 
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/_stdio.h:93:16: note: insert '_Nonnull' if the pointer should never be null
        unsigned char   *_base;
                        ^
                          _Nonnull 
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/_stdio.h:138:32: warning: pointer is missing a nullability type specifier (_Nonnull, _Nullable, or _Null_unspecified) [-Wnullability-completeness]
        int     (* _Nullable _read) (void *, char *, int);
                                          ^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/_stdio.h:138:32: note: insert '_Nullable' if the pointer may be null
        int     (* _Nullable _read) (void *, char *, int);
                                          ^
                                           _Nullable
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/_stdio.h:138:32: note: insert '_Nonnull' if the pointer should never be null
        int     (* _Nullable _read) (void *, char *, int);
                                          ^
                                           _Nonnull
...
int     sigvec(int, struct sigvec *, struct sigvec *);
                                                   ^
                                                    _Nonnull
219 warnings generated.

When trying the first example in the docs:

from jitcdde import jitcdde, y, t
import numpy

τ = 15
n = 10
β = 0.25
γ = 0.1

f = [ β * y(0,t-τ) / (1 + y(0,t-τ)**n) - γ*y(0) ]
DDE = jitcdde(f)

DDE.add_past_point(-1.0, [1.0], [0.0])
DDE.add_past_point( 0.0, [1.0], [0.0])

DDE.step_on_discontinuities()

This is the output I see before the output array([ 1.18037923]) is returned:

Generating, compiling, and loading C code.
Using default integration parameters.

C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\jitcdde\_jitcdde.py:470: UserWarning: Traceback (most recent call last):
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\_msvccompiler.py", line 382, in compile
    self.spawn(args)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\_msvccompiler.py", line 501, in spawn
    return super().spawn(cmd)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\ccompiler.py", line 909, in spawn
    spawn(cmd, dry_run=self.dry_run)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\spawn.py", line 38, in spawn
    _spawn_nt(cmd, search_path, dry_run=dry_run)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\spawn.py", line 81, in _spawn_nt
    "command %r failed with exit status %d" % (cmd, rc))
distutils.errors.DistutilsExecError: command 'C:\\Program Files (x86)\\Microsoft Visual Studio 14.0\\VC\\BIN\\x86_amd64\\cl.exe' failed with exit status 2

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\core.py", line 148, in setup
    dist.run_commands()
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\dist.py", line 955, in run_commands
    self.run_command(cmd)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\dist.py", line 974, in run_command
    cmd_obj.run()
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\setuptools-27.2.0-py3.6.egg\setuptools\command\build_ext.py", line 77, in run
    _build_ext.run(self)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\Cython\Distutils\old_build_ext.py", line 185, in run
    _build_ext.build_ext.run(self)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\command\build_ext.py", line 339, in run
    self.build_extensions()
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\Cython\Distutils\old_build_ext.py", line 193, in build_extensions
    self.build_extension(ext)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\setuptools-27.2.0-py3.6.egg\setuptools\command\build_ext.py", line 198, in build_extension
    _build_ext.build_extension(self, ext)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\command\build_ext.py", line 533, in build_extension
    depends=ext.depends)
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\_msvccompiler.py", line 384, in compile
    raise CompileError(msg)
distutils.errors.CompileError: command 'C:\\Program Files (x86)\\Microsoft Visual Studio 14.0\\VC\\BIN\\x86_amd64\\cl.exe' failed with exit status 2

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\jitcdde\_jitcdde.py", line 468, in _initiate
    self.generate_f_C()
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\jitcdde\_jitcdde.py", line 458, in generate_f_C
    verbose = verbose
  File "C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\distutils\core.py", line 163, in setup
    raise SystemExit("error: " + str(msg))
SystemExit: error: command 'C:\\Program Files (x86)\\Microsoft Visual Studio 14.0\\VC\\BIN\\x86_amd64\\cl.exe' failed with exit status 2

  warn(format_exc())
C:\Users\Gommersr\AppData\Local\Continuum\Anaconda3\lib\site-packages\jitcdde\_jitcdde.py:471: UserWarning: Generating compiled integrator failed; resorting to lambdified functions.
  warn("Generating compiled integrator failed; resorting to lambdified functions.")

That's way too noisy; shouldn't be too hard to silence that - I would expect just a UserWarning.

Hi,

can you please explain where UserWarning: The target time is smaller than the current time. comes from and how we can mute it?

(as I am running a fit, the optimizer spit it out a huge amount of times, making the jupyter notebook very heavy...)

Thanks

ps: I am on mac os, with jitcdde.__version__ = 1.5.0

Hi,
I am working on a system of delay differential equations where the vector field is quite nonlinear. My main interest is calculating the largest Lyapunov Exponent as function of vector field parameters. For this purpose, I can use the algorithm itself because it automatically calculates the Jacobian matrix. Unfortunately, this Jacobian matrix is again nonlinear due to nonlinear character of the vector field and I obtain incorrect values considering the expected Lyapunov Exponents.
My question is: how could I set the Jacobian matrix before hand instead of being dependent on the algorithm?

Thank you,
Edmilson Roque.

After changing the code from the Kuramoto example to this:

rom jitcdde import jitcdde, y, t
from numpy import pi, arange, random, max
from symengine import sin
import numpy as np
import matplotlib.pyplot as plt
from time import gmtime, strftime
import os
import sympy
	
n=2
delay=1
np.random.seed(123)

def orderParam(phis):
		R=0
		for phi in phis:

				R+=np.exp(1j*phi)
		return np.abs(R)/len(phis)

def Simulate(c,initphi,forward=True):
	print ("hier komt de ordeparameter van de initiele toestand")
	print (orderParam(initphi))
	
	omega= random.standard_cauchy(n)
	
	q = 0.10
	A = random.choice( [1,0], size=(n,n), p=[q,1-q] )
	
		
	def kuramotos():
		for i in range(n):
			yield omega[i]+ omega[i]*sum(
						sin(y(j,t-delay)-y(i))
						for j in range(n)
						if A[j,i]
						)
		
	I = jitcdde(kuramotos,n=n,verbose=False)
	I.set_integration_parameters(rtol=0,atol=1e-5)
		
	I.constant_past(initphi, time=0.0 )
	I.integrate_blindly( delay , 0.1 )

	orderparams=[]
	A=[]
	for time in I.t + arange(0,300,0.2):
		#print(*I.integrate(time) % (2*pi))
		print (time)
		A=(I.integrate(time) % (2*pi))
		Rt=orderParam(A)
		orderparams.append(Rt)
	np.savetxt(simulationtitle+"/coupling: "+str(c)+"fw:"+str(forward)+" .txt",orderparams)
	
	print (c)
	
	sync=np.average(orderparams[-1000:])
	return A,sync

The code runs way way slower now, but I don't see why at all, it also keeps running that slow when I decrease the accuracy.

Hi,

First off, this is a really neat package that I think is exactly what I need to run DDE models for my research. Thanks!

I have a question about how to speed up the integration of my code (a minimal working example is below). Very briefly, I am using a variable-time delay model to predict the abundance of an insect species with juvenile (J) and adult (A) life stages. A third DDE equation describes juvenile survival (S) and the fourth equation describes the development time delay (tau). The maturation rate from juveniles to adults is dependent on the habitat temperature, which is modeled via a sine wave. For simplicity, I am using a constant past and I am integrating blindly (following advice in the documentation). To deal with initial discontinuities, I am simply running the model for long enough for transients to dissipate and the system to settle into equilibrium cycles.

The codes seems to be working (awesome!), but it takes 11-12 hours to run on my Mac laptop with Python 3.8.8 and JitCDDE 1.8. For reference, I ran a similar model years ago using the PyDDE package on the same laptop and the code would take ~15-20 minutes to run. I don't have access to a better computer, but I feel that I should be able to write a more efficient code that speeds up the integration. I am a Python newbie, so I am sure that any problems lie with my implementation of the code. I tried my best to follow the Common Mistakes and Questions suggestions, but must admit that some of the information about generators and Clang went over my head.

I hesitate to post this as an issue because it's very likely my own naive coding, but I am wondering if you have any suggests for speeding up the integration. Thanks so much in advance, it would be a life saver!

Best,
Chris

# Delay differential equation model for predicting insect population dynamics
# under seasonal temperature variation and climate change


# IMPORT PACKAGES
from numpy import array, arange, hstack, vstack, savetxt, pi
from jitcdde import jitcdde, y, t
from symengine import exp, sin
from matplotlib.pyplot import subplots, xlabel, ylabel, xlim, ylim #, yscale
from pandas import read_csv


# INPUT TEMPERATURE RESPONSE DATA
tempData = read_csv("Temperature response data.csv")



# DEFINE MODEL PARAMETERS
# Time parameters
yr = 365.0 # days in year
max_years = 10. # how long to run simulations
tstep = 1. # time step = 1 day

# Habitat temperature parameters
meanT = 300.2
amplT = 1.36 
shiftT = 35.39

# Life history and competitive traits
# fecundity
bTopt = 8.9
Toptb = 300.4
sb = 3.64
# maturation
mTR = 0.0334
TR = 298
AmJ = 7000
skew = 1
AL = -70694
TL = 289.8
AH = 146729
TH = 308.7
# mortality
dJTR = 0.013
AdJ = 23770
dATR = 0.0265
AdA = 9710
# competition
qTR = 1
Toptq = Toptb
sq = sb
Aq = AdA
Tmax = meanT + amplT
qTemp = 1
qTopt = qTR*exp(Aq*(1./TR - 1./Tmax))


# FUNCTIONS
# Seasonal temperature variation (K) over time
def T(x):
    return meanT + amplT * sin(2*pi*(x + shiftT)/yr)

# Life history functions
# fecundity
def b(x):
    return bTopt * exp(-(T(x)-Toptb)**2/2./sb**2)

# maturation rates
def mJ(x):
    return mTR * T(x)/TR * exp(AmJ * (1./TR - 1./T(x))) / (1. + skew * (exp(AL*(1./TL-1./T(x)))+exp(AH*(1./TH-1./T(x)))))

# mortality rates
def dJ(x):
    return dJTR * exp(AdJ * (1./TR - 1./T(x)))
def dA(x):
    return dATR * exp(AdA * (1./TR - 1./T(x)))

# density-dependence due to competition
def q(x):
    return (1-qTemp) * qTR + qTemp * qTopt * exp(-(T(x)-Toptq)**2/2./sq**2)



# DDE MODEL
# Define state variables
J,A,S,τ = [y(i) for i in range(4)]


# Model
f = {
    J: b(t)*A*exp(-q(t)*A) - b(t-τ)*y(1,t-τ)*exp(-q(t-τ)*y(1,t-τ))*mJ(t)/mJ(t-τ)*S - dJ(t)*J,
    
    A: b(t-τ)*y(1,t-τ)*exp(-q(t-τ)*y(1,t-τ))*mJ(t)/mJ(t-τ)*S - dA(t)*A,
    
    S: S*(mJ(t)/mJ(t-τ)*dJ(t-τ) - dJ(t)),
    
    τ: 1. -  mJ(t)/mJ(t-τ)
    }


# RUN DDE SOLVER
# Time and initial conditions
times = arange(0., max_years*yr, tstep)
init = array([ 10., 1., exp(-dJ(-1e-3)/mTR), 1./mTR ])


# DDE solver
DDE = jitcdde(f, max_delay=100, verbose=False)
DDE.constant_past(init)
DDE.integrate_blindly(0.01)
DDE.compile_C(simplify=False, do_cse=False, chunk_size=30, verbose=True) # These options are used in an attempt to speed up the compiler


# Save data array containing time and state variables
data = vstack([ hstack([time,DDE.integrate(time)]) for time in times ])
#filename = 'Time_series.csv'  
#savetxt(filename, data, fmt='%s', delimiter=",", header="Time,J,A,S,tau", comments='') 



# PLOT
fig,ax = subplots()
ax.plot(data[:,0], data[:,1], label='J')
ax.plot(data[:,0], data[:,2], label='A')
#ax.plot(data[:,0], data[:,3], label='S')
#ax.plot(data[:,0], data[:,4], label='τ')
ax.legend(loc='best')
xlabel("time (days)")
ylabel("population density")
xlim((max_years-10)*yr,max_years*yr)
ylim(0,40)
#yscale("log")
#ylim(0.1,100)

i want to use jitcdde package using conda install. what is the alternative to pip3 install.[ all my python packages are conda installed instead of pip3 install]

First of all, I have to admit that I don't have much experience with DDEs. That being said I get the concept behind them, and JiTCDDE seems like the best Python package out there for solving them. (I am using version 1.3.3.) I think my code is almost there, but I'm getting an error that I believe is related to step size. Here is a working (or not-quite-working) example. (My actual code uses 6 different equations, but this example with 4 is sufficient to raise the error.)

import numpy as np
from jitcdde import jitcdde, y, t

# Params
a1 = 40.0
a2 = 50.0
u  = 0.69
f1 = 22.0
f2 = 8.4
d  = 3.0
e1 = 3.5
e2 = 5.0
K  = 22.0
τ  = 0.093

func = {
        y(1),
        -a1 * (y(0)**2 - u) * y(1) - f1 * y(0) * (y(0) + d) * (y(0) + e1), 
        y(3), 
        -a2 * (y(2)**2 - u) * y(3) - f2 * y(2) * (y(2) + d) * (y(2) + e2) + K * (y(1, t-τ) - y(3)),
       }

y0 = [-0.1, 0.025, -0.1, 0.025]

I = jitcdde(func)
I.set_integration_parameters()
I.constant_past(y0, time=0.0)
I.step_on_discontinuities()

data = []
for time in np.linspace(0, 2):
    data.append( I.integrate(time).to_list() )

It seemingly works for early iterations. (In this example, the first 35 or so steps take less than 3 seconds collectively.) But then the speed of each iteration rapidly slow down so that a few iterations later, it takes 40 seconds for a single iteration. The iteration after that takes a very long time (a couple minutes) before raising the following error:

UnsuccessfulIntegration Traceback (most recent call last)
in ()
8 data = []
9 for time in np.linspace(0, 2)):
---> 10 data.append( I.integrate(time).to_list() )

~\AppData\Local\Continuum\anaconda3\lib\site-packages\jitcdde_jitcdde.py in integrate(self, target_time)
795 continue
796
--> 797 if self._adjust_step_size():
798 self.DDE.accept_step()
799

~\AppData\Local\Continuum\anaconda3\lib\site-packages\jitcdde_jitcdde.py in _adjust_step_size(self)
728 if p > self.decrease_threshold:
729 self.dt = max(self.safety_factorp**(-1/self.q), self.min_factor)
--> 730 self._control_for_min_step()
731 return False
732 else:

~\AppData\Local\Continuum\anaconda3\lib\site-packages\jitcdde_jitcdde.py in _control_for_min_step(self)
714 if self.atol==0:
715 message += "\n• You did not allow for an absolute error tolerance (atol) though your DDE calls for it. Even a very small absolute tolerance (1e-16) may sometimes help."
--> 716 raise UnsuccessfulIntegration(message)
717
718 def _increase_chance(self, new_dt):

UnsuccessfulIntegration:
Could not integrate with the given tolerance parameters:

atol: 1.000000e-10
rtol: 1.000000e-05
min_step: 1.000000e-10

The most likely reasons for this are:
• You did not sufficiently address initial discontinuities.
• The DDE is ill-posed or stiff.

I've tried playing the a_tol, r_tol, and min_step, but that doesn't seem to help. I'm not sure how else to address the initial discontinuities besides I.step_on_discontinuities(). Finally, I don't believe there is a problem with the DDE itself since it has been successfully solved with MatLab's dde23.

Though I realize this is almost certainly a user issue rather than a package issue, any help and understanding would be greatly appreciated!

Thanks for making this package, it is exactly what I need for my master's thesis on explosive synchronisation in the Kuramoto model.

I have the following problem:
when I use your example code on the Kuramoto model (https://github.com/neurophysik/jitcdde/blob/master/examples/kuramoto_network.py)
and only do the following changes: fix tau to a fixed number (say 5) and take omega from a distribution.

I get the following warning "differential equation does not include a delay term", although in the code there is still written a delay term (tau is a fixed non-zero number, this also happens when I keep tau a random distribution)

from jitcdde import jitcdde, y, t
from numpy import pi, arange, random, max
from symengine import sin
	
n = 100
ω = random.uniform(0,1,n)
c = 0
q = 0.05
A = random.choice( [1,0], size=(n,n), p=[q,1-q] )

	
def kuramotos():
	for i in range(n):
		yield ω [i]+ c/(n-1)*sum(
					sin(y(j,t-5)-y(i))
					for j in range(n)
					if A[j,i]
					)
	
I = jitcdde(kuramotos,n=n,verbose=False)
I.set_integration_parameters(rtol=0,atol=1e-5)
	
I.constant_past( random.uniform(0,2*pi,n), time=0.0 )
I.integrate_blindly( max(τ) , 0.1 )
	
for time in I.t + arange(0,400,0.2):
	print(*I.integrate(time) % (2*pi))