Merge remote-tracking branch 'fork/develop' into develop

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/XYSampleData.cpp

@@ -234,6 +234,21 @@ DVec XYSampleData::getYError(const Index j)
     return data_.getYError(j);
 }
 
+bool XYSampleData::checkFit()
+{
+    return goodFit_;
+}
+
+void XYSampleData::checkChi2PerDof(double Chi2PerDof)
+{
+    if(Chi2PerDof >= 2 or Chi2PerDof < 0 or isnan(Chi2PerDof)) 
+    {
+        goodFit_ = false;
+        cerr << "chi2PerDof = " << Chi2PerDof << ". Aborting fit now." << endl;
+    }
+
+}
+
 // get total fit variance matrix and its pseudo-inverse ////////////////////////
 const DMat & XYSampleData::getFitVarMat(void)
 {
@@ -292,24 +307,34 @@ SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
     result.resize(nSample_);
     result.chi2_.resize(nSample_);
     result.model_.resize(v.size());
+    double chi2PerDof;
+    goodFit_ = true;
     FOR_STAT_ARRAY(result, s)
     {
-        setDataToSample(s);
+        if(goodFit_)
-        if (s == central)
         {
-            sampleResult = data_.fit(minimizer, initCopy, v);
+            setDataToSample(s);
-            initCopy     = sampleResult.segment(0, initCopy.size());
+            if (s == central)
-        }
+            {
-        else
+                sampleResult = data_.fit(minimizer, initCopy, v);
-        {
+                initCopy     = sampleResult.segment(0, initCopy.size());
-            sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
+                chi2PerDof   = sampleResult.getChi2PerDof();
-        }
+                checkChi2PerDof(chi2PerDof);
-        result[s]       = sampleResult;
+            }
-        result.chi2_[s] = sampleResult.getChi2();
+            else
-        for (unsigned int j = 0; j < v.size(); ++j)
+            {
-        {
+                sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
-            result.model_[j].resize(nSample_);
+                chi2PerDof   = sampleResult.getChi2PerDof();
-            result.model_[j][s] = sampleResult.getModel(j);
+                checkChi2PerDof(chi2PerDof);
+            }
+            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);
+
+            }
         }
     }
     result.nPar_    = sampleResult.getNPar();

lib/Statistics/XYSampleData.hpp

@@ -91,6 +91,8 @@ public:
     const DMat &       getXYVar(const Index i, const Index j);
     DVec               getXError(const Index i);
     DVec               getYError(const Index j);
+    bool               checkFit(); // check fit candidate based on chi2PerDof
+    void               checkChi2PerDof(double Chi2PerDof);
     // get total fit variance matrix and its pseudo-inverse
     const DMat & getFitVarMat(void);
     const DMat & getFitVarMatPInv(void);
@@ -133,6 +135,7 @@ private:
     Index                                 nSample_, dataSample_{central};
     bool                                  initData_{true}, computeVarMat_{true};
     bool                                  initXMap_{true};
+    bool                                  goodFit_{true}; // used to break minimisation if central sample chi2PerDof is bad
 };
 
 /******************************************************************************

physics/2pt-fit.cpp

@@ -110,10 +110,6 @@ int main(int argc, char *argv[])
     nt      = corr[central].rows();
     corr    = corr.block(0, 0, nt, 1);
     corr    = CorrelatorUtils::shift(corr, shift);
-    if (fold)
-    {
-        corr = CorrelatorUtils::fold(corr);
-    }
     
     // make model //////////////////////////////////////////////////////////////
     CorrelatorFitter fitter(corr);
@@ -140,6 +136,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;
-}
+}