mixturepdf.hpp

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

#ifndef INDII_ML_AUX_MIXTUREPDF_HPP
#define INDII_ML_AUX_MIXTUREPDF_HPP

#include "Pdf.hpp"

#include <vector>

namespace indii {
  namespace ml {
    namespace aux {
    
/**
 * Mixture probability density.
 *
 * @author Lawrence Murray <[e-mail]>
 * @version $Rev: 406 $
 * @date $Date: 2008-03-05 16:33:45 +0000 (Wed, 05 Mar 2008) $
 *
 * @param T Component type, should be derivative of Pdf.
 *
 * @section MixturePdf_serialization Serialization
 *
 * This class supports serialization through the Boost.Serialization library.
 *
 * @section MixturePdf_parallel Parallelisation
 *
 * This class supports distributed storage. MixturePdf objects may be
 * instantiated on all nodes of a communicator, each with a different
 * set of mixture components, such that the union of all components
 * defines the full distribution.
 *
 * "Regular" methods such as getNumComponents(), getTotalWeight() and
 * getExpectation() use only those components stored on the local
 * node. Special "distributed" methods such as
 * getDistributedNumComponents(), getDistributedTotalWeight() and
 * getDistributedExpectation() use all components stored across the
 * nodes.
 */
template <class T>
class MixturePdf : public Pdf {
public:
  /**
   * Weighted component.
   */
  struct weighted_component {
  public:
    /**
     * Weight.
     */
    double w;

    /**
     * Component
     */
    T x;

    /**
     * Comparison operator for sorting.
     *
     * @return True if this component has greater weight than the
     * argument's component.
     */
    bool operator<(const weighted_component& o) const;
    
  private:
    /**
     * Serialize, or restore from serialization.
     */
    template<class Archive>
    void serialize(Archive& ar, const unsigned int version);

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

  };

  /**
   * Weighted component vector.
   */
  typedef std::vector<weighted_component> weighted_component_vector;

  /**
   * Weighted component iterator.
   */
  typedef typename weighted_component_vector::iterator
      weighted_component_iterator;

  /**
   * Weighted component const iterator.
   */
  typedef typename weighted_component_vector::const_iterator
      weighted_component_const_iterator;

  /**
   * Default constructor.
   *
   * Initialises the mixture with zero dimensions. This should
   * generally only be used when the object is to be restored from a
   * serialization. Indeed, there is no other way to resize the
   * mixture to nonzero dimensionality except by subsequently
   * restoring from a serialization.
   */
  MixturePdf();

  /**
   * Constructor. One or more components should be added with
   * addComponent() after construction.
   *
   * @param N Dimensionality of the distribution.
   */
  MixturePdf(const unsigned int N);

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

  /**
   * Assignment operator. Both sides must have the same dimensionality.
   */
  virtual MixturePdf<T>& operator=(const MixturePdf<T>& o);

  virtual unsigned int getDimensions() const;

  virtual void setDimensions(const unsigned int N,
      const bool preserve = false);

  /**
   * @name Local storage methods
   *
   * Use these methods if:
   *
   * @li You are not working in a parallel environment,
   * @li You are working in a parallel environment but are not using
   * the distributed storage features of this class, or
   * @li You are working in a parallel environment and are using the
   * distributed storage features of this class, but only want to deal
   * with the mixture components stored on the local node.
   */
  //@{

  /**
   * Add a component to the distribution on the local node.
   *
   * @param x The component.
   * @param w Unnormalised weight of the component.
   *
   * If the weight is not-a-number or zero, the sample is ignored for
   * reasons of efficiency, as it will have no effect.
   */
  void addComponent(const T& x, const double w = 1.0);

  /**
   * Add a component to the distribution on the local node.
   *
   * @param xw The weighted component.
   *
   * If the weight is not-a-number or zero, the component is ignored for
   * reasons of efficiency, as it will have no effect.
   */
  virtual void addComponent(const weighted_component& xw);
  
  /**
   * Get weighted components on the local node.
   *
   * @return \f$\{(\mathbf{x}^{(i)},w^{(i)})\}\f$; set of weighted
   * components defining the distribution, as a vector.
   */
  virtual const weighted_component_vector& getComponents() const;

  /**
   * Get weighted components on the local node.
   *
   * @return \f$\{(\mathbf{x}^{(i)},w^{(i)})\}\f$; set of weighted
   * components defining the distribution, as a vector.
   *
   * @note Modifying any of the components in the mixture may have
   * unintended side effects, especially since Pdf classes rely heavily on
   * precalculations which may consequently become out of date. This method
   * is provided only for those situations where precalculations must be
   * made within the component objects themselves, such as within their
   * getExpectation() methods, such that a non-const context is required.
   * For all other situations, favour the const version of this method.
   */
  virtual weighted_component_vector& getComponents();

  /**
   * Sort all components in descending order by weight. This is
   * particularly worthwhile to improve performance before repetitive
   * calls to sample().
   */
  void sort();

  /**
   * Clear all components from the distribution on the local node.
   */
  void clear();

  /**
   * Get the number of components in the distribution on the local
   * node.
   *
   * @return \f$K\f$; the number of components in the distribution.
   */
  unsigned int getNumComponents() const;

  /**
   * Get the total weight of components on the local node.
   *
   * @return \f$W\f$; the total weight of components.
   */
  double getTotalWeight() const;

  /**
   * Normalise weights on the local node to sum to 1.
   *
   * @warning If in a distributed storage situation this is probably
   * not what you want. Consider distributedNormalise()
   * instead. Normalising the weights only on the local node will skew
   * the weighting of mixture components across all nodes.
   */
  void normalise();

  /**
   * Sample from the distribution on the local node.
   *
   * @return A sample from the distribution.
   */
  virtual vector sample();

  /**
   * Calculate the density on the local node at a given point.
   *
   * \f[p(\mathbf{x}) = \frac{1}{W}\sum_{i=1}^{K}w^{(i)}p^{(i)}(\mathbf{x})\f]
   *
   * @param x \f$\mathbf{x}\f$; the point at which to calculate the
   * density.
   *
   * @return \f$p(\mathbf{x})\f$; the density at \f$\mathbf{x}\f$.
   */
  virtual double densityAt(const vector& x);

  /**
   * Get the expected value of the distribution on the local node.
   *
   * @return \f$\mathbf{\mu}\f$; expected value of the distribution.
   */
  virtual const vector& getExpectation();

  //@}

  /**
   * @name Distributed storage methods
   *
   * Use these methods if:
   *
   * @li You are working in a parallel environment and are using the
   * distributed storage features of this class, and want to deal with
   * the mixture components stored across all nodes.
   *
   * These methods are used for distributed storage of mixtures. They
   * require synchronization and communication between all nodes in
   * the communicator, such that if any of these methods is called on
   * one node, the same method should eventually, and preferably as
   * soon as possible, be called on all other nodes in the same
   * communicator.
   */
  //@{

  /**
   * Get the number of components in the distribution across all
   * nodes.
   *
   * @return \f$\sum_i K_i\f$; the number of components in the distribution.
   */
  unsigned int getDistributedNumComponents() const;

  /**
   * Get the total weight of components across all nodes.
   *
   * @return \f$\sum_i W_i\f$; the total weight of components across
   * all nodes.
   */
  double getDistributedTotalWeight() const;

  /**
   * Normalise weights across all nodes to sum to 1.
   */
  void distributedNormalise();

  /**
   * Sample from the full distribution.
   *
   * @return A sample from the full distribution.
   */
  virtual vector distributedSample();

  /**
   * Draw multiple samples from the full distribution. This is
   * significantly more efficient than multiple calls to
   * distributedSample(), as it involves less message passing.
   *
   * @param P \f$P\f$; number of samples to draw.
   *
   * @return Vector of samples drawn on the local node. The number of
   * samples drawn across all nodes will total \f$P\f$.
   */
  virtual std::vector<vector> distributedSample(const unsigned int P);

  /**
   * Calculate the density of the full distribution at a given point.
   *
   * @param x \f$\mathbf{x}\f$; the point at which to calculate the
   * density.
   *
   * @return \f$p(\mathbf{x})\f$; the density at \f$\mathbf{x}\f$.
   */
  virtual double distributedDensityAt(const vector& x);

  /**
   * Get the expected value of the full distribution.
   *
   * @return \f$\mathbf{\mu}\f$; expected value of the full
   * distribution.
   */
  virtual vector getDistributedExpectation();

  /**
   * Redistribute components across nodes by number. This attempts to
   * redistribute the components of the full distribution across
   * nodes, so that each node stores as close to an equal number of
   * components as possible.
   */
  void redistributeByCount();

  /**
   * Redistribute components across nodes by weight. This attempts to
   * redistribute the components of the full distribution across
   * nodes, so that the total weight of components at each node is as
   * close to an equal number as possible.
   */
  void redistributeByWeight();

  //@}

protected:
  /**
   * Called when changes are made to the distribution, such as when a
   * new component is added. This allows pre-calculations to be
   * refreshed.
   */
  virtual void dirty();

private:
  /**
   * Node weight.
   */
  template <class P>
  struct node_property {
    /**
     * Node rank.
     */
    int rank;

    /**
     * Property at node.
     */
    P p;

    /**
     * Comparison operator for sorting.
     *
     * @return True if this node's property is greater than the
     * argument's property.
     */
    bool operator<(const node_property& o) const;

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

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

  };

  typedef struct node_property<int> node_count;
  typedef std::vector<node_count> node_count_vector;
  typedef typename node_count_vector::iterator node_count_iterator;

  typedef struct node_property<double> node_weight;
  typedef std::vector<node_weight> node_weight_vector;
  typedef typename node_weight_vector::iterator node_weight_iterator;

  /**
   * Dimensionality of the distribution.
   */
  unsigned int N;

  /**
   * \f$W\f$; total weight.
   */
  double W;

  /**
   * Weighted components.
   */
  weighted_component_vector xs;

  /**
   * Last calculated expectation.
   */
  vector mu;

  /**
   * Last calculated unnormalised expectation.
   */
  vector Zmu;

  /**
   * Is the last calculated expectation up to date?
   */
  bool haveMu;

  /**
   * Is part of the component list sorted?
   */
  bool havePartialSort;

  /**
   * Points to the end of the sorted portion of the component list.
   */
  weighted_component_iterator sorted;

  /**
   * Calculate expectation from current components on the local node.
   */
  void calculateExpectation();

  /**
   * Given a number \f$u \in [0,1]\f$, returns the component \f$x_i\f$
   * on the local node such that \f$\sum_{j=1}^{i-1} w_j < Wu \le
   * \sum_{j=1}^{i} w_j\f$.
   *
   * @param u \f$u\f$
   *
   * @return \f$x_i\f$
   */
  T& findComponent(const double u);

  /**
   * Serialize.
   */
  template<class Archive>
  void save(Archive& ar, const unsigned int version) const;

  /**
   * Restore from serialization.
   */
  template<class Archive>
  void load(Archive& ar, const unsigned int version);

  /*
   * Boost.Serialization requirements.
   */
  BOOST_SERIALIZATION_SPLIT_MEMBER()
  friend class boost::serialization::access;

};

    }
  }
}

#include "Random.hpp"
#include "Parallel.hpp"

#include "boost/serialization/base_object.hpp"
#include "boost/serialization/vector.hpp"

#include "boost/mpi.hpp"
#include "boost/mpi/environment.hpp"
#include "boost/mpi/communicator.hpp"

#include <algorithm>

template <class T>
indii::ml::aux::MixturePdf<T>::MixturePdf() : N(0), mu(0), Zmu(0) {
  haveMu = false;
  havePartialSort = false;
  W = 0.0;
}

template <class T>
indii::ml::aux::MixturePdf<T>::MixturePdf(const unsigned int N) : N(N), mu(N),
    Zmu(N) {
  haveMu = false;
  havePartialSort = false;
  W = 0.0;
}

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

  xs = o.xs;
  W = o.W;
  mu = o.mu;
  Zmu = o.Zmu;
  haveMu = o.haveMu;
  havePartialSort = o.havePartialSort;
  sorted = o.sorted;

  return *this;
}

template <class T>
indii::ml::aux::MixturePdf<T>::~MixturePdf() {
  //
}

template <class T>
void indii::ml::aux::MixturePdf<T>::addComponent(const T& x, const double w) {
  weighted_component s = { w, x };
  addComponent(s);
}

template <class T>
void indii::ml::aux::MixturePdf<T>::addComponent(const weighted_component& xw){