diff --git a/lib/Core/Utilities.hpp b/lib/Core/Utilities.hpp index 0ca5e4a..4c0ffd8 100644 --- a/lib/Core/Utilities.hpp +++ b/lib/Core/Utilities.hpp @@ -108,6 +108,23 @@ inline std::string strFrom(const T x) } // specialization for vectors +template<> +inline std::vector strTo>(const std::string &str) +{ + std::vector res; + std::vector vbuf; + double buf; + std::istringstream stream(str); + + while (!stream.eof()) + { + stream >> buf; + res.push_back(buf); + } + + return res; +} + template<> inline DVec strTo(const std::string &str) { diff --git a/lib/Physics/CorrelatorFitter.cpp b/lib/Physics/CorrelatorFitter.cpp index af7b4f3..e450ab2 100644 --- a/lib/Physics/CorrelatorFitter.cpp +++ b/lib/Physics/CorrelatorFitter.cpp @@ -253,16 +253,39 @@ DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts) } } -DMatSample CorrelatorUtils::fold(const DMatSample &c) +DMatSample CorrelatorUtils::fold(const DMatSample &c, const CorrelatorModels::ModelPar &par) { const Index nt = c[central].rows(); DMatSample buf = c; - - FOR_STAT_ARRAY(buf, s) + int sign; + bool fold = false; + + switch (par.type) { - for (Index t = 0; t < nt; ++t) + 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) { - buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt)); + for (Index t = 0; t < nt; ++t) + { + buf[s](t) = 0.5*(c[s](t) + sign*c[s]((nt - t) % nt)); + } } } diff --git a/lib/Physics/CorrelatorFitter.hpp b/lib/Physics/CorrelatorFitter.hpp index 694fd9b..816c0ca 100644 --- a/lib/Physics/CorrelatorFitter.hpp +++ b/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); + DMatSample fold(const DMatSample &c, const CorrelatorModels::ModelPar &par); DMatSample fourierTransform(const DMatSample &c, FFT &fft, const unsigned int dir = FFT::Forward); }; diff --git a/lib/Statistics/Histogram.cpp b/lib/Statistics/Histogram.cpp index 4f2cf06..832d4d7 100644 --- a/lib/Statistics/Histogram.cpp +++ b/lib/Statistics/Histogram.cpp @@ -146,6 +146,16 @@ 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.); diff --git a/lib/Statistics/Histogram.hpp b/lib/Statistics/Histogram.hpp index 43c24b1..8a9d126 100644 --- a/lib/Statistics/Histogram.hpp +++ b/lib/Statistics/Histogram.hpp @@ -52,6 +52,8 @@ public: const StatArray & getData(void) const; const StatArray & 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 diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 0676564..413e3b8 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -300,6 +300,67 @@ const XYStatData & XYSampleData::getData(void) } // fit ///////////////////////////////////////////////////////////////////////// + +void XYSampleData::fitSample(std::vector &minimizer, + const std::vector &v, + SampleFitResult &result, + DVec &init, + Index s) +{ + result.resize(nSample_); + result.chi2_.resize(nSample_); + result.model_.resize(v.size()); + FitResult sampleResult; + setDataToSample(s); + if (s == central) + { + sampleResult = data_.fit(minimizer, init, v); + init = sampleResult.segment(0, init.size()); + result.nPar_ = sampleResult.getNPar(); + result.nDof_ = sampleResult.nDof_; + result.parName_ = sampleResult.parName_; + result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); + } + else + { + sampleResult = data_.fit(*(minimizer.back()), init, v); + } + result[s] = sampleResult; + result.chi2_[s] = sampleResult.getChi2(); + for (unsigned int j = 0; j < v.size(); ++j) + { + result.model_[j].resize(nSample_); + result.model_[j][s] = sampleResult.getModel(j); + } + + +} + +SampleFitResult XYSampleData::fit(std::vector &minimizer, + const DVec &init, + const std::vector &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 &v, + Index s) +{ + vector mv{&minimizer}; + + return fit(mv, init, v, s); +} + SampleFitResult XYSampleData::fit(std::vector &minimizer, const DVec &init, const std::vector &v) @@ -307,43 +368,14 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, 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()); FOR_STAT_ARRAY(result, s) { - setDataToSample(s); - if (s == central) - { - sampleResult = data_.fit(minimizer, initCopy, v); - initCopy = sampleResult.segment(0, initCopy.size()); - if (verbCopy != Minimizer::Verbosity::Debug) - { - minimizer.back()->setVerbosity(Minimizer::Verbosity::Silent); - } - } - else - { - - sampleResult = data_.fit(*(minimizer.back()), initCopy, v); - } - result[s] = sampleResult; - result.chi2_[s] = sampleResult.getChi2(); - for (unsigned int j = 0; j < v.size(); ++j) - { - result.model_[j].resize(nSample_); - result.model_[j][s] = sampleResult.getModel(j); - } + 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; } diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index bb82748..963e4ee 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -103,9 +103,16 @@ public: // get internal XYStatData const XYStatData & getData(void); // fit - SampleFitResult fit(std::vector &minimizer, const DVec &init, + void fitSample(std::vector &minimizer, + const std::vector &v, + SampleFitResult &sampleResult, DVec &init, Index s); + SampleFitResult fit(std::vector &minimizer, const DVec &init, + const std::vector &v, Index s); + SampleFitResult fit(Minimizer &minimizer, const DVec &init, + const std::vector &v, Index s); + SampleFitResult fit(std::vector &minimizer, const DVec &init, const std::vector &v); - SampleFitResult fit(Minimizer &minimizer, const DVec &init, + SampleFitResult fit(Minimizer &minimizer, const DVec &init, const std::vector &v); template SampleFitResult fit(std::vector &minimizer, const DVec &init, diff --git a/physics/2pt-fit.cpp b/physics/2pt-fit.cpp index 413ddb2..6ba59be 100644 --- a/physics/2pt-fit.cpp +++ b/physics/2pt-fit.cpp @@ -114,6 +114,7 @@ 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 filter = {1., -2., 1.}; @@ -155,6 +156,11 @@ int main(int argc, char *argv[]) } } + if (fold) + { + corr = CorrelatorUtils::fold(corr,modelPar); + } + // fit ///////////////////////////////////////////////////////////////////// DVec init(nPar); NloptMinimizer globMin(NloptMinimizer::Algorithm::GN_CRS2_LM); diff --git a/utils/sample-fake.cpp b/utils/sample-fake.cpp index f34921b..5f1816c 100644 --- a/utils/sample-fake.cpp +++ b/utils/sample-fake.cpp @@ -18,30 +18,42 @@ */ #include +#include + using namespace std; using namespace Latan; int main(int argc, char *argv[]) { + OptParser opt; Index nSample; double val, err; string outFileName; - if (argc != 5) + opt.addOption("r", "seed" , OptParser::OptType::value, true, + "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 << " " << endl; + cerr << endl << "Possible options:" << endl << opt << endl; return EXIT_FAILURE; } + val = strTo(argv[1]); err = strTo(argv[2]); nSample = strTo(argv[3]); outFileName = argv[4]; random_device rd; - mt19937 gen(rd()); + SeedType seed = (opt.gotOption("r")) ? opt.optionValue("r") : rd(); + mt19937 gen(seed); normal_distribution<> dis(val, err); DSample res(nSample); @@ -59,4 +71,4 @@ int main(int argc, char *argv[]) Io::save(res, outFileName); return EXIT_SUCCESS; -} +} \ No newline at end of file diff --git a/utils/sample-read.cpp b/utils/sample-read.cpp index 52f72af..83e667e 100644 --- a/utils/sample-read.cpp +++ b/utils/sample-read.cpp @@ -38,9 +38,23 @@ int main(int argc, char *argv[]) { DMatSample s = Io::load(fileName); string name = Io::getFirstName(fileName); + Index nRows = s[central].rows(); + Index nCols = s[central].cols(); cout << scientific; - cout << "central value:\n" << s[central] << endl; - cout << "standard deviation:\n" << s.variance().cwiseSqrt() << endl; + cout << "central value +/- standard deviation\n" << endl; + cout << "Re:" << 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); @@ -51,8 +65,8 @@ int main(int argc, char *argv[]) DSample s = Io::load(fileName); string name = Io::getFirstName(fileName); cout << scientific; - cout << "central value:\n" << s[central] << endl; - cout << "standard deviation:\n" << sqrt(s.variance()) << endl; + cout << "central value +/- standard deviation\n" << endl; + cout << s[central] << " +/- " << sqrt(s.variance()) << endl; if (!copy.empty()) { Io::save(s, copy, File::Mode::write, name);