diff --git a/.github/workflows/build-macos.yml b/.github/workflows/build-macos.yml new file mode 100644 index 0000000..587b668 --- /dev/null +++ b/.github/workflows/build-macos.yml @@ -0,0 +1,26 @@ +name: Build macOS + +on: [push] + +jobs: + build: + + runs-on: macos-11 + + steps: + - name: Checkout + uses: actions/checkout@v2 + - name: Install basic dependencies + run: brew install automake autoconf libtool bison flex + - name: Build dependencies + shell: bash + run: | + export PATH=/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:${PATH} + cd ci-scripts + ./install-deps.sh prefix 6 + - name: Build LatAnalyze + shell: bash + run: | + export PATH=/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:${PATH} + cd ci-scripts + ./install-latan.sh prefix 6 diff --git a/.github/workflows/build-ubuntu.yml b/.github/workflows/build-ubuntu.yml new file mode 100644 index 0000000..3442ee5 --- /dev/null +++ b/.github/workflows/build-ubuntu.yml @@ -0,0 +1,26 @@ +name: Build Ubuntu + +on: [push] + +jobs: + build: + + runs-on: ubuntu-20.04 + + steps: + - name: Checkout + uses: actions/checkout@v2 + - name: Install basic dependencies + run: | + sudo bash -c "$(wget -O - https://apt.llvm.org/llvm.sh)" + sudo apt install cmake bison flex + - name: Build dependencies + shell: bash + run: | + cd ci-scripts + CC=clang CXX=clang++ ./install-deps.sh prefix 6 + - name: Build LatAnalyze + shell: bash + run: | + cd ci-scripts + CC=clang CXX=clang++ ./install-latan.sh prefix 6 diff --git a/.gitignore b/.gitignore index 6ed0930..7ef4ca1 100644 --- a/.gitignore +++ b/.gitignore @@ -16,6 +16,7 @@ autom4te.cache/* *.in~ config.h* configure +configure~ .buildutils/* aclocal.m4 diff --git a/Readme.md b/Readme.md index 703747c..9137da0 100644 --- a/Readme.md +++ b/Readme.md @@ -1,21 +1,6 @@ # LatAnalyze -License: GNU General Public License v3 - - - - - - - - - - -
Last stable release - -
Development branch - -
+[![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) ## Description LatAnalyze is a C++11 library for statistical data analysis based on bootstrap @@ -164,4 +149,4 @@ Fixes: #### v3.0 Commit `7b4f2884a5e99bbfab4d4bd7623f609a55403c39`. First 'stable' version of LatAnalyze in C++. The v2.0 refers to the [C version](https://github.com/aportelli/LatAnalyze-legacy) and v1.0 to an old undistributed version. -**This version compiles fine on OS X with clang but does have many portability issues to other platforms/compilers, v3.1 is the first real release.** \ No newline at end of file +**This version compiles fine on OS X with clang but does have many portability issues to other platforms/compilers, v3.1 is the first real release.** diff --git a/ci-scripts/install-deps.sh b/ci-scripts/install-deps.sh index 6faf97b..e08c42a 100755 --- a/ci-scripts/install-deps.sh +++ b/ci-scripts/install-deps.sh @@ -1,15 +1,16 @@ #!/usr/bin/env bash -if (( $# != 1 )); then - echo "usage: `basename $0` " 1>&2 +if (( $# != 2 )); then + echo "usage: `basename $0` " 1>&2 exit 1 fi PREFIX=$1 +NTASKS=$2 set -ex mkdir -p local/build for d in gsl nlopt minuit hdf5; do if [ ! -e local/.built.${d} ]; then - ./install-${d}.sh ${PREFIX} + ./install-${d}.sh ${PREFIX} ${NTASKS} fi done diff --git a/ci-scripts/install-gsl.sh b/ci-scripts/install-gsl.sh index fe39df2..41ed0c2 100755 --- a/ci-scripts/install-gsl.sh +++ b/ci-scripts/install-gsl.sh @@ -2,11 +2,12 @@ NAME='gsl-2.6' -if (( $# != 1 )); then - echo "usage: `basename $0` " 1>&2 +if (( $# != 2 )); then + echo "usage: `basename $0` " 1>&2 exit 1 fi PREFIX=$1 +NTASKS=$2 set -ex INITDIR=$(pwd -P) @@ -19,7 +20,7 @@ tar -xzvf ${NAME}.tar.gz mkdir -p ${NAME}/build cd ${NAME}/build ../configure --prefix=${PREFIX} -make -j4 +make -j${NTASKS} make install cd ${INITDIR}/local touch .built.gsl diff --git a/ci-scripts/install-hdf5.sh b/ci-scripts/install-hdf5.sh index 1788268..77beb3e 100755 --- a/ci-scripts/install-hdf5.sh +++ b/ci-scripts/install-hdf5.sh @@ -1,12 +1,13 @@ #!/usr/bin/env bash -NAME='hdf5-1.10.5' +NAME='hdf5-1.10.8' -if (( $# != 1 )); then - echo "usage: `basename $0` " 1>&2 +if (( $# != 2 )); then + echo "usage: `basename $0` " 1>&2 exit 1 fi PREFIX=$1 +NTASKS=$2 set -ex INITDIR=$(pwd -P) @@ -19,7 +20,7 @@ tar -xzvf ${NAME}.tar.gz mkdir ${NAME}/build cd ${NAME}/build ../configure --prefix=${PREFIX} --enable-cxx -make -j4 +make -j${NTASKS} make install cd ${INITDIR}/local touch .built.hdf5 diff --git a/ci-scripts/install-latan.sh b/ci-scripts/install-latan.sh index 797b07f..3e8a377 100755 --- a/ci-scripts/install-latan.sh +++ b/ci-scripts/install-latan.sh @@ -1,10 +1,11 @@ #!/usr/bin/env bash -if (( $# != 1 )); then - echo "usage: `basename $0` " 1>&2 +if (( $# != 2 )); then + echo "usage: `basename $0` " 1>&2 exit 1 fi PREFIX=$1 +NTASKS=$2 set -ex INITDIR=$(pwd -P) @@ -12,11 +13,11 @@ mkdir -p ${PREFIX} cd ${PREFIX} PREFIX=$(pwd -P) cd ${INITDIR} -./install-deps.sh ${PREFIX} +./install-deps.sh ${PREFIX} ${NTASKS} cd .. ./bootstrap.sh mkdir -p build cd build -../configure --prefix=${PREFIX} --with-minuit=${PREFIX} --with-nlopt=${PREFIX} --with-hdf5=${PREFIX} --with-gsl=${PREFIX} CXXFLAGS="${CXXFLAGS} -O3 -march=haswell -mtune=haswell" -make -j4 +../configure --prefix=${PREFIX} --with-minuit=${PREFIX} --with-nlopt=${PREFIX} --with-hdf5=${PREFIX} --with-gsl=${PREFIX} CXXFLAGS="${CXXFLAGS} -O3 -march=native -mtune=native" +make -j${NTASKS} make install diff --git a/ci-scripts/install-minuit.sh b/ci-scripts/install-minuit.sh index 3c15118..02e96bb 100755 --- a/ci-scripts/install-minuit.sh +++ b/ci-scripts/install-minuit.sh @@ -1,12 +1,12 @@ #!/usr/bin/env bash -NAME='Minuit2-5.34.14' -if (( $# != 1 )); then - echo "usage: `basename $0` " 1>&2 +if (( $# != 2 )); then + echo "usage: `basename $0` " 1>&2 exit 1 fi PREFIX=$1 +NTASKS=$2 set -ex INITDIR=$(pwd -P) @@ -14,12 +14,12 @@ mkdir -p ${PREFIX} cd ${PREFIX} PREFIX=$(pwd -P) cd ${INITDIR}/local/build -wget http://www.cern.ch/mathlibs/sw/5_34_14/Minuit2/${NAME}.tar.gz -tar -xzvf ${NAME}.tar.gz -mkdir -p ${NAME}/build -cd ${NAME}/build -../configure --prefix=${PREFIX} --disable-openmp -make -j4 +rm -rf root +git clone https://github.com/root-project/root.git +cd root/math/minuit2/ +mkdir build; cd build +cmake .. -Dminuit2_standalone=ON -DCMAKE_INSTALL_PREFIX=${PREFIX} +make -j${NTASKS} make install cd ${INITDIR}/local touch .built.minuit diff --git a/ci-scripts/install-nlopt.sh b/ci-scripts/install-nlopt.sh index 3851888..23b1572 100755 --- a/ci-scripts/install-nlopt.sh +++ b/ci-scripts/install-nlopt.sh @@ -2,11 +2,12 @@ NAME='2.6.1' -if (( $# != 1 )); then - echo "usage: `basename $0` " 1>&2 +if (( $# != 2 )); then + echo "usage: `basename $0` " 1>&2 exit 1 fi PREFIX=$1 +NTASKS=$2 set -ex INITDIR=$(pwd -P) @@ -20,7 +21,7 @@ NAME=nlopt-${NAME} mkdir -p ${NAME}/build cd ${NAME}/build cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} -DCMAKE_BUILD_WITH_INSTALL_NAME_DIR=TRUE -DCMAKE_INSTALL_NAME_DIR="${PREFIX}/lib" .. -make -j4 +make -j${NTASKS} make install cd ${INITDIR}/local touch .built.nlopt diff --git a/configure.ac b/configure.ac index 17b3aed..27538e5 100644 --- a/configure.ac +++ b/configure.ac @@ -36,7 +36,7 @@ AC_ARG_WITH([gsl], AC_ARG_WITH([minuit], [AS_HELP_STRING([--with-minuit=prefix], [try this for a non-standard install prefix of the Minuit2 library])], - [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_minuit/include"] + [AM_CXXFLAGS="$AM_CXXFLAGS -I$with_minuit/include -I$with_minuit/include/Minuit2 -I$with_minuit/include/Fit"] [AM_LDFLAGS="$AM_LDFLAGS -L$with_minuit/lib"]) AC_ARG_WITH([nlopt], [AS_HELP_STRING([--with-nlopt=prefix], @@ -74,6 +74,7 @@ CXXFLAGS_CPY=$CXXFLAGS LDFLAGS_CPY=$LDFLAGS CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS" LDFLAGS="$AM_LDFLAGS $LDFLAGS" +AC_CHECK_LIB([pthread],[pthread_create],[],[AC_MSG_ERROR([pthread library not found])]) AC_CHECK_LIB([m],[cos],[],[AC_MSG_ERROR([libm library not found])]) AC_CHECK_LIB([gslcblas],[cblas_dgemm],[], [AC_MSG_ERROR([GSL CBLAS library not found])]) @@ -90,10 +91,10 @@ AC_CHECK_LIB([hdf5_cpp],[H5Fopen], [AC_MSG_ERROR([HDF5 library not found])], [-lhdf5]) SAVED_LDFLAGS=$LDFLAGS LDFLAGS="$LDFLAGS -lMinuit2" -AC_MSG_CHECKING([for ROOT::Minuit2::BasicMinimumError in -lMinuit2]); +AC_MSG_CHECKING([for ROOT::Minuit2::VariableMetricMinimizer in -lMinuit2]); AC_LINK_IFELSE( - [AC_LANG_PROGRAM([#include ], - [ROOT::Minuit2::BasicMinimumError dummy(0)])], + [AC_LANG_PROGRAM([#include ], + [ROOT::Minuit2::VariableMetricMinimizer dummy()])], [LIBS="$LIBS -lMinuit2"] [AC_DEFINE([HAVE_MINUIT2], [1], @@ -103,6 +104,20 @@ AC_LINK_IFELSE( [have_minuit=false] [AC_MSG_RESULT([no])]) AM_CONDITIONAL([HAVE_MINUIT], [test x$have_minuit = xtrue]) +LDFLAGS="$LDFLAGS -lMinuit2Math" +AC_MSG_CHECKING([for ROOT::Math::MinimizerOptions in -lMinuit2Math]); +AC_LINK_IFELSE( + [AC_LANG_PROGRAM([#include ], + [ROOT::Math::MinimizerOptions dummy()])], + [LIBS="$LIBS -lMinuit2Math"] + [AC_DEFINE([HAVE_MINUIT2MATH], + [1], + [Define to 1 if you have the `Minuit2Math' library (-lMinuit2Math).])] + [have_minuitmath=true] + [AC_MSG_RESULT([yes])], + [have_minuitmath=false] + [AC_MSG_RESULT([no])]) +AM_CONDITIONAL([HAVE_MINUITMATH], [test x$have_minuit = xtrue]) LDFLAGS=$SAVED_LDFLAGS CXXFLAGS=$CXXFLAGS_CPY LDFLAGS=$LDFLAGS_CPY diff --git a/examples/Makefile.am b/examples/Makefile.am index dcb5a7c..123bafb 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -19,7 +19,8 @@ noinst_PROGRAMS = \ exPlot \ exPValue \ exRand \ - exRootFinder + exRootFinder \ + exThreadPool exCompiledDoubleFunction_SOURCES = exCompiledDoubleFunction.cpp exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS) @@ -73,4 +74,8 @@ exRootFinder_SOURCES = exRootFinder.cpp exRootFinder_CXXFLAGS = $(COM_CXXFLAGS) exRootFinder_LDFLAGS = -L../lib/.libs -lLatAnalyze +exThreadPool_SOURCES = exThreadPool.cpp +exThreadPool_CXXFLAGS = $(COM_CXXFLAGS) +exThreadPool_LDFLAGS = -L../lib/.libs -lLatAnalyze + ACLOCAL_AMFLAGS = -I .buildutils/m4 diff --git a/examples/exThreadPool.cpp b/examples/exThreadPool.cpp new file mode 100644 index 0000000..f22ea43 --- /dev/null +++ b/examples/exThreadPool.cpp @@ -0,0 +1,29 @@ +#include + +using namespace std; +using namespace Latan; + +int main(void) +{ + ThreadPool pool; + + cout << "Using " << pool.getThreadNum() << " threads" << endl; + for (unsigned int i = 1; i <= 20; ++i) + { + pool.addJob([i, &pool](void) + { + pool.critical([i](void) + { + cout << "job " << i << " wait for " << i*100 << " ms" << endl; + }); + this_thread::sleep_for(chrono::milliseconds(i*100)); + pool.critical([i](void) + { + cout << "job " << i << " done" << endl; + }); + }); + } + pool.terminate(); + + return EXIT_SUCCESS; +} diff --git a/lib/Core/EigenPlugin.hpp b/lib/Core/EigenPlugin.hpp index 48bd820..ee241f2 100644 --- a/lib/Core/EigenPlugin.hpp +++ b/lib/Core/EigenPlugin.hpp @@ -17,7 +17,7 @@ * along with LatAnalyze. If not, see . */ -Derived pInverse(const double tolerance = 1.0e-10) +Derived pInverse(const double tolerance = 1.0e-10) const { auto svd = jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV); const auto u = svd.matrixU(); @@ -52,7 +52,7 @@ Derived pInverse(const double tolerance = 1.0e-10) return v*s.asDiagonal()*u.transpose(); } -Derived singularValues(void) +Derived singularValues(void) const { auto svd = jacobiSvd(); diff --git a/lib/Core/Math.cpp b/lib/Core/Math.cpp index df9e8ba..5278412 100644 --- a/lib/Core/Math.cpp +++ b/lib/Core/Math.cpp @@ -29,7 +29,8 @@ using namespace Latan; ******************************************************************************/ DMat MATH_NAMESPACE::varToCorr(const DMat &var) { - DMat res = var, invDiag = res.diagonal(); + DMat res = var; + DVec invDiag = res.diagonal(); invDiag = invDiag.cwiseInverse().cwiseSqrt(); res = (invDiag*invDiag.transpose()).cwiseProduct(res); @@ -37,6 +38,28 @@ DMat MATH_NAMESPACE::varToCorr(const DMat &var) return res; } +DMat MATH_NAMESPACE::corrToVar(const DMat &corr, const DVec &varDiag) +{ + DMat res = corr; + DVec varSqrtDiag = varDiag.cwiseSqrt(); + + res = (varSqrtDiag*varSqrtDiag.transpose()).cwiseProduct(res); + + return res; +} + +double MATH_NAMESPACE::svdDynamicRange(const DMat &mat) +{ + DVec s = mat.singularValues(); + + return s.maxCoeff()/s.minCoeff(); +} + +double MATH_NAMESPACE::svdDynamicRangeDb(const DMat &mat) +{ + return 10.*log10(svdDynamicRange(mat)); +} + /****************************************************************************** * Standard C functions * ******************************************************************************/ diff --git a/lib/Core/Math.hpp b/lib/Core/Math.hpp index 3802ceb..0ece91a 100644 --- a/lib/Core/Math.hpp +++ b/lib/Core/Math.hpp @@ -70,6 +70,11 @@ namespace MATH_NAMESPACE // convert variance matrix to correlation matrix DMat varToCorr(const DMat &var); + DMat corrToVar(const DMat &corr, const DVec &varDiag); + + // matrix SVD dynamic range + double svdDynamicRange(const DMat &mat); + double svdDynamicRangeDb(const DMat &mat); // Constants constexpr double pi = 3.1415926535897932384626433832795028841970; diff --git a/lib/Core/Plot.cpp b/lib/Core/Plot.cpp index d7b44fa..42848ed 100644 --- a/lib/Core/Plot.cpp +++ b/lib/Core/Plot.cpp @@ -112,7 +112,7 @@ PlotHeadCommand::PlotHeadCommand(const string &command) } // PlotData constructor //////////////////////////////////////////////////////// -PlotData::PlotData(const DMatSample &x, const DMatSample &y) +PlotData::PlotData(const DMatSample &x, const DMatSample &y, const bool abs) { if (x[central].rows() != y[central].rows()) { @@ -128,10 +128,17 @@ PlotData::PlotData(const DMatSample &x, const DMatSample &y) d.col(3) = y.variance().cwiseSqrt(); tmpFileName = dumpToTmpFile(d); pushTmpFile(tmpFileName); - setCommand("'" + tmpFileName + "' u 1:3:2:4 w xyerr"); + if (!abs) + { + setCommand("'" + tmpFileName + "' u 1:3:2:4 w xyerr"); + } + else + { + setCommand("'" + tmpFileName + "' u 1:(abs($3)):2:4 w xyerr"); + } } -PlotData::PlotData(const DVec &x, const DMatSample &y) +PlotData::PlotData(const DVec &x, const DMatSample &y, const bool abs) { if (x.rows() != y[central].rows()) { @@ -146,10 +153,17 @@ PlotData::PlotData(const DVec &x, const DMatSample &y) d.col(2) = y.variance().cwiseSqrt(); tmpFileName = dumpToTmpFile(d); pushTmpFile(tmpFileName); - setCommand("'" + tmpFileName + "' u 1:2:3 w yerr"); + if (!abs) + { + setCommand("'" + tmpFileName + "' u 1:2:3 w yerr"); + } + else + { + setCommand("'" + tmpFileName + "' u 1:(abs($2)):3 w yerr"); + } } -PlotData::PlotData(const DMatSample &x, const DVec &y) +PlotData::PlotData(const DMatSample &x, const DVec &y, const bool abs) { if (x[central].rows() != y.rows()) { @@ -164,14 +178,29 @@ PlotData::PlotData(const DMatSample &x, const DVec &y) d.col(1) = x.variance().cwiseSqrt(); tmpFileName = dumpToTmpFile(d); pushTmpFile(tmpFileName); - setCommand("'" + tmpFileName + "' u 1:3:2 w xerr"); + if (!abs) + { + setCommand("'" + tmpFileName + "' u 1:3:2 w xerr"); + } + else + { + setCommand("'" + tmpFileName + "' u 1:($3):2 w xerr"); + } } -PlotData::PlotData(const XYStatData &data, const Index i, const Index j) +PlotData::PlotData(const XYStatData &data, const Index i, const Index j, const bool abs) { string usingCmd, tmpFileName; - usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr"; + if (!abs) + { + usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr"; + } + else + { + usingCmd = (data.isXExact(i)) ? "u 1:(abs($3)):4 w yerr" : "u 1:(abs($3)):2:4 w xyerr"; + } + tmpFileName = dumpToTmpFile(data.getTable(i, j)); pushTmpFile(tmpFileName); setCommand("'" + tmpFileName + "' " + usingCmd); @@ -195,6 +224,24 @@ PlotLine::PlotLine(const DVec &x, const DVec &y) setCommand("'" + tmpFileName + "' u 1:2 w lines"); } +// PlotPoints constructor //////////////////////////////////////////////////////// +PlotPoints::PlotPoints(const DVec &x, const DVec &y) +{ + if (x.size() != y.size()) + { + LATAN_ERROR(Size, "x and y vectors do not have the same size"); + } + + DMat d(x.size(), 2); + string usingCmd, tmpFileName; + + d.col(0) = x; + d.col(1) = y; + tmpFileName = dumpToTmpFile(d); + pushTmpFile(tmpFileName); + setCommand("'" + tmpFileName + "' u 1:2"); +} + // PlotHLine constructor /////////////////////////////////////////////////////// PlotHLine::PlotHLine(const double y) { @@ -217,7 +264,8 @@ PlotBand::PlotBand(const double xMin, const double xMax, const double yMin, // PlotFunction constructor //////////////////////////////////////////////////// PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin, - const double xMax, const unsigned int nPoint) + const double xMax, const unsigned int nPoint, + const bool abs) { DMat d(nPoint, 2); string tmpFileName; @@ -230,7 +278,14 @@ PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin, } tmpFileName = dumpToTmpFile(d); pushTmpFile(tmpFileName); - setCommand("'" + tmpFileName + "' u 1:2 w lines"); + if (!abs) + { + setCommand("'" + tmpFileName + "' u 1:2 w lines"); + } + else + { + setCommand("'" + tmpFileName + "' u 1:(abs($2)) w lines"); + } } // PlotPredBand constructor //////////////////////////////////////////////////// diff --git a/lib/Core/Plot.hpp b/lib/Core/Plot.hpp index e079906..1c2901c 100644 --- a/lib/Core/Plot.hpp +++ b/lib/Core/Plot.hpp @@ -89,10 +89,11 @@ class PlotData: public PlotObject { public: // constructor - PlotData(const DMatSample &x, const DMatSample &y); - PlotData(const DVec &x, const DMatSample &y); - PlotData(const DMatSample &x, const DVec &y); - PlotData(const XYStatData &data, const Index i = 0, const Index j = 0); + PlotData(const DMatSample &x, const DMatSample &y, const bool abs = false); + PlotData(const DVec &x, const DMatSample &y, const bool abs = false); + PlotData(const DMatSample &x, const DVec &y, const bool abs = false); + PlotData(const XYStatData &data, const Index i = 0, const Index j = 0, + const bool abs = false); // destructor virtual ~PlotData(void) = default; }; @@ -115,6 +116,15 @@ public: virtual ~PlotLine(void) = default; }; +class PlotPoints: public PlotObject +{ +public: + // constructor + PlotPoints(const DVec &x, const DVec &y); + // destructor + virtual ~PlotPoints(void) = default; +}; + class PlotBand: public PlotObject { public: @@ -130,7 +140,8 @@ class PlotFunction: public PlotObject public: // constructor PlotFunction(const DoubleFunction &function, const double xMin, - const double xMax, const unsigned int nPoint = 1000); + const double xMax, const unsigned int nPoint = 1000, + const bool abs = false); // destructor virtual ~PlotFunction(void) = default; }; @@ -182,6 +193,11 @@ PlotRange(Axis::x, -.5, (m).cols() - .5) <<\ PlotRange(Axis::y, (m).rows() - .5, -.5) <<\ PlotMatrixNoRange(m) +#define PlotCorrMatrix(m)\ +PlotHeadCommand("set cbrange [-1:1]") <<\ +PlotHeadCommand("set palette defined (0 'blue', 1 'white', 2 'red')") <<\ +PlotMatrix(m) + /****************************************************************************** * Plot modifiers * ******************************************************************************/ diff --git a/lib/Core/ThreadPool.cpp b/lib/Core/ThreadPool.cpp new file mode 100644 index 0000000..fa8e27e --- /dev/null +++ b/lib/Core/ThreadPool.cpp @@ -0,0 +1,117 @@ +/* + * ThreadPool.cpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2021 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; + +/****************************************************************************** + * ThreadPool implementation * + ******************************************************************************/ +// constructors //////////////////////////////////////////////////////////////// +ThreadPool::ThreadPool(void) +: ThreadPool(std::thread::hardware_concurrency()) +{} + +ThreadPool::ThreadPool(const unsigned int nThreads) +: nThreads_(nThreads) +{ + for (unsigned int t = 0; t < nThreads_; ++t) + { + threads_.push_back(thread(&ThreadPool::workerLoop, this)); + } +} + +// destructor ////////////////////////////////////////////////////////////////// +ThreadPool::~ThreadPool(void) +{ + terminate(); +} + +// get the number of threads /////////////////////////////////////////////////// +unsigned int ThreadPool::getThreadNum(void) const +{ + return nThreads_; +} + +// get the pool mutex for synchronisation ////////////////////////////////////// +std::mutex & ThreadPool::getMutex(void) +{ + return mutex_; +} + +// worker loop ///////////////////////////////////////////////////////////////// +void ThreadPool::workerLoop(void) +{ + while (true) + { + Job job; + { + unique_lock lock(mutex_); + + condition_.wait(lock, [this](){ + return !queue_.empty() || terminatePool_; + }); + if (terminatePool_ and queue_.empty()) + { + return; + } + job = queue_.front(); + queue_.pop(); + } + job(); + } +} + +// add jobs //////////////////////////////////////////////////////////////////// +void ThreadPool::addJob(Job newJob) +{ + { + unique_lock lock(mutex_); + + queue_.push(newJob); + } + condition_.notify_one(); +} + +// critical section //////////////////////////////////////////////////////////// +void ThreadPool::critical(Job fn) +{ + unique_lock lock(mutex_); + + fn(); +} + +// wait for completion ///////////////////////////////////////////////////////// +void ThreadPool::terminate(void) +{ + { + unique_lock lock(mutex_); + + terminatePool_ = true; + } + condition_.notify_all(); + for (auto &thread: threads_) + { + thread.join(); + } + threads_.clear(); +} diff --git a/lib/Core/ThreadPool.hpp b/lib/Core/ThreadPool.hpp new file mode 100644 index 0000000..6594a8e --- /dev/null +++ b/lib/Core/ThreadPool.hpp @@ -0,0 +1,56 @@ +/* + * ThreadPool.hpp, part of LatAnalyze 3 + * + * Copyright (C) 2013 - 2021 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 . + */ + +#ifndef Latan_ThreadPool_hpp_ +#define Latan_ThreadPool_hpp_ + +#include + +class ThreadPool +{ +public: + typedef std::function Job; +public: + // constructors/destructor + ThreadPool(void); + ThreadPool(const unsigned int nThreads); + virtual ~ThreadPool(void); + // get the number of threads + unsigned int getThreadNum(void) const; + // get the pool mutex for synchronisation + std::mutex & getMutex(void); + // add jobs + void addJob(Job newJob); + // critical section + void critical(Job fn); + // wait for completion and terminate + void terminate(void); +private: + // worker loop + void workerLoop(void); +private: + unsigned int nThreads_; + std::condition_variable condition_; + std::vector threads_; + bool terminatePool_{false}; + std::queue queue_; + std::mutex mutex_; +}; + +#endif diff --git a/lib/Core/stdincludes.hpp b/lib/Core/stdincludes.hpp index bec9bf8..e976118 100644 --- a/lib/Core/stdincludes.hpp +++ b/lib/Core/stdincludes.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +41,7 @@ #include #include #include +#include #include #include #include diff --git a/lib/Makefile.am b/lib/Makefile.am index ddb23b6..9fb5e8c 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -32,6 +32,7 @@ libLatAnalyze_la_SOURCES = \ Core/MathParser.ypp \ Core/OptParser.cpp \ Core/Plot.cpp \ + Core/ThreadPool.cpp \ Core/Utilities.cpp \ Functional/CompiledFunction.cpp \ Functional/CompiledModel.cpp \ @@ -75,6 +76,7 @@ HPPFILES = \ Core/OptParser.hpp \ Core/ParserState.hpp \ Core/Plot.hpp \ + Core/ThreadPool.hpp \ Core/stdincludes.hpp \ Core/Utilities.hpp \ Functional/CompiledFunction.hpp \ diff --git a/lib/Numerical/MinuitMinimizer.cpp b/lib/Numerical/MinuitMinimizer.cpp index c01ff19..bab7638 100644 --- a/lib/Numerical/MinuitMinimizer.cpp +++ b/lib/Numerical/MinuitMinimizer.cpp @@ -19,6 +19,20 @@ #include #include + +// forward declaration necessary in the ROOT-based version of Minuit2 +namespace ROOT +{ + namespace Fit + { + class ParameterSettings; + }; +}; + +// macros necessary in the ROOT-based version of Minuit2 +#define ROOT_Math_VecTypes +#define MATHCORE_STANDALONE + #include #include diff --git a/lib/Statistics/MatSample.hpp b/lib/Statistics/MatSample.hpp index a48db9d..a103215 100644 --- a/lib/Statistics/MatSample.hpp +++ b/lib/Statistics/MatSample.hpp @@ -103,6 +103,10 @@ public: const Index nCol); // resize all matrices void resizeMat(const Index nRow, const Index nCol); + // covariance matrix + Mat covarianceMatrix(const MatSample &sample) const; + Mat varianceMatrix(void) const; + Mat correlationMatrix(void) const; }; // non-member operators @@ -379,6 +383,78 @@ void MatSample::resizeMat(const Index nRow, const Index nCol) } } +// covariance matrix /////////////////////////////////////////////////////////// +template +Mat MatSample::covarianceMatrix(const MatSample &sample) const +{ + if (((*this)[central].cols() != 1) or (sample[central].cols() != 1)) + { + LATAN_ERROR(Size, "samples have more than one column"); + } + + Index n1 = (*this)[central].rows(), n2 = sample[central].rows(); + Index nSample = this->size(); + Mat tmp1(n1, nSample), tmp2(n2, nSample), res(n1, n2); + Mat s1(n1, 1), s2(n2, 1), one(nSample, 1); + + one.fill(1.); + s1.fill(0.); + s2.fill(0.); + for (unsigned int s = 0; s < nSample; ++s) + { + s1 += (*this)[s]; + tmp1.col(s) = (*this)[s]; + } + tmp1 -= s1*one.transpose()/static_cast(nSample); + for (unsigned int s = 0; s < nSample; ++s) + { + s2 += sample[s]; + tmp2.col(s) = sample[s]; + } + tmp2 -= s2*one.transpose()/static_cast(nSample); + res = tmp1*tmp2.transpose()/static_cast(nSample - 1); + + return res; +} + +template +Mat MatSample::varianceMatrix(void) const +{ + if ((*this)[central].cols() != 1) + { + LATAN_ERROR(Size, "samples have more than one column"); + } + + Index n1 = (*this)[central].rows(); + Index nSample = this->size(); + Mat tmp1(n1, nSample), res(n1, n1); + Mat s1(n1, 1), one(nSample, 1); + + one.fill(1.); + s1.fill(0.); + for (unsigned int s = 0; s < nSample; ++s) + { + s1 += (*this)[s]; + tmp1.col(s) = (*this)[s]; + } + tmp1 -= s1*one.transpose()/static_cast(nSample); + res = tmp1*tmp1.transpose()/static_cast(nSample - 1); + + return res; +} + +template +Mat MatSample::correlationMatrix(void) const +{ + Mat res = varianceMatrix(); + Mat invDiag(res.rows(), 1); + + invDiag = res.diagonal(); + invDiag = invDiag.cwiseInverse().cwiseSqrt(); + res = (invDiag*invDiag.transpose()).cwiseProduct(res); + + return res; +} END_LATAN_NAMESPACE diff --git a/lib/Statistics/StatArray.hpp b/lib/Statistics/StatArray.hpp index 9637089..4e70dbb 100644 --- a/lib/Statistics/StatArray.hpp +++ b/lib/Statistics/StatArray.hpp @@ -51,14 +51,11 @@ public: const T & operator[](const Index s) const; // statistics void bin(Index binSize); + T sum(const Index pos = 0, const Index n = -1) const; + T meanOld(const Index pos = 0, const Index n = -1) const; T mean(const Index pos = 0, const Index n = -1) const; - T covariance(const StatArray &array, const Index pos = 0, - const Index n = -1) const; - T covarianceMatrix(const StatArray &array, const Index pos = 0, - const Index n = -1) const; - T variance(const Index pos = 0, const Index n = -1) const; - T varianceMatrix(const Index pos = 0, const Index n = -1) const; - T correlationMatrix(const Index pos = 0, const Index n = -1) const; + T covariance(const StatArray &array) const; + T variance(void) const; // IO type virtual IoType getType(void) const; public: @@ -66,7 +63,7 @@ public: }; // reduction operations -namespace ReducOp +namespace StatOp { // general templates template @@ -148,128 +145,67 @@ void StatArray::bin(Index binSize) } } +template +T StatArray::sum(const Index pos, const Index n) const +{ + T result; + const Index m = (n >= 0) ? n : size(); + + result = (*this)[pos]; + for (Index i = pos + 1; i < pos + m; ++i) + { + result += (*this)[i]; + } + + return result; +} + template T StatArray::mean(const Index pos, const Index n) const { - T result = T(); const Index m = (n >= 0) ? n : size(); - if (m) + return sum(pos, n)/static_cast(m); +} + +template +T StatArray::covariance(const StatArray &array) const +{ + T s1, s2, res; + + s1 = array.sum(); + s2 = this->sum(); + res = StatOp::prod(array[0], (*this)[0]); + for (Index i = 1; i < size(); ++i) { - result = this->segment(pos+os, m).redux(&ReducOp::sum); + res += StatOp::prod(array[i], (*this)[i]); } - return result/static_cast(m); -} - -template -T StatArray::covariance(const StatArray &array, const Index pos, - const Index n) const -{ - T s1, s2, prs, res = T(); - const Index m = (n >= 0) ? n : size(); - - if (m) - { - auto arraySeg = array.segment(pos+os, m); - auto thisSeg = this->segment(pos+os, m); - - s1 = thisSeg.redux(&ReducOp::sum); - s2 = arraySeg.redux(&ReducOp::sum); - prs = thisSeg.binaryExpr(arraySeg, &ReducOp::prod) - .redux(&ReducOp::sum); - res = prs - ReducOp::prod(s1, s2)/static_cast(m); - } - - return res/static_cast(m - 1); -} - -template -T StatArray::covarianceMatrix(const StatArray &array, - const Index pos, const Index n) const -{ - T s1, s2, prs, res = T(); - const Index m = (n >= 0) ? n : size(); - - if (m) - { - auto arraySeg = array.segment(pos+os, m); - auto thisSeg = this->segment(pos+os, m); - - s1 = thisSeg.redux(&ReducOp::sum); - s2 = arraySeg.redux(&ReducOp::sum); - prs = thisSeg.binaryExpr(arraySeg, &ReducOp::tensProd) - .redux(&ReducOp::sum); - res = prs - ReducOp::tensProd(s1, s2)/static_cast(m); - } - - return res/static_cast(m - 1); -} - -template -T StatArray::variance(const Index pos, const Index n) const -{ - return covariance(*this, pos, n); -} - -template -T StatArray::varianceMatrix(const Index pos, const Index n) const -{ - return covarianceMatrix(*this, pos, n); -} - -template -T StatArray::correlationMatrix(const Index pos, const Index n) const -{ - T res = varianceMatrix(pos, n); - T invDiag(res.rows(), 1); - - invDiag = res.diagonal(); - invDiag = invDiag.cwiseInverse().cwiseSqrt(); - res = (invDiag*invDiag.transpose()).cwiseProduct(res); + res -= StatOp::prod(s1, s2)/static_cast(size()); + res /= static_cast(size() - 1); return res; } -// reduction operations //////////////////////////////////////////////////////// -namespace ReducOp +template +T StatArray::variance(void) const { - template - inline T sum(const T &a, const T &b) - { - return a + b; - } + return covariance(*this); +} +// reduction operations //////////////////////////////////////////////////////// +namespace StatOp +{ template inline T prod(const T &a, const T &b) { return a*b; } - template - inline T tensProd(const T &v1 __dumb, const T &v2 __dumb) - { - LATAN_ERROR(Implementation, - "tensorial product not implemented for this type"); - } - template <> inline Mat prod(const Mat &a, const Mat &b) { return a.cwiseProduct(b); } - - template <> - inline Mat tensProd(const Mat &v1, - const Mat &v2) - { - if ((v1.cols() != 1) or (v2.cols() != 1)) - { - LATAN_ERROR(Size, - "tensorial product is only valid with column vectors"); - } - - return v1*v2.transpose(); - } } // IO type ///////////////////////////////////////////////////////////////////// diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 48f6d4d..f014c94 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -62,6 +62,11 @@ double SampleFitResult::getPValue(const Index s) const return Math::chi2PValue(getChi2(s), getNDof()); } +double SampleFitResult::getCorrRangeDb(void) const +{ + return corrRangeDb_; +} + double SampleFitResult::getCcdf(const Index s) const { return Math::chi2Ccdf(getChi2(s), getNDof()); @@ -107,9 +112,11 @@ void SampleFitResult::print(const bool printXsi, ostream &out) const getChi2(), static_cast(getNDof()), getChi2PerDof(), getCcdf(), getPValue()); out << buf << endl; + sprintf(buf, "correlation dynamic range= %.1f dB", getCorrRangeDb()); + out << buf << endl; for (Index p = 0; p < pMax; ++p) { - sprintf(buf, "%8s= % e +/- %e", parName_[p].c_str(), + sprintf(buf, "%12s= % e +/- %e", parName_[p].c_str(), (*this)[central](p), err(p)); out << buf << endl; } @@ -249,6 +256,20 @@ const DMat & XYSampleData::getFitVarMatPInv(void) return data_.getFitVarMatPInv(); } +const DMat & XYSampleData::getFitCorrMat(void) +{ + computeVarMat(); + + return data_.getFitCorrMat(); +} + +const DMat & XYSampleData::getFitCorrMatPInv(void) +{ + computeVarMat(); + + return data_.getFitCorrMatPInv(); +} + // set data to a particular sample ///////////////////////////////////////////// void XYSampleData::setDataToSample(const Index s) { @@ -330,9 +351,10 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, { computeVarMat(); - SampleFitResult result; - FitResult sampleResult; - DVec initCopy = init; + SampleFitResult result; + FitResult sampleResult; + DVec initCopy = init; + Minimizer::Verbosity verbCopy = minimizer.back()->getVerbosity(); result.resize(nSample_); result.chi2_.resize(nSample_); @@ -348,9 +370,11 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.model_[j][s] = sampleResult.getModel(j); } } - result.nPar_ = sampleResult.getNPar(); - result.nDof_ = sampleResult.nDof_; - result.parName_ = sampleResult.parName_; + minimizer.back()->setVerbosity(verbCopy); + result.nPar_ = sampleResult.getNPar(); + result.nDof_ = sampleResult.nDof_; + result.parName_ = sampleResult.parName_; + result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); return result; } @@ -382,6 +406,29 @@ XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit) return res; } +XYSampleData XYSampleData::getNormalisedResiduals(const SampleFitResult &fit) +{ + XYSampleData res(*this); + + for (Index j = 0; j < getNYDim(); ++j) + { + const DoubleFunctionSample &f = fit.getModel(_, j); + + for (auto &p: yData_[j]) + { + res.y(p.first, j) -= f(x(p.first)); + } + + const DMat &var = res.getYYVar(j, j); + for (auto &p: yData_[j]) + { + res.y(p.first, j) /= sqrt(var(p.first, p.first)); + } + } + + return res; +} + XYSampleData XYSampleData::getPartialResiduals(const SampleFitResult &fit, const DVec &ref, const Index i) { diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index 59e9bb7..32e4cbe 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -49,6 +49,7 @@ public: double getNDof(void) const; Index getNPar(void) const; double getPValue(const Index s = central) const; + double getCorrRangeDb(void) const; double getCcdf(const Index s = central) const; const DoubleFunction & getModel(const Index s = central, const Index j = 0) const; @@ -60,6 +61,7 @@ public: std::ostream &out = std::cout) const; private: DSample chi2_; + double corrRangeDb_{0.}; Index nDof_{0}, nPar_{0}; std::vector model_; std::vector parName_; @@ -91,9 +93,11 @@ public: const DMat & getXYVar(const Index i, const Index j); DVec getXError(const Index i); DVec getYError(const Index j); - // get total fit variance matrix and its pseudo-inverse + // get total fit variance & correlation matrices and their pseudo-inverse const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); + const DMat & getFitCorrMat(void); + const DMat & getFitCorrMatPInv(void); // set data to a particular sample void setDataToSample(const Index s); // get internal XYStatData @@ -118,6 +122,7 @@ public: const DoubleModel &model, const Ts... models); // residuals XYSampleData getResiduals(const SampleFitResult &fit); + XYSampleData getNormalisedResiduals(const SampleFitResult &fit); XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x, const Index i); private: diff --git a/lib/Statistics/XYStatData.cpp b/lib/Statistics/XYStatData.cpp index d324842..690a3a6 100644 --- a/lib/Statistics/XYStatData.cpp +++ b/lib/Statistics/XYStatData.cpp @@ -60,6 +60,11 @@ double FitResult::getCcdf(void) const return Math::chi2Ccdf(getChi2(), getNDof());; } +double FitResult::getCorrRangeDb(void) const +{ + return corrRangeDb_; +} + const DoubleFunction & FitResult::getModel(const Index j) const { return model_[j]; @@ -75,9 +80,11 @@ void FitResult::print(const bool printXsi, ostream &out) const getChi2(), static_cast(getNDof()), getChi2PerDof(), getCcdf(), getPValue()); out << buf << endl; + sprintf(buf, "correlation dynamic range= %.1f dB", getCorrRangeDb()); + out << buf << endl; for (Index p = 0; p < pMax; ++p) { - sprintf(buf, "%8s= %e", parName_[p].c_str(), (*this)(p)); + sprintf(buf, "%12s= %e", parName_[p].c_str(), (*this)(p)); out << buf << endl; } } @@ -216,7 +223,7 @@ DVec XYStatData::getXError(const Index i) const DVec XYStatData::getYError(const Index j) const { - checkXDim(j); + checkYDim(j); return yyVar_(j, j).diagonal().cwiseSqrt(); } @@ -259,6 +266,20 @@ const DMat & XYStatData::getFitVarMatPInv(void) return fitVarInv_; } +const DMat & XYStatData::getFitCorrMat(void) +{ + updateFitVarMat(); + + return fitCorr_; +} + +const DMat & XYStatData::getFitCorrMatPInv(void) +{ + updateFitVarMat(); + + return fitCorrInv_; +} + // fit ///////////////////////////////////////////////////////////////////////// FitResult XYStatData::fit(vector &minimizer, const DVec &init, const vector &v) @@ -337,9 +358,10 @@ FitResult XYStatData::fit(vector &minimizer, const DVec &init, result = (*m)(chi2); totalInit = result; } - result.chi2_ = chi2(result); - result.nPar_ = nPar; - result.nDof_ = layout.totalYSize - nPar; + result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); + result.chi2_ = chi2(result); + result.nPar_ = nPar; + result.nDof_ = layout.totalYSize - nPar; result.model_.resize(v.size()); for (unsigned int j = 0; j < v.size(); ++j) { @@ -379,6 +401,27 @@ XYStatData XYStatData::getResiduals(const FitResult &fit) return res; } +XYStatData XYStatData::getNormalisedResiduals(const FitResult &fit) +{ + XYStatData res(*this); + + for (Index j = 0; j < getNYDim(); ++j) + { + const DoubleFunction &f = fit.getModel(j); + const DVec err = getYError(j); + Index row = 0; + + for (auto &p: yData_[j]) + { + res.y(p.first, j) -= f(x(p.first)); + res.y(p.first, j) /= err(row); + row++; + } + } + + return res; +} + XYStatData XYStatData::getPartialResiduals(const FitResult &fit, const DVec &ref, const Index i) { @@ -530,8 +573,11 @@ void XYStatData::updateFitVarMat(void) chi2DataVec_.resize(layout.totalSize); chi2ModVec_.resize(layout.totalSize); chi2Vec_.resize(layout.totalSize); - fitVar_ = fitVar_.cwiseProduct(makeCorrFilter()); - fitVarInv_ = fitVar_.pInverse(getSvdTolerance()); + fitVar_ = fitVar_.cwiseProduct(makeCorrFilter()); + fitCorr_ = Math::varToCorr(fitVar_); + fitCorrInv_ = fitCorr_.pInverse(getSvdTolerance()); + fitVarInv_ = Math::corrToVar(fitCorrInv_, fitVar_.diagonal().cwiseInverse()); + scheduleFitVarMatInit(false); } } diff --git a/lib/Statistics/XYStatData.hpp b/lib/Statistics/XYStatData.hpp index e14b15c..49cbd28 100644 --- a/lib/Statistics/XYStatData.hpp +++ b/lib/Statistics/XYStatData.hpp @@ -48,12 +48,13 @@ public: Index getNPar(void) const; double getPValue(void) const; double getCcdf(void) const; + double getCorrRangeDb(void) const; const DoubleFunction & getModel(const Index j = 0) const; // IO void print(const bool printXsi = false, std::ostream &out = std::cout) const; private: - double chi2_{0.}; + double chi2_{0.}, corrRangeDb_{0.}; Index nDof_{0}, nPar_{0}; std::vector model_; std::vector parName_; @@ -88,9 +89,11 @@ public: DVec getXError(const Index i) const; DVec getYError(const Index j) const; DMat getTable(const Index i, const Index j) const; - // get total fit variance matrix and its pseudo-inverse + // get total fit variance & correlation matrices and their pseudo-inverse const DMat & getFitVarMat(void); const DMat & getFitVarMatPInv(void); + const DMat & getFitCorrMat(void); + const DMat & getFitCorrMatPInv(void); // fit FitResult fit(std::vector &minimizer, const DVec &init, const std::vector &v); @@ -104,6 +107,7 @@ public: const DoubleModel &model, const Ts... models); // residuals XYStatData getResiduals(const FitResult &fit); + XYStatData getNormalisedResiduals(const FitResult &fit); XYStatData getPartialResiduals(const FitResult &fit, const DVec &ref, const Index i); protected: @@ -130,7 +134,7 @@ private: std::vector xData_; std::vector xMap_; Mat xxVar_, yyVar_, xyVar_; - DMat fitVar_, fitVarInv_; + DMat fitVar_, fitVarInv_, fitCorr_, fitCorrInv_; DVec chi2DataVec_, chi2ModVec_, chi2Vec_; DVec xBuf_; bool initXMap_{true}; diff --git a/utils/sample-plot-corr.cpp b/utils/sample-plot-corr.cpp index b13136e..02316c0 100644 --- a/utils/sample-plot-corr.cpp +++ b/utils/sample-plot-corr.cpp @@ -67,7 +67,8 @@ int main(int argc, char *argv[]) sample = sample.block(0, 0, sample[central].rows(), 1); var = sample.varianceMatrix(); corr = sample.correlationMatrix(); - p << PlotMatrix(corr); + + p << PlotCorrMatrix(corr); p.display(); if (!outVarName.empty()) {