python类odeint()的实例源码

def likelihood(parameter_vector):

    parameter_vector = 10**np.array(parameter_vector)

    #Solve ODE system given parameter vector
    yout = odeint(odefunc, y0, tspan, args=(parameter_vector,))

    cout = yout[:, 2]

    #Calculate log probability contribution given simulated experimental values.

    logp_ctotal = np.sum(like_ctot.logpdf(cout))

    #If simulation failed due to integrator errors, return a log probability of -inf.
    if np.isnan(logp_ctotal):
        logp_ctotal = -np.inf

    return logp_ctotal


# Add vector of rate parameters to be sampled as unobserved random variables in DREAM with uniform priors.

def velocity(stateVec, t):
    """
    return the velocity field of Lorentz system.
    stateVec : the state vector in the full space. [x, y, z]
    t : time is used since odeint() requires it. 
    """

    x = stateVec[0]
    y = stateVec[1]
    z = stateVec[2]

    # complete the flowing 3 lines.
    vx =  G_sigma*(y - x)
    vy = G_rho*x - y - x*z
    vz = x*y - G_b*z

    return np.array([vx, vy, vz])

def velocity(stateVec, t):      
    """
    velocity in the full state space.

    stateVec: state vector [x1, y1, x2, y2]
    t: just for convention of odeint, not used.
    return: velocity at stateVec. Dimension [1 x 4]
    """
    x1 = stateVec[0]
    y1 = stateVec[1]
    x2 = stateVec[2]
    y2 = stateVec[3]

    r2 = x1**2 + y1**2

    velo = np.array([(G_mu1-r2) * x1 + G_c1 * (x1*x2 + y1*y2),
             (G_mu1-r2) * y1 + G_c1 * (x1*y2 - x2*y1),
             x2 + y2 + x1**2 - y1**2 + G_a2 * x2 * r2,
             -x2 + y2 + 2.0 * x1 * y1 + G_a2 * y2 * r2])

    return velo

def solve_coupled_constecc_solution(F0, e0, phase0, mc, t):
    """
    Compute the solution to the coupled system of equations
    from from Peters (1964) and Barack & Cutler (2004) at
    a given time.

    :param F0: Initial orbital frequency [Hz]
    :param mc: Chirp mass of binary [Solar Mass]
    :param t: Time at which to evaluate solution [s]

    :returns: (F(t), phase(t))
    """

    y0 = np.array([F0, phase0])

    y, infodict = odeint(get_coupled_constecc_eqns, y0, t,
                         args=(mc,e0), full_output=True)

    if infodict['message'] == 'Integration successful.':
        ret = y
    else:
        ret = 0

    return ret

def step(self, action):
        action = self._discrete_actions[action[0]]
        sa = np.append(self._state, action)
        new_state = odeint(self._dpds, sa, [0, self._dt])

        self._state = new_state[-1, :-1]

        if self._state[0] < -self.max_pos or \
                np.abs(self._state[1]) > self.max_velocity:
            reward = -1
            absorbing = True
        elif self._state[0] > self.max_pos and \
                np.abs(self._state[1]) <= self.max_velocity:
            reward = 1
            absorbing = True
        else:
            reward = 0
            absorbing = False

        return self._state, reward, absorbing, {}

def step(self, action):
        action = self._discrete_actions[action[0]]
        action += np.random.uniform(-10., 10.)
        sa = np.append(self._state, action)
        new_state = odeint(self._dpds, sa, [0, self._dt])

        self._state = new_state[-1, :-1]
        self._state[0] = self._range_pi(self._state[0])

        if np.abs(self._state[0]) > np.pi / 2.:
            reward = -1
            absorbing = True
        else:
            reward = 0
            absorbing = False

        return self._state, reward, absorbing, {}

def Wdesired(x):
        ''' x is the augmented variable represented in the network 
            it obeys \tau_s x_\alpha = -x_\alpha + Wdesired_\alpha(x) 
            x is related to the original augmented variable \tilde{x} by x_\alpha = varFactors_\alpha \tilde{x}_\alpha 
            where varFactors_alpha = angleFactor | velocityFactor | torqueFactor
            now, original augmented variable obeys \dot{\tilde{x}}=f(\tilde{x})
            so, we have Wdesired_\alpha(x) = \tau_s * varFactor_alpha * f_\alpha(\tilde{x}) + x
        '''
        # \tilde{x} (two zeroes at x[N:N+N//2] are ignored
        xtilde = x/varFactors
        if XY: angles = armAngles(xtilde[:N])
        else: angles = xtilde[:N//2]
        # f(\tilde{x}), \dot{u} part is not needed
        qdot,dqdot = evolveFns(angles,xtilde[Nobs-N//2:Nobs],xtilde[Nobs:],XY,dt)
                                                                        # returns deltaposn if XY else deltaangles
        # \tau_s * varFactors_alpha * f_\alpha(\tilde{x}) + x
        return np.append(np.append(qdot,dqdot),np.zeros(N//2))*varFactors*tau + x
                                                                        # integral on torque u also
                                                                        # VERY IMP to compensate for synaptic decay on torque
        #return np.append( np.append(qdot,dqdot)*tau*varFactors[:Nobs] + x[:Nobs], np.zeros(N//2) )
                                                                        # normal synaptic decay on torque u

    ##### For the reference, choose EITHER robot simulation rateEvolveProbe above
    #####  OR evolve Wdesired inverted / evolveFns using odeint as below -- both should be exactly same

def reaction(y0, t, k0, k0r):
    """
    Wrapper for the reaction.
    It receives an `np.array` `y0` with the initial concentrations and an
    `np.array` `t` with timestamps (it should also include `0`).
    This function solves the corresponding ODE system and returns an `np.array`
    `Y` in which each column represents a chemical species and each line a
    timestamp.
    """
    def dydt(y, t):
        return np.array([
                         -1*v_0(y[0], y[1]),
                         +1*v_0(y[0], y[1]),
                         ])

    # A <-> B
    def v_0(A, B):
        return k0 * A**1 - k0r * B**1

    return odeint(dydt, y0, t)

# Reaction rates:
# A <-> B

def reaction(y0, t, k0, k0r):
    """
    Wrapper for the reaction.
    It receives an `np.array` `y0` with the initial concentrations and an
    `np.array` `t` with timestamps (it should also include `0`).
    This function solves the corresponding ODE system and returns an `np.array`
    `Y` in which each column represents a chemical species and each line a
    timestamp.
    """
    def dydt(y, t):
        return np.array([
                         -2*v_0(y[0], y[1]),
                         +1*v_0(y[0], y[1]),
                         ])

    # 2*A <-> C
    def v_0(A, C):
        return k0 * A**2 - k0r * C**1

    return odeint(dydt, y0, t)

# Reaction rates:
# 2*A <-> C

def reaction(y0, t, k0, k0r):
    """
    Wrapper for the reaction.
    It receives an `np.array` `y0` with the initial concentrations and an
    `np.array` `t` with timestamps (it should also include `0`).
    This function solves the corresponding ODE system and returns an `np.array`
    `Y` in which each column represents a chemical species and each line a
    timestamp.
    """
    def dydt(y, t):
        return np.array([
                         -1*v_0(y[1], y[0], y[2]),
                         +2*v_0(y[1], y[0], y[2]),
                         +2*v_0(y[1], y[0], y[2]),
                         ])

    # S <-> 2*A + 2*C
    def v_0(A, S, C):
        return k0 * S**1 - k0r * A**2 * C**2

    return odeint(dydt, y0, t)

# Reaction rates:
# S <-> 2*A + 2*C

def raw_sim_one(start, steps, dy_data, settings, include_step_zero=False):
    '''
    simulate from a single point at the given times

    return an nparray of states at those times, possibly excluding time zero
    '''

    start.shape = (dy_data.num_dims,)

    times = np.linspace(0, settings.step * steps, num=steps+1)

    der_func = dy_data.make_der_func()
    jac_func, max_upper, max_lower = dy_data.make_jac_func()
    sim_tol = settings.sim_tol

    result = odeint(der_func, start, times, Dfun=jac_func, col_deriv=True, \
                    mxstep=int(1e8), mu=max_upper, ml=max_lower, \
                    atol=sim_tol, rtol=sim_tol)

    if not include_step_zero:
        result = result[1:]

    return result

def sim(start, der_func, time_amount, num_steps, quick):
    'simulate for some fixed time, and return the resultant (states, times) tuple'

    tol = 1.49012e-8

    if not quick:
        tol /= 1e5 # more accurate simulation

    times = np.linspace(0, time_amount, num=1 + num_steps)
    states = odeint(der_func, start, times, col_deriv=True, rtol=tol, atol=tol, mxstep=50000)

    states = states[1:]
    times = times[1:]
    assert len(states) == num_steps
    assert len(times) == num_steps

    return (states, times)

def run(self,t=1000,ts=100000,plot=False):
        """
        Run it for time t with ts number of time steps
        Outputs the concentration profile for the run
        default values are 100000 and 1000000
        """ 
        t_index = np.linspace(0, t, ts)
        y = self.concentrations
        output = odeint(odes, y, t_index, args= (self.reactions,self.rates))
        self.concentrations = output[-1,:]
        if plot:
            for i in xrange(output.shape[1]):
                label = self.species[i]
                plt.plot(t_index, output[:, i], label=label)
            plt.legend(loc='best')
            plt.xlabel('t')
            plt.title('Concentration Profile')
            plt.grid()
            plt.show()
        return output

def solveUnst(self,t):
    solution = [None]*len(self.mapping)
    # This has the unsteady part
    # Solver, just live and let live  
    for ix, (i,k) in enumerate(self.mapping):
      solution[ix] = self.b[i][k]
    soln = odeint(self.rUnst, solution, t,hmax=(t[-1]-t[0])/len(t), \
      rtol = 1e-4, atol = 1e-4)

    # final update
    self.updateUnst(t[-1])
    self.update(soln[-1,:])

    # Lets collect all the steps
    fullSolution = dict([(b.name + '_'+s,[]) for b in self.b for s in b.state])
    for j in range(0,len(t)):
      for ix, (i,k) in enumerate(self.mapping):
        fullSolution[self.b[i].name+'_'+k].append(soln[j,ix])
    fullSolution['t'] = t
    return fullSolution

def newton_rad_func(E_val,D,m,L0):

   E = E_val

   c2 = D*D*E/4.

   L = newton(newton_ang_func,L0,args=(c2,m),tol=1e-8,maxiter=500)

   slope = -(-L+m*(m+1.)+2.*D+c2)/(2.*(m+1.))

   z0 = [1+step*slope,slope]

   z = odeint(g,z0,x_rad,args=(c2,L,D,m))

   temp=pow(x_rad,2.0)-1.
   temp=pow(temp,m/2.)

   zz=temp*z[:,0]

   first_zz = np.array([1])
   zz=np.append(first_zz, zz)

   return z[:,1][-1]

def plottrajectories(fs, x0, t=np.linspace(1,400,10000), **kw):
    """
    plots trajectory of the solution

    f  -- must accept an array of X and t=0, and return a 2D array of \dot y and \dot x
    x0 -- vector

    Example
    =======

    plottrajectories(lambda X,t=0:array([ X[0] -   X[0]*X[1] ,
                   -X[1] + X[0]*X[1] ]), [ 5,5], color='red')
    """
    x0 = np.array(x0)
    #f = lambda X,t=0: array(fs(X[0],X[1]))
    #fa = lambda X,t=0:array(fs(X[0],X[1]))
    X = integrate.odeint( fs, x0, t)
    plt.plot(X[:,0], X[:,1], **kw)

def solve(self):
        self.psi = odeint(self.RichardsEquation, self.psi0, self.t, args=(
            self.dz, self.n, self.p, self.qTop, self.qBot, self.psiTop, self.psiBot), mxstep=500)
        self.psi0 = self.psi[-1, :]

def calculate_solution(self):
        solution = odeint(func=self.set_up_equation, y0=self.initial_condition, t=self.steps,
                          args=(self.coefficient_matrix, self.steps))
        return solution

def update(self, dt, F, M):
        # limit thrust and Moment
        L = params.arm_length
        r = params.r
        prop_thrusts = params.invA.dot(np.r_[np.array([[F]]), M])
        prop_thrusts_clamped = np.maximum(np.minimum(prop_thrusts, params.maxF/4), params.minF/4)
        F = np.sum(prop_thrusts_clamped)
        M = params.A[1:].dot(prop_thrusts_clamped)
        self.state = integrate.odeint(self.state_dot, self.state, [0,dt], args = (F, M))[1]

def TC_Simulate(Mode,initialCondition,time_bound):
    time_step = 0.05;
    time_bound = float(time_bound)
    initial = [float(tmp)  for tmp in initialCondition]
    number_points = int(np.ceil(time_bound/time_step))
    t = [i*time_step for i in range(0,number_points)]
    if t[-1] != time_step:
        t.append(time_bound)

    y_initial = initial[0]

    if Mode == 'On':
        rate = 0.1
    elif Mode == 'Off':
        rate = -0.1
    else:
        print('Wrong Mode name!')
    sol = odeint(thermo_dynamic,y_initial,t,args=(rate,),hmax = time_step)

    # Construct the final output
    trace = []
    for j in range(len(t)):
        #print t[j], current_psi
        tmp = []
        tmp.append(t[j])
        tmp.append(sol[j,0])
        trace.append(tmp)
    return trace


# sol = TC_Simulate('Off',[60],10)
# print(sol)

def integrator(init_x, dt, nstp):
    """
    The integator of the Lorentz system.
    init_x: the intial condition
    dt : time step
    nstp: number of integration steps.

    return : a [ nstp x 3 ] vector 
    """

    state = odeint(velocity, init_x, np.arange(0, dt*nstp, dt))
    return state

def Jacobian(ssp, dt, nstp):
    """
    Jacobian function for the trajectory started on ssp, evolved for time t

    Inputs:
    ssp: Initial state space point. dx1 NumPy array: ssp = [x, y, z]
    t: Integration time
    Outputs:
    J: Jacobian of trajectory f^t(ssp). dxd NumPy array
    """
    #CONSTRUCT THIS FUNCTION
    #Hint: See the Jacobian calculation in CycleStability.py
    #J = None
    Jacobian0 = np.identity(3)  # COMPLETE THIS LINE. HINT: Use np.identity(DIMENSION)
    #Initial condition for Jacobian integral is a d+d^2 dimensional matrix
    #formed by concatenation of initial condition for state space and the
    #Jacobian:
    sspJacobian0 = np.zeros(3 + 3 ** 2)  # Initiate
    sspJacobian0[0:3] = ssp  # First 3 elemenets
    sspJacobian0[3:] = np.reshape(Jacobian0, 9)  # Remaining 9 elements
    tInitial = 0  # Initial time
    tFinal = dt*nstp  # Final time
    Nt = nstp  # Number of time points to be used in the integration

    tArray = np.linspace(tInitial, tFinal, Nt)  # Time array for solution

    sspJacobianSolution = odeint(JacobianVelocity, sspJacobian0, tArray)

    xt = sspJacobianSolution[:, 0]  # Read x(t)
    yt = sspJacobianSolution[:, 1]  # Read y(t)
    zt = sspJacobianSolution[:, 2]  # Read z(t)

    #Read the Jacobian for the periodic orbit:
    J = sspJacobianSolution[-1, 3:].reshape((3, 3))

    return J

def integrator(init_state, dt, nstp):
    """
    integrate two modes system in the full state sapce.

    init_state: initial state [x1, y1, x2, y2]
    dt: time step 
    nstp: number of time step
    """
    states = odeint(velocity, init_state, np.arange(0, dt*nstp, dt))
    return states

def solve_coupled_ecc_solution(F0, e0, gamma0, phase0, mc, q, t):
    """
    Compute the solution to the coupled system of equations
    from from Peters (1964) and Barack & Cutler (2004) at
    a given time.

    :param F0: Initial orbital frequency [Hz]
    :param e0: Initial orbital eccentricity
    :param gamma0: Initial angle of precession of periastron [rad]
    :param mc: Chirp mass of binary [Solar Mass]
    :param q: Mass ratio of binary
    :param t: Time at which to evaluate solution [s]

    :returns: (F(t), e(t), gamma(t), phase(t))
    """

    y0 = np.array([F0, e0, gamma0, phase0])

    y, infodict = odeint(get_coupled_ecc_eqns, y0, t,
                         args=(mc,q), full_output=True)

    if infodict['message'] == 'Integration successful.':
        ret = y
    else:
        ret = 0

    return ret

def matevolve2(y,t):
        ''' the reference y is only N-dim i.e. (q,dq), not 2N-dim, as inpfn is used directly as reference for torque u
        '''
        ## invert the nengo function transformation with angleFactor, tau, +x, etc. in Wdesired()
        ## -- some BUG, inversion is not working correctly
        #xfull = np.append(y,inpfn(t))*varFactors
        #return ( (Wdesired(xfull)[:N]/tau/varFactors[:N] - xfull[:N]) \
        #                        if (t%Tperiod)<(Tperiod-Tclamp) else -x/tau )
        # instead of above, directly use evolveFns()
        #########  DOESN'T WORK: should only use armAngles() with valid posn, not all posn-s are valid for an arm!!!  ###########
        if XY: angles = armAngles(y[:N])
        else: angles = y[:N//2]        
        # evolveFns returns deltaposn if XY else deltaangles
        if trialClamp:
            return ( evolveFns( angles, y[Nobs-N//2:Nobs], inpfn(t), XY, dt).flatten()\
                        if (t%Tperiod)<(Tperiod-Tclamp) else -y/dt )
        else:
            return evolveFns( angles, y[Nobs-N//2:Nobs], inpfn(t), XY, dt).flatten()

    ##### uncomment below to override rateEvolveProbe by matevolve2-computed Wdesired-inversion / evolveFns-evolution, as reference signal
    #trange = np.arange(0.0,Tmax,dt)
    #y = odeint(matevolve2,0.001*np.ones(N),trange,hmax=dt)  # set hmax=dt, to avoid adaptive step size
    #                                                        # some systems (van der pol) have origin as a fixed pt
    #                                                        # hence start just off-origin
    #rateEvolveProbe = y                                     # only copies pointer, not full array (no need to use np.copy() here)

###
### Reference evolution used when copycat layer is not used for reference ###
###

# scale the output of the robot simulation or odeint to cover the representation range of the network
# here I scale by angle/velocity factors, below at nodeIn I scale by torque factors.

def reaction(y0, t, k0, k0r, k1, k1r):
    """
    Wrapper for the reaction.
    It receives an `np.array` `y0` with the initial concentrations and an
    `np.array` `t` with timestamps (it should also include `0`).
    This function solves the corresponding ODE system and returns an `np.array`
    `Y` in which each column represents a chemical species and each line a
    timestamp.
    """
    def dydt(y, t):
        return np.array([
                         -2*v_0(y[0], y[2], y[1])-1*v_1(y[0], y[2], y[3]),
                         -1*v_0(y[0], y[2], y[1]),
                         +1*v_0(y[0], y[2], y[1])-1*v_1(y[0], y[2], y[3]),
                         +1*v_1(y[0], y[2], y[3]),
                         ])

    # 2*A + B <-> C
    def v_0(A, C, B):
        return k0 * A**2 * B**1 - k0r * C**1

    # A + C <-> D
    def v_1(A, C, D):
        return k1 * A**1 * C**1 - k1r * D**1

    return odeint(dydt, y0, t)

# Reaction rates:
# 2*A + B <-> C

def reaction(y0, t, k0, k0r, k1, k1r, k2, k2r):
    """
    Wrapper for the reaction.
    It receives an `np.array` `y0` with the initial concentrations and an
    `np.array` `t` with timestamps (it should also include `0`).
    This function solves the corresponding ODE system and returns an `np.array`
    `Y` in which each column represents a chemical species and each line a
    timestamp.
    """
    def dydt(y, t):
        return np.array([
                         -1*v_0(y[1], y[0], y[2])+1*v_1(y[3], y[0], y[2]),
                         -1*v_0(y[1], y[0], y[2])-1*v_2(y[3], y[1]),
                         +1*v_0(y[1], y[0], y[2])-1*v_1(y[3], y[0], y[2]),
                         +1*v_1(y[3], y[0], y[2])+1*v_2(y[3], y[1]),
                         ])

    # E + S <-> ES
    def v_0(S, E, ES):
        return k0 * E**1 * S**1 - k0r * ES**1

    # ES <-> E + P
    def v_1(P, E, ES):
        return k1 * ES**1 - k1r * E**1 * P**1

    # S <-> P
    def v_2(P, S):
        return k2 * S**1 - k2r * P**1

    return odeint(dydt, y0, t)

# Reaction rates:
# E + S <-> ES

def reaction(y0, t, k0, k0r, k1, k1r):
    """
    Wrapper for the reaction.
    It receives an `np.array` `y0` with the initial concentrations and an
    `np.array` `t` with timestamps (it should also include `0`).
    This function solves the corresponding ODE system and returns an `np.array`
    `Y` in which each column represents a chemical species and each line a
    timestamp.
    """
    def dydt(y, t):
        return np.array([
                         -1*v_0(y[1], y[0], y[2])+1*v_1(y[3], y[0], y[2]),
                         -1*v_0(y[1], y[0], y[2]),
                         +1*v_0(y[1], y[0], y[2])-1*v_1(y[3], y[0], y[2]),
                         +1*v_1(y[3], y[0], y[2]),
                         ])

    # E + S <-> ES
    def v_0(S, E, ES):
        return k0 * E**1 * S**1 - k0r * ES**1

    # ES <-> E + P
    def v_1(P, E, ES):
        return k1 * ES**1 - k1r * E**1 * P**1

    return odeint(dydt, y0, t)

# Reaction rates:
# E + S <-> ES

def _transform_eigenfunction(self):

        eigenfunction = si.odeint(self._ff, self._init_state_vect, self._domain)

        return [eigenfunction[:, 0], eigenfunction[:, 1], eigenfunction[:, 2], eigenfunction[:, 3]]

def run_till(self,delta_time=0.001,eps = 10**-12,every=10, stop=100000):
        """
        Runs till the gradients are less than eps
        returns time taken with array of concentrations
        default values of time_step and eps are 10**-3 and 10**-12 
        """
        y = self.concentrations
        output = [y]
        t=0
        t_index = np.linspace(0,every,int(every/(delta_time)))
        if every>10*delta_time:
             while True:
                gy = odes(y,t,self.reactions,self.rates)
                if (np.abs(gy) <eps).all():
                    break
                o = odeint(odes, y, t_index, args= (self.reactions,self.rates))
                y = o[-1,:]
                t+=every
                output = output + o[1:]
                if t>stop:
                    break 
        else:
            while True:
                gy = odes(y,t,self.reactions,self.rates)
                if (np.abs(gy) <eps).all():
                    break
                y+= gy*delta_time
                t+=delta_time
                output.append(y)
                if t>stop:
                    break

        output = np.array(output)
        self.concentrations =y
        return output,t