#include <cstdlib>
#include <time.h>
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <sstream>

#include "BPSolver.h"
#include "BpNode.h"


BPSolver* Edge::klass = 0;
StatisticMap Statistics::stats_;
unsigned Statistics::numCreatedNets     = 0;
unsigned Statistics::numSolvedPolyTrees = 0;
unsigned Statistics::numSolvedLoopyNets = 0;
unsigned Statistics::numUnconvergedRuns = 0;
unsigned Statistics::maxIterations      = 0;
unsigned Statistics::totalOfIterations  = 0;


BPSolver::BPSolver (const BayesNet& bn) : Solver (&bn)
{
  bn_ = &bn;
  forceGenericSolver_ = false;
  //forceGenericSolver_ = true;
  schedule_           = S_SEQ_FIXED;
  //schedule_           = S_SEQ_RANDOM;
  //schedule_           = S_PARALLEL;
  //schedule_           = S_MAX_RESIDUAL;
  maxIter_            = 205;
  accuracy_           = 0.000001;
}



BPSolver::~BPSolver (void)
{
  for (unsigned i = 0; i < msgs_.size(); i++) {
    delete msgs_[i];
  }
}



void
BPSolver::runSolver (void)
{
  if (DL >= 1) {
    //bn_->printNetwork();
  }

  clock_t start_ = clock();
  if (bn_->isSingleConnected() && !forceGenericSolver_) {
    runPolyTreeSolver();
    Statistics::numSolvedPolyTrees ++;
  } else {
    runGenericSolver();
    Statistics::numSolvedLoopyNets ++;
    if (nIter_ >= maxIter_) {
      Statistics::numUnconvergedRuns ++;
    } else {
      Statistics::updateIterations (nIter_);
    }
    if (DL >= 1) {
      cout << endl;
      if (nIter_ < maxIter_) {
        cout << "Belief propagation converged in " ; 
        cout << nIter_ << " iterations" << endl;
      } else {
        cout << "The maximum number of iterations was hit, terminating..." ;
        cout << endl;
      }
    }
  }
  double time = (double (clock() - start_)) / CLOCKS_PER_SEC;
  unsigned size = bn_->getNumberOfNodes();
  Statistics::updateStats (size, time);
  //if (size > 30) {
  //   stringstream ss;
  //  ss << size << "." << Statistics::getCounting (size) << ".dot" ;
  //  bn_->exportToDotFile (ss.str().c_str());
  //}
}



ParamSet
BPSolver::getPosterioriOf (const Variable* var) const
{
  assert (var);
  assert (var == bn_->getNode (var->getVarId()));
  assert (var->getIndex() < msgs_.size());
  return msgs_[var->getIndex()]->getBeliefs();
}




ParamSet
BPSolver::getJointDistribution (const NodeSet& jointVars) const
{
  if (DL >= 1) {
    cout << "calculating joint distribuition on: " ;
    for (unsigned i = 0; i < jointVars.size(); i++) {
      cout << jointVars[i]->getLabel() << " " ;
    }
    cout << endl;
  }

  //BayesNet* workingNet = bn_->pruneNetwork (bn_->getNodes());
  //FIXME see if this works:
  BayesNet* workingNet = bn_->pruneNetwork (jointVars);
  BayesNode* node      = workingNet->getNode (jointVars[0]->getVarId());

  BayesNet* tempNet    = workingNet->pruneNetwork (node);
  BPSolver solver (*tempNet);
  solver.runSolver();

  NodeSet observedVars = { jointVars[0] };

  node = tempNet->getNode (jointVars[0]->getVarId());
  ParamSet prevBeliefs = solver.getPosterioriOf (node);

  delete tempNet;

  for (unsigned i = 1; i < jointVars.size(); i++) {
    node = workingNet->getNode (observedVars[i - 1]->getVarId());
    if (!node->hasEvidence()) {
      node->setEvidence (0);
    }
    node = workingNet->getNode (jointVars[i]->getVarId());
    tempNet = workingNet->pruneNetwork (node);

    ParamSet allBeliefs;
    vector<DomainConf> confs = 
        BayesNet::getDomainConfigurationsOf (observedVars);
    for (unsigned j = 0; j < confs.size(); j++) {
      for (unsigned k = 0; k < observedVars.size(); k++) {
        node = tempNet->getNode (observedVars[k]->getVarId());
        if (!observedVars[k]->hasEvidence()) {
          if (node) {
            node->setEvidence (confs[j][k]);
          } else {
            // FIXME try optimize
            //assert (false);
            cout << observedVars[k]->getLabel();
            cout << " is not in temporary net!" ;
            cout << endl;
          }
        } else {
          cout << observedVars[k]->getLabel();
          cout << " already has evidence in original net!" ;
          cout << endl;
        }
      }
      BPSolver solver (*tempNet);
      node = tempNet->getNode (jointVars[i]->getVarId());
      solver.runSolver();
      ParamSet beliefs = solver.getPosterioriOf (node);
      for (unsigned k = 0; k < beliefs.size(); k++) {
        allBeliefs.push_back (beliefs[k]);
      }
    }

    int count = -1;
    for (unsigned j = 0; j < allBeliefs.size(); j++) {
      if (j % jointVars[i]->getDomainSize() == 0) {
        count ++;
      }
      allBeliefs[j] *= prevBeliefs[count];
    }
    prevBeliefs = allBeliefs;
    observedVars.push_back (jointVars[i]);
    delete tempNet;
  }
  delete workingNet;
  return prevBeliefs;
}



void 
BPSolver::initializeSolver (void)
{
  if (DL >= 1) {
    cout << "Initializing solver" << endl;
    cout << "-> schedule        = ";
    if (forceGenericSolver_) {
      switch (schedule_) {
        case S_SEQ_FIXED:    cout << "sequential fixed" ;  break;
        case S_SEQ_RANDOM:   cout << "sequential random" ; break;
        case S_PARALLEL:     cout << "parallel" ;          break;
        case S_MAX_RESIDUAL: cout << "max residual" ;      break;
      }
    } else {
      cout << "polytree solver" ;
    }
    cout << endl;
    cout << "-> max iters       = " << maxIter_ << endl;
    cout << "-> accuracy        = " << accuracy_ << endl;
    cout << endl;
  }

  const NodeSet& nodes = bn_->getNodes();
  for (unsigned i = 0; i < msgs_.size(); i++) {
    delete msgs_[i];
  }
  msgs_.clear();
  msgs_.reserve (nodes.size());
  updateOrder_.clear();
  sortedOrder_.clear();
  edgeMap_.clear();

  for (unsigned i = 0; i < nodes.size(); i++) {
    msgs_.push_back (new BpNode (nodes[i]));
  }

  NodeSet roots = bn_->getRootNodes();
  for (unsigned i = 0; i < roots.size(); i++) {
    const ParamSet& params = roots[i]->getParameters();
    ParamSet& piVals = M(roots[i])->getPiValues();
    for (int ri = 0; ri < roots[i]->getDomainSize(); ri++) {
      piVals[ri] = params[ri];
    }
  }
}



void 
BPSolver::incorporateEvidence (BayesNode* x)
{
  ParamSet& piVals = M(x)->getPiValues();
  ParamSet& ldVals = M(x)->getLambdaValues();
  for (int xi = 0; xi < x->getDomainSize(); xi++) {
    piVals[xi] = 0.0;
    ldVals[xi] = 0.0;
  }
  piVals[x->getEvidence()] = 1.0;
  ldVals[x->getEvidence()] = 1.0;
}



void
BPSolver::runPolyTreeSolver (void)
{
  initializeSolver();
  const NodeSet& nodes = bn_->getNodes();

  // Hack: I need this else this can happen with bayes ball
  // Variable: 174
  // Id:       174
  // Domain:   -1, 0, 1
  // Evidence: 1
  // Parents:  
  // Childs:   176
  // cpt
  // ----------------------------------------------------
  // -1             0           0           0           0 ...
  // 0       0.857143    0.857143    0.857143    0.857143 ...
  // 1       0.142857    0.142857    0.142857    0.142857 ...
  // the cpt for this node would be 0,0,0

  for (unsigned i = 0; i < nodes.size(); i++) {
    if (nodes[i]->hasEvidence()) {
      incorporateEvidence (nodes[i]);
    }
  } 

  // first compute all node marginals ...
  NodeSet roots = bn_->getRootNodes();
  for (unsigned i = 0; i < roots.size(); i++) {
    const NodeSet& childs = roots[i]->getChilds();
    for (unsigned j = 0; j < childs.size(); j++) {
      polyTreePiMessage (roots[i], childs[j]);
    }
  }
  // then propagate the evidence
  for (unsigned i = 0; i < nodes.size(); i++) {
    if (nodes[i]->hasEvidence()) {
      incorporateEvidence (nodes[i]); 
      const NodeSet& parents = nodes[i]->getParents();
      for (unsigned j = 0; j < parents.size(); j++) {
        if (!parents[j]->hasEvidence()) {
          polyTreeLambdaMessage (nodes[i], parents[j]);
        }
      }
      const NodeSet& childs = nodes[i]->getChilds();
      for (unsigned j = 0; j < childs.size(); j++) {
        polyTreePiMessage (nodes[i], childs[j]);
      }
    }
  }
}



void
BPSolver::polyTreePiMessage (BayesNode* z, BayesNode* x)
{
  if (DL >= 1) {
    cout << PI << " (" << z->getLabel();
    cout << " --> " << x->getLabel();
    cout << ")" << endl;
  }
  calculateNextPiMessage (z, x);
  updatePiMessage (z, x);

  if (!x->hasEvidence()) {
    updatePiValues (x);
    const NodeSet& xChilds = x->getChilds();
    for (unsigned i = 0; i < xChilds.size(); i++) {
      polyTreePiMessage (x, xChilds[i]);
    }
  }

  if (M(x)->hasReceivedChildInfluence()) {
    const NodeSet& xParents = x->getParents();
    for (unsigned i = 0; i < xParents.size(); i++) {
      if (xParents[i] != z && !xParents[i]->hasEvidence()) {
        polyTreeLambdaMessage (x, xParents[i]);
      }
    }
  }
}



void 
BPSolver::polyTreeLambdaMessage (BayesNode* y, BayesNode* x)
{
  if (DL >= 1) {
    cout << LD << " (" << y->getLabel();
    cout << " --> " << x->getLabel();
    cout << ")" << endl;
  }
  calculateNextLambdaMessage (y, x);
  updateLambdaMessage (y, x);
  updateLambdaValues (x);

  const NodeSet& xParents = x->getParents();
  for (unsigned i = 0; i < xParents.size(); i++) {
    if (!xParents[i]->hasEvidence()) {
      polyTreeLambdaMessage (x, xParents[i]);
    }
  }

  const NodeSet& xChilds = x->getChilds();
  for (unsigned i = 0; i < xChilds.size(); i++) {
    if (xChilds[i] != y) {
      polyTreePiMessage (x, xChilds[i]);
    }
  }
}



void
BPSolver::runGenericSolver()
{ 
  initializeSolver();
  const NodeSet& nodes = bn_->getNodes();
  for (unsigned i = 0; i < nodes.size(); i++) {
    if (nodes[i]->hasEvidence()) {
      incorporateEvidence (nodes[i]);
    }
  }
 
  for (unsigned i = 0; i < nodes.size(); i++) {
    // pi messages
    const NodeSet& childs = nodes[i]->getChilds();
    for (unsigned j = 0; j < childs.size(); j++) {
      updateOrder_.push_back (Edge (nodes[i], childs[j], PI_MSG));
    }
    // lambda messages
    const NodeSet& parents = nodes[i]->getParents();
    for (unsigned j = 0; j < parents.size(); j++) {
      if (!parents[j]->hasEvidence()) {
        updateOrder_.push_back (Edge (nodes[i], parents[j], LAMBDA_MSG));
      }
    }
  }

  nIter_ = 0;
  while (!converged() && nIter_ < maxIter_) {

    nIter_++;
    if (DL >= 1) {
      cout << endl;
      cout << "****************************************" ;
      cout << "****************************************" ;
      cout << endl;
      cout << " Iteration " << nIter_ << endl;
      cout << "****************************************" ;
      cout << "****************************************" ;
      cout << endl;
    }
 
    switch (schedule_) {

      case S_SEQ_RANDOM:
        random_shuffle (updateOrder_.begin(), updateOrder_.end());
        // no break
  
      case S_SEQ_FIXED:
        for (unsigned i = 0; i < updateOrder_.size(); i++) {
          calculateNextMessage (updateOrder_[i]);
          updateMessage (updateOrder_[i]);
          updateValues (updateOrder_[i]);
        }
        break;

      case S_PARALLEL:
        for (unsigned i = 0; i < updateOrder_.size(); i++) {
          calculateNextMessage (updateOrder_[i]);
        }
        for (unsigned i = 0; i < updateOrder_.size(); i++) {
          updateMessage (updateOrder_[i]);
          updateValues (updateOrder_[i]);
        }
        break;
 
      case S_MAX_RESIDUAL:
        maxResidualSchedule();
        break;

    }
  }
}



void
BPSolver::maxResidualSchedule (void)
{
  if (nIter_ == 1) {
    Edge::klass = this;
    for (unsigned i = 0; i < updateOrder_.size(); i++) {
      calculateNextMessage (updateOrder_[i]);
      updateResidual (updateOrder_[i]);
      SortedOrder::iterator it = sortedOrder_.insert (updateOrder_[i]);
      edgeMap_.insert (make_pair (updateOrder_[i].getId(), it));
    }
    return;
  }

  for (unsigned c = 0; c < sortedOrder_.size(); c++) {
    if (DL >= 1) {
      for (set<Edge, compare>::iterator it = sortedOrder_.begin();
          it != sortedOrder_.end(); it ++) {
        cout << it->toString() << " residual = " ;
        cout << getResidual (*it) << endl;
      }
    }
 
    set<Edge, compare>::iterator it = sortedOrder_.begin();
    Edge e = *it;
    if (getResidual (e) < accuracy_) {
      return;
    }
    updateMessage (e);
    updateValues (e);
    clearResidual (e);
    sortedOrder_.erase (it);
    assert (edgeMap_.find (e.getId()) != edgeMap_.end());
    edgeMap_.find (e.getId())->second = sortedOrder_.insert (e);

     // update the messages that depend on message source --> destination
    const NodeSet& childs = e.destination->getChilds();
    for (unsigned i = 0; i < childs.size(); i++) {
      if (childs[i] != e.source) {
        Edge neighbor (e.destination, childs[i], PI_MSG);
        calculateNextMessage (neighbor);
        updateResidual (neighbor);
        assert (edgeMap_.find (neighbor.getId()) != edgeMap_.end());
        EdgeMap::iterator iter = edgeMap_.find (neighbor.getId());
        sortedOrder_.erase (iter->second);
        iter->second = sortedOrder_.insert (neighbor);
      }
    }
    const NodeSet& parents = e.destination->getParents();
    for (unsigned i = 0; i < parents.size(); i++) {
      if (parents[i] != e.source && !parents[i]->hasEvidence()) {
        Edge neighbor (e.destination, parents[i], LAMBDA_MSG);
        calculateNextMessage (neighbor);
        updateResidual (neighbor);
        assert (edgeMap_.find (neighbor.getId()) != edgeMap_.end());
        EdgeMap::iterator iter = edgeMap_.find (neighbor.getId());
        sortedOrder_.erase (iter->second);
        iter->second = sortedOrder_.insert (neighbor);
      }
    }
  }
}



bool
BPSolver::converged (void) const
{
  bool converged = true;
  if (schedule_ == S_MAX_RESIDUAL) {
    if (nIter_ <= 2) {
      return false;
    }
    // this can happen if every node does not have neighbors
    if (sortedOrder_.size() == 0) {
      return true;
    }
    Param maxResidual = getResidual (*(sortedOrder_.begin()));
    if (maxResidual > accuracy_) {
      return false;
    }
  } else {
    if (nIter_ == 0) {
      return false;
    }
    const NodeSet& nodes = bn_->getNodes();
    for (unsigned i = 0; i < nodes.size(); i++) {
      if (!nodes[i]->hasEvidence()) {
        double change = M(nodes[i])->getBeliefChange();
        if (DL >= 1) {
          cout << nodes[i]->getLabel() + " belief change = " ;
          cout << change << endl;
        }
        if (change > accuracy_) {
          converged = false;
          if (DL == 0) break;
        }
      }
    }
  }

  return converged;
}



void
BPSolver::updatePiValues (BayesNode* x)
{
  // π(Xi)
  const NodeSet& parents = x->getParents();
  const vector<CptEntry>& entries = x->getCptEntries();
  assert (parents.size() != 0);
  stringstream* calcs1;
  stringstream* calcs2;

  ParamSet messageProducts (entries.size());
  for (unsigned k = 0; k < entries.size(); k++) {
    if (DL >= 5) {
      calcs1 = new stringstream;
      calcs2 = new stringstream;
    }
    double messageProduct = 1.0;
    const DomainConf& conf = entries[k].getParentConfigurations();
    for (unsigned i = 0; i < parents.size(); i++) {
      messageProduct *= M(parents[i])->getPiMessageValue(x, conf[i]);
      if (DL >= 5) {
        if (i != 0) *calcs1 << "." ;
        if (i != 0) *calcs2 << "*" ;
        *calcs1 << PI << "(" << x->getLabel() << ")" ;
        *calcs1 << "[" << parents[i]->getDomain()[conf[i]] << "]";
        *calcs2 << M(parents[i])->getPiMessageValue(x, conf[i]);
      }
    }
    messageProducts[k] = messageProduct;
    if (DL >= 5) {
      cout << "    mp" << k;
      cout << " = " << (*calcs1).str();
      if (parents.size() == 1) {
        cout << " = " << messageProduct << endl;
      } else {
        cout << " = " << (*calcs2).str();
        cout << " = " << messageProduct << endl;
      }
      delete calcs1;
      delete calcs2;
    }
  }

  for (int xi = 0; xi < x->getDomainSize(); xi++) {
    double sum = 0.0;
    if (DL >= 5) {
      calcs1 = new stringstream;
      calcs2 = new stringstream;
    }
    for (unsigned k = 0; k < entries.size(); k++) {
      sum += x->getProbability (xi, entries[k]) * messageProducts[k];
      if (DL >= 5) {
        if (k != 0) *calcs1 << " + " ;
        if (k != 0) *calcs2 << " + " ;
        *calcs1 << x->cptEntryToString (xi, entries[k]); 
        *calcs1 << ".mp" << k;
        *calcs2 << x->getProbability (xi, entries[k]);
        *calcs2 << "*" << messageProducts[k];
      }
    }
    M(x)->setPiValue (xi, sum);
    if (DL >= 5) {
      cout << "    " << PI << "(" << x->getLabel() << ")" ;
      cout << "[" << x->getDomain()[xi] << "]" ;
      cout << " = " << (*calcs1).str();
      cout << " = " << (*calcs2).str();
      cout << " = " << sum << endl;
      delete calcs1;
      delete calcs2;
    }
  }
}



void
BPSolver::updateLambdaValues (BayesNode* x)
{
  // λ(Xi)
  const NodeSet& childs = x->getChilds();
  assert (childs.size() != 0);
  stringstream* calcs1;
  stringstream* calcs2;

  for (int xi = 0; xi < x->getDomainSize(); xi++) {
    double product = 1.0;
    if (DL >= 5) {
      calcs1 = new stringstream;
      calcs2 = new stringstream;
    }
    for (unsigned i = 0; i < childs.size(); i++) {
      product *= M(x)->getLambdaMessageValue(childs[i], xi);
      if (DL >= 5) {
        if (i != 0) *calcs1 << "." ;
        if (i != 0) *calcs2 << "*" ;
        *calcs1 << LD << "(" << childs[i]->getLabel();
        *calcs1 << "-->" << x->getLabel() << ")" ;
        *calcs1 << "[" << x->getDomain()[xi] << "]" ;
        *calcs2 << M(x)->getLambdaMessageValue(childs[i], xi);
      }
    }
    M(x)->setLambdaValue (xi, product);
    if (DL >= 5) {
      cout << "    " << LD << "(" << x->getLabel() << ")" ;
      cout << "[" << x->getDomain()[xi] << "]" ;
      cout << " = " << (*calcs1).str();
      if (childs.size() == 1) {
        cout << " = " << product << endl;
      } else {
        cout << " = " << (*calcs2).str();
        cout << " = " << product << endl;
      }
      delete calcs1;
      delete calcs2;
    }
  }
}



void
BPSolver::calculateNextPiMessage (BayesNode* z, BayesNode* x)
{
  // πX(Zi)
  ParamSet& zxPiNextMessage = M(z)->piNextMessageReference (x);
  const NodeSet& zChilds = z->getChilds();
  stringstream* calcs1;
  stringstream* calcs2;

  for (int zi = 0; zi < z->getDomainSize(); zi++) {
    double product = M(z)->getPiValue (zi);
    if (DL >= 5) {
      calcs1 = new stringstream;
      calcs2 = new stringstream;
      *calcs1 << PI << "(" << z->getLabel() << ")";
      *calcs1 << "[" << z->getDomain()[zi] << "]" ;
      *calcs2 << product;
    }
    for (unsigned i = 0; i < zChilds.size(); i++) {
      if (zChilds[i] != x) {
        product *= M(z)->getLambdaMessageValue(zChilds[i], zi);
        if (DL >= 5) {
          *calcs1 << "." << LD << "(" << zChilds[i]->getLabel();
          *calcs1 << "-->" << z->getLabel() << ")";
          *calcs1 << "[" << z->getDomain()[zi] + "]" ;
          *calcs2 << " * " << M(z)->getLambdaMessageValue(zChilds[i], zi);
        }
      }
    }
    zxPiNextMessage[zi] = product;
    if (DL >= 5) {
      cout << "    " << PI << "(" << z->getLabel();
      cout << "-->" << x->getLabel() << ")" ;
      cout << "["  << z->getDomain()[zi] << "]" ;
      cout << " = " << (*calcs1).str();
      if (zChilds.size() == 1) {
        cout << " = " << product << endl;
      } else {
        cout << " = " << (*calcs2).str();
        cout << " = " << product << endl;
      }
      delete calcs1;
      delete calcs2;
    }
  }
}



void
BPSolver::calculateNextLambdaMessage (BayesNode* y, BayesNode* x)
{
  // λY(Xi)
  //if (!y->hasEvidence() && !M(y)->hasReceivedChildInfluence()) {
  //  if (DL >= 5) {
  //    cout << "unnecessary calculation" << endl;
  //  }
  //  return;
  //}
  ParamSet& yxLambdaNextMessage       = M(x)->lambdaNextMessageReference (y);
  const NodeSet& yParents             = y->getParents();
  const vector<CptEntry>& allEntries  = y->getCptEntries();
  int parentIndex                     = y->getIndexOfParent (x);
  stringstream* calcs1;
  stringstream* calcs2;
 
  vector<CptEntry> entries;
  DomainConstr constr = make_pair (parentIndex, 0);
  for (unsigned i = 0; i < allEntries.size(); i++) {
    if (allEntries[i].matchConstraints(constr)) {
      entries.push_back (allEntries[i]);
    }
  }

  ParamSet messageProducts (entries.size());
  for (unsigned k = 0; k < entries.size(); k++) {
    if (DL >= 5) {
      calcs1 = new stringstream;
      calcs2 = new stringstream;
    }
    double messageProduct = 1.0;
    const DomainConf& conf = entries[k].getParentConfigurations();
    for (unsigned i = 0; i < yParents.size(); i++) {
      if (yParents[i] != x) {
        if (DL >= 5) {
          if (messageProduct != 1.0) *calcs1 << "*" ;
          if (messageProduct != 1.0) *calcs2 << "*" ;
          *calcs1 << PI << "(" << yParents[i]->getLabel();
          *calcs1 << "-->" << y->getLabel() << ")" ;
          *calcs1 << "[" << yParents[i]->getDomain()[conf[i]] << "]" ;
          *calcs2 << M(yParents[i])->getPiMessageValue(y, conf[i]);
        }
        messageProduct *= M(yParents[i])->getPiMessageValue(y, conf[i]);
      }
    }
    messageProducts[k] = messageProduct;
    if (DL >= 5) {
      cout << "    mp" << k;
      cout << " = " << (*calcs1).str();
      if  (yParents.size() == 1) {
        cout << 1 << endl;
      } else if (yParents.size() == 2) {
        cout << " = " << messageProduct << endl;
      } else {
        cout << " = " << (*calcs2).str();
        cout << " = " << messageProduct << endl;
      }
      delete calcs1;
      delete calcs2;
    }
  }

  for (int xi = 0; xi < x->getDomainSize(); xi++) {
    if (DL >= 5) {
      calcs1 = new stringstream;
      calcs2 = new stringstream;
    }
    vector<CptEntry> entries;
    DomainConstr constr = make_pair (parentIndex, xi);
    for (unsigned i = 0; i < allEntries.size(); i++) {
      if (allEntries[i].matchConstraints(constr)) {
        entries.push_back (allEntries[i]);
      }
    }
    double outerSum = 0.0;
    for (int yi = 0; yi < y->getDomainSize(); yi++) {
      if (DL >= 5) {
        (yi != 0) ? *calcs1 << " + {" : *calcs1 << "{" ;
        (yi != 0) ? *calcs2 << " + {" : *calcs2 << "{" ;
      } 
      double innerSum = 0.0;
      for (unsigned k = 0; k < entries.size(); k++) {
        if (DL >= 5) {
          if (k != 0) *calcs1 << " + " ;
          if (k != 0) *calcs2 << " + " ;
          *calcs1 << y->cptEntryToString (yi, entries[k]);
          *calcs1 << ".mp" << k;
          *calcs2 << y->getProbability (yi, entries[k]);
          *calcs2 << "*" << messageProducts[k];
        }
        innerSum += y->getProbability (yi, entries[k]) * messageProducts[k];
      }
      outerSum += innerSum * M(y)->getLambdaValue (yi);
      if (DL >= 5) {
        *calcs1 << "}." << LD << "(" << y->getLabel() << ")" ;
        *calcs1 << "["  << y->getDomain()[yi] << "]";
        *calcs2 << "}*" << M(y)->getLambdaValue (yi);
      }
    }
    yxLambdaNextMessage[xi] = outerSum;
    if (DL >= 5) {
      cout << "    " << LD << "(" << y->getLabel();
      cout << "-->" << x->getLabel() << ")" ;
      cout << "[" << x->getDomain()[xi] << "]" ;
      cout << " = " << (*calcs1).str();
      cout << " = " << (*calcs2).str();
      cout << " = " << outerSum << endl;
      delete calcs1;
      delete calcs2;
    }
  }
}



void
BPSolver::printMessageStatusOf (const BayesNode* var) const
{
  cout << left;
  cout << setw (10) << "domain" ;
  cout << setw (20) << PI << "(" + var->getLabel() + ")" ;
  cout << setw (20) << LD << "(" + var->getLabel() + ")" ;
  cout << setw (16) << "belief" ;
  cout << endl;
  cout << "--------------------------------" ;
  cout << "--------------------------------" ;
  cout << endl;

  BpNode*         x        = M(var);
  ParamSet&       piVals   = x->getPiValues();
  ParamSet&       ldVals   = x->getLambdaValues();
  ParamSet        beliefs  = x->getBeliefs();
  const Domain&   domain   = var->getDomain();
  const NodeSet&  childs   = var->getChilds();

  for (int xi = 0; xi < var->getDomainSize(); xi++) {
    cout << setw (10) << domain[xi];
    cout << setw (19) << piVals[xi];
    cout << setw (19) << ldVals[xi];
    cout.precision (PRECISION);
    cout << setw (16) << beliefs[xi];
    cout << endl;
  }
  cout << endl;
  if (childs.size() > 0) {
    string s = "(" + var->getLabel() + ")" ;
    for (unsigned j = 0; j < childs.size(); j++) {
      cout << setw (10) << "domain" ;
      cout << setw (28) << PI + childs[j]->getLabel() + s;
      cout << setw (28) << LD + childs[j]->getLabel() + s;
      cout << endl;
      cout << "--------------------------------" ;
      cout << "--------------------------------" ;
      cout << endl;
      const ParamSet& piMessage     = x->getPiMessage (childs[j]);
      const ParamSet& lambdaMessage = x->getLambdaMessage (childs[j]);
      for (int xi = 0; xi < var->getDomainSize(); xi++) {
        cout << setw (10) << domain[xi];
        cout.precision (PRECISION);
        cout << setw (27) << piMessage[xi];
        cout.precision (PRECISION);
        cout << setw (27) << lambdaMessage[xi];
        cout << endl;
      }
      cout << endl;
    }
  }
}



void 
BPSolver::printAllMessageStatus (void) const
{ 
  const NodeSet& nodes = bn_->getNodes();
  for (unsigned i = 0; i < nodes.size(); i++) {
    printMessageStatusOf (nodes[i]);
  }
}