#include <vector>
#include <gsl/gsl_rng.h>
#include "rtnorm.hpp"

using namespace std;

template <typename T>
struct var_alea {
    typedef T result_type;
    var_alea() : value(0) {};
    var_alea(T value) : value(value) {};
    virtual ~var_alea() {}; 
    virtual T operator()() = 0;
    T current() const { return value; };
    protected: 
    	T value;
};

typedef var_alea<double> var_alea_real;

struct gaussian_truncated : public var_alea_real
{
  gaussian_truncated(double a, double b, double mean=0, double sigma=1)
    :a(a), b(b), mean(mean), sigma(sigma) {
    const gsl_rng_type* type = gsl_rng_default;
    gen = gsl_rng_alloc(type);
  };
  double operator()() {
    pair<double, double> p = rtnorm(gen, a, b, mean, sigma);
    return value = p.first;
  };
  ~gaussian_truncated() { gsl_rng_free(gen); }
private:
  double mean, sigma, a, b;
  gsl_rng *gen;
};

template <typename Gen>
struct stratified_sampling {
  stratified_sampling(vector<double> p, vector<Gen> gen)
    :p(p), gen(gen) {};
  void update(int N);
  //vector<double> get_mean();
  //double estimator();
private:
  vector<double> p;
  vector<int> M;
  vector<double> sigma;
  vector<Gen> gen;
 		      
};