From 82d41d0b5d2f3e0f80496c24aa3de998ff6fb8e5 Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Thu, 24 Jan 2019 12:20:13 +0000 Subject: [PATCH] Now includes functioning const fit; tidying up now --- physics/2pt-fit.cpp | 140 +++++++++++++++++++++++++++++--------------- 1 file changed, 93 insertions(+), 47 deletions(-) diff --git a/physics/2pt-fit.cpp b/physics/2pt-fit.cpp index 8a6cd84..68b81ec 100644 --- a/physics/2pt-fit.cpp +++ b/physics/2pt-fit.cpp @@ -117,7 +117,7 @@ int main(int argc, char *argv[]) // make model ////////////////////////////////////////////////////////////// DoubleModel mod; - bool coshModel = false, linearModel = false; + bool coshModel = false, linearModel = false, constModel = false; if ((model == "exp") or (model == "exp1")) { @@ -185,11 +185,12 @@ int main(int argc, char *argv[]) } else if (model == "const") { - nPar = 1; - mod.setFunction([](const double *x, const double *p) + constModel = true; + nPar = 1; + mod.setFunction([](const double *x __dumb, const double *p) { return p[0]; - }, 0, nPar); + }, 1, nPar); } else { @@ -222,21 +223,35 @@ 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, "delta_m"); } + 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)); } + for (Index p = 2; p < nPar; p += 2) { init(p) = 2*init(p - 2); @@ -249,14 +264,24 @@ 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, -1.1*fabs(init(0))); + locMin.setLowLimit(p, 0.); + globMin.setHighLimit(p, 1.1*fabs(init(0))); + } else { globMin.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) + { + 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); @@ -286,69 +311,89 @@ 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; + // if (!constModel) + // { + // p << PlotRange(Axis::x, 0, nt - 1); + // if (!linearModel) + // { + // p << LogScale(Axis::y); + // } + // p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1); //<<<< problematic line for const fit + // p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1); + // p << Color("rgb 'red'") << PlotData(data.getData()); + // p.display(); + // } 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 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1); //<<<< problematic line for const fit + 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) + + // effective mass plot ////////////////////////////////////////////////////// + if (!constModel) { - FOR_STAT_ARRAY(effMass, s) + 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) { - 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(); } - 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; @@ -367,6 +412,7 @@ int main(int argc, char *argv[]) } } + // output ////////////////////////////////////////////////////////////////// if (!outFileName.empty()) {