From 68d22eca11aafa3673bffe7cbbee29b2e5253ea2 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Fri, 27 Nov 2020 17:09:31 +0000 Subject: [PATCH 01/23] Badges update --- Readme.md | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) 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.** From b4b6bd22fa9aae8926c13c71d371df8769e8c88b Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 16 Nov 2021 21:34:25 +0000 Subject: [PATCH 02/23] compatible with Minuit2 in ROOT --- ci-scripts/install-minuit.sh | 11 +++++------ configure.ac | 22 ++++++++++++++++++---- lib/Numerical/MinuitMinimizer.cpp | 14 ++++++++++++++ 3 files changed, 37 insertions(+), 10 deletions(-) diff --git a/ci-scripts/install-minuit.sh b/ci-scripts/install-minuit.sh index 3c15118..0f6c716 100755 --- a/ci-scripts/install-minuit.sh +++ b/ci-scripts/install-minuit.sh @@ -1,6 +1,5 @@ #!/usr/bin/env bash -NAME='Minuit2-5.34.14' if (( $# != 1 )); then echo "usage: `basename $0` " 1>&2 @@ -14,11 +13,11 @@ 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 +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 -j4 make install cd ${INITDIR}/local diff --git a/configure.ac b/configure.ac index 17b3aed..d1f2fab 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], @@ -90,10 +90,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 +103,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/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 From 9341a31cf4200c2ede69eff7a714ebb4fb934eaa Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 16 Nov 2021 21:34:39 +0000 Subject: [PATCH 03/23] gitignore update --- .gitignore | 1 + 1 file changed, 1 insertion(+) 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 From d6e5ba724d7767a5a97b4085d29d5258f90b76d9 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 16 Nov 2021 21:35:35 +0000 Subject: [PATCH 04/23] CI scripts expose number of build tasks --- ci-scripts/install-deps.sh | 7 ++++--- ci-scripts/install-gsl.sh | 7 ++++--- ci-scripts/install-hdf5.sh | 7 ++++--- ci-scripts/install-latan.sh | 11 ++++++----- ci-scripts/install-minuit.sh | 7 ++++--- ci-scripts/install-nlopt.sh | 7 ++++--- 6 files changed, 26 insertions(+), 20 deletions(-) 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..c0617e3 100755 --- a/ci-scripts/install-hdf5.sh +++ b/ci-scripts/install-hdf5.sh @@ -2,11 +2,12 @@ NAME='hdf5-1.10.5' -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 0f6c716..02e96bb 100755 --- a/ci-scripts/install-minuit.sh +++ b/ci-scripts/install-minuit.sh @@ -1,11 +1,12 @@ #!/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) @@ -18,7 +19,7 @@ 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 -j4 +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 From a3054a0f449ea669efecce3f8ca8748f54e935b4 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 Nov 2021 16:06:26 +0000 Subject: [PATCH 05/23] first GitHub CI --- .github/workflows/build-macos.yml | 25 +++++++++++++++++++++++++ .github/workflows/build-ubuntu.yml | 26 ++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) create mode 100644 .github/workflows/build-macos.yml create mode 100644 .github/workflows/build-ubuntu.yml diff --git a/.github/workflows/build-macos.yml b/.github/workflows/build-macos.yml new file mode 100644 index 0000000..b3068e3 --- /dev/null +++ b/.github/workflows/build-macos.yml @@ -0,0 +1,25 @@ +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: | + 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 From 75485219d8066eb95f5893f201f819b14f205a0e Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 Nov 2021 16:17:45 +0000 Subject: [PATCH 06/23] CI uses HDF5 1.10.8 hoping to solve compilation issues on Big Sur --- ci-scripts/install-hdf5.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ci-scripts/install-hdf5.sh b/ci-scripts/install-hdf5.sh index c0617e3..77beb3e 100755 --- a/ci-scripts/install-hdf5.sh +++ b/ci-scripts/install-hdf5.sh @@ -1,6 +1,6 @@ #!/usr/bin/env bash -NAME='hdf5-1.10.5' +NAME='hdf5-1.10.8' if (( $# != 2 )); then echo "usage: `basename $0` " 1>&2 From ccb837a244c18424b7d45a196074ad1983e1a018 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 17 Nov 2021 16:37:05 +0000 Subject: [PATCH 07/23] Update build-macos.yml --- .github/workflows/build-macos.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/build-macos.yml b/.github/workflows/build-macos.yml index b3068e3..587b668 100644 --- a/.github/workflows/build-macos.yml +++ b/.github/workflows/build-macos.yml @@ -21,5 +21,6 @@ jobs: - 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 From 5f192ad30fcc8a6a2c1718ce566cccf320163759 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Sun, 28 Nov 2021 23:26:17 +0000 Subject: [PATCH 08/23] Thread pool implementation --- examples/Makefile.am | 7 ++- examples/exThreadPool.cpp | 29 ++++++++++ lib/Core/ThreadPool.cpp | 109 ++++++++++++++++++++++++++++++++++++++ lib/Core/ThreadPool.hpp | 54 +++++++++++++++++++ lib/Core/stdincludes.hpp | 2 + lib/Makefile.am | 2 + 6 files changed, 202 insertions(+), 1 deletion(-) create mode 100644 examples/exThreadPool.cpp create mode 100644 lib/Core/ThreadPool.cpp create mode 100644 lib/Core/ThreadPool.hpp 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..85c714d --- /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) + { + { + unique_lock lock(pool.getMutex()); + cout << "job " << i << " wait for " << i*100 << " ms" << endl; + } + this_thread::sleep_for(chrono::milliseconds(i*100)); + { + unique_lock lock(pool.getMutex()); + cout << "job " << i << " done" << endl; + } + }); + } + pool.terminate(); + + return EXIT_SUCCESS; +} diff --git a/lib/Core/ThreadPool.cpp b/lib/Core/ThreadPool.cpp new file mode 100644 index 0000000..de4f575 --- /dev/null +++ b/lib/Core/ThreadPool.cpp @@ -0,0 +1,109 @@ +/* + * 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(); +} + +// 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..5a504e5 --- /dev/null +++ b/lib/Core/ThreadPool.hpp @@ -0,0 +1,54 @@ +/* + * 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); + // 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 \ From 78351a9b76499d8c12056ee96e0fc6cf66fdb501 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Sun, 28 Nov 2021 23:51:56 +0000 Subject: [PATCH 09/23] better interface for critical threaded sections --- examples/exThreadPool.cpp | 8 ++++---- lib/Core/ThreadPool.cpp | 8 ++++++++ lib/Core/ThreadPool.hpp | 2 ++ 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/examples/exThreadPool.cpp b/examples/exThreadPool.cpp index 85c714d..f22ea43 100644 --- a/examples/exThreadPool.cpp +++ b/examples/exThreadPool.cpp @@ -12,15 +12,15 @@ int main(void) { pool.addJob([i, &pool](void) { + pool.critical([i](void) { - unique_lock lock(pool.getMutex()); cout << "job " << i << " wait for " << i*100 << " ms" << endl; - } + }); this_thread::sleep_for(chrono::milliseconds(i*100)); + pool.critical([i](void) { - unique_lock lock(pool.getMutex()); cout << "job " << i << " done" << endl; - } + }); }); } pool.terminate(); diff --git a/lib/Core/ThreadPool.cpp b/lib/Core/ThreadPool.cpp index de4f575..fa8e27e 100644 --- a/lib/Core/ThreadPool.cpp +++ b/lib/Core/ThreadPool.cpp @@ -92,6 +92,14 @@ void ThreadPool::addJob(Job newJob) condition_.notify_one(); } +// critical section //////////////////////////////////////////////////////////// +void ThreadPool::critical(Job fn) +{ + unique_lock lock(mutex_); + + fn(); +} + // wait for completion ///////////////////////////////////////////////////////// void ThreadPool::terminate(void) { diff --git a/lib/Core/ThreadPool.hpp b/lib/Core/ThreadPool.hpp index 5a504e5..6594a8e 100644 --- a/lib/Core/ThreadPool.hpp +++ b/lib/Core/ThreadPool.hpp @@ -37,6 +37,8 @@ public: 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: From b92fb84e9dbbee5a34748d74d772b52d0bc9c921 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 29 Nov 2021 00:03:19 +0000 Subject: [PATCH 10/23] adding pthread link --- configure.ac | 1 + 1 file changed, 1 insertion(+) diff --git a/configure.ac b/configure.ac index d1f2fab..27538e5 100644 --- a/configure.ac +++ b/configure.ac @@ -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])]) From c796187d1ef60eb11bb1621d3f18f5ade3461c91 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 14 Dec 2021 13:13:36 +0000 Subject: [PATCH 11/23] Absolute value plot (useful for log scale) --- lib/Core/Plot.cpp | 57 ++++++++++++++++++++++++++++++++++++++--------- lib/Core/Plot.hpp | 12 +++++----- 2 files changed, 54 insertions(+), 15 deletions(-) diff --git a/lib/Core/Plot.cpp b/lib/Core/Plot.cpp index d7b44fa..df21eba 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); @@ -217,7 +246,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 +260,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..23835d1 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; }; @@ -130,7 +131,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; }; From 57c600479711bdfb8715c7b07f230acc6b307950 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 20 Dec 2021 01:25:13 +0100 Subject: [PATCH 12/23] significant optimisation of covariance routines + checks --- examples/Makefile.am | 7 ++- examples/exVarBenchmark.cpp | 119 +++++++++++++++++++++++++++++++++++ lib/Statistics/MatSample.hpp | 76 ++++++++++++++++++++++ lib/Statistics/StatArray.hpp | 119 +++++++++++++++++++++++++---------- 4 files changed, 288 insertions(+), 33 deletions(-) create mode 100644 examples/exVarBenchmark.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index 123bafb..8a51549 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -20,7 +20,8 @@ noinst_PROGRAMS = \ exPValue \ exRand \ exRootFinder \ - exThreadPool + exThreadPool \ + exVarBenchmark exCompiledDoubleFunction_SOURCES = exCompiledDoubleFunction.cpp exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS) @@ -78,4 +79,8 @@ exThreadPool_SOURCES = exThreadPool.cpp exThreadPool_CXXFLAGS = $(COM_CXXFLAGS) exThreadPool_LDFLAGS = -L../lib/.libs -lLatAnalyze +exVarBenchmark_SOURCES = exVarBenchmark.cpp +exVarBenchmark_CXXFLAGS = $(COM_CXXFLAGS) +exVarBenchmark_LDFLAGS = -L../lib/.libs -lLatAnalyze + ACLOCAL_AMFLAGS = -I .buildutils/m4 diff --git a/examples/exVarBenchmark.cpp b/examples/exVarBenchmark.cpp new file mode 100644 index 0000000..4c0ea9c --- /dev/null +++ b/examples/exVarBenchmark.cpp @@ -0,0 +1,119 @@ +#include +#include +#include +#include +#include + +using namespace std; +using namespace Latan; + +constexpr Index size = 1000; +constexpr Index nSample = 2000; + +int main(void) +{ + random_device rd; + DMat var(size, size); + DVec mean(size); + DMatSample sample(nSample, size, 1), sample2(nSample, size, 1); + + cout << "-- generating " << nSample << " Gaussian random vectors..." << endl; + var = DMat::Random(size, size); + var *= var.adjoint(); + mean = DVec::Random(size); + RandomNormal mgauss(mean, var, rd()); + sample[central] = mgauss(); + FOR_STAT_ARRAY(sample, s) + { + sample[s] = mgauss(); + } + sample2[central] = mgauss(); + FOR_STAT_ARRAY(sample, s) + { + sample2[s] = mgauss(); + } + + cout << "-- check new routines" << endl; + + DMat v, vo; + + cout << "var" << endl; + auto start = chrono::high_resolution_clock::now(); + vo = sample.varianceOld(); + auto end = chrono::high_resolution_clock::now(); + chrono::duration diff = end - start; + cout << "time " << diff.count() << endl; + start = chrono::high_resolution_clock::now(); + v = sample.variance(); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + cout << "diff " << (v - vo).norm() << endl; + + cout << "cov" << endl; + start = chrono::high_resolution_clock::now(); + vo = sample.covarianceOld(sample2); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + start = chrono::high_resolution_clock::now(); + v = sample.covariance(sample2); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + cout << "diff " << (v - vo).norm() << endl; + + cout << "mean" << endl; + start = chrono::high_resolution_clock::now(); + vo = sample.meanOld(3, 5); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + start = chrono::high_resolution_clock::now(); + v = sample.mean(3, 5); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + cout << "diff " << (v - vo).norm() << endl; + + cout << "varmat" << endl; + start = chrono::high_resolution_clock::now(); + vo = sample.varianceMatrixOld(); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + start = chrono::high_resolution_clock::now(); + v = sample.varianceMatrix(); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + cout << "diff " << (v - vo).norm() << endl; + + cout << "covarmat" << endl; + start = chrono::high_resolution_clock::now(); + vo = sample.covarianceMatrixOld(sample2); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + start = chrono::high_resolution_clock::now(); + v = sample.covarianceMatrix(sample2); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + cout << "diff " << (v - vo).norm() << endl; + + cout << "corrmat" << endl; + start = chrono::high_resolution_clock::now(); + vo = sample.correlationMatrixOld(); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + start = chrono::high_resolution_clock::now(); + v = sample.correlationMatrix(); + end = chrono::high_resolution_clock::now(); + diff = end - start; + cout << "time " << diff.count() << endl; + cout << "diff " << (v - vo).norm() << endl; + + return EXIT_SUCCESS; +} 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..15179e1 100644 --- a/lib/Statistics/StatArray.hpp +++ b/lib/Statistics/StatArray.hpp @@ -51,14 +51,17 @@ 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, + T covarianceOld(const StatArray &array) const; + T covariance(const StatArray &array) const; + T covarianceMatrixOld(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 varianceOld(void) const; + T variance(void) const; + T varianceMatrixOld(const Index pos = 0, const Index n = -1) const; + T correlationMatrixOld(const Index pos = 0, const Index n = -1) const; // IO type virtual IoType getType(void) const; public: @@ -66,7 +69,7 @@ public: }; // reduction operations -namespace ReducOp +namespace StatOp { // general templates template @@ -149,42 +152,82 @@ void StatArray::bin(Index binSize) } template -T StatArray::mean(const Index pos, const Index n) const +T StatArray::meanOld(const Index pos, const Index n) const { T result = T(); const Index m = (n >= 0) ? n : size(); if (m) { - result = this->segment(pos+os, m).redux(&ReducOp::sum); + result = this->segment(pos+os, m).redux(&StatOp::sum); } return result/static_cast(m); } template -T StatArray::covariance(const StatArray &array, const Index pos, - const Index n) const +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 +{ + const Index m = (n >= 0) ? n : size(); + + return sum(pos, n)/static_cast(m); +} + +template +T StatArray::covarianceOld(const StatArray &array) const { T s1, s2, prs, res = T(); - const Index m = (n >= 0) ? n : size(); + const Index m = size(); if (m) { - auto arraySeg = array.segment(pos+os, m); - auto thisSeg = this->segment(pos+os, m); + auto arraySeg = array.segment(os, m); + auto thisSeg = this->segment(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); + s1 = thisSeg.redux(&StatOp::sum); + s2 = arraySeg.redux(&StatOp::sum); + prs = thisSeg.binaryExpr(arraySeg, &StatOp::prod) + .redux(&StatOp::sum); + res = prs - StatOp::prod(s1, s2)/static_cast(m); } return res/static_cast(m - 1); } template -T StatArray::covarianceMatrix(const StatArray &array, +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) + { + res += StatOp::prod(array[i], (*this)[i]); + } + res -= StatOp::prod(s1, s2)/static_cast(size()); + res /= static_cast(size() - 1); + + return res; +} + +template +T StatArray::covarianceMatrixOld(const StatArray &array, const Index pos, const Index n) const { T s1, s2, prs, res = T(); @@ -195,32 +238,38 @@ T StatArray::covarianceMatrix(const StatArray &array, 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); + s1 = thisSeg.redux(&StatOp::sum); + s2 = arraySeg.redux(&StatOp::sum); + prs = thisSeg.binaryExpr(arraySeg, &StatOp::tensProd) + .redux(&StatOp::sum); + res = prs - StatOp::tensProd(s1, s2)/static_cast(m); } return res/static_cast(m - 1); } template -T StatArray::variance(const Index pos, const Index n) const +T StatArray::variance(void) const { - return covariance(*this, pos, n); + return covariance(*this); } template -T StatArray::varianceMatrix(const Index pos, const Index n) const +T StatArray::varianceOld(void) const { - return covarianceMatrix(*this, pos, n); + return covarianceOld(*this); } template -T StatArray::correlationMatrix(const Index pos, const Index n) const +T StatArray::varianceMatrixOld(const Index pos, const Index n) const { - T res = varianceMatrix(pos, n); + return covarianceMatrixOld(*this, pos, n); +} + +template +T StatArray::correlationMatrixOld(const Index pos, const Index n) const +{ + T res = varianceMatrixOld(pos, n); T invDiag(res.rows(), 1); invDiag = res.diagonal(); @@ -231,8 +280,14 @@ T StatArray::correlationMatrix(const Index pos, const Index n) const } // reduction operations //////////////////////////////////////////////////////// -namespace ReducOp +namespace StatOp { + template + inline void zero(T &a) + { + a = 0.; + } + template inline T sum(const T &a, const T &b) { From 24a7b9c2033ecb85ee3e7064b3997a7256eb2fd9 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 20 Dec 2021 01:29:29 +0100 Subject: [PATCH 13/23] remove new covariance routine regression code --- examples/Makefile.am | 7 +-- examples/exVarBenchmark.cpp | 119 ----------------------------------- lib/Statistics/StatArray.hpp | 119 ----------------------------------- 3 files changed, 1 insertion(+), 244 deletions(-) delete mode 100644 examples/exVarBenchmark.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index 8a51549..123bafb 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -20,8 +20,7 @@ noinst_PROGRAMS = \ exPValue \ exRand \ exRootFinder \ - exThreadPool \ - exVarBenchmark + exThreadPool exCompiledDoubleFunction_SOURCES = exCompiledDoubleFunction.cpp exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS) @@ -79,8 +78,4 @@ exThreadPool_SOURCES = exThreadPool.cpp exThreadPool_CXXFLAGS = $(COM_CXXFLAGS) exThreadPool_LDFLAGS = -L../lib/.libs -lLatAnalyze -exVarBenchmark_SOURCES = exVarBenchmark.cpp -exVarBenchmark_CXXFLAGS = $(COM_CXXFLAGS) -exVarBenchmark_LDFLAGS = -L../lib/.libs -lLatAnalyze - ACLOCAL_AMFLAGS = -I .buildutils/m4 diff --git a/examples/exVarBenchmark.cpp b/examples/exVarBenchmark.cpp deleted file mode 100644 index 4c0ea9c..0000000 --- a/examples/exVarBenchmark.cpp +++ /dev/null @@ -1,119 +0,0 @@ -#include -#include -#include -#include -#include - -using namespace std; -using namespace Latan; - -constexpr Index size = 1000; -constexpr Index nSample = 2000; - -int main(void) -{ - random_device rd; - DMat var(size, size); - DVec mean(size); - DMatSample sample(nSample, size, 1), sample2(nSample, size, 1); - - cout << "-- generating " << nSample << " Gaussian random vectors..." << endl; - var = DMat::Random(size, size); - var *= var.adjoint(); - mean = DVec::Random(size); - RandomNormal mgauss(mean, var, rd()); - sample[central] = mgauss(); - FOR_STAT_ARRAY(sample, s) - { - sample[s] = mgauss(); - } - sample2[central] = mgauss(); - FOR_STAT_ARRAY(sample, s) - { - sample2[s] = mgauss(); - } - - cout << "-- check new routines" << endl; - - DMat v, vo; - - cout << "var" << endl; - auto start = chrono::high_resolution_clock::now(); - vo = sample.varianceOld(); - auto end = chrono::high_resolution_clock::now(); - chrono::duration diff = end - start; - cout << "time " << diff.count() << endl; - start = chrono::high_resolution_clock::now(); - v = sample.variance(); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - cout << "diff " << (v - vo).norm() << endl; - - cout << "cov" << endl; - start = chrono::high_resolution_clock::now(); - vo = sample.covarianceOld(sample2); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - start = chrono::high_resolution_clock::now(); - v = sample.covariance(sample2); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - cout << "diff " << (v - vo).norm() << endl; - - cout << "mean" << endl; - start = chrono::high_resolution_clock::now(); - vo = sample.meanOld(3, 5); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - start = chrono::high_resolution_clock::now(); - v = sample.mean(3, 5); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - cout << "diff " << (v - vo).norm() << endl; - - cout << "varmat" << endl; - start = chrono::high_resolution_clock::now(); - vo = sample.varianceMatrixOld(); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - start = chrono::high_resolution_clock::now(); - v = sample.varianceMatrix(); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - cout << "diff " << (v - vo).norm() << endl; - - cout << "covarmat" << endl; - start = chrono::high_resolution_clock::now(); - vo = sample.covarianceMatrixOld(sample2); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - start = chrono::high_resolution_clock::now(); - v = sample.covarianceMatrix(sample2); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - cout << "diff " << (v - vo).norm() << endl; - - cout << "corrmat" << endl; - start = chrono::high_resolution_clock::now(); - vo = sample.correlationMatrixOld(); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - start = chrono::high_resolution_clock::now(); - v = sample.correlationMatrix(); - end = chrono::high_resolution_clock::now(); - diff = end - start; - cout << "time " << diff.count() << endl; - cout << "diff " << (v - vo).norm() << endl; - - return EXIT_SUCCESS; -} diff --git a/lib/Statistics/StatArray.hpp b/lib/Statistics/StatArray.hpp index 15179e1..4e70dbb 100644 --- a/lib/Statistics/StatArray.hpp +++ b/lib/Statistics/StatArray.hpp @@ -54,14 +54,8 @@ public: 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 covarianceOld(const StatArray &array) const; T covariance(const StatArray &array) const; - T covarianceMatrixOld(const StatArray &array, const Index pos = 0, - const Index n = -1) const; - T varianceOld(void) const; T variance(void) const; - T varianceMatrixOld(const Index pos = 0, const Index n = -1) const; - T correlationMatrixOld(const Index pos = 0, const Index n = -1) const; // IO type virtual IoType getType(void) const; public: @@ -151,19 +145,6 @@ void StatArray::bin(Index binSize) } } -template -T StatArray::meanOld(const Index pos, const Index n) const -{ - T result = T(); - const Index m = (n >= 0) ? n : size(); - - if (m) - { - result = this->segment(pos+os, m).redux(&StatOp::sum); - } - return result/static_cast(m); -} - template T StatArray::sum(const Index pos, const Index n) const { @@ -187,27 +168,6 @@ T StatArray::mean(const Index pos, const Index n) const return sum(pos, n)/static_cast(m); } -template -T StatArray::covarianceOld(const StatArray &array) const -{ - T s1, s2, prs, res = T(); - const Index m = size(); - - if (m) - { - auto arraySeg = array.segment(os, m); - auto thisSeg = this->segment(os, m); - - s1 = thisSeg.redux(&StatOp::sum); - s2 = arraySeg.redux(&StatOp::sum); - prs = thisSeg.binaryExpr(arraySeg, &StatOp::prod) - .redux(&StatOp::sum); - res = prs - StatOp::prod(s1, s2)/static_cast(m); - } - - return res/static_cast(m - 1); -} - template T StatArray::covariance(const StatArray &array) const { @@ -226,105 +186,26 @@ T StatArray::covariance(const StatArray &array) const return res; } -template -T StatArray::covarianceMatrixOld(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(&StatOp::sum); - s2 = arraySeg.redux(&StatOp::sum); - prs = thisSeg.binaryExpr(arraySeg, &StatOp::tensProd) - .redux(&StatOp::sum); - res = prs - StatOp::tensProd(s1, s2)/static_cast(m); - } - - return res/static_cast(m - 1); -} - template T StatArray::variance(void) const { return covariance(*this); } -template -T StatArray::varianceOld(void) const -{ - return covarianceOld(*this); -} - -template -T StatArray::varianceMatrixOld(const Index pos, const Index n) const -{ - return covarianceMatrixOld(*this, pos, n); -} - -template -T StatArray::correlationMatrixOld(const Index pos, const Index n) const -{ - T res = varianceMatrixOld(pos, n); - T invDiag(res.rows(), 1); - - invDiag = res.diagonal(); - invDiag = invDiag.cwiseInverse().cwiseSqrt(); - res = (invDiag*invDiag.transpose()).cwiseProduct(res); - - return res; -} - // reduction operations //////////////////////////////////////////////////////// namespace StatOp { - template - inline void zero(T &a) - { - a = 0.; - } - - template - inline T sum(const T &a, const T &b) - { - return a + b; - } - 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 ///////////////////////////////////////////////////////////////////// From adf2c9cc69910085aacda39378badb46627eefdd Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 20 Dec 2021 01:30:03 +0100 Subject: [PATCH 14/23] normalised residuals routines --- lib/Statistics/XYSampleData.cpp | 23 +++++++++++++++++++++++ lib/Statistics/XYSampleData.hpp | 1 + lib/Statistics/XYStatData.cpp | 18 ++++++++++++++++++ lib/Statistics/XYStatData.hpp | 1 + 4 files changed, 43 insertions(+) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 6e2d2d6..b2fb4a4 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -346,6 +346,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 6cba855..a629ece 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -111,6 +111,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..2a88072 100644 --- a/lib/Statistics/XYStatData.cpp +++ b/lib/Statistics/XYStatData.cpp @@ -379,6 +379,24 @@ 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); + + for (auto &p: yData_[j]) + { + res.y(p.first, j) -= f(x(p.first)); + res.y(p.first, j) /= sqrt(getYYVar(j, j)(p.first, p.first)); + } + } + + return res; +} + XYStatData XYStatData::getPartialResiduals(const FitResult &fit, const DVec &ref, const Index i) { diff --git a/lib/Statistics/XYStatData.hpp b/lib/Statistics/XYStatData.hpp index e14b15c..b733dee 100644 --- a/lib/Statistics/XYStatData.hpp +++ b/lib/Statistics/XYStatData.hpp @@ -104,6 +104,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: From 60d91cbff535a6e16469c537ed8f58d75b94b46f Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Mon, 20 Dec 2021 01:30:26 +0100 Subject: [PATCH 15/23] plot data with points --- lib/Core/Plot.cpp | 18 ++++++++++++++++++ lib/Core/Plot.hpp | 9 +++++++++ 2 files changed, 27 insertions(+) diff --git a/lib/Core/Plot.cpp b/lib/Core/Plot.cpp index df21eba..42848ed 100644 --- a/lib/Core/Plot.cpp +++ b/lib/Core/Plot.cpp @@ -224,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) { diff --git a/lib/Core/Plot.hpp b/lib/Core/Plot.hpp index 23835d1..f84db3d 100644 --- a/lib/Core/Plot.hpp +++ b/lib/Core/Plot.hpp @@ -116,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: From bac8356de57ec12717e978cba046b93ebfee2328 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 21 Dec 2021 18:18:14 +0100 Subject: [PATCH 16/23] Ignore verbose minimiser for samples --- lib/Statistics/XYSampleData.cpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index b2fb4a4..99e4996 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -285,9 +285,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_); @@ -299,9 +300,14 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, { sampleResult = data_.fit(minimizer, initCopy, v); initCopy = sampleResult.segment(0, initCopy.size()); + if (verbCopy != Minimizer::Verbosity::Debug) + { + minimizer.back()->setVerbosity(Minimizer::Verbosity::Silent); + } } else { + sampleResult = data_.fit(*(minimizer.back()), initCopy, v); } result[s] = sampleResult; @@ -312,6 +318,7 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.model_[j][s] = sampleResult.getModel(j); } } + minimizer.back()->setVerbosity(verbCopy); result.nPar_ = sampleResult.getNPar(); result.nDof_ = sampleResult.nDof_; result.parName_ = sampleResult.parName_; From 8cd29c2bee61daa10ada741f80f8464509355abb Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Sun, 26 Dec 2021 22:21:03 +0100 Subject: [PATCH 17/23] XYStatData fixes --- lib/Statistics/XYStatData.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/lib/Statistics/XYStatData.cpp b/lib/Statistics/XYStatData.cpp index 2a88072..7407998 100644 --- a/lib/Statistics/XYStatData.cpp +++ b/lib/Statistics/XYStatData.cpp @@ -216,7 +216,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(); } @@ -385,12 +385,15 @@ XYStatData XYStatData::getNormalisedResiduals(const FitResult &fit) for (Index j = 0; j < getNYDim(); ++j) { - const DoubleFunction &f = fit.getModel(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) /= sqrt(getYYVar(j, j)(p.first, p.first)); + res.y(p.first, j) /= err(row); + row++; } } From e4cefae515d9bf1f99b6a4c6c1f0ca1081e245b7 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Tue, 28 Dec 2021 12:17:02 +0100 Subject: [PATCH 18/23] special macro to plot correlation matrices --- lib/Core/Plot.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/lib/Core/Plot.hpp b/lib/Core/Plot.hpp index f84db3d..1c2901c 100644 --- a/lib/Core/Plot.hpp +++ b/lib/Core/Plot.hpp @@ -193,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 * ******************************************************************************/ From 0de8091f3c47a194e449e3ad03f73596627ca320 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 16 Feb 2022 18:53:45 +0000 Subject: [PATCH 19/23] Eigen plugin const --- lib/Core/EigenPlugin.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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(); From 857a8e59c9aac045c1814b6055cb4d64081dce17 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 16 Feb 2022 18:54:16 +0000 Subject: [PATCH 20/23] corr to var and SVD dynamic range --- lib/Core/Math.cpp | 25 ++++++++++++++++++++++++- lib/Core/Math.hpp | 5 +++++ 2 files changed, 29 insertions(+), 1 deletion(-) 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; From ebc1bd4c2e9b59bfad2f69db501f22c3d847a587 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 16 Feb 2022 18:55:08 +0000 Subject: [PATCH 21/23] fit: stable variance inversion and SVD dynamic range --- lib/Statistics/XYSampleData.cpp | 30 ++++++++++++++++++++++---- lib/Statistics/XYSampleData.hpp | 6 +++++- lib/Statistics/XYStatData.cpp | 37 +++++++++++++++++++++++++++------ lib/Statistics/XYStatData.hpp | 9 +++++--- 4 files changed, 68 insertions(+), 14 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index 99e4996..b2aa4bd 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::getSvdRangeDb(void) const +{ + return svdRangeDb_; +} + 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", getSvdRangeDb()); + 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) { @@ -319,9 +340,10 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, } } minimizer.back()->setVerbosity(verbCopy); - result.nPar_ = sampleResult.getNPar(); - result.nDof_ = sampleResult.nDof_; - result.parName_ = sampleResult.parName_; + result.nPar_ = sampleResult.getNPar(); + result.nDof_ = sampleResult.nDof_; + result.parName_ = sampleResult.parName_; + result.svdRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); return result; } diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index a629ece..b5c7a45 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 getSvdRangeDb(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 svdRangeDb_{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 diff --git a/lib/Statistics/XYStatData.cpp b/lib/Statistics/XYStatData.cpp index 7407998..2052be8 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::getSvdRangeDb(void) const +{ + return svdRangeDb_; +} + 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", getSvdRangeDb()); + 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; } } @@ -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.svdRangeDb_ = 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) { @@ -551,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 b733dee..45cd2be 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 getSvdRangeDb(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.}, svdRangeDb_{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); @@ -131,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}; From 35f67332923688ecaa28634a858d2592d68ce43a Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 16 Feb 2022 18:55:24 +0000 Subject: [PATCH 22/23] sample-plot-corr better palette --- utils/sample-plot-corr.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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()) { From 47d0b3f0404b9e37d75a358dde7ada4162d9b8c7 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Wed, 16 Feb 2022 19:03:19 +0000 Subject: [PATCH 23/23] correlation dynamic range renaming --- lib/Statistics/XYSampleData.cpp | 8 ++++---- lib/Statistics/XYSampleData.hpp | 4 ++-- lib/Statistics/XYStatData.cpp | 8 ++++---- lib/Statistics/XYStatData.hpp | 4 ++-- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/lib/Statistics/XYSampleData.cpp b/lib/Statistics/XYSampleData.cpp index b2aa4bd..0676564 100644 --- a/lib/Statistics/XYSampleData.cpp +++ b/lib/Statistics/XYSampleData.cpp @@ -62,9 +62,9 @@ double SampleFitResult::getPValue(const Index s) const return Math::chi2PValue(getChi2(s), getNDof()); } -double SampleFitResult::getSvdRangeDb(void) const +double SampleFitResult::getCorrRangeDb(void) const { - return svdRangeDb_; + return corrRangeDb_; } double SampleFitResult::getCcdf(const Index s) const @@ -112,7 +112,7 @@ 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", getSvdRangeDb()); + sprintf(buf, "correlation dynamic range= %.1f dB", getCorrRangeDb()); out << buf << endl; for (Index p = 0; p < pMax; ++p) { @@ -343,7 +343,7 @@ SampleFitResult XYSampleData::fit(std::vector &minimizer, result.nPar_ = sampleResult.getNPar(); result.nDof_ = sampleResult.nDof_; result.parName_ = sampleResult.parName_; - result.svdRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); + result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); return result; } diff --git a/lib/Statistics/XYSampleData.hpp b/lib/Statistics/XYSampleData.hpp index b5c7a45..bb82748 100644 --- a/lib/Statistics/XYSampleData.hpp +++ b/lib/Statistics/XYSampleData.hpp @@ -49,7 +49,7 @@ public: double getNDof(void) const; Index getNPar(void) const; double getPValue(const Index s = central) const; - double getSvdRangeDb(void) const; + double getCorrRangeDb(void) const; double getCcdf(const Index s = central) const; const DoubleFunction & getModel(const Index s = central, const Index j = 0) const; @@ -61,7 +61,7 @@ public: std::ostream &out = std::cout) const; private: DSample chi2_; - double svdRangeDb_{0.}; + double corrRangeDb_{0.}; Index nDof_{0}, nPar_{0}; std::vector model_; std::vector parName_; diff --git a/lib/Statistics/XYStatData.cpp b/lib/Statistics/XYStatData.cpp index 2052be8..690a3a6 100644 --- a/lib/Statistics/XYStatData.cpp +++ b/lib/Statistics/XYStatData.cpp @@ -60,9 +60,9 @@ double FitResult::getCcdf(void) const return Math::chi2Ccdf(getChi2(), getNDof());; } -double FitResult::getSvdRangeDb(void) const +double FitResult::getCorrRangeDb(void) const { - return svdRangeDb_; + return corrRangeDb_; } const DoubleFunction & FitResult::getModel(const Index j) const @@ -80,7 +80,7 @@ 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", getSvdRangeDb()); + sprintf(buf, "correlation dynamic range= %.1f dB", getCorrRangeDb()); out << buf << endl; for (Index p = 0; p < pMax; ++p) { @@ -358,7 +358,7 @@ FitResult XYStatData::fit(vector &minimizer, const DVec &init, result = (*m)(chi2); totalInit = result; } - result.svdRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); + result.corrRangeDb_ = Math::svdDynamicRangeDb(getFitCorrMat()); result.chi2_ = chi2(result); result.nPar_ = nPar; result.nDof_ = layout.totalYSize - nPar; diff --git a/lib/Statistics/XYStatData.hpp b/lib/Statistics/XYStatData.hpp index 45cd2be..49cbd28 100644 --- a/lib/Statistics/XYStatData.hpp +++ b/lib/Statistics/XYStatData.hpp @@ -48,13 +48,13 @@ public: Index getNPar(void) const; double getPValue(void) const; double getCcdf(void) const; - double getSvdRangeDb(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.}, svdRangeDb_{0.}; + double chi2_{0.}, corrRangeDb_{0.}; Index nDof_{0}, nPar_{0}; std::vector model_; std::vector parName_;