wienerprocess.hpp

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

}

template <class T>
WienerProcess<T>::WienerProcess(const vector& mu,
    const symmetric_matrix& sigma) : N(mu.size()), mu(N), sigma(N), L(N,N),
    sigmaI(N), sigmaIDiag(N), haveCholesky(false), haveInverse(false),
    haveDeterminant(false), haveZ(false), haveDensityOptimisations(false) {
  /* pre-condition */
  assert(mu.size() == sigma.size1());

  setDrift(mu);
  setDiffusion(sigma);
}

template <class T>
WienerProcess<T>::WienerProcess(unsigned int N) : N(N), mu(N), sigma(N),
    L(N,N), sigmaI(N), sigmaIDiag(N), haveCholesky(false), haveInverse(false),
    haveDeterminant(false), haveZ(false), haveDensityOptimisations(false) {
  setDrift(aux::zero_vector(N));
  setDiffusion(aux::identity_matrix(N));
}

template <class T>
WienerProcess<T>& WienerProcess<T>::operator=(const WienerProcess<T>& o) {
  /* pre-condition */
  assert (o.N == N);

  haveCholesky = o.haveCholesky;
  haveInverse = o.haveInverse;
  haveDeterminant = o.haveDeterminant;
  haveZ = o.haveZ;
  haveDensityOptimisations = o.haveDensityOptimisations;

  noalias(mu) = o.mu;
  noalias(sigma) = o.sigma;
  if (haveCholesky) {
    noalias(L) = o.L;
  }
  if (haveInverse) {
    noalias(sigmaI) = o.sigmaI;
  }
  if (haveDeterminant) {
    sigmaDet = o.sigmaDet;
  }
  if (haveZ) {
    ZI = o.ZI;
  }
  if (haveDensityOptimisations) {
    isMuZero = o.isMuZero;
    isSigmaIDiagonal = o.isSigmaIDiagonal;
    if (isSigmaIDiagonal) {
      sigmaIDiag = o.sigmaIDiag;
    }
  }  

  return *this;
}

template <class T>
WienerProcess<T>::~WienerProcess() {
  //
}

template <class T>
unsigned int WienerProcess<T>::getDimensions() const {
  return N;
}

template <class T>
void WienerProcess<T>::setDimensions(const unsigned int N,
    const bool preserve) {
  this->N = N;

  mu.resize(N, preserve);
  sigma.resize(N, preserve);
  L.resize(N, N, preserve);
  sigmaI.resize(N, preserve);
  sigmaIDiag.resize(N, preserve);

  dirty();
}

template <class T>
vector WienerProcess<T>::getExpectation(const T delta) const {
  vector result(mu * delta);
  return result;
}  

template <class T>
symmetric_matrix WienerProcess<T>::getCovariance(const T delta) const {
  symmetric_matrix result(sigma * delta);
  return result;
}

template <class T>
vector WienerProcess<T>::getExpectation(const T delta) {
  vector result(mu * delta);
  return result;
}  

template <class T>
symmetric_matrix WienerProcess<T>::getCovariance(const T delta) {
  symmetric_matrix result(sigma * delta);
  return result;
}

template <class T>
const vector& WienerProcess<T>::getDrift() {
  return mu;
}  

template <class T>
const symmetric_matrix& WienerProcess<T>::getDiffusion() {
  return sigma;
}

template <class T>
const vector& WienerProcess<T>::getDrift() const {
  return mu;
}  

template <class T>
const symmetric_matrix& WienerProcess<T>::getDiffusion() const {
  return sigma;
}

template <class T>
void WienerProcess<T>::setDrift(const vector& mu) {
  /* pre-condition */
  assert(mu.size() == N); // new same size as old

  this->mu = mu;

  dirty();
}

template <class T>
void WienerProcess<T>::setDiffusion(const symmetric_matrix& sigma) {
  /* pre-condition */
  assert(sigma.size1() == N); // new same size as old

  this->sigma = sigma;

  dirty();
}

template <class T>
vector WienerProcess<T>::sample(const T delta) {
  vector z(N);
  unsigned int i;

  if (!haveCholesky) {
    calculateCholesky();
    assert(haveCholesky);
  }

  /* \f$z\f$ sampled from \f$\mathcal{N}(\mathbf{0}_N,\sqrt{\Delta t}I_N)\f$ */
  for (i = 0; i < N; i++) {
    z(i) = Random::gaussian(0.0, sqrt(delta));
  }

  /* \f$\mathbf{x} = \Delta t \mathbf{\mu} + L\mathbf{z})\f$ */
  return delta*mu + prod(L,z);
}

template <class T>
double WienerProcess<T>::densityAt(const T delta, const vector& x) {
  if (!haveDensityOptimisations) {
    calculateDensityOptimisations();
    assert(haveDensityOptimisations);
  }

  double deltaI = 1 / delta;
  double exponent;
  if (isMuZero && isSigmaIDiagonal) {
    exponent = inner_prod(x, element_prod(sigmaIDiag, x));
  } else {
    if (isMuZero) {
      exponent = inner_prod(x, prod(sigmaI, x));
    } else {
      aux::vector d(x - mu);
      if (isSigmaIDiagonal) {
	exponent = inner_prod(d, element_prod(sigmaIDiag, d));
      } else {
	exponent = inner_prod(d, prod(sigmaI, d));
      }
    }
  }

  double z = sqrt(std::pow(deltaI, static_cast<int>(N))) * ZI;
  return z * exp(-0.5 * deltaI * exponent);
}

template <class T>
void WienerProcess<T>::dirty() {
  haveCholesky = false;
  haveInverse = false;
  haveDeterminant = false;
  haveZ = false;
  haveDensityOptimisations = false;
}

template <class T>
double WienerProcess<T>::calculateDensity(const T delta, const vector& x) {
  return densityAt(delta, x);
}

template <class T>
void WienerProcess<T>::calculateCholesky() {
  matrix tmp(sigma);
  int err;

  err = lapack::potrf('L', tmp);
  assert (err == 0);
  noalias(L) = ublas::triangular_adaptor<matrix, ublas::lower>(tmp);
  haveCholesky = true;
  
  /* post-condition */
  assert(haveCholesky);
}

template <class T>
void WienerProcess<T>::calculateInverse() {
  if (!haveCholesky) {
    calculateCholesky();
    assert(haveCholesky);
  }

  matrix X(L), I(N,N);
  I = identity_matrix(N);

  int err;
  err = lapack::potrs('L', X, I);
  assert (err == 0);
  noalias(sigmaI) = ublas::symmetric_adaptor<matrix, ublas::lower>(I);
  haveInverse = true;

  /* post-condition */
  assert(haveInverse);
}

template <class T>
void WienerProcess<T>::calculateDeterminant() {
  if (!haveCholesky) {
    calculateCholesky();
    assert(haveCholesky);
  }

  ublas::matrix_vector_range<lower_triangular_matrix> d = diag(L);
  assert(d.size() > 0);
  unsigned int i;
  double LDet = d(0);
  for (i = 1; i < d.size(); i++) {
    LDet *= d(i);
  }
  sigmaDet = LDet*LDet;
  haveDeterminant = true;

  /* post-condition */
  assert(haveDeterminant);
}

template <class T>
void WienerProcess<T>::calculateZ() {
  if (!haveDeterminant) {
    calculateDeterminant();
    assert(haveDeterminant);
  }
  ZI = 1.0 / sqrt(std::pow(2.0*M_PI, static_cast<int>(N)) * sigmaDet);
  haveZ = true;

  /* post-condition */
  assert(haveZ);
}

template <class T>
void WienerProcess<T>::calculateDensityOptimisations() {
  unsigned int i, j;

  if (!haveDeterminant) {
    calculateDeterminant();
    assert(haveDeterminant);
  }
  if (!haveInverse) {
    calculateInverse();
    assert(haveInverse);
  }
  if (!haveZ) {
    calculateZ();
    assert(haveZ);
  }

  /* is expectation zero? */
  isMuZero = true;
  for (i = 0; i < N && isMuZero; i++) {
    isMuZero = isMuZero && mu(i) == 0.0;
  }

  /* is covariance diagonal? */
  isSigmaIDiagonal = true;
  for (j = 0; j < N && isSigmaIDiagonal; j++) {
    for (i = 0; i < j && isSigmaIDiagonal; i++) {
      isSigmaIDiagonal = isSigmaIDiagonal && sigma(i,j) == 0.0;
      /* ^ better to check sigma rather than sigmaI elements here --
	 probably supplied by user, making floating point equality
	 comparisons reliable. */
    }
  }
  if (isSigmaIDiagonal) {
    noalias(sigmaIDiag) = diag(sigmaI);
  }

  haveDensityOptimisations = true;

  /* post-condition */
  assert (haveDensityOptimisations == true);
}

template <class T>
template <class Archive>
void indii::ml::aux::WienerProcess<T>::save(Archive& ar,
    const unsigned int version) const {
  ar & boost::serialization::base_object<StochasticProcess<T> >(*this);
  ar & N;
  ar & mu;

  /* serialization of symmetric_matrix not yet supported in ublas */
  const matrix tmp(sigma);
  ar & tmp;
}

template <class T>
template <class Archive>
void indii::ml::aux::WienerProcess<T>::load(Archive& ar,
    const unsigned int version) {
  ar & boost::serialization::base_object<StochasticProcess<T> >(*this);
  ar & N;
  ar & mu;

  sigma.resize(N, false);
  L.resize(N,N, false);
  sigmaI.resize(N, false);
  sigmaIDiag.resize(N, false);

  /* serialization of symmetric_matrix not yet supported in ublas */
  matrix tmp(N,N);
  ar & tmp;
  noalias(sigma) = boost::numeric::ublas::symmetric_adaptor<matrix>(tmp);

  haveCholesky = false;
  haveInverse = false;
  haveDeterminant = false;
  haveZ = false;
  haveDensityOptimisations = false;
}

#endif