encapsulage du stratified sampling

4 files changed, 123 insertions, 97 deletions

diff --git a/Makefile b/Makefile
index c412ee3..62e9c4e 100644
--- a/Makefile
+++ b/Makefile
@@ -4,11 +4,11 @@ RM = rm
 CXXFLAGS=-std=c++11 -O3
 GSL_FLAGS:=$(shell pkg-config --libs gsl)
 
-stratified_sampling: rtnorm.o stratified_sampling.o
+main: main.o rtnorm.o stratified_sampling.o
 	$(CXX) $^ -o $@ $(GSL_FLAGS)
 
 test: test.o mt19937.o
 	$(CXX) $^ -o $@
 
 clean:
-	-$(RM) -f *.o test stratified_sampling
+	-$(RM) -f *.o test stratified_sampling main
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..c1eb3ec
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,35 @@
+#include <iostream>
+#include <gsl/gsl_rng.h>
+#include <vector>
+//#include "rtnorm.hpp"
+#include <gsl/gsl_cdf.h>
+#include <gsl/gsl_math.h>
+#include "stratified_sampling.hpp"
+#include <cmath>
+#include <algorithm>
+
+//--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;
+}
+
+int main()
+{
+  //--- GSL random init ---
+  gsl_rng_env_setup();
+  std::vector<double> q;
+  q = quantile_norm(10, 1);
+  vector<double> p(10, 0.1);
+  vector<gaussian_truncated> gen;
+  gen.push_back(gaussian_truncated(GSL_NEGINF,q[0]));
+  for (int i=1; i<10; i++){
+    gen.push_back(gaussian_truncated(q[i-1],q[i]));
+  }
+  stratified_sampling<gaussian_truncated> S(p, gen);
+  //S.update(1000);
+  return 0;
+}
diff --git a/stratified_sampling.cpp b/stratified_sampling.cpp
index 8f408cc..c03d447 100644
--- a/stratified_sampling.cpp
+++ b/stratified_sampling.cpp
@@ -1,31 +1,6 @@
-#include <iostream>
-#include <gsl/gsl_rng.h>
-#include <vector>
-#include "rtnorm.hpp"
-#include <gsl/gsl_cdf.h>
-#include <gsl/gsl_math.h>
-
-#include <cmath>
+#include "stratified_sampling.hpp"
 #include <algorithm>
-
-//--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;
- }
+#include <iostream>
 
 std::pair<double, double> mean_var( std::vector<double> r){
     std::pair<double, double> p;
@@ -39,75 +14,40 @@ std::pair<double, double> mean_var( std::vector<double> r){
     return p;
 }
 //actualisation du nombre de tirages à faire par strates
-std::vector<int> update_sampling (std::vector<double> p, 
-				  std::vector<double> sigma, int Nk) {
-    int I = p.size();
-	std::vector<int> M(I, 1); // notre vecteur final à retourner
-	std::vector<double> m(I, 0); //le vecteur des m_i idéals
-
-    if (sigma.empty()) {
-		for (int i=0; i<I; i++) {
-			m[i] = (Nk-I)*p[i];
-		}
-	}
-	else {
-		double scal = std::inner_product(p.begin(), p.end(), sigma.begin(), (double) 0);
-		for (int i=0; i < I; i++) {
-			m[i] = (Nk-I)*p[i]*sigma[i]/scal;
-			std::cout<<m[i]<<std::endl;
-		}
-	}
-    M[0]+=floor(m[0]);
-    double current = m[0];
-	for (int i=1; i<I; i++){
-        M[i] += floor(current+m[i]) - floor(current);
-        current += m[i];
-	}
-	return M;
-}
-			
-
-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;
+template <typename Gen>
+void stratified_sampling<Gen>::update(int Nk) {
+  int I = p.size();
+  //reinitialistation du vecteur M du nombre de tirages par strates
+  if (M.empty()) {
+    M.resize(I,1);
+  }
+  else {
+    for(int i=0; i<I; i++){
+      M[i]=1;
+    }
   }
-  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);
+  
+  std::vector<double> m(I, 0); //le vecteur des m_i idéals
 
-		
+  if (sigma.empty()) {
+    for (int i=0; i<I; i++) {
+      m[i] = (Nk-I)*p[i];
+    }
+  }
+  else {
+    double scal = std::inner_product(p.begin(), p.end(), sigma.begin(), (double) 0);
+    for (int i=0; i < I; i++) {
+      m[i] = (Nk-I)*p[i]*sigma[i]/scal;
+      std::cout<<m[i]<<std::endl;
+    }
+  }
+  M[0]+=floor(m[0]);
+  double current = m[0];
+  for (int i=1; i<I; i++){
+    M[i] += floor(current+m[i]) - floor(current);
+    current += m[i];
+  }
   
-  return 0;
 }
+			
+
diff --git a/stratified_sampling.hpp b/stratified_sampling.hpp
new file mode 100644
index 0000000..b63cd52
--- /dev/null
+++ b/stratified_sampling.hpp
@@ -0,0 +1,51 @@
+#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;
+ 		      
+};