Merge pull request #40 from pbakalov/master

bug fix: setting sigma_ from covariance matrix eigenvalues
CQMP · May 27, 2018 · f5e50d3 · f5e50d3
2 parents d0108da + affb036
commit f5e50d3
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 7 deletions.
diff --git a/src/maxent_params.cpp b/src/maxent_params.cpp
@@ -238,9 +238,11 @@ void ContiParameters::decompose_covariance_matrix(const alps::params& p){
     // We drop eigenvalues of the covariance matrix which are smaller than 1e-10
     // as they represent bad data directions (usually there is a steep drop
     // below that value)
+    double cutoff = p["CM_EIGENVALUE_CUTOFF"];
+    std::cout << "Cutoff for singular eigenvalues of the covariance: " << cutoff << "\n";
     int new_ndat_,old_ndat_=ndat();
     for (new_ndat_ =0;new_ndat_<ndat();new_ndat_++)
-      if (var[new_ndat_]>1e-10) break;
+      if (var[new_ndat_]>cutoff) break; 
     // This is the number of good data
     ndat_ = old_ndat_ - new_ndat_;
     if (new_ndat_ ==0)
@@ -259,7 +261,7 @@ void ContiParameters::decompose_covariance_matrix(const alps::params& p){
       }
     }
     for (int i=0; i<ndat(); ++i) {
-      sigma_[i] = std::abs(var(new_ndat_+i))/static_cast<double>(p["NORM"]);
+      sigma_[i] = std::sqrt(std::abs(var(new_ndat_+i)))/static_cast<double>(p["NORM"]);
       if (p["VERBOSE"])
         std::cout << "# " << var(new_ndat_+i) << "\n";
     }

diff --git a/src/maxent_simulation.cpp b/src/maxent_simulation.cpp
@@ -57,6 +57,7 @@ void MaxEntSimulation::define_parameters(alps::params &p){
   p.define<int>("NDAT","# of input points");
   p.define<std::string>("DATA","","data file input");
   p.define<std::string>("COVARIANCE_MATRIX","","name of covariance matrix file");
+  p.define<double>("CM_EIGENVALUE_CUTOFF",1e-10, "Cutoff for eigenvalues of the covariance matrix");
   p.define<std::string>("BASENAME","","Specified output name (generated if not given)");
   p.define<int>("MODEL_RUNS","How many default model runs");
   p.define<double>("X_0","G input for param file entry");
@@ -309,15 +310,20 @@ void MaxEntSimulation::evaluate(){
     // for the self energy: use Im Sigma(omega)=-A(omega)*pi
     ofstream_ maxspec_self_str;maxspec_self_str.open((name+"maxspec_self.dat").c_str());
     ofstream_ avspec_self_str; avspec_self_str.open((name+"avspec_self.dat").c_str());
+    ofstream_ chispec_self_str; chispec_self_str.open((name+"chispec_self.dat").c_str());
     for (std::size_t  i=0; i<avspec.size(); ++i){ 
       avspec_self_str << omega_coord(i) << " " << -avspec[i]*M_PI<< " " << -def[i]*norm*M_PI<<std::endl;
     }
     for (std::size_t i=0; i<spectra[0].size(); ++i){
       maxspec_self_str << omega_coord(i) << " " << -spectra[max_a][i]*norm*M_PI<< " " << -def[i]*norm*M_PI << std::endl;
     }
+    for (std::size_t  i=0; i<specchi.size(); ++i){ 
+      chispec_self_str << omega_coord(i) << " " << -specchi[i]*M_PI<< " " << -def[i]*norm*M_PI<<std::endl;
+    }
     //for public facing variables
     avspec*=-M_PI;
     maxspec*=-M_PI;
+    specchi*=-M_PI;
   }
 
   if(make_back){
@@ -337,7 +343,7 @@ void MaxEntSimulation::evaluate(){
 
     std::cerr << "spectra"<<sp<< " max backcont diff" <<sp<<  "chi^2 value " <<std::endl;
     std::cerr << "======="<< sp<<" ================="<< sp<< "=========== " <<std::endl;
-    backcontinue(chispec_back_file,specchi,norm,"chispec",chispec_back);
+    backcontinue(chispec_back_file,specchi,norm_fix,"chispec",chispec_back);
     backcontinue(avspec_back_file,avspec,norm_fix,"avspec ",avspec_back);
     backcontinue(maxspec_back_file,maxspec,norm_fix,"maxspec",maxspec_back);
   }

diff --git a/test/paramsTest.cpp b/test/paramsTest.cpp
@@ -328,8 +328,8 @@ TEST(Parameters,CovarianceDataInFile){
   EXPECT_EQ(c.T(),0.5);
 
   for(int i=0;i<c.ndat();i++){
-    EXPECT_NEAR(c.y(i),(i+1)*.1/0.3,1e-10);
-    EXPECT_EQ(c.sigma(i),0.3);
+    EXPECT_NEAR(c.y(i),(i+1)*.1/std::sqrt(0.3),1e-10);
+    EXPECT_EQ(c.sigma(i),std::sqrt(0.3));
   }
 
   std::remove(pf.c_str());
@@ -392,8 +392,8 @@ TEST(Parameters,CovarianceHDF5Params){
 
 	//check input values
   for(int i=0;i<c.ndat();i++){
-    EXPECT_NEAR(c.y(i),(i+1)*.1/0.3,1e-10);
-    EXPECT_EQ(c.sigma(i),0.3);
+    EXPECT_NEAR(c.y(i),(i+1)*.1/std::sqrt(0.3),1e-10);
+    EXPECT_EQ(c.sigma(i),std::sqrt(0.3));
   }
   std::remove(tf.c_str());
 }