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

Merge branch 'release/3.1'

Browse Source
This commit is contained in:
Antonin Portelli 2016-04-13 12:47:07 +01:00
commit 4b2d070e02
102 changed files with 4009 additions and 2758 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -27,3 +27,6 @@ lib/*Parser.hpp
# Eigen headers
lib/Eigen/*
lib/eigen_files.mk
# CI builds
ci-scripts/local/*

84
.travis.yml Normal file
View File

@ -0,0 +1,84 @@
language: cpp
notifications:
email: false
slack: ukqcd:mQLXCtz8D2cg89xT8j1a4wku
cache:
directories:
- ci-scripts/local
matrix:
include:
- os: osx
osx_image: xcode7.2
compiler: clang
- os: osx
osx_image: xcode7.2
compiler: gcc
env: VERSION=-5
- compiler: gcc
addons:
apt:
sources:
- ubuntu-toolchain-r-test
packages:
- g++-4.9
- libgsl0-dev
- flex
- bison
env: VERSION=-4.9
- compiler: gcc
addons:
apt:
sources:
- ubuntu-toolchain-r-test
packages:
- g++-5
- libgsl0-dev
- flex
- bison
env: VERSION=-5
- compiler: clang
addons:
apt:
sources:
- ubuntu-toolchain-r-test
- llvm-toolchain-precise-3.6
packages:
- clang-3.6
- libgsl0-dev
- flex
- bison
env: VERSION=-3.6
- compiler: clang
addons:
apt:
sources:
- ubuntu-toolchain-r-test
- llvm-toolchain-precise-3.7
packages:
- clang-3.7
- libgsl0-dev
- flex
- bison
env: VERSION=-3.7
before_install:
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew update; fi
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install gsl; fi
- if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install bison flex; export PATH="/usr/local/opt/flex/bin:/usr/local/opt/bison/bin:$PATH"; fi
- if [[ "$TRAVIS_OS_NAME" == "osx" ]] && [[ "$CC" == "gcc" ]]; then brew install gcc5; fi
install:
- export LATDIR=`pwd`
- export CC=$CC$VERSION
- export CXX=$CXX$VERSION
- cd ci-scripts
- ./install-deps.sh `pwd`/local
- cd ..
- if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then export LD_LIBRARY_PATH=${LATDIR}/ci-scripts/local/lib:$LD_LIBRARY_PATH; fi
script:
- cd ci-scripts
- ./install-latan.sh `pwd`/local $TRAVIS_OS_NAME

74
Readme.md Normal file
View File

@ -0,0 +1,74 @@
# LatAnalyze
Contributors: Antonin Portelli, Matt Spraggs
License: GNU General Public License v3
<table>
<tr>
<td>Last stable release</td>
<td><a href="https://travis-ci.org/aportelli/LatAnalyze">
<img src="https://travis-ci.org/aportelli/LatAnalyze.svg?branch=master"></a>
</td>
</tr>
<tr>
<td>Development branch</td>
<td><a href="https://travis-ci.org/aportelli/LatAnalyze">
<img src="https://travis-ci.org/aportelli/LatAnalyze.svg?branch=develop"></a>
</td>
</tr>
</table>
## Description
LatAnalyze is a C++11 library for statistical data analysis based on bootstrap
resampling. It has been written with lattice QCD data analysis in mind (hence
the name), but its features are not lattice specific and can be used more general statistical context.
Sadly a proper documentation was never written, but some comprehensive examples covering most features can be found in the `examples` directory.
The main features are the following:
* Array and matrix types with fast arithmetic operations based on [Eigen](http://eigen.tuxfamily.org).
* High-level types for bootstrap sample manipulation (including various statistical estimators and histogramming).
* Mathematical expression parser for runtime defined functions.
* Data I/O in ASCII and HDF5 (optional).
* High-level wrappers to minimisation routines from the [GSL](http://www.gnu.org/software/gsl/), [Minuit](http://seal.web.cern.ch/seal/snapshot/work-packages/mathlibs/minuit/) (optional) and [NLopt](http://ab-initio.mit.edu/wiki/index.php/NLopt).
* Non-linear regression with error on independent variables (through total least squares).
* High-level wrappers to numerical integrator and non-linear solver from the [GSL](http://www.gnu.org/software/gsl/).
* High-level functional types for function of model. General functions can be defined from C pointers, C++ objects, strings of mathematical expressions or tabulated data.
* High-level plotting functions based on [gnuplot](http://www.gnuplot.info), with the possibility of generating and saving plot scripts.
## Installation
The head of the `master` branch always points to the latest stable release. The `develop` branch is the main unstable branch of LatAnalyze.
LatAnalyze is written in C++11 and requires a rather recent C++ compiler to be built. It has been successfully built on various Linux and OS X platforms using clang (from 3.7), GCC (from 4.9) and the Intel C++ compiler (2016).
The strict dependencies are the [GSL](http://www.gnu.org/software/gsl/) and [LatCore](https://github.com/aportelli/LatCore).
Additionally, autoconf, automake (from 1.11), libtool, bison (from 3.0) and flex are necessary to build the library. Unless you use a very exotic system, these tools are standard on any Unix platform and should be already present or easy to install through a package manager.
Optional dependencies are [HDF5](https://www.hdfgroup.org/HDF5/) (built with C++ support), [Minuit](http://seal.web.cern.ch/seal/snapshot/work-packages/mathlibs/minuit/) and [NLopt](http://ab-initio.mit.edu/wiki/index.php/NLopt).
For a quick installation with all possible extensions execute:
```
./install-latan.sh <prefix> {osx|linux}
```
in the `ci-scripts` directory where `<prefix>` is where you want LatAnalyze (and its dependencies) to be installed and `{osx|linux}` should be the name of your OS. This script will automatically download, build and install LatCore, HDF5, Minuit and NLopt. It is assumed that the GSL can be found in a standard location (_e.g._ it has been installed through a package manager for Linux and is present in `/usr/local` for OS X).
For a more customised installation, one first needs to generate the build system by running `./bootstrap.sh` in the root directory. Then the library can be built and installed through the usual GNU mantra `./configure <options> && make && make install`. Use `./configure --help` to obtain a list of possible options for `./configure`. Because Eigen expressions rely a lot on inlining and compiler optimisations it is strongly recommended to set the `CXXFLAGS` variable to `-O3 -march=native -mtune=native`.
## History
#### v3.1
Additions:
* Wrappers to NLopt and GSL minimisers.
* Command-line tool to plot the correlation heatmap of a boostrap sample file.
* I/O functions for `DSample` type.
Changes:
* Internal random generator removed (obsolete because of C++11 pseudo-random generators).
* Fit interface and `XY*Data` classes rewritten from scratch for improved flexibility and performance.
Fixes:
* **Loads** of portability and compatibility fixes and [CI with Travis](https://travis-ci.org/aportelli/LatAnalyze).
#### 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.**

16
ci-scripts/install-deps.sh Executable file
View File

@ -0,0 +1,16 @@
#!/usr/bin/env bash
if (( $# != 1 )); then
echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
exit 1
fi
PREFIX=$1
set -ex
mkdir -p local/build
for d in nlopt minuit hdf5; do
if [ ! -e local/.built.${d} ]; then
./install-${d}.sh ${PREFIX}
fi
done
./install-latcore.sh ${PREFIX}

23
ci-scripts/install-hdf5.sh Executable file
View File

@ -0,0 +1,23 @@
#!/usr/bin/env bash
NAME='hdf5-1.8.16'
if (( $# != 1 )); then
echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
exit 1
fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd local/build
wget http://www.hdfgroup.org/ftp/HDF5/current/src/${NAME}.tar.gz
tar -xzvf ${NAME}.tar.gz
mkdir ${NAME}/build
cd ${NAME}/build
../configure --prefix=${PREFIX} --enable-cxx
make -j4
make install
cd ${INITDIR}/local
touch .built.hdf5
cd ${INITDIR}

19
ci-scripts/install-latan.sh Executable file
View File

@ -0,0 +1,19 @@
#!/usr/bin/env bash
if (( $# != 2 )); then
echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
exit 1
fi
PREFIX=$1
OS=$2
set -ex
./install-deps.sh ${PREFIX}
cd ..
./bootstrap.sh
mkdir -p build
cd build
if [[ "$OS" == "osx" ]]; then EXTRA_FLAGS='--with-gsl=/usr/local'; fi
../configure --prefix=$PREFIX --with-minuit=$PREFIX --with-nlopt=$PREFIX --with-latcore=$PREFIX --with-hdf5=$PREFIX $EXTRA_FLAGS CXXFLAGS='-O3 -march=native -mtune=native'
make -j4
make install

25
ci-scripts/install-latcore.sh Executable file
View File

@ -0,0 +1,25 @@
#!/usr/bin/env bash
if (( $# != 1 )); then
echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
exit 1
fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd local/build
if [ -d LatCore ]; then
cd LatCore
git pull origin master
else
git clone https://github.com/aportelli/LatCore.git
mkdir LatCore/build
cd LatCore
./bootstrap.sh
fi
cd build
../configure --prefix=${PREFIX}
make -j4
make install
cd ${INITDIR}

23
ci-scripts/install-minuit.sh Executable file
View File

@ -0,0 +1,23 @@
#!/usr/bin/env bash
NAME='Minuit2-5.34.14'
if (( $# != 1 )); then
echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
exit 1
fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd 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
make install
cd ${INITDIR}/local
touch .built.minuit
cd ${INITDIR}

23
ci-scripts/install-nlopt.sh Executable file
View File

@ -0,0 +1,23 @@
#!/usr/bin/env bash
NAME='nlopt-2.4.2'
if (( $# != 1 )); then
echo "usage: `basename $0` <prefix> {osx|linux}" 1>&2
exit 1
fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd local/build
wget http://ab-initio.mit.edu/nlopt/${NAME}.tar.gz
tar -xzvf ${NAME}.tar.gz
mkdir -p ${NAME}/build
cd ${NAME}/build
../configure --prefix=${PREFIX} --with-cxx --without-guile --without-python --without-octave --without-matlab --with-pic
make -j4
make install
cd ${INITDIR}/local
touch .built.nlopt
cd ${INITDIR}

94
configure.ac
View File

@ -2,16 +2,17 @@
# Initialization
AC_PREREQ([2.63])
AC_INIT([LatAnalyze],[3.0],[antonin.portelli@me.com],[LatAnalyze])
AC_INIT([LatAnalyze],[3.1-rc],[antonin.portelli@me.com],[LatAnalyze])
AC_CONFIG_AUX_DIR([.buildutils])
AC_CONFIG_SRCDIR([lib/Global.cpp])
AC_CONFIG_SRCDIR([utils/sample_read.cpp])
AC_CONFIG_SRCDIR([examples/exMathInterpreter.cpp])
AC_CONFIG_MACRO_DIR([.buildutils/m4])
AM_INIT_AUTOMAKE([-Wall -Werror foreign])
AM_INIT_AUTOMAKE([1.11 -Wall -Werror foreign])
AM_SILENT_RULES([yes])
AC_CONFIG_HEADERS([config.h])
AM_CONDITIONAL([HAVE_AM_MINOR_LE_11],
[test `automake --version | grep automake | awk -F '.' '{print $2}'` -le 11])
# Checks for programs
AC_PROG_CXX
AC_PROG_AWK
@ -27,37 +28,34 @@ m4_ifdef([AM_PROG_AR],[AM_PROG_AR])
LT_INIT
# Options
AC_ARG_ENABLE([SSE],
[AS_HELP_STRING([--enable-SSE],
[compiles SSE version of ranlux random generator])],
[AC_DEFINE([HAVE_SSE],
[1],
[Define to 1 if your CPU support SSE instructions.])])
AC_ARG_WITH([gsl],
[AS_HELP_STRING([--with-gsl=prefix],
[try this for a non-standard install prefix of the GSL library])],
[AM_CFLAGS="$AM_CFLAGS -I$with_gsl/include"]
[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_gsl/include"]
[AM_LDFLAGS="$AM_LDFLAGS -L$with_gsl/lib"])
AC_ARG_WITH([Minuit2],
[AS_HELP_STRING([--with-Minuit2=prefix],
AC_ARG_WITH([minuit],
[AS_HELP_STRING([--with-minuit=prefix],
[try this for a non-standard install prefix of the Minuit2 library])],
[AM_CFLAGS="$AM_CFLAGS -I$with_Minuit2/include"]
[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_Minuit2/include"]
[AM_LDFLAGS="$AM_LDFLAGS -L$with_Minuit2/lib"])
AC_ARG_WITH([LatCore],
[AS_HELP_STRING([--with-LatCore=prefix],
[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_minuit/include"]
[AM_LDFLAGS="$AM_LDFLAGS -L$with_minuit/lib"])
AC_ARG_WITH([nlopt],
[AS_HELP_STRING([--with-nlopt=prefix],
[try this for a non-standard install prefix of the NLopt library])],
[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_nlopt/include"]
[AM_LDFLAGS="$AM_LDFLAGS -L$with_nlopt/lib"])
AC_ARG_WITH([hdf5],
[AS_HELP_STRING([--with-hdf5=prefix],
[try this for a non-standard install prefix of the HDF5 library])],
[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_hdf5/include"]
[AM_LDFLAGS="$AM_LDFLAGS -L$with_hdf5/lib"])
AC_ARG_WITH([latcore],
[AS_HELP_STRING([--with-latcore=prefix],
[use this option for a non-standard install prefix of the LatCore library])],
[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_LatCore/include"]
[AM_LDFLAGS="$AM_LDFLAGS -L$with_LatCore/lib"])
AX_LIB_HDF5()
if test x$with_hdf5 = xno; then
AC_MSG_ERROR([HDF5 library not found])
fi
CFLAGS="$AM_CFLAGS $HDF5_CFLAGS $CFLAGS"
CXXFLAGS="$AM_CXXFLAGS $HDF5_CPPFLAGS $CXXFLAGS"
LDFLAGS="$AM_LDFLAGS $HDF5_LDFLAGS $LDFLAGS"
LIBS="$LIBS $HDF5_LIBS -lhdf5_cpp"
[AM_CXXFLAGS="$AM_CXXFLAGS -I$with_latcore/include"]
[AM_LDFLAGS="$AM_LDFLAGS -L$with_latcore/lib"])
CFLAGS="$AM_CFLAGS $CFLAGS"
CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
LDFLAGS="$AM_LDFLAGS $LDFLAGS"
# Get compilers informations
AX_COMPILER_VENDOR
@ -84,8 +82,20 @@ 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])])
AC_CHECK_LIB([gsl],[gsl_blas_dgemm],[],[AC_MSG_ERROR([GSL library not found])])
AC_CHECK_LIB([LatCore],[_ZN7LatCore12testFunctionEv],[],
[AC_MSG_ERROR([LatCore library not found])])
AC_CHECK_LIB([nlopt_cxx],[nlopt_create],
[AC_DEFINE([HAVE_NLOPT],
[1],
[Define to 1 if you have the `NLopt' library (-lnlopt_cxx).])]
[have_nlopt=true]
[LIBS="$LIBS -lnlopt_cxx"],[])
AM_CONDITIONAL([HAVE_NLOPT], [test x$have_nlopt = xtrue])
AC_CHECK_LIB([hdf5_cpp],[H5Fopen],
[AC_DEFINE([HAVE_HDF5],
[1],
[Define to 1 if you have the `HDF5' library (-lhdf5_cpp).])]
[have_hdf5=true]
[LIBS="$LIBS -lhdf5_cpp -lhdf5"],[],[-lhdf5])
AM_CONDITIONAL([HAVE_HDF5], [test x$have_hdf5 = xtrue])
SAVED_LDFLAGS=$LDFLAGS
LDFLAGS="$LDFLAGS -lMinuit2"
AC_MSG_CHECKING([for ROOT::Minuit2::BasicMinimumError in -lMinuit2]);
@ -102,6 +112,17 @@ AC_LINK_IFELSE(
[AC_MSG_RESULT([no])])
AM_CONDITIONAL([HAVE_MINUIT], [test x$have_minuit = xtrue])
LDFLAGS=$SAVED_LDFLAGS
SAVED_LDFLAGS=$LDFLAGS
LDFLAGS="$LDFLAGS -lLatCore"
AC_MSG_CHECKING([for LatCore::XmlReader in -lLatCore]);
AC_LINK_IFELSE(
[AC_LANG_PROGRAM([#include <LatCore/XmlReader.hpp>],
[LatCore::XmlReader dummy()])],
[LIBS="$LIBS -lLatCore"]
[AC_MSG_RESULT([yes])],
[AC_MSG_RESULT([no])]
[AC_MSG_ERROR([LatCore library not found])])
LDFLAGS=$SAVED_LDFLAGS
# Checks for header files.
AC_HEADER_STDC
@ -112,3 +133,18 @@ AC_SUBST([AM_LDFLAGS])
AC_CONFIG_FILES([Makefile lib/Makefile utils/Makefile examples/Makefile])
AC_OUTPUT
echo "*********************************************"
echo "* ${PACKAGE_NAME} v${VERSION}" build options
echo "*********************************************"
echo "* C++ compiler : $CXX"
echo "* HDF5 support : `test x$HAVE_HDF5_TRUE = 'x' && echo yes \
|| echo no`"
echo "* Minuit minimizers: `test x$HAVE_MINUIT_TRUE = 'x' && echo yes \
|| echo no`"
echo "* NLopt minimizers : `test x$HAVE_NLOPT_TRUE = 'x' && echo yes \
|| echo no`"
echo "*********************************************"

91
examples/Makefile.am
View File

@ -1,80 +1,71 @@
if CC_GNU
COM_CFLAGS = -Wall -W -pedantic
else
if CC_INTEL
COM_CFLAGS = -Wall
endif
endif
if CXX_GNU
COM_CXXFLAGS = -Wall -W -pedantic
COM_CXXFLAGS = -Wall -W -pedantic -Wno-deprecated-declarations
else
if CXX_INTEL
COM_CXXFLAGS = -Wall
COM_CXXFLAGS = -wd1682 -Wall
endif
endif
noinst_PROGRAMS = \
exCompiledDoubleFunction\
exDerivative \
exFit \
exFitSample \
exIntegrator \
exInterp \
exMat \
exMathInterpreter \
exMin \
exPlot \
exRand \
exRootFinder
if HAVE_MINUIT
noinst_PROGRAMS += exFit exMin
endif
exCompiledDoubleFunction_SOURCES = exCompiledDoubleFunction.cpp
exCompiledDoubleFunction_CXXFLAGS = $(COM_CXXFLAGS)
exCompiledDoubleFunction_LDFLAGS = -L../lib/.libs -lLatAnalyze
exCompiledDoubleFunction_SOURCES = exCompiledDoubleFunction.cpp
exCompiledDoubleFunction_CFLAGS = -g -O2
exCompiledDoubleFunction_LDFLAGS = -L../lib/.libs -lLatAnalyze
exDerivative_SOURCES = exDerivative.cpp
exDerivative_CXXFLAGS = $(COM_CXXFLAGS)
exDerivative_LDFLAGS = -L../lib/.libs -lLatAnalyze
exDerivative_SOURCES = exDerivative.cpp
exDerivative_CFLAGS = -g -O2
exDerivative_LDFLAGS = -L../lib/.libs -lLatAnalyze
exFit_SOURCES = exFit.cpp
exFit_CXXFLAGS = $(COM_CXXFLAGS)
exFit_LDFLAGS = -L../lib/.libs -lLatAnalyze
if HAVE_MINUIT
exFit_SOURCES = exFit.cpp
exFit_CFLAGS = -g -O2
exFit_LDFLAGS = -L../lib/.libs -lLatAnalyze
endif
exFitSample_SOURCES = exFitSample.cpp
exFitSample_CXXFLAGS = $(COM_CXXFLAGS)
exFitSample_LDFLAGS = -L../lib/.libs -lLatAnalyze
exInterp_SOURCES = exInterp.cpp
exInterp_CFLAGS = -g -O2
exInterp_LDFLAGS = -L../lib/.libs -lLatAnalyze
exInterp_SOURCES = exInterp.cpp
exInterp_CXXFLAGS = $(COM_CXXFLAGS)
exInterp_LDFLAGS = -L../lib/.libs -lLatAnalyze
exIntegrator_SOURCES = exIntegrator.cpp
exIntegrator_CFLAGS = -g -O2
exIntegrator_LDFLAGS = -L../lib/.libs -lLatAnalyze
exIntegrator_SOURCES = exIntegrator.cpp
exIntegrator_CXXFLAGS = $(COM_CXXFLAGS)
exIntegrator_LDFLAGS = -L../lib/.libs -lLatAnalyze
exMat_SOURCES = exMat.cpp
exMat_CFLAGS = -g -O2
exMat_LDFLAGS = -L../lib/.libs -lLatAnalyze
exMat_SOURCES = exMat.cpp
exMat_CXXFLAGS = $(COM_CXXFLAGS)
exMat_LDFLAGS = -L../lib/.libs -lLatAnalyze
if HAVE_MINUIT
exMin_SOURCES = exMin.cpp
exMin_CFLAGS = -g -O2
exMin_LDFLAGS = -L../lib/.libs -lLatAnalyze
endif
exMin_SOURCES = exMin.cpp
exMin_CXXFLAGS = $(COM_CXXFLAGS)
exMin_LDFLAGS = -L../lib/.libs -lLatAnalyze
exMathInterpreter_SOURCES = exMathInterpreter.cpp
exMathInterpreter_CFLAGS = -g -O2
exMathInterpreter_LDFLAGS = -L../lib/.libs -lLatAnalyze
exMathInterpreter_SOURCES = exMathInterpreter.cpp
exMathInterpreter_CXXFLAGS = $(COM_CXXFLAGS)
exMathInterpreter_LDFLAGS = -L../lib/.libs -lLatAnalyze
exPlot_SOURCES = exPlot.cpp
exPlot_CFLAGS = -g -O2
exPlot_LDFLAGS = -L../lib/.libs -lLatAnalyze
exPlot_SOURCES = exPlot.cpp
exPlot_CXXFLAGS = $(COM_CXXFLAGS)
exPlot_LDFLAGS = -L../lib/.libs -lLatAnalyze
exRand_SOURCES = exRand.cpp
exRand_CFLAGS = -g -O2
exRand_LDFLAGS = -L../lib/.libs -lLatAnalyze
exRand_SOURCES = exRand.cpp
exRand_CXXFLAGS = $(COM_CXXFLAGS)
exRand_LDFLAGS = -L../lib/.libs -lLatAnalyze
exRootFinder_SOURCES = exRootFinder.cpp
exRootFinder_CFLAGS = -g -O2
exRootFinder_LDFLAGS = -L../lib/.libs -lLatAnalyze
exRootFinder_SOURCES = exRootFinder.cpp
exRootFinder_CXXFLAGS = $(COM_CXXFLAGS)
exRootFinder_LDFLAGS = -L../lib/.libs -lLatAnalyze
ACLOCAL_AMFLAGS = -I .buildutils/m4

2
examples/exCompiledDoubleFunction.cpp
View File

@ -1,5 +1,3 @@
#include <iostream>
#include <iomanip>
#include <LatAnalyze/CompiledFunction.hpp>
using namespace std;

1
examples/exDerivative.cpp
View File

@ -1,4 +1,3 @@
#include <iostream>
#include <LatAnalyze/Derivative.hpp>
#include <LatAnalyze/CompiledFunction.hpp>
#include <LatAnalyze/Math.hpp>

88
examples/exFit.cpp
View File

@ -1,53 +1,77 @@
#include <iostream>
#include <cmath>
#include <LatAnalyze/CompiledModel.hpp>
#include <LatAnalyze/MinuitMinimizer.hpp>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/GslMinimizer.hpp>
#include <LatAnalyze/Plot.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <LatAnalyze/XYStatData.hpp>
using namespace std;
using namespace Latan;
const Index nPoint = 20;
const double exactPar[2] = {0.5,5.0}, dx = 10.0/static_cast<double>(nPoint);
const Index nPoint1 = 10, nPoint2 = 10;
const double xErr = .1, yErr = .3;
const double exactPar[2] = {0.5,5.};
const double dx1 = 10.0/static_cast<double>(nPoint1);
const double dx2 = 5.0/static_cast<double>(nPoint2);
int main(void)
{
// generate fake data
XYStatData data(nPoint, 1, 1);
RandGen rg;
double x_k, y_k;
DoubleModel f = compile("return p_1*exp(-x_0*p_0);", 1, 2);
for (Index k = 0; k < nPoint; ++k)
XYStatData data;
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis;
double xBuf[2];
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
cout << "-- generating fake data..." << endl;
data.addXDim(nPoint1);
data.addXDim(nPoint2);
data.addYDim();
for (Index i1 = 0; i1 < nPoint1; ++i1)
{
x_k = k*dx;
y_k = f(&x_k, exactPar) + rg.gaussian(0.0, 0.1);
cout << x_k << " " << y_k << " " << 0.1 << endl;
data.x(0, k) = x_k;
data.y(0, k) = y_k;
xBuf[0] = i1*dx1;
data.x(i1, 0) = xErr*dis(gen) + xBuf[0];
for (Index i2 = 0; i2 < nPoint2; ++i2)
{
xBuf[1] = i2*dx2;
data.x(i2, 1) = xBuf[1];
data.y(data.dataIndex(i1, i2), 0) = yErr*dis(gen)
+ f(xBuf, exactPar);
}
}
data.yyVar(0, 0).diagonal() = DMat::Constant(nPoint, 1, 0.1*0.1);
data.assumeXExact(0);
data.setXError(0, DVec::Constant(data.getXSize(0), xErr));
data.assumeXExact(true, 1);
data.setYError(0, DVec::Constant(data.getYSize(), yErr));
// set minimizers
DVec init = DVec::Constant(2, 0.1);
FitResult p;
GslMinimizer min(GslMinimizer::Algorithm::bfgs2);
// fit
DVec init = DVec::Constant(2, 0.5);
FitResult p;
MinuitMinimizer minimizer;
min.setVerbosity(Minimizer::Verbosity::Debug);
cout << "-- fit..." << endl;
f.parName().setName(0, "m");
f.parName().setName(1, "A");
p = data.fit(min, init, f);
p.print();
data.fitAllPoints();
p = data.fit(minimizer, init, f);
// plot
Plot plot;
DVec ref(2);
XYStatData res;
cout << "a= " << p(0) << " b= " << p(1);
cout << " chi^2/ndof= " << p.getChi2PerDof();
cout << " p-value= " << p.getPValue() <<endl;
// plot result
Plot plot;
plot << LogScale(Axis::y) << PlotData(data);
plot << Color("rgb 'blue'") << PlotFunction(p.getModel(), 0.0, 10.0);
cout << "-- generating plots..." << endl;
ref(1) = 0.;
res = data.getPartialResiduals(p, ref, 0);
plot << PlotRange(Axis::x, 0., 10.);
plot << Color("rgb 'blue'");
plot << PlotFunction(p.getModel().bind(0, ref), 0., 10.);
plot << Color("rgb 'red'");
plot << PlotData(res);
plot.display();
return EXIT_SUCCESS;

78
examples/exFitSample.cpp Normal file
View File

@ -0,0 +1,78 @@
#include <LatAnalyze/CompiledModel.hpp>
#include <LatAnalyze/GslMinimizer.hpp>
#include <LatAnalyze/Plot.hpp>
#include <LatAnalyze/XYSampleData.hpp>
using namespace std;
using namespace Latan;
const Index nPoint1 = 10, nPoint2 = 10;
const Index nSample = 1000;
const double xErr = .1, yErr = .3;
const double exactPar[2] = {0.5,5.};
const double dx1 = 10.0/static_cast<double>(nPoint1);
const double dx2 = 5.0/static_cast<double>(nPoint2);
int main(void)
{
// generate fake data
XYSampleData data(nSample);
double xBuf[2];
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis;
DoubleModel f([](const double *x, const double *p)
{return p[1]*exp(-x[0]*p[0])+x[1];}, 2, 2);
cout << "-- generating fake data..." << endl;
data.addXDim(nPoint1);
data.addXDim(nPoint2);
data.addYDim();
for (Index s = central; s < nSample; ++s)
{
for (Index i1 = 0; i1 < nPoint1; ++i1)
{
xBuf[0] = i1*dx1;
data.x(i1, 0)[s] = xErr*dis(gen) + xBuf[0];
for (Index i2 = 0; i2 < nPoint2; ++i2)
{
xBuf[1] = i2*dx2;
data.x(i2, 1)[s] = xBuf[1];
data.y(data.dataIndex(i1, i2), 0)[s] = yErr*dis(gen)
+ f(xBuf, exactPar);
}
}
}
data.assumeXExact(true, 1);
// set minimizers
DVec init = DVec::Constant(2, 0.1);
SampleFitResult p;
GslMinimizer min(GslMinimizer::Algorithm::bfgs2);
// fit
cout << "-- fit..." << endl;
f.parName().setName(0, "m");
f.parName().setName(1, "A");
p = data.fit(min, init, f);
p.print();
// plot
Plot plot;
DVec ref(2);
XYStatData res;
cout << "-- generating plots..." << endl;
ref(1) = 0.;
res = data.getPartialResiduals(p, ref, 0).getData();
plot << PlotRange(Axis::x, 0., 10.);
plot << Color("rgb 'blue'");
plot << PlotPredBand(p.getModel(_).bind(0, ref), 0., 10.);
plot << Color("rgb 'blue'");
plot << PlotFunction(p.getModel().bind(0, ref), 0., 10.);
plot << Color("rgb 'red'");
plot << PlotData(res);
plot.display();
return EXIT_SUCCESS;
}

1
examples/exIntegrator.cpp
View File

@ -1,4 +1,3 @@
#include <iostream>
#include <LatAnalyze/CompiledFunction.hpp>
#include <LatAnalyze/GslQagsIntegrator.hpp>

19
examples/exInterp.cpp
View File

@ -1,22 +1,3 @@
/*
* TabFunction.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli, Matt Spraggs
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/TabFunction.hpp>
int main(void)

3
examples/exMat.cpp
View File

@ -1,4 +1,3 @@
#include <iostream>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/Math.hpp>
@ -9,7 +8,7 @@ using namespace Latan;
int main(void)
{
DMat A(2, 3), B(3, 2), C = DMat::Random(6, 6);
const string fileName = "exMat.h5";
const string fileName = "exMat.dat";
A << 1, 2, 3,
4, 5, 6;

1
examples/exMathInterpreter.cpp
View File

@ -1,4 +1,3 @@
#include <iostream>
#include <LatAnalyze/Math.hpp>
#include <LatAnalyze/MathInterpreter.hpp>

7
examples/exMin.cpp
View File

@ -1,6 +1,5 @@
#include <iostream>
#include <LatAnalyze/CompiledFunction.hpp>
#include <LatAnalyze/MinuitMinimizer.hpp>
#include <LatAnalyze/GslMinimizer.hpp>
using namespace std;
using namespace Latan;
@ -17,8 +16,8 @@ int main(int argc, char* argv[])
}
source = argv[1];
DoubleFunction f = compile(source, 1);
MinuitMinimizer minimize;
DoubleFunction f = compile(source, 1);
GslMinimizer minimize;
DVec init(1);
double min;

1
examples/exPlot.cpp
View File

@ -1,4 +1,3 @@
#include <iostream>
#include <LatAnalyze/CompiledFunction.hpp>
#include <LatAnalyze/Math.hpp>
#include <LatAnalyze/Plot.hpp>

41
examples/exRand.cpp
View File

@ -1,53 +1,26 @@
#include <iostream>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/CompiledFunction.hpp>
#include <LatAnalyze/Plot.hpp>
#include <LatAnalyze/RandGen.hpp>
using namespace std;
using namespace Latan;
constexpr int seqLength = 25;
constexpr int saveStep = 9;
constexpr Index nDraw = 20000;
const string stateFileName = "exRand.seed";
int main(void)
{
RandGenState state;
RandGen gen[2];
AsciiFile stateFile(stateFileName, File::Mode::write|File::Mode::read);
DVec gauss(nDraw);
Plot p;
Histogram h;
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis;
DVec gauss(nDraw);
Plot p;
Histogram h;
cout << "- GENERATOR STATE I/O TESTS" << endl;
cout << "-- generating a " << seqLength << " steps random sequence..."
<< endl;
for (int i = 0; i < seqLength; ++i)
{
if (i == saveStep)
{
state = gen[0].getState();
stateFile.save(state, "exRand");
cout << "generator state after step " << saveStep - 1
<< " saved in '" << stateFileName << "'" << endl;
}
cout << "step " << i << "\t: " << gen[0].uniform() <<endl;
}
cout << "-- setting up another generator from '" << stateFileName << "'..."
<< endl;
gen[1].setState(stateFile.read<RandGenState>("exRand"));
cout << "-- generating a " << seqLength << " steps random sequence..."
<< endl;
for (int i = 0; i < seqLength; ++i)
{
cout << "step " << i << "\t: " << gen[1].uniform() <<endl;
}
cout << "-- generating " << nDraw << " Gaussian random numbers..." << endl;
FOR_VEC(gauss, i)
{
gauss(i) = gen[0].gaussian();
gauss(i) = dis(gen);
}
h.setFromData(gauss, -5., 5., 40);
h.normalize();

1
examples/exRootFinder.cpp
View File

@ -1,4 +1,3 @@
#include <iostream>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/GslHybridRootFinder.hpp>

41
lib/AsciiFile.cpp
View File

@ -1,7 +1,7 @@
/*
* AsciiFile.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -71,7 +71,22 @@ void AsciiFile::save(const DMat &m, const std::string &name)
fileStream_.precision(defaultPrec);
}
void AsciiFile::save(const DMatSample &s, const std::string &name)
void AsciiFile::save(const DSample &ds, const std::string &name)
{
if (name.empty())
{
LATAN_ERROR(Io, "trying to save data with an empty name");
}
checkWritability();
isParsed_ = false;
fileStream_ << "#L latan_begin rs_sample " << name << endl;
fileStream_ << ds.size() << endl;
save(ds.matrix(), name + "_data");
fileStream_ << "#L latan_end rs_sample " << endl;
}
void AsciiFile::save(const DMatSample &ms, const std::string &name)
{
if (name.empty())
{
@ -81,29 +96,15 @@ void AsciiFile::save(const DMatSample &s, const std::string &name)
checkWritability();
isParsed_ = false;
fileStream_ << "#L latan_begin rs_sample " << name << endl;
fileStream_ << s.size() << endl;
save(s[central], name + "_C");
for (Index i = 0; i < s.size(); ++i)
fileStream_ << ms.size() << endl;
save(ms[central], name + "_C");
for (Index i = 0; i < ms.size(); ++i)
{
save(s[i], name + "_S_" + strFrom(i));
save(ms[i], name + "_S_" + strFrom(i));
}
fileStream_ << "#L latan_end rs_sample " << endl;
}
void AsciiFile::save(const RandGenState &state, const std::string &name)
{
if (name.empty())
{
LATAN_ERROR(Io, "trying to save data with an empty name");
}
checkWritability();
isParsed_ = false;
fileStream_ << "#L latan_begin rg_state " << name << endl;
fileStream_ << state << endl;
fileStream_ << "#L latan_end rg_state " << endl;
}
// read first name ////////////////////////////////////////////////////////////
string AsciiFile::getFirstName(void)
{

12
lib/AsciiFile.hpp
View File

@ -1,7 +1,7 @@
/*
* AsciiFile.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -25,8 +25,6 @@
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/ParserState.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <fstream>
BEGIN_LATAN_NAMESPACE
@ -48,11 +46,11 @@ public:
bool isFirst;
std::string first;
// parsing buffers
RandGenState stateBuf;
int intBuf;
DSample dSampleBuf;
DMatSample dMatSampleBuf;
std::queue<DMat> dMatQueue;
std::queue<double> doubleQueue;
std::queue<int> intQueue;
private:
// allocation/deallocation functions defined in IoAsciiLexer.lpp
virtual void initScanner(void);
@ -66,8 +64,8 @@ public:
virtual ~AsciiFile(void);
// access
virtual void save(const DMat &m, const std::string &name);
virtual void save(const DMatSample &s, const std::string &name);
virtual void save(const RandGenState &state, const std::string &name);
virtual void save(const DSample &ds, const std::string &name);
virtual void save(const DMatSample &ms, const std::string &name);
// read first name
virtual std::string getFirstName(void);
// tests

17
lib/AsciiLexer.lpp
View File

@ -1,7 +1,7 @@
/*
* AsciiLexer.lpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -25,20 +25,9 @@
%option yylineno
%{
#include <iostream>
#include <LatAnalyze/AsciiFile.hpp>
#include "AsciiParser.hpp"
#if (defined __INTEL_COMPILER)
#pragma warning disable 1682 2259
#elif (defined __GNUC__)||(defined __clang__)
#pragma GCC diagnostic ignored "-Wconversion"
#pragma GCC diagnostic ignored "-Wsign-conversion"
#pragma GCC diagnostic ignored "-Wsign-compare"
#pragma GCC diagnostic ignored "-Wunused-function"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#endif
using namespace std;
using namespace Latan;
@ -76,7 +65,7 @@ BLANK [ \t]
%%
{LMARK} {BEGIN(MARK);}
{INT} {yylval->val_int = strTo<int>(yytext); RETTOK(INT);}
{INT} {yylval->val_int = strTo<long int>(yytext); RETTOK(INT);}
{FLOAT} {yylval->val_double = strTo<double>(yytext); RETTOK(FLOAT);}
({ALPHA}|{DIGIT})+ {strcpy(yylval->val_str,yytext); RETTOK(ID);}
<MARK>latan_begin {BEGIN(TYPE); RETTOK(OPEN);}

89
lib/AsciiParser.ypp
View File

@ -1,7 +1,7 @@
/*
* AsciiParser.ypp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -22,18 +22,6 @@
#include <LatAnalyze/AsciiFile.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <iostream>
#include <sstream>
#include <utility>
#include <cstring>
#if (defined __INTEL_COMPILER)
#pragma warning disable 1682 2259
#elif (defined __GNUC__)||(defined __clang__)
#pragma GCC diagnostic ignored "-Wconversion"
#pragma GCC diagnostic ignored "-Wsign-conversion"
#endif
using namespace std;
using namespace Latan;
@ -47,7 +35,7 @@
%}
%pure-parser
%name-prefix = "_Ascii_"
%name-prefix "_Ascii_"
%locations
%defines
%error-verbose
@ -57,10 +45,10 @@
%union
{
int val_int;
double val_double;
char val_char;
char val_str[256];
long int val_int;
double val_double;
char val_char;
char val_str[256];
}
%token <val_char> ERR
@ -69,7 +57,7 @@
%token <val_str> ID
%token OPEN CLOSE MAT SAMPLE RG_STATE
%type <val_str> mat sample rg_state
%type <val_str> mat matsample dsample
%{
int _Ascii_lex(YYSTYPE* lvalp, YYLTYPE* llocp, void* scanner);
@ -101,16 +89,16 @@ data:
(*state->data)[$1].reset(new DMat(state->dMatQueue.front()));
state->dMatQueue.pop();
}
| sample
| dsample
{
TEST_FIRST($1);
(*state->data)[$1].reset(new DSample(state->dSampleBuf));
}
| matsample
{
TEST_FIRST($1);
(*state->data)[$1].reset(new DMatSample(state->dMatSampleBuf));
}
| rg_state
{
TEST_FIRST($1);
(*state->data)[$1].reset(new RandGenState(state->stateBuf));
}
;
mat:
@ -138,15 +126,37 @@ mat:
}
;
sample:
OPEN SAMPLE ID INT mats CLOSE SAMPLE
dsample:
OPEN SAMPLE ID INT mat CLOSE SAMPLE
{
const unsigned int nSample = $4, os = DMatSample::offset;
DMat &m = state->dMatQueue.front();
if (m.rows() != nSample + os)
{
LATAN_ERROR(Size, "double sample '" + *state->streamName + ":"
+ $3 + "' has a wrong size");
}
if (m.cols() != 1)
{
LATAN_ERROR(Size, "double sample '" + *state->streamName + ":"
+ $3 + "' is not stored as a column vector");
}
state->dSampleBuf = m.array();
state->dMatQueue.pop();
strcpy($$, $3);
}
;
matsample:
OPEN SAMPLE ID INT mat mats CLOSE SAMPLE
{
const unsigned int nSample = $4, os = DMatSample::offset;
if (state->dMatQueue.size() != nSample + os)
{
LATAN_ERROR(Size, "sample '" + *state->streamName + ":" + $3 +
"' has a wrong size");
LATAN_ERROR(Size, "matrix sample '" + *state->streamName + ":"
+ $3 + "' has a wrong size");
}
state->dMatSampleBuf.resize(nSample);
state->dMatSampleBuf[central] = state->dMatQueue.front();
@ -158,22 +168,6 @@ sample:
}
strcpy($$, $3);
}
rg_state:
OPEN RG_STATE ID ints CLOSE RG_STATE
{
if (state->intQueue.size() != RLXG_STATE_SIZE)
{
LATAN_ERROR(Size, "random generator state '" + *state->streamName
+ ":" + $3 + "' has a wrong size");
}
for (Index i = 0; i < RLXG_STATE_SIZE; ++i)
{
state->stateBuf[i] = state->intQueue.front();
state->intQueue.pop();
}
strcpy($$, $3);
}
;
mats:
@ -187,8 +181,3 @@ floats:
| FLOAT {state->doubleQueue.push($1);}
| INT {state->doubleQueue.push(static_cast<double>($1));}
;
ints:
ints INT {state->intQueue.push($2);}
| INT {state->intQueue.push($1);}
;

325
lib/Chi2Function.cpp
View File

@ -1,325 +0,0 @@
/*
* Chi2Function.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/Chi2Function.hpp>
#include <LatAnalyze/includes.hpp>
#include <LatAnalyze/XYStatData.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* Chi2Function implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
Chi2Function::Chi2Function(const XYStatData &data)
: data_(&data)
, buffer_(new Chi2FunctionBuffer)
{
resizeBuffer();
}
Chi2Function::Chi2Function(const XYStatData &data,
const vector<const DoubleModel *> &modelVector)
: Chi2Function(data)
{
setModel(modelVector);
}
// access //////////////////////////////////////////////////////////////////////
Index Chi2Function::getNArg(void) const
{
if (nPar_ < 0)
{
LATAN_ERROR(Memory, "no model set");
}
return nPar_ + data_->getStatXDim()*data_->getNFitPoint();
}
Index Chi2Function::getNDof(void) const
{
if (nPar_ < 0)
{
LATAN_ERROR(Memory, "no model set");
}
return data_->getYDim()*data_->getNFitPoint() - nPar_;
}
Index Chi2Function::getNPar(void) const
{
if (nPar_ < 0)
{
LATAN_ERROR(Memory, "no model set");
}
return nPar_;
}
void Chi2Function::setModel(const DoubleModel &model, const Index j)
{
typedef decltype(model_.size()) size_type;
if (static_cast<Index>(model_.size()) != data_->getYDim())
{
model_.resize(static_cast<size_type>(data_->getYDim()));
}
if (model.getNArg() != data_->getXDim())
{
LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
}
for (unsigned int l = 0; l < data_->getYDim(); ++l)
{
if (model_[l]&&(l != j))
{
if (model_[l]->getNPar() != model.getNPar())
{
LATAN_ERROR(Size, "model number of parameter mismatch");
}
}
}
model_[static_cast<size_type>(j)] = &model;
nPar_ = model.getNPar();
}
void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
{
typedef decltype(model_.size()) size_type;
if (static_cast<Index>(model_.size()) != data_->getYDim())
{
model_.resize(static_cast<size_type>(data_->getYDim()));
}
if (modelVector.size() != static_cast<size_type>(data_->getYDim()))
{
LATAN_ERROR(Size, "number of models and y-dimension mismatch");
}
for (unsigned int j = 0; j < model_.size(); ++j)
{
if (!modelVector[j])
{
LATAN_ERROR(Memory, "trying to set a null model");
}
if (modelVector[j]->getNArg() != data_->getXDim())
{
LATAN_ERROR(Size, "model number of arguments and x-dimension mismatch");
}
model_[static_cast<size_type>(j)] = modelVector[j];
if (modelVector[j]->getNPar() != modelVector[0]->getNPar())
{
LATAN_ERROR(Size, "model number of parameter mismatch");
}
}
nPar_ = modelVector[0]->getNPar();
}
void Chi2Function::requestInit(void) const
{
buffer_->isInit = false;
}
void Chi2Function::resizeBuffer(void) const
{
Index size;
size = (data_->getYDim() + data_->getStatXDim())*data_->getNFitPoint();
buffer_->v.setConstant(size, 0.0);
buffer_->x.setConstant(data_->getXDim(), 0.0);
buffer_->invVar.setConstant(size, size, 0.0);
buffer_->xInd.setConstant(data_->getStatXDim(), 0);
buffer_->dInd.setConstant(data_->getNFitPoint(), 0);
}
// compute variance matrix inverse /////////////////////////////////////////////
void Chi2Function::setVarianceBlock(const Index l1, const Index l2,
ConstBlock<MatBase<double>> m) const
{
const Index nPoint = data_->getNFitPoint();
FOR_VEC(buffer_->dInd, k2)
FOR_VEC(buffer_->dInd, k1)
{
if (data_->isDataCorrelated(buffer_->dInd(k1), buffer_->dInd(k2)))
{
buffer_->invVar(l1*nPoint + k1, l2*nPoint + k2) =
m(buffer_->dInd(k1), buffer_->dInd(k2));
}
}
}
void Chi2Function::initBuffer(void) const
{
const Index xDim = data_->getXDim();
const Index statXDim = data_->getStatXDim();
const Index yDim = data_->getYDim();
const Index nData = data_->getNData();
const Index nPoint = data_->getNFitPoint();
Index is, kf;
// resize buffer
resizeBuffer();
// build index tables
is = 0;
for (Index i = 0; i < xDim; ++i)
{
if (!data_->isXExact(i))
{
buffer_->xInd(is) = i;
is++;
}
}
kf = 0;
for (Index k = 0; k < nData; ++k)
{
if (data_->isFitPoint(k))
{
buffer_->dInd(kf) = k;
kf++;
}
}
// set y/y variance matrix
for (Index j2 = 0; j2 < yDim; ++j2)
for (Index j1 = 0; j1 < yDim; ++j1)
{
if (data_->isYYCorrelated(j1, j2))
{
setVarianceBlock(j1, j2, data_->yyVar(j1, j2));
}
}
// set x/x variance matrix
FOR_VEC(buffer_->xInd, i2)
FOR_VEC(buffer_->xInd, i1)
{
if (data_->isXXCorrelated(buffer_->xInd(i1), buffer_->xInd(i2)))
{
setVarianceBlock(i1 + yDim, i2 + yDim,
data_->xxVar(buffer_->xInd(i1), buffer_->xInd(i2)));
}
}
// set y/x variance matrix
FOR_VEC(buffer_->xInd, i)
for (Index j = 0; j < yDim; ++j)
{
if (data_->isYXCorrelated(j, buffer_->xInd(i)))
{
setVarianceBlock(j, i + yDim, data_->yxVar(j, buffer_->xInd(i)));
}
}
auto lowerYX = buffer_->invVar.block(yDim*nPoint, 0, yDim*statXDim,
yDim*nPoint);
auto upperYX = buffer_->invVar.block(0, yDim*nPoint, yDim*nPoint,
yDim*statXDim);
lowerYX = upperYX.transpose();
// inversion
buffer_->invVar = buffer_->invVar.pInverse(data_->getSvdTolerance());
buffer_->isInit = true;
}
// function call ///////////////////////////////////////////////////////////////
double Chi2Function::operator()(const double *arg) const
{
typedef decltype(model_.size()) size_type;
if (!model_[0])
{
LATAN_ERROR(Memory, "null model");
}
const Index xDim = data_->getXDim();
const Index yDim = data_->getYDim();
const Index nPoint = data_->getNFitPoint();
Index is;
ConstMap<DVec> xi(arg + nPar_, getNArg() - nPar_, 1);
double res;
// initialize buffers if necessary
if (!buffer_->isInit)
{
initBuffer();
}
// set the upper part of v: f_j(xi, p) - y_j
for (Index j = 0; j < yDim; ++j)
FOR_VEC(buffer_->dInd, k)
{
const DoubleModel *f = model_[static_cast<size_type>(j)];
double f_jk, y_jk = data_->y(j, buffer_->dInd(k));
if (!f)
{
LATAN_ERROR(Memory, "null model");
}
is = 0;
for (Index i = 0; i < xDim; ++i)
{
if (data_->isXExact(i))
{
buffer_->x(i) = data_->x(i, buffer_->dInd(k));
}
else
{
buffer_->x(i) = xi(is*nPoint + k);
is++;
}
}
// call model on double pointers to avoid any copy
f_jk = (*f)(buffer_->x.data(), arg);
buffer_->v(j*nPoint + k) = f_jk - y_jk;
}
// set the down part of v: xi_ik - x_ik
FOR_VEC(buffer_->xInd, i)
FOR_VEC(buffer_->dInd, k)
{
double x_ik = data_->x(buffer_->xInd(i), buffer_->dInd(k));
double xi_ik = xi(i*nPoint + k);
buffer_->v(yDim*nPoint + i*nPoint + k) = xi_ik - x_ik;
}
// compute result
res = buffer_->v.dot(buffer_->invVar*buffer_->v);
return res;
}
// DoubleFunction factory //////////////////////////////////////////////////////
DoubleFunction Chi2Function::makeFunction(const bool makeHardCopy) const
{
DoubleFunction res;
if (makeHardCopy)
{
Chi2Function copy(*this);
res.setFunction([copy](const double *p){return copy(p);}, getNArg());
}
else
{
res.setFunction([this](const double *p){return (*this)(p);}, getNArg());
}
return res;
}

80
lib/Chi2Function.hpp
View File

@ -1,80 +0,0 @@
/*
* Chi2Function.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_Chi2Function_hpp_
#define Latan_Chi2Function_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Model.hpp>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* <Class> *
******************************************************************************/
// forward declaration of XYStatData
class XYStatData;
class Chi2Function: public DoubleFunctionFactory
{
private:
struct Chi2FunctionBuffer
{
DVec v, x;
DMat invVar;
bool isInit{false};
Vec<Index> xInd, dInd;
};
public:
// constructor
explicit Chi2Function(const XYStatData &data);
Chi2Function(const XYStatData &data,
const std::vector<const DoubleModel *> &modelVector);
// destructor
virtual ~Chi2Function(void) = default;
// access
virtual Index getNArg(void) const;
Index getNDof(void) const;
Index getNPar(void) const;
void setModel(const DoubleModel &model, const Index j = 0);
void setModel(const std::vector<const DoubleModel *> &modelVector);
void requestInit(void) const;
// function call
double operator()(const double *arg) const;
// factory
virtual DoubleFunction makeFunction(const bool makeHardCopy = true) const;
private:
// access
void resizeBuffer(void) const;
// compute variance matrix inverse
void setVarianceBlock(const Index l1, const Index l2,
ConstBlock<MatBase<double>> m) const;
void initBuffer(void) const;
private:
const XYStatData *data_;
std::shared_ptr<Chi2FunctionBuffer> buffer_;
std::vector<const DoubleModel *> model_;
Index nPar_{-1};
};
END_LATAN_NAMESPACE
#endif // Latan_Chi2Function_hpp_

2
lib/CompiledFunction.cpp
View File

@ -1,7 +1,7 @@
/*
* CompiledFunction.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/CompiledFunction.hpp
View File

@ -1,7 +1,7 @@
/*
* CompiledFunction.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/CompiledModel.cpp
View File

@ -1,7 +1,7 @@
/*
* CompiledModel.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/CompiledModel.hpp
View File

@ -1,7 +1,7 @@
/*
* CompiledModel.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

37
lib/Dataset.hpp
View File

@ -1,7 +1,7 @@
/*
* Dataset.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,9 +23,6 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/File.hpp>
#include <LatAnalyze/StatArray.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <fstream>
#include <vector>
BEGIN_LATAN_NAMESPACE
@ -35,6 +32,8 @@ BEGIN_LATAN_NAMESPACE
template <typename T>
class Dataset: public StatArray<T>
{
public:
typedef std::random_device::result_type SeedType;
public:
// constructors
Dataset(void) = default;
@ -46,7 +45,8 @@ public:
template <typename FileType>
void load(const std::string &listFileName, const std::string &dataName);
// resampling
Sample<T> bootstrapMean(const Index nSample, RandGen& generator);
Sample<T> bootstrapMean(const Index nSample, const SeedType seed);
Sample<T> bootstrapMean(const Index nSample);
private:
// mean from pointer vector for resampling
void ptVectorMean(T &m, const std::vector<const T *> &v);
@ -82,26 +82,23 @@ void Dataset<T>::load(const std::string &listFileName,
// resampling //////////////////////////////////////////////////////////////////
template <typename T>
Sample<T> Dataset<T>::bootstrapMean(const Index nSample, RandGen& generator)
Sample<T> Dataset<T>::bootstrapMean(const Index nSample, const SeedType seed)
{
typedef typename std::vector<const T *>::size_type size_type;
size_type nData = static_cast<size_type>(this->size());
std::vector<const T *> data(nData);
Sample<T> s(nSample);
std::vector<const T *> data(this->size());
Sample<T> s(nSample);
std::mt19937 gen(seed);
std::uniform_int_distribution<Index> dis(0, this->size() - 1);
for (size_type j = 0; j < nData; ++j)
for (unsigned int j = 0; j < this->size(); ++j)
{
data[j] = &((*this)[static_cast<Index>(j)]);
}
ptVectorMean(s[central], data);
for (Index i = 0; i < nSample; ++i)
{
for (size_type j = 0; j < nData; ++j)
for (unsigned int j = 0; j < this->size(); ++j)
{
Index k= static_cast<Index>(generator.discreteUniform(static_cast<unsigned int>(nData)));
data[j] = &((*this)[k]);
data[j] = &((*this)[dis(gen)]);
}
ptVectorMean(s[i], data);
}
@ -109,6 +106,14 @@ Sample<T> Dataset<T>::bootstrapMean(const Index nSample, RandGen& generator)
return s;
}
template <typename T>
Sample<T> Dataset<T>::bootstrapMean(const Index nSample)
{
std::random_device rd;
return bootstrapMean(nSample, rd());
}
template <typename T>
void Dataset<T>::ptVectorMean(T &m, const std::vector<const T *> &v)
{

23
lib/Derivative.cpp
View File

@ -1,7 +1,7 @@
/*
* Derivative.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -31,11 +31,12 @@ using namespace Math;
// constructor /////////////////////////////////////////////////////////////////
Derivative::Derivative(const DoubleFunction &f, const Index dir,
const double step)
: f_(f)
, dir_(dir)
, step_(step)
, buffer_(new DVec(f.getNArg()))
{}
: buffer_(new DVec(f.getNArg()))
{
setFunction(f);
setDir(dir);
setStep(step);
}
Derivative::Derivative(const DoubleFunction &f, const Index dir,
const Index order, const DVec &point, const double step)
@ -45,6 +46,11 @@ Derivative::Derivative(const DoubleFunction &f, const Index dir,
}
// access //////////////////////////////////////////////////////////////////////
Index Derivative::getDir(void) const
{
return dir_;
}
Index Derivative::getOrder(void) const
{
return order_;
@ -60,6 +66,11 @@ double Derivative::getStep(void) const
return step_;
}
void Derivative::setDir(const Index dir)
{
dir_ = dir;
}
void Derivative::setFunction(const DoubleFunction &f)
{
f_ = f;

11
lib/Derivative.hpp
View File

@ -1,7 +1,7 @@
/*
* Derivative.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -39,9 +39,11 @@ public:
// destructor
virtual ~Derivative(void) = default;
// access
Index getDir(void) const;
Index getNPoint(void) const;
Index getOrder(void) const;
double getStep(void) const;
void setDir(const Index dir);
void setFunction(const DoubleFunction &f);
void setOrderAndPoint(const Index order, const DVec &point);
void setStep(const double step);
@ -73,16 +75,13 @@ public:
static const Index defaultPrecOrder = 2;
public:
// constructor
CentralDerivative(const DoubleFunction &f, const Index dir = 0,
CentralDerivative(const DoubleFunction &f = DoubleFunction(),
const Index dir = 0,
const Index order = 1,
const Index precOrder = defaultPrecOrder);
// destructor
virtual ~CentralDerivative(void) = default;
// access
using Derivative::getNPoint;
using Derivative::getStep;
using Derivative::getOrder;
using Derivative::setStep;
Index getPrecOrder(void) const;
void setOrder(const Index order, const Index precOrder = defaultPrecOrder);
// function call

3
lib/Exceptions.cpp
View File

@ -1,7 +1,7 @@
/*
* Exceptions.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -41,6 +41,7 @@ CONST_EXC(Range, Logic("range error: " + msg, loc))
CONST_EXC(Size, Logic("size error: " + msg, loc))
// runtime errors
CONST_EXC(Runtime, runtime_error(Env::msgPrefix + msg + ERR_SUFF))
CONST_EXC(Argument, Runtime("argument error: " + msg, loc))
CONST_EXC(Compilation, Runtime("compilation error: " + msg, loc))
CONST_EXC(Io, Runtime("IO error: " + msg, loc))
CONST_EXC(Memory, Runtime("memory error: " + msg, loc))

3
lib/Exceptions.hpp
View File

@ -1,7 +1,7 @@
/*
* Exceptions.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -51,6 +51,7 @@ namespace Exceptions
DECL_EXC(Size, Logic);
// runtime errors
DECL_EXC(Runtime, std::runtime_error);
DECL_EXC(Argument, Runtime);
DECL_EXC(Compilation, Runtime);
DECL_EXC(Io, Runtime);
DECL_EXC(Memory, Runtime);

6
lib/File.cpp
View File

@ -1,7 +1,7 @@
/*
* File.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -52,10 +52,6 @@ unsigned int File::getMode(void) const
// internal functions //////////////////////////////////////////////////////////
void File::deleteData(void)
{
for (auto &i : data_)
{
i.second.reset();
}
data_.clear();
}

9
lib/File.hpp
View File

@ -1,7 +1,7 @@
/*
* File.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -24,9 +24,6 @@
#include <LatAnalyze/IoObject.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <queue>
#include <unordered_map>
BEGIN_LATAN_NAMESPACE
@ -61,8 +58,8 @@ public:
template <typename IoT>
const IoT & read(const std::string &name = "");
virtual void save(const DMat &m, const std::string &name) = 0;
virtual void save(const DMatSample &state, const std::string &name) = 0;
virtual void save(const RandGenState &state, const std::string &name) = 0;
virtual void save(const DSample &ds, const std::string &name) = 0;
virtual void save(const DMatSample &ms, const std::string &name) = 0;
// read first name
virtual std::string getFirstName(void) = 0;
// tests

796
lib/FitInterface.cpp
View File

@ -1,7 +1,7 @@
/*
* FitInterface.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -24,186 +24,770 @@ using namespace std;
using namespace Latan;
/******************************************************************************
* FitInterface implementation *
* FitInterface implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
FitInterface::FitInterface(const Index nData, const Index xDim,
const Index yDim)
// constructor /////////////////////////////////////////////////////////////////
FitInterface::FitInterface(void)
: xName_("x")
, yName_("y")
{}
// copy object (not as a constructor to be accessed from derived class) ////////
void FitInterface::copyInterface(const FitInterface &d)
{
resize(nData, xDim, yDim);
*this = d;
scheduleFitVarMatInit();
}
// access //////////////////////////////////////////////////////////////////////
void FitInterface::assumeXExact(const Index i, const bool isExact)
// add dimensions //////////////////////////////////////////////////////////////
void FitInterface::addXDim(const Index nData, const string name,
const bool isExact)
{
isXExact_(i) = (isExact) ? 1 : 0;
}
void FitInterface::assumeXXCorrelated(const Index i1, const Index i2,
const bool isCorrelated)
{
int val = (isCorrelated) ? 1 : 0;
isXXCorr_(i1, i2) = val;
isXXCorr_(i2, i1) = val;
}
void FitInterface::assumeYYCorrelated(const Index j1, const Index j2,
const bool isCorrelated)
{
int val = (isCorrelated) ? 1 : 0;
isYYCorr_(j1, j2) = val;
isYYCorr_(j2, j1) = val;
}
void FitInterface::assumeYXCorrelated(const Index j, const Index i,
const bool isCorrelated)
{
int val = (isCorrelated) ? 1 : 0;
isYXCorr_(j, i) = val;
}
void FitInterface::assumeDataCorrelated(const Index k1, const Index k2,
const bool isCorrelated)
{
int val = (isCorrelated) ? 1 : 0;
isDataCorr_(k1, k2) = val;
isDataCorr_(k2, k1) = val;
}
void FitInterface::assumeDataCorrelated(const bool isCorrelated)
{
FOR_MAT(isDataCorr_, k1, k2)
if (getYSize() != 0)
{
if (k1 != k2)
LATAN_ERROR(Logic, "cannot add an X dimension if fit data is "
"not empty");
}
else
{
xSize_.push_back(nData);
xIsExact_.push_back(isExact);
maxDataIndex_ *= nData;
createXData(name, nData);
scheduleLayoutInit();
scheduleDataCoordInit();
if (!name.empty())
{
assumeDataCorrelated(k1, k2, isCorrelated);
xName().setName(getNXDim(), name);
}
}
}
void FitInterface::fitPoint(const Index i, const bool isFitPoint)
void FitInterface::addYDim(const string name)
{
isFitPoint_(i) = (isFitPoint) ? 1 : 0;
yDataIndex_.push_back(map<Index, bool>());
createYData(name);
scheduleLayoutInit();
if (!name.empty())
{
yName().setName(getNYDim(), name);
}
}
void FitInterface::fitPointRange(const Index k1, const Index k2,
const bool isFitPoint)
// access //////////////////////////////////////////////////////////////////////
Index FitInterface::getNXDim(void) const
{
int size = static_cast<int>(k2-k1+1);
return xSize_.size();
}
Index FitInterface::getNYDim(void) const
{
return yDataIndex_.size();
}
Index FitInterface::getXSize(void) const
{
Index size = 0;
isFitPoint_.segment(k1, size) = IVec::Constant(size, (isFitPoint) ? 1 : 0);
for (Index i = 0; i < getNXDim(); ++i)
{
size += getXSize(i);
}
return size;
}
void FitInterface::fitAllPoints(const bool isFitPoint)
Index FitInterface::getXSize(const Index i) const
{
fitPointRange(0, getNData()-1, isFitPoint);
checkXDim(i);
return xSize_[i];
}
Index FitInterface::getNData(void) const
Index FitInterface::getYSize(void) const
{
return isFitPoint_.size();
Index size = 0;
for (Index j = 0; j < getNYDim(); ++j)
{
size += getYSize(j);
}
return size;
}
Index FitInterface::getNFitPoint(void) const
Index FitInterface::getYSize(const Index j) const
{
return isFitPoint_.sum();
checkYDim(j);
return static_cast<Index>(yDataIndex_[j].size());
}
Index FitInterface::getXDim(void) const
Index FitInterface::getXFitSize(void) const
{
return isXXCorr_.rows();
Index size = 0;
for (Index i = 0; i < getNXDim(); ++i)
{
size += getXFitSize(i);
}
return size;
}
Index FitInterface::getYDim(void) const
Index FitInterface::getXFitSize(const Index i) const
{
return isYYCorr_.rows();
set<Index> fitCoord;
vector<Index> v;
checkXDim(i);
for (Index j = 0; j < getNYDim(); ++j)
{
for (auto &p: yDataIndex_[j])
{
if (p.second)
{
v = dataCoord(p.first);
fitCoord.insert(v[i]);
}
}
}
return fitCoord.size();
}
Index FitInterface::getStatXDim(void) const
Index FitInterface::getYFitSize(void) const
{
return isXExact_.size() - isXExact_.sum();
Index size = 0;
for (Index j = 0; j < getNYDim(); ++j)
{
size += getYFitSize(j);
}
return size;
}
Index FitInterface::getYFitSize(const Index j) const
{
Index size;
auto pred = [](const pair<Index, bool> &p)
{
return p.second;
};
checkYDim(j);
size = count_if(yDataIndex_[j].begin(), yDataIndex_[j].end(), pred);
return size;
}
Index FitInterface::getMaxDataIndex(void) const
{
return maxDataIndex_;
}
const set<Index> & FitInterface::getDataIndexSet(void) const
{
return dataIndexSet_;
}
double FitInterface::getSvdTolerance(void) const
{
return svdTolerance_;
return svdTol_;
}
void FitInterface::setFitInterface(const FitInterface &fitInterface)
void FitInterface::setSvdTolerance(const double &tol)
{
if (&fitInterface != this)
svdTol_ = tol;
}
VarName & FitInterface::xName(void)
{
return xName_;
}
const VarName & FitInterface::xName(void) const
{
return xName_;
}
VarName & FitInterface::yName(void)
{
return yName_;
}
const VarName & FitInterface::yName(void) const
{
return yName_;
}
// Y dimension index helper ////////////////////////////////////////////////////
Index FitInterface::dataIndex(const vector<Index> &v) const
{
Index k, n = v.size();
checkDataCoord(v);
k = xSize_[1]*v[0];
for (unsigned int d = 1; d < n-1; ++d)
{
isXExact_ = fitInterface.isXExact_;
isFitPoint_ = fitInterface.isFitPoint_;
isXXCorr_ = fitInterface.isXXCorr_;
isYYCorr_ = fitInterface.isYYCorr_;
isYXCorr_ = fitInterface.isYXCorr_;
isDataCorr_ = fitInterface.isDataCorr_;
svdTolerance_ = fitInterface.svdTolerance_;
k = xSize_[d+1]*(v[d] + k);
}
k += v[n-1];
return k;
}
const vector<Index> & FitInterface::dataCoord(const Index k) const
{
checkDataIndex(k);
updateDataCoord();
return dataCoord_.at(k);
}
// enable fit points ///////////////////////////////////////////////////////////
void FitInterface::fitPoint(const bool isFitPoint, const Index k, const Index j)
{
checkPoint(k, j);
yDataIndex_[j][k] = isFitPoint;
scheduleLayoutInit();
}
// variance interface //////////////////////////////////////////////////////////
void FitInterface::assumeXExact(const bool isExact, const Index i)
{
checkXDim(i);
xIsExact_[i] = isExact;
scheduleLayoutInit();
}
void FitInterface::addCorr(set<array<Index, 4>> &s, const bool isCorr,
const array<Index, 4> &c)
{
if (isCorr)
{
s.insert(c);
}
else
{
auto it = s.find(c);
if (it != s.end())
{
s.erase(it);
}
}
}
void FitInterface::setNData(const Index nData)
void FitInterface::assumeXXCorrelated(const bool isCorr, const Index r1,
const Index i1, const Index r2,
const Index i2)
{
resize(nData, getXDim(), getYDim());
array<Index, 4> c{{r1, i1, r2, i2}};
checkXIndex(r1, i1);
checkXIndex(r2, i2);
if ((i1 != i2) or (r1 != r2))
{
addCorr(xxCorr_, isCorr, c);
}
scheduleFitVarMatInit();
}
void FitInterface::setXDim(const Index xDim)
void FitInterface::assumeXXCorrelated(const bool isCorr, const Index i1,
const Index i2)
{
resize(getNData(), xDim, getYDim());
for (Index r1 = 0; r1 < getXSize(i1); ++r1)
for (Index r2 = 0; r2 < getXSize(i2); ++r2)
{
assumeXXCorrelated(isCorr, r1, i1, r2, i2);
}
}
void FitInterface::setYDim(const Index yDim)
void FitInterface::assumeYYCorrelated(const bool isCorr, const Index k1,
const Index j1, const Index k2,
const Index j2)
{
resize(getNData(), getXDim(), yDim);
array<Index, 4> c{{k1, j1, k2, j2}};
checkPoint(k1, j1);
checkPoint(k2, j2);
if ((j1 != j2) or (k1 != k2))
{
addCorr(yyCorr_, isCorr, c);
}
scheduleFitVarMatInit();
}
void FitInterface::setSvdTolerance(const double tolerance)
void FitInterface::assumeYYCorrelated(const bool isCorr, const Index j1,
const Index j2)
{
svdTolerance_ = tolerance;
checkYDim(j1);
checkYDim(j2);
for (auto &p1: yDataIndex_[j1])
for (auto &p2: yDataIndex_[j2])
{
assumeYYCorrelated(isCorr, p1.first, j1, p2.first, j2);
}
}
void FitInterface::resize(const Index nData, const Index xDim, const Index yDim)
void FitInterface::assumeXYCorrelated(const bool isCorr, const Index r,
const Index i, const Index k,
const Index j)
{
isXExact_.setConstant(xDim, 0);
isFitPoint_.setConstant(nData, 0);
isXXCorr_.setIdentity(xDim, xDim);
isYYCorr_.setIdentity(yDim, yDim);
isYXCorr_.setConstant(yDim, xDim, 0);
isDataCorr_.setIdentity(nData, nData);
array<Index, 4> c{{r, i, k, j}};
checkXIndex(r, i);
checkPoint(k, j);
addCorr(xyCorr_, isCorr, c);
scheduleFitVarMatInit();
}
// test ////////////////////////////////////////////////////////////////////////
bool FitInterface::isFitPoint(const Index k) const
void FitInterface::assumeXYCorrelated(const bool isCorr, const Index i,
const Index j)
{
return (isFitPoint_[k] == 1);
checkYDim(j);
for (Index r = 0; r < getXSize(i); ++r)
for (auto &p: yDataIndex_[j])
{
assumeXYCorrelated(isCorr, r, i, p.first, j);
}
}
// tests ///////////////////////////////////////////////////////////////////////
bool FitInterface::pointExists(const Index k) const
{
bool isUsed = false;
for (Index j = 0; j < getNYDim(); ++j)
{
isUsed = isUsed or pointExists(k, j);
}
return isUsed;
}
bool FitInterface::pointExists(const Index k, const Index j) const
{
checkDataIndex(k);
checkYDim(j);
return !(yDataIndex_[j].find(k) == yDataIndex_[j].end());
}
bool FitInterface::isXExact(const Index i) const
{
return (isXExact_[i] == 1);
checkXDim(i);
return xIsExact_[i];
}
bool FitInterface::isXXCorrelated(const Index i1, const Index i2) const
bool FitInterface::isXUsed(const Index r, const Index i, const bool inFit) const
{
return (isXXCorr_(i1, i2) == 1);
vector<Index> v;
checkXDim(i);
for (Index j = 0; j < getNYDim(); ++j)
{
for (auto &p: yDataIndex_[j])
{
if (p.second or !inFit)
{
v = dataCoord(p.first);
if (v[i] == r)
{
return true;
}
}
}
}
return false;
}
bool FitInterface::isYYCorrelated(const Index j1, const Index j2) const
bool FitInterface::isFitPoint(const Index k, const Index j) const
{
return (isYYCorr_(j1, j2) == 1);
checkPoint(k, j);
return yDataIndex_[j].at(k);
}
bool FitInterface::isYXCorrelated(const Index j, const Index i) const
bool FitInterface::isXXCorrelated(const Index r1, const Index i1,
const Index r2, const Index i2) const
{
return (isYXCorr_(j, i) == 1);
array<Index, 4> c{{r1, i1, r2, i2}};
auto it = xxCorr_.find(c);
return (it != xxCorr_.end());
}
bool FitInterface::isDataCorrelated(const Index k1, const Index k2) const
bool FitInterface::isYYCorrelated(const Index k1, const Index j1,
const Index k2, const Index j2) const
{
return (isDataCorr_(k1, k2) == 1);
array<Index, 4> c{{k1, j1, k2, j2}};
auto it = yyCorr_.find(c);
return (it != yyCorr_.end());
}
bool FitInterface::isXYCorrelated(const Index r, const Index i,
const Index k, const Index j) const
{
array<Index, 4> c{{r, i, k, j}};
auto it = xyCorr_.find(c);
return (it != xyCorr_.end());
}
bool FitInterface::hasCorrelations(void) const
{
return ((xxCorr_.size() != 0) or (yyCorr_.size() != 0)
or (xyCorr_.size() != 0));
}
// make correlation filter for fit variance matrix /////////////////////////////
DMat FitInterface::makeCorrFilter(void)
{
updateLayout();
DMat f = DMat::Identity(layout.totalSize, layout.totalSize);
Index row, col;
for (auto &c: xxCorr_)
{
row = indX(c[0], c[1]);
col = indX(c[2], c[3]);
if ((row != -1) and (col != -1))
{
f(row, col) = 1.;
f(col, row) = 1.;
}
}
for (auto &c: yyCorr_)
{
row = indY(c[0], c[1]);
col = indY(c[2], c[3]);
if ((row != -1) and (col != -1))
{
f(row, col) = 1.;
f(col, row) = 1.;
}
}
for (auto &c: xyCorr_)
{
row = indX(c[0], c[1]);
col = indY(c[2], c[3]);
if ((row != -1) and (col != -1))
{
f(row, col) = 1.;
f(col, row) = 1.;
}
}
return f;
}
// schedule variance matrix initialization /////////////////////////////////////
void FitInterface::scheduleFitVarMatInit(const bool init)
{
initVarMat_ = init;
}
// register a data point ///////////////////////////////////////////////////////
void FitInterface::registerDataPoint(const Index k, const Index j)
{
checkYDim(j);
yDataIndex_[j][k] = true;
dataIndexSet_.insert(k);
scheduleLayoutInit();
}
// coordinate buffering ////////////////////////////////////////////////////////
void FitInterface::scheduleDataCoordInit(void)
{
initDataCoord_ = true;
scheduleFitVarMatInit();
}
void FitInterface::updateDataCoord(void) const
{
FitInterface * modThis = const_cast<FitInterface *>(this);
if (initDataCoord_)
{
modThis->dataCoord_.clear();
for (auto k: getDataIndexSet())
{
modThis->dataCoord_[k] = rowMajToCoord(k);
}
modThis->initDataCoord_ = false;
}
}
// global layout management ////////////////////////////////////////////////////
void FitInterface::scheduleLayoutInit(void)
{
initLayout_ = true;
scheduleFitVarMatInit();
}
bool FitInterface::initVarMat(void) const
{
return initVarMat_;
}
void FitInterface::updateLayout(void) const
{
if (initLayout_)
{
FitInterface * modThis = const_cast<FitInterface *>(this);
Layout & l = modThis->layout;
Index size, ifit;
vector<Index> v;
l.nXFitDim = 0;
l.nYFitDim = 0;
l.totalSize = 0;
l.totalXSize = 0;
l.totalYSize = 0;
l.xSize.clear();
l.ySize.clear();
l.dataIndexSet.clear();
l.xDim.clear();
l.yDim.clear();
l.xFitDim.clear();
l.yFitDim.clear();
l.x.clear();
l.y.clear();
l.xFit.clear();
l.yFit.clear();
ifit = 0;
for (Index i = 0; i < getNXDim(); ++i)
{
if (!xIsExact_[i])
{
l.nXFitDim++;
size = getXFitSize(i);
l.xSize.push_back(size);
l.totalXSize += size;
l.xDim.push_back(i);
l.xFitDim.push_back(layout.xDim.size() - 1);
l.x.push_back(vector<Index>());
l.xFit.push_back(vector<Index>());
for (Index r = 0; r < getXSize(i); ++r)
{
if (isXUsed(r, i))
{
l.x[ifit].push_back(r);
l.xFit[i].push_back(layout.x[ifit].size() - 1);
}
else
{
l.xFit[i].push_back(-1);
}
}
ifit++;
}
else
{
l.xFitDim.push_back(-1);
l.xFit.push_back(vector<Index>());
for (Index r = 0; r < getXSize(i); ++r)
{
l.xFit[i].push_back(-1);
}
}
}
for (Index j = 0; j < getNYDim(); ++j)
{
Index s = 0;
l.nYFitDim++;
size = getYFitSize(j);
l.ySize.push_back(size);
l.totalYSize += size;
l.yDim.push_back(j);
l.yFitDim.push_back(layout.yDim.size() - 1);
l.y.push_back(vector<Index>());
l.yFit.push_back(vector<Index>());
l.data.push_back(vector<Index>());
l.yFitFromData.push_back(map<Index, Index>());
for (auto &p: yDataIndex_[j])
{
if (p.second)
{
l.dataIndexSet.insert(p.first);
l.y[j].push_back(s);
l.yFit[j].push_back(layout.y[j].size() - 1);
l.data[j].push_back(p.first);
l.yFitFromData[j][p.first] = layout.y[j].size() - 1;
}
else
{
l.yFit[j].push_back(-1);
l.yFitFromData[j][p.first] = -1;
}
s++;
}
}
l.totalSize = layout.totalXSize + layout.totalYSize;
l.nXFitDim = static_cast<Index>(layout.xSize.size());
l.nYFitDim = static_cast<Index>(layout.ySize.size());
l.xIndFromData.resize(getMaxDataIndex());
for (Index k: layout.dataIndexSet)
{
v = dataCoord(k);
for (Index i = 0; i < getNXDim(); ++i)
{
l.xIndFromData[k].push_back(indX(v[i], i));
}
}
modThis->initLayout_ = false;
}
}
Index FitInterface::indX(const Index r, const Index i) const
{
Index ind = -1;
if (layout.xFit[i][r] != -1)
{
Index ifit = layout.xFitDim[i], rfit = layout.xFit[i][r];
ind = layout.totalYSize;
for (Index a = 0; a < ifit; ++a)
{
ind += layout.xSize[a];
}
ind += rfit;
}
return ind;
}
Index FitInterface::indY(const Index k, const Index j) const
{
Index ind = -1;
if (layout.yFitFromData[j].at(k) != -1)
{
Index jfit = layout.yFitDim[j], sfit = layout.yFitFromData[j].at(k);
ind = 0;
for (Index b = 0; b < jfit; ++b)
{
ind += layout.ySize[b];
}
ind += sfit;
}
return ind;
}
// function to convert an row-major index into coordinates /////////////////////
vector<Index> FitInterface::rowMajToCoord(const Index k) const
{
vector<Index> v(getNXDim());
Index buf, dimProd;
checkDataIndex(k);
buf = k;
dimProd = 1;
for (Index d = getNXDim() - 1; d >= 0; --d)
{
v[d] = (buf/dimProd)%xSize_[d];
buf -= dimProd*v[d];
dimProd *= xSize_[d];
}
return v;
}
// IO //////////////////////////////////////////////////////////////////////////
ostream & Latan::operator<<(ostream &out, FitInterface &f)
{
out << "X dimensions: " << f.getNXDim() << endl;
for (Index i = 0; i < f.getNXDim(); ++i)
{
out << " * " << i << " \"" << f.xName().getName(i) << "\": ";
out << f.getXSize(i) << " value(s)";
if (f.isXExact(i))
{
out << " (assumed exact)";
}
out << endl;
}
out << "Y dimensions: " << f.getNYDim() << endl;
for (Index j = 0; j < f.getNYDim(); ++j)
{
out << " * " << j << " \"" << f.yName().getName(j) << "\": ";
out << f.getYSize(j) << " value(s)" << endl;
for (auto &p: f.yDataIndex_[j])
{
out << " " << setw(3) << p.first << " (";
for (auto vi: f.dataCoord(p.first))
{
out << vi << ",";
}
out << "\b) fit: " << (p.second ? "true" : "false") << endl;
}
}
out << "X/X correlations (r1 i1 r2 i2): ";
if (f.xxCorr_.empty())
{
out << "no" << endl;
}
else
{
out << endl;
for (auto &c: f.xxCorr_)
{
out << " * ";
for (auto i: c)
{
out << i << " ";
}
out << endl;
}
}
out << "Y/Y correlations (k1 j1 k2 j2): ";
if (f.yyCorr_.empty())
{
out << "no" << endl;
}
else
{
out << endl;
for (auto &c: f.yyCorr_)
{
out << " * ";
for (auto i: c)
{
out << i << " ";
}
out << endl;
}
}
out << "X/Y correlations (r i k j): ";
if (f.xyCorr_.empty())
{
out << "no";
}
else
{
out << endl;
for (auto &c: f.xyCorr_)
{
out << " * ";
for (auto i: c)
{
out << i << " ";
}
out << endl;
}
}
return out;
}

231
lib/FitInterface.hpp
View File

@ -1,7 +1,7 @@
/*
* FitInterface.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -21,63 +21,208 @@
#define Latan_FitInterface_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <LatAnalyze/Mat.hpp>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* base class for data fit *
* FitInterface *
******************************************************************************/
class FitInterface
{
private:
typedef struct
{
Index nXFitDim, nYFitDim;
// X/Y block sizes
Index totalSize, totalXSize, totalYSize;
// size of each X/Y dimension
std::vector<Index> xSize, ySize;
// set of active data indices
std::set<Index> dataIndexSet;
// lookup tables
// xDim : x fit dim ifit -> x dim i
// x : x fit point ifit,rfit -> x point r
// xFitDim : x dim i -> x fit dim ifit (-1 if empty)
// xFit : x point i,r -> x fit point rfit (-1 if empty)
// data : y fit point jfit,sfit -> y point index k
// yFitFromData: y point index k,j -> y fit point sfit (-1 if empty)
// xIndFromData: data index k -> index of coordinates of associated x
std::vector<Index> xDim, yDim, xFitDim, yFitDim;
std::vector<std::vector<Index>> x, y, data, xFit, yFit;
std::vector<std::map<Index, Index>> yFitFromData;
// no map here for fit performance
std::vector<std::vector<Index>> xIndFromData;
} Layout;
public:
typedef Minimizer::Verbosity FitVerbosity;
public:
// constructors
FitInterface(void) = default;
FitInterface(const Index nData, const Index xDim, const Index yDim);
// constructor
FitInterface(void);
// destructor
virtual ~FitInterface(void) = default;
// copy object (not as a constructor to be accessed from derived class)
void copyInterface(const FitInterface &d);
// add dimensions
void addXDim(const Index nData, const std::string name = "",
const bool isExact = false);
void addYDim(const std::string name = "");
// access
void assumeXExact(const Index i, const bool isExact = true);
void assumeXXCorrelated(const Index i1, const Index i2,
const bool isCorrelated = true);
void assumeYYCorrelated(const Index j1, const Index j2,
const bool isCorrelated = true);
void assumeYXCorrelated(const Index j, const Index i,
const bool isCorrelated = true);
void assumeDataCorrelated(const Index k1, const Index k2,
const bool isCorrelated = true);
void assumeDataCorrelated(const bool isCorrelated = true);
void fitPoint(const Index k, const bool isFitPoint = true);
void fitPointRange(const Index k1, const Index k2,
const bool isFitPoint = true);
void fitAllPoints(const bool isFitPoint = true);
Index getNData(void) const;
Index getNFitPoint(void) const;
Index getXDim(void) const;
Index getYDim(void) const;
Index getStatXDim(void) const;
double getSvdTolerance(void) const;
void setFitInterface(const FitInterface &fitInterface);
void setNData(const Index nData);
void setXDim(const Index xDim);
void setYDim(const Index yDim);
void setSvdTolerance(const double tolerance);
void resize(const Index nData, const Index xDim, const Index yDim);
// test
bool isFitPoint(const Index k) const;
Index getNXDim(void) const;
Index getNYDim(void) const;
Index getXSize(void) const;
Index getXSize(const Index i) const;
Index getYSize(void) const;
Index getYSize(const Index j) const;
Index getXFitSize(void) const;
Index getXFitSize(const Index i) const;
Index getYFitSize(void) const;
Index getYFitSize(const Index j) const;
Index getMaxDataIndex(void) const;
const std::set<Index> & getDataIndexSet(void) const;
double getSvdTolerance(void) const;
void setSvdTolerance(const double &tol);
VarName & xName(void);
const VarName & xName(void) const;
VarName & yName(void);
const VarName & yName(void) const;
// Y dimension index helper
template <typename... Ts>
Index dataIndex(const Ts... is) const;
Index dataIndex(const std::vector<Index> &v) const;
const std::vector<Index> & dataCoord(const Index k) const;
// enable fit points
void fitPoint(const bool isFitPoint, const Index k, const Index j = 0);
// variance interface
void assumeXExact(const bool isExact, const Index i);
void assumeXXCorrelated(const bool isCorr, const Index r1, const Index i1,
const Index r2, const Index i2);
void assumeXXCorrelated(const bool isCorr, const Index i1, const Index i2);
void assumeYYCorrelated(const bool isCorr, const Index k1, const Index j1,
const Index k2, const Index j2);
void assumeYYCorrelated(const bool isCorr, const Index j1, const Index j2);
void assumeXYCorrelated(const bool isCorr, const Index r, const Index i,
const Index k, const Index j);
void assumeXYCorrelated(const bool isCorr, const Index i, const Index j);
// tests
bool pointExists(const Index k) const;
bool pointExists(const Index k, const Index j) const;
bool isXExact(const Index i) const;
bool isXXCorrelated(const Index i1, const Index i2) const;
bool isYYCorrelated(const Index j1, const Index j2) const;
bool isYXCorrelated(const Index j, const Index i) const;
bool isDataCorrelated(const Index k1, const Index k2) const;
bool isXUsed(const Index r, const Index i, const bool inFit = true) const;
bool isFitPoint(const Index k, const Index j) const;
bool isXXCorrelated(const Index r1, const Index i1, const Index r2,
const Index i2) const;
bool isYYCorrelated(const Index k1, const Index j1, const Index k2,
const Index j2) const;
bool isXYCorrelated(const Index r, const Index i, const Index k,
const Index j) const;
bool hasCorrelations(void) const;
// make correlation filter for fit variance matrix
DMat makeCorrFilter(void);
// schedule variance matrix initialization
void scheduleFitVarMatInit(const bool init = true);
// IO
friend std::ostream & operator<<(std::ostream &out, FitInterface &f);
protected:
// register a data point
void registerDataPoint(const Index k, const Index j = 0);
// add correlation to a set
static void addCorr(std::set<std::array<Index, 4>> &s, const bool isCorr,
const std::array<Index, 4> &c);
// abstract methods to create data containers
virtual void createXData(const std::string name, const Index nData) = 0;
virtual void createYData(const std::string name) = 0;
// coordinate buffering
void scheduleDataCoordInit(void);
void updateDataCoord(void) const;
// global layout management
void scheduleLayoutInit(void);
bool initVarMat(void) const;
void updateLayout(void) const;
Index indX(const Index r, const Index i) const;
Index indY(const Index k, const Index j) const;
private:
IVec isXExact_, isFitPoint_;
IMat isXXCorr_, isYYCorr_, isYXCorr_, isDataCorr_;
double svdTolerance_{1.0e-10};
// function to convert an row-major index into coordinates
std::vector<Index> rowMajToCoord(const Index k) const;
protected:
Layout layout;
private:
VarName xName_, yName_;
std::vector<Index> xSize_;
std::vector<bool> xIsExact_;
std::map<Index, std::vector<Index>> dataCoord_;
std::set<Index> dataIndexSet_;
std::vector<std::map<Index, bool>> yDataIndex_;
std::set<std::array<Index, 4>> xxCorr_, yyCorr_, xyCorr_;
Index maxDataIndex_{1};
bool initLayout_{true};
bool initVarMat_{true};
bool initDataCoord_{true};
double svdTol_{1.e-10};
};
std::ostream & operator<<(std::ostream &out, FitInterface &f);
/******************************************************************************
* FitInterface template implementation *
******************************************************************************/
// Y dimension index helper ////////////////////////////////////////////////////
template <typename... Ts>
Index FitInterface::dataIndex(const Ts... coords) const
{
static_assert(static_or<std::is_convertible<Index, Ts>::value...>::value,
"fitPoint arguments are not compatible with Index");
const std::vector<Index> coord = {coords...};
return dataIndex(coord);
}
/******************************************************************************
* error check macros *
******************************************************************************/
#define checkXDim(i)\
if ((i) >= getNXDim())\
{\
LATAN_ERROR(Range, "X dimension " + strFrom(i) + " out of range");\
}
#define checkXIndex(vi, i)\
if ((vi) >= getXSize(i))\
{\
LATAN_ERROR(Range, "index " + strFrom(vi) + " in X dimension "\
+ strFrom(i) + " out of range");\
}
#define checkYDim(j)\
if ((j) >= getNYDim())\
{\
LATAN_ERROR(Range, "Y dimension " + strFrom(j) + " out of range");\
}
#define checkDataIndex(k)\
if ((k) >= getMaxDataIndex())\
{\
LATAN_ERROR(Range, "data point index " + strFrom(k) + " invalid");\
}
#define checkDataCoord(v)\
if (static_cast<Index>((v).size()) != getNXDim())\
{\
LATAN_ERROR(Size, "number of coordinates and number of X dimensions "\
"mismatch");\
}\
for (unsigned int i_ = 0; i_ < (v).size(); ++i_)\
{\
checkXIndex((v)[i_], i_);\
}
#define checkPoint(k, j)\
if (!pointExists(k, j))\
{\
LATAN_ERROR(Range, "no data point in Y dimension " + strFrom(j)\
+ " with index " + strFrom(k));\
}
END_LATAN_NAMESPACE
#endif // Latan_FitInterface_hpp_

13
lib/Function.cpp
View File

@ -1,7 +1,7 @@
/*
* Function.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -29,6 +29,7 @@ using namespace Latan;
// constructor /////////////////////////////////////////////////////////////////
DoubleFunction::DoubleFunction(const vecFunc &f, const Index nArg)
: buffer_(new DVec)
, varName_("x")
{
setFunction(f, nArg);
}
@ -45,6 +46,16 @@ void DoubleFunction::setFunction(const vecFunc &f, const Index nArg)
f_ = f;
}
VarName & DoubleFunction::varName(void)
{
return varName_;
}
const VarName & DoubleFunction::varName(void) const
{
return varName_;
}
// error checking //////////////////////////////////////////////////////////////
void DoubleFunction::checkSize(const Index nPar) const
{

12
lib/Function.hpp
View File

@ -1,7 +1,7 @@
/*
* Function.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,9 +23,6 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <functional>
#include <stack>
#include <vector>
BEGIN_LATAN_NAMESPACE
@ -43,8 +40,10 @@ public:
// destructor
virtual ~DoubleFunction(void) = default;
// access
virtual Index getNArg(void) const;
void setFunction(const vecFunc &f, const Index nArg);
virtual Index getNArg(void) const;
void setFunction(const vecFunc &f, const Index nArg);
VarName & varName(void);
const VarName & varName(void) const;
// function call
double operator()(const double *arg) const;
double operator()(const DVec &arg) const;
@ -75,6 +74,7 @@ private:
void checkSize(const Index nPar) const;
private:
std::shared_ptr<DVec> buffer_{nullptr};
VarName varName_;
vecFunc f_;
};

4
lib/Global.cpp
View File

@ -1,7 +1,7 @@
/*
* Global.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,6 +23,8 @@
using namespace std;
using namespace Latan;
PlaceHolder Latan::_;
const string Env::fullName = PACKAGE_STRING;
const string Env::name = PACKAGE_NAME;
const string Env::version = PACKAGE_VERSION;

2
lib/Global.hpp
View File

@ -1,7 +1,7 @@
/*
* Global.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/GslHybridRootFinder.cpp
View File

@ -1,7 +1,7 @@
/*
* GslHybridRootFinder.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/GslHybridRootFinder.hpp
View File

@ -1,7 +1,7 @@
/*
* GslHybridRootFinder.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

361
lib/GslMinimizer.cpp Normal file
View File

@ -0,0 +1,361 @@
/*
* GslMinimizer.cpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2016 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/GslMinimizer.hpp>
#include <LatAnalyze/includes.hpp>
#include <LatAnalyze/Math.hpp>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_blas.h>
using namespace std;
using namespace Latan;
/******************************************************************************
* GslMinimizer implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
GslMinimizer::GslMinimizer(const Algorithm algorithm)
{
setAlgorithm(algorithm);
der_.setOrder(1, 1);
}
// access //////////////////////////////////////////////////////////////////////
GslMinimizer::Algorithm GslMinimizer::getAlgorithm(void) const
{
return algorithm_;
}
void GslMinimizer::setAlgorithm(const Algorithm algorithm)
{
algorithm_ = algorithm;
}
bool GslMinimizer::supportLimits(void) const
{
return false;
}
// test ////////////////////////////////////////////////////////////////////////
bool GslMinimizer::isDerAlgorithm(const Algorithm algorithm)
{
return (algorithm <= Algorithm::lastDerAlg);
}
// minimization ////////////////////////////////////////////////////////////////
const DVec & GslMinimizer::operator()(const DoubleFunction &f)
{
DVec &x = getState();
// resize minimizer state to match function number of arguments
if (f.getNArg() != x.size())
{
resize(f.getNArg());
}
// set function data
GslFuncData data;
der_.setFunction(f);
data.f = &f;
data.d = &der_;
// set initial position
gsl_vector *gslX = gsl_vector_alloc(getDim());
for (Index i = 0; i < getDim(); ++i)
{
gsl_vector_set(gslX, i, x(i));
}
// minimization
int status;
if (isDerAlgorithm(getAlgorithm()))
{
// set function
gsl_multimin_function_fdf gslFunc;
gslFunc.n = getDim();
gslFunc.f = &fWrapper;
gslFunc.df = &dfWrapper;
gslFunc.fdf = &fdfWrapper;
gslFunc.params = &data;
// create and set minimizer
const gsl_multimin_fdfminimizer_type *gslAlg;
gsl_multimin_fdfminimizer *gslMin;
switch (getAlgorithm())
{
case Algorithm::cgFR:
gslAlg = gsl_multimin_fdfminimizer_conjugate_fr;
break;
case Algorithm::cgPR:
gslAlg = gsl_multimin_fdfminimizer_conjugate_pr;
break;
case Algorithm::bfgs:
gslAlg = gsl_multimin_fdfminimizer_vector_bfgs;
break;
case Algorithm::bfgs2:
gslAlg = gsl_multimin_fdfminimizer_vector_bfgs2;
break;
case Algorithm::steepDesc:
gslAlg = gsl_multimin_fdfminimizer_vector_bfgs2;
break;
default:
LATAN_ERROR(Argument, "unknow GSL minization algorithm "
+ strFrom(static_cast<int>(getAlgorithm())));
break;
}
gslMin = gsl_multimin_fdfminimizer_alloc(gslAlg, getDim());
// minimize
unsigned int pass = 0, it;
double dxRel;
do
{
pass++;
gsl_multimin_fdfminimizer_set(gslMin, &gslFunc, gslX, 0.01, 0.001);
if (getVerbosity() >= Verbosity::Normal)
{
cout << "========== GSL minimization, pass #" << pass;
cout << " ==========" << endl;
cout << "Algorithm: " << getAlgorithmName(getAlgorithm());
cout << endl;
cout << "Max eval.= " << getMaxIteration();
cout << " -- Precision= " << getPrecision() << endl;
printf("Starting f(x)= %.10e\n", f(x));
}
it = 0;
do
{
it++;
gsl_multimin_fdfminimizer_iterate(gslMin);
dxRel = gsl_blas_dnrm2(gslMin->dx)/gsl_blas_dnrm2(gslMin->x);
status = (dxRel < getPrecision()) ? GSL_SUCCESS : GSL_CONTINUE;
if (getVerbosity() >= Verbosity::Debug)
{
printf("iteration %4d: f= %.10e dxrel= %.10e eval= %d\n",
it, gslMin->f, dxRel, data.evalCount);
}
} while (status == GSL_CONTINUE and
(data.evalCount < getMaxIteration()));
if (getVerbosity() >= Verbosity::Normal)
{
printf("Found minimum %.10e at:\n", gslMin->f);
for (Index i = 0; i < x.size(); ++i)
{
printf("%8s= %.10e\n", f.varName().getName(i).c_str(),
gsl_vector_get(gslMin->x, i));
}
cout << "after " << data.evalCount << " evaluations" << endl;
cout << "Minimization ended with code " << status;
cout << endl;
}
data.evalCount = 0;
for (Index i = 0; i < getDim(); ++i)
{
gsl_vector_set(gslX, i, gsl_vector_get(gslMin->x, i));
}
} while (status != GSL_SUCCESS and (pass < getMaxPass()));
// deallocate GSL minimizer
gsl_multimin_fdfminimizer_free(gslMin);
}
else
{
// set function
gsl_multimin_function gslFunc;
gslFunc.n = getDim();
gslFunc.f = &fWrapper;
gslFunc.params = &data;
// create and set minimizer
const gsl_multimin_fminimizer_type *gslAlg;
gsl_multimin_fminimizer *gslMin;
switch (getAlgorithm())
{
case Algorithm::simplex:
gslAlg = gsl_multimin_fminimizer_nmsimplex;
break;
case Algorithm::simplex2:
gslAlg = gsl_multimin_fminimizer_nmsimplex2;
break;
case Algorithm::simplex2R:
gslAlg = gsl_multimin_fminimizer_nmsimplex2rand;
break;
default:
LATAN_ERROR(Argument, "unknow GSL minization algorithm "
+ strFrom(static_cast<int>(getAlgorithm())));
break;
}
gslMin = gsl_multimin_fminimizer_alloc(gslAlg, getDim());
// minimize
unsigned int pass = 0, it;
gsl_vector *step = gsl_vector_alloc(getDim());
double relSize;
gsl_vector_set_all(step, 0.01);
do
{
pass++;
gsl_multimin_fminimizer_set(gslMin, &gslFunc, gslX, step);
if (getVerbosity() >= Verbosity::Normal)
{
cout << "========== GSL minimization, pass #" << pass;
cout << " ==========" << endl;
cout << "Algorithm: " << getAlgorithmName(getAlgorithm());
cout << endl;
cout << "Max eval.= " << getMaxIteration();
cout << " -- Precision= " << getPrecision() << endl;
printf("Starting f(x)= %.10e\n", f(x));
}
it = 0;
do
{
it++;
gsl_multimin_fminimizer_iterate(gslMin);
relSize = Math::pow<2>(gslMin->size)/gsl_blas_dnrm2(gslMin->x);
status = (relSize < getPrecision()) ? GSL_SUCCESS
: GSL_CONTINUE;
if (getVerbosity() >= Verbosity::Debug)
{
printf("iteration %4d: f= %.10e relSize= %.10e eval= %d\n",
it, gslMin->fval, relSize, data.evalCount);
}
} while (status == GSL_CONTINUE and
(data.evalCount < getMaxIteration()));
if (getVerbosity() >= Verbosity::Normal)
{
printf("Found minimum %.10e at:\n", gslMin->fval);
for (Index i = 0; i < x.size(); ++i)
{
printf("%8s= %.10e\n", f.varName().getName(i).c_str(),
gsl_vector_get(gslMin->x, i));
}
cout << "after " << data.evalCount << " evaluations" << endl;
cout << "Minimization ended with code " << status;
cout << endl;
}
data.evalCount = 0;
for (Index i = 0; i < getDim(); ++i)
{
gsl_vector_set(gslX, i, gsl_vector_get(gslMin->x, i));
}
} while (status != GSL_SUCCESS and (pass < getMaxPass()));
// deallocate GSL minimizer
gsl_multimin_fminimizer_free(gslMin);
gsl_vector_free(step);
}
if (status != GSL_SUCCESS)
{
LATAN_WARNING("invalid minimum: maximum number of call reached");
}
// save final result
for (Index i = 0; i < getDim(); ++i)
{
x(i) = gsl_vector_get(gslX, i);
}
// deallocate GSL state and return
gsl_vector_free(gslX);
return x;
}
// function wrappers ///////////////////////////////////////////////////////////
double GslMinimizer::fWrapper(const gsl_vector *x, void *vdata)
{
GslFuncData &data = *static_cast<GslFuncData *>(vdata);
data.evalCount++;
return (*data.f)(x->data);
}
void GslMinimizer::dfWrapper(const gsl_vector *x, void *vdata, gsl_vector * df)
{
GslFuncData &data = *static_cast<GslFuncData *>(vdata);
const unsigned int n = data.f->getNArg();
for (unsigned int i = 0; i < n; ++i)
{
data.d->setDir(i);
gsl_vector_set(df, i, (*(data.d))(x->data));
}
data.evalCount += data.d->getNPoint()*n;
}
void GslMinimizer::fdfWrapper(const gsl_vector *x, void *vdata, double *f,
gsl_vector * df)
{
GslFuncData &data = *static_cast<GslFuncData *>(vdata);
const unsigned int n = data.f->getNArg();
for (unsigned int i = 0; i < n; ++i)
{
data.d->setDir(i);
gsl_vector_set(df, i, (*(data.d))(x->data));
}
*f = (*data.f)(x->data);
data.evalCount += data.d->getNPoint()*n + 1;
}
// algorithm names /////////////////////////////////////////////////////////////
string GslMinimizer::getAlgorithmName(const Algorithm algorithm)
{
switch (algorithm)
{
case Algorithm::cgFR:
return "Fletcher-Reeves conjugate gradient";
break;
case Algorithm::cgPR:
return "Polak-Ribiere conjugate gradient";
break;
case Algorithm::bfgs:
return "Broyden-Fletcher-Goldfarb-Shanno";
break;
case Algorithm::bfgs2:
return "improved Broyden-Fletcher-Goldfarb-Shanno";
break;
case Algorithm::steepDesc:
return "steepest descent";
break;
case Algorithm::simplex:
return "Nelder-Mead simplex";
break;
case Algorithm::simplex2:
return "improved Nelder-Mead simplex";
break;
case Algorithm::simplex2R:
return "improved Nelder-Mead simplex with random start";
break;
}
return "";
}

86
lib/GslMinimizer.hpp Normal file
View File

@ -0,0 +1,86 @@
/*
* GslMinimizer.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2016 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_GslMinimizer_hpp_
#define Latan_GslMinimizer_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Derivative.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <gsl/gsl_vector.h>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* interface to the GSL minimizers *
******************************************************************************/
class GslMinimizer: public Minimizer
{
public:
enum class Algorithm
{
cgFR = 1,
cgPR = 2,
bfgs = 3,
bfgs2 = 4,
steepDesc = 5,
lastDerAlg = 5,
simplex = 6,
simplex2 = 7,
simplex2R = 8
};
private:
struct GslFuncData
{
const DoubleFunction *f{nullptr};
Derivative *d{nullptr};
unsigned int evalCount{0};
};
public:
// constructor
explicit GslMinimizer(const Algorithm algorithm = defaultAlg_);
// destructor
virtual ~GslMinimizer(void) = default;
// access
Algorithm getAlgorithm(void) const;
void setAlgorithm(const Algorithm algorithm);
virtual bool supportLimits(void) const;
// minimization
virtual const DVec & operator()(const DoubleFunction &f);
private:
// test
static bool isDerAlgorithm(const Algorithm algorithm);
// function wrappers
static double fWrapper(const gsl_vector *x, void * params);
static void dfWrapper(const gsl_vector *x, void * params,
gsl_vector * df);
static void fdfWrapper(const gsl_vector *x, void *params, double *f,
gsl_vector * df);
// algorithm names
std::string getAlgorithmName(const Algorithm algorithm);
private:
Algorithm algorithm_;
static constexpr Algorithm defaultAlg_ = Algorithm::simplex2;
CentralDerivative der_;
};
END_LATAN_NAMESPACE
#endif // Latan_GslMinimizer_hpp_

2
lib/GslQagsIntegrator.cpp
View File

@ -1,7 +1,7 @@
/*
* GslQagsIntegrator.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/GslQagsIntegrator.hpp
View File

@ -1,7 +1,7 @@
/*
* GslQagsIntegrator.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

91
lib/Hdf5File.cpp
View File

@ -1,7 +1,7 @@
/*
* Hdf5File.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli, Matt Spraggs
* Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@ -30,8 +30,8 @@ using namespace H5NS;
constexpr unsigned int maxGroupNameSize = 1024u;
const short dMatType = static_cast<short>(IoObject::IoType::dMat);
const short dSampleType = static_cast<short>(IoObject::IoType::dSample);
const short dMatSampleType = static_cast<short>(IoObject::IoType::dMatSample);
const short rgStateType = static_cast<short>(IoObject::IoType::rgState);
/******************************************************************************
* Hdf5File implementation *
@ -65,7 +65,6 @@ void Hdf5File::save(const DMat &m, const string &name)
hsize_t dim[2] = {static_cast<hsize_t>(m.rows()),
static_cast<hsize_t>(m.cols())};
hsize_t attrDim = 1;
DataSpace dataSpace(2, dim), attrSpace(1, &attrDim);
group = h5File_->createGroup(name.c_str() + nameOffset(name));
@ -75,7 +74,31 @@ void Hdf5File::save(const DMat &m, const string &name)
dataset.write(m.data(), PredType::NATIVE_DOUBLE);
}
void Hdf5File::save(const DMatSample &sample, const string &name)
void Hdf5File::save(const DSample &ds, const string &name)
{
if (name.empty())
{
LATAN_ERROR(Io, "trying to save data with an empty name");
}
Group group;
Attribute attr;
DataSet dataset;
hsize_t dim = static_cast<hsize_t>(ds.size() + 1);
hsize_t attrDim = 1;
DataSpace dataSpace(1, &dim), attrSpace(1, &attrDim);
const long int nSample = ds.size();
group = h5File_->createGroup(name.c_str() + nameOffset(name));
attr = group.createAttribute("type", PredType::NATIVE_SHORT, attrSpace);
attr.write(PredType::NATIVE_SHORT, &dSampleType);
attr = group.createAttribute("nSample", PredType::NATIVE_LONG, attrSpace);
attr.write(PredType::NATIVE_LONG, &nSample);
dataset = group.createDataSet("data", PredType::NATIVE_DOUBLE, dataSpace);
dataset.write(ds.data(), PredType::NATIVE_DOUBLE);
}
void Hdf5File::save(const DMatSample &ms, const string &name)
{
if (name.empty())
{
@ -85,11 +108,11 @@ void Hdf5File::save(const DMatSample &sample, const string &name)
Group group;
Attribute attr;
DataSet dataset;
hsize_t dim[2] = {static_cast<hsize_t>(sample[central].rows()),
static_cast<hsize_t>(sample[central].cols())};
hsize_t dim[2] = {static_cast<hsize_t>(ms[central].rows()),
static_cast<hsize_t>(ms[central].cols())};
hsize_t attrDim = 1;
DataSpace dataSpace(2, dim), attrSpace(1, &attrDim);
const long int nSample = sample.size();
const long int nSample = ms.size();
string datasetName;
group = h5File_->createGroup(name.c_str() + nameOffset(name));
@ -97,37 +120,16 @@ void Hdf5File::save(const DMatSample &sample, const string &name)
attr.write(PredType::NATIVE_SHORT, &dMatSampleType);
attr = group.createAttribute("nSample", PredType::NATIVE_LONG, attrSpace);
attr.write(PredType::NATIVE_LONG, &nSample);
FOR_STAT_ARRAY(sample, s)
FOR_STAT_ARRAY(ms, s)
{
datasetName = (s == central) ? "data_C" : ("data_S_" + strFrom(s));
dataset = group.createDataSet(datasetName.c_str(),
PredType::NATIVE_DOUBLE,
dataSpace);
dataset.write(sample[s].data(), PredType::NATIVE_DOUBLE);
dataset.write(ms[s].data(), PredType::NATIVE_DOUBLE);
}
}
void Hdf5File::save(const RandGenState &state, const string &name)
{
if (name.empty())
{
LATAN_ERROR(Io, "trying to save data with an empty name");
}
Group group;
Attribute attr;
DataSet dataset;
hsize_t dim = RLXG_STATE_SIZE;
hsize_t attrDim = 1;
DataSpace dataSpace(1, &dim), attrSpace(1, &attrDim);
group = h5File_->createGroup(name.c_str() + nameOffset(name));
attr = group.createAttribute("type", PredType::NATIVE_SHORT, attrSpace);
attr.write(PredType::NATIVE_SHORT, &rgStateType);
dataset = group.createDataSet("data", PredType::NATIVE_INT, dataSpace);
dataset.write(state.data(), PredType::NATIVE_INT);
}
// read first name ////////////////////////////////////////////////////////////
string Hdf5File::getFirstName(void)
{
@ -253,18 +255,15 @@ void Hdf5File::load(DMat &m, const DataSet &d)
d.read(m.data(), PredType::NATIVE_DOUBLE);
}
void Hdf5File::load(RandGenState &state, const DataSet &d)
void Hdf5File::load(DSample &ds, const DataSet &d)
{
DataSpace dataspace;
hsize_t dim[1];
dataspace = d.getSpace();
dataspace.getSimpleExtentDims(dim);
if (dim[0] != RLXG_STATE_SIZE)
{
LATAN_ERROR(Io, "random generator state has a wrong length");
}
d.read(state.data(), PredType::NATIVE_INT);
ds.resize(dim[0] - 1);
d.read(ds.data(), PredType::NATIVE_DOUBLE);
}
string Hdf5File::load(const string &name)
@ -296,6 +295,15 @@ string Hdf5File::load(const string &name)
load(*pt, dataset);
break;
}
case IoObject::IoType::dSample:
{
DSample *pt = new DSample;
data_[groupName].reset(pt);
dataset = group.openDataSet("data");
load(*pt, dataset);
break;
}
case IoObject::IoType::dMatSample:
{
DMatSample *pt = new DMatSample;
@ -320,15 +328,6 @@ string Hdf5File::load(const string &name)
}
break;
}
case IoObject::IoType::rgState:
{
RandGenState *pt = new RandGenState;
data_[groupName].reset(pt);
dataset = group.openDataSet("data");
load(*pt, dataset);
break;
}
default:
{
LATAN_ERROR(Io, "unknown data type ("

9
lib/Hdf5File.hpp
View File

@ -1,7 +1,7 @@
/*
* Hdf5File.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli, Matt Spraggs
* Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@ -24,7 +24,6 @@
#include <LatAnalyze/File.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <H5Cpp.h>
BEGIN_LATAN_NAMESPACE
@ -46,8 +45,8 @@ public:
virtual ~Hdf5File(void);
// access
virtual void save(const DMat &m, const std::string &name);
virtual void save(const DMatSample &s, const std::string &name);
virtual void save(const RandGenState &state, const std::string &name);
virtual void save(const DSample &ds, const std::string &name);
virtual void save(const DMatSample &ms, const std::string &name);
// read first name
virtual std::string getFirstName(void);
// tests
@ -60,8 +59,8 @@ private:
std::string getFirstGroupName(void);
virtual std::string load(const std::string &name = "");
void load(DMat &m, const H5NS::DataSet &d);
void load(DSample &ds, const H5NS::DataSet &d);
void load(DMatSample &s, const H5NS::DataSet &d);
void load(RandGenState &state, const H5NS::DataSet &d);
// check name for forbidden characters
static size_t nameOffset(const std::string &name);
private:

2
lib/Histogram.cpp
View File

@ -1,7 +1,7 @@
/*
* Histogram.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/Histogram.hpp
View File

@ -1,7 +1,7 @@
/*
* Histogram.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/Integrator.hpp
View File

@ -1,7 +1,7 @@
/*
* Integrator.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

18
lib/Io.cpp
View File

@ -1,7 +1,7 @@
/*
* Io.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -19,6 +19,10 @@
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/includes.hpp>
#include <LatAnalyze/AsciiFile.hpp>
#ifdef HAVE_HDF5
#include <LatAnalyze/Hdf5File.hpp>
#endif
using namespace std;
using namespace Latan;
@ -34,14 +38,16 @@ unique_ptr<File> Io::open(const std::string &fileName, const unsigned int mode)
{
string ext = extension(fileName);
if (ext == "h5")
{
return unique_ptr<File>(new Hdf5File(fileName, mode));
}
else if ((ext == "dat")||(ext == "sample")||(ext == "seed"))
if ((ext == "dat")||(ext == "sample")||(ext == "seed"))
{
return unique_ptr<File>(new AsciiFile(fileName, mode));
}
#ifdef HAVE_HDF5
else if (ext == "h5")
{
return unique_ptr<File>(new Hdf5File(fileName, mode));
}
#endif
else
{
LATAN_ERROR(Io, "unknown file extension '" + ext + "'");

5
lib/Io.hpp
View File

@ -1,7 +1,7 @@
/*
* Io.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -21,8 +21,7 @@
#define Latan_Io_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/AsciiFile.hpp>
#include <LatAnalyze/Hdf5File.hpp>
#include <LatAnalyze/File.hpp>
BEGIN_LATAN_NAMESPACE

5
lib/IoObject.hpp
View File

@ -1,7 +1,7 @@
/*
* IoObject.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -30,12 +30,13 @@ BEGIN_LATAN_NAMESPACE
class IoObject
{
public:
// conserve order for datafile retro-compatibility!
enum class IoType: short int
{
noType = 0,
dMat = 1,
dMatSample = 2,
rgState = 3
dSample = 3
};
public:
// destructor

57
lib/Makefile.am
View File

@ -1,24 +1,28 @@
COM_CXXFLAGS = -Wall -Wconversion
COM_CXXFLAGS = -Wall
if CXX_GNU
COM_CXXFLAGS += -W -pedantic
COM_CXXFLAGS += -W -pedantic -Wno-deprecated-declarations
else
if CXX_INTEL
COM_CXXFLAGS +=
COM_CXXFLAGS += -wd1682
endif
endif
AM_LFLAGS = -olex.yy.c
AM_YFLAGS = -d
BUILT_SOURCES = AsciiParser.hpp MathParser.hpp
lib_LTLIBRARIES = libLatAnalyze.la
noinst_LTLIBRARIES = libLexers.la
libLexers_la_SOURCES = AsciiLexer.lpp MathLexer.lpp
if CXX_GNU
libLexers_la_CXXFLAGS = $(COM_CXXFLAGS) -Wno-unused-parameter -Wno-unused-function -Wno-deprecated-register
else
libLexers_la_CXXFLAGS = $(COM_CXXFLAGS)
endif
libLatAnalyze_la_SOURCES = \
AsciiFile.cpp \
AsciiParser.ypp \
AsciiLexer.lpp \
Chi2Function.cpp \
CompiledFunction.cpp \
CompiledModel.cpp \
Derivative.cpp \
@ -28,8 +32,8 @@ libLatAnalyze_la_SOURCES = \
Function.cpp \
Global.cpp \
GslHybridRootFinder.cpp\
GslMinimizer.cpp \
GslQagsIntegrator.cpp \
Hdf5File.cpp \
Histogram.cpp \
includes.hpp \
Io.cpp \
@ -37,14 +41,12 @@ libLatAnalyze_la_SOURCES = \
Math.cpp \
MathInterpreter.cpp \
MathParser.ypp \
MathLexer.lpp \
MatSample.cpp \
Minimizer.cpp \
Model.cpp \
Plot.cpp \
RandGen.cpp \
RootFinder.cpp \
Solver.cpp \
StatArray.cpp \
TabFunction.cpp \
XYSampleData.cpp \
XYStatData.cpp \
@ -52,19 +54,18 @@ libLatAnalyze_la_SOURCES = \
libLatAnalyze_ladir = $(pkgincludedir)
libLatAnalyze_la_HEADERS = \
AsciiFile.hpp \
Chi2Function.hpp \
CompiledFunction.hpp \
CompiledModel.hpp \
Dataset.hpp \
Derivative.hpp \
Exceptions.hpp \
FitInterface.hpp \
Function.hpp \
File.hpp \
FitInterface.hpp \
Global.hpp \
GslHybridRootFinder.hpp\
GslMinimizer.hpp \
GslQagsIntegrator.hpp \
Hdf5File.hpp \
Histogram.hpp \
Integrator.hpp \
Io.hpp \
@ -77,17 +78,43 @@ libLatAnalyze_la_HEADERS = \
Model.hpp \
ParserState.hpp \
Plot.hpp \
RandGen.hpp \
RootFinder.hpp \
TabFunction.hpp \
Solver.hpp \
StatArray.hpp \
XYSampleData.hpp \
XYStatData.hpp
if HAVE_HDF5
libLatAnalyze_la_SOURCES += Hdf5File.cpp
libLatAnalyze_la_HEADERS += Hdf5File.hpp
endif
if HAVE_MINUIT
libLatAnalyze_la_SOURCES += MinuitMinimizer.cpp
libLatAnalyze_la_HEADERS += MinuitMinimizer.hpp
endif
if HAVE_NLOPT
libLatAnalyze_la_SOURCES += NloptMinimizer.cpp
libLatAnalyze_la_HEADERS += NloptMinimizer.hpp
endif
libLatAnalyze_la_CXXFLAGS = $(COM_CXXFLAGS)
libLatAnalyze_la_LIBADD = libLexers.la
if HAVE_AM_MINOR_LE_11
AsciiParser.hpp: AsciiParser.ypp
$(AM_V_YACC) $(YACC) -o AsciiParser.cpp --defines=AsciiParser.hpp $<
MathParser.hpp: MathParser.ypp
$(AM_V_YACC) $(YACC) -o MathParser.cpp --defines=MathParser.hpp $<
endif
BUILT_SOURCES = AsciiParser.hpp MathParser.hpp
CLEANFILES = \
MathLexer.cpp \
AsciiLexer.cpp \
AsciiParser.cpp\
AsciiParser.hpp\
MathParser.cpp \
MathParser.hpp
ACLOCAL_AMFLAGS = -I .buildutils/m4

2
lib/Mat.cpp
View File

@ -1,7 +1,7 @@
/*
* Mat.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/Mat.hpp
View File

@ -1,7 +1,7 @@
/*
* Mat.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

14
lib/MatSample.hpp
View File

@ -1,7 +1,7 @@
/*
* MatSample.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,7 +23,6 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/StatArray.hpp>
#include <functional>
BEGIN_LATAN_NAMESPACE
@ -35,7 +34,7 @@ BEGIN_LATAN_NAMESPACE
std::placeholders::_1, x))
template <typename T>
class MatSample: public Sample<Mat<T>>, public IoObject
class MatSample: public Sample<Mat<T>>
{
public:
// block type template
@ -104,9 +103,6 @@ public:
const Index nCol);
// resize all matrices
void resizeMat(const Index nRow, const Index nCol);
// IO type
virtual IoType getType(void) const;
};
// non-member operators
@ -383,12 +379,6 @@ void MatSample<T>::resizeMat(const Index nRow, const Index nCol)
}
}
// IO type /////////////////////////////////////////////////////////////////////
template <typename T>
IoObject::IoType MatSample<T>::getType(void) const
{
return IoType::noType;
}
END_LATAN_NAMESPACE

15
lib/Math.cpp
View File

@ -1,7 +1,7 @@
/*
* Math.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -24,6 +24,19 @@
using namespace std;
using namespace Latan;
/******************************************************************************
* Custom math functions *
******************************************************************************/
DMat MATH_NAMESPACE::varToCorr(const DMat &var)
{
DMat res = var, invDiag = res.diagonal();
invDiag = invDiag.cwiseInverse().cwiseSqrt();
res = (invDiag*invDiag.transpose()).cwiseProduct(res);
return res;
}
/******************************************************************************
* Standard C functions *
******************************************************************************/

5
lib/Math.hpp
View File

@ -1,7 +1,7 @@
/*
* Math.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -68,6 +68,9 @@ namespace MATH_NAMESPACE
return res;
}
// convert variance matrix to correlation matrix
DMat varToCorr(const DMat &var);
// Constants
const double pi = 3.1415926535897932384626433832795028841970;
const double e = 2.7182818284590452353602874713526624977572;

8
lib/MathInterpreter.cpp
View File

@ -1,7 +1,7 @@
/*
* MathInterpreter.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -411,7 +411,7 @@ void ExprNode::pushArg(ExprNode *node)
// ExprNode operators //////////////////////////////////////////////////////////
const ExprNode &ExprNode::operator[](const Index i) const
{
return *arg_[static_cast<unsigned int>(i)];
return *arg_[i];
}
ostream &Latan::operator<<(ostream &out, const ExprNode &n)
@ -451,7 +451,7 @@ try\
}\
catch (out_of_range)\
{\
address = static_cast<unsigned int>((table).size());\
address = (table).size();\
(table)[(name)] = address;\
}\
@ -582,7 +582,7 @@ MathInterpreter::MathInterpreter(const std::string &code)
// access //////////////////////////////////////////////////////////////////////
const Instruction * MathInterpreter::operator[](const Index i) const
{
return program_[static_cast<unsigned int>(i)].get();
return program_[i].get();
}
const ExprNode * MathInterpreter::getAST(void) const

7
lib/MathInterpreter.hpp
View File

@ -1,7 +1,7 @@
/*
* MathInterpreter.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -20,11 +20,6 @@
#ifndef Latan_MathInterpreter_hpp_
#define Latan_MathInterpreter_hpp_
#include <iostream>
#include <map>
#include <queue>
#include <string>
#include <stack>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/ParserState.hpp>

15
lib/MathLexer.lpp
View File

@ -1,7 +1,7 @@
/*
* MathLexer.lpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -25,20 +25,9 @@
%option yylineno
%{
#include <iostream>
#include <LatAnalyze/MathInterpreter.hpp>
#include "MathParser.hpp"
#if (defined __INTEL_COMPILER)
#pragma warning disable 1682 2259
#elif (defined __GNUC__)||(defined __clang__)
#pragma GCC diagnostic ignored "-Wconversion"
#pragma GCC diagnostic ignored "-Wsign-conversion"
#pragma GCC diagnostic ignored "-Wsign-compare"
#pragma GCC diagnostic ignored "-Wunused-function"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#endif
using namespace std;
using namespace Latan;

14
lib/MathParser.ypp
View File

@ -1,7 +1,7 @@
/*
* MathParser.ypp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -18,25 +18,15 @@
*/
%{
#include <iostream>
#include <sstream>
#include <cstring>
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/MathInterpreter.hpp>
#if (defined __INTEL_COMPILER)
#pragma warning disable 1682 2259
#elif (defined __GNUC__)||(defined __clang__)
#pragma GCC diagnostic ignored "-Wconversion"
#pragma GCC diagnostic ignored "-Wsign-conversion"
#endif
using namespace std;
using namespace Latan;
%}
%pure-parser
%name-prefix = "_math_"
%name-prefix "_math_"
%locations
%defines
%error-verbose

111
lib/Minimizer.cpp
View File

@ -1,7 +1,7 @@
/*
* Minimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,12 +23,6 @@
using namespace std;
using namespace Latan;
// constructor /////////////////////////////////////////////////////////////////
Minimizer::Minimizer(const Index dim)
{
resize(dim);
}
// access //////////////////////////////////////////////////////////////////////
void Minimizer::resize(const Index dim)
{
@ -49,119 +43,152 @@ void Minimizer::resize(const Index dim)
}
}
// limits //////////////////////////////////////////////////////////////////////
#define checkSupport \
if (!supportLimits())\
{\
LATAN_ERROR(Implementation, "minimizer does not support limits");\
}
double Minimizer::getHighLimit(const Index i) const
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid variable index");
}
return highLimit_(i);
}
const DVec & Minimizer::getHighLimit(const PlaceHolder ph __dumb) const
{
checkSupport;
return highLimit_;
}
double Minimizer::getLowLimit(const Index i) const
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid variable index");
}
return lowLimit_(i);
}
const DVec & Minimizer::getLowLimit(const PlaceHolder ph __dumb) const
{
checkSupport;
return lowLimit_;
}
bool Minimizer::hasHighLimit(const Index i) const
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid variable index");
}
return hasHighLimit_(i);
}
bool Minimizer::hasLowLimit(const Index i) const
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid variable index");
}
return hasLowLimit_(i);
}
void Minimizer::setHighLimit(const Index i, const double l)
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid limit index");
}
else
{
highLimit_(i) = l;
useHighLimit(i);
resize(i + 1);
}
highLimit_(i) = l;
useHighLimit(i);
}
void Minimizer::setHighLimit(const PlaceHolder ph __dumb, const DVec &l)
{
checkSupport;
if (l.size() != getDim())
{
LATAN_ERROR(Size, "invalid limit vector size");
}
else
{
highLimit_ = l;
useHighLimit(_);
resize(l.size());
}
highLimit_ = l;
useHighLimit(_);
}
void Minimizer::setLowLimit(const Index i, const double l)
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid limit index");
}
else
{
lowLimit_(i) = l;
useLowLimit(i);
resize(i + 1);
}
lowLimit_(i) = l;
useLowLimit(i);
}
void Minimizer::setLowLimit(const PlaceHolder ph __dumb, const DVec &l)
{
checkSupport;
if (l.size() != getDim())
{
LATAN_ERROR(Size, "invalid limit vector size");
}
else
{
lowLimit_ = l;
useLowLimit(_);
resize(l.size());
}
lowLimit_ = l;
useLowLimit(_);
}
void Minimizer::useHighLimit(const Index i, const bool use)
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid limit index");
}
else
{
hasHighLimit_(i) = use;
resize(i + 1);
}
hasHighLimit_(i) = use;
}
void Minimizer::useHighLimit(const PlaceHolder ph __dumb, const bool use)
{
checkSupport;
hasHighLimit_.fill(use);
}
void Minimizer::useLowLimit(const Index i, const bool use)
{
checkSupport;
if (i >= getDim())
{
LATAN_ERROR(Size, "invalid limit index");
}
else
{
hasLowLimit_(i) = use;
resize(i + 1);
}
hasLowLimit_(i) = use;
}
void Minimizer::useLowLimit(const PlaceHolder ph __dumb, const bool use)
{
checkSupport;
hasLowLimit_.fill(use);
}
unsigned int Minimizer::getMaxPass(void) const
{
return maxPass_;
}
void Minimizer::setMaxPass(const unsigned int maxPass)
{
maxPass_ = maxPass;
}

14
lib/Minimizer.hpp
View File

@ -1,7 +1,7 @@
/*
* Minimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -36,12 +36,10 @@ class Minimizer: public Solver
public:
// constructor
Minimizer(void) = default;
explicit Minimizer(const Index dim);
// destructor
virtual ~Minimizer(void) = default;
// access
virtual void resize(const Index dim);
// limits
virtual void resize(const Index dim);
virtual double getHighLimit(const Index i) const ;
virtual const DVec & getHighLimit(const PlaceHolder ph = _) const;
virtual double getLowLimit(const Index i) const;
@ -58,11 +56,15 @@ public:
virtual void useLowLimit(const Index i, const bool use = true);
virtual void useLowLimit(const PlaceHolder ph = _,
const bool use = true);
virtual bool supportLimits(void) const = 0;
virtual unsigned int getMaxPass(void) const;
virtual void setMaxPass(const unsigned int maxPass);
// minimization
virtual const DVec & operator()(const DoubleFunction &f) = 0;
private:
DVec highLimit_, lowLimit_;
Vec<bool> hasHighLimit_, hasLowLimit_;
DVec highLimit_, lowLimit_;
Vec<bool> hasHighLimit_, hasLowLimit_;
unsigned int maxPass_{5u};
};
END_LATAN_NAMESPACE

230
lib/MinuitMinimizer.cpp
View File

@ -1,7 +1,7 @@
/*
* MinuitMinimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -19,186 +19,158 @@
#include <LatAnalyze/MinuitMinimizer.hpp>
#include <LatAnalyze/includes.hpp>
#pragma GCC diagnostic ignored "-Wconversion"
#include <Minuit2/FCNBase.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/MnPrint.h>
#include <Minuit2/MnPlot.h>
#include <Minuit2/MnScan.h>
#include <Minuit2/MnSimplex.h>
#include <Minuit2/ScanMinimizer.h>
#include <Minuit2/SimplexMinimizer.h>
#include <Minuit2/VariableMetricMinimizer.h>
#include <Minuit2/Minuit2Minimizer.h>
#include <Math/Functor.h>
using namespace std;
using namespace ROOT;
using namespace Minuit2;
using namespace Latan;
static constexpr double initErr = 0.5;
static constexpr unsigned int maxTry = 10u;
static constexpr double initErr = 0.1;
/******************************************************************************
* MinuitMinimizer implementation *
******************************************************************************/
// MinuitFunction constructor //////////////////////////////////////////////////
MinuitMinimizer::MinuitFunction::MinuitFunction(const DoubleFunction &f)
: f_(&f)
{}
// MinuitFunction minuit members ///////////////////////////////////////////////
double MinuitMinimizer::MinuitFunction::operator()
(const vector<double> &x) const
{
return (*f_)(x);
}
double MinuitMinimizer::MinuitFunction::Up(void) const
{
return 1.;
}
// constructors ////////////////////////////////////////////////////////////////
MinuitMinimizer::MinuitMinimizer(const Algorithm algorithm)
{
setAlgorithm(algorithm);
}
MinuitMinimizer::MinuitMinimizer(const Index dim, const Algorithm algorithm)
: Minimizer(dim)
{
setAlgorithm(algorithm);
}
// access //////////////////////////////////////////////////////////////////////
MinuitMinimizer::Algorithm MinuitMinimizer::getAlgorithm(void) const
{
return algorithm_;
}
double MinuitMinimizer::getPrecision(void) const
{
LATAN_ERROR(Implementation,
"Minuit minimizer precision cannot be accessed");
return 0.;
}
void MinuitMinimizer::setAlgorithm(const Algorithm algorithm)
{
algorithm_ = algorithm;
}
void MinuitMinimizer::setPrecision(const double precision __dumb)
bool MinuitMinimizer::supportLimits(void) const
{
LATAN_ERROR(Implementation,
"Minuit minimizer precision cannot be accessed");
return true;
}
// minimization ////////////////////////////////////////////////////////////////
const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
{
DVec &x = getState();
Verbosity verbosity = getVerbosity();
using namespace ROOT;
using namespace Minuit2;
DVec &x = getState();
int printLevel = 0;
EMinimizerType minuitAlg = kCombined;
// convert Latan parameters to Minuit parameters
switch (getVerbosity())
{
case Verbosity::Silent:
printLevel = 0;
break;
case Verbosity::Normal:
printLevel = 2;
break;
case Verbosity::Debug:
printLevel = 3;
break;
}
switch (getAlgorithm())
{
case Algorithm::migrad:
minuitAlg = kMigrad;
break;
case Algorithm::simplex:
minuitAlg = kSimplex;
break;
case Algorithm::combined:
minuitAlg = kCombined;
break;
}
// resize minimizer state to match function number of arguments
if (f.getNArg() != x.size())
{
resize(f.getNArg());
}
// set parameters
MinuitFunction minuitF(f);
MnUserParameters parameters;
// create and set minimizer
Minuit2Minimizer min(minuitAlg);
min.SetStrategy(2);
min.SetMaxFunctionCalls(getMaxIteration());
min.SetTolerance(getPrecision());
min.SetPrintLevel(printLevel);
// set function and variables
Math::Functor minuitF(f, x.size());
string name;
double val, step;
min.SetFunction(minuitF);
for (Index i = 0; i < x.size(); ++i)
{
parameters.Add("x_" + strFrom(i), x(i), initErr*fabs(x(i)));
if (hasLowLimit(i)&&hasHighLimit(i))
name = f.varName().getName(i);
val = x(i);
step = (fabs(x(i)) != 0.) ? initErr*fabs(x(i)) : 1.;
if (hasHighLimit(i) and !hasLowLimit(i))
{
parameters.SetLimits(static_cast<unsigned int>(i),
getLowLimit(i), getHighLimit(i));
min.SetUpperLimitedVariable(i, name, val, step, getHighLimit(i));
}
else if (hasLowLimit(i))
else if (!hasHighLimit(i) and hasLowLimit(i))
{
parameters.SetLowerLimit(static_cast<unsigned int>(i),
getLowLimit(i));
min.SetLowerLimitedVariable(i, name, val, step, getLowLimit(i));
}
else if (hasHighLimit(i))
else if (hasHighLimit(i) and hasLowLimit(i))
{
parameters.SetUpperLimit(static_cast<unsigned int>(i),
getHighLimit(i));
min.SetLimitedVariable(i, name, val, step, getLowLimit(i),
getHighLimit(i));
}
else
{
min.SetVariable(i, name, val, step);
}
}
// pre-minimization
MnSimplex preMinimizer(minuitF, parameters, 1);
FunctionMinimum preMin = preMinimizer();
// minimize
int status;
unsigned int n = 0;
if (verbosity >= Verbosity::Debug)
do
{
cout << "-- MINUIT pre-minimization -----------------------------";
cout << preMin;
cout << "--------------------------------------------------------";
cout << endl;
}
for (unsigned int i = 0; i < x.size(); ++i)
{
parameters.SetValue(i, preMin.UserParameters().Value(i));
parameters.SetError(i, 2.*preMin.UserParameters().Error(i));
}
// minimization and output
unique_ptr<MnApplication> minimizer(nullptr);
if (algorithm_ == Algorithm::Migrad)
{
minimizer.reset(new MnMigrad(minuitF, parameters, 2));
}
else if (algorithm_ == Algorithm::Simplex)
{
minimizer.reset(new MnSimplex(minuitF, parameters, 2));
}
unsigned int iTry = 0;
FunctionMinimum min = (*minimizer)();
while ((!min.IsValid())&&(iTry < maxTry))
{
min = (*minimizer)();
iTry++;
}
if (!min.IsValid())
{
LATAN_WARNING("MINUIT library reported that minimization result is not valid");
}
for (unsigned int i = 0; i < x.size(); ++i)
{
x(i) = min.UserParameters().Value(i);
}
if (verbosity >= Verbosity::Normal)
{
cout << "-- MINUIT minimization ---------------------------------";
cout << min;
cout << "--------------------------------------------------------";
cout << endl;
}
if (verbosity >= Verbosity::Debug)
{
vector<pair<double, double>> scanRes;
MnPlot plot;
MnScan scanner(minuitF, preMin.UserParameters(), 2);
cout << "-- MINUIT scan -----------------------------------------";
cout << endl;
for (unsigned int i = 0; i < x.size(); i++)
if (getVerbosity() >= Verbosity::Normal)
{
cout << "Parameter x_" << i << endl;
scanRes = scanner.Scan(i);
plot(scanRes);
cout << "========== Minuit minimization, pass #" << n + 1;
cout << " =========" << endl;
}
cout << "--------------------------------------------------------";
cout << endl;
min.Minimize();
status = min.Status();
n++;
} while (status and (n < getMaxPass()));
if (getVerbosity() >= Verbosity::Normal)
{
cout << "=================================================" << endl;
}
switch (status)
{
case 1:
LATAN_WARNING("invalid minimum: covariance matrix was made positive");
break;
case 2:
LATAN_WARNING("invalid minimum: Hesse analysis is not valid");
break;
case 3:
LATAN_WARNING("invalid minimum: requested precision not reached");
break;
case 4:
LATAN_WARNING("invalid minimum: iteration limit reached");
break;
}
// save and return result
for (Index i = 0; i < x.size(); ++i)
{
x(i) = min.X()[i];
}
return x;

42
lib/MinuitMinimizer.hpp
View File

@ -1,7 +1,7 @@
/*
* MinuitMinimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,52 +23,36 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <Minuit2/FCNBase.h>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* Minuit minimizer *
* interface to CERN Minuit minimizer *
* ( http://www.cern.ch/minuit ) *
******************************************************************************/
class MinuitMinimizer: public Minimizer
{
public:
enum class Algorithm
{
Migrad = 1,
Simplex = 2
};
private:
class MinuitFunction: public ROOT::Minuit2::FCNBase
{
public:
// constructor
explicit MinuitFunction(const DoubleFunction &f);
// destructor
virtual ~MinuitFunction(void) = default;
// minuit members
virtual double operator()(const std::vector<double> &x) const;
virtual double Up(void) const;
private:
const DoubleFunction *f_;
migrad = 1,
simplex = 2,
combined = 3
};
public:
// constructors
MinuitMinimizer(const Algorithm algorithm = Algorithm::Migrad);
explicit MinuitMinimizer(const Index dim,
const Algorithm algorithm = Algorithm::Migrad);
// constructor
explicit MinuitMinimizer(const Algorithm algorithm = defaultAlg_);
// destructor
virtual ~MinuitMinimizer(void) = default;
// access
virtual double getPrecision(void) const;
Algorithm getAlgorithm(void) const;
virtual void setPrecision(const double precision);
void setAlgorithm(const Algorithm algorithm);
Algorithm getAlgorithm(void) const;
void setAlgorithm(const Algorithm algorithm);
virtual bool supportLimits(void) const;
// minimization
virtual const DVec & operator()(const DoubleFunction &f);
private:
Algorithm algorithm_;
Algorithm algorithm_;
static constexpr Algorithm defaultAlg_ = Algorithm::combined;
};
END_LATAN_NAMESPACE

25
lib/Model.cpp
View File

@ -1,7 +1,7 @@
/*
* Model.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -19,7 +19,6 @@
#include <LatAnalyze/Model.hpp>
#include <LatAnalyze/includes.hpp>
#include <functional>
using namespace std;
using namespace std::placeholders;
@ -31,6 +30,8 @@ using namespace Latan;
// constructor /////////////////////////////////////////////////////////////////
DoubleModel::DoubleModel(const vecFunc &f, const Index nArg, const Index nPar)
: size_(new ModelSize)
, varName_("x")
, parName_("p")
{
setFunction(f, nArg, nPar);
}
@ -54,6 +55,26 @@ void DoubleModel::setFunction(const vecFunc &f, const Index nArg,
f_ = f;
}
VarName & DoubleModel::varName(void)
{
return varName_;
}
const VarName & DoubleModel::varName(void) const
{
return varName_;
}
VarName & DoubleModel::parName(void)
{
return parName_;
}
const VarName & DoubleModel::parName(void) const
{
return parName_;
}
// error checking //////////////////////////////////////////////////////////////
void DoubleModel::checkSize(const Index nArg, const Index nPar) const
{

16
lib/Model.hpp
View File

@ -1,7 +1,7 @@
/*
* Model.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,7 +23,6 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Mat.hpp>
#include <vector>
BEGIN_LATAN_NAMESPACE
@ -43,10 +42,14 @@ public:
// destructor
virtual ~DoubleModel(void) = default;
// access
virtual Index getNArg(void) const;
virtual Index getNPar(void) const;
void setFunction(const vecFunc &f, const Index nArg,
const Index nPar);
virtual Index getNArg(void) const;
virtual Index getNPar(void) const;
void setFunction(const vecFunc &f, const Index nArg,
const Index nPar);
VarName & varName(void);
const VarName & varName(void) const;
VarName & parName(void);
const VarName & parName(void) const;
// function call
double operator()(const DVec &data, const DVec &par) const;
double operator()(const std::vector<double> &data,
@ -61,6 +64,7 @@ private:
void checkSize(const Index nArg, const Index nPar) const;
private:
std::shared_ptr<ModelSize> size_;
VarName varName_, parName_;
vecFunc f_;
};

201
lib/NloptMinimizer.cpp Normal file
View File

@ -0,0 +1,201 @@
/*
* NloptMinimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 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 <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/NloptMinimizer.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* NloptMinimizer implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
NloptMinimizer::NloptMinimizer(const Algorithm algorithm)
{
setAlgorithm(algorithm);
der_.setOrder(1, 1);
}
// access //////////////////////////////////////////////////////////////////////
NloptMinimizer::Algorithm NloptMinimizer::getAlgorithm(void) const
{
return algorithm_;
}
void NloptMinimizer::setAlgorithm(const Algorithm algorithm)
{
algorithm_ = algorithm;
}
bool NloptMinimizer::supportLimits(void) const
{
return true;
}
// minimization ////////////////////////////////////////////////////////////////
const DVec & NloptMinimizer::operator()(const DoubleFunction &f)
{
DVec &x = getState();
// resize minimizer state to match function number of arguments
if (f.getNArg() != x.size())
{
resize(f.getNArg());
}
// create and set minimizer
nlopt::opt min(getAlgorithm(), x.size());
NloptFuncData data;
vector<double> lb(x.size()), hb(x.size());
min.set_maxeval(getMaxIteration());
min.set_xtol_rel(getPrecision());
min.set_ftol_rel(-1.);
der_.setFunction(f);
data.f = &f;
data.d = &der_;
min.set_min_objective(&funcWrapper, &data);
for (Index i = 0; i < x.size(); ++i)
{
lb[i] = hasLowLimit(i) ? getLowLimit(i) : -HUGE_VAL;
hb[i] = hasHighLimit(i) ? getHighLimit(i) : HUGE_VAL;
}
min.set_lower_bounds(lb);
min.set_upper_bounds(hb);
// minimize
double res;
vector<double> vx(x.size());
nlopt::result status;
unsigned int n = 0;
for (Index i = 0; i < x.size(); ++i)
{
vx[i] = x(i);
}
do
{
if (getVerbosity() >= Verbosity::Normal)
{
cout << "========== NLopt minimization, pass #" << n + 1;
cout << " ==========" << endl;
cout << "Algorithm: " << min.get_algorithm_name() << endl;
cout << "Max eval.= " << min.get_maxeval();
cout << " -- Precision= " << min.get_xtol_rel() << endl;
printf("Starting f(x)= %.10e\n", f(x));
}
try
{
status = min.optimize(vx, res);
}
catch (invalid_argument &e)
{
LATAN_ERROR(Runtime, "NLopt has reported receving invalid "
"arguments (if you are using a global minimizer, did "
"you specify limits for all variables?)");
}
if (getVerbosity() >= Verbosity::Normal)
{
printf("Found minimum %.10e at:\n", res);
for (Index i = 0; i < x.size(); ++i)
{
printf("%8s= %.10e\n", f.varName().getName(i).c_str(), vx[i]);
}
cout << "after " << data.evalCount << " evaluations" << endl;
cout << "Minimization ended with code " << status;
cout << " (" << returnMessage(status) << ")";
cout << endl;
}
data.evalCount = 0;
for (Index i = 0; i < x.size(); ++i)
{
x(i) = vx[i];
}
n++;
} while (!minSuccess(status) and (n < getMaxPass()));
if (getVerbosity() >= Verbosity::Normal)
{
cout << "=================================================" << endl;
}
if (!minSuccess(status))
{
LATAN_WARNING("invalid minimum: " + returnMessage(status));
}
return x;
}
// NLopt return code parser ////////////////////////////////////////////////////
string NloptMinimizer::returnMessage(const nlopt::result status)
{
switch (status)
{
case nlopt::SUCCESS:
return "success";
case nlopt::STOPVAL_REACHED:
return "stopping value reached";
case nlopt::FTOL_REACHED:
return "tolerance on function reached";
case nlopt::XTOL_REACHED:
return "tolerance on variable reached";
case nlopt::MAXEVAL_REACHED:
return "maximum function evaluation reached";
case nlopt::MAXTIME_REACHED:
return "maximum time reached";
default:
return "";
}
}
// NLopt function wrapper //////////////////////////////////////////////////////
double NloptMinimizer::funcWrapper(unsigned int n, const double *arg,
double *grad , void *vdata)
{
NloptFuncData &data = *static_cast<NloptFuncData *>(vdata);
if (grad)
{
for (unsigned int i = 0; i < n; ++i)
{
data.d->setDir(i);
grad[i] = (*(data.d))(arg);
}
data.evalCount += data.d->getNPoint()*n;
}
data.evalCount++;
return (*data.f)(arg);
}
// NLopt return status parser //////////////////////////////////////////////////
bool NloptMinimizer::minSuccess(const nlopt::result status)
{
switch (status)
{
case nlopt::SUCCESS:
case nlopt::FTOL_REACHED:
case nlopt::XTOL_REACHED:
return true;
break;
default:
return false;
break;
}
}

76
lib/NloptMinimizer.hpp Normal file
View File

@ -0,0 +1,76 @@
/*
* NloptMinimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 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 <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_NloptMinimizer_hpp_
#define Latan_NloptMinimizer_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Derivative.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <nlopt.hpp>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* interface to NLOpt minimizers *
* ( http://ab-initio.mit.edu/wiki/index.php/NLopt ) *
* -------------------------------------------------------------------------- *
* cf. http://ab-initio.mit.edu/wiki/index.php/NLopt_Algorithms for algorithm *
* references and naming conventions *
******************************************************************************/
class NloptMinimizer: public Minimizer
{
public:
typedef nlopt::algorithm Algorithm;
private:
struct NloptFuncData
{
const DoubleFunction *f{nullptr};
Derivative *d{nullptr};
unsigned int evalCount{0};
};
public:
// constructor
explicit NloptMinimizer(const Algorithm algorithm = defaultAlg_);
// destructor
virtual ~NloptMinimizer(void) = default;
// access
Algorithm getAlgorithm(void) const;
void setAlgorithm(const Algorithm algorithm);
virtual bool supportLimits(void) const;
// minimization
virtual const DVec & operator()(const DoubleFunction &f);
private:
// NLopt return code parser
static std::string returnMessage(const nlopt::result status);
// NLopt function wrapper
static double funcWrapper(unsigned int n, const double *arg,
double *grad , void *vdata);
// NLopt return status parser
static bool minSuccess(const nlopt::result status);
private:
Algorithm algorithm_;
static constexpr Algorithm defaultAlg_ = Algorithm::LN_NELDERMEAD;
CentralDerivative der_;
};
END_LATAN_NAMESPACE
#endif // Latan_NloptMinimizer_hpp_

3
lib/ParserState.hpp
View File

@ -1,7 +1,7 @@
/*
* ParserState.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -21,7 +21,6 @@
#define Latan_ParserState_hpp_
#include <LatAnalyze/Global.hpp>
#include <iostream>
BEGIN_LATAN_NAMESPACE

17
lib/Plot.cpp
View File

@ -1,7 +1,7 @@
/*
* Plot.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -169,15 +169,10 @@ PlotData::PlotData(const DMatSample &x, const DVec &y)
PlotData::PlotData(const XYStatData &data, const Index i, const Index j)
{
DMat d(data.getNData(), 4);
string usingCmd, tmpFileName;
d.col(0) = data.x(i);
d.col(2) = data.y(j);
d.col(1) = data.xxVar(i, i).diagonal().array().sqrt();
d.col(3) = data.yyVar(j, j).diagonal().array().sqrt();
usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:3:2:4 w xyerr";
tmpFileName = dumpToTmpFile(d);
tmpFileName = dumpToTmpFile(data.getTable(i, j));
pushTmpFile(tmpFileName);
setCommand("'" + tmpFileName + "' " + usingCmd);
}
@ -252,6 +247,14 @@ PlotHistogram::PlotHistogram(const Histogram &h)
setCommand("'" + tmpFileName + "' u 1:2 w steps");
}
// PlotMatrixNoRange constructor ///////////////////////////////////////////////
PlotMatrixNoRange::PlotMatrixNoRange(const DMat &m)
{
string tmpFileName = dumpToTmpFile(m);
setCommand("'" + tmpFileName + "' matrix w image");
}
/******************************************************************************
* Plot modifiers *
******************************************************************************/

28
lib/Plot.hpp
View File

@ -1,7 +1,7 @@
/*
* Plot.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -21,11 +21,11 @@
#define Latan_Plot_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/Histogram.hpp>
#include <LatAnalyze/XYStatData.hpp>
#include <sstream>
#include <vector>
// gnuplot default parameters
#ifndef GNUPLOT_BIN
@ -47,11 +47,11 @@ public:
// destructor
virtual ~PlotObject(void) = default;
// access
std::string popTmpFile(void);
const std::string & getCommand(void) const;
const std::string & getHeadCommand(void) const;
std::string popTmpFile(void);
const std::string & getCommand(void) const;
const std::string & getHeadCommand(void) const;
// test
bool gotTmpFile(void) const;
bool gotTmpFile(void) const;
protected:
// access
void pushTmpFile(const std::string &fileName);
@ -136,6 +136,20 @@ public:
virtual ~PlotHistogram(void) = default;
};
class PlotMatrixNoRange: public PlotObject
{
public:
// constructor
PlotMatrixNoRange(const DMat &m);
// destructor
virtual ~PlotMatrixNoRange(void) = default;
};
#define PlotMatrix(m)\
PlotRange(Axis::x, -.5, (m).cols() - .5) <<\
PlotRange(Axis::y, (m).rows() - .5, -.5) <<\
PlotMatrixNoRange(m)
/******************************************************************************
* Plot modifiers *
******************************************************************************/

690
lib/RandGen.cpp
View File

@ -1,690 +0,0 @@
/*
* RandGen.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/RandGen.hpp>
#include <LatAnalyze/includes.hpp>
#pragma GCC diagnostic ignored "-Wfloat-conversion"
#ifndef RLXD_LEVEL
#define RLXD_LEVEL 1
#endif
using namespace std;
using namespace Latan;
/******************************************************************************
* RandGenState implementation *
******************************************************************************/
// IO type ///////////////////////////////////////////////////////////////
IoObject::IoType RandGenState::getType(void) const
{
return IoType::rgState;
}
/******************************************************************************
* RandGen implementation *
******************************************************************************/
// RanLxd implementation ///////////////////////////////////////////////////////
RandGen::RanLxd::RanLxd(void)
{
// avoid a warning in the SSE case
one_bit = 0.0;
}
/****************************************************************************/
/* Copyright (C) 2005 Martin Luescher (GPL)
* This software is distributed under the terms of the GNU General Public
* License (GPL)
*
* Random number generator "ranlxd". See the notes
*
* "User's guide for ranlxs and ranlxd v3.2" (December 2005)
*
* "Algorithms used in ranlux v3.0" (May 2001)
*
* for a detailed description
*
* The functions are
*
* void ranlxd(double r[],int n)
* Computes the next n double-precision random numbers and
* assigns them to the elements r[0],...,r[n-1] of the array r[]
*
* void rlxd_init(int level,int seed)
* Initialization of the generator
*
* int rlxd_size(void)
* Returns the number of integers required to save the state of
* the generator
*
* void rlxd_get(int state[])
* Extracts the current state of the generator and stores the
* information in the array state[N] where N>=rlxd_size()
*
* void rlxd_reset(int state[])
* Resets the generator to the state defined by the array state[N]
*/
#ifdef HAVE_SSE
#define STEP(pi,pj) \
__asm__ __volatile__ ("movaps %4, %%xmm4 \n\t" \
"movaps %%xmm2, %%xmm3 \n\t" \
"subps %2, %%xmm4 \n\t" \
"movaps %%xmm1, %%xmm5 \n\t" \
"cmpps $0x6, %%xmm4, %%xmm2 \n\t" \
"andps %%xmm2, %%xmm5 \n\t" \
"subps %%xmm3, %%xmm4 \n\t" \
"andps %%xmm0, %%xmm2 \n\t" \
"addps %%xmm4, %%xmm5 \n\t" \
"movaps %%xmm5, %0 \n\t" \
"movaps %5, %%xmm6 \n\t" \
"movaps %%xmm2, %%xmm3 \n\t" \
"subps %3, %%xmm6 \n\t" \
"movaps %%xmm1, %%xmm7 \n\t" \
"cmpps $0x6, %%xmm6, %%xmm2 \n\t" \
"andps %%xmm2, %%xmm7 \n\t" \
"subps %%xmm3, %%xmm6 \n\t" \
"andps %%xmm0, %%xmm2 \n\t" \
"addps %%xmm6, %%xmm7 \n\t" \
"movaps %%xmm7, %1" \
: \
"=m" ((*pi).c1), \
"=m" ((*pi).c2) \
: \
"m" ((*pi).c1), \
"m" ((*pi).c2), \
"m" ((*pj).c1), \
"m" ((*pj).c2) \
: \
"xmm2", "xmm3", "xmm4", "xmm5", "xmm6", "xmm7")
void RandGen::RanLxd::error(int no) const
{
switch(no)
{
case 1:
LATAN_ERROR(Range, "Bad choice of luxury level (should be 1 or 2)");
break;
case 2:
LATAN_ERROR(Range, "Bad choice of seed (should be between 1 and 2^31-1)");
break;
case 3:
LATAN_ERROR(Runtime, "Undefined state (ranlxd is not initialized)");
break;
case 5:
LATAN_ERROR(Logic, "Unexpected input data");
break;
}
}
void RandGen::RanLxd::update(void)
{
int k,kmax;
rlxd_dble_vec_t *pmin,*pmax,*pi,*pj;
kmax=rlxd_pr;
pmin=&rlxd_x.vec[0];
pmax=pmin+12;
pi=&rlxd_x.vec[ir];
pj=&rlxd_x.vec[jr];
__asm__ __volatile__ ("movaps %0, %%xmm0 \n\t"
"movaps %1, %%xmm1 \n\t"
"movaps %2, %%xmm2"
:
:
"m" (one_bit_sse),
"m" (one_sse),
"m" (carry)
:
"xmm0", "xmm1", "xmm2");
for (k=0;k<kmax;k++)
{
STEP(pi,pj);
pi+=1;
pj+=1;
if (pi==pmax)
pi=pmin;
if (pj==pmax)
pj=pmin;
}
__asm__ __volatile__ ("movaps %%xmm2, %0"
:
"=m" (carry));
ir+=prm;
jr+=prm;
if (ir>=12)
ir-=12;
if (jr>=12)
jr-=12;
is=8*ir;
is_old=is;
}
void RandGen::RanLxd::define_constants(void)
{
int k;
float b;
one_sse.c1=1.0f;
one_sse.c2=1.0f;
one_sse.c3=1.0f;
one_sse.c4=1.0f;
b=(float)(ldexp(1.0,-24));
one_bit_sse.c1=b;
one_bit_sse.c2=b;
one_bit_sse.c3=b;
one_bit_sse.c4=b;
for (k=0;k<96;k++)
{
next[k]=(k+1)%96;
if ((k%4)==3)
next[k]=(k+5)%96;
}
}
void RandGen::RanLxd::rlxd_init(int level,int seed)
{
int i,k,l;
int ibit,jbit,xbit[31];
int ix,iy;
define_constants();
if (level==1)
rlxd_pr=202;
else if (level==2)
rlxd_pr=397;
else
error(1);
i=seed;
for (k=0;k<31;k++)
{
xbit[k]=i%2;
i/=2;
}
if ((seed<=0)||(i!=0))
error(2);
ibit=0;
jbit=18;
for (i=0;i<4;i++)
{
for (k=0;k<24;k++)
{
ix=0;
for (l=0;l<24;l++)
{
iy=xbit[ibit];
ix=2*ix+iy;
xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
ibit=(ibit+1)%31;
jbit=(jbit+1)%31;
}
if ((k%4)!=i)
ix=16777215-ix;
rlxd_x.num[4*k+i]=(float)(ldexp((double)(ix),-24));
}
}
carry.c1=0.0f;
carry.c2=0.0f;
carry.c3=0.0f;
carry.c4=0.0f;
ir=0;
jr=7;
is=91;
is_old=0;
prm=rlxd_pr%12;
init=1;
}
void RandGen::RanLxd::ranlxd(double r[],int n)
{
int k;
if (init==0)
rlxd_init(1,1);
for (k=0;k<n;k++)
{
is=next[is];
if (is==is_old)
update();
r[k]=(double)(rlxd_x.num[is+4])+(double)(one_bit_sse.c1*rlxd_x.num[is]);
}
}
int RandGen::RanLxd::rlxd_size(void) const
{
return(RLXG_STATE_SIZE);
}
void RandGen::RanLxd::rlxd_get(int state[]) const
{
int k;
float base;
if (init==0)
error(3);
base=(float)(ldexp(1.0,24));
state[0]=rlxd_size();
for (k=0;k<96;k++)
state[k+1]=(int)(base*rlxd_x.num[k]);
state[97]=(int)(base*carry.c1);
state[98]=(int)(base*carry.c2);
state[99]=(int)(base*carry.c3);
state[100]=(int)(base*carry.c4);
state[101]=rlxd_pr;
state[102]=ir;
state[103]=jr;
state[104]=is;
}
void RandGen::RanLxd::rlxd_reset(const int state[])
{
int k;
define_constants();
if (state[0]!=rlxd_size())
error(5);
for (k=0;k<96;k++)
{
if ((state[k+1]<0)||(state[k+1]>=167777216))
error(5);
rlxd_x.num[k]=(float)(ldexp((double)(state[k+1]),-24));
}
if (((state[97]!=0)&&(state[97]!=1))||
((state[98]!=0)&&(state[98]!=1))||
((state[99]!=0)&&(state[99]!=1))||
((state[100]!=0)&&(state[100]!=1)))
error(5);
carry.c1=(float)(ldexp((double)(state[97]),-24));
carry.c2=(float)(ldexp((double)(state[98]),-24));
carry.c3=(float)(ldexp((double)(state[99]),-24));
carry.c4=(float)(ldexp((double)(state[100]),-24));
rlxd_pr=state[101];
ir=state[102];
jr=state[103];
is=state[104];
is_old=8*ir;
prm=rlxd_pr%12;
init=1;
if (((rlxd_pr!=202)&&(rlxd_pr!=397))||
(ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))||
(is<0)||(is>91))
error(5);
}
#else
#define BASE 0x1000000
#define MASK 0xffffff
#define STEP(pi,pj) \
d=(*pj).c1.c1-(*pi).c1.c1-carry.c1; \
(*pi).c2.c1+=(d<0); \
d+=BASE; \
(*pi).c1.c1=d&MASK; \
d=(*pj).c1.c2-(*pi).c1.c2-carry.c2; \
(*pi).c2.c2+=(d<0); \
d+=BASE; \
(*pi).c1.c2=d&MASK; \
d=(*pj).c1.c3-(*pi).c1.c3-carry.c3; \
(*pi).c2.c3+=(d<0); \
d+=BASE; \
(*pi).c1.c3=d&MASK; \
d=(*pj).c1.c4-(*pi).c1.c4-carry.c4; \
(*pi).c2.c4+=(d<0); \
d+=BASE; \
(*pi).c1.c4=d&MASK; \
d=(*pj).c2.c1-(*pi).c2.c1; \
carry.c1=(d<0); \
d+=BASE; \
(*pi).c2.c1=d&MASK; \
d=(*pj).c2.c2-(*pi).c2.c2; \
carry.c2=(d<0); \
d+=BASE; \
(*pi).c2.c2=d&MASK; \
d=(*pj).c2.c3-(*pi).c2.c3; \
carry.c3=(d<0); \
d+=BASE; \
(*pi).c2.c3=d&MASK; \
d=(*pj).c2.c4-(*pi).c2.c4; \
carry.c4=(d<0); \
d+=BASE; \
(*pi).c2.c4=d&MASK
void RandGen::RanLxd::error(int no) const
{
switch(no)
{
case 0:
LATAN_ERROR(System, "Arithmetic on this machine is not suitable for ranlxd");
break;
case 1:
LATAN_ERROR(Range, "Bad choice of luxury level (should be 1 or 2)");
break;
case 2:
LATAN_ERROR(Range, "Bad choice of seed (should be between 1 and 2^31-1)");
break;
case 3:
LATAN_ERROR(Runtime, "Undefined state (ranlxd is not initialized)");
case 4:
LATAN_ERROR(System, "Arithmetic on this machine is not suitable for ranlxd");
break;
case 5:
LATAN_ERROR(Logic, "Unexpected input data");
break;
}
}
void RandGen::RanLxd::update(void)
{
int k,kmax,d;
rlxd_dble_vec_t *pmin,*pmax,*pi,*pj;
kmax=rlxd_pr;
pmin=&rlxd_x.vec[0];
pmax=pmin+12;
pi=&rlxd_x.vec[ir];
pj=&rlxd_x.vec[jr];
for (k=0;k<kmax;k++)
{
STEP(pi,pj);
pi+=1;
pj+=1;
if (pi==pmax)
pi=pmin;
if (pj==pmax)
pj=pmin;
}
ir+=prm;
jr+=prm;
if (ir>=12)
ir-=12;
if (jr>=12)
jr-=12;
is=8*ir;
is_old=is;
}
void RandGen::RanLxd::define_constants(void)
{
int k;
one_bit=ldexp(1.0,-24);
for (k=0;k<96;k++)
{
next[k]=(k+1)%96;
if ((k%4)==3)
next[k]=(k+5)%96;
}
}
void RandGen::RanLxd::rlxd_init(int level,int seed)
{
int i,k,l;
int ibit,jbit,xbit[31];
int ix,iy;
if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)||
(DBL_MANT_DIG<48))
error(0);
define_constants();
if (level==1)
rlxd_pr=202;
else if (level==2)
rlxd_pr=397;
else
error(1);
i=seed;
for (k=0;k<31;k++)
{
xbit[k]=i%2;
i/=2;
}
if ((seed<=0)||(i!=0))
error(2);
ibit=0;
jbit=18;
for (i=0;i<4;i++)
{
for (k=0;k<24;k++)
{
ix=0;
for (l=0;l<24;l++)
{
iy=xbit[ibit];
ix=2*ix+iy;
xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
ibit=(ibit+1)%31;
jbit=(jbit+1)%31;
}
if ((k%4)!=i)
ix=16777215-ix;
rlxd_x.num[4*k+i]=ix;
}
}
carry.c1=0;
carry.c2=0;
carry.c3=0;
carry.c4=0;
ir=0;
jr=7;
is=91;
is_old=0;
prm=rlxd_pr%12;
init=1;
}
void RandGen::RanLxd::ranlxd(double r[],int n)
{
int k;
if (init==0)
rlxd_init(1,1);
for (k=0;k<n;k++)
{
is=next[is];
if (is==is_old)
update();
r[k]=one_bit*((double)(rlxd_x.num[is+4])+one_bit*(double)(rlxd_x.num[is]));
}
}
int RandGen::RanLxd::rlxd_size(void) const
{
return(RLXG_STATE_SIZE);
}
void RandGen::RanLxd::rlxd_get(int state[]) const
{
int k;
if (init==0)
error(3);
state[0]=rlxd_size();
for (k=0;k<96;k++)
state[k+1]=rlxd_x.num[k];
state[97]=carry.c1;
state[98]=carry.c2;
state[99]=carry.c3;
state[100]=carry.c4;
state[101]=rlxd_pr;
state[102]=ir;
state[103]=jr;
state[104]=is;
}
void RandGen::RanLxd::rlxd_reset(const int state[])
{
int k;
if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)||
(DBL_MANT_DIG<48))
error(4);
define_constants();
if (state[0]!=rlxd_size())
error(5);
for (k=0;k<96;k++)
{
if ((state[k+1]<0)||(state[k+1]>=167777216))
error(5);
rlxd_x.num[k]=state[k+1];
}
if (((state[97]!=0)&&(state[97]!=1))||
((state[98]!=0)&&(state[98]!=1))||
((state[99]!=0)&&(state[99]!=1))||
((state[100]!=0)&&(state[100]!=1)))
error(5);
carry.c1=state[97];
carry.c2=state[98];
carry.c3=state[99];
carry.c4=state[100];
rlxd_pr=state[101];
ir=state[102];
jr=state[103];
is=state[104];
is_old=8*ir;
prm=rlxd_pr%12;
init=1;
if (((rlxd_pr!=202)&&(rlxd_pr!=397))||
(ir<0)||(ir>11)||(jr<0)||(jr>11)||(jr!=((ir+7)%12))||
(is<0)||(is>91))
error(5);
}
#endif
// constructors ////////////////////////////////////////////////////////////////
RandGen::RandGen(void)
{
generator_.rlxd_init(RLXD_LEVEL, (int)time(NULL));
}
RandGen::RandGen(const int seed)
{
generator_.rlxd_init(RLXD_LEVEL, seed);
}
RandGen::RandGen(const RandGenState &state)
{
setState(state);
}
// state management ////////////////////////////////////////////////////////////
RandGenState RandGen::getState(void) const
{
RandGenState state;
generator_.rlxd_get(state.data());
return state;
}
void RandGen::setState(const RandGenState &state)
{
generator_.rlxd_reset(state.data());
}
// generators //////////////////////////////////////////////////////////////////
double RandGen::uniform(const double a, const double b)
{
double rx;
generator_.ranlxd(&rx, 1);
return (b-a)*rx + a;
}
unsigned int RandGen::discreteUniform(const unsigned int n)
{
return ((unsigned int)(uniform()*(double)(n)));
}
double RandGen::gaussian(const double mean, const double sigma)
{
double rx, ry, sqNrm;
do
{
rx = uniform(-1.0, 1.0);
ry = uniform(-1.0, 1.0);
sqNrm = rx*rx + ry*ry;
} while ((sqNrm > 1.0)||(sqNrm == 0.0));
return sigma*rx*sqrt(-2.0*log(sqNrm)/sqNrm) + mean;
}

103
lib/RandGen.hpp
View File

@ -1,103 +0,0 @@
/*
* RandGen.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 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 <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_RandGen_hpp_
#define Latan_RandGen_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/IoObject.hpp>
#define RLXG_STATE_SIZE 105u
BEGIN_LATAN_NAMESPACE
class RandGenState: public Array<int, RLXG_STATE_SIZE, 1>,
public IoObject
{
private:
typedef Array<int, RLXG_STATE_SIZE, 1> Base;
public:
// destructor
virtual ~RandGenState(void) = default;
// IO type
IoType getType(void) const;
};
/******************************************************************************
* Random generator class *
******************************************************************************/
class RandGen
{
private:
// Martin Luescher's ranlxd generator interface
class RanLxd
{
private:
typedef struct alignas(16)
{
float c1,c2,c3,c4;
} rlxd_vec_t;
typedef struct alignas(16)
{
rlxd_vec_t c1,c2;
} rlxd_dble_vec_t;
public:
RanLxd(void);
~RanLxd(void) = default;
void ranlxd(double r[],int n);
void rlxd_init(int level,int seed);
int rlxd_size(void) const;
void rlxd_get(int state[]) const;
void rlxd_reset(const int state[]);
private:
void error(int no) const;
void update(void);
void define_constants(void);
private:
int init{0}, rlxd_pr, prm, ir, jr, is, is_old, next[96];
rlxd_vec_t one_sse, one_bit_sse, carry;
double one_bit;
union alignas(16)
{
rlxd_dble_vec_t vec[12];
float num[96];
} rlxd_x;
};
public:
// constructors
RandGen(void);
explicit RandGen(const int seed);
explicit RandGen(const RandGenState &state);
// destructor
virtual ~RandGen(void) = default;
// state management
RandGenState getState(void) const;
void setState(const RandGenState &state);
// generators
double uniform(const double a = 0.0, const double b = 1.0);
unsigned int discreteUniform(const unsigned int n);
double gaussian(const double mean = 0.0, const double sigma = 1.0);
private:
RanLxd generator_;
};
END_LATAN_NAMESPACE
#endif // Latan_RandGen_hpp_

2
lib/RootFinder.cpp
View File

@ -1,7 +1,7 @@
/*
* RootFinder.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/RootFinder.hpp
View File

@ -1,7 +1,7 @@
/*
* RootFinder.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

2
lib/Solver.cpp
View File

@ -1,7 +1,7 @@
/*
* Solver.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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

4
lib/Solver.hpp
View File

@ -1,7 +1,7 @@
/*
* Solver.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -31,7 +31,7 @@ BEGIN_LATAN_NAMESPACE
class Solver
{
public:
static const unsigned int defaultMaxIteration = 1000u;
static const unsigned int defaultMaxIteration = 10000u;
static constexpr double defaultPrec = 1.0e-7;
public:
enum class Verbosity

8
lib/MatSample.cpp → lib/StatArray.cpp
View File

@ -1,7 +1,7 @@
/*
* MatSample.cpp, part of LatAnalyze 3
* StatArray.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -17,7 +17,7 @@
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/StatArray.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
@ -25,7 +25,7 @@ using namespace std;
namespace Latan
{
template <>
IoObject::IoType MatSample<double>::getType(void) const
IoObject::IoType StatArray<Mat<double>, -1>::getType(void) const
{
return IoType::dMatSample;
}

22
lib/StatArray.hpp
View File

@ -1,7 +1,7 @@
/*
* StatArray.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -22,7 +22,6 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Mat.hpp>
#include <iostream>
#define FOR_STAT_ARRAY(ar, i) \
for (Latan::Index i = -(ar).offset; i < (ar).size(); ++i)
@ -33,7 +32,7 @@ BEGIN_LATAN_NAMESPACE
* Array class with statistics *
******************************************************************************/
template <typename T, Index os = 0>
class StatArray: public Array<T, dynamic, 1>
class StatArray: public Array<T, dynamic, 1>, public IoObject
{
protected:
typedef Array<T, dynamic, 1> Base;
@ -60,8 +59,10 @@ public:
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;
// IO type
virtual IoType getType(void) const;
public:
static const Index offset = os;
static constexpr Index offset = os;
};
// reduction operations
@ -77,12 +78,10 @@ namespace ReducOp
}
// Sample types
#define SAMPLE_OFFSET 1
const int central = -SAMPLE_OFFSET;
const int central = -1;
template <typename T>
using Sample = StatArray<T, SAMPLE_OFFSET>;
using Sample = StatArray<T, 1>;
typedef Sample<double> DSample;
typedef Sample<std::complex<double>> CSample;
@ -273,6 +272,13 @@ namespace ReducOp
}
}
// IO type /////////////////////////////////////////////////////////////////////
template <typename T, Index os>
IoObject::IoType StatArray<T, os>::getType(void) const
{
return IoType::noType;
}
END_LATAN_NAMESPACE
#endif // Latan_StatArray_hpp_

32
lib/TabFunction.cpp
View File

@ -1,7 +1,7 @@
/*
* TabFunction.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -34,12 +34,12 @@ TabFunction::TabFunction(const DVec &x, const DVec &y,
setData(x, y);
}
TabFunction::TabFunction(const XYStatData &data, const Index i, const Index j,
const InterpType interpType)
: interpType_(interpType)
{
setData(data, i, j);
}
//TabFunction::TabFunction(const XYStatData &data, const Index i, const Index j,
// const InterpType interpType)
//: interpType_(interpType)
//{
// setData(data, i, j);
//}
// access //////////////////////////////////////////////////////////////////////
void TabFunction::setData(const DVec &x, const DVec &y)
@ -54,10 +54,10 @@ void TabFunction::setData(const DVec &x, const DVec &y)
}
}
void TabFunction::setData(const XYStatData &data, const Index i, const Index j)
{
setData(data.x(i), data.y(j));
}
//void TabFunction::setData(const XYStatData &data, const Index i, const Index j)
//{
// setData(data.x(i), data.y(j));
//}
// function call ///////////////////////////////////////////////////////////////
double TabFunction::operator()(const double *arg) const
@ -154,11 +154,11 @@ DoubleFunction Latan::interpolate(const DVec &x, const DVec &y,
return TabFunction(x, y, interpType).makeFunction();
}
DoubleFunction Latan::interpolate(const XYStatData &data, const Index i,
const Index j, const InterpType interpType)
{
return TabFunction(data, i, j, interpType).makeFunction();
}
//DoubleFunction Latan::interpolate(const XYStatData &data, const Index i,
// const Index j, const InterpType interpType)
//{
// return TabFunction(data, i, j, interpType).makeFunction();
//}
map<double, double>::const_iterator TabFunction::nearest(const double x) const
{

16
lib/TabFunction.hpp
View File

@ -1,7 +1,7 @@
/*
* TabFunction.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -20,8 +20,6 @@
#ifndef Latan_TabFunction_hpp_
#define Latan_TabFunction_hpp_
#include <algorithm>
#include <map>
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Math.hpp>
@ -47,13 +45,13 @@ public:
TabFunction(void) = default;
TabFunction(const DVec &x, const DVec &y,
const InterpType interpType = InterpType::LINEAR);
TabFunction(const XYStatData &data, const Index i = 0, const Index j = 0,
const InterpType interpType = InterpType::LINEAR);
//TabFunction(const XYStatData &data, const Index i = 0, const Index j = 0,
// const InterpType interpType = InterpType::LINEAR);
// destructor
virtual ~TabFunction(void) = default;
// access
void setData(const DVec &x, const DVec &y);
void setData(const XYStatData &data, const Index i = 0, const Index j = 0);
// void setData(const XYStatData &data, const Index i = 0, const Index j = 0);
// function call
double operator()(const double *arg) const;
// factory
@ -67,9 +65,9 @@ private:
DoubleFunction interpolate(const DVec &x, const DVec &y,
const InterpType interpType = InterpType::LINEAR);
DoubleFunction interpolate(const XYStatData &data, const Index i = 0,
const Index j = 0,
const InterpType interpType = InterpType::LINEAR);
//DoubleFunction interpolate(const XYStatData &data, const Index i = 0,
// const Index j = 0,
// const InterpType interpType = InterpType::LINEAR);
END_LATAN_NAMESPACE

596
lib/XYSampleData.cpp
View File

@ -1,7 +1,7 @@
/*
* XYSampleData.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -18,12 +18,11 @@
*/
#include <LatAnalyze/XYSampleData.hpp>
#include <LatAnalyze/Math.hpp>
#include <LatAnalyze/includes.hpp>
#include <LatAnalyze/Math.hpp>
using namespace std;
using namespace Latan;
using namespace Math;
/******************************************************************************
* SampleFitResult implementation *
@ -53,28 +52,33 @@ double SampleFitResult::getNDof(void) const
return static_cast<double>(nDof_);
}
double SampleFitResult::getPValue(const Index s) const
Index SampleFitResult::getNPar(void) const
{
return chi2PValue(getChi2(s), getNDof());
return nPar_;
}
const DoubleFunction & SampleFitResult::getModel(const Index s,
double SampleFitResult::getPValue(const Index s) const
{
return Math::chi2PValue(getChi2(s), getNDof());
}
const DoubleFunction & SampleFitResult::getModel(const Index s,
const Index j) const
{
return model_[static_cast<unsigned int>(j)][s];
return model_[j][s];
}
const DoubleFunctionSample & SampleFitResult::getModel(
const PlaceHolder ph __dumb,
const Index j) const
const PlaceHolder ph __dumb,
const Index j) const
{
return model_[static_cast<unsigned int>(j)];
return model_[j];
}
FitResult SampleFitResult::getFitResult(const Index s) const
{
FitResult fit;
fit = (*this)[s];
fit.chi2_ = getChi2();
fit.nDof_ = static_cast<Index>(getNDof());
@ -83,292 +87,406 @@ FitResult SampleFitResult::getFitResult(const Index s) const
{
fit.model_[k] = model_[k][s];
}
return fit;
}
// IO //////////////////////////////////////////////////////////////////////////
void SampleFitResult::print(const bool printXsi, ostream &out) const
{
char buf[256];
Index pMax = printXsi ? size() : nPar_;
DMat err = this->variance().cwiseSqrt();
sprintf(buf, "chi^2/dof= %.1f/%d= %.2f -- p-value= %.2e", getChi2(),
static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
out << buf << endl;
for (Index p = 0; p < pMax; ++p)
{
sprintf(buf, "%8s= % e +/- %e", parName_[p].c_str(),
(*this)[central](p), err(p));
out << buf << endl;
}
}
/******************************************************************************
* XYSampleData implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
XYSampleData::XYSampleData(const Index nData, const Index xDim,
const Index yDim, const Index nSample)
// constructor /////////////////////////////////////////////////////////////////
XYSampleData::XYSampleData(const Index nSample)
: nSample_(nSample)
{}
// data access /////////////////////////////////////////////////////////////////
DSample & XYSampleData::x(const Index r, const Index i)
{
resize(nData, xDim, yDim, nSample);
checkXIndex(r, i);
scheduleDataInit();
scheduleComputeVarMat();
if (xData_[i][r].size() == 0)
{
xData_[i][r].resize(nSample_);
}
return xData_[i][r];
}
// access //////////////////////////////////////////////////////////////////////
const XYStatData & XYSampleData::getData(const Index s)
const DSample & XYSampleData::x(const Index r, const Index i) const
{
setDataToSample(s);
checkXIndex(r, i);
return xData_[i][r];
}
const DMatSample & XYSampleData::x(const Index k)
{
checkDataIndex(k);
updateXMap();
return xMap_.at(k);
}
DSample & XYSampleData::y(const Index k, const Index j)
{
checkYDim(j);
if (!pointExists(k, j))
{
registerDataPoint(k, j);
}
scheduleDataInit();
scheduleComputeVarMat();
if (yData_[j][k].size() == 0)
{
yData_[j][k].resize(nSample_);
}
return yData_[j][k];
}
const DSample & XYSampleData::y(const Index k, const Index j) const
{
checkPoint(k, j);
return yData_[j].at(k);
}
const DMat & XYSampleData::getXXVar(const Index i1, const Index i2)
{
checkXDim(i1);
checkXDim(i2);
computeVarMat();
return data_.getXXVar(i1, i2);
}
const DMat & XYSampleData::getYYVar(const Index j1, const Index j2)
{
checkYDim(j1);
checkYDim(j2);
computeVarMat();
return data_.getYYVar(j1, j2);
}
const DMat & XYSampleData::getXYVar(const Index i, const Index j)
{
checkXDim(i);
checkYDim(j);
computeVarMat();
return data_.getXYVar(i, j);
}
DVec XYSampleData::getXError(const Index i)
{
checkXDim(i);
computeVarMat();
return data_.getXError(i);
}
DVec XYSampleData::getYError(const Index j)
{
checkYDim(j);
computeVarMat();
return data_.getYError(j);
}
// get total fit variance matrix and its pseudo-inverse ////////////////////////
const DMat & XYSampleData::getFitVarMat(void)
{
computeVarMat();
return data_.getFitVarMat();
}
const DMat & XYSampleData::getFitVarMatPInv(void)
{
computeVarMat();
return data_.getFitVarMatPInv();
}
// set data to a particular sample /////////////////////////////////////////////
void XYSampleData::setDataToSample(const Index s)
{
if (initData_ or (s != dataSample_))
{
for (Index i = 0; i < getNXDim(); ++i)
for (Index r = 0; r < getXSize(i); ++r)
{
data_.x(r, i) = xData_[i][r][s];
}
for (Index j = 0; j < getNYDim(); ++j)
for (auto &p: yData_[j])
{
data_.y(p.first, j) = p.second[s];
}
dataSample_ = s;
initData_ = false;
}
}
// get internal XYStatData /////////////////////////////////////////////////////
const XYStatData & XYSampleData::getData(void)
{
setDataToSample(central);
return data_;
}
void XYSampleData::resize(const Index nData, const Index xDim,
const Index yDim, const Index nSample)
{
FitInterface::resize(nData, xDim, yDim);
x_.resize(nSample);
x_.resizeMat(nData, xDim);
y_.resize(nSample);
y_.resizeMat(nData, yDim);
data_.resize(nData, xDim, yDim);
isCovarianceInit_ = false;
}
XYSampleData::SampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
{
isCovarianceInit_ = false;
return x_.block(0, 0, getNData(), getXDim());
}
XYSampleData::ConstSampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
const
{
return x_.block(0, 0, getNData(), getXDim());
}
XYSampleData::SampleBlock XYSampleData::x(const Index i,
const PlaceHolder ph2 __dumb)
{
isCovarianceInit_ = false;
return x_.block(0, i, getNData(), 1);
}
XYSampleData::ConstSampleBlock XYSampleData::x(const Index i,
const PlaceHolder ph2 __dumb)
const
{
return x_.block(0, i, getNData(), 1);
}
XYSampleData::SampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
const Index k)
{
isCovarianceInit_ = false;
return x_.block(k, 0, 1, getXDim());
}
XYSampleData::ConstSampleBlock XYSampleData::x(const PlaceHolder ph1 __dumb,
const Index k) const
{
return x_.block(k, 0, 1, getXDim());
}
XYSampleData::SampleBlock XYSampleData::x(const Index i, const Index k)
{
isCovarianceInit_ = false;
return x_.block(k, i, 1, 1);
}
XYSampleData::ConstSampleBlock XYSampleData::x(const Index i, const Index k)
const
{
return x_.block(k, i, 1, 1);
}
XYSampleData::SampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
{
isCovarianceInit_ = false;
return y_.block(0, 0, getNData(), getYDim());
}
XYSampleData::ConstSampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
const
{
return y_.block(0, 0, getNData(), getYDim());
}
XYSampleData::SampleBlock XYSampleData::y(const Index j,
const PlaceHolder ph2 __dumb)
{
isCovarianceInit_ = false;
return y_.block(0, j, getNData(), 1);
}
XYSampleData::ConstSampleBlock XYSampleData::y(const Index j,
const PlaceHolder ph2 __dumb)
const
{
return y_.block(0, j, getNData(), 1);
}
XYSampleData::SampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
const Index k)
{
isCovarianceInit_ = false;
return y_.block(k, 0, 1, getYDim());
}
XYSampleData::ConstSampleBlock XYSampleData::y(const PlaceHolder ph1 __dumb,
const Index k) const
{
return y_.block(k, 0, 1, getYDim());
}
XYSampleData::SampleBlock XYSampleData::y(const Index j, const Index k)
{
isCovarianceInit_ = false;
return y_.block(k, j, 1, 1);
}
XYSampleData::ConstSampleBlock XYSampleData::y(const Index j, const Index k)
const
{
return y_.block(k, j, 1, 1);
}
// fit /////////////////////////////////////////////////////////////////////////
SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &modelVector)
SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
const DVec &init,
const std::vector<const DoubleModel *> &v)
{
const Index nSample = x_.size();
FitResult sampleResult;
SampleFitResult result;
DVec initBuf = init;
computeVarMat();
result.resize(nSample);
result.chi2_.resize(nSample);
FOR_STAT_ARRAY(x_, s)
SampleFitResult result;
FitResult sampleResult;
DVec initCopy = init;
result.resize(nSample_);
result.chi2_.resize(nSample_);
result.model_.resize(v.size());
FOR_STAT_ARRAY(result, s)
{
// reinit chi^2 for central value only
setDataToSample(s);
if (s == central)
{
data_.reinitChi2(true);
sampleResult = data_.fit(minimizer, initCopy, v);
}
else
{
data_.reinitChi2(false);
sampleResult = data_.fit(*(minimizer.back()), initCopy, v);
}
// set data
setDataToSample(s);
// fit
sampleResult = data_.fit(minimizer, initBuf, modelVector);
if (s == central)
{
initBuf = sampleResult;
}
// store result
initCopy = sampleResult.segment(0, initCopy.size());
result[s] = sampleResult;
result.chi2_[s] = sampleResult.getChi2();
result.nDof_ = sampleResult.getNDof();
result.model_.resize(modelVector.size());
for (unsigned int j = 0; j < modelVector.size(); ++j)
for (unsigned int j = 0; j < v.size(); ++j)
{
result.model_[j].resize(nSample);
result.model_[j].resize(nSample_);
result.model_[j][s] = sampleResult.getModel(j);
}
}
result.nPar_ = sampleResult.getNPar();
result.nDof_ = sampleResult.nDof_;
result.parName_ = sampleResult.parName_;
return result;
}
void XYSampleData::setDataToSample(const Index s)
SampleFitResult XYSampleData::fit(Minimizer &minimizer,
const DVec &init,
const std::vector<const DoubleModel *> &v)
{
// compute covariance matrices if necessary
if (!isCovarianceInit_)
{
DMatSample buf1, buf2;
for (Index i2 = 0; i2 < getXDim(); ++i2)
for (Index i1 = 0; i1 < getXDim(); ++i1)
{
buf1 = x(i1);
buf2 = x(i2);
data_.xxVar(i1, i2) = buf1.covarianceMatrix(buf2);
}
for (Index j2 = 0; j2 < getYDim(); ++j2)
for (Index j1 = 0; j1 < getYDim(); ++j1)
{
buf1 = y(j1);
buf2 = y(j2);
data_.yyVar(j1, j2) = buf1.covarianceMatrix(buf2);
}
for (Index i = 0; i < getXDim(); ++i)
for (Index j = 0; j < getYDim(); ++j)
{
buf1 = y(j);
buf2 = x(i);
data_.yxVar(j, i) = buf1.covarianceMatrix(buf2);
}
isCovarianceInit_ = true;
}
vector<Minimizer *> mv{&minimizer};
// copy interface to sample data
data_.setFitInterface(*this);
// set data
data_.x() = x_[s];
data_.y() = y_[s];
return fit(mv, init, v);
}
// residuals ///////////////////////////////////////////////////////////////////
XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit) const
XYSampleData XYSampleData::getResiduals(const SampleFitResult &fit)
{
const Index nSample = x_.size();
XYSampleData res(*this);
DMatSample xBuf(nSample, getXDim(), 1), tmp(nSample, 1, 1);
for (Index j = 0; j < res.getYDim(); ++j)
for (Index j = 0; j < getNYDim(); ++j)
{
const DoubleFunctionSample &f = fit.getModel(_, j);
for (Index k = 0; k < res.getNData(); ++k)
for (auto &p: yData_[j])
{
xBuf = this->x(_, k);
tmp = this->y(j, k);
FOR_STAT_ARRAY(xBuf, s)
{
tmp[s](0) -= f[s](xBuf[s].transpose());
}
res.y(j, k) = tmp;
res.y(p.first, j) -= f(x(p.first));
}
}
return res;
}
XYSampleData XYSampleData::getPartialResiduals(const SampleFitResult &fit,
const DVec &x,
const Index i) const
const DVec &ref, const Index i)
{
const Index nSample = x_.size();
XYSampleData res(*this);
DMatSample xBuf(nSample, getXDim(), 1), tmp(nSample, 1, 1);
DVec buf(x);
for (Index j = 0; j < res.getYDim(); ++j)
DMatSample buf(nSample_);
buf.fill(ref);
for (Index j = 0; j < getNYDim(); ++j)
{
const DoubleFunctionSample &f = fit.getModel(_, j);
for (Index k = 0; k < res.getNData(); ++k)
for (auto &p: yData_[j])
{
xBuf = this->x(_, k);
tmp = this->y(j, k);
FOR_STAT_ARRAY(xBuf, s)
FOR_STAT_ARRAY(buf, s)
{
buf(i) = xBuf[s](i);
tmp[s](0) -= f[s](xBuf[s].transpose()) - f[s](buf);
buf[s](i) = x(p.first)[s](i);
}
res.y(j, k) = tmp;
res.y(p.first, j) -= f(x(p.first)) - f(buf);
}
}
return res;
}
// buffer list of x vectors ////////////////////////////////////////////////////
void XYSampleData::scheduleXMapInit(void)
{
initXMap_ = true;
}
void XYSampleData::updateXMap(void)
{
if (initXMap_)
{
for (Index s = central; s < nSample_; ++s)
{
setDataToSample(s);
for (auto k: getDataIndexSet())
{
if (s == central)
{
xMap_[k].resize(nSample_);
}
xMap_[k][s] = data_.x(k);
}
}
initXMap_ = false;
}
}
// schedule data initilization from samples ////////////////////////////////////
void XYSampleData::scheduleDataInit(void)
{
initData_ = true;
}
// variance matrix computation /////////////////////////////////////////////////
void XYSampleData::scheduleComputeVarMat(void)
{
computeVarMat_ = true;
}
void XYSampleData::computeVarMat(void)
{
if (computeVarMat_)
{
// initialize data if necessary
setDataToSample(central);
// compute relevant sizes
Index size = 0, ySize = 0;
for (Index j = 0; j < getNYDim(); ++j)
{
size += getYSize(j);
}
ySize = size;
for (Index i = 0; i < getNXDim(); ++i)
{
size += getXSize(i);
}
// compute total matrix
DMatSample z(nSample_, size, 1);
DMat var;
Index a;
FOR_STAT_ARRAY(z, s)
{
a = 0;
for (Index j = 0; j < getNYDim(); ++j)
for (auto &p: yData_[j])
{
z[s](a, 0) = p.second[s];
a++;
}
for (Index i = 0; i < getNXDim(); ++i)
for (Index r = 0; r < getXSize(i); ++r)
{
z[s](a, 0) = xData_[i][r][s];
a++;
}
}
var = z.varianceMatrix();
// assign blocks to data
Index a1, a2;
a1 = ySize;
for (Index i1 = 0; i1 < getNXDim(); ++i1)
{
a2 = ySize;
for (Index i2 = 0; i2 < getNXDim(); ++i2)
{
data_.setXXVar(i1, i2,
var.block(a1, a2, getXSize(i1), getXSize(i2)));
a2 += getXSize(i2);
}
a1 += getXSize(i1);
}
a1 = 0;
for (Index j1 = 0; j1 < getNYDim(); ++j1)
{
a2 = 0;
for (Index j2 = 0; j2 < getNYDim(); ++j2)
{
data_.setYYVar(j1, j2,
var.block(a1, a2, getYSize(j1), getYSize(j2)));
a2 += getYSize(j2);
}
a1 += getYSize(j1);
}
a1 = ySize;
for (Index i = 0; i < getNXDim(); ++i)
{
a2 = 0;
for (Index j = 0; j < getNYDim(); ++j)
{
data_.setXYVar(i, j,
var.block(a1, a2, getXSize(i), getYSize(j)));
a2 += getYSize(j);
}
a1 += getXSize(i);
}
computeVarMat_ = false;
}
if (initVarMat())
{
data_.copyInterface(*this);
scheduleFitVarMatInit(false);
}
}
// create data /////////////////////////////////////////////////////////////////
void XYSampleData::createXData(const string name, const Index nData)
{
data_.addXDim(nData, name);
xData_.push_back(vector<DSample>(nData));
}
void XYSampleData::createYData(const string name)
{
data_.addYDim(name);
yData_.push_back(map<Index, DSample>());
}

154
lib/XYSampleData.hpp
View File

@ -1,7 +1,7 @@
/*
* XYSampleData.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -22,9 +22,9 @@
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/FitInterface.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/Model.hpp>
#include <LatAnalyze/XYStatData.hpp>
BEGIN_LATAN_NAMESPACE
@ -47,90 +47,138 @@ public:
double getChi2PerDof(const Index s = central) const;
DSample getChi2PerDof(const PlaceHolder ph) const;
double getNDof(void) const;
Index getNPar(void) const;
double getPValue(const Index s = central) const;
const DoubleFunction & getModel(const Index s = central,
const Index j = 0) const;
const DoubleFunctionSample & getModel(const PlaceHolder ph,
const Index j = 0) const;
FitResult getFitResult(const Index s = central) const;
// IO
void print(const bool printXsi = false,
std::ostream &out = std::cout) const;
private:
DSample chi2_;
double nDof_{0.};
Index nDof_{0}, nPar_{0};
std::vector<DoubleFunctionSample> model_;
std::vector<std::string> parName_;
};
/******************************************************************************
* XYSampleData *
******************************************************************************
* index convention: i: X, j: Y, k: data
*/
* XYSampleData *
******************************************************************************/
class XYSampleData: public FitInterface
{
public:
typedef DMatSample::Block SampleBlock;
typedef DMatSample::ConstBlock ConstSampleBlock;
public:
// constructors
XYSampleData(void) = default;
XYSampleData(const Index nData, const Index xDim, const Index yDim,
const Index nSample);
// constructor
explicit XYSampleData(const Index nSample);
// destructor
virtual ~XYSampleData(void) = default;
// access
const XYStatData & getData(const Index s = central);
void resize(const Index nData, const Index xDim,
const Index yDim, const Index nSample);
SampleBlock x(const PlaceHolder ph1 = _, const PlaceHolder ph2 = _);
ConstSampleBlock x(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _) const;
SampleBlock x(const Index i, const PlaceHolder ph2 = _);
ConstSampleBlock x(const Index i, const PlaceHolder ph2 = _) const;
SampleBlock x(const PlaceHolder ph1, const Index k);
ConstSampleBlock x(const PlaceHolder ph1, const Index k) const;
SampleBlock x(const Index i, const Index k);
ConstSampleBlock x(const Index i, const Index k) const;
SampleBlock y(const PlaceHolder ph1 = _, const PlaceHolder ph2 = _);
ConstSampleBlock y(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _) const;
SampleBlock y(const Index j, const PlaceHolder ph2 = _);
ConstSampleBlock y(const Index j, const PlaceHolder ph2 = _) const;
SampleBlock y(const PlaceHolder ph1, const Index k);
ConstSampleBlock y(const PlaceHolder ph1, const Index k) const;
SampleBlock y(const Index j, const Index k);
ConstSampleBlock y(const Index j, const Index k) const;
// fit
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &modelVector);
template <typename... Mods>
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Mods... models);
// residuals
XYSampleData getResiduals(const SampleFitResult &fit) const;
XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x,
const Index j) const;
private:
// data access
DSample & x(const Index r, const Index i);
const DSample & x(const Index r, const Index i) const;
const DMatSample & x(const Index k);
DSample & y(const Index k, const Index j);
const DSample & y(const Index k, const Index j) const;
template <typename... Ts>
void setUnidimData(const DMatSample &xData,
const Ts & ...yDatas);
const DMat & getXXVar(const Index i1, const Index i2);
const DMat & getYYVar(const Index j1, const Index j2);
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
const DMat & getFitVarMat(void);
const DMat & getFitVarMatPInv(void);
// set data to a particular sample
void setDataToSample(const Index s);
// get internal XYStatData
const XYStatData & getData(void);
// fit
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &v);
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &v);
template <typename... Ts>
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models);
template <typename... Ts>
SampleFitResult fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models);
// residuals
XYSampleData getResiduals(const SampleFitResult &fit);
XYSampleData getPartialResiduals(const SampleFitResult &fit, const DVec &x,
const Index i);
private:
bool isCovarianceInit_;
DMatSample x_, y_;
XYStatData data_;
// buffer list of x vectors
void scheduleXMapInit(void);
void updateXMap(void);
// schedule data initilization from samples
void scheduleDataInit(void);
// variance matrix computation
void scheduleComputeVarMat(void);
void computeVarMat(void);
// create data
virtual void createXData(const std::string name, const Index nData);
virtual void createYData(const std::string name);
private:
std::vector<std::map<Index, DSample>> yData_;
std::vector<std::vector<DSample>> xData_;
std::map<Index, DMatSample> xMap_;
XYStatData data_;
Index nSample_, dataSample_{central};
bool initData_{true}, computeVarMat_{true};
bool initXMap_{true};
};
/******************************************************************************
* XYSampleData template implementation *
******************************************************************************/
template <typename... Ts>
SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
void XYSampleData::setUnidimData(const DMatSample &xData, const Ts & ...yDatas)
{
static_assert(static_or<std::is_assignable<DMatSample, Ts>::value...>::value,
"y data arguments are not compatible with DMatSample");
std::vector<const DMatSample *> yData{&yDatas...};
FOR_STAT_ARRAY(xData, s)
FOR_VEC(xData[central], r)
{
x(r, 0)[s] = xData[s](r);
for (unsigned int j = 0; j < yData.size(); ++j)
{
y(r, j)[s] = (*(yData[j]))[s](r);
}
}
}
template <typename... Ts>
SampleFitResult XYSampleData::fit(std::vector<Minimizer *> &minimizer,
const DVec &init,
const DoubleModel &model, const Ts... models)
{
static_assert(static_or<std::is_assignable<DoubleModel &, Ts>::value...>::value,
"model arguments are not compatible with DoubleModel &");
"model arguments are not compatible with DoubleModel");
std::vector<const DoubleModel *> modelVector{&model, &models...};
return fit(minimizer, init, modelVector);
}
template <typename... Ts>
SampleFitResult XYSampleData::fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models)
{
static_assert(static_or<std::is_assignable<DoubleModel &, Ts>::value...>::value,
"model arguments are not compatible with DoubleModel");
std::vector<Minimizer *> mv{&minimizer};
return fit(mv, init, model, models...);
}
END_LATAN_NAMESPACE
#endif // Latan_XYSampleData_hpp_

671
lib/XYStatData.cpp
View File

@ -1,7 +1,7 @@
/*
* XYStatData.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -23,7 +23,8 @@
using namespace std;
using namespace Latan;
using namespace Math;
static constexpr double maxXsiDev = 10.;
/******************************************************************************
* FitResult implementation *
@ -44,272 +45,566 @@ double FitResult::getNDof(void) const
return static_cast<double>(nDof_);
}
Index FitResult::getNPar(void) const
{
return nPar_;
}
double FitResult::getPValue(void) const
{
return chi2PValue(getChi2(), getNDof());;
return Math::chi2PValue(getChi2(), getNDof());;
}
const DoubleFunction & FitResult::getModel(const Index j) const
{
return model_[static_cast<unsigned int>(j)];
return model_[j];
}
// IO //////////////////////////////////////////////////////////////////////////
void FitResult::print(const bool printXsi, ostream &out) const
{
char buf[256];
Index pMax = printXsi ? size() : nPar_;
sprintf(buf, "chi^2/dof= %.1f/%d= %.2f -- p-value= %.2e", getChi2(),
static_cast<int>(getNDof()), getChi2PerDof(), getPValue());
out << buf << endl;
for (Index p = 0; p < pMax; ++p)
{
sprintf(buf, "%8s= %e", parName_[p].c_str(), (*this)(p));
out << buf << endl;
}
}
/******************************************************************************
* XYStatData implementation *
* XYStatData implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
XYStatData::XYStatData(void)
: chi2_(*this)
{}
XYStatData::XYStatData(const Index nData, const Index xDim, const Index yDim)
: XYStatData()
// data access /////////////////////////////////////////////////////////////////
double & XYStatData::x(const Index r, const Index i)
{
resize(nData, xDim, yDim);
checkXIndex(r, i);
scheduleXMapInit();
scheduleChi2DataVecInit();
return xData_[i](r);
}
// access //////////////////////////////////////////////////////////////////////
void XYStatData::resize(const Index nData, const Index xDim, const Index yDim)
const double & XYStatData::x(const Index r, const Index i) const
{
FitInterface::resize(nData, xDim, yDim);
x_.resize(nData, xDim);
y_.resize(nData, yDim);
var_[xx].resize(xDim, xDim);
var_[yy].resize(yDim, yDim);
var_[yx].resize(yDim, xDim);
FOR_MAT(var_[xx], i1, i2)
checkXIndex(r, i);
return xData_[i](r);
}
const DVec & XYStatData::x(const Index k) const
{
checkDataIndex(k);
updateXMap();
return xMap_.at(k);
}
double & XYStatData::y(const Index k, const Index j)
{
checkYDim(j);
if (!pointExists(k, j))
{
var_[xx](i1, i2).resize(nData, nData);
registerDataPoint(k, j);
resizeVarMat();
}
FOR_MAT(var_[yy], j1, j2)
scheduleXMapInit();
scheduleChi2DataVecInit();
return yData_[j][k];
}
const double & XYStatData::y(const Index k, const Index j) const
{
checkPoint(k, j);
return yData_[j].at(k);
}
void XYStatData::setXXVar(const Index i1, const Index i2, const DMat &m)
{
checkXDim(i1);
checkXDim(i2);
checkVarMat(m, xxVar_(i1, i2));
xxVar_(i1, i2) = m;
if (i1 != i2)
{
var_[yy](j1, j2).resize(nData, nData);
xxVar_(i2, i1) = m.transpose();
}
FOR_MAT(var_[yx], j, i)
scheduleFitVarMatInit();
}
void XYStatData::setYYVar(const Index j1, const Index j2, const DMat &m)
{
checkYDim(j1);
checkYDim(j2);
checkVarMat(m, yyVar_(j1, j2));
yyVar_(j1, j2) = m;
if (j1 != j2)
{
var_[yx](j, i).resize(nData, nData);
yyVar_(j2, j1) = m.transpose();
}
scheduleFitVarMatInit();
}
void XYStatData::reinitChi2(const bool doReinit)
void XYStatData::setXYVar(const Index i, const Index j, const DMat &m)
{
reinitChi2_ = doReinit;
checkXDim(i);
checkYDim(j);
checkVarMat(m, xyVar_(i, j));
xyVar_(i, j) = m;
scheduleFitVarMatInit();
}
Block<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
void XYStatData::setXError(const Index i, const DVec &err)
{
return x_.block(0, 0, getNData(), getXDim());
checkXDim(i);
checkErrVec(err, xxVar_(i, i));
xxVar_(i, i).diagonal() = err.cwiseProduct(err);
scheduleFitVarMatInit();
}
ConstBlock<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
const
void XYStatData::setYError(const Index j, const DVec &err)
{
return x_.block(0, 0, getNData(), getXDim());
checkXDim(j);
checkErrVec(err, yyVar_(j, j));
yyVar_(j, j).diagonal() = err.cwiseProduct(err);
scheduleFitVarMatInit();
}
Block<MatBase<double>> XYStatData::x(const Index i,
const PlaceHolder ph2 __dumb)
const DMat & XYStatData::getXXVar(const Index i1, const Index i2) const
{
return x_.block(0, i, getNData(), 1);
checkXDim(i1);
checkXDim(i2);
return xxVar_(i1, i2);
}
ConstBlock<MatBase<double>> XYStatData::x(const Index i,
const PlaceHolder ph2 __dumb)
const
const DMat & XYStatData::getYYVar(const Index j1, const Index j2) const
{
return x_.block(0, i, getNData(), 1);
checkYDim(j1);
checkYDim(j2);
return yyVar_(j1, j2);
}
Block<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
const Index k)
const DMat & XYStatData::getXYVar(const Index i, const Index j) const
{
return x_.block(k, 0, 1, getXDim());
checkXDim(i);
checkYDim(j);
return xyVar_(i, j);
}
ConstBlock<MatBase<double>> XYStatData::x(const PlaceHolder ph1 __dumb,
const Index k) const
DVec XYStatData::getXError(const Index i) const
{
return x_.block(k, 0, 1, getXDim());
checkXDim(i);
return xxVar_(i, i).diagonal().cwiseSqrt();
}
double & XYStatData::x(const Index i, const Index k)
DVec XYStatData::getYError(const Index j) const
{
return x_(k, i);
checkXDim(j);
return yyVar_(j, j).diagonal().cwiseSqrt();
}
const double & XYStatData::x(const Index i, const Index k) const
DMat XYStatData::getTable(const Index i, const Index j) const
{
return x_(k, i);
checkXDim(i);
checkYDim(j);
DMat table(getYSize(j), 4);
Index row = 0;
for (auto &p: yData_[j])
{
Index k = p.first;
Index r = dataCoord(k)[i];
table(row, 0) = x(k)(i);
table(row, 2) = p.second;
table(row, 1) = xxVar_(i, i).diagonal().cwiseSqrt()(r);
table(row, 3) = yyVar_(j, j).diagonal().cwiseSqrt()(row);
row++;
}
return table;
}
Block<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
// get total fit variance matrix ///////////////////////////////////////////////
const DMat & XYStatData::getFitVarMat(void)
{
return y_.block(0, 0, getNData(), getYDim());
updateFitVarMat();
return fitVar_;
}
ConstBlock<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb,
const PlaceHolder ph2 __dumb)
const
const DMat & XYStatData::getFitVarMatPInv(void)
{
return y_.block(0, 0, getNData(), getYDim());
}
Block<MatBase<double>> XYStatData::y(const Index j,
const PlaceHolder ph2 __dumb)
{
return y_.block(0, j, getNData(), 1);
}
ConstBlock<MatBase<double>> XYStatData::y(const Index j,
const PlaceHolder ph2 __dumb)
const
{
return y_.block(0, j, getNData(), 1);
}
Block<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb, const Index k)
{
return y_.block(k, 0, 1, getYDim());
}
ConstBlock<MatBase<double>> XYStatData::y(const PlaceHolder ph1 __dumb,
const Index k) const
{
return y_.block(k, 0, 1, getYDim());
}
double & XYStatData::y(const Index j, const Index k)
{
return y_(k, j);
}
const double & XYStatData::y(const Index j, const Index k) const
{
return y_(k, j);
}
#define FULL_BLOCK(m) (m).block(0, 0, (m).rows(), (m).cols())
Block<MatBase<double>> XYStatData::xxVar(const Index i1, const Index i2)
{
return FULL_BLOCK(var_[xx](i1, i2));
}
ConstBlock<MatBase<double>> XYStatData::xxVar(const Index i1,
const Index i2) const
{
return FULL_BLOCK(var_[xx](i1, i2));
}
Block<MatBase<double>> XYStatData::yyVar(const Index j1, const Index j2)
{
return FULL_BLOCK(var_[yy](j1, j2));
}
ConstBlock<MatBase<double>> XYStatData::yyVar(const Index j1,
const Index j2) const
{
return FULL_BLOCK(var_[yy](j1, j2));
}
Block<MatBase<double>> XYStatData::yxVar(const Index j, const Index i)
{
return FULL_BLOCK(var_[yx](j, i));
}
ConstBlock<MatBase<double>> XYStatData::yxVar(const Index j,
const Index i) const
{
return FULL_BLOCK(var_[yx](j, i));
updateFitVarMat();
return fitVarInv_;
}
// fit /////////////////////////////////////////////////////////////////////////
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
const vector<const DoubleModel *> &modelVector)
FitResult XYStatData::fit(vector<Minimizer *> &minimizer, const DVec &init,
const vector<const DoubleModel *> &v)
{
// initialization
chi2_.setModel(modelVector);
if (reinitChi2_)
// check model consistency
checkModelVec(v);
// buffering
updateLayout();
updateFitVarMat();
updateChi2DataVec();
// get number of parameters
Index nPar = v[0]->getNPar();
Index nXDim = getNXDim();
Index totalNPar = nPar + layout.totalXSize;
// chi^2 functions
auto corrChi2Func = [this, nPar, nXDim, totalNPar, &v](const double *x)->double
{
chi2_.requestInit();
ConstMap<DVec> p(x, totalNPar);
updateChi2ModVec(p, v, nPar, nXDim);
chi2Vec_ = (chi2ModVec_ - chi2DataVec_);
return chi2Vec_.dot(fitVarInv_*chi2Vec_);
};
DoubleFunction corrChi2(corrChi2Func, totalNPar);
auto uncorrChi2Func = [this, nPar, nXDim, totalNPar, &v](const double *x)->double
{
ConstMap<DVec> p(x, totalNPar);
updateChi2ModVec(p, v, nPar, nXDim);
chi2Vec_ = (chi2ModVec_ - chi2DataVec_);
return chi2Vec_.dot(chi2Vec_.cwiseQuotient(fitVar_.diagonal()));
};
DoubleFunction uncorrChi2(uncorrChi2Func, totalNPar);
DoubleFunction &chi2 = hasCorrelations() ? corrChi2 : uncorrChi2;
for (Index p = 0; p < nPar; ++p)
{
chi2.varName().setName(p, v[0]->parName().getName(p));
}
// initial parameters
const Index nPoint = getNFitPoint();
DVec fullInit = init;
Index is = 0, kf = 0;
fullInit.conservativeResize(chi2_.getNArg());
for (Index i = 0; i < getXDim(); ++i)
for (Index p = 0; p < totalNPar - nPar; ++p)
{
if (!isXExact(i))
chi2.varName().setName(p + nPar, "xsi_" + strFrom(p));
}
// minimization
FitResult result;
DVec totalInit(totalNPar);
//// set total init vector
totalInit.segment(0, nPar) = init;
totalInit.segment(nPar, layout.totalXSize) =
chi2DataVec_.segment(layout.totalYSize, layout.totalXSize);
for (auto &m: minimizer)
{
m->setInit(totalInit);
if (m->supportLimits())
{
kf = 0;
for (Index k = 0; k < getNData(); ++k)
//// do not allow more than maxXsiDev std. deviations on the x-axis
for (Index p = nPar; p < totalNPar; ++p)
{
if (isFitPoint(k))
{
fullInit(chi2_.getNPar() + nPoint*is + kf) = x(i, k);
kf++;
}
double err;
err = sqrt(fitVar_.diagonal()(layout.totalYSize + p - nPar));
m->useLowLimit(p);
m->useHighLimit(p);
m->setLowLimit(p, totalInit(p) - maxXsiDev*err);
m->setHighLimit(p, totalInit(p) + maxXsiDev*err);
}
is++;
}
//// minimize and store results
result = (*m)(chi2);
totalInit = result;
}
minimizer.setInit(fullInit);
// fit
DoubleFunction chi2 = chi2_.makeFunction(false);
FitResult result;
result = minimizer(chi2);
result.chi2_ = chi2(result);
result.nDof_ = chi2_.getNDof();
result.model_.resize(modelVector.size());
for (unsigned int j = 0; j < modelVector.size(); ++j)
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)
{
result.model_[j] = modelVector[j]->fixPar(result);
result.model_[j] = v[j]->fixPar(result);
}
for (Index p = 0; p < totalNPar; ++p)
{
result.parName_.push_back(chi2.varName().getName(p));
}
return result;
}
// residuals ///////////////////////////////////////////////////////////////////
XYStatData XYStatData::getResiduals(const FitResult &fit) const
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
const vector<const DoubleModel *> &v)
{
XYStatData res(*this);
for (Index j = 0; j < res.getYDim(); ++j)
{
const DoubleFunction &f = fit.getModel(j);
for (Index k = 0; k < res.getNData(); ++k)
{
res.y(j, k) -= f(res.x(_, k).transpose());
}
}
return res;
vector<Minimizer *> mv{&minimizer};
return fit(mv, init, v);
}
XYStatData XYStatData::getPartialResiduals(const FitResult &fit, const DVec &x,
const Index i) const
// residuals ///////////////////////////////////////////////////////////////////
XYStatData XYStatData::getResiduals(const FitResult &fit)
{
XYStatData res(*this);
DVec buf(x), xk;
for (Index j = 0; j < res.getYDim(); ++j)
for (Index j = 0; j < getNYDim(); ++j)
{
const DoubleFunction &f = fit.getModel(j);
for (Index k = 0; k < res.getNData(); ++k)
for (auto &p: yData_[j])
{
buf(i) = res.x(i, k);
res.y(j, k) -= f(res.x(_, k).transpose()) - f(buf);
res.y(p.first, j) -= f(x(p.first));
}
}
return res;
}
XYStatData XYStatData::getPartialResiduals(const FitResult &fit,
const DVec &ref, const Index i)
{
XYStatData res(*this);
DVec buf(ref);
for (Index j = 0; j < res.getNYDim(); ++j)
{
const DoubleFunction &f = fit.getModel(j);
for (auto &p: yData_[j])
{
buf(i) = x(p.first)(i);
res.y(p.first, j) -= f(x(p.first)) - f(buf);
}
}
return res;
}
// create data /////////////////////////////////////////////////////////////////
void XYStatData::createXData(const std::string name __dumb, const Index nData)
{
xData_.push_back(DVec::Zero(nData));
xBuf_.resize(xData_.size());
resizeVarMat();
}
void XYStatData::createYData(const std::string name __dumb)
{
yData_.push_back(map<Index, double>());
resizeVarMat();
}
void XYStatData::resizeVarMat(void)
{
xxVar_.conservativeResize(getNXDim(), getNXDim());
for (Index i1 = 0; i1 < getNXDim(); ++i1)
for (Index i2 = 0; i2 < getNXDim(); ++i2)
{
xxVar_(i1, i2).conservativeResize(getXSize(i1), getXSize(i2));
}
yyVar_.conservativeResize(getNYDim(), getNYDim());
for (Index j1 = 0; j1 < getNYDim(); ++j1)
for (Index j2 = 0; j2 < getNYDim(); ++j2)
{
yyVar_(j1, j2).conservativeResize(getYSize(j1), getYSize(j2));
}
xyVar_.conservativeResize(getNXDim(), getNYDim());
for (Index i = 0; i < getNXDim(); ++i)
for (Index j = 0; j < getNYDim(); ++j)
{
xyVar_(i, j).conservativeResize(getXSize(i), getYSize(j));
}
scheduleFitVarMatInit();
}
// schedule buffer computation /////////////////////////////////////////////////
void XYStatData::scheduleXMapInit(void)
{
initXMap_ = true;
}
void XYStatData::scheduleChi2DataVecInit(void)
{
initChi2DataVec_ = true;
}
// buffer total fit variance matrix ////////////////////////////////////////////
void XYStatData::updateFitVarMat(void)
{
if (initVarMat())
{
updateLayout();
DMat &v = fitVar_;
Index roffs, coffs;
v.resize(layout.totalSize, layout.totalSize);
roffs = layout.totalYSize;
for (Index ifit1 = 0; ifit1 < layout.nXFitDim; ++ifit1)
{
coffs = layout.totalYSize;
for (Index ifit2 = 0; ifit2 < layout.nXFitDim; ++ifit2)
{
for (Index rfit1 = 0; rfit1 < layout.xSize[ifit1]; ++rfit1)
for (Index rfit2 = 0; rfit2 < layout.xSize[ifit2]; ++rfit2)
{
Index i1, i2, r1, r2;
i1 = layout.xDim[ifit1];
i2 = layout.xDim[ifit2];
r1 = layout.x[ifit1][rfit1];
r2 = layout.x[ifit2][rfit2];
v(roffs+rfit1, coffs+rfit2) = xxVar_(i1, i2)(r1, r2);
v(coffs+rfit2, roffs+rfit1) = v(roffs+rfit1, coffs+rfit2);
}
coffs += layout.xSize[ifit2];
}
roffs += layout.xSize[ifit1];
}
roffs = 0;
for (Index jfit1 = 0; jfit1 < layout.nYFitDim; ++jfit1)
{
coffs = 0;
for (Index jfit2 = 0; jfit2 < layout.nYFitDim; ++jfit2)
{
for (Index sfit1 = 0; sfit1 < layout.ySize[jfit1]; ++sfit1)
for (Index sfit2 = 0; sfit2 < layout.ySize[jfit2]; ++sfit2)
{
Index j1, j2, s1, s2;
j1 = layout.yDim[jfit1];
j2 = layout.yDim[jfit2];
s1 = layout.y[jfit1][sfit1];
s2 = layout.y[jfit2][sfit2];
v(roffs+sfit1, coffs+sfit2) = yyVar_(j1, j2)(s1, s2);
v(coffs+sfit2, roffs+sfit1) = v(roffs+sfit1, coffs+sfit2);
}
coffs += layout.ySize[jfit2];
}
roffs += layout.ySize[jfit1];
}
roffs = layout.totalYSize;
for (Index ifit = 0; ifit < layout.nXFitDim; ++ifit)
{
coffs = 0;
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
{
for (Index rfit = 0; rfit < layout.xSize[ifit]; ++rfit)
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
{
Index i, j, r, s;
i = layout.xDim[ifit];
j = layout.yDim[jfit];
r = layout.x[ifit][rfit];
s = layout.y[jfit][sfit];
v(roffs+rfit, coffs+sfit) = xyVar_(i, j)(r, s);
v(coffs+sfit, roffs+rfit) = v(roffs+rfit, coffs+sfit);
}
coffs += layout.ySize[jfit];
}
roffs += layout.xSize[ifit];
}
chi2DataVec_.resize(layout.totalSize);
chi2ModVec_.resize(layout.totalSize);
chi2Vec_.resize(layout.totalSize);
fitVar_ = fitVar_.cwiseProduct(makeCorrFilter());
fitVarInv_ = fitVar_.pInverse(getSvdTolerance());
scheduleFitVarMatInit(false);
}
}
// buffer list of x vectors ////////////////////////////////////////////////////
void XYStatData::updateXMap(void) const
{
if (initXMap_)
{
XYStatData * modThis = const_cast<XYStatData *>(this);
modThis->xMap_.clear();
modThis->xMap_.resize(getMaxDataIndex());
for (auto k: getDataIndexSet())
{
modThis->xMap_[k] = DVec(getNXDim());
for (Index i = 0; i < getNXDim(); ++i)
{
modThis->xMap_[k](i) = xData_[i](dataCoord(k)[i]);
}
}
modThis->initXMap_ = false;
}
}
// buffer chi^2 vectors ////////////////////////////////////////////////////////
void XYStatData::updateChi2DataVec(void)
{
if (initChi2DataVec_)
{
Index a = 0, j, k, i, r;
updateLayout();
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
{
j = layout.yDim[jfit];
k = layout.data[jfit][sfit];
chi2DataVec_(a) = yData_[j][k];
a++;
}
for (Index ifit = 0; ifit < layout.nXFitDim; ++ifit)
for (Index rfit = 0; rfit < layout.xSize[ifit]; ++rfit)
{
i = layout.xDim[ifit];
r = layout.x[ifit][rfit];
chi2DataVec_(a) = xData_[i](r);
a++;
}
initChi2DataVec_ = false;
}
}
// WARNING: updateChi2ModVec is heavily called by fit
void XYStatData::updateChi2ModVec(const DVec p,
const vector<const DoubleModel *> &v,
const Index nPar, const Index nXDim)
{
updateLayout();
updateXMap();
Index a = 0, j, k, ind;
auto &par = p.segment(0, nPar), &xsi = p.segment(nPar, layout.totalXSize);
for (Index jfit = 0; jfit < layout.nYFitDim; ++jfit)
{
j = layout.yDim[jfit];
for (Index sfit = 0; sfit < layout.ySize[jfit]; ++sfit)
{
k = layout.data[jfit][sfit];
for (Index i = 0; i < nXDim; ++i)
{
ind = layout.xIndFromData[k][i] - layout.totalYSize;
xBuf_(i) = (ind >= 0) ? xsi(ind) : xMap_[k](i);
}
chi2ModVec_(a) = (*v[j])(xBuf_.data(), par.data());
a++;
}
}
chi2ModVec_.segment(a, layout.totalXSize) = xsi;
}

228
lib/XYStatData.hpp
View File

@ -1,7 +1,7 @@
/*
* XYData.hpp, part of LatAnalyze 3
* XYStatData.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -17,17 +17,13 @@
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_XYData_hpp_
#define Latan_XYData_hpp_
#ifndef Latan_XYStatData_hpp_
#define Latan_XYStatData_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Chi2Function.hpp>
#include <LatAnalyze/FitInterface.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <LatAnalyze/Model.hpp>
#include <vector>
BEGIN_LATAN_NAMESPACE
@ -37,6 +33,7 @@ BEGIN_LATAN_NAMESPACE
class FitResult: public DVec
{
friend class XYStatData;
friend class XYSampleData;
friend class SampleFitResult;
public:
// constructors
@ -48,103 +45,192 @@ public:
double getChi2(void) const;
double getChi2PerDof(void) const;
double getNDof(void) const;
Index getNPar(void) const;
double getPValue(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.0};
Index nDof_{0};
double chi2_{0.};
Index nDof_{0}, nPar_{0};
std::vector<DoubleFunction> model_;
std::vector<std::string> parName_;
};
/******************************************************************************
* object for X vs. Y statistical data *
******************************************************************************
* index convention: i: X, j: Y, k: data
*/
* class for X vs. Y statistical data *
******************************************************************************/
class XYStatData: public FitInterface
{
public:
enum
{
xx = 0,
yy = 1,
yx = 2
};
public:
// constructors
XYStatData(void);
XYStatData(const Index nData, const Index nXDim, const Index nYDim);
// constructor
XYStatData(void) = default;
// destructor
virtual ~XYStatData(void) = default;
// access
void resize(const Index nData, const Index xDim,
const Index yDim);
void reinitChi2(const bool doReinit = true);
Block<MatBase<double>> x(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _);
ConstBlock<MatBase<double>> x(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _) const;
Block<MatBase<double>> x(const Index i, const PlaceHolder ph2 = _);
ConstBlock<MatBase<double>> x(const Index i,
const PlaceHolder ph2 = _) const;
Block<MatBase<double>> x(const PlaceHolder ph1, const Index k);
ConstBlock<MatBase<double>> x(const PlaceHolder ph1,
const Index k) const;
double & x(const Index i, const Index k);
const double & x(const Index i, const Index k) const;
Block<MatBase<double>> y(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _);
ConstBlock<MatBase<double>> y(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _) const;
Block<MatBase<double>> y(const Index i, const PlaceHolder ph2 = _);
ConstBlock<MatBase<double>> y(const Index i,
const PlaceHolder ph2 = _) const;
Block<MatBase<double>> y(const PlaceHolder ph1, const Index k);
ConstBlock<MatBase<double>> y(const PlaceHolder ph1,
const Index k) const;
double & y(const Index i, const Index k);
const double & y(const Index i, const Index k) const;
Block<MatBase<double>> xxVar(const Index i1, const Index i2);
ConstBlock<MatBase<double>> xxVar(const Index i1, const Index i2) const;
Block<MatBase<double>> yyVar(const Index j1, const Index j2);
ConstBlock<MatBase<double>> yyVar(const Index j1, const Index j2) const;
Block<MatBase<double>> yxVar(const Index j, const Index i);
ConstBlock<MatBase<double>> yxVar(const Index j, const Index i) const;
// data access
double & x(const Index r, const Index i);
const double & x(const Index r, const Index i) const;
const DVec & x(const Index k) const;
double & y(const Index k, const Index j);
const double & y(const Index k, const Index j) const;
void setXXVar(const Index i1, const Index i2, const DMat &m);
void setYYVar(const Index j1, const Index j2, const DMat &m);
void setXYVar(const Index i, const Index j, const DMat &m);
void setXError(const Index i, const DVec &err);
void setYError(const Index j, const DVec &err);
template <typename... Ts>
void setUnidimData(const DMat &xData, const Ts & ...yDatas);
const DMat & getXXVar(const Index i1, const Index i2) const;
const DMat & getYYVar(const Index j1, const Index j2) const;
const DMat & getXYVar(const Index i, const Index j) const;
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
const DMat & getFitVarMat(void);
const DMat & getFitVarMatPInv(void);
// fit
FitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &v);
FitResult fit(Minimizer &minimizer, const DVec &init,
const std::vector<const DoubleModel *> &modelVector);
const std::vector<const DoubleModel *> &v);
template <typename... Ts>
FitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models);
template <typename... Ts>
FitResult fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models);
// residuals
XYStatData getResiduals(const FitResult &fit) const;
XYStatData getPartialResiduals(const FitResult &fit, const DVec &x,
const Index j) const;
XYStatData getResiduals(const FitResult &fit);
XYStatData getPartialResiduals(const FitResult &fit, const DVec &ref,
const Index i);
protected:
// create data
virtual void createXData(const std::string name, const Index nData);
virtual void createYData(const std::string name);
void resizeVarMat(void);
private:
DMat x_, y_;
Mat<DMat> var_[3];
IVec isXExact_, isFitPoint_;
Chi2Function chi2_;
bool reinitChi2_{true};
// schedule buffer computation
void scheduleXMapInit(void);
void scheduleChi2DataVecInit(void);
// buffer total fit variance matrix
void updateFitVarMat(void);
// buffer list of x vectors
void updateXMap(void) const;
// buffer chi^2 vectors
void updateChi2DataVec(void);
void updateChi2ModVec(const DVec p,
const std::vector<const DoubleModel *> &v,
const Index nPar, const Index nXDim);
private:
std::vector<std::map<Index, double>> yData_;
// no map here for fit performance
std::vector<DVec> xData_;
std::vector<DVec> xMap_;
Mat<DMat> xxVar_, yyVar_, xyVar_;
DMat fitVar_, fitVarInv_;
DVec chi2DataVec_, chi2ModVec_, chi2Vec_;
DVec xBuf_;
bool initXMap_{true};
bool initChi2DataVec_{true};
};
/******************************************************************************
* XYStatData template implementation *
******************************************************************************/
template <typename... Ts>
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
void XYStatData::setUnidimData(const DMat &xData, const Ts & ...yDatas)
{
static_assert(static_or<std::is_assignable<DMat, Ts>::value...>::value,
"y data arguments are not compatible with DMat");
std::vector<const DMat *> yData{&yDatas...};
FOR_VEC(xData, r)
{
x(r, 0) = xData(r);
for (unsigned int j = 0; j < yData.size(); ++j)
{
y(r, j) = (*(yData[j]))(r);
}
}
}
template <typename... Ts>
FitResult XYStatData::fit(std::vector<Minimizer *> &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models)
{
static_assert(static_or<std::is_assignable<DoubleModel &, Ts>::value...>::value,
"model arguments are not compatible with DoubleModel &");
"model arguments are not compatible with DoubleModel");
std::vector<const DoubleModel *> modelVector{&model, &models...};
return fit(minimizer, init, modelVector);
}
template <typename... Ts>
FitResult XYStatData::fit(Minimizer &minimizer, const DVec &init,
const DoubleModel &model, const Ts... models)
{
static_assert(static_or<std::is_assignable<DoubleModel &, Ts>::value...>::value,
"model arguments are not compatible with DoubleModel");
std::vector<Minimizer *> mv{&minimizer};
return fit(mv, init, model, models...);
}
/******************************************************************************
* error check macros *
******************************************************************************/
#define checkVarMat(m, var)\
if (((m).rows() != (var).rows()) or ((m).cols() != (var).cols()))\
{\
LATAN_ERROR(Size, "provided variance matrix has a wrong size"\
" (expected " + strFrom((var).rows()) + "x"\
+ strFrom((var).cols()) + ", got " + strFrom((m).rows())\
+ "x" + strFrom((m).cols()) + ")");\
}
#define checkErrVec(err, var)\
if ((err).size() != (var).rows())\
{\
LATAN_ERROR(Size, "provided error vector has a wrong size"\
" (expected " + strFrom((var).rows()) + ", got "\
+ strFrom((err).size()) + ")");\
}
#define checkModelVec(v)\
if (static_cast<Index>((v).size()) != getNYDim())\
{\
LATAN_ERROR(Size, "provided model vector has a wrong size"\
" (expected " + strFrom(getNYDim()) + ", got "\
+ strFrom((v).size()) + ")");\
}\
for (unsigned int _i = 1; _i < (v).size(); ++_i)\
{\
if ((v)[_i]->getNArg() != getNXDim())\
{\
LATAN_ERROR(Size, "model " + strFrom(_i) + " has a wrong"\
+ " number of argument (expected " + strFrom(getNXDim())\
+ ", got " + strFrom((v)[_i]->getNArg()));\
}\
}\
{\
Index _nPar = (v)[0]->getNPar();\
for (unsigned int _i = 1; _i < (v).size(); ++_i)\
{\
if ((v)[_i]->getNPar() != _nPar)\
{\
LATAN_ERROR(Size, "model " + strFrom(_i) + " has a wrong"\
+ " number of parameter (expected " + strFrom(_nPar)\
+ ", got " + strFrom((v)[_i]->getNPar()));\
}\
}\
}
END_LATAN_NAMESPACE
#endif // Latan_XYData_hpp_
#endif // Latan_XYStatData_hpp_

19
lib/includes.hpp
View File

@ -1,7 +1,7 @@
/*
* includes.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -20,23 +20,6 @@
#ifndef Latan_includes_hpp_
#define Latan_includes_hpp_
#include <algorithm>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <limits>
#include <sstream>
#include <utility>
#include <cfloat>
#include <climits>
#include <cmath>
#include <cstdarg>
#include <cstdio>
#include <cstdlib>
#include <sys/stat.h>
#include <unistd.h>
#include <config.h>
#endif // Latan_includes_hpp_

22
utils/Makefile.am
View File

@ -1,30 +1,18 @@
if CC_GNU
COM_CFLAGS = -Wall -W -pedantic
else
if CC_INTEL
COM_CFLAGS = -Wall
endif
endif
if CXX_GNU
COM_CXXFLAGS = -Wall -W -pedantic
COM_CXXFLAGS = -Wall -W -pedantic -Wno-deprecated-declarations
else
if CXX_INTEL
COM_CXXFLAGS = -Wall
COM_CXXFLAGS = -wd1682 -Wall
endif
endif
bin_PROGRAMS = \
latan_create_rg_state \
latan_make_fake_sample\
latan_sample_combine \
latan_sample_plot_corr\
latan_sample_read \
latan_resample
latan_create_rg_state_SOURCES = create_rg_state.cpp
latan_create_rg_state_CXXFLAGS = $(COM_CXXFLAGS)
latan_create_rg_state_LDFLAGS = -L../lib/.libs -lLatAnalyze
latan_make_fake_sample_SOURCES = make_fake_sample.cpp
latan_make_fake_sample_CXXFLAGS = $(COM_CXXFLAGS)
latan_make_fake_sample_LDFLAGS = -L../lib/.libs -lLatAnalyze
@ -33,6 +21,10 @@ latan_sample_combine_SOURCES = sample_combine.cpp
latan_sample_combine_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_combine_LDFLAGS = -L../lib/.libs -lLatAnalyze
latan_sample_plot_corr_SOURCES = sample_plot_corr.cpp
latan_sample_plot_corr_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_plot_corr_LDFLAGS = -L../lib/.libs -lLatAnalyze
latan_sample_read_SOURCES = sample_read.cpp
latan_sample_read_CXXFLAGS = $(COM_CXXFLAGS)
latan_sample_read_LDFLAGS = -L../lib/.libs -lLatAnalyze

16
utils/make_fake_sample.cpp
View File

@ -1,7 +1,7 @@
/*
* make_fake_sample.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -17,9 +17,7 @@
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/RandGen.hpp>
using namespace std;
using namespace Latan;
@ -42,21 +40,23 @@ int main(int argc, char *argv[])
nSample = strTo<Index>(argv[3]);
outFileName = argv[4];
RandGen gen;
DMatSample res(nSample, 1, 1);
random_device rd;
mt19937 gen(rd());
normal_distribution<> dis(val, err);
DSample res(nSample);
FOR_STAT_ARRAY(res, s)
{
if (s == central)
{
res[s](0, 0) = val;
res[s] = val;
}
else
{
res[s](0, 0) = gen.gaussian(val, err);
res[s] = dis(gen);
}
}
Io::save<DMatSample>(res, outFileName);
Io::save<DSample>(res, outFileName);
return EXIT_SUCCESS;
}

32
utils/resample.cpp
View File

@ -1,7 +1,7 @@
/*
* resample.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -17,10 +17,6 @@
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <map>
#include <libgen.h>
#include <unistd.h>
#include <LatAnalyze/Dataset.hpp>
#include <LatAnalyze/Io.hpp>
@ -34,7 +30,7 @@ using namespace Latan;
static void usage(const string &cmdName)
{
cerr << "usage: " << cmdName;
cerr << " [-n <nsample> -b <bin size> -r <state> -o <output dir> -f {h5|sample}]";
cerr << " [-n <nsample> -b <bin size> -r <seed> -o <output dir> -f {h5|sample}]";
cerr << " <data list> <name list>";
cerr << endl;
exit(EXIT_FAILURE);
@ -43,10 +39,12 @@ static void usage(const string &cmdName)
int main(int argc, char *argv[])
{
// argument parsing ////////////////////////////////////////////////////////
int c;
string manFileName, nameFileName, stateFileName, cmdName, outDirName = ".";
string ext = "h5";
Index binSize = 1, nSample = DEF_NSAMPLE;
int c;
random_device rd;
SeedType seed = rd();
string manFileName, nameFileName, cmdName, outDirName = ".";
string ext = "h5";
Index binSize = 1, nSample = DEF_NSAMPLE;
opterr = 0;
cmdName = basename(argv[0]);
@ -64,7 +62,7 @@ int main(int argc, char *argv[])
outDirName = optarg;
break;
case 'r':
stateFileName = optarg;
seed = strTo<SeedType>(optarg);
break;
case 'f':
ext = optarg;
@ -126,25 +124,15 @@ int main(int argc, char *argv[])
// data resampling /////////////////////////////////////////////////////////
DMatSample s(nSample);
RandGen g;
RandGenState state;
cout << "-- resampling data..." << endl;
if (!stateFileName.empty())
{
state = Io::load<RandGenState>(stateFileName);
}
for (unsigned int i = 0; i < name.size(); ++i)
{
const string outFileName = name[i] + "_" + manFileName + "." + ext;
cout << '\r' << ProgressBar(i + 1, name.size());
data[name[i]].bin(binSize);
if (!stateFileName.empty())
{
g.setState(state);
}
s = data[name[i]].bootstrapMean(nSample, g);
s = data[name[i]].bootstrapMean(nSample, seed);
Io::save<DMatSample>(s, outDirName + "/" + outFileName,
File::Mode::write, outFileName);
}

159
utils/sample_combine.cpp
View File

@ -1,7 +1,7 @@
/*
* sample_combine.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2016 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
@ -17,9 +17,6 @@
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <libgen.h>
#include <unistd.h>
#include <LatAnalyze/Io.hpp>
#include <LatAnalyze/CompiledFunction.hpp>
@ -35,6 +32,108 @@ static void usage(const string &cmdName)
exit(EXIT_FAILURE);
}
template <typename T>
static void loadAndCheck(vector<T> &sample, const vector<string> &fileName)
{
const unsigned int n = sample.size();
Index nSample = 0;
cout << "-- loading data..." << endl;
for (unsigned int i = 0; i < n; ++i)
{
sample[i] = Io::load<T>(fileName[i]);
if (i == 0)
{
nSample = sample[i].size();
}
else if (sample[i].size() != nSample)
{
cerr << "error: number of sample mismatch (between '";
cerr << fileName[0] << "' and '" << fileName[i] << "')" << endl;
abort();
}
}
}
template <typename T>
static void combine(const string &outFileName __dumb,
const vector<T> &sample __dumb, const string &code __dumb)
{
abort();
}
template <>
void combine(const string &outFileName, const vector<DSample> &sample,
const string &code)
{
const unsigned int n = sample.size();
DoubleFunction f = compile(code, n);
DSample result(sample[0]);
DVec buf(n);
cout << "-- combining data..." << endl;
result = sample[0];
FOR_STAT_ARRAY(result, s)
{
for (unsigned int k = 0; k < n; ++k)
{
buf[k] = sample[k][s];
}
result[s] = f(buf);
}
cout << scientific;
cout << "central value:\n" << result[central];
cout << endl;
cout << "standard deviation:\n" << sqrt(result.variance());
cout << endl;
if (!outFileName.empty())
{
Io::save<DSample>(result, outFileName);
}
}
template <>
void combine(const string &outFileName, const vector<DMatSample> &sample,
const string &code)
{
const unsigned int n = sample.size();
DoubleFunction f = compile(code, n);
DVec buf(n);
DMatSample result(sample[0]);
cout << "-- combining data..." << endl;
FOR_STAT_ARRAY(result, s)
{
FOR_MAT(result[s], i, j)
{
for (unsigned int k = 0; k < n; ++k)
{
buf[k] = sample[k][s](i,j);
}
result[s](i, j) = f(buf);
}
}
cout << scientific;
cout << "central value:\n" << result[central];
cout << endl;
cout << "standard deviation:\n" << result.variance().cwiseSqrt();
cout << endl;
if (!outFileName.empty())
{
Io::save<DMatSample>(result, outFileName);
}
}
template <typename T>
void process(const string &outFileName, const vector<string> &fileName,
const string &code)
{
vector<T> sample(fileName.size());
loadAndCheck(sample, fileName);
combine(outFileName, sample, code);
}
int main(int argc, char *argv[])
{
// argument parsing ////////////////////////////////////////////////////////
@ -82,54 +181,16 @@ int main(int argc, char *argv[])
{
usage(cmdName);
}
// data loading ////////////////////////////////////////////////////////////
vector<DMatSample> sample(n);
Index nSample = 0;
cout << "-- loading data..." << endl;
for (unsigned int i = 0; i < n; ++i)
// process data ////////////////////////////////////////////////////////////
try
{
sample[i] = Io::load<DMatSample>(fileName[i]);
if (i == 0)
{
nSample = sample[i].size();
}
else if (sample[i].size() != nSample)
{
cerr << "error: number of sample mismatch (between '";
cerr << fileName[0] << "' and '" << fileName[i] << "')" << endl;
return EXIT_FAILURE;
}
process<DSample>(outFileName, fileName, code);
}
// combine data ////////////////////////////////////////////////////////////
DoubleFunction f = compile(code, n);
DVec buf(n);
DMatSample result(sample[0]);
cout << "-- combining data..." << endl;
FOR_STAT_ARRAY(result, s)
catch (bad_cast &e)
{
FOR_MAT(result[s], i, j)
{
for (unsigned int k = 0; k < n; ++k)
{
buf[k] = sample[k][s](i,j);
}
result[s](i, j) = f(buf);
}
process<DMatSample>(outFileName, fileName, code);
}
// output //////////////////////////////////////////////////////////////////
cout << scientific;
cout << "central value:\n" << result[central] << endl;
cout << "standard deviation:\n" << result.variance().cwiseSqrt() << endl;
if (!outFileName.empty())
{
Io::save<DMatSample>(result, outFileName);
}
return EXIT_SUCCESS;
}

Some files were not shown because too many files have changed in this diff Show More