#include <iostream>
#include <gsl/gsl_rng.h>
#include <vector>
#include "rtnorm.hpp"
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>
#define LOW_DISCREPANCY_HPP
#include <cmath>
#include <climits>
#include <list>
#include <algorithm>
#include <numeric>


//--génération quantiles--
std::vector<double> quantile_norm(int n, double sigma){
  std::vector<double> q(n);
  for (int i=0; i<n; i++) {
    q[i] = gsl_cdf_gaussian_Pinv ((double)(i+1)/n, sigma);
  }
  return q;
}

//--tirage de normale tronquée entre les quantile de taille 1/n i et i+1--
double quantile_truncate_normal (int i, int n, double mu, 
 double sigma, gsl_rng *gen) {
	 std::vector<double> q;
	 q = quantile_norm(n, sigma);
	 std::pair<double, double> p;
	 p = rtnorm (gen, q[i], q[i+1], mu, sigma);
	 return p.first;
 }

std::pair<double, double> mean_var( std::vector<double> r){
    std::pair<double, double> p;
    for(auto &x: r){
        p.first += x;
        p.second += x*x;
    }
    p.first /= r.size();
    p.second /= r.size();
    p.second -= p.first * p.first;
    return p;
}
//actualisation du nombre de tirages à faire par strates
std::vector<int> update_sampling (std::vector<double> p, 
std::vector<double> sigma, int n) {
	std::vector<int> r; // notre vecteur final à retourner
	std::vector<double> m; //le vecteur des Mi idéals 
	int s = p.size();
	if (sigma.empty()) {
		for (int i = 0; i<s ; i++) {
			m.push_back(double(n)*p[i]);
		}
	}
	else {
		double scal = std::inner_product(p.begin(), p.end(), sigma.begin(), (double) 0);
		for (int i=0; i<s; i++) {
			double v = p[i]*sigma[i]/scal;
			m.push_back(v*n);
			std::cout<<m[i]<<std::endl;
		}
	}
	r.push_back((int)(m[0]));
	double v1 = 0;
	for (int i=1; i<s-1; i++){
		for (int j=0; j<i; j++) {
			v1=v1+m[j];
		}
		r.push_back((int)v1 - (int)(v1-m[i]));
	}
	r.push_back(n - std::accumulate( r.begin(), r.end(), 0 ));
	return r;
}
			

int main()
{
  //--- GSL random init ---
  gsl_rng_env_setup();                          // Read variable environnement
  const gsl_rng_type* type = gsl_rng_default;   // Default algorithm 'twister'
  gsl_rng *gen = gsl_rng_alloc (type);          // Rand generator allocation
  std::vector<double> q;
  q = quantile_norm(10, 1);


  std::pair<double, double> p;
  std::pair<double, double> mv;
  //number of classes
  int I = 10;
  //number of samples
  int N = 10000;
  std::vector<double> r(N);
  double a;
  for (int i=0; i<I; i++){
      if(i==0){
          a = GSL_NEGINF;
      }else{
          a = q[i-1];
      }
      for(int j=0; j<N; j++){
          p = rtnorm (gen, a, q[i], 0, 1);
          r[j] = p.first;
      }
      mv = mean_var(r);
      //std::cout<<"mean :"<<mv.first<<" var :"<<mv.second<<std::endl;
  }
  std::vector<int> k;
  std::vector<double> z = {(double)1/3,(double)1/3,(double)1/3};
	std::vector<double> sigma = {0.1, 0.4, 0.3};
	k = update_sampling(z, sigma, 10000);
	for (int j=0; j<k.size(); j++){
		std::cout<<k[j]<<std::endl;
	}
  gsl_rng_free(gen);

		
  
  return 0;
}