tuning fixes and update comparing with Andrew

lib/Core/Plot.cpp

@@ -184,7 +184,7 @@ PlotData::PlotData(const DMatSample &x, const DVec &y, const bool abs)
     }
     else
     {
-        setCommand("'" + tmpFileName + "' u 1:(abs($3)):2 w xerr");
+        setCommand("'" + tmpFileName + "' u 1:($3):2 w xerr");
     }
 }
 
@@ -206,60 +206,6 @@ PlotData::PlotData(const XYStatData &data, const Index i, const Index j, const b
     setCommand("'" + tmpFileName + "' " + usingCmd);
 }
 
-// PlotPoint constructor ///////////////////////////////////////////////////////
-PlotPoint::PlotPoint(const double x, const double y)
-{
-    DMat d(1, 2);
-    string usingCmd, tmpFileName;
-
-    d(0, 0)     = x;
-    d(0, 1)     = y;
-    tmpFileName = dumpToTmpFile(d);
-    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' u 1:2");
-}
-
-PlotPoint::PlotPoint(const DSample &x, const double y)
-{
-    DMat d(1, 3);
-    string usingCmd, tmpFileName;
-
-    d(0, 0)     = x[central];
-    d(0, 2)     = y;
-    d(0, 1)     = sqrt(x.variance());
-    tmpFileName = dumpToTmpFile(d);
-    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' u 1:3:2 w xerr");
-}
-
-PlotPoint::PlotPoint(const double x, const DSample &y)
-{
-    DMat d(1, 3);
-    string usingCmd, tmpFileName;
-
-    d(0, 0)     = x;
-    d(0, 1)     = y[central];
-    d(0, 2)     = sqrt(y.variance());
-    tmpFileName = dumpToTmpFile(d);
-    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' u 1:2:3 w yerr");
-}
-
-PlotPoint::PlotPoint(const DSample &x, const DSample &y)
-{
-    DMat d(1, 4);
-    string usingCmd, tmpFileName;
-
-    d(0, 0)     = x[central];
-    d(0, 2)     = y[central];
-    d(0, 1)     = sqrt(x.variance());
-    d(0, 3)     = sqrt(y.variance());
-    tmpFileName = dumpToTmpFile(d);
-    pushTmpFile(tmpFileName);
-    setCommand("'" + tmpFileName + "' u 1:3:2:4 w xyerr");
-}
-
-
 // PlotLine constructor ////////////////////////////////////////////////////////
 PlotLine::PlotLine(const DVec &x, const DVec &y)
 {

lib/Core/Plot.hpp

@@ -98,18 +98,6 @@ public:
     virtual ~PlotData(void) = default;
 };
 
-class PlotPoint: public PlotObject
-{
-public:
-    // constructor
-    PlotPoint(const double  x,  const double  y);
-    PlotPoint(const DSample &x, const double  y);
-    PlotPoint(const double  x,  const DSample &y);
-    PlotPoint(const DSample &x, const DSample &y);
-    // destructor
-    virtual ~PlotPoint(void) = default;
-};
-
 class PlotHLine: public PlotObject
 {
 public:

lib/Core/Utilities.hpp

@@ -108,23 +108,6 @@ 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,39 +253,16 @@ DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
     }
 }
 
-DMatSample CorrelatorUtils::fold(const DMatSample &c, const CorrelatorModels::ModelPar &par)
+DMatSample CorrelatorUtils::fold(const DMatSample &c)
 {
     const Index nt  = c[central].rows();
     DMatSample  buf = c;
-    int         sign;
-    bool        fold = false;
 
-    switch (par.type)
-    {
-        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)
     {
         for (Index t = 0; t < nt; ++t)
         {
-                buf[s](t) = 0.5*(c[s](t) + sign*c[s]((nt - t) % nt));
+            buf[s](t) = 0.5*(c[s](t) + 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, const CorrelatorModels::ModelPar &par);
+    DMatSample fold(const DMatSample &c);
     DMatSample fourierTransform(const DMatSample &c, FFT &fft, 
                                 const unsigned int dir = FFT::Forward);
 };

lib/Statistics/Histogram.cpp

@@ -146,16 +146,6 @@ 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,8 +52,6 @@ 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,30 +300,36 @@ const XYStatData & XYSampleData::getData(void)
 }
 
 // fit /////////////////////////////////////////////////////////////////////////
-
+SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
-void XYSampleData::fitSample(std::vector<Minimizer *> &minimizer,
+                                  const DVec &init,
-                             const std::vector<const DoubleModel *> &v,
+                                  const std::vector<const DoubleModel *> &v)
-                             SampleFitResult &result,
-                             DVec &init,
-                             Index s)
 {
+    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());
-    FitResult sampleResult;
+    FOR_STAT_ARRAY(result, s)
+    {
         setDataToSample(s);
         if (s == central)
         {
-        sampleResult        = data_.fit(minimizer, init, v);
+            sampleResult = data_.fit(minimizer, initCopy, v);
-        init                = sampleResult.segment(0, init.size());
+            initCopy     = sampleResult.segment(0, initCopy.size());
-        result.nPar_        = sampleResult.getNPar();
+            if (verbCopy != Minimizer::Verbosity::Debug)
-        result.nDof_        = sampleResult.nDof_;
+            {
-        result.parName_     = sampleResult.parName_;
+                minimizer.back()->setVerbosity(Minimizer::Verbosity::Silent);
-        result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat());
+            }
         }
         else
         {
-        sampleResult = data_.fit(*(minimizer.back()), init, v);
+            
+            sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
         }
         result[s]       = sampleResult;
         result.chi2_[s] = sampleResult.getChi2();
@@ -332,50 +338,12 @@ void XYSampleData::fitSample(std::vector<Minimizer *> &minimizer,
             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)
-{
-    computeVarMat();
-    
-    SampleFitResult      result;
-    DVec                 initCopy = init;
-    Minimizer::Verbosity verbCopy = minimizer.back()->getVerbosity();
-    
-    FOR_STAT_ARRAY(result, s)
-    {
-        fitSample(minimizer, v, result, initCopy, s);
     }
     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,13 +103,6 @@ public:
     // get internal XYStatData
     const XYStatData & getData(void);
     // fit
-    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,

physics/2pt-fit.cpp

@@ -114,7 +114,6 @@ 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.};
@@ -156,11 +155,6 @@ 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/Makefile.am

@@ -14,6 +14,7 @@ bin_PROGRAMS =            \
     latan-sample-fake          \
     latan-sample-ft            \
     latan-sample-merge         \
+    latan-sample-noise-analysis\
     latan-sample-plot          \
     latan-sample-plot-corr     \
     latan-sample-read          \
@@ -47,6 +48,10 @@ latan_sample_merge_SOURCES  = sample-merge.cpp
 latan_sample_merge_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_merge_LDFLAGS  = -L../lib/.libs -lLatAnalyze
 
+latan_sample_noise_analysis_SOURCES  = sample-noise-analysis.cpp
+latan_sample_noise_analysis_CXXFLAGS = $(COM_CXXFLAGS)
+latan_sample_noise_analysis_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
 latan_sample_plot_corr_SOURCES  = sample-plot-corr.cpp
 latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_plot_corr_LDFLAGS  = -L../lib/.libs -lLatAnalyze

utils/sample-fake.cpp

@@ -18,42 +18,30 @@
  */
 
 #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;
 
-    opt.addOption("r", "seed"      , OptParser::OptType::value,   true,
+    if (argc != 5)
-                    "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;
-    SeedType              seed = (opt.gotOption("r")) ? opt.optionValue<SeedType>("r") : rd();
+    mt19937               gen(rd());
-    mt19937               gen(seed);
     normal_distribution<> dis(val, err);
     DSample               res(nSample);
 

utils/sample-noise-analysis.cpp

@@ -0,0 +1,193 @@
+/*
+ * sample-noise-analysis.cpp, part of LatAnalyze 3
+ *
+ * Copyright (C) 2013 - 2020 Antonin Portelli
+ *
+ * LatAnalyze 3 is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LatAnalyze 3 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with LatAnalyze 3.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <LatAnalyze/Core/OptParser.hpp>
+#include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Core/Math.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
+#include <LatAnalyze/Numerical/GslFFT.hpp>
+#include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
+#include <LatAnalyze/Statistics/XYSampleData.hpp>
+
+using namespace std;
+using namespace Latan;
+using namespace Math;
+
+int main(int argc, char *argv[])
+{
+    // argument parsing ////////////////////////////////////////////////////////
+    OptParser      opt;
+    bool           parsed;
+    string         filename;
+
+    opt.addOption("" , "help"  , OptParser::OptType::trigger, true,
+                  "show this help message and exit");
+    parsed = opt.parse(argc, argv);
+    if (!parsed or (opt.getArgs().size() != 1) or opt.gotOption("help"))
+    {
+        cerr << "usage: " << argv[0];
+        cerr << "<options> <sample file>" << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
+    }
+    filename = opt.getArgs()[0];
+
+    // load data ///////////////////////////////////////////////////////////////
+    DMatSample sample;
+
+    cout << "-- load data" << endl;
+    sample = Io::load<DMatSample>(filename);
+
+    // compute power spectrum //////////////////////////////////////////////////
+    DMat       av, err;
+    double     l0;
+    Index      nSample = sample.size(), n = sample[central].rows();
+    DMatSample noise(nSample), pow(nSample, n, 1);
+    CMatSample ftBuf(nSample, n, 1);
+    GslFFT     fft(n);
+
+    cout << "-- compute power spectrum" << endl;
+    FOR_STAT_ARRAY(sample, s)
+    {
+        sample[s].conservativeResize(n, 1);
+    }
+    av  = sample.mean();
+    err = sample.variance().cwiseSqrt();
+    FOR_STAT_ARRAY(sample, s)
+    {
+        noise[s] = sample[s] - av;
+        ftBuf[s].real() = noise[s];
+        ftBuf[s].imag().fill(0.);
+        fft(ftBuf[s]);
+        pow[s] = ftBuf[s].cwiseAbs2().unaryExpr([](const double x){return 10.*log10(x);});
+        //pow[s] = ftBuf[s].cwiseAbs2();
+        pow[s].conservativeResize(n/2, 1);
+    }
+    pow[central] = pow.mean();
+    // {
+    //     Plot p;
+    //     DVec x;
+
+    //     x.setLinSpaced(n/2, 0., n/2 - 1.);
+    //     p << LogScale(Axis::x);
+    //     p << PlotData(x, pow);
+    //     p.display();
+    // }
+    // l0 = pow.mean()(1);
+    // FOR_STAT_ARRAY(sample, s)
+    // {
+    //     pow[s] = pow[s].unaryExpr([l0](const double x){return x - l0;});
+    // }
+
+    // fit decay ///////////////////////////////////////////////////////////////
+    DVec            x, init(2);
+    DMat            fitErr;
+    DMatSample      xs(nSample, n/2, 1);
+    DSample         beta(nSample);
+    XYSampleData    data(nSample);
+    MinuitMinimizer min;
+    DoubleModel     lin([](const double *x, const double *p){return x[0]*p[0] + p[1];}, 1, 2);
+
+    cout << "-- fit decay" << endl;
+    x.setLinSpaced(n/2, 0., n/2 - 1.);
+    FOR_VEC(x, i)
+    {
+        x(i) = log2(x(i));
+    }
+    xs.fill(x);
+    data.addXDim(n/2, "f", true);
+    data.addYDim("pow");
+    data.setUnidimData(xs, pow);
+    data.assumeYYCorrelated(true, 0, 0);
+    for (unsigned int i = 0; i < n/2; ++i)
+    {
+        data.fitPoint((x(i) >= 2.) and (x(i) <= log2(n/2) - 0.5), i);
+    }
+    init(0) = -0.1; init(1) = -0.1;
+    auto fit = data.fit(min, init, lin);
+    fitErr   = fit.variance().cwiseSqrt();
+    FOR_STAT_ARRAY(beta, s)
+    {
+        beta[s] = fit[s](0)/(10.*log10(2.));
+    }
+    printf("chi^2/dof = %.1e/%d= %.2e -- chi^2 CCDF = %.2e -- p-value = %.2e -- CDR = %.1f dB\n",
+        fit.getChi2(), static_cast<int>(fit.getNDof()), fit.getChi2PerDof(), 
+        fit.getCcdf(), fit.getPValue(), fit.getCorrRangeDb());
+    printf("    decay = %.2f +/- %.2f dB/oct\n", fit[central](0), fitErr(0));
+    printf(" exponent = %.2f +/- %.2f\n", beta[central], sqrt(beta.variance()));
+
+    Plot p;
+
+    p << Caption("noise power spectrum");
+    p << PlotRange(Axis::x, -0.5, log2(n/2) + 0.5) 
+      << Label("frequency (oct)", Axis::x) << Label("power (dB)", Axis::y);
+    p << Color("1") << PlotPredBand(fit.getModel(_), 0., log2(n/2) + 0.5);
+    p << Color("1") << PlotFunction(fit.getModel(), 0., log2(n/2) + 0.5);
+    p << Color("2") << PlotData(x, pow);
+
+    p.display();
+
+    // p.reset();
+    // p << PlotCorrMatrix(data.getFitCorrMat());
+    // p.display();
+
+    // filter correlator ///////////////////////////////////////////////////////
+    DVec       filter(n);
+    DMatSample fsample(nSample, n, 1);
+
+    FOR_VEC(filter, i)
+    {
+        filter(i) = -std::pow(2.*sin(pi/n*i), 2);//-beta[central]*.5);
+    }
+    FOR_STAT_ARRAY(sample, s)
+    {
+        ftBuf[s].real() = sample[s].col(0);
+        ftBuf[s].imag().fill(0.);
+        fft(ftBuf[s], FFT::Forward);
+        ftBuf[s] = ftBuf[s].cwiseProduct(filter);
+        fft(ftBuf[s], FFT::Backward);
+        fsample[s] = ftBuf[s].real();
+    }
+
+    // p.reset();
+    x.setLinSpaced(n, 0., n - 1.);
+    // p << PlotData(x, sample);
+    // p << PlotData(x, fsample);
+    // p.display();
+
+    p.reset();
+    FOR_VEC(x, i)
+    {
+        x(i) = log2(x(i));
+    }
+    p << PlotRange(Axis::x, -0.5, log2(n/2) + 0.5);
+    p << PlotPoints(x, -filter);
+    p.display();
+    p.reset();
+    p << PlotCorrMatrix(sample.correlationMatrix());
+    p.display();
+    p.reset();
+    p << PlotCorrMatrix(fsample.correlationMatrix());
+    p.display();
+    Io::save(fsample, "test.h5");
+
+
+    return EXIT_SUCCESS;
+}

utils/sample-read.cpp

@@ -38,23 +38,9 @@ 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 +/-  standard deviation\n" << endl;
+        cout << "central value:\n"      << s[central]               << endl;
-        cout << "Re:" << endl;
+        cout << "standard deviation:\n" << s.variance().cwiseSqrt() << 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);
@@ -65,8 +51,8 @@ int main(int argc, char *argv[])
         DSample s    = Io::load<DSample>(fileName);
         string  name = Io::getFirstName(fileName);
         cout << scientific;
-        cout << "central value +/-  standard deviation\n" << endl;
+        cout << "central value:\n"      << s[central]         << endl;
-        cout << s[central] << " +/- " << sqrt(s.variance()) << endl;
+        cout << "standard deviation:\n" << sqrt(s.variance()) << endl;
         if (!copy.empty())
         {
             Io::save(s, copy, File::Mode::write, name);