Commit a92d3b0c authored by Julien David's avatar Julien David
Browse files

First state revisited

parent ffe85080
Loading
Loading
Loading
Loading
+49 −33
Changes for src/RandomDistributionGenerator.hpp: 49 added lines, 33 removed lines.
Original line number Diff line number Diff line
@@ -4,6 +4,8 @@
#include"ConcaveFunction.hpp"
#include<assert.h>
#include<unistd.h>
#include<stdexcept>
#include <iomanip>

template<typename F = ConcaveFunction>
class RandomDistributionGenerator{
@@ -19,21 +21,6 @@ private:

  F cf;

  /**
   Initializes all the metrics of the generator's benchmarks.
   @param nb_xp is the number of conducted experiments
   @ensures all metrics are set to zero, except for the number of experiments
  */
  void init_benchmarks(){
    __AIM_STEPS__ = 0;
    __FIRST_STEP_FAILURES__ = 0;
    __INC_DIST_FAILURES__ = 0;
    __MARKOV_STEP_FAILURES__ = 0;
    __FIRST_STATE_ITERATIONS__ = 0;
    __NB_XP__ = 0;
    __NEWTON_STEPS__ = 0;
    __NEWTON_CALLS__ = 0;
  }

  void distribution_with_max_concave(long double * res, int n){
    int i;
@@ -78,9 +65,10 @@ private:
    // If the target is unreachable, exit function, only happens when precision error occurs
    //assert(sum<=1);
    if(aim_p < target || cf.aim_double(EPSILON, sum-EPSILON) > target) {
      //std::cout<<"FAILURE "<<aim_p<<" "<<target<<" "<<cf.aim_double(EPSILON, sum-EPSILON)<<std::endl;
      return -1.;
    }
    while (step > EPSILON && (double) aim_p != (double) target && p > 0){
    while ((double)step > (double)EPSILON && (double) aim_p != (double) target && p > 0){
      if(aim_p < target)
	p += step;
      else
@@ -102,52 +90,64 @@ private:
  void distribution_markov_chain(long double * res, long double target, int size, int steps){
    int i;

    while( reach_first_state(res, target, size) < 0)
    while( reach_first_state(res, target, size) < 0){
      __FIRST_STEP_FAILURES__++;
      if(__FIRST_STEP_FAILURES__ == 10)
	throw (target);
    }

    for (i = 0; i < steps; i++)
      markov_chain_step(res, size);        
  }

  
  
  int reach_first_state(long double * distribution, long double target, int size){
    long double incomplete_sum, contribution;
    long double min_aim, max_aim;
    long double min_aim, max_aim, update;
    long double delta = 0.5/(long double)size;
    distribution_with_max_concave(distribution, size);
    incomplete_sum = distribution[0]*(size-2);
    contribution = cf.contribution(distribution, size-2);
    min_aim = cf.minimal_of_two_values(1-incomplete_sum);
    max_aim = cf.maximal_of_two_values(1-incomplete_sum);
    
    long double update = cf.update_target(target, contribution, size, size-2);

    while( max_aim < update || min_aim > update ){
    distribution[size-2] = (1-incomplete_sum)/2;
    distribution[size-1] = (1-incomplete_sum)/2;
    max_aim = cf.h(cf.contribution(distribution, size), size);
    distribution[size-2] = (1-incomplete_sum-EPSILON);
    distribution[size-1] = EPSILON;
    min_aim = cf.h(cf.contribution(distribution, size), size);
    
    while( (max_aim < target || min_aim > target) && delta > EPSILON ){
      __FIRST_STATE_ITERATIONS__++;

      if( min_aim > target )
	add_delta(distribution, size-2, cf.sign()*delta);
	add_delta(distribution, size-2, -delta);
      else
	add_delta(distribution, size-2, -cf.sign()*delta);
	add_delta(distribution, size-2, delta);
      
      incomplete_sum = distribution[0]*(size-2);
      contribution = cf.contribution(distribution, size-2);
      min_aim = cf.minimal_of_two_values(1-incomplete_sum);
      max_aim = cf.maximal_of_two_values(1-incomplete_sum);
      distribution[size-2] = (1-incomplete_sum)/2;
      distribution[size-1] = (1-incomplete_sum)/2;
      max_aim = cf.h(cf.contribution(distribution, size),size);
      distribution[size-2] = (1-incomplete_sum-EPSILON);
      distribution[size-1] = EPSILON;
      min_aim = cf.h(cf.contribution(distribution, size),size);

      delta/=2;
      update = cf.update_target(target, contribution, size, size-2);
      //std::cout<<distribution[0]<<" "<<min_aim<<" "<<max_aim<<" IT:"<<__FIRST_STATE_ITERATIONS__<<" "<<delta<<std::endl;
      //std::cout<<"x:"<<distribution[0]<<" min:"<<min_aim<<" max:"<<max_aim<<" IT:"<<__FIRST_STATE_ITERATIONS__<<" delta:"<<delta<<std::endl;
    }
    contribution = cf.contribution(distribution, size-2);
    update = cf.update_target(target, contribution, size, size-2);
    distribution[size-2] = aim_concave_value(1-incomplete_sum, update);
    
    if(distribution[size-2] < 0)
      return -1;
  										
    distribution[size-1] = 1-incomplete_sum - distribution[size-2];
    //print_distribution(distribution, size);
    return 0; //If equal to -1, then a problem occured.
  }


  

  void markov_chain_step(long double * dist, int k){
    int triplet[3], i;
    long double target;
@@ -206,6 +206,22 @@ public:
      std::cerr<<"The target is too big to be reached"<<std::endl;
  }

  /**
   Initializes all the metrics of the generator's benchmarks.
   @param nb_xp is the number of conducted experiments
   @ensures all metrics are set to zero, except for the number of experiments
  */
  void init_benchmarks(){
    __AIM_STEPS__ = 0;
    __FIRST_STEP_FAILURES__ = 0;
    __INC_DIST_FAILURES__ = 0;
    __MARKOV_STEP_FAILURES__ = 0;
    __FIRST_STATE_ITERATIONS__ = 0;
    __NB_XP__ = 0;
    __NEWTON_STEPS__ = 0;
    __NEWTON_CALLS__ = 0;
  }


  void print_benchmarks(){
    print_benchmarks(__NB_XP__);