---

integrator.h

/***************************************************************************
                          integrator.h  -  description
                             -------------------
    begin                : Tue Aug 12 2000
    copyright            : (C) 2000 by Michael Peeters
    email                : Michael.Peeters@vub.ac.be
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef INTEGRATOR_H
#define INTEGRATOR_H 

#include "numerictypes.h"
#include "numerictraits.h"
#include "vectorfunction.h"
// #include "odesystem.h"
#include "random.h"
#include "jacobian.h"

namespace MODEL {

  // predef
  template<integer dims, typename nelem, class NT>
  class ODESystem;

  enum Integration {Euler, RungeKutta, Milshtein,Heun};

  template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
  class Integrator
  {
        friend class ODESystem<dims,nelem,NT>;

  public:
        typedef typename NT::number                                     numT;
        typedef typename NT::vect                                       vect;
        typedef typename NT::vf                                         vf;
        typedef typename NT::matrix                                     matrix;

        Integrator(vf& d, vf& s)
          :  deter(&d), stoch(&s){}
        
        virtual ~Integrator(){};

        virtual vect& step(vect& current, time& dt)=0;

        vect&   currentpoint(void) {return current;}
                
  private:
        vf* get_deterministic(void)
        {
          return deter;
        }

        vf* get_stochastic(void)
        {
          return stoch;
        }

  protected:
        vf* deter;
        vf* stoch;
        
  };    
  //------------------------------------------------------------

  template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
  class IEuler : public Integrator<dims,nelem,NT>
  {
  public:
        IEuler(vf& d, vf& s) 
          : Integrator<dims,nelem,NT>(d,s) {}

        virtual vect& step(vect& current, time& dt)
        {
          vect noise;
          for(counter i=0;i<dims;i++) noise[i]=rnd();
          vect g=(*stoch)(current);
          for(counter i=0;i<dims;i++) g[i]*=noise[i];

          return current+=dt*(*deter)(current)+sqrt(dt)*g;
        }

  private:
        Normal rnd; // will be seeded with time()

  };

        
  //-------------------------------------------------------------------------
  template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
  class IRungeKutta : public Integrator<dims,nelem,NT>
  {
  public:
        IRungeKutta(vf& d, vf& s) 
          : Integrator<dims,nelem,NT>(d,s) {}

        virtual vect& step(vect& current, time& dt)
        {
          // Runge-Kutta needs 4 evaluations:
          // 1) at the current point
          // 2) at a halfway point, using the values of 1)
          // 3) at a halfway point, using the values of 2)
          // 4) at the end point, using the values of 3)
          
          vect k1,k2,k3,k4;
          number dt2=dt/2.;

          // Step 1) 

          k1=(*deter)(current);
          k1*=dt2;

          // Step 2)
          // halfway in t

          vect pos(k1);
          pos+=current; //halfway in u
 
          k2=(*deter)(pos);
          k2*=dt2;

          // Step 3)
          pos=k2;
          pos+=current;

          k3=(*deter)(pos);
          k3*=dt2;
          // Step 4)
          //another half step

          pos=k3; pos+=current; // BUG FIX: Forgot to add current point
          k4=(*deter)(pos);
          k4*=dt2;

          // Now calculate next point. Caveat
          // emptor, this is very slow (compared to what it could be).
          return current+=(k1+2.*k2+2.*k3+k4)/3.;
        }

  };

  //----------------------------------------------------------------------              
  template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
  class IMilshtein : public Integrator<dims,nelem,NT> 
  {
  public:
        IMilshtein(vf& d, vf& s) 
          : Integrator<dims,nelem,NT>(d,s){}

        virtual vect& step(vect& current, time& dt)
        {
          // The Milshtein method is the basic algo for noisy problems
          // => Needs a random number generator !

          vect oldcurrent(current);

          // Deterministix q
          vect q=(*deter)(oldcurrent);
          current+=dt*q;
                  
          // Stochastix g
          // Generate numbers: noise is _NOT_ correlated between modes.
          vect noise;
          for(counter i=0;i<dims;i++) noise[i]=rnd();
                  
          // Additive component (Euler)
          number sqdt=sqrt(dt);
          vect g=(*stoch)(oldcurrent);
                  
          // Multiply each component with its noise
          vect ng(g);
          for(counter i=0;i<dims;i++) ng[i]*=noise[i]; 
                  
          // Result
          current+=sqdt*ng;

          // Extra Millstein 
          Jacobian<dims,NT> J(*stoch);
          matrix dgdy=J.calculate(oldcurrent,g);
                  
          vect dgs(0.);
          for(counter i=0;i<dims;i++)
                for (counter j=0;j<dims;j++)
                  dgs[i]+=0.5*dgdy[i][j]*g[j]*noise[i]*noise[j];
        
          current+=dt*dgs;
          return current;
        }

  private:              
        Normal rnd; // will be seeded with time()
  };

  template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
  class IHeun : public Integrator<dims,nelem,NT> 
  {
  public:
        IHeun(vf& d, vf& s) 
          : Integrator<dims,nelem,NT>(d,s){}

        virtual vect& step(vect& current, time& dt)
        {
          // The Heun method is the best easy algo for noisy problems
          // => Needs a random number generator !
          
          // u is N(0,1)

          // k=h.q(t,x)
          // l=sqrth.u.g(t,x)
          //
          // qavg = (1/2)*( q(t,x) + q(t+h,x+k+l) )
          // gavg = (1/2)*( g(t,x) + g(t+h,x+k+l) )
          //
          // Watch out: we might have to take a step for the modulators
          //
          // xnew = x + h*qavg + sqrth*u*gavg

          vect oldcurrent(current);

          // Some to be optimized out temps
          const number h=dt;
          const number sh=sqrt(dt);

          vect q = (*deter)(oldcurrent);
          vect g = (*stoch)(oldcurrent);

          // Generate numbers: noise is _NOT_ correlated between modes.
          vect u;
          for(counter i=0;i<dims;i++) u[i]=rnd();
          
          // uneff: vect k=h*q;
          // uneff: vect l=sh*g;

          // generate next point
          for(counter i=0;i<dims;i++)
                {
                  oldcurrent[i]+=h*q[i]+sh*g[i]*u[i];
                  // The stochastic part should be sqrt(2D) when the autocorrelation
                  // is given by <F,F>=2D.delta; 2D is the variance (in fact)
                }

          // uneff:               l[i]*=u[i];

          // uneff: oldcurrent+=k;
          // uneff: oldcurrent+=l;
          
          // vect qavg=0.5* h * ( q + (*deter)(oldcurrent) );
          // vect gavg=0.5* sh * ( g + (*stoch)(oldcurrent) );

          vect qn=(*deter)(oldcurrent);
          vect gn=(*stoch)(oldcurrent);

          for(counter i=0;i<dims;i++) 
                {
                  current[i]+=0.5*( h*( q[i] + qn[i]) + sh*u[i]*(g[i]+gn[i]));
                  
                  // uneff:  gavg[i]*=u[i];
                }

          // uneff: current += qavg;
          // uneff: current += gavg;
          
          return current;
        }

  private:              
        Normal rnd; // will be seeded with time()
  };

  template<integer dims, typename nelem=number, class NT = NumericTraits<nelem,dims> >
  class IHeunSimpleCorr : public Integrator<dims,nelem,NT> 
  {
  public:
        IHeunSimpleCorr(vf& d, vf& s, vf& transform) 
          : Integrator<dims,nelem,NT>(d,s), trans(&transform){}

        virtual vect& step(vect& current, time& dt)       
        {         
        // The Heun method is the best easy algo for noisy problems
          // => Needs a random number generator !
          
          // u is N(0,1)

          // k=h.q(t,x)
          // l=sqrth.u.g(t,x)
          //
          // qavg = (1/2)*( q(t,x) + q(t+h,x+k+l) )
          // gavg = (1/2)*( g(t,x) + g(t+h,x+k+l) )
          //
          // Watch out: we might have to take a step for the modulators
          //
          // xnew = x + h*qavg + sqrth*u*gavg

          vect oldcurrent(current);

          // Some to be optimized out temps
          const number h=dt;
          const number sh=sqrt(dt);

          vect q = (*deter)(oldcurrent);
          vect g = (*stoch)(oldcurrent);

          // Generate numbers: noise is _NOT_ correlated between modes.
          vect u,u0;
          for(counter i=0;i<dims;i++) u0[i]=rnd();

          // Do transform for eventual correlation
          trans->function(u,u0);
          
          // uneff: vect k=h*q;
          // uneff: vect l=sh*g;

          // generate next point
          for(counter i=0;i<dims;i++)
                {
                  oldcurrent[i]+=h*q[i]+sh*g[i]*u[i];
                }

          // uneff:               l[i]*=u[i];

          // uneff: oldcurrent+=k;
          // uneff: oldcurrent+=l;
          
          // vect qavg=0.5* h * ( q + (*deter)(oldcurrent) );
          // vect gavg=0.5* sh * ( g + (*stoch)(oldcurrent) );

          vect qn=(*deter)(oldcurrent);
          vect gn=(*stoch)(oldcurrent);

          for(counter i=0;i<dims;i++) 
                {
                  current[i]+=0.5*( h*( q[i] + qn[i]) + sh*u[i]*(g[i]+gn[i]));
                  
                  // uneff:  gavg[i]*=u[i];
                }

          // uneff: current += qavg;
          // uneff: current += gavg;
          
          return current;
        }

  private:              
        Normal rnd; // will be seeded with time()

        vf* trans;
  };

} // end namespace
#endif

/*********************************************************************
$Id: integrator.h,v 1.1.2.1 2002/06/19 13:15:47 mpeeters Exp $
**********************************************************************

$Log: integrator.h,v $
Revision 1.1.2.1  2002/06/19 13:15:47  mpeeters
New integrator manipulators, Euler-Maryuama, transforms.

Revision 1.1  2001/05/22 10:54:55  mpeeters
Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter".

Revision 1.4  2001/05/21 11:53:16  mpeeters
Removed Makefile.in, which is automatically generated anyway.

Revision 1.3  2000/09/22 08:50:03  mpeeters
Changed dependencies causing trouble.

Revision 1.2  2000/09/15 10:26:31  mpeeters
Cleaned out header and added CVS tails to files, removed superfuous
@author comments, inserted dates

*********************************************************************/

---

To get the sources or tarballs, please go to SourceForge or you can use the CVS repository.

More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be

Last update: June 2002.

---

Looking for Open Source? Check out !