diff --git a/physics/2pt-fit.cpp b/physics/2pt-fit.cpp index 977d3f5..df67ca3 100644 --- a/physics/2pt-fit.cpp +++ b/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|)", "cosh"); + "fit model (exp|exp2|exp3|sinh|cosh|cosh2|cosh3|explin|const|)", "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("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(p + 1, "Z_" + strFrom(p/2)); + mod.parName().setName(0, "const"); } + 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))); - globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1))); + if(!constModel) + { + 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); - effMassT.setLinSpaced(maxT, 1, maxT); - fitErr = fit.variance().cwiseSqrt(); - e0 = fit[central](0); - e0Err = fitErr(0); - if (coshModel) + if(savePlot != "") { - 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)) - /(2.*corr[s](t))); + for (Index t = 1; t < nt - 1; ++t) + { + effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1)) + /(2.*corr[s](t))); + } } } - } - else if (linearModel) - { - FOR_STAT_ARRAY(effMass, s) + else if (linearModel) { - 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 - { - FOR_STAT_ARRAY(effMass, s) + else { - 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()) { diff --git a/utils/Makefile.am b/utils/Makefile.am index a5747cd..660649d 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -7,6 +7,7 @@ endif endif bin_PROGRAMS = \ + latan-plot \ latan-sample-combine \ latan-sample-element \ latan-sample-fake \ @@ -14,7 +15,11 @@ bin_PROGRAMS = \ latan-sample-plot \ latan-sample-plot-corr\ latan-sample-read \ - latan-resample + 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) diff --git a/utils/plot.cpp b/utils/plot.cpp new file mode 100644 index 0000000..0238741 --- /dev/null +++ b/utils/plot.cpp @@ -0,0 +1,142 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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 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> " << 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("l"); + yName = opt.optionValue("y"); + spacing = opt.optionValue("spacing"); + save = opt.optionValue("s"); + title = opt.optionValue("t"); + outFileName = opt.optionValue("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(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(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(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; +}