Python scipy.integrate.ode() Examples

def _run_solout_test(self, integrator):
        # Check correct usage of solout
        ts = []
        ys = []
        t0 = 0.0
        tend = 10.0
        y0 = [1.0, 2.0]
        def solout(t, y):
            ts.append(t)
            ys.append(y.copy())
        def rhs(t, y):
            return [y[0] + y[1], -y[1]**2]
        ig = ode(rhs).set_integrator(integrator)
        ig.set_solout(solout)
        ig.set_initial_value(y0, t0)
        ret = ig.integrate(tend)
        assert_array_equal(ys[0], y0)
        assert_array_equal(ys[-1], ret)
        assert_equal(ts[0], t0)
        assert_equal(ts[-1], tend) 

def _run_solout_test(self, integrator):
        # Check correct usage of solout
        ts = []
        ys = []
        t0 = 0.0
        tend = 10.0
        y0 = [1.0, 2.0]

        def solout(t, y):
            ts.append(t)
            ys.append(y.copy())

        def rhs(t, y):
            return [y[0] + y[1], -y[1]**2]

        ig = ode(rhs).set_integrator(integrator)
        ig.set_solout(solout)
        ig.set_initial_value(y0, t0)
        ret = ig.integrate(tend)
        assert_array_equal(ys[0], y0)
        assert_array_equal(ys[-1], ret)
        assert_equal(ts[0], t0)
        assert_equal(ts[-1], tend) 

def _do_problem(self, problem, integrator, method='adams'):

        # ode has callback arguments in different order than odeint
        f = lambda t, z: problem.f(z, t)
        jac = None
        if hasattr(problem, 'jac'):
            jac = lambda t, z: problem.jac(z, t)

        ig = ode(f, jac)
        ig.set_integrator(integrator,
                          atol=problem.atol/10,
                          rtol=problem.rtol/10,
                          method=method)
        ig.set_initial_value(problem.z0, t=0.0)
        z = ig.integrate(problem.stop_t)

        assert_(ig.successful(), (problem, method))
        assert_(problem.verify(array([z]), problem.stop_t), (problem, method)) 

def __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0,
                 max_step=np.inf, rtol=1e-3, atol=1e-6, jac=None, lband=None,
                 uband=None, vectorized=False, **extraneous):
        warn_extraneous(extraneous)
        super(LSODA, self).__init__(fun, t0, y0, t_bound, vectorized)

        if first_step is None:
            first_step = 0  # LSODA value for automatic selection.
        elif first_step <= 0:
            raise ValueError("`first_step` must be positive or None.")

        if max_step == np.inf:
            max_step = 0  # LSODA value for infinity.
        elif max_step <= 0:
            raise ValueError("`max_step` must be positive.")

        if min_step < 0:
            raise ValueError("`min_step` must be nonnegative.")

        rtol, atol = validate_tol(rtol, atol, self.n)

        if jac is None:  # No lambda as PEP8 insists.
            def jac():
                return None

        solver = ode(self.fun, jac)
        solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,
                              min_step=min_step, first_step=first_step,
                              lband=lband, uband=uband)
        solver.set_initial_value(y0, t0)

        # Inject t_bound into rwork array as needed for itask=5.
        solver._integrator.rwork[0] = self.t_bound
        solver._integrator.call_args[4] = solver._integrator.rwork

        self._lsoda_solver = solver 

def _do_problem(self, problem, integrator, method='adams'):

        # ode has callback arguments in different order than odeint
        f = lambda t, z: problem.f(z, t)
        jac = None
        if hasattr(problem, 'jac'):
            jac = lambda t, z: problem.jac(z, t)
        ig = complex_ode(f, jac)
        ig.set_integrator(integrator,
                          atol=problem.atol/10,
                          rtol=problem.rtol/10,
                          method=method)
        ig.set_initial_value(problem.z0, t=0.0)
        z = ig.integrate(problem.stop_t)
        z2 = ig.y

        assert_array_equal(z, z2)
        assert_(ig.successful(), (problem, method))
        assert_(problem.verify(array([z]), problem.stop_t), (problem, method)) 

def reset(self):
        '''
        Reset the patient state to default intial state
        '''
        self._last_Qsto = self.init_state[0] + self.init_state[1]
        self._last_foodtaken = 0
        self.name = self._params.Name

        self._odesolver = ode(self.model).set_integrator('dopri5')
        self._odesolver.set_initial_value(self.init_state, self.t0)

        self._last_action = Action(CHO=0, insulin=0)
        self.is_eating = False
        self.planned_meal = 0 

def _my_odeint_(dfunc, V0, times, args=()):
    r'''For some of the systems odeint will switch to the BDF solver.
    In large enough systems, it then gets stuck trying to estimate the 
    Jacobian.  This will even cause segmentation faults.

    This routine has identical inputs to integrate.odeint, but relies on 
    integrate.ode.  It avoids BDF.

    In particular, this seems to be important for SIS heterogeneous 
    pairwise where the number of equations is very large.  I have found 
    that near equilibrium, this often is interpreted as being a stiff 
    system and it switches to bdf, which requires calculating a 
    Jacobian.  In some systems this is impractically large.
    
    See this question: 
        http://stackoverflow.com/q/40317096/2966723,
    with the answer by 
        Phillip: http://stackoverflow.com/users/1881610/phillip
        
    INPUT and OUTPUT are as integrate.odeint
    
    :Arguments: 
    Identical to integrate.odeint
        
    :Returns: 
    Identical to integrate.odeint
    '''

    r = integrate.ode(lambda t, X: dfunc(X, t, *args))
    r.set_integrator('vode', method='adams')
    r.set_initial_value(V0,times[0])
    V=[V0]
    for time in times[1:]:
        V.append(r.integrate(time))
    V = np.array(V)
    return V 

def set_state(self, state):
        if self.state is None or\
           np.linalg.norm(np.array(state)-np.array(self.state)) > 1e-12:
            # float64 required for the ode integrator
            self.state = np.array(state, dtype=np.float64).flatten()
        # set solver internal state
        self.solver = self.solver.set_initial_value(self.state)
        # get time from solver
        self.t = self.solver.t 

def __init__(self, name='ODEPlant', integrator='dopri5',
                 atol=1e-12, rtol=1e-12,
                 *args, **kwargs):
        super(ODEPlant, self).__init__(name=name, *args, **kwargs)
        integrator = kwargs.get('integrator', 'dopri5')
        atol = kwargs.get('atol', 1e-12)
        rtol = kwargs.get('rtol', 1e-12)

        # initialize ode solver
        self.solver = ode(self.dynamics).set_integrator(integrator,
                                                        atol=atol,
                                                        rtol=rtol) 

def _dynamics(self, t, state, action):
        """ Evaluate the system dynamics

        Parameters
        ----------
        t: float
            Input Parameter required for the odesolver for time-dependent
            odes. Has no influence in this system.
        state: 2x1 array[float]
            The current state of the system
        action: 1x1 array[float]
            The action to be applied at the current time step

        Returns
        -------
        dz: 2x1 array[float]
            The ode evaluated at the given inputs.
        """

        inertia = self.m * self.l ** 2
        dz = np.zeros((2, 1))
        dz[0] = self.g / self.l * np.sin(
            state[1]) + action / inertia - self.b / inertia * state[0]
        dz[1] = state[0]

        return dz 

def set_integrator(self, name, **kwargs):
		self.f = ode(self.model.integration_loop)
		self.f.set_integrator(name, **kwargs)
		self.f.set_initial_value(self.init_state, self.t0) 

def _get_solver(self, f, jac):
        solver = ode(f, jac)
        if self.solver_uses_jac:
            solver.set_integrator(self.solver_name, atol=1e-9, rtol=1e-7,
                                  with_jacobian=self.solver_uses_jac)
        else:
            # XXX Shouldn't set_integrator *always* accept the keyword arg
            # 'with_jacobian', and perhaps raise an exception if it is set
            # to True if the solver can't actually use it?
            solver.set_integrator(self.solver_name, atol=1e-9, rtol=1e-7)
        return solver 

def _run_solout_break_test(self, integrator):
        # Check correct usage of stopping via solout
        ts = []
        ys = []
        t0 = 0.0
        tend = 10.0
        y0 = [1.0, 2.0]

        def solout(t, y):
            ts.append(t)
            ys.append(y.copy())
            if t > tend/2.0:
                return -1

        def rhs(t, y):
            return [y[0] + y[1], -y[1]**2]

        ig = ode(rhs).set_integrator(integrator)
        ig.set_solout(solout)
        ig.set_initial_value(y0, t0)
        ret = ig.integrate(tend)
        assert_array_equal(ys[0], y0)
        assert_array_equal(ys[-1], ret)
        assert_equal(ts[0], t0)
        assert_(ts[-1] > tend/2.0)
        assert_(ts[-1] < tend) 

def test_concurrent_ok(self):
        f = lambda t, y: 1.0

        for k in xrange(3):
            for sol in ('vode', 'zvode', 'lsoda', 'dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.1)
                assert_allclose(r2.y, 0.2)

            for sol in ('dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.3)
                assert_allclose(r2.y, 0.2) 

def test_concurrent_fail(self):
        for sol in ('vode', 'zvode', 'lsoda'):
            f = lambda t, y: 1.0

            r = ode(f).set_integrator(sol)
            r.set_initial_value(0, 0)

            r2 = ode(f).set_integrator(sol)
            r2.set_initial_value(0, 0)

            r.integrate(r.t + 0.1)
            r2.integrate(r2.t + 0.1)

            assert_raises(RuntimeError, r.integrate, r.t + 0.1) 

def _do_problem(self, problem, integrator, method='adams'):

        # ode has callback arguments in different order than odeint
        f = lambda t, z: problem.f(z, t)
        jac = None
        if hasattr(problem, 'jac'):
            jac = lambda t, z: problem.jac(z, t)

        integrator_params = {}
        if problem.lband is not None or problem.uband is not None:
            integrator_params['uband'] = problem.uband
            integrator_params['lband'] = problem.lband

        ig = self.ode_class(f, jac)
        ig.set_integrator(integrator,
                          atol=problem.atol/10,
                          rtol=problem.rtol/10,
                          method=method,
                          **integrator_params)

        ig.set_initial_value(problem.z0, t=0.0)
        z = ig.integrate(problem.stop_t)

        assert_array_equal(z, ig.y)
        assert_(ig.successful(), (problem, method))
        assert_(ig.get_return_code() > 0, (problem, method))
        assert_(problem.verify(array([z]), problem.stop_t), (problem, method)) 

def set_system_equation(self, system_equation, jac=None):
        # Docstring of superclass
        super().set_system_equation(system_equation, jac)
        self._ode = ode(system_equation, jac).set_integrator(self._integrator, **self._solver_args) 

def __init__(self, integrator='dopri5', **kwargs):
        """
        Args:
            integrator(str): String to choose the integrator from the scipy.integrate.ode
            kwargs(dict): All parameters that can be set in the "set_integrator"-method of scipy.integrate.ode
        """
        self._solver = None
        self._solver_args = kwargs
        self._integrator = integrator 

def set_rhs(self, rhs=None, reset=True):
        """This set_rhs function fully instantiates the scipy ode object behind the scenes."""

        if rhs is None:
            rhs = {'rhs': None}

        if callable(rhs):
            rhs = {'rhs': rhs}

        if 'rhs' not in rhs:
            raise Exception('ODE_Method requires at minimum a specification of an rhs function.')

        self.rhs = self._state_type_converter.transform_rhs_funcs(rhs)

        self._ODE = ode(self.rhs['rhs'])

        self._ODE._integrator = qiskit_zvode(method=self.options.method,
                                             order=self.options.order,
                                             atol=self.options.atol,
                                             rtol=self.options.rtol,
                                             nsteps=self.options.nsteps,
                                             first_step=self.options.first_step,
                                             min_step=self.options.min_step,
                                             max_step=self.options.max_step
                                             )

        # Forces complex ODE solving
        if not self._ODE._y:
            self._ODE.t = 0.0
            self._ODE._y = np.array([0.0], complex)
        self._ODE._integrator.reset(len(self._ODE._y), self._ODE.jac is not None)

        self._ODE.set_initial_value(self._y, self._t)

        self._reset_method(reset) 

def __init__(self, t0=None, y0=None, rhs=None, options=None):

        # all de specification arguments are necessary to instantiate scipy ode object
        if (t0 is None) or (y0 is None) or (rhs is None):
            raise Exception('QiskitZVODE solver requires t0, y0, and rhs at instantiation.')

        # initialize internal attribute for storing scipy ode object
        self._ODE = None

        super().__init__(t0, y0, rhs, options) 

def _get_solver(self, f, jac):
        solver = ode(f, jac)
        if self.solver_uses_jac:
            solver.set_integrator(self.solver_name, atol=1e-9, rtol=1e-7,
                                  with_jacobian=self.solver_uses_jac)
        else:
            # XXX Shouldn't set_integrator *always* accept the keyword arg
            # 'with_jacobian', and perhaps raise an exception if it is set
            # to True if the solver can't actually use it?
            solver.set_integrator(self.solver_name, atol=1e-9, rtol=1e-7)
        return solver 

def test_concurrent_ok(self):
        f = lambda t, y: 1.0

        for k in xrange(3):
            for sol in ('vode', 'zvode', 'lsoda', 'dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.1)
                assert_allclose(r2.y, 0.2)

            for sol in ('dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.3)
                assert_allclose(r2.y, 0.2) 

def test_concurrent_fail(self):
        for sol in ('vode', 'zvode', 'lsoda'):
            f = lambda t, y: 1.0

            r = ode(f).set_integrator(sol)
            r.set_initial_value(0, 0)

            r2 = ode(f).set_integrator(sol)
            r2.set_initial_value(0, 0)

            r.integrate(r.t + 0.1)
            r2.integrate(r2.t + 0.1)

            assert_raises(RuntimeError, r.integrate, r.t + 0.1) 

def __init__(self, fun, t0, y0, t_bound, first_step=None, min_step=0.0,
                 max_step=np.inf, rtol=1e-3, atol=1e-6, jac=None, lband=None,
                 uband=None, vectorized=False, **extraneous):
        warn_extraneous(extraneous)
        super(LSODA, self).__init__(fun, t0, y0, t_bound, vectorized)

        if first_step is None:
            first_step = 0  # LSODA value for automatic selection.
        elif first_step <= 0:
            raise ValueError("`first_step` must be positive or None.")

        if max_step == np.inf:
            max_step = 0  # LSODA value for infinity.
        elif max_step <= 0:
            raise ValueError("`max_step` must be positive.")

        if min_step < 0:
            raise ValueError("`min_step` must be nonnegative.")

        rtol, atol = validate_tol(rtol, atol, self.n)

        if jac is None:  # No lambda as PEP8 insists.
            def jac():
                return None

        solver = ode(self.fun, jac)
        solver.set_integrator('lsoda', rtol=rtol, atol=atol, max_step=max_step,
                              min_step=min_step, first_step=first_step,
                              lband=lband, uband=uband)
        solver.set_initial_value(y0, t0)

        # Inject t_bound into rwork array as needed for itask=5.
        solver._integrator.rwork[0] = self.t_bound
        solver._integrator.call_args[4] = solver._integrator.rwork

        self._lsoda_solver = solver 

def line(self, x0):
        """Returns the field line passing through x0.
        Refs: http://folk.uib.no/fcihh/seminar/lec1.pdf and lect2.pdf
              http://numbercrunch.de/blog/2013/05/visualizing-streamlines/
        and especially: "Electric field lines don't work",
        http://scitation.aip.org/content/aapt/journal/ajp/64/6/10.1119/1.18237
        """

        if None in [XMIN, XMAX, YMIN, YMAX]:
            raise ValueError('Domain must be set using init().')

        # Set up integrator for the field line
        streamline = lambda t, y: list(self.direction(y))
        solver = ode(streamline).set_integrator('vode')

        # Initialize the coordinate lists
        x = [x0]

        # Integrate in both the forward and backward directions
        dt = 0.008

        # Solve in both the forward and reverse directions
        for sign in [1, -1]:

            # Set the starting coordinates and time
            solver.set_initial_value(x0, 0)

            # Integrate field line over successive time steps
            while solver.successful():

                # Find the next step
                solver.integrate(solver.t + sign*dt)

                # Save the coordinates
                if sign > 0:
                    x.append(solver.y)
                else:
                    x.insert(0, solver.y)

                # Check if line connects to a charge
                flag = False
                for c in self.charges:
                    if c.is_close(solver.y):
                        flag = True
                        break

                # Terminate line at charge or if it leaves the area of interest
                if flag or not (XMIN < solver.y[0] < XMAX) or \
                  not YMIN < solver.y[1] < YMAX:
                    break

        return FieldLine(x) 

def __init__(self, name="InvertedPendulum", l=.5, m=.15, g=9.82, b=0.,
                 dt=.05, init_m=0., init_std=.01,
                 plant_noise=np.array([0.01, 0.01]) ** 2,
                 u_min=np.array([-1.]), u_max=np.array([1.]),
                 target=np.array([0.0, 0.0]),
                 verbosity=1, norm_x=None, norm_u=None):
        """
        Parameters
        ----------
        name: str, optional
            The name of the system
        l: float, optional
            The length of the pendulum
        m: float, optional
            The mass of the pendulum
        g: float, optional
            The gravitation constant
        b: float, optional
            The friction coefficient of the system
        init_m: 2x0 1darray[float], optional
            The initial state mean
        init_std: float, optional
            The standard deviation of the start state sample distribution.
            Note: This is not(!) the uncertainty of the state but merely allows
            for variation in the initial (deterministic) state.
        u_min: 1x0 1darray[float], optional
            The maximum torque applied to the system
        u_max: 1x0 1darray[float], optional
            The maximum torquie applied to the system
        target: 2x0 1darray[float], optional
            The target state
        """
        super(InvertedPendulum, self).__init__(name, 2, 1, dt, init_m, init_std,
                                               plant_noise, u_min, u_max, target,
                                               verbosity)
        self.odesolver = ode(self._dynamics)
        self.l = l
        self.m = m
        self.g = g
        self.b = b
        self.p_origin = np.array([0.0, 0.0])
        self.l_mu = np.array([0.05, .02])  # TODO: This should be somewhere else
        self.l_sigm = np.array([0.05, .02])
        self.target = target
        self.target_ilqr = init_m

        warnings.warn("Normalization turned off for now. Need to look into it")
        max_deg = 30
        if norm_x is None:
            norm_x = np.array([1., 1.])
            # norm_x = np.array([np.sqrt(g/l), np.deg2rad(max_deg)])

        if norm_u is None:
            norm_u = np.array([1.])
            # norm_u = np.array([g*m*l*np.sin(np.deg2rad(max_deg))])

        self.norm = [norm_x, norm_u]
        self.inv_norm = [arr ** -1 for arr in self.norm]

        self._init_safety_constraints() 

def __init__(self, name='CartPole', dt=0.1, l=0.5, m=0.5, M=0.5, b=0.1, g=9.82,
                 init_m=np.array([0.0, 0.0, 0.0, 0.0]), visualize=True, init_std=0.0,
                 u_min=np.array([-4.0]), u_max=np.array([4.0]),
                 norm_x=None, norm_u=None, plant_noise=np.array([0.02, 0.05, 0.02, 0.05]) ** 2, verbosity=1):
        super(CartPole, self).__init__(name, 4, 1, dt, init_m, init_std,
                                       plant_noise,
                                       u_min, u_max,
                                       np.array([0.0, l, 0.0]), verbosity)

        self.ns_ode = 4

        self.current_state = None

        self.T = 15
        # initialize the physical properties
        self.l = l
        self.m = m
        self.M = M
        self.b = b
        self.g = g
        self.visualize = visualize

        self.l_mu = np.array(
            [.05, .05, .05, .05])  # TODO: This should be somewhere else
        self.l_sigm = np.array([.05, .05, .05, .05])

        self.idx_angles = np.array([2])
        self.obs_angles_sin = np.array([3])
        self.obs_angles_cos = np.array([4])

        self.target_ilqr = init_m

        self.delay = 20.0  # fps

        self.D_cost = np.array([40, 20, 40])
        self.R_cost = np.array([1.0])

        if norm_x is None:
            x_max = 3.0
            x_dot_max = 5.0
            max_deg = 30
            norm_x = np.array([x_max, x_dot_max, np.deg2rad(max_deg), 2 * np.sqrt(g / l)])

        if norm_u is None:
            norm_u = self.u_max - self.u_min

        self.norm = [norm_x, norm_u]
        self.inv_norm = [arr ** -1 for arr in self.norm]

        self._init_safety_constraints()

        self.odesolver = ode(self._dynamics)
        self.name = name 

def imageProcessing(self, inputLocation, outputLocation,
                        tempA, tempB, initialCondition, Ib, t):
        """Process the image with the input arguments.

        Args:
            inputLocation (str): The string path for the input image.
            outputLocation (str): The string path for the output processed
                image.
            tempA (:obj:`list` of :obj:`list`of :obj:`float`): Feedback
                template.
            tempB (:obj:`list` of :obj:`list`of :obj:`float`): Control
                template.
            initialCondition (float): The initial state.
            Ib (float): System bias.
            t (numpy.ndarray): A numpy array with evenly spaced numbers
                representing time points.
        """
        gray = img.open(inputLocation).convert('RGB')
        self.m, self.n = gray.size
        u = np.array(gray)
        u = u[:, :, 0]
        z0 = u * initialCondition
        Bu = sig.convolve2d(u, tempB, 'same')
        z0 = z0.flatten()
        tFinal = t.max()
        tInitial = t.min()
        if t.size > 1:
            dt = t[1] - t[0]
        else:
            dt = t[0]
        ode = sint.ode(self.f) \
            .set_integrator('vode') \
            .set_initial_value(z0, tInitial) \
            .set_f_params(Ib, Bu, tempA)
        while ode.successful() and ode.t < tFinal + 0.1:
            ode_result = ode.integrate(ode.t + dt)
        z = self.cnn(ode_result)
        out_l = z[:].reshape((self.n, self.m))
        out_l = out_l / (255.0)
        out_l = np.uint8(np.round(out_l * 255))
        # The direct vectorization was causing problems on Raspberry Pi.
        # In case anyone face a similar issue, use the below
        # loops rather than the above direct vectorization.
        # for i in range(out_l.shape[0]):
        #     for j in range(out_l.shape[1]):
        #         out_l[i][j] = np.uint8(round(out_l[i][j] * 255))
        out_l = img.fromarray(out_l).convert('RGB')
        out_l.save(outputLocation)

    # general image processing for given templates