Python scipy.integrate() Examples

def _integrateOneStep(r, t, func, jac, args=(), full_output=False):
    '''
    Perform integration with just one step
    '''
    r.integrate(t)
    # wish to do checking at each iteration?  A feature to be
    # considered in the future
    if r.successful():
        if full_output:
            e = np.linalg.eig(jac(r.t, r.y, *args))[0]
            return r.y, r.successful(), e, max(e), min(e)
        else:
            return r.y
    else:
        try:
            np.linalg.eig(jac(r.t, r.y, *args))
        except:
            raise IntegrationError("Failed integration with x =  " + str(r.y))
        else:
            a = np.linalg.eig(jac(r.t, r.y, *args))[0]
            raise IntegrationError("Failed integration with x =  " + str(r.y) +
                                   " and max/min eigenvalues of the jacobian is "
                                   + str(max(a)) + " and " + str(min(a))) 

def integrate(ode, x0, t, full_output=False):
    '''
    A wrapper on top of :mod:`odeint <scipy.integrate.odeint>` using
    :class:`DeterministicOde <pygom.model.DeterministicOde>`.

    Parameters
    ----------
    ode: object
        of type :class:`DeterministicOde <pygom.model.DeterministicOde>`
    t: array like
        the time points including initial time
    full_output: bool, optional
        If the additional information from the integration is required

    '''

    # INTEGRATE!!! (shout it out loud, in Dalek voice)
    # determine the number of output we want
    solution, output = scipy.integrate.odeint(ode.ode,
                                              x0, t,
                                              Dfun=ode.jacobian,
                                              mu=None, ml=None,
                                              col_deriv=False,
                                              mxstep=10000,
                                              full_output=True)

    if full_output == True:
        # have both
        return solution, output
    else:
        return solution 

def _setupIntegrator(func, jac, x0, t0, args=(), method=None, nsteps=10000):
    if method == 'dopri5':
        # if we are going to use rk5, then one thing for sure is that we
        # know for sure that the set of equations are not stiff.
        # Furthermore, the jacobian information will never be used as
        # it evaluate f(x) directly
        r = scipy.integrate.ode(func).set_integrator('dopri5', nsteps=nsteps,
                                                     atol=atol, rtol=rtol)
    elif method == 'dop853':
        r = scipy.integrate.ode(func).set_integrator('dop853', nsteps=nsteps,
                                                     atol=atol, rtol=rtol)
    elif method == 'vode':
        r = scipy.integrate.ode(func, jac).set_integrator('vode',
                                                          with_jacobian=True,
                                                          lband=None, uband=None,
                                                          nsteps=nsteps,
                                                          atol=atol, rtol=rtol)
    elif method == 'ivode':
        r = scipy.integrate.ode(func, jac).set_integrator('vode', method='bdf',
                                                          with_jacobian=True,
                                                          lband=None, uband=None,
                                                          nsteps=nsteps,
                                                          atol=atol, rtol=rtol)
    elif method == 'lsoda':
        r = scipy.integrate.ode(func, jac).set_integrator('lsoda',
                                                          with_jacobian=True,
                                                          lband=None, uband=None,
                                                          nsteps=nsteps,
                                                          atol=atol, rtol=rtol)
    else:
        r = scipy.integrate.ode(func, jac).set_integrator('lsoda',
                                                          with_jacobian=True,
                                                          lband=None, uband=None,
                                                          nsteps=nsteps,
                                                          atol=atol, rtol=rtol)

    r.set_f_params(*args).set_jac_params(*args)
    r.set_initial_value(x0, t0)
    return r 

def test_SIR(self):
        '''
        Test the SIR model from the set of pre-defined models in common_models
        '''
        # We we wish to test another (simpler) model
        ode = common_models.SIR()

        # define the parameters
        param_eval = [
                      ('beta', 0.5),
                      ('gamma', 1.0/3.0)
                     ]

        ode.parameters = param_eval
        # the initial state, normalized to zero one
        initial_state = [1, 1.27e-6, 0]

        # evaluating the ode
        ode.ode(initial_state, 1)
        ode.jacobian(initial_state, 1)
        ode.grad(initial_state, 1)
        # b.sensitivity(sensitivity, t, state)
        ode.sensitivity(numpy.zeros(6), 1, initial_state)

        ode.linear_ode()
        # set the time sequence that we would like to observe
        t = numpy.linspace(1, 150, 100)
        # now find the solution
        _solution, output = scipy.integrate.odeint(ode.ode,
                                                   initial_state,
                                                   t,
                                                   full_output=True)

        self.assertTrue(output['message'] == 'Integration successful.')
        # Happy! :) 

def test_SEIR_periodic(self):
        '''
        Test the SEIR model from the set of pre-defined models in common_models
        '''
        ode = common_models.SEIR_Birth_Death_Periodic()
        t = numpy.linspace(0, 100, 1001)
        x0 = [0.0658, 0.0007, 0.0002, 0.]
        ode.initial_values = (x0,0)
        ode.parameters = [0.02, 35.84, 100, 1800, 0.27]
        # try to integrate to see if there is any problem
        _solution, _output = ode.integrate(t[1::], True) 

def test_SensJacobian(self):
        """
        Analytic Jacobian for the forward sensitivity equations against
        the forward differencing numeric Jacobian
        """
        # initial conditions
        s0 = np.zeros(self.d*self.p)
        ffParam = np.append(self.x0, s0)

        # integrate without using the analytical Jacobian
        solution_sens, _out = scipy.integrate.odeint(self.ode.ode_and_sensitivity,
                                                     ffParam, self.t,
                                                     full_output=True)

        # the Jacobian of the ode itself
        ff0 = solution_sens[self.index, :]
        J0 = self.ode.ode_and_sensitivity(ff0, self.t[self.index])
        p1 = self.p + 1
        dp1 = self.d*p1
        J = np.zeros((dp1, dp1))
        for i in range(dp1):
            for j in range(dp1):
                ff_temp = copy.deepcopy(ff0)
                ff_temp[j] += self.h
                J[i,j] = (self.ode.ode_and_sensitivity(ff_temp, self.t[self.index])[i] - J0[i])/self.h

        JAnalytic = self.ode.ode_and_sensitivity_jacobian(ff0, self.t[self.index])
        # JAnalytic = ode.odeAndSensitivityJacobian(ff0,t[index])

        self.assertTrue(np.allclose(J, JAnalytic)) 

def test_HessianJacobian(self):
        """
        Analytic Jacobian for the forward forward sensitivity equations
        i.e. the Hessian of the objective function against
        the forward differencing numeric Jacobian
        """
        ff0 = np.zeros(self.d*self.p*self.p)
        s0 = np.zeros(self.d*self.p)
        ffParam = np.append(np.append(self.x0, s0), ff0)
        # our integration
        sol_hess, _o = scipy.integrate.odeint(self.ode.ode_and_forwardforward,
                                              ffParam, self.t, full_output=True)

        numFF = len(ffParam)
        J = np.zeros((numFF, numFF))
        # get the info
        ff0 = sol_hess[self.index, :]
        # evaluate at target point
        J0 = self.ode.ode_and_forwardforward(ff0, self.t[self.index])
        # J0 = ode.odeAndForwardforward(ff0, t[index])
        # the Analytical solution is
        JAnalytic = self.ode.ode_and_forwardforward_jacobian(ff0, self.t[self.index])
        # JAnalytic = ode.odeAndForwardforwardJacobian(ff0, t[index])
        # now we go and find the finite difference Jacobian
        for i in range(numFF):
            for j in range(numFF):
                ff_temp = copy.deepcopy(ff0)
                ff_temp[j] += self.h
                J[i,j] = (self.ode.ode_and_forwardforward(ff_temp, self.t[self.index])[i] - J0[i])/self.h 

def cdf(self, x):
		self._check_dimension(x)
		tv = np.asarray([ scipy.stats.t.ppf(u, df=self.df) for u in x ])

		def fun(a, b):
			return multivariate_t_distribution(np.asarray([a, b]), np.asarray([0, 0]), self.R, self.df, self.dim)
			
		lim_0 = lambda x: -10
		lim_1 = lambda x: tv[1]
		return fun(x[0], x[1])
		#return scipy.integrate.dblquad(fun, -10, tv[0], lim_0, lim_1)[0] 

def lazy_quad(f, a, b, args=(), epsrel=1.49e-08, epsabs=1.49e-8, **kwargs):
    if not IS_PYPY:
        from scipy.integrate import quad
        return quad(f, a, b, args, epsrel=epsrel, epsabs=epsabs, **kwargs)
    else:
        return quad_adaptive(f, a, b, ags=args, epsrel=epsrel, epsabs=epsabs)
#        n = 300
#        return fixed_quad_Gauss_Kronrod(f, a, b, kronrod_points[n], kronrod_weights[n], legendre_weights[n], args) 

def dblquad(func, a, b, hfun, gfun, args=(), epsrel=1.48e-12, epsabs=1.48e-15):
    '''Nominally working, but trying to use it has exposed the expected bugs in 
    `quad_adaptive`.
    '''
    def inner_func(y, *args):
        full_args = (y,)+args
        quad_fluids = quad_adaptive(func, hfun(y, *args), gfun(y, *args), args=full_args, epsrel=epsrel, epsabs=epsabs)[0]
        # from scipy.integrate import quad as quad2
        # quad_sp = quad2(func, hfun(y), gfun(y), args=full_args, epsrel=epsrel, epsabs=epsabs)[0]
        # print(quad_fluids, quad_sp, hfun(y), gfun(y), full_args, )
        return quad_fluids
#         return quad(func, hfun(y), gfun(y), args=(y,)+args, epsrel=epsrel, epsabs=epsabs)[0]
    return quad_adaptive(inner_func, a, b, args=args, epsrel=epsrel, epsabs=epsabs) 

def integrate(self, axes=None, method=scipy.integrate.simps):
        '''
        Calculates the definite integral along the given axes.

        Parameters
        ----------
        method: callable
            Choose the method to use. Available options:

            * 'constant'
            * any function with the same signature as scipy.integrate.simps (default).
        '''
        if not callable(method) and method not in ['constant', 'fast']:
            raise ValueError("Requested method {} is not supported".format(method))

        if axes is None:
            axes = range(self.dimensions)

        if not isinstance(axes, Iterable):
            axes = (axes,)

        if method == 'constant':
            return self._integrate_constant(axes)
        elif method == 'fast':
            return self._integrate_fast(axes)
        else:
            return self._integrate_scipy(axes, method) 

def controllability_gramian(A, B, T = np.inf):
    '''Compute the causal controllability Gramian of the continuous time system.
    
    The system is described as
     dx = A*x + B*u
     
    T is the horizon over which to compute the Gramian. If not specified, the 
    infinite horizon Gramian is computed. Note that the infinite horizon Gramian
    only exists for asymptotically stable systems.
    
    If T is specified, we compute the Gramian as
     Wc = integrate exp(A*t)*B*B.H*exp(A.H*t) dt 
    
    Returns the matrix Wc.
    '''
    
    assert A.shape[0]==A.shape[1], "Matrix A is not square"
    assert A.shape[0]==B.shape[0], "Matrix A and B do not align"

    if not np.isfinite(T):
        #Infinite time Gramian:
        assert is_hurwitz(A), "Can only compute infinite horizon Gramian for a stable system."
        
        Wc = scipy.linalg.solve_lyapunov(A, -B.dot(B.T))
        return Wc
    
    # We need to solve the finite time Gramian
    # Boils down to solving an ODE:
    A = np.array(A,dtype=float)
    B = np.array(B,dtype=float)
    T = np.float(T)
    
    def gramian_ode(y, t0, A, B):
        temp = np.dot(scipy.linalg.expm(A*t0),B)
        dQ = np.dot(temp,np.conj(temp.T))
         
        return dQ.reshape((A.shape[0]**2,1))[:,0]
     
    y0 = np.zeros([A.shape[0]**2,1])[:,0]
    out = scipy.integrate.odeint(gramian_ode, y0, [0,T], args=(A,B))
    Q = out[1,:].reshape([A.shape[0], A.shape[0]])
    return Q 

def is_detectable(C, A):
    '''Compute whether the pair (C,A) is detectable.

    Returns True if detectable, False otherwise.
    '''

    return is_stabilisable(A.conj().T, C.conj().T)


#TODO
# def observability_gramian(A, B, T = np.inf):
#     '''Compute the observability Gramian of the continuous time system.
#     
#     The system is described as
#      dx = A*x + B*u
#      
#     T is the horizon over which to compute the Gramian. If not specified, the 
#     infinite horizon Gramian is computed. Note that the infinite horizon Gramian
#     only exists for asymptotically stable systems.
#     
#     If T is specified, we compute the Gramian as
#      Wc = integrate exp(A*t)*B*B.H*exp(A.H*t) dt 
#     
#     Returns the matrix Wc.
#     '''
#     
#     assert A.shape[0]==A.shape[1], "Matrix A is not square"
#     assert A.shape[0]==B.shape[0], "Matrix A and B do not align"
# 
#     if not np.isfinite(T):
#         #Infinite time Gramian:
#         eigVals, eigVecs = scipy.linalg.eig(A)
#         assert np.max(np.real(eigVals)) < 0, "Can only compute infinite horizon Gramian for a stable system."
#         
#         Wc = scipy.linalg.solve_lyapunov(A, -B*B.T)
#         return Wc
#     
#     # We need to solve the finite time Gramian
#     # Boils down to solving an ODE:
#     A = np.array(A,dtype=float)
#     B = np.array(B,dtype=float)
#     T = np.float(T)
#     
#     def gramian_ode(y, t0, A, B):
#         temp = np.dot(scipy.linalg.expm(A*t0),B)
#         dQ = np.dot(temp,np.conj(temp.T))
#          
#         return dQ.reshape((A.shape[0]**2,1))[:,0]
#      
#     y0 = np.zeros([A.shape[0]**2,1])[:,0]
#     out = scipy.integrate.odeint(gramian_ode, y0, [0,T], args=(A,B))
#     Q = out[1,:].reshape([A.shape[0], A.shape[0]])
#     return Q 

def PlanToEndEffectorPose(self, robot, goal_pose, timelimit=5.0,
                              pose_error_tol=0.01,
                              integration_interval=10.0,
                              **kw_args):
        """
        Plan to an end effector pose by following a geodesic loss function
        in SE(3) via an optimized Jacobian.

        @param robot
        @param goal_pose desired end-effector pose
        @param timelimit time limit before giving up
        @param pose_error_tol in meters
        @param integration_interval The time interval to integrate over
        @return traj
        """
        manip = robot.GetActiveManipulator()

        def vf_geodesic():
            """
            Define a joint-space vector field, that moves along the
            geodesic (shortest path) from the start pose to the goal pose.
            """
            twist = util.GeodesicTwist(manip.GetEndEffectorTransform(),
                                       goal_pose)
            dqout, tout = util.ComputeJointVelocityFromTwist(
                robot, twist, joint_velocity_limits=numpy.PINF)

            # Go as fast as possible
            vlimits = robot.GetDOFVelocityLimits(robot.GetActiveDOFIndices())
            return min(abs(vlimits[i] / dqout[i])
                       if dqout[i] != 0. else 1.
                       for i in xrange(vlimits.shape[0])) * dqout

        def CloseEnough():
            """
            The termination condition.
            At each integration step, the geodesic error between the
            start and goal poses is compared. If within threshold,
            the integration will terminate.
            """
            pose_error = util.GetGeodesicDistanceBetweenTransforms(
                manip.GetEndEffectorTransform(), goal_pose)
            if pose_error < pose_error_tol:
                return Status.TERMINATE
            return Status.CONTINUE

        traj = self.FollowVectorField(robot, vf_geodesic, CloseEnough,
                                      integration_interval,
                                      timelimit,
                                      **kw_args)

        # Flag this trajectory as unconstrained. This overwrites the
        # constrained flag set by FollowVectorField.
        util.SetTrajectoryTags(traj, {Tags.CONSTRAINED: False}, append=True)
        return traj