Merge pull request #14 from AndrewYongZhenNing/develop

Develop

physics/2pt-fit.cpp

@@ -16,7 +16,7 @@ int main(int argc, char *argv[])
     // parse arguments /////////////////////////////////////////////////////////
     OptParser            opt;
     bool                 parsed, doPlot, doHeatmap, doCorr, fold;
-    string               corrFileName, model, outFileName, outFmt;
+    string               corrFileName, model, outFileName, outFmt, savePlot;
     Index                ti, tf, shift, nPar, thinning;
     double               svdTol;
     Minimizer::Verbosity verbosity;
@@ -30,7 +30,7 @@ int main(int argc, char *argv[])
     opt.addOption("s", "shift"    , OptParser::OptType::value  , true,
                   "time variable shift", "0");
     opt.addOption("m", "model"    , OptParser::OptType::value  , true,
-                  "fit model (exp|exp2|exp3|cosh|cosh2|cosh3|explin|<interpreter code>)", "cosh");
+                  "fit model (exp|exp2|exp3|sinh|cosh|cosh2|cosh3|explin|const|<interpreter code>)", "cosh");
     opt.addOption("" , "nPar"     , OptParser::OptType::value  , true,
                   "number of model parameters for custom models "
                   "(-1 if irrelevant)", "-1");
@@ -48,6 +48,8 @@ int main(int argc, char *argv[])
                   "show the fit plot");
     opt.addOption("h", "heatmap"  , OptParser::OptType::trigger, true,
                   "show the fit correlation heatmap");
+    opt.addOption("s", "save-plot", OptParser::OptType::value, true,
+                    "saves the source and .pdf", "");
     opt.addOption("", "help"      , OptParser::OptType::trigger, true,
                   "show this help message and exit");
     parsed = opt.parse(argc, argv);
@@ -71,6 +73,7 @@ int main(int argc, char *argv[])
     fold         = opt.gotOption("fold");
     doPlot       = opt.gotOption("p");
     doHeatmap    = opt.gotOption("h");
+    savePlot     = opt.optionValue("s");
     switch (opt.optionValue<unsigned int>("v"))
     {
         case 0:
@@ -117,7 +120,7 @@ int main(int argc, char *argv[])
     
     // make model //////////////////////////////////////////////////////////////
     DoubleModel mod;
-    bool        coshModel = false, linearModel = false;
+    bool        sinhModel = false, coshModel = false, linearModel = false, constModel = false;
     
     if ((model == "exp") or (model == "exp1"))
     {
@@ -144,6 +147,15 @@ int main(int argc, char *argv[])
                                    + p[5]*exp(-p[4]*x[0]);
                         }, 1, nPar);
     }
+    else if (model == "sinh")
+    {
+        sinhModel = true;
+        nPar      = 2;
+        mod.setFunction([nt](const double *x, const double *p)
+                        {
+                            return p[1]*(exp(-p[0]*x[0])-exp(-p[0]*(nt-x[0])));
+                        }, 1, nPar);
+    }
     else if ((model == "cosh") or (model == "cosh1"))
     {
         coshModel = true;
@@ -183,6 +195,15 @@ int main(int argc, char *argv[])
                             return p[1] - p[0]*x[0];
                         }, 1, nPar);
     }
+    else if (model == "const")
+    {
+        constModel = true;
+        nPar       = 1;
+        mod.setFunction([](const double *x __dumb, const double *p)
+                        {
+                            return p[0];
+                        }, 1, nPar);
+    }
     else
     {
         if (nPar > 0)
@@ -214,26 +235,43 @@ int main(int argc, char *argv[])
     data.addXDim(nt, "t/a", true);
     data.addYDim("C(t)");
     data.setUnidimData(tvec, corr);
-    for (Index p = 0; p < nPar; p += 2)
+    // set parameter name /////////////
+    if(constModel)
     {
-        mod.parName().setName(p, "E_" + strFrom(p/2));
+        mod.parName().setName(0, "const");
-        mod.parName().setName(p + 1, "Z_" + strFrom(p/2));
     }
+    else
+    {
+        for (Index p = 0; p < nPar; p += 2)
+        {
+            mod.parName().setName(p, "E_" + strFrom(p/2));
+            mod.parName().setName(p + 1, "Z_" + strFrom(p/2));   
+        }
+    }
+    //set initial values ////////////////
     if (linearModel)
     {
         init(0) = data.y(nt/4, 0)[central] - data.y(nt/4 + 1, 0)[central];
         init(1) = data.y(nt/4, 0)[central] + nt/4*init(0);
     }
+    else if(constModel)
+    {
+        init(0) = data.y(nt/4, 0)[central];
+
+    }
     else
     {
         init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]);
         init(1) = data.y(nt/4, 0)[central]/(exp(-init(0)*nt/4));
+        // cout << init(0) << "\t" << init(1) << endl;
     }
     for (Index p = 2; p < nPar; p += 2)
     {
         init(p)     = 2*init(p - 2);
         init(p + 1) = init(p - 1)/2.;
+        
     }
+    // set limits for minimiser //////////////
     for (Index p = 0; p < nPar; p += 2)
     {
         if (linearModel)
@@ -241,20 +279,32 @@ int main(int argc, char *argv[])
             globMin.setLowLimit(p, -10.*fabs(init(p)));
             globMin.setHighLimit(p, 10.*fabs(init(p)));
         }
+        else if(constModel)
+        {
+            globMin.setLowLimit(p, -10*fabs(init(0)));
+            locMin.setLowLimit(p, -10*fabs(init(0)));
+            globMin.setHighLimit(p, 10*fabs(init(0)));
+        }
         else
         {
             globMin.setLowLimit(p, 0.);
-            locMin.setLowLimit(p, 0.);
+            // locMin.setLowLimit(p, 0.);
             globMin.setHighLimit(p, 10.*init(p));
         }
-        globMin.setLowLimit(p + 1, -10.*fabs(init(p + 1)));
+        if(!constModel)
-        globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1)));
+        {
+            locMin.setLowLimit(p+1, -1);
+            globMin.setLowLimit(p + 1, -10.*fabs(init(p + 1)));
+            globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1)));
+        }
+        
     }
     globMin.setPrecision(0.001);
     globMin.setMaxIteration(100000);
     globMin.setVerbosity(verbosity);
     locMin.setMaxIteration(1000000);
     locMin.setVerbosity(verbosity);
+    // fit /////////////////////////////////
     for (Index t = 0; t < nt; ++t)
     {
         data.fitPoint((t >= ti) and (t <= tf)
@@ -278,69 +328,86 @@ int main(int argc, char *argv[])
         fit = data.fit(locMin, init, mod);
         fit.print();
     }
-    
     // plots ///////////////////////////////////////////////////////////////////
     if (doPlot)
     {
         Plot       p;
-        DMatSample effMass(nSample);
-        DVec       effMassT, fitErr;
-        Index      maxT = (coshModel) ? (nt - 2) : (nt - 1);
-        double     e0, e0Err;
 
         p << PlotRange(Axis::x, 0, nt - 1);
-        if (!linearModel)
+        if (!linearModel and !constModel)
         {
             p << LogScale(Axis::y);
         }
         p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1);
         p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1);
-        p << Color("rgb 'red'") << PlotData(data.getData());
+        p << Color("rgb 'red'")  << PlotData(data.getData());
         p.display();
-        effMass.resizeMat(maxT, 1);
+        if(savePlot != "")
-        effMassT.setLinSpaced(maxT, 1, maxT);
-        fitErr = fit.variance().cwiseSqrt();
-        e0     = fit[central](0);
-        e0Err  = fitErr(0);
-        if (coshModel)
         {
-            FOR_STAT_ARRAY(effMass, s)
+            cout << "Saving plot and source code to " << savePlot << endl;
+            p.save(savePlot);
+        }
+        // effective mass plot //////////////////////////////////////////////////////
+        if (!constModel)
+        {
+            DMatSample effMass(nSample);
+            DVec       effMassT, fitErr;
+            Index      maxT = (coshModel) ? (nt - 2) : (nt - 1);
+            double     e0, e0Err;
+
+            effMass.resizeMat(maxT, 1);
+            effMassT.setLinSpaced(maxT, 0, maxT-1);
+            fitErr = fit.variance().cwiseSqrt();
+            e0     = fit[central](0);
+            e0Err  = fitErr(0);
+            if (coshModel or sinhModel)
             {
-                for (Index t = 1; t < nt - 1; ++t)
+                FOR_STAT_ARRAY(effMass, s)
                 {
-                    effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1))
+                    for (Index t = 1; t < nt - 1; ++t)
-                                              /(2.*corr[s](t)));
+                    {
+                        effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1))
+                                                /(2.*corr[s](t)));
+                    }
                 }
             }
-        }
+            else if (linearModel)
-        else if (linearModel)
-        {
-            FOR_STAT_ARRAY(effMass, s)
             {
-                for (Index t = 0; t < nt - 1; ++t)
+                FOR_STAT_ARRAY(effMass, s)
                 {
-                    effMass[s](t) = corr[s](t) - corr[s](t+1);
+                    for (Index t = 0; t < nt - 1; ++t)
+                    {
+                        effMass[s](t) = corr[s](t) - corr[s](t+1);
+                    }
                 }
             }
-        }
+            else
-        else
-        {
-            FOR_STAT_ARRAY(effMass, s)
             {
-                for (Index t = 1; t < nt; ++t)
+                FOR_STAT_ARRAY(effMass, s)
                 {
-                    effMass[s](t - 1) = log(corr[s](t-1)/corr[s](t));
+                    for (Index t = 1; t < nt; ++t)
+                    {
+                        effMass[s](t - 1) = log(corr[s](t-1)/corr[s](t));
+                    }
                 }
             }
+            p.reset();
+            p << PlotRange(Axis::x, 0, maxT);
+            p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err);
+            p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err);
+            p << Color("rgb 'blue'") << PlotHLine(e0);
+            p << Color("rgb 'red'") << PlotData(effMassT, effMass);
+            p << Caption("Effective Mass");
+            p.display();
+            if(savePlot != "")
+            {
+                string savename = savePlot + "_effMass";
+                cout << "Saving effective mass plot and source code to " << savename << endl;
+            p.save(savename);
+            }
         }
-        p.reset();
-        p << PlotRange(Axis::x, 1, maxT);
-        p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err);
-        p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err);
-        p << Color("rgb 'blue'") << PlotHLine(e0);
-        p << Color("rgb 'red'") << PlotData(effMassT, effMass);
-        p.display();
     }         
+ 
     if (doHeatmap)
     {
         Plot  p;
@@ -359,6 +426,7 @@ int main(int argc, char *argv[])
         }
     }
     
+    
     // output //////////////////////////////////////////////////////////////////
     if (!outFileName.empty())
     {

utils/Makefile.am

@@ -7,6 +7,7 @@ endif
 endif
 
 bin_PROGRAMS =            \
+    latan-plot            \
     latan-sample-combine  \
     latan-sample-element  \
     latan-sample-fake     \
@@ -16,6 +17,10 @@ bin_PROGRAMS =            \
     latan-sample-read     \
     latan-resample        
 
+latan_plot_SOURCES  = plot.cpp
+latan_plot_CXXFLAGS = $(COM_CXXFLAGS)
+latan_plot_LDFLAGS  = -L../lib/.libs -lLatAnalyze
+
 latan_sample_combine_SOURCES  = sample-combine.cpp
 latan_sample_combine_CXXFLAGS = $(COM_CXXFLAGS)
 latan_sample_combine_LDFLAGS  = -L../lib/.libs -lLatAnalyze

utils/plot.cpp

@@ -0,0 +1,142 @@
+#include <LatAnalyze/Core/OptParser.hpp>
+#include <LatAnalyze/Functional/CompiledModel.hpp>
+#include <LatAnalyze/Io/Io.hpp>
+#include <LatAnalyze/Statistics/MatSample.hpp>
+#include <LatAnalyze/Core/Math.hpp>
+#include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
+#include <LatAnalyze/Numerical/NloptMinimizer.hpp>
+#include <LatAnalyze/Core/Plot.hpp>
+#include <LatAnalyze/Statistics/XYSampleData.hpp>
+
+using namespace std;
+using namespace Latan;
+
+int main(int argc, char *argv[])
+{
+    // parse arguments /////////////////////////////////////////////////////////
+    OptParser            opt;
+    bool                 parsed, imag;
+    string               plotFileName, outFileName, xName, yName, title, save;
+    vector<string>       inFileName;
+    double               xLow, xHigh, spacing;
+
+    opt.addOption("i" , "imag"   , OptParser::OptType::trigger, true,
+                "plot imaginary");
+    opt.addOption("o", "output", OptParser::OptType::value  , true,
+                  "output file", "");
+    opt.addOption("x", "xAxis", OptParser::OptType::value  , true,
+                  "x-axis name", "");
+    opt.addOption("l", "xLow", OptParser::OptType::value  , true,
+                  "x-axis lower bound", "0");
+    opt.addOption("y", "yAxis", OptParser::OptType::value  , true,
+                  "y-axis name", "");
+    opt.addOption("", "spacing", OptParser::OptType::value  , true,
+                  "spacing between points", "1");
+    opt.addOption("s", "save", OptParser::OptType::value, true,
+                    "saves the source and .pdf", "");
+    opt.addOption("t", "title", OptParser::OptType::value  , true,
+                  "plot title", "");
+    opt.addOption("", "help"      , OptParser::OptType::trigger, true,
+                   "show this help message and exit");
+    parsed = opt.parse(argc, argv);
+    if (!parsed or opt.gotOption("help") or opt.getArgs().size() != 1)
+    {
+        cerr << "usage: " << argv[0] << " <.h5/manifest file> <options> " << endl;
+        cerr << endl << "Possible options:" << endl << opt << endl;
+        
+        return EXIT_FAILURE;
+    }
+    
+    plotFileName = opt.getArgs().front();
+    imag         = opt.gotOption("i");
+    xName        = opt.optionValue("x");
+    xLow         = opt.optionValue<double>("l");
+    yName        = opt.optionValue("y");
+    spacing      = opt.optionValue<double>("spacing");
+    save         = opt.optionValue("s");
+    title        = opt.optionValue("t");
+    outFileName  = opt.optionValue<string>("o");
+
+
+    if(plotFileName.find(".h5") == string::npos)
+    {
+        inFileName = readManifest(plotFileName);
+    }
+    
+    // load and plot file(s) /////////////////////////////////////////////////////////
+    DMatSample tmp;
+    Index      nt;
+    Plot       p;
+    DVec       tAxis;
+
+
+    if(inFileName.size() == 0)
+    {
+        tmp     = Io::load<DMatSample>(plotFileName);
+        nt      = tmp[central].rows();
+        if(imag)
+        {
+            tmp     = tmp.block(0, 1, nt, 1);
+        }
+        else
+        {
+            tmp     = tmp.block(0, 0, nt, 1);
+        }
+        xHigh= xLow+spacing*(nt-1);
+        tAxis.setLinSpaced(nt, xLow, xHigh);
+        p  << PlotData(tAxis, tmp); 
+    
+    }
+    else
+    {
+        tmp = Io::load<DMatSample>(inFileName[0]);
+        nt      = tmp[central].rows();
+        xHigh= xLow+spacing*(nt-1);
+        tAxis.setLinSpaced(nt, xLow, xHigh);
+
+        for(unsigned long i = 0; i < inFileName.size(); i++)
+        {
+            plotFileName = inFileName[i];
+            tmp = Io::load<DMatSample>(plotFileName);
+        if(imag)
+        {
+            tmp     = tmp.block(0, 1, nt, 1);
+        }
+        else
+        {
+            tmp     = tmp.block(0, 0, nt, 1);
+        }   
+            p  << PlotData(tAxis, tmp);
+            
+        }
+        
+    }
+
+    p << Label(xName, Axis::x);
+    p << Label(yName, Axis::y);
+    p << PlotRange(Axis::x, xLow, xHigh);
+    p << Caption(title);
+
+    if(save != "")
+    {
+        cout << "Saving plot and source code to " << save << endl;
+        p.save(save + "/" + title);
+    }
+    
+    else
+    {
+        cout << "Displaying plot..." << endl;
+        p.display();    
+    }
+    
+    
+
+    // output //////////////////////////////////////////////////////////////////
+    if (!outFileName.empty())
+    {
+        Io::save(tmp, outFileName);
+        cout << "File saved as: " << outFileName << endl;
+    }
+
+    return EXIT_SUCCESS;
+}