Difference between revisions of "Team:EPF Lausanne/Software"

"""
Copyright (C) 2015 iGEM Team EPF_Lausanne

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>
"""

import numpy as np
from scipy.integrate import *

class Solver:
    """
    Class that allows the solution of a system of non-linear ODEs. The system is specified by the function fun

        dy/dt = fun(t,y)

    where t is a number and y and dy/dt are numpy arrays or lists.

    The solution is performed with the dopri5 method, an explicit Runge-Kutta method of order (4)5.
    The method is due to Dormand & Prince, and is implemented by E. Hairer and G. Wanner.

    See
        http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.ode.html
    for more details.

    NOTE:
    Our Solver can take a function of the form

        f(t,y,pars)

    where PARS are parameters. PARS can be eventually passed to the constructor of the Solver.
    """

    def __init__(self,dt,fun,t0,T,y0,pars=[]):
        self.dt = dt # Time step
        self.fun = fun # Function representing the ODE
        self.t0 = t0 # Initial time
        self.T = T # Final time
        self.y0 = y0 # Initial condition

        self.pars = pars # Parameters of the system of ODEs

    def solve(self):
        """
        Solve the system of ODEs

            dy/dt = fun(t,y)

        on the interval [t0,T], with the initial condition y(0)=y0.

        Returns two lists, time and solution, containing time points and the solution at these time points.
        """

        # Choose integrator type
        r = ode(self.fun).set_integrator('dopri5')

        # Initialize the integrator
        r.set_initial_value(self.y0, self.t0)

        # Set parameters for the ODE function
        r.set_f_params(*self.pars)

        # Initialize solution list and time points list
        solution = np.asarray(self.y0)
        time = np.asarray(self.t0)

        while r.successful() and r.t < self.T:
            r.integrate(r.t + self.dt) # Perform one integration step, i.e. obtain the solution y at time t+dt

            time = np.append(time,r.t) # Append the new time
            solution = np.vstack((solution,r.y)) # Append the new solution

        return time, solution # Return time and solution vectors

    def solve_for_t(self,t):
        """
        Solve the system of ODEs

            dy/dt = fun(t,y)

        on the interval [t0,T], with the initial condition y(0)=y0.

        Returns two lists, solution and time, containing time points and the solution at these time points.

        The solution is computed at the points specified in t, i.e. the time step dt is ignored.
        """

        # Choose integrator type: dopri5 in this case
        r = ode(self.fun).set_integrator('dopri5')

        # Initialize the integrator
        r.set_initial_value(self.y0, self.t0)

        # Set parameters for the ODE function
        r.set_f_params(*self.pars)

        # Initialize solution list and time points list
        solution = []
        time = []

        for tt in t:
            r.integrate(tt) # Perform one integration step

            time.append(tt) # Append the new time
            solution.append(r.y) # Append the new solution

        return np.asarray(time), np.asarray(solution) # Return time and solution vectors


if __name__ == "__main__":
    """
    Our test functions:
       rapid_equilibrium (standard function)
       rapid_equilibrium_from_string() (returns a function compiled from a string)
    """

    import matplotlib.pylab as plt
    from test import * # Import test functions for the ODE integrator

    dt = 0.1

    t0 = 0
    T = 100
    y0 = [1,0,0]

    # Store the funtion compiled from a string
    rapid_equilibrium_s = rapid_equilibrium_from_string()

    mysolver = Solver(dt,rapid_equilibrium,t0,T,y0)
    mysolver_string = Solver(dt,rapid_equilibrium_s,t0,T,y0)

    t,y = mysolver.solve()
    tt,yt = mysolver.solve_for_t(np.linspace(t0,T,10))
    ts,ys = mysolver_string.solve()

    plt.plot(t,y)
    plt.plot(tt,yt,'x')
    plt.plot(ts,ys)
    plt.show()