Merge pull request #19 from AndrewYongZhenNing/develop

Corrected Sinh model's definition, init fit parameter & fold definition.

lib/Core/Utilities.hpp

@@ -108,6 +108,23 @@ inline std::string strFrom(const T x)
 }
 
 // specialization for vectors
+template<>
+inline std::vector<Index> strTo<std::vector<Index>>(const std::string &str)
+{
+    std::vector<Index>  res;
+    std::vector<double> vbuf;
+    double              buf;
+    std::istringstream  stream(str);
+    
+    while (!stream.eof())
+    {
+        stream >> buf;
+        res.push_back(buf);
+    }
+    
+    return res;
+}
+
 template<>
 inline DVec strTo<DVec>(const std::string &str)
 {

lib/Physics/CorrelatorFitter.cpp

@@ -253,16 +253,39 @@ DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
     }
 }
 
-DMatSample CorrelatorUtils::fold(const DMatSample &c)
+DMatSample CorrelatorUtils::fold(const DMatSample &c, const CorrelatorModels::ModelPar &par)
 {
     const Index nt  = c[central].rows();
     DMatSample  buf = c;
-
+    int         sign;
-    FOR_STAT_ARRAY(buf, s)
+    bool        fold = false;
+    
+    switch (par.type)
     {
-        for (Index t = 0; t < nt; ++t)
+        case CorrelatorType::cosh:
+        case CorrelatorType::cst:
+            sign = 1;
+            fold = true;
+            break;
+        case CorrelatorType::sinh:
+            sign = -1;
+            fold = true;
+            break;
+        case CorrelatorType::linear:
+            cout << "Linear model is asymmetric: will not fold." << endl;
+            break;
+        default:
+            break;
+    }
+
+    if (fold)
+    {
+        FOR_STAT_ARRAY(buf, s)
         {
-            buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt));
+            for (Index t = 0; t < nt; ++t)
+            {
+                buf[s](t) = 0.5*(c[s](t) + sign*c[s]((nt - t) % nt));
+            }
         }
     }
 

lib/Physics/CorrelatorFitter.hpp

@@ -56,7 +56,7 @@ namespace CorrelatorModels
 namespace CorrelatorUtils
 {
     DMatSample shift(const DMatSample &c, const Index ts);
-    DMatSample fold(const DMatSample &c);
+    DMatSample fold(const DMatSample &c, const CorrelatorModels::ModelPar &par);
     DMatSample fourierTransform(const DMatSample &c, FFT &fft, 
                                 const unsigned int dir = FFT::Forward);
 };

lib/Statistics/Histogram.cpp

@@ -146,6 +146,16 @@ double Histogram::getX(const Index i) const
     return x_(i);
 }
 
+double Histogram::getXMin(void) const
+{
+    return xMin_;
+}
+
+double Histogram::getXMax(void) const
+{
+    return xMax_;
+}
+
 double Histogram::operator[](const Index i) const
 {
     return bin_(i)*(isNormalized() ? norm_ : 1.);

lib/Statistics/Histogram.hpp

@@ -52,6 +52,8 @@ public:
     const StatArray<double> & getData(void) const;
     const StatArray<double> & getWeight(void) const;
     double                    getX(const Index i) const;
+    double                    getXMin(void) const;
+    double                    getXMax(void) const;
     double                    operator[](const Index i) const;
     double                    operator()(const double x) const;
     // percentiles & confidence interval

lib/Statistics/XYSampleData.cpp

@@ -300,6 +300,67 @@ const XYStatData & XYSampleData::getData(void)
 }
 
 // fit /////////////////////////////////////////////////////////////////////////
+
+void XYSampleData::fitSample(std::vector<Minimizer *> &minimizer,
+                             const std::vector<const DoubleModel *> &v,
+                             SampleFitResult &result,
+                             DVec &init,
+                             Index s)
+{
+    result.resize(nSample_);
+    result.chi2_.resize(nSample_);
+    result.model_.resize(v.size());
+    FitResult sampleResult;
+    setDataToSample(s);
+    if (s == central)
+    {
+        sampleResult        = data_.fit(minimizer, init, v);
+        init                = sampleResult.segment(0, init.size());
+        result.nPar_        = sampleResult.getNPar();
+        result.nDof_        = sampleResult.nDof_;
+        result.parName_     = sampleResult.parName_;
+        result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
+    }
+    else
+    {
+        sampleResult = data_.fit(*(minimizer.back()), init, v);
+    }
+    result[s]       = sampleResult;
+    result.chi2_[s] = sampleResult.getChi2();
+    for (unsigned int j = 0; j < v.size(); ++j)
+    {
+        result.model_[j].resize(nSample_);
+        result.model_[j][s] = sampleResult.getModel(j);
+    }
+
+    
+}
+
+SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
+                                  const DVec &init,
+                                  const std::vector<const DoubleModel *> &v,
+                                  Index s)
+{
+    computeVarMat();
+    
+    SampleFitResult result;
+    DVec      initCopy = init;
+
+    fitSample(minimizer, v, result, initCopy, s);
+    
+    return result;
+}
+
+SampleFitResult XYSampleData::fit(Minimizer &minimizer,
+                                  const DVec &init,
+                                  const std::vector<const DoubleModel *> &v,
+                                  Index s)
+{
+    vector<Minimizer *> mv{&minimizer};
+    
+    return fit(mv, init, v, s);
+}
+
 SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
                                   const DVec &init,
                                   const std::vector<const DoubleModel *> &v)
@@ -307,43 +368,14 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
     computeVarMat();
     
     SampleFitResult      result;
-    FitResult            sampleResult;
     DVec                 initCopy = init;
     Minimizer::Verbosity verbCopy = minimizer.back()->getVerbosity();
     
-    result.resize(nSample_);
-    result.chi2_.resize(nSample_);
-    result.model_.resize(v.size());
     FOR_STAT_ARRAY(result, s)
     {
-        setDataToSample(s);
+        fitSample(minimizer, v, result, initCopy, s);
-        if (s == central)
-        {
-            sampleResult = data_.fit(minimizer, initCopy, v);
-            initCopy     = sampleResult.segment(0, initCopy.size());
-            if (verbCopy != Minimizer::Verbosity::Debug)
-            {
-                minimizer.back()->setVerbosity(Minimizer::Verbosity::Silent);
-            }
-        }
-        else
-        {
-            
-            sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
-        }
-        result[s]       = sampleResult;
-        result.chi2_[s] = sampleResult.getChi2();
-        for (unsigned int j = 0; j < v.size(); ++j)
-        {
-            result.model_[j].resize(nSample_);
-            result.model_[j][s] = sampleResult.getModel(j);
-        }
     }
     minimizer.back()->setVerbosity(verbCopy);
-    result.nPar_       = sampleResult.getNPar();
-    result.nDof_       = sampleResult.nDof_;
-    result.parName_    = sampleResult.parName_;
-    result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
     
     return result;
 }

lib/Statistics/XYSampleData.hpp

@@ -103,9 +103,16 @@ public:
     // get internal XYStatData
     const XYStatData & getData(void);
     // fit
-    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
+    void            fitSample(std::vector<Minimizer *> &minimizer, 
+                              const std::vector<const DoubleModel *> &v, 
+                              SampleFitResult &sampleResult, DVec &init, Index s);
+    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init, 
+                        const std::vector<const DoubleModel *> &v, Index s);
+    SampleFitResult fit(Minimizer &minimizer, const DVec &init, 
+                        const std::vector<const DoubleModel *> &v, Index s);
+    SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init, 
                         const std::vector<const DoubleModel *> &v);
-    SampleFitResult fit(Minimizer &minimizer, const DVec &init,
+    SampleFitResult fit(Minimizer &minimizer, const DVec &init, 
                         const std::vector<const DoubleModel *> &v);
     template <typename... Ts>
     SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,

physics/2pt-fit.cpp

@@ -114,6 +114,7 @@ int main(int argc, char *argv[])
     nt      = corr[central].rows();
     corr    = corr.block(0, 0, nt, 1);
     corr    = CorrelatorUtils::shift(corr, shift);
+
     if (doLaplace)
     {
         vector<double> filter = {1., -2., 1.};
@@ -155,6 +156,11 @@ int main(int argc, char *argv[])
         }
     }
     
+    if (fold)
+    {
+        corr = CorrelatorUtils::fold(corr,modelPar);
+    }
+
     // fit /////////////////////////////////////////////////////////////////////
     DVec                init(nPar);
     NloptMinimizer      globMin(NloptMinimizer::Algorithm::GN_CRS2_LM);

utils/sample-fake.cpp

@@ -18,30 +18,42 @@
  */
 
 #include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Core/OptParser.hpp>
+
 
 using namespace std;
 using namespace Latan;
 
 int main(int argc, char *argv[])
 {
+    OptParser            opt;
     Index  nSample;
     double val, err;
     string outFileName;
 
-    if (argc != 5)
+    opt.addOption("r", "seed"      , OptParser::OptType::value,   true,
+                    "random generator seed (default: random)"); 
+    opt.addOption("", "help"      , OptParser::OptType::trigger, true,
+                  "show this help message and exit");
+
+    bool parsed = opt.parse(argc, argv);
+    if (!parsed or (opt.getArgs().size() != 4) or opt.gotOption("help"))
     {
         cerr << "usage: " << argv[0];
         cerr << " <central value> <error> <nSample> <output file>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
 
         return EXIT_FAILURE;
     }
+
     val         = strTo<double>(argv[1]);
     err         = strTo<double>(argv[2]);
     nSample     = strTo<Index>(argv[3]);
     outFileName = argv[4];
 
     random_device         rd;
-    mt19937               gen(rd());
+    SeedType              seed = (opt.gotOption("r")) ? opt.optionValue<SeedType>("r") : rd();
+    mt19937               gen(seed);
     normal_distribution<> dis(val, err);
     DSample               res(nSample);
 
@@ -59,4 +71,4 @@ int main(int argc, char *argv[])
     Io::save<DSample>(res, outFileName);
 
     return EXIT_SUCCESS;
-}
+}

utils/sample-read.cpp

@@ -38,9 +38,23 @@ int main(int argc, char *argv[])
     {
         DMatSample s    = Io::load<DMatSample>(fileName);
         string     name = Io::getFirstName(fileName);
+        Index nRows     = s[central].rows();
+        Index nCols     = s[central].cols();
         cout << scientific;
-        cout << "central value:\n"      << s[central]               << endl;
+        cout << "central value +/-  standard deviation\n" << endl;
-        cout << "standard deviation:\n" << s.variance().cwiseSqrt() << endl;
+        cout << "Re:" << endl;
+        for(Index i = 0; i < nRows; i++)
+        {
+            cout << s[central](i,0) << " +/- " << s.variance().cwiseSqrt()(i,0) << endl;
+        }
+        if(nCols == 2)
+        {
+            cout << "\nIm:" << endl;
+            for(Index i = 0; i < nRows; i++)
+            {
+                cout << s[central](i,1) << " +/- " << s.variance().cwiseSqrt()(i,1) << endl;
+            }   
+        }
         if (!copy.empty())
         {
             Io::save(s, copy, File::Mode::write, name);
@@ -51,8 +65,8 @@ int main(int argc, char *argv[])
         DSample s    = Io::load<DSample>(fileName);
         string  name = Io::getFirstName(fileName);
         cout << scientific;
-        cout << "central value:\n"      << s[central]         << endl;
+        cout << "central value +/-  standard deviation\n" << endl;
-        cout << "standard deviation:\n" << sqrt(s.variance()) << endl;
+        cout << s[central] << " +/- " << sqrt(s.variance()) << endl;
         if (!copy.empty())
         {
             Io::save(s, copy, File::Mode::write, name);