#ifndef CHAIN_H
#define CHAIN_H

#define _USE_MATH_DEFINES
#include <cmath>
#include "../treedef.h"
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/exponential_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <fstream>

class ChainManager;

class Chain
{
public:

	Chain(ChainManager * man, int id, bool b, bool bstrict, const std::string & treePath, const std::string & dataPath, const std::string & outputBasePath, double igammashape, double igammascale, long seed);
	~Chain();
	
	void setCycleLength(int i) {cycleLength = i;}
	double igamShape() const {return igamma_shape;}
	double igamScale() const {return igamma_scale;}
	std::string errorMessage() const {return errorMsg;}
	void setAlpha(double a) {stable_alpha = a;}
    void setScale(double c) {stable_scale = c;}
    void setBrownian (bool b) {brownian = b; if (brownian) stable_alpha=2;}
	Tree * getTree() {return tree;}
	Taxa * getTaxa() {return taxa;}
	
	double logL() const;
	double logCondL(double x, int numCons);
	
	double mcmcSampleAlpha(double logP);
    double mcmcSampleScale(double logP);
	double sliceSampleTraitValue(double curVal, double logP, int numCons);
	
	void run(int pos, int iters);
	
private:

    int myid;
    Tree * tree;
    Taxa * taxa;
    double igamma_shape;
    double igamma_scale;
    double stable_alpha;
    double stable_scale;
    double igamma_mode;
    boost::mt19937 rng;
    boost::variate_generator<boost::mt19937, boost::uniform_real<double> > unif01;
    boost::variate_generator<boost::mt19937, boost::exponential_distribution<double> > expdist1;
  	double lgamshape;
  	bool brownian;
  	bool brownianstrict;
  	double * cval;
  	double * clen;
  	int cycleLength;
  	std::string errorMsg;
  	std::ofstream output;
  	ChainManager * manager;
};

#endif // CHAIN_H
