portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2024-11-10 08:55:37 +00:00
Code Releases Activity

Last commit before adding const fit

Browse Source
This commit is contained in:
Andrew Zhen Ning Yong 2019-01-22 11:11:04 +00:00
parent 5dfe91ddb4
commit 1919497b49
1 changed files with 48 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

91
physics/2pt-fit.cpp
View File

@ -13,8 +13,6 @@ using namespace Latan;
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
cout << "Version edited by: Andrew Yong, 4/01/19\n" << endl;
// parse arguments ///////////////////////////////////////////////////////// // parse arguments /////////////////////////////////////////////////////////
OptParser opt; OptParser opt;
bool parsed, doPlot, doHeatmap, doCorr, fold; bool parsed, doPlot, doHeatmap, doCorr, fold;
@ -32,7 +30,7 @@ int main(int argc, char *argv[])
opt.addOption("s", "shift" , OptParser::OptType::value , true, opt.addOption("s", "shift" , OptParser::OptType::value , true,
"time variable shift", "0"); "time variable shift", "0");
opt.addOption("m", "model" , OptParser::OptType::value , true, opt.addOption("m", "model" , OptParser::OptType::value , true,
"fit model (exp|exp2|exp3|cosh|cosh2|cosh3|linear|<interpreter code>)", "cosh"); "fit model (exp|exp2|exp3|cosh|cosh2|cosh3|explin|const|<interpreter code>)", "cosh");
opt.addOption("" , "nPar" , OptParser::OptType::value , true, opt.addOption("" , "nPar" , OptParser::OptType::value , true,
"number of model parameters for custom models " "number of model parameters for custom models "
"(-1 if irrelevant)", "-1"); "(-1 if irrelevant)", "-1");
@ -119,9 +117,7 @@ int main(int argc, char *argv[])
// make model ////////////////////////////////////////////////////////////// // make model //////////////////////////////////////////////////////////////
DoubleModel mod; DoubleModel mod;
bool coshModel = false; bool coshModel = false, linearModel = false;
bool linearModel = false;
if ((model == "exp") or (model == "exp1")) if ((model == "exp") or (model == "exp1"))
{ {
@ -178,8 +174,7 @@ int main(int argc, char *argv[])
+ p[5]*(exp(-p[2]*x[0])+exp(-p[4]*(nt-x[0]))); + p[5]*(exp(-p[2]*x[0])+exp(-p[4]*(nt-x[0])));
}, 1, nPar); }, 1, nPar);
} }
else if (model == "explin")
else if (model == "linear")
{ {
linearModel = true; linearModel = true;
nPar = 2; nPar = 2;
@ -188,6 +183,14 @@ int main(int argc, char *argv[])
return p[1] - p[0]*x[0]; return p[1] - p[0]*x[0];
}, 1, nPar); }, 1, nPar);
} }
else if (model == "const")
{
nPar = 1;
mod.setFunction([](const double *x, const double *p)
{
return p[0];
}, 0, nPar);
}
else else
{ {
if (nPar > 0) if (nPar > 0)
@ -219,64 +222,41 @@ int main(int argc, char *argv[])
data.addXDim(nt, "t/a", true); data.addXDim(nt, "t/a", true);
data.addYDim("C(t)"); data.addYDim("C(t)");
data.setUnidimData(tvec, corr); data.setUnidimData(tvec, corr);
for (Index p = 0; p < nPar; p += 2) // naming parameters for (Index p = 0; p < nPar; p += 2)
{
if((model == "cosh") or (model =="cosh1") or (model == "cosh2") or (model == "cosh3"))
{ {
mod.parName().setName(p, "E_" + strFrom(p/2)); mod.parName().setName(p, "E_" + strFrom(p/2));
mod.parName().setName(p + 1, "Z_" + strFrom(p/2)); mod.parName().setName(p + 1, "Z_" + strFrom(p/2));
} }
if (linearModel)
else if(model == "linear")
{
mod.parName().setName(p, "dm");
mod.parName().setName(p + 1, "dA/A_0");
}
else // to edit when necessary
{
mod.parName().setName(p, "E_" + strFrom(p/2));
mod.parName().setName(p + 1, "Z_" + strFrom(p/2));
}
}
// initial values & limits//////////////////////////////////////////////////////////
if(model == "linear")
{ {
init(0) = data.y(nt/4, 0)[central] - data.y(nt/4 + 1, 0)[central]; init(0) = data.y(nt/4, 0)[central] - data.y(nt/4 + 1, 0)[central];
init(1) = data.y(nt/4,0)[central] + init(0)*nt/4; init(1) = data.y(nt/4, 0)[central] + nt/4*init(0);
cout << "init(0) = " << init(0) << "\tinit(1) = " << init(1) << endl;
double bound = 30.;
for (Index p = 0; p < nPar; p += 2) // setting appropriate limits for global min
{
globMin.setLowLimit(p, -bound*fabs(init(p)));
globMin.setHighLimit(p, bound*fabs(init(p)));
globMin.setLowLimit(p + 1, -bound*fabs(init(p + 1)));
globMin.setHighLimit(p + 1, bound*fabs(init(p + 1)));
} }
}
else else
{ {
init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]); 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)); init(1) = data.y(nt/4, 0)[central]/(exp(-init(0)*nt/4));
cout << "init(0) = " << init(0) << "\tinit(1) = " << init(1) << endl; }
for (Index p = 2; p < nPar; p += 2) for (Index p = 2; p < nPar; p += 2)
{ {
init(p) = 2*init(p - 2); init(p) = 2*init(p - 2);
init(p + 1) = init(p - 1)/2.; init(p + 1) = init(p - 1)/2.;
} }
for (Index p = 0; p < nPar; p += 2) for (Index p = 0; p < nPar; p += 2)
{ {
cout << "p: " << p << endl; if (linearModel)
globMin.setLowLimit(p, 0.); {
globMin.setHighLimit(p, 10.*init(p)); globMin.setLowLimit(p, -10.*fabs(init(p)));
globMin.setLowLimit(p + 1, -10.*init(p + 1)); globMin.setHighLimit(p, 10.*fabs(init(p)));
globMin.setHighLimit(p + 1, 10.*init(p + 1));
locMin.setLowLimit(p, 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)));
} }
globMin.setPrecision(0.001); globMin.setPrecision(0.001);
globMin.setMaxIteration(100000); globMin.setMaxIteration(100000);
@ -293,7 +273,6 @@ int main(int argc, char *argv[])
cout << "-- uncorrelated fit..." << endl; cout << "-- uncorrelated fit..." << endl;
} }
cout << "using model '" << model << "'" << endl; cout << "using model '" << model << "'" << endl;
cout << "svdTol: " << svdTol << endl;
data.setSvdTolerance(svdTol); data.setSvdTolerance(svdTol);
data.assumeYYCorrelated(false, 0, 0); data.assumeYYCorrelated(false, 0, 0);
fit = data.fit(unCorrMin, init, mod); fit = data.fit(unCorrMin, init, mod);
@ -317,9 +296,11 @@ int main(int argc, char *argv[])
Index maxT = (coshModel) ? (nt - 2) : (nt - 1); Index maxT = (coshModel) ? (nt - 2) : (nt - 1);
double e0, e0Err; double e0, e0Err;
p << Title("Correlated Fit");
p << PlotRange(Axis::x, 0, nt - 1); p << PlotRange(Axis::x, 0, nt - 1);
if (!linearModel)
{
p << LogScale(Axis::y); p << LogScale(Axis::y);
}
p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1); p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1);
p << Color("rgb 'blue'") << PlotFunction(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());
@ -340,6 +321,16 @@ int main(int argc, char *argv[])
} }
} }
} }
else if (linearModel)
{
FOR_STAT_ARRAY(effMass, s)
{
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_STAT_ARRAY(effMass, s)
@ -351,14 +342,12 @@ int main(int argc, char *argv[])
} }
} }
p.reset(); p.reset();
p << Title("Effective Mass");
p << PlotRange(Axis::x, 1, maxT); p << PlotRange(Axis::x, 1, maxT);
p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err); 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'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err);
p << Color("rgb 'blue'") << PlotHLine(e0); p << Color("rgb 'blue'") << PlotHLine(e0);
p << Color("rgb 'red'") << PlotData(effMassT, effMass); p << Color("rgb 'red'") << PlotData(effMassT, effMass);
p.display(); p.display();
p.save("test");
} }
if (doHeatmap) if (doHeatmap)
{ {