diracmixturepdf.cpp

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

//#if defined(__GNUC__) && defined(GCC_PCH)
//  #include "aux.hpp"
//#else
  #include "DiracMixturePdf.hpp"
//#endif

using namespace indii::ml::aux;

namespace ublas = boost::numeric::ublas;

DiracMixturePdf::DiracMixturePdf() : MixturePdf<DiracPdf>(0), sigma(0),
    Zsigma(0) {
  haveSigma = false;
}

DiracMixturePdf::DiracMixturePdf(Pdf& o, const unsigned int P) :
    MixturePdf<DiracPdf>(o.getDimensions()), sigma(o.getDimensions()),
    Zsigma(o.getDimensions()), haveSigma(false) {
  unsigned int i;

  for (i = 0; i < P; i++) {
    addComponent(o.sample());
  }
}

DiracMixturePdf::DiracMixturePdf(const unsigned int N) :
  MixturePdf<DiracPdf>(N), sigma(N), Zsigma(N) {
  haveSigma = false;
}

DiracMixturePdf::~DiracMixturePdf() {
  //
}

void DiracMixturePdf::setDimensions(const unsigned int N,
    const bool preserve) {
  MixturePdf<DiracPdf>::setDimensions(N, preserve);

  Zsigma.resize(N, false);
  sigma.resize(N, false);
}

const symmetric_matrix& DiracMixturePdf::getCovariance() {
  if (!haveSigma) {
    calculateCovariance();
  }
  return sigma;
}

symmetric_matrix DiracMixturePdf::getDistributedCovariance() {
  const unsigned int N = getDimensions();

  boost::mpi::communicator world;
  matrix Zsigma_d(N,N), sigma_d(N,N);
  vector mu_d(N);

  if (getTotalWeight() > 0.0) {
    if (!haveSigma) {
      calculateCovariance();
    }
  } else {
    Zsigma.clear();
  }

  noalias(mu_d) = getDistributedExpectation();
  noalias(Zsigma_d) = boost::mpi::all_reduce(world, matrix(Zsigma),
      std::plus<matrix>());
  noalias(sigma_d) = Zsigma_d / getDistributedTotalWeight() -
      outer_prod(mu_d, mu_d);

  return ublas::symmetric_adaptor<matrix, ublas::lower>(sigma_d);
}

void DiracMixturePdf::dirty() {
  MixturePdf<DiracPdf>::dirty();
  haveSigma = false;
}

void DiracMixturePdf::calculateCovariance() {
  /* pre-condition */
  assert (getTotalWeight() > 0.0);

  weighted_component_const_iterator iter, end;
  const vector& mu = getExpectation();

  Zsigma.clear();
  iter = getComponents().begin();
  end = getComponents().end();
  while (iter != end) {
    noalias(Zsigma) += iter->w * outer_prod(iter->x, iter->x);
    iter++;
  }
  noalias(sigma) = Zsigma / getTotalWeight() - outer_prod(mu, mu);
  haveSigma = true;
}