From b6c2efa666a4b678a40df4dde25512a8ca96282a Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Tue, 7 Jul 2020 10:07:05 +0100 Subject: [PATCH 01/26] 1) fold now requires second argument: modelPar to properly fold correlators with cosh/sinh form. 2) makeSinhModel has the appropriate alternate sign s.t the prefactor remains positive. 3) parameterGuess uses correct analytic form for sinh-type correlator's init(1), corresponding to the prefactor. --- lib/Physics/CorrelatorFitter.cpp | 57 +++++++++++++++++++++++++++----- lib/Physics/CorrelatorFitter.hpp | 2 +- 2 files changed, 50 insertions(+), 9 deletions(-) diff --git a/lib/Physics/CorrelatorFitter.cpp b/lib/Physics/CorrelatorFitter.cpp index 7784f5d..5061a02 100644 --- a/lib/Physics/CorrelatorFitter.cpp +++ b/lib/Physics/CorrelatorFitter.cpp @@ -84,7 +84,7 @@ DoubleModel CorrelatorModels::makeSinhModel(const Index nState, const Index nt) for (unsigned int i = 0; i < nState; ++i) { - res += p[2*i + 1]*(exp(-p[2*i]*x[0]) - exp(-p[2*i]*(nt - x[0]))); + res += p[2*i + 1]*(- exp(-p[2*i]*x[0]) + exp(-p[2*i]*(nt - x[0]))); } return res; @@ -200,8 +200,6 @@ DVec CorrelatorModels::parameterGuess(const DMatSample &corr, LATAN_ERROR(Argument, "correlator type is undefined"); break; case CorrelatorType::exp: - case CorrelatorType::cosh: - case CorrelatorType::sinh: init.resize(2*par.nState); init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1)); init(1) = corr[central](nt/4)/(exp(-init(0)*nt/4)); @@ -211,6 +209,26 @@ DVec CorrelatorModels::parameterGuess(const DMatSample &corr, init(p + 1) = init(p - 1)/2.; } break; + case CorrelatorType::cosh: + init.resize(2*par.nState); + init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1)); + init(1) = corr[central](nt/4)/(exp(-init(0)*nt/4) + exp(-init(0)*3*nt/4)); + for (Index p = 2; p < init.size(); p += 2) + { + init(p) = 2*init(p - 2); + init(p + 1) = init(p - 1)/2.; + } + break; + case CorrelatorType::sinh: + init.resize(2*par.nState); + init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1)); + init(1) = corr[central](nt/4)/( -exp(-init(0)*nt/4) +exp(-init(0)*3*nt/4)); + for (Index p = 2; p < init.size(); p += 2) + { + init(p) = 2*init(p - 2); + init(p + 1) = init(p - 1)/2.; + } + break; case CorrelatorType::linear: init.resize(2); init(0) = corr[central](nt/4) - corr[central](nt/4 + 1, 0); @@ -253,16 +271,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); }; From 3e3cdf2d69dfd8b04d77f7604cce7999b6b41f75 Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Tue, 7 Jul 2020 10:07:48 +0100 Subject: [PATCH 02/26] 2pt-fit.cpp now has correct number of arguments for fold function if --fold is called. --- physics/2pt-fit.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/physics/2pt-fit.cpp b/physics/2pt-fit.cpp index 2b19b4f..739f18f 100644 --- a/physics/2pt-fit.cpp +++ b/physics/2pt-fit.cpp @@ -110,10 +110,6 @@ int main(int argc, char *argv[]) nt = corr[central].rows(); corr = corr.block(0, 0, nt, 1); corr = CorrelatorUtils::shift(corr, shift); - if (fold) - { - corr = CorrelatorUtils::fold(corr); - } // make model ////////////////////////////////////////////////////////////// CorrelatorFitter fitter(corr); @@ -140,6 +136,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); From ddee922e72ed9375d37b3c39aa6d9de4a7f8feac Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Fri, 10 Jul 2020 08:33:28 +0100 Subject: [PATCH 03/26] 1) restored original definition of sinh in makeSinhModel() and its init in parameterGuess(). 2) fold definition includes sign factor due to cosh/sinh model. --- lib/Physics/CorrelatorFitter.cpp | 22 ++-------------------- 1 file changed, 2 insertions(+), 20 deletions(-) diff --git a/lib/Physics/CorrelatorFitter.cpp b/lib/Physics/CorrelatorFitter.cpp index 5061a02..02cf5eb 100644 --- a/lib/Physics/CorrelatorFitter.cpp +++ b/lib/Physics/CorrelatorFitter.cpp @@ -84,7 +84,7 @@ DoubleModel CorrelatorModels::makeSinhModel(const Index nState, const Index nt) for (unsigned int i = 0; i < nState; ++i) { - res += p[2*i + 1]*(- exp(-p[2*i]*x[0]) + exp(-p[2*i]*(nt - x[0]))); + res += p[2*i + 1]*(exp(-p[2*i]*x[0]) - exp(-p[2*i]*(nt - x[0]))); } return res; @@ -200,29 +200,11 @@ DVec CorrelatorModels::parameterGuess(const DMatSample &corr, LATAN_ERROR(Argument, "correlator type is undefined"); break; case CorrelatorType::exp: - init.resize(2*par.nState); - init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1)); - init(1) = corr[central](nt/4)/(exp(-init(0)*nt/4)); - for (Index p = 2; p < init.size(); p += 2) - { - init(p) = 2*init(p - 2); - init(p + 1) = init(p - 1)/2.; - } - break; case CorrelatorType::cosh: - init.resize(2*par.nState); - init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1)); - init(1) = corr[central](nt/4)/(exp(-init(0)*nt/4) + exp(-init(0)*3*nt/4)); - for (Index p = 2; p < init.size(); p += 2) - { - init(p) = 2*init(p - 2); - init(p + 1) = init(p - 1)/2.; - } - break; case CorrelatorType::sinh: init.resize(2*par.nState); init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1)); - init(1) = corr[central](nt/4)/( -exp(-init(0)*nt/4) +exp(-init(0)*3*nt/4)); + init(1) = corr[central](nt/4)/(exp(-init(0)*nt/4)); for (Index p = 2; p < init.size(); p += 2) { init(p) = 2*init(p - 2); From 1b12f2ca2d5bd9cb6f55a10fd41536840cf26e2f Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Fri, 13 Nov 2020 07:56:58 +0000 Subject: [PATCH 04/26] latan-sample-fake now takes option -r/--seed where seed can be specified if fake bootstrap samples with 100% correlation are needed. --- utils/sample-fake.cpp | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) 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 From 9e8d53463581fb917bc506d6de73d561fc01b59b Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Thu, 20 May 2021 11:04:21 +0100 Subject: [PATCH 05/26] New methods in XYSampleData class to skip fits that are non-converging/large chi2PerDof. This is tracked by boolean goodFit_. --- lib/Statistics/XYSampleData.cpp | 55 ++++++++++++++++++++++++--------- lib/Statistics/XYSampleData.hpp | 3 ++ 2 files changed, 43 insertions(+), 15 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 6e2d2d6..19342c9 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -234,6 +234,21 @@ DVec XYSampleData::getYError(const Index j) return data_.getYError(j); } +bool XYSampleData::checkFit() +{ + return goodFit_; +} + +void XYSampleData::checkChi2PerDof(double Chi2PerDof) +{ + if(Chi2PerDof >= 2 or Chi2PerDof < 0 or isnan(Chi2PerDof)) + { + goodFit_ = false; + cerr << "chi2PerDof = " << Chi2PerDof << ". Aborting fit now." << endl; + } + +} + // get total fit variance matrix and its pseudo-inverse //////////////////////// const DMat & XYSampleData::getFitVarMat(void) { @@ -292,24 +307,34 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.resize(nSample_); result.chi2_.resize(nSample_); result.model_.resize(v.size()); + double chi2PerDof; + goodFit_ = true; FOR_STAT_ARRAY(result, s) { - setDataToSample(s); - if (s == central) + if(goodFit_) { - sampleResult = data_.fit(minimizer, initCopy, v); - initCopy = sampleResult.segment(0, initCopy.size()); - } - 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); + setDataToSample(s); + if (s == central) + { + sampleResult = data_.fit(minimizer, initCopy, v); + initCopy = sampleResult.segment(0, initCopy.size()); + chi2PerDof = sampleResult.getChi2PerDof(); + checkChi2PerDof(chi2PerDof); + } + else + { + sampleResult = data_.fit(*(minimizer.back()), initCopy, v); + chi2PerDof = sampleResult.getChi2PerDof(); + checkChi2PerDof(chi2PerDof); + } + 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); + + } } } result.nPar_ = sampleResult.getNPar(); diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index 6cba855..a66c9dd 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -91,6 +91,8 @@ public: const DMat & getXYVar(const Index i, const Index j); DVec getXError(const Index i); DVec getYError(const Index j); + bool checkFit(); // check fit candidate based on chi2PerDof + void checkChi2PerDof(double Chi2PerDof); // get total fit variance matrix and its pseudo-inverse const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); @@ -133,6 +135,7 @@ private: Index nSample_, dataSample_{central}; bool initData_{true}, computeVarMat_{true}; bool initXMap_{true}; + bool goodFit_{true}; // used to break minimisation if central sample chi2PerDof is bad }; /****************************************************************************** From c48e2be20b34844d23a213c9284fc8802887efe6 Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Thu, 20 May 2021 11:49:03 +0100 Subject: [PATCH 06/26] Added feature to reject fit if 25% of sample fits are bad (with chi2perdof > 2 or is a NaN) --- lib/Statistics/XYSampleData.cpp | 16 ++++++++++++++-- lib/Statistics/XYSampleData.hpp | 1 + 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 19342c9..b336b62 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -244,9 +244,15 @@ void XYSampleData::checkChi2PerDof(double Chi2PerDof) if(Chi2PerDof >= 2 or Chi2PerDof < 0 or isnan(Chi2PerDof)) { goodFit_ = false; - cerr << "chi2PerDof = " << Chi2PerDof << ". Aborting fit now." << endl; } +} +void XYSampleData::checkChi2PerDof(double Chi2PerDof, unsigned int &counter) +{ + if(Chi2PerDof >= 2 or Chi2PerDof < 0 or isnan(Chi2PerDof)) + { + counter++; + } } // get total fit variance matrix and its pseudo-inverse //////////////////////// @@ -308,6 +314,7 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.chi2_.resize(nSample_); result.model_.resize(v.size()); double chi2PerDof; + unsigned int badSampleFits = 0; goodFit_ = true; FOR_STAT_ARRAY(result, s) { @@ -325,7 +332,7 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, { sampleResult = data_.fit(*(minimizer.back()), initCopy, v); chi2PerDof = sampleResult.getChi2PerDof(); - checkChi2PerDof(chi2PerDof); + checkChi2PerDof(chi2PerDof, badSampleFits); } result[s] = sampleResult; result.chi2_[s] = sampleResult.getChi2(); @@ -335,6 +342,11 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.model_[j][s] = sampleResult.getModel(j); } + if(badSampleFits > 0.25*nSample_) + { + goodFit_ = false; + cerr << "At least " << 0.25*nSample_ << " of the sample fits have bad chi2/dof. Aborting fit." << endl; + } } } result.nPar_ = sampleResult.getNPar(); diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index a66c9dd..c5c61b5 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -93,6 +93,7 @@ public: DVec getYError(const Index j); bool checkFit(); // check fit candidate based on chi2PerDof void checkChi2PerDof(double Chi2PerDof); + void checkChi2PerDof(double Chi2PerDof, unsigned int &counter); // get total fit variance matrix and its pseudo-inverse const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); From 375b8fd038261403de535de76c5e9820d1dcb06f Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Fri, 21 May 2021 15:29:24 +0100 Subject: [PATCH 07/26] Removed criteria of rejecting fits if sample fit chi2PerDof is bad. --- lib/Statistics/XYSampleData.cpp | 14 -------------- lib/Statistics/XYSampleData.hpp | 1 - 2 files changed, 15 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index b336b62..1d9b0cb 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -247,14 +247,6 @@ void XYSampleData::checkChi2PerDof(double Chi2PerDof) } } -void XYSampleData::checkChi2PerDof(double Chi2PerDof, unsigned int &counter) -{ - if(Chi2PerDof >= 2 or Chi2PerDof < 0 or isnan(Chi2PerDof)) - { - counter++; - } -} - // get total fit variance matrix and its pseudo-inverse //////////////////////// const DMat & XYSampleData::getFitVarMat(void) { @@ -314,7 +306,6 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.chi2_.resize(nSample_); result.model_.resize(v.size()); double chi2PerDof; - unsigned int badSampleFits = 0; goodFit_ = true; FOR_STAT_ARRAY(result, s) { @@ -342,11 +333,6 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.model_[j][s] = sampleResult.getModel(j); } - if(badSampleFits > 0.25*nSample_) - { - goodFit_ = false; - cerr << "At least " << 0.25*nSample_ << " of the sample fits have bad chi2/dof. Aborting fit." << endl; - } } } result.nPar_ = sampleResult.getNPar(); diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index c5c61b5..a66c9dd 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -93,7 +93,6 @@ public: DVec getYError(const Index j); bool checkFit(); // check fit candidate based on chi2PerDof void checkChi2PerDof(double Chi2PerDof); - void checkChi2PerDof(double Chi2PerDof, unsigned int &counter); // get total fit variance matrix and its pseudo-inverse const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); From af31d1564dd0ed5b56327cfaf48f16be9e73f0ca Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Fri, 21 May 2021 15:33:34 +0100 Subject: [PATCH 08/26] Previous commit failed to compile: needed to remove leftover lines. --- lib/Statistics/XYSampleData.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 1d9b0cb..ae5edba 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -323,7 +323,6 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, { sampleResult = data_.fit(*(minimizer.back()), initCopy, v); chi2PerDof = sampleResult.getChi2PerDof(); - checkChi2PerDof(chi2PerDof, badSampleFits); } result[s] = sampleResult; result.chi2_[s] = sampleResult.getChi2(); @@ -331,7 +330,6 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, { result.model_[j].resize(nSample_); result.model_[j][s] = sampleResult.getModel(j); - } } } From 80e3c27d8e3991379a75dedd250f26a8b851a9ab Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Wed, 25 Aug 2021 14:53:52 +0100 Subject: [PATCH 09/26] Added private variables to set bound on chi2perDof. This is used tas criteria for fit scans. TODO: modify to exception handling. --- lib/Statistics/XYSampleData.cpp | 8 +++++++- lib/Statistics/XYSampleData.hpp | 3 +++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index ae5edba..73e5e51 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -234,6 +234,12 @@ DVec XYSampleData::getYError(const Index j) return data_.getYError(j); } +void XYSampleData::setChi2PerDofBound(double upper_lim, double lower_lim) +{ + chi2PerDofu_ = upper_lim; + chi2PerDofl_ = lower_lim; +} + bool XYSampleData::checkFit() { return goodFit_; @@ -241,7 +247,7 @@ bool XYSampleData::checkFit() void XYSampleData::checkChi2PerDof(double Chi2PerDof) { - if(Chi2PerDof >= 2 or Chi2PerDof < 0 or isnan(Chi2PerDof)) + if(Chi2PerDof >= chi2PerDofu_ or Chi2PerDof < chi2PerDofl_ or isnan(Chi2PerDof)) { goodFit_ = false; } diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index a66c9dd..a18103a 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -91,6 +91,8 @@ public: const DMat & getXYVar(const Index i, const Index j); DVec getXError(const Index i); DVec getYError(const Index j); + // fit criteria + void setChi2PerDofBound(double upper_lim, double lower_lim); bool checkFit(); // check fit candidate based on chi2PerDof void checkChi2PerDof(double Chi2PerDof); // get total fit variance matrix and its pseudo-inverse @@ -136,6 +138,7 @@ private: bool initData_{true}, computeVarMat_{true}; bool initXMap_{true}; bool goodFit_{true}; // used to break minimisation if central sample chi2PerDof is bad + double chi2PerDofu_{1e10}, chi2PerDofl_{0}; }; /****************************************************************************** From f0739047c363f9980e8780efb14a11dbce577540 Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Wed, 25 Aug 2021 15:59:00 +0100 Subject: [PATCH 10/26] XYSampleData throws Runtime exception error if central fit's chi2PerDof exceeds user-set/default bounds for chi2PerDof. This can be used to skip fits with bad fitranges. --- lib/Statistics/XYSampleData.cpp | 19 ++++++------------- lib/Statistics/XYSampleData.hpp | 4 +--- 2 files changed, 7 insertions(+), 16 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 73e5e51..e455518 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -20,6 +20,7 @@ #include #include #include +// #include using namespace std; using namespace Latan; @@ -234,22 +235,18 @@ DVec XYSampleData::getYError(const Index j) return data_.getYError(j); } -void XYSampleData::setChi2PerDofBound(double upper_lim, double lower_lim) +void XYSampleData::setChi2PerDofBound(double uLim, double lLim) { - chi2PerDofu_ = upper_lim; - chi2PerDofl_ = lower_lim; -} - -bool XYSampleData::checkFit() -{ - return goodFit_; + chi2PerDofu_ = uLim; + chi2PerDofl_ = lLim; } void XYSampleData::checkChi2PerDof(double Chi2PerDof) { if(Chi2PerDof >= chi2PerDofu_ or Chi2PerDof < chi2PerDofl_ or isnan(Chi2PerDof)) { - goodFit_ = false; + string msg = "central chi2perDof = " + to_string(Chi2PerDof) + " not within user-set bound of " + to_string(chi2PerDofl_) + " < chi2PerDof < " + to_string(chi2PerDofu_) ; + LATAN_ERROR(Runtime, msg); } } @@ -312,11 +309,8 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.chi2_.resize(nSample_); result.model_.resize(v.size()); double chi2PerDof; - goodFit_ = true; FOR_STAT_ARRAY(result, s) { - if(goodFit_) - { setDataToSample(s); if (s == central) { @@ -337,7 +331,6 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.model_[j].resize(nSample_); result.model_[j][s] = sampleResult.getModel(j); } - } } result.nPar_ = sampleResult.getNPar(); result.nDof_ = sampleResult.nDof_; diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index a18103a..5109b6c 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -92,8 +92,7 @@ public: DVec getXError(const Index i); DVec getYError(const Index j); // fit criteria - void setChi2PerDofBound(double upper_lim, double lower_lim); - bool checkFit(); // check fit candidate based on chi2PerDof + void setChi2PerDofBound(double uLim, double lLim); void checkChi2PerDof(double Chi2PerDof); // get total fit variance matrix and its pseudo-inverse const DMat & getFitVarMat(void); @@ -137,7 +136,6 @@ private: Index nSample_, dataSample_{central}; bool initData_{true}, computeVarMat_{true}; bool initXMap_{true}; - bool goodFit_{true}; // used to break minimisation if central sample chi2PerDof is bad double chi2PerDofu_{1e10}, chi2PerDofl_{0}; }; From 9c5ade4989650a375ab539c84d6cbeac81ec8cca Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Wed, 22 Sep 2021 12:35:41 +0100 Subject: [PATCH 11/26] Changed fit criteria from chi2PerDof to pValue. --- lib/Statistics/XYSampleData.cpp | 13 ++++++------- lib/Statistics/XYSampleData.hpp | 6 +++--- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index e455518..bdb4f19 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -20,7 +20,6 @@ #include #include #include -// #include using namespace std; using namespace Latan; @@ -235,17 +234,17 @@ DVec XYSampleData::getYError(const Index j) return data_.getYError(j); } -void XYSampleData::setChi2PerDofBound(double uLim, double lLim) +void XYSampleData::setPValueBound(double uLim, double lLim) { - chi2PerDofu_ = uLim; - chi2PerDofl_ = lLim; + pValMax = uLim; + pValMin = lLim; } -void XYSampleData::checkChi2PerDof(double Chi2PerDof) +void XYSampleData::checkPValue(double pValue) { - if(Chi2PerDof >= chi2PerDofu_ or Chi2PerDof < chi2PerDofl_ or isnan(Chi2PerDof)) + if(pValue >= pValMax or pValue <= pValMin or isnan(pValue)) { - string msg = "central chi2perDof = " + to_string(Chi2PerDof) + " not within user-set bound of " + to_string(chi2PerDofl_) + " < chi2PerDof < " + to_string(chi2PerDofu_) ; + string msg = "central pValue = " + to_string(pValue) + " not within user-set bound of " + to_string(pValMin) + " < pValue < " + to_string(pValMax) ; LATAN_ERROR(Runtime, msg); } } diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index 5109b6c..f6bcb5f 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -92,8 +92,8 @@ public: DVec getXError(const Index i); DVec getYError(const Index j); // fit criteria - void setChi2PerDofBound(double uLim, double lLim); - void checkChi2PerDof(double Chi2PerDof); + void setPValueBound(double uLim, double lLim); + void checkPValue(double Chi2PerDof); // get total fit variance matrix and its pseudo-inverse const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); @@ -136,7 +136,7 @@ private: Index nSample_, dataSample_{central}; bool initData_{true}, computeVarMat_{true}; bool initXMap_{true}; - double chi2PerDofu_{1e10}, chi2PerDofl_{0}; + double pValMax{1}, pValMin{0}; }; /****************************************************************************** From 28aff209c4e399ce86e8a4802d828c5d520e86fb Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Wed, 22 Sep 2021 12:38:54 +0100 Subject: [PATCH 12/26] New function fitSample to make XYSampleData::fit more modular: Eg application: In fit scans, one may only want central fit for the uncorr fit to speed up process. --- lib/Statistics/XYSampleData.cpp | 69 +++++++++++++++++++++------------ lib/Statistics/XYSampleData.hpp | 6 ++- 2 files changed, 48 insertions(+), 27 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index bdb4f19..0bad9ce 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -294,9 +294,40 @@ const XYStatData & XYSampleData::getData(void) } // fit ///////////////////////////////////////////////////////////////////////// + +void XYSampleData::fitSample(std::vector &minimizer, + const std::vector &v, + SampleFitResult &result, + FitResult &sampleResult, + DVec &init, + Index s) +{ + double pValue = 0; + setDataToSample(s); + if (s == central) + { + sampleResult = data_.fit(minimizer, init, v); + init = sampleResult.segment(0, init.size()); + pValue = sampleResult.getPValue(); + checkPValue(pValue); + } + 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) + const std::vector &v, + bool centralSample) { computeVarMat(); @@ -307,29 +338,16 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.resize(nSample_); result.chi2_.resize(nSample_); result.model_.resize(v.size()); - double chi2PerDof; - FOR_STAT_ARRAY(result, s) + if(centralSample) { - setDataToSample(s); - if (s == central) - { - sampleResult = data_.fit(minimizer, initCopy, v); - initCopy = sampleResult.segment(0, initCopy.size()); - chi2PerDof = sampleResult.getChi2PerDof(); - checkChi2PerDof(chi2PerDof); - } - else - { - sampleResult = data_.fit(*(minimizer.back()), initCopy, v); - chi2PerDof = sampleResult.getChi2PerDof(); - } - 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, sampleResult, initCopy, central); + } + else + { + FOR_STAT_ARRAY(result, s) + { + fitSample(minimizer, v, result, sampleResult, initCopy, s); + } } result.nPar_ = sampleResult.getNPar(); result.nDof_ = sampleResult.nDof_; @@ -340,11 +358,12 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, - const std::vector &v) + const std::vector &v, + bool centralSample) { vector mv{&minimizer}; - return fit(mv, init, v); + return fit(mv, init, v, centralSample); } // residuals /////////////////////////////////////////////////////////////////// diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index f6bcb5f..ac3a464 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -102,10 +102,12 @@ public: // get internal XYStatData const XYStatData & getData(void); // fit + void fitSample(std::vector &minimizer, const std::vector &v, + SampleFitResult &result, FitResult &sampleResult, DVec &init, Index s); SampleFitResult fit(std::vector &minimizer, const DVec &init, - const std::vector &v); + const std::vector &v, bool centralSample = false); SampleFitResult fit(Minimizer &minimizer, const DVec &init, - const std::vector &v); + const std::vector &v, bool centralSample = false); template SampleFitResult fit(std::vector &minimizer, const DVec &init, const DoubleModel &model, const Ts... models); From 6eca1e6fc6d01cf5d8689788523ab8f190feceaa Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Thu, 7 Oct 2021 17:41:36 +0100 Subject: [PATCH 13/26] i) removed pvalue fit criteria ii) overloaded fit() to output FitResult object of a particular sample --- lib/Statistics/XYSampleData.cpp | 74 ++++++++++++++++----------------- lib/Statistics/XYSampleData.hpp | 20 +++++---- 2 files changed, 47 insertions(+), 47 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 0bad9ce..48f6d4d 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -234,21 +234,6 @@ DVec XYSampleData::getYError(const Index j) return data_.getYError(j); } -void XYSampleData::setPValueBound(double uLim, double lLim) -{ - pValMax = uLim; - pValMin = lLim; -} - -void XYSampleData::checkPValue(double pValue) -{ - if(pValue >= pValMax or pValue <= pValMin or isnan(pValue)) - { - string msg = "central pValue = " + to_string(pValue) + " not within user-set bound of " + to_string(pValMin) + " < pValue < " + to_string(pValMax) ; - LATAN_ERROR(Runtime, msg); - } -} - // get total fit variance matrix and its pseudo-inverse //////////////////////// const DMat & XYSampleData::getFitVarMat(void) { @@ -297,37 +282,51 @@ const XYStatData & XYSampleData::getData(void) void XYSampleData::fitSample(std::vector &minimizer, const std::vector &v, - SampleFitResult &result, FitResult &sampleResult, DVec &init, Index s) { - double pValue = 0; setDataToSample(s); if (s == central) { sampleResult = data_.fit(minimizer, init, v); init = sampleResult.segment(0, init.size()); - pValue = sampleResult.getPValue(); - checkPValue(pValue); } 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); - } + +} + +FitResult XYSampleData::fit(std::vector &minimizer, + const DVec &init, + const std::vector &v, + Index s) +{ + computeVarMat(); + + FitResult sampleResult; + DVec initCopy = init; + + fitSample(minimizer, v, sampleResult, initCopy, s); + + return sampleResult; +} + +FitResult 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, - bool centralSample) + const std::vector &v) { computeVarMat(); @@ -338,15 +337,15 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.resize(nSample_); result.chi2_.resize(nSample_); result.model_.resize(v.size()); - if(centralSample) + FOR_STAT_ARRAY(result, s) { - fitSample(minimizer, v, result, sampleResult, initCopy, central); - } - else - { - FOR_STAT_ARRAY(result, s) + fitSample(minimizer, v, sampleResult, initCopy, s); + result[s] = sampleResult; + result.chi2_[s] = sampleResult.getChi2(); + for (unsigned int j = 0; j < v.size(); ++j) { - fitSample(minimizer, v, result, sampleResult, initCopy, s); + result.model_[j].resize(nSample_); + result.model_[j][s] = sampleResult.getModel(j); } } result.nPar_ = sampleResult.getNPar(); @@ -358,12 +357,11 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init, - const std::vector &v, - bool centralSample) + const std::vector &v) { vector mv{&minimizer}; - return fit(mv, init, v, centralSample); + return fit(mv, init, v); } // residuals /////////////////////////////////////////////////////////////////// diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index ac3a464..8e9def6 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -91,9 +91,6 @@ public: const DMat & getXYVar(const Index i, const Index j); DVec getXError(const Index i); DVec getYError(const Index j); - // fit criteria - void setPValueBound(double uLim, double lLim); - void checkPValue(double Chi2PerDof); // get total fit variance matrix and its pseudo-inverse const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); @@ -102,12 +99,17 @@ public: // get internal XYStatData const XYStatData & getData(void); // fit - void fitSample(std::vector &minimizer, const std::vector &v, - SampleFitResult &result, FitResult &sampleResult, DVec &init, Index s); - SampleFitResult fit(std::vector &minimizer, const DVec &init, - const std::vector &v, bool centralSample = false); - SampleFitResult fit(Minimizer &minimizer, const DVec &init, - const std::vector &v, bool centralSample = false); + void fitSample(std::vector &minimizer, + const std::vector &v, + FitResult &sampleResult, DVec &init, Index s); + FitResult fit(std::vector &minimizer, const DVec &init, + const std::vector &v, Index s); + FitResult 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, + const std::vector &v); template SampleFitResult fit(std::vector &minimizer, const DVec &init, const DoubleModel &model, const Ts... models); From 15a5471bef0222b4c6179971df4a5b172f0f3742 Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Fri, 14 Jan 2022 16:03:52 +0000 Subject: [PATCH 14/26] strTo<> now converts string of numbers to vector --- lib/Core/Utilities.hpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/lib/Core/Utilities.hpp b/lib/Core/Utilities.hpp index 0ca5e4a..448fb46 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; + vbuf.push_back(buf); + } + + return res; +} + template<> inline DVec strTo(const std::string &str) { From a389e01aa0b333e84ccd031cc315a4320665976d Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Fri, 14 Jan 2022 16:58:24 +0000 Subject: [PATCH 15/26] corrected own bug in strTo --- lib/Core/Utilities.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/Core/Utilities.hpp b/lib/Core/Utilities.hpp index 448fb46..4c0ffd8 100644 --- a/lib/Core/Utilities.hpp +++ b/lib/Core/Utilities.hpp @@ -119,7 +119,7 @@ inline std::vector strTo>(const std::string &str) while (!stream.eof()) { stream >> buf; - vbuf.push_back(buf); + res.push_back(buf); } return res; From db99951dd448a6d1b8b9005ab1397169ac8c6537 Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Fri, 14 Jan 2022 16:58:44 +0000 Subject: [PATCH 16/26] removed redundant pValue private member --- lib/Statistics/XYSampleData.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index 8e9def6..59e9bb7 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -140,7 +140,6 @@ private: Index nSample_, dataSample_{central}; bool initData_{true}, computeVarMat_{true}; bool initXMap_{true}; - double pValMax{1}, pValMin{0}; }; /****************************************************************************** From 493d641e2fe46cccf1d65827c8ac8fe75b87b1fc Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Fri, 4 Mar 2022 10:36:30 +0000 Subject: [PATCH 17/26] Modified fitSample routines to use SampleFitResult objects. --- lib/Statistics/XYSampleData.cpp | 63 +++++++++++++++++---------------- lib/Statistics/XYSampleData.hpp | 6 ++-- 2 files changed, 35 insertions(+), 34 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index f014c94..413e3b8 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -303,42 +303,58 @@ const XYStatData & XYSampleData::getData(void) void XYSampleData::fitSample(std::vector &minimizer, const std::vector &v, - FitResult &sampleResult, + 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()); + 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); + } + } -FitResult XYSampleData::fit(std::vector &minimizer, - const DVec &init, - const std::vector &v, - Index s) +SampleFitResult XYSampleData::fit(std::vector &minimizer, + const DVec &init, + const std::vector &v, + Index s) { computeVarMat(); - FitResult sampleResult; - DVec initCopy = init; + SampleFitResult result; + DVec initCopy = init; - fitSample(minimizer, v, sampleResult, initCopy, s); + fitSample(minimizer, v, result, initCopy, s); - return sampleResult; + return result; } -FitResult XYSampleData::fit(Minimizer &minimizer, - const DVec &init, - const std::vector &v, - Index s) +SampleFitResult XYSampleData::fit(Minimizer &minimizer, + const DVec &init, + const std::vector &v, + Index s) { vector mv{&minimizer}; @@ -352,29 +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) { - fitSample(minimizer, v, sampleResult, initCopy, s); - 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 32e4cbe..963e4ee 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -105,10 +105,10 @@ public: // fit void fitSample(std::vector &minimizer, const std::vector &v, - FitResult &sampleResult, DVec &init, Index s); - FitResult fit(std::vector &minimizer, const DVec &init, + SampleFitResult &sampleResult, DVec &init, Index s); + SampleFitResult fit(std::vector &minimizer, const DVec &init, const std::vector &v, Index s); - FitResult fit(Minimizer &minimizer, const DVec &init, + 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); From c442a437e59a5ad0ddaceac892f13117a7d81155 Mon Sep 17 00:00:00 2001 From: AndrewYongZhenNing Date: Fri, 25 Mar 2022 09:38:32 +0000 Subject: [PATCH 18/26] Restructured sample-read to print stat err next to central value. --- utils/sample-read.cpp | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) 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); From e02a4bf30d0f854f94f1a5f5e23387b529012839 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:20:38 +0000 Subject: [PATCH 19/26] data plots compatible with multi-column arrays --- lib/Core/Plot.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/lib/Core/Plot.cpp b/lib/Core/Plot.cpp index 42848ed..3a236cd 100644 --- a/lib/Core/Plot.cpp +++ b/lib/Core/Plot.cpp @@ -122,10 +122,10 @@ PlotData::PlotData(const DMatSample &x, const DMatSample &y, const bool abs) DMat d(x[central].rows(), 4); string usingCmd, tmpFileName; - d.col(0) = x[central]; - d.col(2) = y[central]; - d.col(1) = x.variance().cwiseSqrt(); - d.col(3) = y.variance().cwiseSqrt(); + d.col(0) = x[central].col(0); + d.col(2) = y[central].col(0); + d.col(1) = x.variance().cwiseSqrt().col(0); + d.col(3) = y.variance().cwiseSqrt().col(0); tmpFileName = dumpToTmpFile(d); pushTmpFile(tmpFileName); if (!abs) @@ -149,8 +149,8 @@ PlotData::PlotData(const DVec &x, const DMatSample &y, const bool abs) string usingCmd, tmpFileName; d.col(0) = x; - d.col(1) = y[central]; - d.col(2) = y.variance().cwiseSqrt(); + d.col(1) = y[central].col(0); + d.col(2) = y.variance().cwiseSqrt().col(0); tmpFileName = dumpToTmpFile(d); pushTmpFile(tmpFileName); if (!abs) @@ -173,9 +173,9 @@ PlotData::PlotData(const DMatSample &x, const DVec &y, const bool abs) DMat d(x[central].rows(), 3), xerr, yerr; string usingCmd, tmpFileName; - d.col(0) = x[central]; + d.col(0) = x[central].col(0); d.col(2) = y; - d.col(1) = x.variance().cwiseSqrt(); + d.col(1) = x.variance().cwiseSqrt().col(0); tmpFileName = dumpToTmpFile(d); pushTmpFile(tmpFileName); if (!abs) From 4436c575e676d815a60b52fa0d87eb1c92f6c661 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:21:10 +0000 Subject: [PATCH 20/26] correlation matrix plot compute dynamic range --- utils/sample-plot-corr.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/utils/sample-plot-corr.cpp b/utils/sample-plot-corr.cpp index 02316c0..c900dc0 100644 --- a/utils/sample-plot-corr.cpp +++ b/utils/sample-plot-corr.cpp @@ -68,6 +68,7 @@ int main(int argc, char *argv[]) var = sample.varianceMatrix(); corr = sample.correlationMatrix(); + cout << "dynamic range " << Math::svdDynamicRangeDb(corr) << " dB" << endl; p << PlotCorrMatrix(corr); p.display(); if (!outVarName.empty()) From 51a46edb2750e0afd7bf899c535a2a8ae1f4edfe Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:21:26 +0000 Subject: [PATCH 21/26] tool to merge samples --- utils/Makefile.am | 5 +++ utils/sample-merge.cpp | 94 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 99 insertions(+) create mode 100644 utils/sample-merge.cpp diff --git a/utils/Makefile.am b/utils/Makefile.am index 660649d..70c1b87 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -12,6 +12,7 @@ bin_PROGRAMS = \ latan-sample-element \ latan-sample-fake \ latan-sample-ft \ + latan-sample-merge \ latan-sample-plot \ latan-sample-plot-corr\ latan-sample-read \ @@ -37,6 +38,10 @@ latan_sample_ft_SOURCES = sample-ft.cpp latan_sample_ft_CXXFLAGS = $(COM_CXXFLAGS) latan_sample_ft_LDFLAGS = -L../lib/.libs -lLatAnalyze +latan_sample_merge_SOURCES = sample-merge.cpp +latan_sample_merge_CXXFLAGS = $(COM_CXXFLAGS) +latan_sample_merge_LDFLAGS = -L../lib/.libs -lLatAnalyze + latan_sample_plot_corr_SOURCES = sample-plot-corr.cpp latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS) latan_sample_plot_corr_LDFLAGS = -L../lib/.libs -lLatAnalyze diff --git a/utils/sample-merge.cpp b/utils/sample-merge.cpp new file mode 100644 index 0000000..11c2e87 --- /dev/null +++ b/utils/sample-merge.cpp @@ -0,0 +1,94 @@ +/* + * sample-merge.cpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2020 Antonin Portelli + * + * LatAnalyze 3 is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze 3 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with LatAnalyze 3. If not, see . + */ + +#include +#include + +using namespace std; +using namespace Latan; + +int main(int argc, char *argv[]) +{ + // argument parsing //////////////////////////////////////////////////////// + OptParser opt; + bool parsed; + string outFileName = ""; + vector fileName; + unsigned int n = 0; + + opt.addOption("o", "output", OptParser::OptType::value , true, + "output file name (default: result not saved)"); + opt.addOption("" , "help" , OptParser::OptType::trigger, true, + "show this help message and exit"); + parsed = opt.parse(argc, argv); + if (!parsed or (opt.getArgs().size() < 1) or opt.gotOption("help")) + { + cerr << "usage: " << argv[0]; + cerr << " ... " << endl; + cerr << endl << "Possible options:" << endl << opt << endl; + + return EXIT_FAILURE; + } + n = opt.getArgs().size(); + outFileName = opt.optionValue("o"); + for (unsigned int i = 0; i < n; ++i) + { + fileName.push_back(opt.getArgs()[i]); + } + + // figure out dimensions /////////////////////////////////////////////////// + Index nCol, nSample, totRows; + DMatSample buf; + + buf = Io::load(fileName[0]); + nSample = buf.size(); + totRows = buf[central].rows(); + nCol = buf[central].cols(); + for (unsigned int i = 1; i < n; ++i) + { + buf = Io::load(fileName[i]); + if (buf.size() != nSample) + { + LATAN_ERROR(Size, "sample size mismatch"); + } + if (buf[central].cols() != nCol) + { + LATAN_ERROR(Size, "column size mismatch"); + } + totRows += buf[central].rows(); + } + + // merge along rows //////////////////////////////////////////////////////// + DMatSample out(nSample, totRows, nCol); + Index rowo = 0, r; + + for (unsigned int i = 0; i < n; ++i) + { + buf = Io::load(fileName[i]); + r = buf[central].rows(); + out.block(rowo, 0, r, nCol) = buf; + rowo += r; + } + if (!outFileName.empty()) + { + Io::save(out, outFileName); + } + + return EXIT_SUCCESS; +} \ No newline at end of file From b9f61d8c17523dd9e299c396194e410842db66b5 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 17 Feb 2022 19:24:34 +0000 Subject: [PATCH 22/26] first skeleton for DWT --- lib/Makefile.am | 4 + lib/Numerical/DWT.cpp | 57 ++++ lib/Numerical/DWT.hpp | 52 ++++ lib/Numerical/DWTFilters.cpp | 528 +++++++++++++++++++++++++++++++++++ lib/Numerical/DWTFilters.hpp | 53 ++++ 5 files changed, 694 insertions(+) create mode 100644 lib/Numerical/DWT.cpp create mode 100644 lib/Numerical/DWT.hpp create mode 100644 lib/Numerical/DWTFilters.cpp create mode 100644 lib/Numerical/DWTFilters.hpp diff --git a/lib/Makefile.am b/lib/Makefile.am index 9fb5e8c..8f33d71 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -48,6 +48,8 @@ libLatAnalyze_la_SOURCES = \ Io/XmlReader.cpp \ Io/Xml/tinyxml2.cpp \ Numerical/Derivative.cpp \ + Numerical/DWT.cpp \ + Numerical/DWTFilters.cpp \ Numerical/GslFFT.cpp \ Numerical/GslHybridRootFinder.cpp\ Numerical/GslMinimizer.cpp \ @@ -92,6 +94,8 @@ HPPFILES = \ Io/IoObject.hpp \ Io/XmlReader.hpp \ Numerical/Derivative.hpp \ + Numerical/DWT.hpp \ + Numerical/DWTFilters.hpp \ Numerical/FFT.hpp \ Numerical/GslFFT.hpp \ Numerical/GslHybridRootFinder.hpp\ diff --git a/lib/Numerical/DWT.cpp b/lib/Numerical/DWT.cpp new file mode 100644 index 0000000..b5c49a7 --- /dev/null +++ b/lib/Numerical/DWT.cpp @@ -0,0 +1,57 @@ +/* + * DWT.cpp, part of LatAnalyze + * + * Copyright (C) 2013 - 2020 Antonin Portelli + * + * LatAnalyze is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with LatAnalyze. If not, see . + */ + +#include +#include + +using namespace std; +using namespace Latan; + +/****************************************************************************** + * DWT implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +DWT::DWT(const DWTFilter &filter) +: filter_(filter) +{} + +// convolution primitive /////////////////////////////////////////////////////// +DVec DWT::filterConvolution(const DVec &data, const DWTFilter &filter, + const Index offset) +{ + DVec res(data.size()); + + return res; +} + +// DWT ///////////////////////////////////////////////////////////////////////// +std::vector +DWT::forward(const DVec &data, const unsigned int level) const +{ + std::vector dwt(level); + + return dwt; +} + +DVec DWT::backward(const std::vector& dwt) const +{ + DVec res; + + return res; +} diff --git a/lib/Numerical/DWT.hpp b/lib/Numerical/DWT.hpp new file mode 100644 index 0000000..f1ee646 --- /dev/null +++ b/lib/Numerical/DWT.hpp @@ -0,0 +1,52 @@ +/* + * DWT.hpp, part of LatAnalyze + * + * Copyright (C) 2013 - 2020 Antonin Portelli + * + * LatAnalyze is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with LatAnalyze. If not, see . + */ + +#ifndef Latan_DWT_hpp_ +#define Latan_DWT_hpp_ + +#include +#include + +BEGIN_LATAN_NAMESPACE + +/****************************************************************************** + * Discrete wavelet transform class * + ******************************************************************************/ +class DWT +{ +public: + typedef std::pair DWTLevel; +public: + // constructor + DWT(const DWTFilter &filter); + // destructor + virtual ~DWT(void) = default; + // convolution primitive + static DVec filterConvolution(const DVec &data, const DWTFilter &filter, + const Index offset); + // DWT + std::vector forward(const DVec &data, const unsigned int level) const; + DVec backward(const std::vector& dwt) const; +private: + DWTFilter filter_; +}; + +END_LATAN_NAMESPACE + +#endif // Latan_DWT_hpp_ diff --git a/lib/Numerical/DWTFilters.cpp b/lib/Numerical/DWTFilters.cpp new file mode 100644 index 0000000..be3accf --- /dev/null +++ b/lib/Numerical/DWTFilters.cpp @@ -0,0 +1,528 @@ +/* + * DWTFilters.cpp, part of LatAnalyze + * + * Copyright (C) 2013 - 2020 Antonin Portelli + * + * LatAnalyze is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with LatAnalyze. If not, see . + */ + +#include +#include + +// cf. http://wavelets.pybytes.com +// *here we implement the reverse filters more convenient for convolutions* + +using namespace std; +using namespace Latan; + +#define FILTDICT(x) {#x, &DWTFilters::x} + +std::map DWTFilters::fromName = { + FILTDICT(haar), + FILTDICT(db2), + FILTDICT(db3), + FILTDICT(db4), + FILTDICT(db5), + FILTDICT(db6), + FILTDICT(bior13), + FILTDICT(bior15), + FILTDICT(bior22), + FILTDICT(bior24), + FILTDICT(bior31), + FILTDICT(bior33), + FILTDICT(bior35) +}; + +DWTFilter DWTFilters::haar = { + // fwdL + {0.7071067811865476, + 0.7071067811865476}, + // fwdH + {0.7071067811865476, + -0.7071067811865476}, + // bwdL + {0.7071067811865476, + 0.7071067811865476}, + // bwdH + {-0.7071067811865476, + 0.7071067811865476} +}; + +DWTFilter DWTFilters::db2 = { + // fwdL + {0.48296291314469025, + 0.836516303737469, + 0.22414386804185735, + -0.12940952255092145}, + // fwdH + {-0.12940952255092145, + -0.22414386804185735, + 0.836516303737469, + -0.48296291314469025}, + // bwdL + {-0.12940952255092145, + 0.22414386804185735, + 0.836516303737469, + 0.48296291314469025}, + // bwdH + {-0.48296291314469025, + 0.836516303737469, + -0.22414386804185735, + -0.12940952255092145} +}; + +DWTFilter DWTFilters::db3 = { + // fwdL + {0.3326705529509569, + 0.8068915093133388, + 0.4598775021193313, + -0.13501102001039084, + -0.08544127388224149, + 0.035226291882100656}, + // fwdH + {0.035226291882100656, + 0.08544127388224149, + -0.13501102001039084, + -0.4598775021193313, + 0.8068915093133388, + -0.3326705529509569}, + // bwdL + {0.035226291882100656, + -0.08544127388224149, + -0.13501102001039084, + 0.4598775021193313, + 0.8068915093133388, + 0.3326705529509569}, + // bwdH + {-0.3326705529509569, + 0.8068915093133388, + -0.4598775021193313, + -0.13501102001039084, + 0.08544127388224149, + 0.035226291882100656} +}; + +DWTFilter DWTFilters::db4 = { + // fwdL + {0.23037781330885523, + 0.7148465705525415, + 0.6308807679295904, + -0.02798376941698385, + -0.18703481171888114, + 0.030841381835986965, + 0.032883011666982945, + -0.010597401784997278}, + // fwdH + {-0.010597401784997278, + -0.032883011666982945, + 0.030841381835986965, + 0.18703481171888114, + -0.02798376941698385, + -0.6308807679295904, + 0.7148465705525415, + -0.23037781330885523}, + // bwdL + {-0.010597401784997278, + 0.032883011666982945, + 0.030841381835986965, + -0.18703481171888114, + -0.02798376941698385, + 0.6308807679295904, + 0.7148465705525415, + 0.23037781330885523}, + // bwdH + {-0.23037781330885523, + 0.7148465705525415, + -0.6308807679295904, + -0.02798376941698385, + 0.18703481171888114, + 0.030841381835986965, + -0.032883011666982945, + -0.010597401784997278} +}; + +DWTFilter DWTFilters::db5 = { + // fwdL + {0.160102397974125, + 0.6038292697974729, + 0.7243085284385744, + 0.13842814590110342, + -0.24229488706619015, + -0.03224486958502952, + 0.07757149384006515, + -0.006241490213011705, + -0.012580751999015526, + 0.003335725285001549}, + // fwdH + {0.003335725285001549, + 0.012580751999015526, + -0.006241490213011705, + -0.07757149384006515, + -0.03224486958502952, + 0.24229488706619015, + 0.13842814590110342, + -0.7243085284385744, + 0.6038292697974729, + -0.160102397974125}, + // bwdL + {0.003335725285001549, + -0.012580751999015526, + -0.006241490213011705, + 0.07757149384006515, + -0.03224486958502952, + -0.24229488706619015, + 0.13842814590110342, + 0.7243085284385744, + 0.6038292697974729, + 0.160102397974125}, + // bwdH + {-0.160102397974125, + 0.6038292697974729, + -0.7243085284385744, + 0.13842814590110342, + 0.24229488706619015, + -0.03224486958502952, + -0.07757149384006515, + -0.006241490213011705, + 0.012580751999015526, + 0.003335725285001549} +}; + +DWTFilter DWTFilters::db6 = { + // fwdL + {0.11154074335008017, + 0.4946238903983854, + 0.7511339080215775, + 0.3152503517092432, + -0.22626469396516913, + -0.12976686756709563, + 0.09750160558707936, + 0.02752286553001629, + -0.031582039318031156, + 0.0005538422009938016, + 0.004777257511010651, + -0.00107730108499558}, + // fwdH + {-0.00107730108499558, + -0.004777257511010651, + 0.0005538422009938016, + 0.031582039318031156, + 0.02752286553001629, + -0.09750160558707936, + -0.12976686756709563, + 0.22626469396516913, + 0.3152503517092432, + -0.7511339080215775, + 0.4946238903983854, + -0.11154074335008017}, + // bwdL + {-0.00107730108499558, + 0.004777257511010651, + 0.0005538422009938016, + -0.031582039318031156, + 0.02752286553001629, + 0.09750160558707936, + -0.12976686756709563, + -0.22626469396516913, + 0.3152503517092432, + 0.7511339080215775, + 0.4946238903983854, + 0.11154074335008017}, + // bwdH + {-0.11154074335008017, + 0.4946238903983854, + -0.7511339080215775, + 0.3152503517092432, + 0.22626469396516913, + -0.12976686756709563, + -0.09750160558707936, + 0.02752286553001629, + 0.031582039318031156, + 0.0005538422009938016, + -0.004777257511010651, + -0.00107730108499558} +}; + +DWTFilter DWTFilters::bior13 = { + // fwdL + {-0.08838834764831845, + 0.08838834764831845, + 0.7071067811865476, + 0.7071067811865476, + 0.08838834764831845, + -0.08838834764831845}, + // fwdH + {0.0, + 0.0, + 0.7071067811865476, + -0.7071067811865476, + 0.0, + 0.0}, + // bwdL + {0.0, + 0.0, + 0.7071067811865476, + 0.7071067811865476, + 0.0, + 0.0}, + // bwdH + {0.08838834764831845, + 0.08838834764831845, + -0.7071067811865476, + 0.7071067811865476, + -0.08838834764831845, + -0.08838834764831845} +}; + +DWTFilter DWTFilters::bior15 = { + // fwdL + {0.01657281518405971, + -0.01657281518405971, + -0.12153397801643787, + 0.12153397801643787, + 0.7071067811865476, + 0.7071067811865476, + 0.12153397801643787, + -0.12153397801643787, + -0.01657281518405971, + 0.01657281518405971}, + // fwdH + {0.0, + 0.0, + 0.0, + 0.0, + 0.7071067811865476, + -0.7071067811865476, + 0.0, + 0.0, + 0.0, + 0.0}, + // bwdL + {0.0, + 0.0, + 0.0, + 0.0, + 0.7071067811865476, + 0.7071067811865476, + 0.0, + 0.0, + 0.0, + 0.0}, + // bwdH + {-0.01657281518405971, + -0.01657281518405971, + 0.12153397801643787, + 0.12153397801643787, + -0.7071067811865476, + 0.7071067811865476, + -0.12153397801643787, + -0.12153397801643787, + 0.01657281518405971, + 0.01657281518405971} +}; + +DWTFilter DWTFilters::bior22 = { + // fwdL + {-0.1767766952966369, + 0.3535533905932738, + 1.0606601717798214, + 0.3535533905932738, + -0.1767766952966369, + 0.0}, + // fwdH + {0.0, + 0.0, + 0.3535533905932738, + -0.7071067811865476, + 0.3535533905932738, + 0.0}, + // bwdL + {0.0, + 0.0, + 0.3535533905932738, + 0.7071067811865476, + 0.3535533905932738, + 0.0}, + // bwdH + {0.1767766952966369, + 0.3535533905932738, + -1.0606601717798214, + 0.3535533905932738, + 0.1767766952966369, + 0.0} +}; + +DWTFilter DWTFilters::bior24 = { + // fwdL + {0.03314563036811942, + -0.06629126073623884, + -0.1767766952966369, + 0.4198446513295126, + 0.9943689110435825, + 0.4198446513295126, + -0.1767766952966369, + -0.06629126073623884, + 0.03314563036811942, + 0.0}, + // fwdH + {0.0, + 0.0, + 0.0, + 0.0, + 0.3535533905932738, + -0.7071067811865476, + 0.3535533905932738, + 0.0, + 0.0, + 0.0}, + // bwdL + {0.0, + 0.0, + 0.0, + 0.0, + 0.3535533905932738, + 0.7071067811865476, + 0.3535533905932738, + 0.0, + 0.0, + 0.0}, + // bwdH + {-0.03314563036811942, + -0.06629126073623884, + 0.1767766952966369, + 0.4198446513295126, + -0.9943689110435825, + 0.4198446513295126, + 0.1767766952966369, + -0.06629126073623884, + -0.03314563036811942, + 0.0} +}; + +DWTFilter DWTFilters::bior31 = { + // fwdL + {-0.3535533905932738, + 1.0606601717798214, + 1.0606601717798214, + -0.3535533905932738}, + // fwdH + {0.1767766952966369, + -0.5303300858899107, + 0.5303300858899107, + -0.1767766952966369}, + // bwdL + {0.1767766952966369, + 0.5303300858899107, + 0.5303300858899107, + 0.1767766952966369}, + // bwdH + {0.3535533905932738, + 1.0606601717798214, + -1.0606601717798214, + -0.3535533905932738} +}; + +DWTFilter DWTFilters::bior33 = { + // fwdL + {0.06629126073623884, + -0.19887378220871652, + -0.15467960838455727, + 0.9943689110435825, + 0.9943689110435825, + -0.15467960838455727, + -0.19887378220871652, + 0.06629126073623884}, + // fwdH + {0.0, + 0.0, + 0.1767766952966369, + -0.5303300858899107, + 0.5303300858899107, + -0.1767766952966369, + 0.0, + 0.0}, + // bwdL + {0.0, + 0.0, + 0.1767766952966369, + 0.5303300858899107, + 0.5303300858899107, + 0.1767766952966369, + 0.0, + 0.0}, + // bwdH + {-0.06629126073623884, + -0.19887378220871652, + 0.15467960838455727, + 0.9943689110435825, + -0.9943689110435825, + -0.15467960838455727, + 0.19887378220871652, + 0.06629126073623884} +}; + +DWTFilter DWTFilters::bior35 = { + // fwdL + {-0.013810679320049757, + 0.04143203796014927, + 0.052480581416189075, + -0.26792717880896527, + -0.07181553246425874, + 0.966747552403483, + 0.966747552403483, + -0.07181553246425874, + -0.26792717880896527, + 0.052480581416189075, + 0.04143203796014927, + -0.013810679320049757}, + // fwdH + {0.0, + 0.0, + 0.0, + 0.0, + 0.1767766952966369, + -0.5303300858899107, + 0.5303300858899107, + -0.1767766952966369, + 0.0, + 0.0, + 0.0, + 0.0}, + // bwdL + {0.0, + 0.0, + 0.0, + 0.0, + 0.1767766952966369, + 0.5303300858899107, + 0.5303300858899107, + 0.1767766952966369, + 0.0, + 0.0, + 0.0, + 0.0}, + // bwdH + {0.013810679320049757, + 0.04143203796014927, + -0.052480581416189075, + -0.26792717880896527, + 0.07181553246425874, + 0.966747552403483, + -0.966747552403483, + -0.07181553246425874, + 0.26792717880896527, + 0.052480581416189075, + -0.04143203796014927, + -0.013810679320049757} +}; diff --git a/lib/Numerical/DWTFilters.hpp b/lib/Numerical/DWTFilters.hpp new file mode 100644 index 0000000..8ceab9a --- /dev/null +++ b/lib/Numerical/DWTFilters.hpp @@ -0,0 +1,53 @@ +/* + * DWTFilters.hpp, part of LatAnalyze + * + * Copyright (C) 2013 - 2020 Antonin Portelli + * + * LatAnalyze is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with LatAnalyze. If not, see . + */ + +#ifndef Latan_DWTFilters_hpp_ +#define Latan_DWTFilters_hpp_ + +#include + +BEGIN_LATAN_NAMESPACE + +struct DWTFilter +{ + const std::vector fwdL, fwdH, bwdL, bwdH; +}; + +namespace DWTFilters +{ + extern DWTFilter haar; + extern DWTFilter db2; + extern DWTFilter db3; + extern DWTFilter db4; + extern DWTFilter db5; + extern DWTFilter db6; + extern DWTFilter bior13; + extern DWTFilter bior15; + extern DWTFilter bior22; + extern DWTFilter bior24; + extern DWTFilter bior31; + extern DWTFilter bior33; + extern DWTFilter bior35; + + extern std::map fromName; +} + +END_LATAN_NAMESPACE + +#endif // Latan_DWTFilters_hpp_ From 500210a2eb5249562373f56164b76e245cac70a7 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 18 Feb 2022 14:06:52 +0000 Subject: [PATCH 23/26] DWT working and tested --- examples/Makefile.am | 5 +++ examples/exDWT.cpp | 28 +++++++++++++ lib/Numerical/DWT.cpp | 92 ++++++++++++++++++++++++++++++++++++++++--- lib/Numerical/DWT.hpp | 7 +++- 4 files changed, 124 insertions(+), 8 deletions(-) create mode 100644 examples/exDWT.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index 123bafb..bf23b0d 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -9,6 +9,7 @@ endif noinst_PROGRAMS = \ exCompiledDoubleFunction\ exDerivative \ + exDWT \ exFit \ exFitSample \ exIntegrator \ @@ -30,6 +31,10 @@ exDerivative_SOURCES = exDerivative.cpp exDerivative_CXXFLAGS = $(COM_CXXFLAGS) exDerivative_LDFLAGS = -L../lib/.libs -lLatAnalyze +exDWT_SOURCES = exDWT.cpp +exDWT_CXXFLAGS = $(COM_CXXFLAGS) +exDWT_LDFLAGS = -L../lib/.libs -lLatAnalyze + exFit_SOURCES = exFit.cpp exFit_CXXFLAGS = $(COM_CXXFLAGS) exFit_LDFLAGS = -L../lib/.libs -lLatAnalyze diff --git a/examples/exDWT.cpp b/examples/exDWT.cpp new file mode 100644 index 0000000..082d8b5 --- /dev/null +++ b/examples/exDWT.cpp @@ -0,0 +1,28 @@ +#include + +using namespace std; +using namespace Latan; + +int main(void) +{ + DVec data, dataRec; + vector dataDWT; + DWT dwt(DWTFilters::db3); + + cout << "-- random data" << endl; + data.setRandom(16); + cout << data.transpose() << endl; + cout << "-- compute Daubechies 3 DWT" << endl; + dataDWT = dwt.forward(data, 4); + for (unsigned int l = 0; l < dataDWT.size(); ++l) + { + cout << "* level " << l << endl; + cout << "L= " << dataDWT[l].first.transpose() << endl; + cout << "H= " << dataDWT[l].second.transpose() << endl; + } + cout << "-- check inverse DWT" << endl; + dataRec = dwt.backward(dataDWT); + cout << "rel diff = " << 2.*(data - dataRec).norm()/(data + dataRec).norm() << endl; + + return EXIT_SUCCESS; +} diff --git a/lib/Numerical/DWT.cpp b/lib/Numerical/DWT.cpp index b5c49a7..bc0e03a 100644 --- a/lib/Numerical/DWT.cpp +++ b/lib/Numerical/DWT.cpp @@ -32,26 +32,106 @@ DWT::DWT(const DWTFilter &filter) {} // convolution primitive /////////////////////////////////////////////////////// -DVec DWT::filterConvolution(const DVec &data, const DWTFilter &filter, - const Index offset) +void DWT::filterConvolution(DVec &out, const DVec &data, + const std::vector &filter, const Index offset) { - DVec res(data.size()); + Index n = data.size(), nf = n*filter.size(); - return res; + out.resize(n); + out.fill(0.); + for (unsigned int i = 0; i < filter.size(); ++i) + { + FOR_VEC(out, j) + { + out(j) += filter[i]*data((j + i + nf - offset) % n); + } + } +} + +// downsampling/upsampling primitives ////////////////////////////////////////// +void DWT::downsample(DVec &out, const DVec &in) +{ + if (out.size() < in.size()/2) + { + LATAN_ERROR(Size, "output vector smaller than half the input vector size"); + } + for (Index i = 0; i < in.size(); i += 2) + { + out(i/2) = in(i); + } +} + +void DWT::upsample(DVec &out, const DVec &in) +{ + if (out.size() < 2*in.size()) + { + LATAN_ERROR(Size, "output vector smaller than twice the input vector size"); + } + out.segment(0, 2*in.size()).fill(0.); + for (Index i = 0; i < in.size(); i ++) + { + out(2*i) = in(i); + } } // DWT ///////////////////////////////////////////////////////////////////////// std::vector DWT::forward(const DVec &data, const unsigned int level) const { - std::vector dwt(level); + std::vector dwt(level); + DVec *finePt = const_cast(&data); + DVec tmp; + Index n = data.size(), o = filter_.fwdL.size()/2, minSize; + + minSize = 1; + for (unsigned int l = 0; l < level; ++l) minSize *= 2; + if (n < minSize) + { + LATAN_ERROR(Size, "data vector too small for a " + strFrom(level) + + "-level DWT (data size is " + strFrom(n) + ")"); + } + for (unsigned int l = 0; l < level; ++l) + { + n /= 2; + dwt[l].first.resize(n); + dwt[l].second.resize(n); + filterConvolution(tmp, *finePt, filter_.fwdL, o); + downsample(dwt[l].first, tmp); + filterConvolution(tmp, *finePt, filter_.fwdH, o); + downsample(dwt[l].second, tmp); + finePt = &dwt[l].first; + } return dwt; } DVec DWT::backward(const std::vector& dwt) const { - DVec res; + unsigned int level = dwt.size(); + Index n = dwt.back().second.size(), o = filter_.bwdL.size()/2 - 1; + DVec res, tmp, conv; + + res = dwt.back().first; + for (int l = level - 2; l >= 0; --l) + { + n *= 2; + if (dwt[l].second.size() != n) + { + LATAN_ERROR(Size, "DWT result size mismatch"); + } + } + n = dwt.back().second.size(); + for (int l = level - 1; l >= 0; --l) + { + n *= 2; + tmp.resize(n); + upsample(tmp, res); + filterConvolution(conv, tmp, filter_.bwdL, o); + res = conv; + upsample(tmp, dwt[l].second); + filterConvolution(conv, tmp, filter_.bwdH, o); + res += conv; + } return res; } diff --git a/lib/Numerical/DWT.hpp b/lib/Numerical/DWT.hpp index f1ee646..5745245 100644 --- a/lib/Numerical/DWT.hpp +++ b/lib/Numerical/DWT.hpp @@ -38,8 +38,11 @@ public: // destructor virtual ~DWT(void) = default; // convolution primitive - static DVec filterConvolution(const DVec &data, const DWTFilter &filter, - const Index offset); + static void filterConvolution(DVec &out, const DVec &data, + const std::vector &filter, const Index offset); + // downsampling/upsampling primitives + static void downsample(DVec &out, const DVec &in); + static void upsample(DVec &out, const DVec &in); // DWT std::vector forward(const DVec &data, const unsigned int level) const; DVec backward(const std::vector& dwt) const; From 8b259879fff7ff4a363ecbf8d6b492548a88172c Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 10 Mar 2022 08:16:40 +0000 Subject: [PATCH 24/26] utility to compute sample DWT --- utils/Makefile.am | 5 ++ utils/sample-dwt.cpp | 167 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 172 insertions(+) create mode 100644 utils/sample-dwt.cpp diff --git a/utils/Makefile.am b/utils/Makefile.am index 70c1b87..de68500 100644 --- a/utils/Makefile.am +++ b/utils/Makefile.am @@ -9,6 +9,7 @@ endif bin_PROGRAMS = \ latan-plot \ latan-sample-combine \ + latan-sample-dwt \ latan-sample-element \ latan-sample-fake \ latan-sample-ft \ @@ -26,6 +27,10 @@ latan_sample_combine_SOURCES = sample-combine.cpp latan_sample_combine_CXXFLAGS = $(COM_CXXFLAGS) latan_sample_combine_LDFLAGS = -L../lib/.libs -lLatAnalyze +latan_sample_dwt_SOURCES = sample-dwt.cpp +latan_sample_dwt_CXXFLAGS = $(COM_CXXFLAGS) +latan_sample_dwt_LDFLAGS = -L../lib/.libs -lLatAnalyze + latan_sample_element_SOURCES = sample-element.cpp latan_sample_element_CXXFLAGS = $(COM_CXXFLAGS) latan_sample_element_LDFLAGS = -L../lib/.libs -lLatAnalyze diff --git a/utils/sample-dwt.cpp b/utils/sample-dwt.cpp new file mode 100644 index 0000000..958eea6 --- /dev/null +++ b/utils/sample-dwt.cpp @@ -0,0 +1,167 @@ +/* + * sample-dwt.cpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs + * + * LatAnalyze 3 is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * LatAnalyze 3 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with LatAnalyze 3. If not, see . + */ + +#include +#include +#include +#include + +using namespace std; +using namespace Latan; + +int main(int argc, char *argv[]) +{ + // argument parsing //////////////////////////////////////////////////////// + OptParser opt; + bool parsed, doPlot, ss, saveAll; + unsigned int level; + string inFilename, outFilename, filterName; + + opt.addOption("l", "level", OptParser::OptType::value, true, + "number of levels", "1"); + opt.addOption("f", "filter", OptParser::OptType::value, true, + "filter name", "haar"); + opt.addOption("s", "ss", OptParser::OptType::trigger, true, + "super-sampling (inverse DWT on data)"); + opt.addOption("a", "all", OptParser::OptType::trigger, true, + "save all-levels"); + opt.addOption("o", "output", OptParser::OptType::value, true, + "output file name, or directory name if --all is used (default: result not saved)", ""); + opt.addOption("p", "plot", OptParser::OptType::trigger, true, + "show plot"); + opt.addOption("" , "help" , OptParser::OptType::trigger, true, + "show this help message and exit"); + parsed = opt.parse(argc, argv); + if (!parsed or (opt.getArgs().size() != 1) or opt.gotOption("help")) + { + cerr << "usage: " << argv[0]; + cerr << " " << endl; + cerr << endl << "Possible options:" << endl << opt << endl; + + return EXIT_FAILURE; + } + inFilename = opt.getArgs()[0]; + level = opt.optionValue("l"); + filterName = opt.optionValue("f"); + ss = opt.gotOption("s"); + saveAll = opt.gotOption("a"); + outFilename = opt.optionValue("o"); + doPlot = opt.gotOption("p"); + + // DWT ///////////////////////////////////////////////////////////////////// + DMatSample in = Io::load(inFilename), res; + Index nSample = in.size(), n = in[central].rows(); + vector out(ss ? 1 : level, DMatSample(nSample)), + outh(ss ? 0 : level, DMatSample(nSample)); + DWT dwt(*DWTFilters::fromName.at(filterName)); + vector dataDWT(level); + + if (!ss) + { + cout << "-- compute discrete wavelet transform" << endl; + cout << "filter '" << filterName << "' / " << level << " level(s)" << endl; + FOR_STAT_ARRAY(in, s) + { + dataDWT = dwt.forward(in[s].col(0), level); + for (unsigned int l = 0; l < level; ++l) + { + out[l][s] = dataDWT[l].first; + outh[l][s] = dataDWT[l].second; + } + } + } + else + { + Index ssn = n; + + cout << "-- compute inverse discrete wavelet transform" << endl; + cout << "filter '" << filterName << "' / " << level << " level(s)" << endl; + for (int l = level - 1; l >= 0; --l) + { + dataDWT[l].first.resize(ssn); + dataDWT[l].second.resize(ssn); + dataDWT[l].first.fill(0.); + dataDWT[l].second.fill(0.); + ssn *= 2; + } + FOR_STAT_ARRAY(in, s) + { + dataDWT.back().first = in[s].col(0); + out[0][s] = dwt.backward(dataDWT); + } + } + if (!outFilename.empty()) + { + if (!ss and saveAll) + { + Latan::mkdir(outFilename); + for (unsigned int l = 0; l < level; ++l) + { + Io::save(out[l], outFilename + "/L" + strFrom(l) + ".h5"); + Io::save(outh[l], outFilename + "/H" + strFrom(l) + ".h5"); + } + } + else + { + Io::save(out.back(), outFilename); + } + } + + // plot //////////////////////////////////////////////////////////////////// + if (doPlot) + { + Plot p; + DVec x; + + x.setLinSpaced(n, 0., n - 1.); + p << PlotRange(Axis::x, 0., n); + p << Title("original") << PlotData(x, in); + if (!ss) + { + Index ln = n, step = 1; + + for (unsigned int l = 0; l < level; ++l) + { + ln /= 2; + step *= 2; + x.setLinSpaced(ln, 0., n - step); + p << Title("level " + strFrom(l + 1) + " L") << PlotData(x, out[l]); + p << Title("level " + strFrom(l + 1) + " H") << PlotData(x, outh[l]); + } + p.display(); + } + else + { + double step = 1.; + DVec err; + + for (unsigned int l = 0; l < level; ++l) + { + step /= 2.; + } + x.setLinSpaced(out[0][central].size(), 0., n - step); + err = out[0].variance().cwiseSqrt(); + p << Title("supersampled") << Color("3") << PlotPredBand(x, out[0][central], err); + p << Color("3") << PlotLine(x, out[0][central]); + p.display(); + } + } + + return EXIT_SUCCESS; +} From db0855963242ac7ee927f433c93d98214bddcc57 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 Nov 2021 16:20:56 +0000 Subject: [PATCH 25/26] Update Readme.md --- Readme.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Readme.md b/Readme.md index 9137da0..c7079c3 100644 --- a/Readme.md +++ b/Readme.md @@ -1,6 +1,6 @@ # LatAnalyze -[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![DOI](https://zenodo.org/badge/10201777.svg)](https://zenodo.org/badge/latestdoi/10201777) +[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![DOI](https://zenodo.org/badge/10201777.svg)](https://zenodo.org/badge/latestdoi/10201777) [![Build Ubuntu](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-ubuntu.yml) [![Build macOS](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml/badge.svg)](https://github.com/aportelli/LatAnalyze/actions/workflows/build-macos.yml) ## Description LatAnalyze is a C++11 library for statistical data analysis based on bootstrap From be72d313642a1c6bdff8c486695ffdba635c3aee Mon Sep 17 00:00:00 2001 From: Andrew Zhen Ning Yong Date: Thu, 21 Apr 2022 07:33:25 +0100 Subject: [PATCH 26/26] Added histogram feature to return xMax and xMin. --- lib/Statistics/Histogram.cpp | 10 ++++++++++ lib/Statistics/Histogram.hpp | 2 ++ 2 files changed, 12 insertions(+) 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