neurophysik / jitcdde

Just-in-time compilation for delay differential equations
Other
57 stars 14 forks source link

Integration error when doing optimization #5

Closed hangyao closed 7 years ago

hangyao commented 7 years ago

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!

Wrzlprmft commented 7 years ago

At a first glance, you are using the Mackey–Glass equation with very atypical parameters and initial conditions that do not correspond to them. I am not surprised that the integrator cannot handle this; probably this is not a well-posed problem.

Is this really what you want?

Wrzlprmft commented 7 years ago

That being said, if your optimizer tests values that lead to ill-posed problems and you want to handle this, the best way is probably to use Python’s exception handling and catch UnsuccessfulIntegration (which you can import from jitcdde).

However, this only works if not all problems are ill-posed. In your case, it will just further obfuscate the actual problem.

hangyao commented 7 years ago

Thank you @Wrzlprmft ! Yes you are right. My equation and parameters were ill-defined. Now it is working properly.