mixturepdf.hpp

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

  /* pre-condition */
  assert (xw.x.getDimensions() == N);

  if (xw.w > 0.0) {
    xs.push_back(xw);
    W += xw.w;
    
    /**
     * @todo Store cumulative weights for efficient sampling using binary
     * search.
     */

    dirty();
  }
}

template <class T>
unsigned int indii::ml::aux::MixturePdf<T>::getNumComponents() const {
  return xs.size();
}

template <class T>
void indii::ml::aux::MixturePdf<T>::normalise() {
  if (W != 1.0) {
    double WI = 1.0 / W;
    weighted_component_iterator iter, end;
    iter = xs.begin();
    end = xs.end();
    while (iter != end) {
      iter->w *= WI;
      iter++;
    }
    W = 1.0;
  }
}

template <class T>
void indii::ml::aux::MixturePdf<T>::sort() {
  if (havePartialSort) {
    if (sorted != xs.end()) { // whole list not sorted
      std::inplace_merge(xs.begin(), sorted, xs.end());
      sorted = xs.end();
    }
  } else {
    std::sort(xs.begin(), xs.end());
    sorted = xs.end();
    havePartialSort = true;
  }
}

template <class T>
void indii::ml::aux::MixturePdf<T>::clear() {
  xs.clear();
  W = 0.0;

  dirty();
}

template <class T>
const typename indii::ml::aux::MixturePdf<T>::weighted_component_vector&
    indii::ml::aux::MixturePdf<T>::getComponents() const {
  return xs;
}

template <class T>
typename indii::ml::aux::MixturePdf<T>::weighted_component_vector&
    indii::ml::aux::MixturePdf<T>::getComponents() {
  return xs;
}

template <class T>
double indii::ml::aux::MixturePdf<T>::getTotalWeight() const {
  return W;
}

template <class T>
unsigned int indii::ml::aux::MixturePdf<T>::getDimensions() const {
  return N;
}

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

  mu.resize(getDimensions());
  Zmu.resize(getDimensions());

  if (preserve) {
    weighted_component_iterator iter, end;
    iter = xs.begin();
    end = xs.end();
    while (iter != end) {
      iter->x.setDimensions(N, preserve);
      iter++;
    }
  } else {
    xs.empty();
  }

  dirty();
}

template <class T>
const indii::ml::aux::vector& indii::ml::aux::MixturePdf<T>::getExpectation() {
  if (!haveMu) {
    calculateExpectation();
  }
  return mu;
}

template <class T>
indii::ml::aux::vector indii::ml::aux::MixturePdf<T>::sample() {
  return findComponent(Random::uniform(0.0, 1.0)).sample();
}

template <class T>
double indii::ml::aux::MixturePdf<T>::densityAt(const vector& x) {
  double p = 0.0;
  weighted_component_iterator iter, end;
  iter = xs.begin();
  end = xs.end();
  while (iter != end) {
    p += iter->w * iter->x.densityAt(x);
    iter++;
  }
  p /= W;

  /* post-condition */
  assert (p >= 0.0);

  return p;
}

template <class T>
unsigned int indii::ml::aux::MixturePdf<T>::getDistributedNumComponents()
    const {
  boost::mpi::communicator world;
  return boost::mpi::all_reduce(world, getNumComponents(),
      std::plus<unsigned int>());
}

template <class T>
double indii::ml::aux::MixturePdf<T>::getDistributedTotalWeight() const {
  boost::mpi::communicator world;
  return boost::mpi::all_reduce(world, getTotalWeight(),
      std::plus<double>()); 
}

template <class T>
void indii::ml::aux::MixturePdf<T>::distributedNormalise() {
  double W_d = getDistributedTotalWeight();

  if (W_d != 1.0) {
    double WI_d = 1.0 / W_d;
    weighted_component_iterator iter, end;
    iter = xs.begin();
    end = xs.end();
    while (iter != end) {
      iter->w *= WI_d;
      iter++;
    }
    W *= WI_d;
    dirty();
    /**
     * @todo This need only dirty unnormalised precalcs.
     */
  }
}

template <class T>
indii::ml::aux::vector indii::ml::aux::MixturePdf<T>::distributedSample() {
  boost::mpi::communicator world;
  int rank = world.rank();
  int size = world.size();
  int i;

  double W_s;
  double u;
  vector x(N);

  /* scan sum weights across nodes */
  W_s = boost::mpi::scan(world, W, std::plus<double>());

  /* final node has total weight, so it can generate random number up
     to total weight and broadcast to other nodes */
  if (rank == size - 1) {
    u = Random::uniform(0.0, W_s);
  }
  boost::mpi::broadcast(world, u, size - 1);

  if (u >= W_s - W && u < W_s) {
    /* this node draws sample and sends*/
    std::vector<boost::mpi::request> reqs;
    reqs.reserve(size - 1);
    x = sample();
    for (i = 0; i < size; i++) {
      if (i != rank) {
        reqs.push_back(world.isend(i, 0, x));
      }
    }

    // boost::mpi::wait_all seems exceptionally slow here, workaround:
    std::vector<boost::mpi::request>::iterator iter, end;
    iter = reqs.begin();
    end = reqs.end();
    while (iter != end) {
      iter->wait();
      iter++;
    }
    // end workaround

  } else {
    /* this node receives sample */
    world.recv(boost::mpi::any_source, boost::mpi::any_tag, x);
  }

  return x;
}

template <class T>
std::vector<indii::ml::aux::vector>
    indii::ml::aux::MixturePdf<T>::distributedSample(const unsigned int P) {
  boost::mpi::communicator world;
  int rank = world.rank();
  int size = world.size();
  unsigned int i;

  double W_s;
  double u;
  std::vector<double> us;
  std::vector<double>::iterator iter, end;
  std::vector<vector> xs;

  /* scan sum weights across nodes */
  W_s = boost::mpi::scan(world, W, std::plus<double>());

  /* final node has total weight, so it can generate random numbers up
     to total weight and broadcast to other nodes */
  if (rank == size - 1) {
    for (i = 0; i < P; i++) {
      us.push_back(Random::uniform(0.0, W_s));
    }
  }
  boost::mpi::broadcast(world, us, size - 1);

  iter = us.begin();
  end = us.end();
  while (iter != end) {
    u = *iter;
    if (u >= W_s - W && u < W_s) {
      /* this node draws sample */
      xs.push_back(sample());
    }
    iter++;
  }

  return xs;
}

template <class T>
double indii::ml::aux::MixturePdf<T>::distributedDensityAt(const vector& x) {
  boost::mpi::communicator world;  
  double p = boost::mpi::all_reduce(world, densityAt(x), std::plus<double>());
  
  /* post-condition */
  assert (p >= 0.0);

  return p;
}

template <class T>
indii::ml::aux::vector
    indii::ml::aux::MixturePdf<T>::getDistributedExpectation() {
  boost::mpi::communicator world;
  vector mu_d(N);

  if (W > 0.0) {
    if (!haveMu) {
      calculateExpectation();
    }
  } else {
    Zmu.clear();
  }

  noalias(mu_d) = boost::mpi::all_reduce(world, Zmu, std::plus<vector>());
  mu_d /= getDistributedTotalWeight();

  return mu_d;
}

template <class T>
void indii::ml::aux::MixturePdf<T>::redistributeByCount() {
  boost::mpi::communicator world;
  int rank = world.rank();
  int size = world.size();

  weighted_component_vector buffer;

  unsigned int K_total;
  unsigned int K_target;  // target no. components per node
  unsigned int K_change;

  node_count excess, from, to;
  node_count_vector excesses;  // no. excess components at each node

  unsigned int i;

  /* work out target number of components for this node */
  K_total = getDistributedNumComponents();
  K_target = K_total / size;
  if (rank < static_cast<int>(K_total % size)) {
    K_target++; // take a leftover
  }

  /* gather excess components at each node */
  excess.rank = rank;
  excess.p = xs.size() - K_target;

  boost::mpi::all_gather(world, excess, excesses);
  std::sort(excesses.begin(), excesses.end());

  /* while distribution of components is not even */
  while (excesses.front().p > 0) {
    from = excesses.front();  // node with largest excess
    to = excesses.back();  // node with smallest (largest -ve) excess
    excesses.erase(excesses.begin());
    excesses.pop_back();

    K_change = std::min(abs(to.p), from.p);  // no. components to transfer

    if (rank == from.rank) {
      /* send K_change no. components */
      buffer.reserve(K_change);
      for (i = 0; i < K_change; i++) {
        buffer.push_back(xs.back());
        W -= xs.back().w;
        xs.pop_back();
      }
      world.send(to.rank, 0, buffer);
      buffer.clear();
    } else if (rank == to.rank) {
      /* receive K_change no. components */
      world.recv(from.rank, 0, buffer);
      while (buffer.size() > 0) {
        addComponent(buffer.back());
        buffer.pop_back();
      }
    }

    /* update excesses */
    from.p -= K_change;
    to.p += K_change;
  
    /* re-sort */
    excesses.push_back(from); // 'from' must have >= no. components of 'to'
    excesses.push_back(to);
    inplace_merge(excesses.begin(), excesses.end() - 2, excesses.end());
  }

  dirty();
}

template <class T>
void indii::ml::aux::MixturePdf<T>::redistributeByWeight() {
  boost::mpi::communicator world;
  int rank = world.rank();
  int size = world.size();

  weighted_component_vector buffer;

  double W_target = getDistributedTotalWeight() / size;
      // target weight per node
  double W_change, W_actual;

  node_weight excess, from, to;
  node_weight_vector excesses; // excess weight at each node

  /* gather excess weight at each node */
  excess.rank = rank;
  excess.p = W - W_target;
  boost::mpi::all_gather(world, excess, excesses);
  std::sort(excesses.begin(), excesses.end());

  /* while distribution of components is not even */
  from = excesses.front();
  to = from;  // will never enter the loop if only one node
  while ((from.rank != excesses.front().rank ||
      to.rank != excesses.back().rank) &&
      excesses.front().p != excesses.back().p) {
    from = excesses.front();  // node with largest excess
    to = excesses.back();  // node with smallest (largest -ve) excess
    excesses.erase(excesses.begin());
    excesses.pop_back();

    W_change = std::min(fabs(to.p), from.p);  // maximum weight to transfer

    if (rank == from.rank) {
      /* send components of up to W_change weight */
      W_actual = 0.0;
      while (W_actual + xs.back().w < W_change) {
        buffer.push_back(xs.back());
        W -= xs.back().w;
        W_actual += xs.back().w;
        xs.pop_back();
      }
      world.send(to.rank, 0, buffer);
      buffer.clear();
    } else if (rank == to.rank) {
      /* receive components */
      world.recv(from.rank, 0, buffer);
      while (buffer.size() > 0) {
        addComponent(buffer.back());
        buffer.pop_back();
      }
    }

    /* broadcast actual weight transfer to all nodes */
    boost::mpi::broadcast(world, W_actual, from.rank);

    /* update excesses */
    from.p -= W_actual;
    to.p += W_actual;

    /* re-sort */
    excesses.push_back(from); // 'from' must have >= weight of 'to'
    excesses.push_back(to);
    inplace_merge(excesses.begin(), excesses.end() - 2, excesses.end());
  }

  dirty();
}

template <class T>
void indii::ml::aux::MixturePdf<T>::dirty() {
  haveMu = false;
  havePartialSort = false;
}

template <class T>
void indii::ml::aux::MixturePdf<T>::calculateExpectation() {
  /* pre-condition */
  assert (W > 0.0);

  weighted_component_iterator iter, end;

  Zmu.clear();
  iter = xs.begin();
  end = xs.end();
  while (iter != end) {
    noalias(Zmu) += iter->w * iter->x.getExpectation();
    iter++;
  }
  noalias(mu) = Zmu / W;
  haveMu = true;
}

template <class T>
T& indii::ml::aux::MixturePdf<T>::findComponent(const double u) {
  /* pre-condition */
  assert (u >= 0.0 && u <= 1.0);

  double w = W*u;
  weighted_component_iterator iter, end;

  iter = xs.begin();
  end = xs.end();
  while (iter != end) {
    w -= iter->w;
    if (w <= 0.0) {
      return iter->x;
    }
    iter++;
  }
  assert (false);
  return iter->x;
}

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

  ar & N;
  ar & W;
  ar & xs;
}

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

  ar & N;
  ar & W;
  ar & xs;

  mu.resize(N, false);
  Zmu.resize(N, false);
  haveMu = false;
  havePartialSort = false;
}

template <class T>
bool indii::ml::aux::MixturePdf<T>::weighted_component::operator<(
    const weighted_component& o) const {
  return w > o.w;
}

template <class T>
template <class Archive>
void indii::ml::aux::MixturePdf<T>::weighted_component::serialize(Archive& ar,
    const unsigned int version) {
  ar & w;
  ar & x;
}

template <class T>
template <class P>
bool indii::ml::aux::MixturePdf<T>::node_property<P>::operator<(
    const node_property& o) const {
  return p > o.p;
}

template <class T>
template <class P>
template <class Archive>
void indii::ml::aux::MixturePdf<T>::node_property<P>::serialize(Archive& ar,
    const unsigned int version) {
  ar & rank;
  ar & p;
}

#endif