yap-6.3/packages/CLPBN/horus/Util.cpp

#include <limits>

#include <sstream>
#include <fstream>

#include "Util.h"
#include "Indexer.h"
#include "ElimGraph.h"


namespace Globals {
bool logDomain = false;

unsigned verbosity = 0;

InfAlgorithms infAlgorithm = InfAlgorithms::VE;
};




namespace BpOptions {
Schedule schedule = BpOptions::Schedule::SEQ_FIXED;
//Schedule schedule = BpOptions::Schedule::SEQ_RANDOM;
//Schedule schedule = BpOptions::Schedule::PARALLEL;
//Schedule schedule = BpOptions::Schedule::MAX_RESIDUAL;
double    accuracy  = 0.0001;
unsigned  maxIter   = 1000;
}



namespace Util {


template <> std::string
toString (const bool& b)
{
  std::stringstream ss;
  ss << std::boolalpha << b;
  return ss.str();
}



unsigned
stringToUnsigned (string str)
{
  int val;
  stringstream ss;
  ss << str;
  ss >> val;
  if (val < 0) {
    cerr << "error: the readed number is negative" << endl;
    abort();
  }
  return static_cast<unsigned> (val);
}



double
stringToDouble (string str)
{
  double val;
  stringstream ss;
  ss << str;
  ss >> val;
  return val;
}




double
factorial (unsigned num)
{
  double result = 1.0;
  for (unsigned i = 1; i <= num; i++) {
    result *= i;
  }
  return result;
}



double
logFactorial (unsigned num)
{
  double result = 0.0;
  if (num < 150) {
    result = std::log (factorial (num));
  } else {
    for (unsigned i = 1; i <= num; i++) {
      result += std::log (i);
    }
  }
  return result;
}



unsigned
nrCombinations (unsigned n, unsigned k)
{
  assert (n >= k);
  int diff = n - k;
  unsigned result = 0;
  if (n < 150) {
    unsigned prod = 1;
    for (int i = n; i > diff; i--) {
      prod *= i;
    }
    result = prod / factorial (k);
  } else {
    double prod = 0.0;
    for (int i = n; i > diff; i--) {
      prod += std::log (i);
    }
    prod -= logFactorial (k);
    result = static_cast<unsigned> (std::exp (prod));
  }
  return result;
}



size_t
sizeExpected (const Ranges& ranges)
{
  return std::accumulate (
      ranges.begin(), ranges.end(), 1, multiplies<unsigned>());
}



unsigned
nrDigits (int num)
{
  unsigned count = 1;
  while (num >= 10) {
    num /= 10; 
    count ++;
  }
  return count;
}



bool
isInteger (const string& s)
{
  stringstream ss1 (s);
  stringstream ss2;
  int integer;
  ss1 >> integer;
  ss2 << integer;
  return (ss1.str() == ss2.str());
}



string
parametersToString (const Params& v, unsigned precision)
{
  stringstream ss;
  ss.precision (precision);
  ss << "[" ; 
  for (size_t i = 0; i < v.size(); i++) {
    if (i != 0) ss << ", " ;
    ss << v[i];
  }
  ss << "]" ;
  return ss.str();
}



vector<string>
getStateLines (const Vars& vars)
{
  Ranges ranges;
  for (size_t i = 0; i < vars.size(); i++) {
    ranges.push_back (vars[i]->range());
  }
  Indexer indexer (ranges);
  vector<string> jointStrings;
  while (indexer.valid()) {
    stringstream ss;
    for (size_t i = 0; i < vars.size(); i++) {
      if (i != 0) ss << ", " ;
      ss << vars[i]->label() << "=" ;
      ss << vars[i]->states()[(indexer[i])];
    }
    jointStrings.push_back (ss.str());
    ++ indexer;
  }
  return jointStrings;
}



bool
setHorusFlag (string key, string value)
{
  bool returnVal = true;
  if (key == "verbosity") {
    stringstream ss;
    ss << value;
    ss >> Globals::verbosity;
  } else if (key == "inf_alg") {
    if (       value == "ve") {
      Globals::infAlgorithm = InfAlgorithms::VE;
    } else if (value == "bp") {
      Globals::infAlgorithm = InfAlgorithms::BP;
    } else if (value == "cbp") {
      Globals::infAlgorithm = InfAlgorithms::CBP;
    } else {
      cerr << "warning: invalid value `" << value << "' " ;
      cerr << "for `" << key << "'" << endl;
      returnVal = false;
    }
  } else if (key == "elim_heuristic") {
    if (       value == "min_neighbors") {
      ElimGraph::elimHeuristic = ElimHeuristic::MIN_NEIGHBORS;
    } else if (value == "min_weight") {
      ElimGraph::elimHeuristic = ElimHeuristic::MIN_WEIGHT;
    } else if (value == "min_fill") {
      ElimGraph::elimHeuristic = ElimHeuristic::MIN_FILL;
    } else if (value == "weighted_min_fill") {
      ElimGraph::elimHeuristic = ElimHeuristic::WEIGHTED_MIN_FILL;
    } else {
      cerr << "warning: invalid value `" << value << "' " ;
      cerr << "for `" << key << "'" << endl;
      returnVal = false;
    }
  } else if (key == "schedule") {
    if (       value == "seq_fixed") {
      BpOptions::schedule = BpOptions::Schedule::SEQ_FIXED;
    } else if (value == "seq_random") {
      BpOptions::schedule = BpOptions::Schedule::SEQ_RANDOM;
    } else if (value == "parallel") {
      BpOptions::schedule = BpOptions::Schedule::PARALLEL;
    } else if (value == "max_residual") {
      BpOptions::schedule = BpOptions::Schedule::MAX_RESIDUAL;
    } else {
      cerr << "warning: invalid value `" << value << "' " ;
      cerr << "for `" << key << "'" << endl;
      returnVal = false;
    }
  } else if (key == "accuracy") {
    stringstream ss;
    ss << value;
    ss >> BpOptions::accuracy;
  } else if (key == "max_iter") {
    stringstream ss;
    ss << value;
    ss >> BpOptions::maxIter;
  } else if (key == "use_logarithms") {
    if (       value == "true") {
      Globals::logDomain = true;
    } else if (value == "false") {
      Globals::logDomain = false;
    } else {
      cerr << "warning: invalid value `" << value << "' " ;
      cerr << "for `" << key << "'" << endl;
      returnVal = false;
    }
  } else {
    cerr << "warning: invalid key `" << key << "'" << endl;
    returnVal = false;
  }
  return returnVal;
}



void
printHeader (string header, std::ostream& os)
{
  printAsteriskLine (os);
  os << header << endl;
  printAsteriskLine (os);
}



void
printSubHeader (string header, std::ostream& os)
{
  printDashedLine (os);
  os << header << endl;
  printDashedLine (os);
}



void
printAsteriskLine (std::ostream& os)
{
  os << "********************************" ;
  os << "********************************" ;
  os << endl;
}



void
printDashedLine (std::ostream& os)
{
  os << "--------------------------------" ;
  os << "--------------------------------" ;
  os << endl;
}


}



namespace LogAware {

void
normalize (Params& v)
{
  double sum = LogAware::addIdenty();
  if (Globals::logDomain) {
    for (size_t i = 0; i < v.size(); i++) {
      sum = Util::logSum (sum, v[i]);
    }
    assert (sum != -numeric_limits<double>::infinity());
    for (size_t i = 0; i < v.size(); i++) {
      v[i] -= sum;
    }
  } else {
    for (size_t i = 0; i < v.size(); i++) {
      sum += v[i];
    }
    assert (sum != 0.0);
    for (size_t i = 0; i < v.size(); i++) {
      v[i] /= sum;
    }
  }
}



double
getL1Distance (const Params& v1, const Params& v2)
{
  assert (v1.size() == v2.size());
  double dist = 0.0;
  if (Globals::logDomain) {
    for (size_t i = 0; i < v1.size(); i++) {
      dist += abs (exp(v1[i]) - exp(v2[i]));
    }
  } else {
    for (size_t i = 0; i < v1.size(); i++) {
      dist += abs (v1[i] - v2[i]);
    }
  }
  return dist;
}



double
getMaxNorm (const Params& v1, const Params& v2)
{
  assert (v1.size() == v2.size());
  double max = 0.0;
  if (Globals::logDomain) {
    for (size_t i = 0; i < v1.size(); i++) {
      double diff = abs (exp(v1[i]) - exp(v2[i]));
      if (diff > max) {
        max = diff;
      }
    }
  } else {
    for (size_t i = 0; i < v1.size(); i++) {
      double diff = abs (v1[i] - v2[i]);
      if (diff > max) {
        max = diff;
      }
    }
  }
  return max;
}



double
pow (double base, unsigned iexp)
{
  return Globals::logDomain ? base * iexp : std::pow (base, iexp);
}



double
pow (double base, double exp)
{
  // assumes that `expoent' is never in log domain
  return Globals::logDomain ? base * exp : std::pow (base, exp);
}



void
pow (Params& v, unsigned iexp)
{
  if (iexp == 1) {
    return;
  }
  Globals::logDomain ? v *= iexp : v ^= (int)iexp;
}



void
pow (Params& v, double exp)
{
  // `expoent' should not be in log domain
  Globals::logDomain ? v *= exp : v ^= exp;
}

}