mixturepdf.hpp

dysii is a C++ library for distributed probabilistic inference and learning in large-scale dynamical

    j++;
  }
  for (; j < getSize(); j++) {
    Ws[j] = Ws[j - 1] + ws(j);
  }
  
  dirty();
}

template <class T>
inline const indii::ml::aux::vector&
    indii::ml::aux::MixturePdf<T>::getWeights() const {
  return ws;
}

template <class T>
void indii::ml::aux::MixturePdf<T>::setWeights(const vector& ws) {
  /* pre-condition */
  assert (this->ws.size() == ws.size());

  this->ws = ws;

  /* recalculate cumulative weights */
  unsigned int i = 0;
  if (getSize() > 0) {
    Ws[0] = ws(0);
    i++;
  }
  for (; i < getSize(); i++) {
    Ws[i] = Ws[i - 1] + ws(i);
  }
  
  dirty();
}

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

template <class T>
void indii::ml::aux::MixturePdf<T>::normalise() {
  if (!Ws.empty()) {
    const double W = getTotalWeight();
  
    if (W != 1.0) {
      double WI = 1.0 / W;
    
      /* update weights */
      ws *= WI;
    
      /* update cumulative weights */
      std::vector<double>::iterator iter, end;
      iter = Ws.begin();
      end = Ws.end();
      while (iter != end) {
        *iter *= WI;
        iter++;
      }
      
      dirty(); // unnormalised mean out of date
    }
  }
}

template <class T>
void indii::ml::aux::MixturePdf<T>::clear() {
  xs.clear();
  ws.resize(0, false);
  Ws.clear();

  dirty();
}

template <class T>
double indii::ml::aux::MixturePdf<T>::getTotalWeight() const {
  return (Ws.empty() ? 0.0 : Ws.back());
}

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

  mu.resize(N, preserve);
  Zmu.resize(N, preserve);

  if (preserve) {
    unsigned int i;
    for (i = 0; i < xs.size(); i++) {
      xs[i].setDimensions(N, preserve);
    }
  } else {
    clear();
  }

  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() {
  /* pre-condition */
  assert (!xs.empty());

  double u = Random::uniform(0.0, getTotalWeight());
  unsigned int i = std::distance(Ws.begin(), lower_bound(Ws.begin(),
      Ws.end(), u));
  
  return xs[i].sample();
}

template <class T>
double indii::ml::aux::MixturePdf<T>::densityAt(const vector& x) {
  double p = 0.0;
  if (getTotalWeight() > 0.0) {
    unsigned int i;
    for (i = 0; i < xs.size(); i++) {
      p += ws(i) * xs[i].densityAt(x);
    }
    p /= getTotalWeight();
  }

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

  return p;
}

template <class T>
unsigned int indii::ml::aux::MixturePdf<T>::getDistributedSize()
    const {
  boost::mpi::communicator world;

  return boost::mpi::all_reduce(world, getSize(), 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 > 0.0 && W_d != 1.0) {
    double WI_d = 1.0 / W_d;

    /* update weights */
    ws *= WI_d;
  
    /* update cumulative weights */
    std::vector<double>::iterator iter, end;
    iter = Ws.begin();
    end = Ws.end();
    while (iter != end) {
      *iter *= WI_d;
      iter++;
    }

    dirty(); // unnormalised mean out of date
  }
}

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

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

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

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

  if (u >= W_s - getTotalWeight() && 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, getTotalWeight(), 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 - getTotalWeight() && 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>::distributedDensityAt(
    std::vector<vector>& xs) {
  boost::mpi::communicator world;  
  unsigned int size = world.size();

  unsigned int i, k;
  vector p(xs.size());

  p.clear();
  for (k = 0; k < size; k++) {
    for (i = 0; i < xs.size(); i++) {
      p(i) += getTotalWeight() * densityAt(xs[i]);
    }
    rotate(p);
    rotate(xs);
  }
  p /= getDistributedTotalWeight();

  return p;
}

template <class T>
indii::ml::aux::vector
    indii::ml::aux::MixturePdf<T>::getDistributedExpectation() {
  boost::mpi::communicator world;
  const unsigned int size = world.size();
  
  if (size == 0) {
    return getExpectation();
  } else {
    vector mu_d(N);

    if (getTotalWeight() > 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>::redistributeBySize() {
  namespace aux = indii::ml::aux;
  namespace ublas = boost::numeric::ublas;

  boost::mpi::communicator world;
  unsigned int rank = world.rank();
  unsigned int size = world.size();

  std::vector<T> xsBuffer;
  aux::vector wsBuffer;

  boost::mpi::request reqSendXs, reqSendWs;
  boost::mpi::request reqRecvXs, reqRecvWs;

  unsigned int P_total;
  unsigned int P_target;  // target no. components per node
  unsigned int P_change;  // no. components to transfer

  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 */
  P_total = getDistributedSize();
  P_target = P_total / size;
  if (rank < P_total % size) {
    P_target++; // take a leftover
  }

  /* gather excess components at each node */
  excess.rank = rank;
  excess.prop = getSize() - P_target;

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

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

    P_change = std::min(abs(to.prop), from.prop);

    if (rank == from.rank) {
      /* send components */
      xsBuffer.clear();
      xsBuffer.insert(xsBuffer.end(), xs.end() - P_change, xs.end());
      xs.resize(xs.size() - P_change);
      reqSendXs = world.isend(to.rank, 0, xsBuffer);
      
      /* send weights */
      wsBuffer.resize(P_change, false);
      noalias(wsBuffer) = project(ws, ublas::range(ws.size() - P_change,
          ws.size()));
      ws.resize(ws.size() - P_change, true);
      reqSendWs = world.isend(to.rank, 1, wsBuffer);
      
      /* update cumulative weights */
      Ws.resize(Ws.size() - P_change);
      
      /* wait */
      reqSendWs.wait();
      reqSendXs.wait();
    } else if (rank == to.rank) {
      /* receive components */
      reqRecvXs = world.irecv(from.rank, 0, xsBuffer);
      reqRecvWs = world.irecv(from.rank, 1, wsBuffer);
      reqRecvWs.wait();
      reqRecvXs.wait();
      
      assert (xsBuffer.size() == wsBuffer.size());
      
      for (i = 0; i < xsBuffer.size(); i++) {
        add(xsBuffer[i], wsBuffer(i));
      }
    }

    /* update excesses */
    from.prop -= P_change;
    to.prop += P_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();
  
  /* post-conditions */
  assert (xs.size() == ws.size());
  assert (xs.size() == Ws.size());
}

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

  std::vector<T> xsBuffer;
  aux::vector wsBuffer;

  boost::mpi::request reqSendXs, reqSendWs;
  boost::mpi::request reqRecvXs, reqRecvWs;

  double W_target = getDistributedTotalWeight() / size;
       // target weight per node
  double W_change;  // weight to transfer
  double W_actual;  // actual weight transferred

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

  unsigned int i;

  /* gather excess weight at each node */
  excess.rank = rank;
  excess.prop = getTotalWeight() - W_target;
  boost::mpi::all_gather(world, excess, excesses);
  std::stable_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().prop != excesses.back().prop) {
    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.prop), from.prop);

    if (rank == from.rank) {
      /* send components of up to W_change weight */
      W_actual = 0.0;
      xsBuffer.clear();
      wsBuffer.resize(0, false);
      
      while (W_actual + ws(ws.size() - 1) < W_change) {
        /* component */
        xsBuffer.push_back(xs.back());
        xs.pop_back();
        
        /* weight */
        wsBuffer.resize(wsBuffer.size() + 1, true);
        wsBuffer(wsBuffer.size() - 1) = ws(ws.size() - 1);
        ws.resize(ws.size() - 1, true);
        
        /* cumulative weights */
        Ws.pop_back();
        
        W_actual += wsBuffer(wsBuffer.size() - 1);
      }
      
      reqSendXs = world.isend(to.rank, 0, xsBuffer);
      reqSendWs = world.isend(to.rank, 1, wsBuffer);
      reqSendWs.wait();
      reqSendXs.wait();      

    } else if (rank == to.rank) {
      /* receive components */
      reqRecvXs = world.irecv(from.rank, 0, xsBuffer);
      reqRecvWs = world.irecv(from.rank, 1, wsBuffer);
      reqRecvWs.wait();
      reqRecvXs.wait();
      
      assert (xsBuffer.size() == wsBuffer.size());
      
      for (i = 0; i < xsBuffer.size(); i++) {
        add(xsBuffer[i], wsBuffer(i));
      }
    }

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

    /* update excesses */
    from.prop -= W_actual;
    to.prop += 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();

  /* post-conditions */
  assert (xs.size() == ws.size());
  assert (xs.size() == Ws.size());
}

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

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

  unsigned int i;

  Zmu.clear();
  for (i = 0; i < xs.size(); i++) {
    noalias(Zmu) += ws(i) * xs[i].getExpectation();
  }
  noalias(mu) = Zmu / getTotalWeight();
  haveMu = true;
}

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 & xs;
  ar & ws;
  ar & Ws;
}

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 & xs;
  ar & ws;
  ar & Ws;

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

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

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 & prop;
}

#endif