#include <config.h>
#include <distributions/DistReal.h>
#include <sarray/SArray.h>

#include <stdexcept>
#include <cfloat>
#include <cmath>
#include <algorithm>

#include <Rmath.h>

using std::max;
using std::min;
using std::string;
using std::vector;
using std::length_error;
using std::logic_error;


DistReal::DistReal(string const &name, unsigned int npar, Support support, 
		   bool canbound)
  : DistScalar(name, npar, support, canbound, false)
{
}

double DistReal::logLikelihood(SArray const &x,
			       vector<SArray const *> const &parameters) const
{
  SArray const *lb = lowerBound(parameters);
  SArray const *ub = upperBound(parameters);

  double y = x.value()[0];
  double loglik = d(y, parameters, true);
  if (lb && ub) {
    double ll = lowerSupport(0, parameters);
    if (y < ll)
      return -DBL_MAX;
    double uu = upperSupport(0, parameters); 
    if (y > uu) {
      return -DBL_MAX;
    }
    loglik -= log(p(uu, parameters, true, false) - 
		  p(ll, parameters, true, false));
  }
  else if (lb) {
    double ll = lowerSupport(0, parameters);
    if (y < ll)
      return -DBL_MAX;
    loglik -= p(ll, parameters, false, true);
  }
  else if (ub) {
    double uu = upperSupport(0, parameters); 
    if (y > uu) {
      return -DBL_MAX;
    }
    loglik -= p(uu, parameters, true, true);
  }
  return loglik;
}

void 
DistReal::randomSample(SArray &x,
		       vector<SArray const *> const &parameters) const
{

    SArray const *bb = lowerBound(parameters);
    SArray const *ba = upperBound(parameters);

    double y;
    if (!ba && !bb) {
	y = r(parameters);
    }
    else if (bb && ba) {
	double lower = lowerSupport(0, parameters);
	double plower = p(lower, parameters, false, false);
	double upper = upperSupport(0, parameters);
	double pupper = p(upper, parameters, false, false);
	double u = runif(plower, pupper);
	y = q(u, parameters, true, false);
    }
    else if (bb) {
	double lower = lowerSupport(0, parameters);
	double logplower = p(lower, parameters, false, true);
	double logpx = logplower - rexp(1);
	y = q(logpx, parameters, false, true);
    }
    else if (ba) {
	double upper = upperSupport(0, parameters);
	double logpupper = p(upper, parameters, true, true);
	double logpx = logpupper - rexp(1);
	y = q(logpx, parameters, true, true);
    }
    x.setValue(&y,1);
}

bool 
DistReal::checkParameterDim (vector<SArray const *> const &parameters) const
{
  for (unsigned int i = 0; i < parameters.size() - 2; i++) {
    if (parameters[i]->length() != 1)
      return false;
  }
  return true;
}