stochasticdifferentialmodel.hpp

dysii是一款非常出色的滤波函数库

#ifndef INDII_ML_ODE_STOCHASTICDIFFERENTIALMODEL_HPP
#define INDII_ML_ODE_STOCHASTICDIFFERENTIALMODEL_HPP

#include "../aux/vector.hpp"
#include "../aux/matrix.hpp"

#include "boost/serialization/serialization.hpp"

#include <vector>

namespace indii {
  namespace ml {
    namespace sde {

/**
 * StochasticAdaptiveRungeKutta compatible model.
 *
 * @author Lawrence Murray <[e-mail]>
 * @version $Rev: 405 $
 * @date $Date: 2008-03-05 15:05:55 +0000 (Wed, 05 Mar 2008) $
 *
 * The model is of the form:
 *
 * \f[
 * d\mathbf{y} = \mathbf{a}(\mathbf{y},t)\,dt + B(\mathbf{y},t)\,d\mathbf{W},
 * \f]
 *
 * where \f$\mathbf{a}(\mathbf{y},t)\f$ is the drift term and
 * \f$B(\mathbf{y},t)\f$ the diffusion term. At any time \f$t\f$, the
 * time derivatives may be calculated as:
 *
 * \f[
 *   \frac{dy_k}{dt} = a_k(\mathbf{y},t) -
 *   \frac{1}{2}\sum_{ij}B_{ij}(\mathbf{y},t)\frac{\partial
 *   B_{kj}(\mathbf{y},t)}{\partial y_i} + \frac{1}{\Delta
 *   t}\sum_{i}B_{ki}(\mathbf{y},t)\Delta W_i,
 * \f]
 *
 * or in matrix form:
 *
 * \f[
 *   \frac{d\mathbf{y}}{dt} = \mathbf{a}(\mathbf{y},t) -
 *   \frac{1}{2}\sum_i \frac{\partial B(\mathbf{y},t)}{\partial
 *   y_i}B_{i,*}(\mathbf{y},t)^T + \frac{1}{\Delta t}
 *   B(\mathbf{y},t)\Delta \mathbf{W},
 * \f]
 *
 * where \f$\Delta t\f$ is the time step, \f$\Delta \mathbf{W}\f$ the
 * Wiener increment and \f$B_{i,*}(\mathbf{y},t)\f$ the \f$i\f$th row of
 * \f$B(\mathbf{y},t)\f$.
 *
 * This is all calculated by StochasticAdaptiveRungeKutta, although
 * the model must provide \f$\mathbf{a}(\mathbf{y},t)\f$,
 * \f$B(\mathbf{y},t)\f$ and optionally \f$\frac{\partial
 * B_{kj}(\mathbf{y},t)}{\partial y_i}\f$ through implementations of
 * the calculateDrift(), calculateDiffusion() and
 * calculateDiffusionDerivatives() functions, respectively.
 */
class StochasticDifferentialModel {
public:
  /**
   * Default constructor for restoring from serialization.
   */
  StochasticDifferentialModel();

  /**
   * Constructor.
   *
   * @param dimensions Number of state variables in the system.
   *
   * The Wiener process noise is assumed to have the same
   * dimensionality as the state.
   */
  StochasticDifferentialModel(const unsigned int dimensions);

  /**
   * Constructor.
   *
   * @param dimensions Number of state variables in the system.
   * @param noiseDimensions Number of noise variables in the system.
   */
  StochasticDifferentialModel(const unsigned int dimensions,
      const unsigned int noiseDimensions);

  /**
   * Destructor.
   */
  virtual ~StochasticDifferentialModel();

  /**
   * Number of state variables in the system.
   *
   * @return the number of state variables in the system.
   */
  unsigned int getDimensions();

  /**
   * Number of noise variables in the system.
   *
   * @return the number of noise variables in the system.
   */
  unsigned int getNoiseDimensions();

  /**
   * Calculate drift term of the system at a given time.
   *
   * @param ts \f$t_s\f$; the proposed new time.
   * @param y \f$\mathbf{y}(t)\f$; the values of all state variables
   * at time \f$t\f$.
   *
   * @return \f$\mathbf{a}(\mathbf{y}, t)\f$; the drift term at time
   * \f$t\f$.
   */
  virtual indii::ml::aux::vector calculateDrift(double ts,
      const indii::ml::aux::vector& y) = 0;

  /**
   * Calculate diffusion term of the system at a given time.
   *
   * @param ts \f$t_s\f$; the proposed new time.
   * @param y \f$\mathbf{y}(t)\f$; the values of all state variables
   * at time \f$t\f$.
   *
   * @return \f$B(\mathbf{y}, t)\f$; the diffusion term at time
   * \f$t\f$.
   */
  virtual indii::ml::aux::matrix calculateDiffusion(double ts,
      const indii::ml::aux::vector& y) = 0;

  /**
   * Calculate the partial derivatives of the diffusion term with
   * respect to each state variable at a given time.
   *
   * @param ts \f$t_s\f$; the proposed new time.
   * @param y \f$\mathbf{y}(t)\f$; the values of all state variables
   * at time \f$t\f$.
   *
   * @return A vector of \f$N\f$ matrices, where matrix \f$i\f$ is
   * \f$\frac{\partial B(\mathbf{y},t)}{\partial y_i}\f$ at time
   * \f$t\f$. In the special case that \f$B(\mathbf{y}, t)\f$ is not
   * dependent on \f$\mathbf{y}\f$, may return an empty
   * std::vector<aux::matrix>. This is the behaviour of the default
   * implementation and has the same effect as returning a vector of
   * \f$N\f$ zero matrices.
   */
  virtual std::vector<indii::ml::aux::matrix> calculateDiffusionDerivatives(
      double ts, const indii::ml::aux::vector& y);

private:
  /**
   * Number of state variables in the system.
   */
  unsigned int dimensions;

  /**
   * Number of noise variables in the system.
   */
  unsigned int noiseDimensions;

  /**
   * Serialize, or restore from serialization.
   */
  template<class Archive>
  void serialize(Archive& ar, const unsigned int version);

  /*
   * Boost.Serialization requirements.
   */
  friend class boost::serialization::access;
  
};

    }
  }
}

template<class Archive>
void indii::ml::sde::StochasticDifferentialModel::serialize(Archive& ar,
    const unsigned int version) {
  ar & dimensions;
  ar & noiseDimensions;
}

#endif