problems about time-dependent delay and jitcdde parameter f_sym

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

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, re2+im2]) ) 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

First, some clarification: f_sym is not called; it is translated to C code, which is then called. It is a list of (or generator function yielding) SymEngine expressions. You cannot replace these expressions with functions or something similar as SymEngine won’t know what to do with them.

I cannot run your example code because several things are undefined. If you want to define something piecewise, you have to build it with a sharp sigmoidal for two reasons:
- SymEngine doesn’t support an analogue to SymPy’s Piecewise.
- Integrators do not like discontinuities.
Creating multiple jitcdde objects hardly ever makes sense. To access the current time, just use t (imported with from jitcdde import t) like you do in a normal delay term.
Integrating over the past is not natively supported (see Issue #2). What you can try is to implement a reasonable discretisation of your integral into f_sym.
There is an example of a state-dependent delay here. Time-dependent delays are not that different. That being said, if you can provide a good example (in form of equations), ideally with some expectations about the outcome, I might implement it as an example.

Thank you for your reply!

I know f_sym can't be called and it needs a SymEngine expressions. Maybe I should use an example to explain.

If I want to implement following piecewise function (just an example, not my real model),

$f(t,\,y(t-\tau_{1}),...)=\begin{cases} y(t) & t<t_{1}\\ y(t-t_{1}) & t_{1}<t<t_{2}\\ (t-t_{2})y(t-t_{1})+y(t-t_{2}) & t>t_{2} \end{cases}$

As you see, I want the expression to be changed with current time t, how can I implement it? I tried following code, but it doesn't work.

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

t1=5.0
t2=10.0

def f():    # I just use this function to define the delay.
    if t<t1:        # Use t or DDE.t?
        return y(t)
    elif t>t1 and t<t2:
        return y(t-t1)
    else:
        return (t-t2)*y(t-t1)+y(t-t2)

fdelay = [f()]    
DDE = jitcdde(fdelay)
DDE.constant_past([1.0])
DDE.step_on_discontinuities()

data = []
tlist = []
for time in np.arange(DDE.t, DDE.t+1000, 10):
    data.append( DDE.integrate(time)[0] )
    tlist.append(time)

As you said in your last reply, I use t to get current time. However, I don't understand how can translated C codes know the changed expression of delay term if I don't update list fdelay.

1.I just put a part of my previous code so it indeed lost some definition. Just want to illustrate my question.

2.As you see above, I doubted that whether jitcdde could get the changed expression from a single list fdelay, so I try to create multiple jitcdde objects in my previous codes.

3.This Issue #2 is similar to my real model ( $f=\int_{t-\delta}^{t}T(t,\tau)y(t,\tau)d\tau$ ), and I have done discretisation on my integral using scipy.integrate.simps. Although it seems awkward, it does give right symbolic expression. Thank you for your explanation.

4.I know this state-dependent example, but it just add a y as a augment of delay function. Sometimes we need to change the whole function form of delay, as I show in above example. Besides, I think an integral delay like my real model and Issue #2 could be common, maybe they are good examples.

5.How can I add some extra parameters? I may have to sweep those parameters to get multiple solutions.

Thank you again! Wance

The problem with your code is that the conditional is evaluated in the line fdelay = [f()], i.e., before you even call jitcdde. The best way to implement your condition are probably sharp sigmoidals. Also note that you have to give an argument to indicate which component of y you want to access – even if it is one-dimensional.


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

eps = Rational(1,100)
sigmoid = lambda x: (tanh(x/eps)+1)/2

t1 = 5.0
t2 = 10.0

f = [
          sigmoid(-t+t1)*y(0,t)
        + sigmoid( t-t1)*sigmoid(-t+t2)*y(0,t-t1)
        + sigmoid( t-t2)*( (t-t2)*y(0,t-t1)+y(0,t-t2) )
    ]

DDE = jitcdde(f,verbose=False)
DDE.constant_past([1.0])
DDE.step_on_discontinuities()

times = np.arange(DDE.t,DDE.t+20,1)
data = [ DDE.integrate(time)[0] for time in times ]

An alternative to the above is to integrate up to t1, use get_state to obtain the past, initialise a new integrator with it, integrate up to t2 and so on. However, after every re-intialisation, you would have to address initial discontinuities again (unless the cases are joint in a sufficiently smooth fashion). This approach may be computationally more feasible if you dwell in one case long enough, as you can avoid computations of complicated and irrelevant zeros (or almost zeros) like sigmoid(-100).
simps seems to handle symbolic input well, if you operate it correctly (e.g., simps( [ y(0,t-τ) for τ in range(-10,0) ])). The only drawback is that it uses floats for factors, which prevents certain optimisations (which should not matter much though). Anyway, right now, the best I can do is to supply a utility function that is similar to simps. The only better alternative is to implement a semi-symbolic integration on a low level, which is a lot of work.
First note that neither of your examples so far requires time-dependent delays. Still, if you want such a thing, just write a respective term like you would on paper, e.g., y(0,t-2-sin(t)) would give you a time-dependent delay.
Use the argument control_pars of jitcdde (look at its documentation for details).

I forgot to reply since I am busy recently. Thank you for your solution!

Your "sharp sigmoidals" solution is quite clever! But it's impractical for my case: the delay in my model is (I didn't show all details in the last comment.) codecogseqn When t>delta, the delay is easy to handle. The question is when 0<t<delta, it requires a dynamic discretization of the integral, which can't be described by sharp sigmoidals. And I also want to make the discretization more precise by adding more points at first, because the integral region is small when t is small. I did think up a solution: I added a list of anchors whose value is 0 when t<0, and one anchor valued 1 at t=0.

anchors = []
for tim in np.linspace(-delta+tstart,tstart,n+1):
    anchors.append((tim,[0.0,0.0],[0.0,0.0]))
anchors[-1] = (tstart,[re0,im0],[func.we(0.0)*im0,-func.we(0.0)*re0])

Then I only use this delay form, $\int_{t-\delta}^{t}T(t,\tau)y(t,\tau)d\tau$

But in this way I introduce a discontinuity at t=0, the solution will lose a part after I step_on_discontinuty. What's worse, the output result doesn't show what the approximated theory expected. I doubt if the discontinuity cause some mistakes.

So I still hope to find a dynamic discretation which can be accepted by jitcdde. I will try get_state() when I have time. Maybe it will work. Or I will try to use other packages to solve 0<t<delta first, then input it into jitcdde.

Thank you!

Before I answer this:

This equation contains τ outside the integral as well as a the integration variable. This is at the very least bad notation, but it also suggests that your actual conditions are something like 0<t<δ and t>δ. This would make the entire problem much easier to handle.

Oh, you're right. The "tau" outside the integral should be "delta". I made mistakes when I was writing the comment.

Okay, here are some ways to do it:

Preamble

from symengine import cos, tanh, Symbol, sympify, Rational
from jitcdde import t, y, jitcdde, past_y, anchors
import numpy as np

n = 1

def trapezoid(integrand,variable,lower,upper,steps=21):
    step_size = (sympify(upper-lower)/sympify(steps-1)).simplify()

    def evaluations():
        yield integrand.subs(variable,lower)/2
        for i in range(1,steps-1):
            position = lower + i*step_size
            yield integrand.subs(variable,position)
        yield integrand.subs(variable,upper)/2

    result = sum(x for x in evaluations())*step_size
    result = result.subs({
            past_y(t,i,anchors(t)): y(i)
            for i in range(n)
        })
    return result

eps = Rational(1,100)
sigmoid = lambda x: (tanh(x/eps)+1)/2

τ = Symbol("τ")
T = cos(τ)
δ = 2

The main part of this is a symbolic implementation of the trapezoidal method. The complex substitution at the end is to avoid unnecessary overlap. (I will implement a similar function soon.)

The rest are imports, definitions, and an example problem.

Variant 1: One Integral

This is based on the insight that the first case of your integral definition is not necessary if T(t,τ) = 0 for τ<0. Thus we only need apply the integral once and can avoid time-dependent delays altogether. We need a sigmoid to ensure T(t,τ) = 0 for τ<0 though:

f = [ trapezoid( sigmoid(τ)*T*y(0,τ), τ, t-δ, t ) ]
DDE = jitcdde(f)
DDE.compile_C(simplify=False)
DDE.constant_past([1.0])

for time in np.arange(DDE.t,20,0.05):
    print(time,DDE.integrate(time)[0])

Note that we do not need to address initial discontinuities as there are none per construction: f(t=0)=0 since the integral is over a vanishing interval. The line DDE.compile_C(simplify=False) exists to avoid a lengthy and pointless simplification attempt of the integral.

This variant means that the integral is less precise at the beginning. This in turn only is a problem if the integral is large with respect to rest of the derivative and the initial value. The former is the case here; the latter isn’t.

Variant 2: Swapping Integrals with a Sigmoid

Here, we use an adaptive integral for the beginning and swap it with a fixed one (like in the previous variant) via sigmoidals later.

f = [
          sigmoid(-t+δ)*trapezoid(T*y(0,τ),τ, 0 ,t)
        + sigmoid( t-δ)*trapezoid(T*y(0,τ),τ,t-δ,t)
    ]
DDE = jitcdde(f,max_delay=1e8)
DDE.compile_C(simplify=False)
DDE.constant_past([1.0])

for time in np.arange(DDE.t,20,0.05):
    print(time,DDE.integrate(time)[0])

The problem with this is that we have to calculate the first integral all the time, even if it contributes (almost) nothing. Since the first integral always starts at 0, this also means that we have to drag a lot of the past data around.

Variant 3: Reinitialising the Integrator

This is like the above, just that we do not use sigmoids to switch, but instead reinitialise the integrator with a different f:

f_1 = [trapezoid(T*y(0,τ),τ, 0 ,t)]
f_2 = [trapezoid(T*y(0,τ),τ,t-δ,t)]

DDE_1 = jitcdde(f_1,max_delay=δ,verbose=False)
DDE_1.compile_C(simplify=False)
DDE_1.constant_past([1.0])

for time in np.linspace(0.05,δ,21):
    print(time,DDE_1.integrate_blindly(time)[0])

DDE_2 = jitcdde(f_2,verbose=False)
DDE_2.compile_C(simplify=False)
DDE_2.add_past_points(DDE_1.get_state())

for time in np.arange(DDE_2.t,20,0.05):
    print(time,DDE_2.integrate(time)[0])

This is arguably the cleanest variant and comes with the fewest problems. Note that we need to integrate blindly in the first part to make sure we finish the last integration step precisely at δ.

Note that all three variants produce comparable results for this example.

I now implemented a utility function for integration called quadrature – which is similar to trapezoid from the above examples. You can now do something like:


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

δ = 2

τ = Symbol("τ")
f_1 = [quadrature(cos(τ)*y(0,τ),τ, 0 ,t)]
f_2 = [quadrature(cos(τ)*y(0,τ),τ,t-δ,t)]

DDE_1 = jitcdde(f_1,max_delay=δ,verbose=False)
DDE_1.compile_C(simplify=False)
DDE_1.constant_past([1.0])

for time in np.linspace(0.05,δ,21):
    print(time,DDE_1.integrate_blindly(time)[0])

DDE_2 = jitcdde(f_2,verbose=False)
DDE_2.compile_C(simplify=False)
DDE_2.add_past_points(DDE_1.get_state())

for time in np.arange(DDE_2.t,20,0.05):
    print(time,DDE_2.integrate(time)[0])

I would be interested in any feedback you have on this.

Hi,

I am sorry to not reply you for a long time, and really thanks for your solutions! It must cost you a lot of time.

In fact, I have had tried your Variant 1 last month, and meanwhile, I also implemented my own solution (First compute the region (0, δ) using scipy.integrate.solve_ivp, then import the data as anchors into jitcdde. In this way I don't have to change the expression of integrand.) Except for some regions that jitcdde step on discontinuity, the results are perfectly matched (although my method is quite slow.) Your Variant 1: c 0 001 n 200 0 5_sr_60 100s_cl506 My solution:

I haven't tried your quadrature, now I am still using a modified variant 1. If I have time (probably after this semester), I will try it and give you feedback. Thank you again!

Best regards!

neurophysik / jitcdde