final with jade optimiazer in c++

src/generator-rf.cc

@@ -6,6 +6,49 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  void GeneratorRF::EvaluateApproximateDemodulated(std::vector<double> &mag_sin,
+                                        std::vector<double> &mag_cos) {
+    std::random_device rd;
+    std::mt19937_64 om;
+    om.seed(rd());
+    double lbound = 0.0, ubound = 1.0;
+    std::uniform_real_distribution<double> rand(lbound, ubound);
+
+    mag_sin.resize(total_samples_);
+    mag_cos.resize(total_samples_);
+    for (unsigned long i = 0; i < total_samples_; ++i) {
+      mag_sin[i] = 0.0;
+      mag_cos[i] = 0.0;
+      double time = timesteps_[i];
+      for (unsigned int voxel = 0; voxel < voxel_num_; ++voxel) {
+        mag_sin[i] +=
+          amplitudes_[voxel]
+          * std::sin(phases_[voxel]
+                     )
+          * ApproxExp(1,-time/
+                     (T2s_decays_norm_[voxel]*T2s_scale_
+                      )
+                     )
+          + (rand(om)-0.5)*noise_amplitude_
+
+          ;
+        mag_cos[i] +=
+          amplitudes_[voxel]
+          * std::cos(phases_[voxel]
+                     )
+          * ApproxExp(1,-time/
+                     (T2s_decays_norm_[voxel]*T2s_scale_
+                      )
+                     )
+          + (rand(om)-0.5)*noise_amplitude_
+          
+          ;
+       }
+    }
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   void GeneratorRF::EvaluateDemodulated(std::vector<double> &mag_sin,
                                         std::vector<double> &mag_cos) {
     std::random_device rd;
@@ -150,6 +193,45 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  double GeneratorRF::ApproxExp(int terms, double x){
+    if (terms % 2  == 0) 
+      throw std::invalid_argument
+        ("Exponential approximation is only available for even terms!");
+    
+    double result = 1.0;
+    for (int i = 1; i < terms+1; ++i) {
+      result += 1.0*std::pow(x,i)/std::tgamma(i+1);
+      //      std::cout << " %% " << std::tgamma(i+1)<< " "<< i <<" -> "<<result;
+    }
+    if (result < 0) return 0.0;
+    return result;
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double GeneratorRF::TryApproximateFit(const std::vector<double> &values) {
+    SetVoxels(values);
+    double diff = 0.0;
+    std::vector<double> mag_sin, mag_cos;
+    EvaluateApproximateDemodulated(mag_sin, mag_cos);
+    
+    for (unsigned long i = 0; i < mag_sin.size(); ++i) {
+      diff += std::abs((base_mag_sin_[i] - mag_sin[i]));
+      diff += std::abs((base_mag_cos_[i] - mag_cos[i]));
+
+      // diff += std::abs((base_mag_sin_[i] - mag_sin[i])
+      //                  / base_mag_sin_[i]);
+      // diff += std::abs((base_mag_cos_[i] - mag_cos[i])
+      //                  / base_mag_cos_[i]);
+      // diff += std::abs(
+      //                  (base_mag_sin_[i]/base_mag_cos_[i] - mag_sin[i]/mag_cos[i])
+      //                  );
+    }
+    return diff;
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   double GeneratorRF::TryFit(const std::vector<double> &values) {
     SetVoxels(values);
     std::vector<double> mag_sin, mag_cos;
@@ -159,10 +241,26 @@ namespace mri {
       // diff += std::abs((base_mag_sin_[i] - mag_sin[i]));
       // diff += std::abs((base_mag_cos_[i] - mag_cos[i]));
 
-      diff += std::abs((base_mag_sin_[i] - mag_sin[i])/mag_sin[i]);
-      diff += std::abs((base_mag_cos_[i] - mag_cos[i])/mag_cos[i]);
-
-      diff += std::abs(base_mag_sin_[i]/base_mag_cos_[i] - mag_sin[i]/mag_cos[i]);
+      diff += std::abs((base_mag_sin_[i] - mag_sin[i])
+                       / //(pow2(
+                              base_mag_sin_[i]
+                         //      )+
+                         // pow2(base_mag_cos_[i])
+                         // )
+                       );
+      diff += std::abs((base_mag_cos_[i] - mag_cos[i])
+                       / //(pow2(base_mag_sin_[i])+
+                         //pow2(
+                       base_mag_cos_[i]
+                       // )
+                       //   )
+                       );
+      diff += std::abs(
+                       (base_mag_sin_[i]/base_mag_cos_[i] - mag_sin[i]/mag_cos[i])
+                       // /((base_mag_sin_[i])+
+                       //   (base_mag_cos_[i])
+                       //   )                       
+                       );
     }
     return diff;
   };
@@ -189,6 +287,15 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  void GeneratorRF::InitApproximation() {
+    exp_terms_.resize(voxel_num_);
+    for (auto &term:exp_terms_) {
+      term = 1;
+    }
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   void GeneratorRF::InitVoxels(unsigned int voxel_num, double phase_init,
                               double phase_range, double T2s_scale) {
     T2s_scale_ = T2s_scale;
@@ -202,8 +309,8 @@ namespace mri {
     double phase_step = phase_range_/(voxel_num_ - 1);
     phases_.resize(voxel_num);
     for (unsigned int i =0; i< voxel_num_; ++i) {
-      //phases_[i] = phase_init_ + i*phase_step;
-      phases_[i] = phase_init_ + std::sin(1.0*i/voxel_num_);
+      phases_[i] = phase_init_ + i*phase_step;
+      //phases_[i] = phase_init_ + std::sin(1.0*i/voxel_num_);
       std::cout<<phases_[i]<<std::endl;
     }
   };

src/generator-rf.h

@@ -7,6 +7,7 @@ namespace mri {
   public:
     template<class T> inline T pow2(const T value) {return value*value;}
     const double PI_=3.14159265358979323846;
+    void InitApproximation();
     void InitTimer(unsigned long total_periods, double period_ms,
                    int samples_per_period);
     void InitVoxels(unsigned int voxel_num, double phase_init_,
@@ -14,6 +15,8 @@ namespace mri {
     void EvaluateSignal(std::vector<double> &signal_rf);
     void EvaluateDemodulated(std::vector<double> &mag_sin,
                              std::vector<double> &mag_cos);
+    void EvaluateApproximateDemodulated(std::vector<double> &mag_sin,
+                                        std::vector<double> &mag_cos);
     void RandomizeVoxels();
     void SetAmplitudes(std::vector<double> amplitudes);
     void SetBase(const std::vector<double> &mag_sin,
@@ -26,14 +29,17 @@ namespace mri {
     void SetVoxels(const std::vector<double> &values);
     void GetVoxels(std::vector<double> &values);
     
+    double TryApproximateFit(const std::vector<double> &values);
     double TryFit(const std::vector<double> &values);
     const std::vector<double>& GetTimesteps(){return timesteps_;};
 
   private:
+    double ApproxExp(int terms, double x);
     std::vector<double> base_mag_sin_, base_mag_cos_;
     std::vector<double> amplitudes_;
     std::vector<double> phases_;
     std::vector<double> T2s_decays_norm_;
+    std::vector<int> exp_terms_;
     double T2s_scale_ = 0.01; // #ms # need to be 10ms
     /* double T2s_min_ = T2s_scale_/1000.0; */
     

src/jade.cc

@@ -38,6 +38,8 @@
 #include <iterator>
 #include <string>
 #include <vector>
+
+
 namespace jade {
   /// @todo Replace all simple kError returns with something meangfull. TODO change kError to throw exception
   // ********************************************************************** //
@@ -168,6 +170,7 @@ namespace jade {
     evaluated_fitness_for_next_generation_ =
       evaluated_fitness_for_current_vectors_;
     for (long g = 0; g < total_generations_max_; ++g) {
+      //if (flag__ == 1) break;
       //if (process_rank_ == kOutput && g%100 == 0) {
       if (g%10 == 0) {
         double fitness = evaluated_fitness_for_current_vectors_.front().first;

src/main.cc

@@ -12,9 +12,11 @@
 #include <stdexcept>
 #include <string>
 
+
+
 const double PI_=3.14159265358979323846;
 const double eps_=1e-11;
-long total_generations_ = 15000;
+long total_generations_ = 10000;
 double noise_amp = 1e-8;
 
 // int voxel_num = 10;
@@ -30,7 +32,7 @@ double vox[] = {0.822628,0.691376,0.282906,0.226013,0.90703,0.144985,0.328563,0.
 // //double vox[] = {0.822628,0.691376,0.282906,0.226013,0.90703,0.144985};
 // double vox[] = {0.822628,0.691376,0.282906,0.226013,0.30703,0.244985};
 
-double total_periods = 20;
+double total_periods = 40;
 int rf_samples_per_period = 1;
 
 std::vector<double> init_vox (vox, vox + sizeof(vox) / sizeof(double) );
@@ -56,16 +58,19 @@ jade::SubPopulation sub_population_;  // Optimizer.
 
 void InitGenerators();
 void PrintVector(const std::vector<double> &values);
-void SetOptimizer();
+void SetOptimizer(long dimensions, std::vector<std::vector<double> > x_feed_vectors);
 
 double EvaluateFitness(const std::vector<double> &input){
-  return fitness_rf_.TryFit(input);
+  //  return fitness_rf_.TryFit(input);
+  return fitness_rf_.TryApproximateFit(input);
 };
 
 int main(int argc, char *argv[]) {
+  
   MPI_Init(&argc, &argv);
   int rank;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  
   try {
     
     InitGenerators();    
@@ -78,21 +83,23 @@ int main(int argc, char *argv[]) {
     generator_rf_.EvaluateDemodulated(mag_sin,mag_cos);
     fitness_rf_.SetBase(mag_sin, mag_cos);
     
-    SetOptimizer();
+    std::vector<std::vector<double> > x_feed_vectors;
+    x_feed_vectors.push_back(std::vector<double> (2*voxel_num, 0.5));
+    SetOptimizer(2*voxel_num, x_feed_vectors);
     sub_population_.FitnessFunction = &EvaluateFitness;
-
-    // double fit = fitness_rf_.TryFit(origin_values);
-    // std::cout << fit << std::endl;
-      
-    sub_population_.RunOptimization();
     double fit;
+        
+   
+    // fit = fitness_rf_.TryApproximateFit(origin_values);
+    // std::cout<< std::endl << "=========================================================  "<< fit << std::endl;      
+    sub_population_.RunOptimization();
     auto best_local_x = sub_population_.GetBest(&fit);
     std::cout << "origin values:\t";
     PrintVector(origin_values);
     std::cout << "   fit values:\t";
     PrintVector(best_local_x);
     std::vector<double> diff_vec(origin_values.size());
-    for (int i = 0; i < origin_values.size(); ++i)
+    for (unsigned int i = 0; i < origin_values.size(); ++i)
       diff_vec[i] = std::abs(origin_values[i]-best_local_x[i])*100;
     std::cout << "         diff:\t";
     PrintVector(diff_vec);
@@ -121,6 +128,7 @@ void InitGenerators() {
     generator_rf_.SetNoiseAmplitude(noise_amp);
     fitness_rf_.InitVoxels(voxel_num, phase_init, phase_range, T2s_scale);
     fitness_rf_.SetNoiseAmplitude(0.0);
+    fitness_rf_.InitApproximation();
     
     generator_rf_.InitTimer(total_periods, period_ms, rf_samples_per_period);
     fitness_rf_.InitTimer(total_periods, period_ms, rf_samples_per_period);
@@ -136,21 +144,21 @@ void PrintVector(const std::vector<double> &values) {
 // ********************************************************************** //
 // ********************************************************************** //
 // ********************************************************************** //
-void SetOptimizer() {
-  long dimension = 2*voxel_num;
-  long total_population = dimension * 20;
+void SetOptimizer(long dimension, std::vector<std::vector<double> > x_feed_vectors) {
+  long total_population = dimension * 30;
   sub_population_.Init(total_population, dimension);
   /// Low and upper bound for all dimenstions;
   sub_population_.SetAllBounds(eps_, 1.0-eps_);
   sub_population_.SetTargetToMinimum();
   sub_population_.SetTotalGenerationsMax(total_generations_);
+  sub_population_.SetFeed(x_feed_vectors);
   //sub_population_.SwitchOffPMCRADE();
 
   //sub_population_.SetBestShareP(0.1);
   //sub_population_.SetAdapitonFrequencyC(1.0/20.0);
   sub_population_.SetBestShareP(0.1);
   sub_population_.SetAdapitonFrequencyC(1.0/20.0);
-
+  
 }
 // ********************************************************************** //
 // ********************************************************************** //