From a86b1d918a4c31f74024aa6b42fa0021919585de Mon Sep 17 00:00:00 2001
From: Bertrand <bertrand.horel@gmail.com>
Date: Mon, 8 Feb 2016 14:40:20 +0000
Subject: débuggage et seeding
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 main.cpp                |  28 +++++++----
 stratified_sampling.cpp |  38 ---------------
 stratified_sampling.hpp | 126 ++++++++++++++++++++++++++++++++++++++----------
 3 files changed, 119 insertions(+), 73 deletions(-)

diff --git a/main.cpp b/main.cpp
index 88617c4..8692967 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,16 +1,16 @@
 #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>
 
+using namespace std;
 //--génération quantiles--
-std::vector<double> quantile_norm(int n, double sigma){
-    std::vector<double> q(n);
+vector<double> quantile_norm(int n, double sigma){
+    vector<double> q(n);
     for (int i=0; i<n; i++) {
         q[i] = gsl_cdf_gaussian_Pinv ((double)(i+1)/n, sigma);
     }
@@ -21,15 +21,23 @@ int main()
 {
     //--- GSL random init ---
     gsl_rng_env_setup();
-    std::vector<double> q;
-    q = quantile_norm(10, 1);
+    vector<double> q = quantile_norm(10, 1);
     vector<double> p(10, 0.1);
-    vector<gaussian_truncated> gen;
-    gen.push_back(gaussian_truncated(GSL_NEGINF,q[0]));
+    vector<gaussian_truncated> rvar;
+    rvar.push_back(gaussian_truncated(GSL_NEGINF, q[0],0,1,0));
     for (int i=1; i<10; i++){
-        gen.push_back(gaussian_truncated(q[i-1],q[i]));
+        rvar.push_back(gaussian_truncated(q[i-1], q[i],0,1,i));
+    }
+    stratified_sampling<gaussian_truncated> S(p,rvar);
+    S.update(100);
+    S.draw();
+    for(int i=0;i<10;i++){
+        cout<<S.get_mean()[i]<<endl;
+    }
+    S.update(500);
+    S.draw();
+    for(int i=0;i<10;i++){
+        cout<<S.get_mean()[i]<<endl;
     }
-    stratified_sampling<gaussian_truncated> S(p, gen);
-    //S.update(1000);
     return 0;
 }
diff --git a/stratified_sampling.cpp b/stratified_sampling.cpp
index b959776..c1c8dc1 100644
--- a/stratified_sampling.cpp
+++ b/stratified_sampling.cpp
@@ -1,6 +1,4 @@
 #include "stratified_sampling.hpp"
-#include <algorithm>
-#include <iostream>
 
 std::pair<double, double> mean_var( std::vector<double> r){
     std::pair<double, double> p;
@@ -13,39 +11,3 @@ std::pair<double, double> mean_var( std::vector<double> r){
     p.second -= p.first * p.first;
     return p;
 }
-
-//actualisation du nombre de tirages à faire par strates
-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<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];
-    }
-}
diff --git a/stratified_sampling.hpp b/stratified_sampling.hpp
index b63cd52..0b0b4ef 100644
--- a/stratified_sampling.hpp
+++ b/stratified_sampling.hpp
@@ -1,4 +1,6 @@
 #include <vector>
+#include <algorithm>
+#include <iostream>
 #include <gsl/gsl_rng.h>
 #include "rtnorm.hpp"
 
@@ -9,43 +11,117 @@ struct var_alea {
     typedef T result_type;
     var_alea() : value(0) {};
     var_alea(T value) : value(value) {};
-    virtual ~var_alea() {}; 
+    virtual ~var_alea() {};
     virtual T operator()() = 0;
     T current() const { return value; };
-    protected: 
-    	T 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); }
+    gaussian_truncated(double a, double b, double mean=0, double sigma=1, int seed=0)
+        :a(a), b(b), mean(mean), sigma(sigma), seed(seed) {
+        const gsl_rng_type* type = gsl_rng_default;
+        gen = gsl_rng_alloc(type);
+        gsl_rng_set(gen, seed);
+    };
+    gaussian_truncated(gaussian_truncated const &other)
+        :a(other.a),b(other.b), mean(other.mean),sigma(other.sigma){
+        gen = gsl_rng_clone(other.gen);
+    };
+
+    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;
+    int seed;
+    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();
+    stratified_sampling(vector<double> p, vector<Gen> gen)
+        :p(p), gen(gen), mean(p.size(), 0), sigma(p.size(), 0){};
+    void update(int N);
+    void draw();
+    vector<double> get_mean();
+    //double estimator();
 private:
-  vector<double> p;
-  vector<int> M;
-  vector<double> sigma;
-  vector<Gen> gen;
- 		      
+    vector<double> p;
+    vector<int> M;
+    vector<int> cumM;
+    vector<double> mean;
+    vector<double> sigma;
+    vector<Gen> gen;
+
+};
+
+//actualisation du nombre de tirages à faire par strates
+template <typename Gen>
+void stratified_sampling<Gen>::update(int Nk) {
+    int I = p.size();
+    bool first_step = M.empty();
+    //reinitialistation du vecteur M du nombre de tirages par strates
+    if (first_step) {
+        M.resize(I,1);
+        cumM.resize(I,0);
+    }
+    else {
+        for(int i=0; i<I; i++){
+            cumM[i]=cumM[i]+M[i];
+            M[i]=1;
+        }
+    }
+
+    std::vector<double> m(I, 0); //le vecteur des m_i idéals
+
+    if (first_step) {
+        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];
+        cout<<M[i]<<"\t";
+    }
+    cout<<endl;
+}
+
+template <typename Gen>
+void stratified_sampling<Gen>::draw() {
+    int I = p.size();
+    double m, s, oldmean;
+    for(int i=0;i<I;i++){
+        m=0;
+        s=0;
+        for(int j=0;j<M[i];j++){
+            m=m+gen[i]();
+            s=s+gen[i].current()*gen[i].current();
+        }
+        oldmean=mean[i];
+        mean[i]=(mean[i]*cumM[i]+m)/(cumM[i]+M[i]);
+        sigma[i]=((sigma[i]+oldmean*oldmean)*cumM[i] + s)/(cumM[i]+M[i]) - mean[i]*mean[i];
+    }
+};
+
+template <typename Gen>
+vector<double> stratified_sampling<Gen>::get_mean() {
+    return mean;
 };
-- 
cgit v1.2.3-70-g09d2