#include <config.h>
#include <model/BUGSModel.h>
#include <model/TraceMonitor.h>
#include <model/NodeArray.h>
#include <graph/DevianceNode.h>
#include <graph/StochasticNode.h>
#include <graph/Node.h>

#include <list>
#include <utility>
#include <stdexcept>
#include <fstream>

#include <Rmath.h>

using std::vector;
using std::ofstream;
using std::list;
using std::pair;
using std::string;
using std::runtime_error;

BUGSModel::BUGSModel()
  : Model()
{
}

SymTab &BUGSModel::symtab() 
{
  return _symtab;
}

TraceMonitor const *BUGSModel::setMonitor(string const &name, Range const &range, unsigned int thin)
{
  NodeArray *array = _symtab.getVariable(name);
  if (!array) {
     throw runtime_error(string("Unknown variable ") + name);
  }
  Node *node = array->getSubset(range);
  if (node == 0) {
    throw runtime_error(string("Cannot monitor undefined node ") +
			array->name() + print(range));
  }
  if (!graph().contains(node)) {
    addExtraNode(node); 
  }
  TraceMonitor const *monitor = Model::setMonitor(node, thin);
  BUGSMonitor bugs_monitor;
  bugs_monitor.array = array;
  bugs_monitor.range = new Range(range);
  bugs_monitor.monitor = monitor;
  _bugs_monitors.push_back(bugs_monitor);
  return monitor;
}

void BUGSModel::clearMonitor(string const &name, Range const &range)
{
    NodeArray *array = _symtab.getVariable(name);
    if (!array) {
	throw runtime_error(string("Unknown variable ") + print(range));
    }
    if (array->range().contains(range)) {
	Node const *node = array->getSubset(range);
	if(node) {
	    vector<BUGSMonitor>::iterator p = _bugs_monitors.begin();
	    for (; p != _bugs_monitors.end(); ++p) {
		if (p->monitor->node() == node) {
		    delete p->range;
		    _bugs_monitors.erase(p);
		    break;
		}
	    }
	    Model::clearMonitor(node);	  
	}
    }
    else {
	throw runtime_error(string("Requested invalid subset ") + name
			    + print(range));
    }
}

static void writeDouble(double x, std::ostream &out)
{
  if (x == JAGS_NA) {
    out << "NA";
  }
  /* FIXME
     This shouldn't happen. On the other hand, I'm only commenting this
     out to avoid problems on MacOs X
  
  else if (ISNAN(x)) {
    out << "NaN";
  }
  */
  else if (!R_FINITE(x)) {
    if (x > 0) {
      out << "Inf";
    }
    else {
      out << "-Inf";
    }
  }
  else {
    out << x;
  }
}

void BUGSModel::coda(vector<Node const*> const &nodes, ofstream &output, 
		     ofstream &index)
{
  /* 
     Dump monitor for given node in CODA format.

     CODA requires two files, the output file (output) and the index
     file (index)
  */
  long lineno = 0;
  for (unsigned int i = 0; i < nodes.size(); i++) {
    for (vector<BUGSMonitor>::const_iterator j = _bugs_monitors.begin();
	 j != _bugs_monitors.end(); ++j) 
      {
	Node const *node = nodes[i];
	if (j->monitor->node() == node) {
	  TraceMonitor const *monitor = j->monitor;
	  int nvar = node->data.length();
	  if (nvar != 1) {
	    /* Multivariate node */
	    for (long offset = 0; offset < nvar; ++offset) {
	      index << j->array->name() << print(j->range->leftIndex(offset))
		    << " "  
		    << lineno + 1 << "  "  << lineno + monitor->size() << '\n';
	      long iter = monitor->start();
	      for (int k = 0; k < monitor->size(); k++) {
		output << iter << "  ";
		writeDouble(monitor->values()[k * nvar + offset], output);
		output << '\n';
		iter += monitor->thin();
	      }
	      lineno += monitor->size();
	    }
	  }
	  else {
	    /* Univariate node */
	    index << _symtab.getName(node) << "  " << lineno + 1 << "  " 
		  << lineno + monitor->size() << '\n';
	    long iter = monitor->start();
	    for (int k = 0; k < monitor->size(); k++) {
	      output << iter << "  ";
	      writeDouble(monitor->values()[k], output);
	      output << '\n';
	      iter += monitor->thin();
	    }
	    lineno += monitor->size();
	  }
	}
      }
  }
}

void BUGSModel::addDevianceNode()
{
  NodeArray const *array = _symtab.getVariable("deviance");
  if (array)
    return; //Deviance already defined by user

  _symtab.addVariable("deviance", Index(1));
  NodeArray *deviance = _symtab.getVariable("deviance");
  vector<Node*> nodes;
  graph().getNodes(nodes);
  std::set<StochasticNode*> parameters;
  for (vector<Node*>::const_iterator p = nodes.begin(); p != nodes.end(); ++p)
    {
      if (isObserved(*p)) {
        StochasticNode *snode = dynamic_cast<StochasticNode*>(*p);
        if (snode)
          parameters.insert(snode);
      }
    }
  DevianceNode *dnode = new DevianceNode(parameters);
  deviance->insert(dnode, Index(1));
  addExtraNode(dnode);
}