mag sin cos fit 5 point tolerance 1e-14

go.sh

@@ -82,6 +82,6 @@ BuildMRI
 echo "Executing $mri_bin on host -> $HOST <-"
 MPI_options=
 #MPI_options="--bind-to-core"
-MPI_size=1
+MPI_size=2
 cd $path_bin
-time mpirun -np $MPI_size $MPI_options ./$mri_bin > out.txt 
+time mpirun -np $MPI_size $MPI_options ./$mri_bin 

src/generator-rf.cc

@@ -10,7 +10,7 @@ namespace mri {
                                         std::vector<double> &mag_cos) {
     mag_sin.resize(total_samples_);
     mag_cos.resize(total_samples_);
-    for (long i = 0; i < total_samples_; ++i) {
+    for (unsigned long i = 0; i < total_samples_; ++i) {
       mag_sin[i] = 0.0;
       mag_cos[i] = 0.0;
       double time = timesteps_[i];
@@ -41,7 +41,7 @@ namespace mri {
   // ********************************************************************** //
   void GeneratorRF::EvaluateSignal(std::vector<double> &signal_rf) {
     signal_rf.resize(total_samples_);
-    for (long i = 0; i < total_samples_; ++i) {
+    for (unsigned long i = 0; i < total_samples_; ++i) {
       signal_rf[i] = 0.0;
       double time = timesteps_[i];
       for (unsigned int voxel = 0; voxel < voxel_num_; ++voxel) {
@@ -74,37 +74,93 @@ namespace mri {
     T2s_decays_norm_.clear();
     amplitudes_.resize(voxel_num_);
     T2s_decays_norm_.resize(voxel_num_);
-    std::cout << "#";
     for (double &amp:amplitudes_){
       amp = rand(om);
-      std::cout << amp << "\t";
     }
     for (double &decay:T2s_decays_norm_){
       decay = rand(om);
-      std::cout << decay << "\t";
     }
-    std::cout<<std::endl;
   };
-
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
   void GeneratorRF::SetAmplitudes(std::vector<double> amplitudes){
     if (amplitudes.size() != voxel_num_)
-      throw std::invalid_argument("Number of amplitudes do not fit number of voxels");    
+      throw std::invalid_argument
+        ("Number of amplitudes do not fit number of voxels");
     amplitudes_ = amplitudes;
   };
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  void GeneratorRF::SetBase(const std::vector<double> &mag_sin,
+                            const std::vector<double> &mag_cos){
+    if ((mag_sin.size() != mag_cos.size())
+      ||
+        (mag_sin.size() != total_samples_)) {
+      throw std::invalid_argument
+        ("Base signal size does not fit settings of the fitness signal.");
+    }
+    base_mag_sin_.resize(total_samples_);
+    base_mag_cos_.resize(total_samples_);
+    for (unsigned long i = 0; i<total_samples_; ++i) {
+      base_mag_sin_[i] = mag_sin[i];
+      base_mag_cos_[i] = mag_cos[i];
+    }
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   void GeneratorRF::SetT2sDecaysNorm(std::vector<double> T2s_decays_norm){
     T2s_decays_norm_ = T2s_decays_norm;
   };
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
+  // amplitudes -> values[:size/2]
+  // T2s_decays_norm -> values[size/2:]
+  void GeneratorRF::SetVoxels(const std::vector<double> &values) {
+    for (unsigned int i = 0; i < voxel_num_; ++i) {
+      amplitudes_[i] = values[i];
+      T2s_decays_norm_[i] = values[i+voxel_num_];
+    }
+  }; 
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // amplitudes -> values[:size/2]
+  // T2s_decays_norm -> values[size/2:]
+  void GeneratorRF::GetVoxels(std::vector<double> &values) {
+    values.resize(2 * voxel_num_);
+    for (unsigned int i = 0; i < voxel_num_; ++i) {
+      values[i] = amplitudes_[i];
+      values[i+voxel_num_] = T2s_decays_norm_[i];
+    }
+  }; 
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
+  double GeneratorRF::TryFit(const std::vector<double> &values) {
+    SetVoxels(values);
+    std::vector<double> mag_sin, mag_cos;
+    EvaluateDemodulated(mag_sin, mag_cos);
+    double diff = 0.0;
+    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])/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]);
+    }
+    return diff;
+  };
+  // ********************************************************************** //
+  // ********************************************************************** //
+  // ********************************************************************** //
   // period is set in milliseconds (ms)
-  void GeneratorRF::SetTimer(long total_periods, double period_ms,
+  void GeneratorRF::InitTimer(unsigned long total_periods, double period_ms,
                              int samples_per_period) {
     total_periods_ = total_periods;
     period_ms_ = period_ms;
@@ -113,7 +169,7 @@ namespace mri {
     total_samples_ = total_periods * samples_per_period;
     timesteps_.resize(total_samples_);
     double a_timestep = period_ms/samples_per_period;
-    for (long i = 0; i < total_samples_; ++i) {
+    for (unsigned long i = 0; i < total_samples_; ++i) {
       timesteps_[i] = i*a_timestep;      
     }
     // for (double step: timesteps_) std::cout << step<<' ';
@@ -123,12 +179,14 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   // ********************************************************************** //
-  void GeneratorRF::SetVoxels(double voxel_num, double phase_init,
+  void GeneratorRF::InitVoxels(unsigned int voxel_num, double phase_init,
                               double phase_range, double T2s_scale) {
     T2s_scale_ = T2s_scale;
     if (voxel_num < 2)
       throw std::invalid_argument("We need at least two voxels for phase gradient");    
     voxel_num_ = voxel_num;
+    amplitudes_.resize(voxel_num_);
+    T2s_decays_norm_.resize(voxel_num_);
     phase_init_ = phase_init;
     phase_range_ = phase_range;
     double phase_step = phase_range_/(voxel_num_ - 1);

src/generator-rf.h

@@ -7,16 +7,29 @@ namespace mri {
   public:
     template<class T> inline T pow2(const T value) {return value*value;}
     const double PI_=3.14159265358979323846;
+    void InitTimer(unsigned long total_periods, double period_ms,
+                   int samples_per_period);
+    void InitVoxels(unsigned int voxel_num, double phase_init_,
+                    double phase_range_, double T2s_scale);
     void EvaluateSignal(std::vector<double> &signal_rf);
     void EvaluateDemodulated(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,
+                 const std::vector<double> &mag_cos);
     void SetT2sDecaysNorm(std::vector<double> T2s_decays_norm);
-    void SetTimer(long total_periods, double period_ms, int samples_per_period);
-    void SetVoxels(double voxel_num, double phase_init_, double phase_range_, double T2s_scale);
+
+    // amplitudes -> values[:size/2]
+    // T2s_decays_norm -> values[size/2:]
+    void SetVoxels(const std::vector<double> &values);
+    void GetVoxels(std::vector<double> &values);
+    
+    double TryFit(const std::vector<double> &values);
     const std::vector<double>& GetTimesteps(){return timesteps_;};
+
   private:
+    std::vector<double> base_mag_sin_, base_mag_cos_;
     std::vector<double> amplitudes_;
     std::vector<double> phases_;
     std::vector<double> T2s_decays_norm_;
@@ -32,9 +45,9 @@ namespace mri {
     // B0=1.5T freq=64Mhz, period = 15.6 ns
     double period_ms_ = 15.6/1000/1000; //ms
     double omega_ = 2.0*PI_/period_ms_; 
-    long total_periods_ = 1000;
+    unsigned long total_periods_ = 1000;
     int samples_per_period_ = 10;
-    long total_samples_ = total_periods_ * samples_per_period_;
+    unsigned long total_samples_ = total_periods_ * samples_per_period_;
 
     std::vector<double> timesteps_;
   };  // end of class GeneratorRF

src/jade.cc

@@ -168,7 +168,11 @@ namespace jade {
     evaluated_fitness_for_next_generation_ =
       evaluated_fitness_for_current_vectors_;
     for (long g = 0; g < total_generations_max_; ++g) {
-      if (process_rank_ == kOutput && g%100 == 0) printf("%li\n",g);      
+      //if (process_rank_ == kOutput && g%100 == 0) {
+      if (g%10 == 0) {
+        double fitness = evaluated_fitness_for_current_vectors_.front().first;
+        printf("%li -> %g\n",g,fitness);
+      }
       to_be_archived_best_A_.clear();
       successful_mutation_parameters_S_F_.clear();
       successful_crossover_parameters_S_CR_.clear();        
@@ -177,7 +181,7 @@ namespace jade {
       //   printf("==============  Generation %li =============\n", g);
       // PrintPopulation();      
       // PrintEvaluated();
-      //end of debug section
+      // //end of debug section
       for (long i = 0; i < subpopulation_; ++i) {
         SetCRiFi(i);
         std::vector<double> mutated_v, crossover_u;

src/jade.h

@@ -45,7 +45,7 @@ namespace jade {
   class SubPopulation {
    public:
     /// @brief Externaly defined fitness function, used by pointer.
-    double (*FitnessFunction)(std::vector<double> x) = nullptr;
+    double (*FitnessFunction)(const std::vector<double> &x) = nullptr;
     /// @brief Class initialization.
     int Init(long total_population, long dimension);              // NOLINT
     /// @brief Vizualize used random distributions (to do manual check).

src/lock-in.cc

@@ -21,8 +21,8 @@ namespace mri {
     // Stage 2: Low-pas filter:
     //FilterLowPassFirstOrder();
     //FilterLowPassForthOrder();
-    //FilterLowPassPeriodAverage();
-    FilterLowPassWindowedSinc();
+    FilterLowPassPeriodAverage();
+    //FilterLowPassWindowedSinc();
 
 
     unsigned long filtered_size = mag_sin_.size();
@@ -31,8 +31,6 @@ namespace mri {
     for (unsigned long i = 0; i<filtered_size; ++i) {
       out_sin[i] = mag_sin_[i];
       out_cos[i] = mag_cos_[i];
-      // out_sin[i] = mul_sin_[i];
-      // out_cos[i] = mul_cos_[i];
     }
     
     
@@ -44,10 +42,11 @@ namespace mri {
     mul_sin_.resize(total_samples_);
     mul_cos_.resize(total_samples_);
     for (unsigned long i = 0; i < total_samples_; ++i) {
-      // mul_sin_[i] = signal[i]*lockin_sin_[i];
-      // mul_cos_[i] = signal[i]*lockin_cos_[i];
-      mul_sin_[i] = lockin_sin_[i]+1;
-      mul_cos_[i] = lockin_cos_[i]+2;
+      mul_sin_[i] = signal[i]*lockin_sin_[i];
+      mul_cos_[i] = signal[i]*lockin_cos_[i];
+      //Test
+      // mul_sin_[i] = lockin_sin_[i]+1;
+      // mul_cos_[i] = lockin_cos_[i]+2;
     }
   };
   // ********************************************************************** //
@@ -133,12 +132,13 @@ namespace mri {
   // ********************************************************************** //
   // ********************************************************************** //
   void LockIn::FilterLowPassWindowedSinc() {
-    double fc = 0.05; // Cutoff frequency as a fraction of the sampling rate (in (0, 0.5)).
-    int N = 1001; // Kernel size, should be odd.
+    double fc = 0.005; // Cutoff frequency as a fraction of the sampling rate (in (0, 0.5)).
+    int N = 5001; // Kernel size, should be odd.
     int M = N-1;
     int mid = static_cast<int>(M/2);
     //Evaluate kernel
-    std::vector<double> h(N);
+    std::vector<double> h;
+    h.resize(N);
     for (int i = 0; i<N; ++i) {
       double x = static_cast<double>(i-mid);
       if (i == mid) {h[i] = 2.0*PI_*fc;}
@@ -153,12 +153,12 @@ namespace mri {
     for (auto &kern:h) kern/=sum;
     // for (int i = 0; i<N; ++i) std::cout<< h[i] <<std::endl;
     //Convolution
-    for (long j = 0; j < N; ++j) {
-      mag_sin_[j] = 0;
-      mag_cos_[j] = 0;
-      mul_sin_[j] = 0;
-      mul_cos_[j] = 0;
-    }
+    // for (long j = 0; j < N; ++j) {
+    //   mag_sin_[j] = 0;
+    //   mag_cos_[j] = 0;
+    //   mul_sin_[j] = 0;
+    //   mul_cos_[j] = 0;
+    // }
     for (unsigned long j = N-1; j < total_samples_; ++j) {
       mag_sin_[j] = 0;
       mag_cos_[j] = 0;

src/main.cc

@@ -13,8 +13,20 @@
 #include <string>
 
 const double PI_=3.14159265358979323846;
+const double eps_=1e-11;
+long total_generations_ = 5000;
+
+// int voxel_num = 10;
+// double vox[] = {0.652282,0.286431,0.727236,0.271936,0.526648,0.47129,0.631995,0.808049,0.322672,0.0329494,0.283877,0.930383,0.0681125,0.269938,0.205617,0.810485,0.615148,0.717859,0.846448,0.264384};
+
+int voxel_num = 5;
+double vox[] = {0.822628,0.691376,0.282906,0.226013,0.90703,0.144985,0.328563,0.440353,0.662462,0.720518};
+
+double total_periods = 20;
+int rf_samples_per_period = 1;
+
+std::vector<double> init_vox (vox, vox + sizeof(vox) / sizeof(double) );
 
-int voxel_num = 30;
 //double phase_range = 0;
 double phase_range = PI_/2;
 double phase_init = 0.0;
@@ -24,8 +36,6 @@ double noise_ratio = 1e8;
 
 // B0=1.5T freq=64Mhz, period = 15.6 ns
 double period_ms = 15.6/1000/1000; //ms
-double total_periods = 500;
-int rf_samples_per_period = 10;
 // double period_ms = 10; //ms
 // double total_periods = 2;
 // int rf_samples_per_period = 3;
@@ -34,47 +44,58 @@ int rf_samples_per_period = 10;
 //double T2s_min = T2s_scale/1000.0; //
 double T2s_scale = period_ms*total_periods; //ms # need to be 10ms
 template<class T> inline T pow2(const T value) {return value*value;}
+mri::GeneratorRF generator_rf_, fitness_rf_;
+jade::SubPopulation sub_population_;  // Optimizer.
+
+void InitGenerators();
+void PrintVector(const std::vector<double> &values);
+void SetOptimizer();
+
+double EvaluateFitness(const std::vector<double> &input){
+  return fitness_rf_.TryFit(input);
+};
 
-  
 int main(int argc, char *argv[]) {
   MPI_Init(&argc, &argv);
   int rank;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
   try {
-    //std::cout<<"Runing .... "<<std::endl;
-    jade::SubPopulation sub_population;  // Optimizer.
-    mri::GeneratorRF generator_rf;
-    generator_rf.SetVoxels(voxel_num, phase_init, phase_range, T2s_scale);
-    generator_rf.SetTimer(total_periods, period_ms, rf_samples_per_period);
-    generator_rf.RandomizeVoxels();
-    std::vector<double> signal_rf, magnitude_rf, phase_rf;
+    
+    InitGenerators();    
+    generator_rf_.SetVoxels(init_vox);
+    //generator_rf_.RandomizeVoxels(); // Set signal to fit
+    std::vector<double> origin_values;
+    generator_rf_.GetVoxels(origin_values);
+    PrintVector(origin_values);
     std::vector<double> mag_sin, mag_cos;
-    generator_rf.EvaluateSignal(signal_rf);
-    generator_rf.EvaluateDemodulated(mag_sin,mag_cos);
-    auto timesteps = generator_rf.GetTimesteps();
-    mri::LockIn lock_in;
-    std::vector<double> out_sin, out_cos;
-    lock_in.Analyse(signal_rf, period_ms, timesteps,
-                    rf_samples_per_period, out_sin, out_cos);
-    for (unsigned long i = 0; i<timesteps.size(); ++i) {
-      //      if (i > static_cast<unsigned long>(rf_samples_per_period)*2) 
-        std::cout << timesteps[i] << "\t"
-                  << signal_rf[i] << "\t"
-        // << mul_sin[i] << "\t"
-        // << mul_cos[i] << "\t"
-                  << std::sqrt(pow2(mag_sin[i])+pow2(mag_cos[i]))
-          /std::sqrt(pow2(out_sin[i])+pow2(out_cos[i])) << "\t"
-                  << std::sqrt(pow2(mag_sin[i])+pow2(mag_cos[i])) << "\t"
-                  << std::sqrt(pow2(out_sin[i])+pow2(out_cos[i])) << "\t"
-                  << out_sin[i] << "\t"
-                  << out_cos[i] << "\t"
-                  << mag_sin[i] << "\t"
-                  << mag_cos[i]
-                  << std::endl;
-    }
-    //std::cout<< timesteps.size() <<std::endl;
-    // for (auto value:signal_rf){
-    //   std::cout << value <<std::endl;
+    generator_rf_.EvaluateDemodulated(mag_sin,mag_cos);
+    fitness_rf_.SetBase(mag_sin, mag_cos);
+    
+    SetOptimizer();
+    sub_population_.FitnessFunction = &EvaluateFitness;
+
+    // double fit = fitness_rf_.TryFit(origin_values);
+    // std::cout << fit << std::endl;
+      
+    sub_population_.RunOptimization();
+    double fit;
+    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)
+      diff_vec[i] = std::abs(origin_values[i]-best_local_x[i])*100;
+    std::cout << "         diff:\t";
+    PrintVector(diff_vec);
+    std::cout << "fit:" << fit <<std::endl;
+    //auto timesteps = generator_rf_.GetTimesteps();
+    // for (unsigned long i = 0; i<timesteps.size(); ++i) {
+    //     std::cout << timesteps[i] << "\t"
+    //               << mag_sin[i] << "\t"
+    //               << mag_cos[i]
+    //               << std::endl;
     // }
     
   } catch( const std::invalid_argument& ia ) {
@@ -88,3 +109,40 @@ int main(int argc, char *argv[]) {
 // ********************************************************************** //
 // ********************************************************************** //
 // ********************************************************************** //
+void InitGenerators() {
+    generator_rf_.InitVoxels(voxel_num, phase_init, phase_range, T2s_scale);
+    fitness_rf_.InitVoxels(voxel_num, phase_init, phase_range, T2s_scale);
+    
+    generator_rf_.InitTimer(total_periods, period_ms, rf_samples_per_period);
+    fitness_rf_.InitTimer(total_periods, period_ms, rf_samples_per_period);
+};  
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+void PrintVector(const std::vector<double> &values) {
+  for (const double &value : values)
+    std::cout << value << " ";
+  std::cout << std::endl;
+};
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //
+void SetOptimizer() {
+  long dimension = 2*voxel_num;
+  long total_population = dimension * 20;
+  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_.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);
+
+}
+// ********************************************************************** //
+// ********************************************************************** //
+// ********************************************************************** //