portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2025-10-31 14:54:33 +00:00
Code Releases Activity

Merge branch 'release/v3.5.1'

Browse Source
This commit is contained in:
Antonin Portelli
2020-11-27 13:34:26 +00:00
parent f356700a4f cfb2c0c5e8
commit 448f34477b
104 changed files with 1012 additions and 428 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
Readme.md
View File

@@ -71,6 +71,9 @@ make install
where `<prefix>` should be replaced by the desired prefix for LatAnalyze installation, and `<n>` is the number of parallel build processes (typically twice your number of cores).
## History
#### v3.5.1
Various fixes and cleaning of outdated code.
#### v3.5
Additions:
* 'Impulse' & line type plots

2
bootstrap.sh
View File

@@ -2,5 +2,5 @@
rm -rf .buildutils
mkdir -p .buildutils/m4
./update_eigen.sh eigen-3.3.7.tar.bz2
./update_eigen.sh eigen-3.3.8.tar.bz2
autoreconf -fvi

7
ci-scripts/install-gsl.sh
View File

@@ -9,8 +9,11 @@ fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd local/build
INITDIR=$(pwd -P)
mkdir -p ${PREFIX}
cd ${PREFIX}
PREFIX=$(pwd -P)
cd ${INITDIR}/local/build
wget http://ftpmirror.gnu.org/gsl/${NAME}.tar.gz
tar -xzvf ${NAME}.tar.gz
mkdir -p ${NAME}/build

7
ci-scripts/install-hdf5.sh
View File

@@ -9,8 +9,11 @@ fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd local/build
INITDIR=$(pwd -P)
mkdir -p ${PREFIX}
cd ${PREFIX}
PREFIX=$(pwd -P)
cd ${INITDIR}/local/build
wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/${NAME}/src/${NAME}.tar.gz
tar -xzvf ${NAME}.tar.gz
mkdir ${NAME}/build

8
ci-scripts/install-latan.sh
View File

@@ -5,14 +5,18 @@ if (( $# != 1 )); then
exit 1
fi
PREFIX=$1
OS=$2
set -ex
INITDIR=$(pwd -P)
mkdir -p ${PREFIX}
cd ${PREFIX}
PREFIX=$(pwd -P)
cd ${INITDIR}
./install-deps.sh ${PREFIX}
cd ..
./bootstrap.sh
mkdir -p build
cd build
../configure --prefix=$PREFIX --with-minuit=$PREFIX --with-nlopt=$PREFIX --with-latcore=$PREFIX --with-hdf5=$PREFIX --with-gsl=$PREFIX CXXFLAGS="${CXXFLAGS} -O3 -march=native -mtune=native"
../configure --prefix=${PREFIX} --with-minuit=${PREFIX} --with-nlopt=${PREFIX} --with-hdf5=${PREFIX} --with-gsl=${PREFIX} CXXFLAGS="${CXXFLAGS} -O3 -march=haswell -mtune=haswell"
make -j4
make install

7
ci-scripts/install-minuit.sh
View File

@@ -9,8 +9,11 @@ fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd local/build
INITDIR=$(pwd -P)
mkdir -p ${PREFIX}
cd ${PREFIX}
PREFIX=$(pwd -P)
cd ${INITDIR}/local/build
wget http://www.cern.ch/mathlibs/sw/5_34_14/Minuit2/${NAME}.tar.gz
tar -xzvf ${NAME}.tar.gz
mkdir -p ${NAME}/build

9
ci-scripts/install-nlopt.sh
View File

@@ -9,14 +9,17 @@ fi
PREFIX=$1
set -ex
INITDIR=`pwd`
cd local/build
INITDIR=$(pwd -P)
mkdir -p ${PREFIX}
cd ${PREFIX}
PREFIX=$(pwd -P)
cd ${INITDIR}/local/build
wget https://github.com/stevengj/nlopt/archive/v${NAME}.tar.gz
tar -xzvf v${NAME}.tar.gz
NAME=nlopt-${NAME}
mkdir -p ${NAME}/build
cd ${NAME}/build
cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} ..
cmake -DCMAKE_INSTALL_PREFIX=${PREFIX} -DCMAKE_BUILD_WITH_INSTALL_NAME_DIR=TRUE -DCMAKE_INSTALL_NAME_DIR="${PREFIX}/lib" ..
make -j4
make install
cd ${INITDIR}/local

2
configure.ac
View File

@@ -2,7 +2,7 @@
# Initialization
AC_PREREQ([2.63])
AC_INIT([LatAnalyze],[3.5],[antonin.portelli@me.com],[LatAnalyze])
AC_INIT([LatAnalyze],[3.5.1],[antonin.portelli@me.com],[LatAnalyze])
AC_CONFIG_AUX_DIR([.buildutils])
AC_CONFIG_SRCDIR([lib/Global.cpp])
AC_CONFIG_SRCDIR([utils/sample_read.cpp])

BIN
eigen-3.3.7.tar.bz2
View File

Binary file not shown.

BIN
eigen-3.3.8.tar.bz2 Normal file
View File

Binary file not shown.

2
lib/Core/Eigen.hpp
View File

@@ -1,7 +1,7 @@
/*
* Eigen.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/Exceptions.cpp
View File

@@ -1,7 +1,7 @@
/*
* Exceptions.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/Exceptions.hpp
View File

@@ -1,7 +1,7 @@
/*
* Exceptions.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/Mat.cpp
View File

@@ -1,7 +1,7 @@
/*
* Mat.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/Mat.hpp
View File

@@ -1,7 +1,7 @@
/*
* Mat.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/Math.cpp
View File

@@ -1,7 +1,7 @@
/*
* Math.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/Math.hpp
View File

@@ -1,7 +1,7 @@
/*
* Math.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/MathInterpreter.cpp
View File

@@ -1,7 +1,7 @@
/*
* MathInterpreter.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/MathInterpreter.hpp
View File

@@ -1,7 +1,7 @@
/*
* MathInterpreter.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/MathLexer.lpp
View File

@@ -1,7 +1,7 @@
/*
* MathLexer.lpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/MathParser.ypp
View File

@@ -1,7 +1,7 @@
/*
* MathParser.ypp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/ParserState.hpp
View File

@@ -1,7 +1,7 @@
/*
* ParserState.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

152
lib/Core/Plot.cpp
View File

@@ -1,7 +1,7 @@
/*
* Plot.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -236,14 +236,21 @@ PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin,
// PlotPredBand constructor ////////////////////////////////////////////////////
void PlotPredBand::makePredBand(const DMat &low, const DMat &high, const double opacity)
{
string lowFileName, highFileName;
string lowFileName, highFileName, contFileName;
DMat contour(low.rows() + high.rows() + 1, 2);
lowFileName = dumpToTmpFile(low);
highFileName = dumpToTmpFile(high);
pushTmpFile(lowFileName);
pushTmpFile(highFileName);
setCommand("'< (cat " + lowFileName + "; tac " + highFileName +
"; head -n1 " + lowFileName + ")' u 1:2 w filledcurves closed" +
FOR_MAT(low, i, j)
{
contour(i, j) = low(i, j);
}
FOR_MAT(high, i, j)
{
contour(low.rows() + i, j) = high(high.rows() - i - 1, j);
}
contour.row(low.rows() + high.rows()) = low.row(0);
contFileName = dumpToTmpFile(contour);
pushTmpFile(contFileName);
setCommand("'" + contFileName + "' u 1:2 w filledcurves closed" +
" fs solid " + strFrom(opacity) + " noborder");
}
@@ -500,7 +507,7 @@ Plot::Plot(void)
// default options /////////////////////////////////////////////////////////////
void Plot::initOptions(void)
{
options_.terminal = "qt";
options_.terminal = "qt noenhanced font 'Arial,12'";
options_.output = "";
options_.caption = "";
options_.title = "";
@@ -580,57 +587,48 @@ Plot & Plot::operator<<(PlotModifier &&modifier)
}
// find gnuplot ////////////////////////////////////////////////////////////////
void Plot::getProgramPath(void)
#define SEARCH_DIR(dir) \
sprintf(buf, "%s/%s", dir, gnuplotBin_.c_str());\
if (access(buf, X_OK) == 0)\
{\
return dir;\
}
std::string Plot::getProgramPath(void)
{
int i, j, lg;
char *path;
static char buf[MAX_PATH_LENGTH];
// trivial case: try in CWD
sprintf(buf,"./%s", gnuplotBin_.c_str()) ;
if (access(buf, X_OK) == 0)
{
sprintf(buf,".");
gnuplotPath_ = buf;
}
// try out in all paths given in the PATH variable
else
buf[0] = 0;
path = getenv("PATH") ;
if (path)
{
buf[0] = 0;
path = getenv("PATH") ;
if (path)
for (i=0;path[i];)
{
for (i=0;path[i];)
for (j=i;(path[j]) and (path[j]!=':');j++);
lg = j - i;
strncpy(buf,path + i,(size_t)(lg));
if (lg == 0)
{
for (j=i;(path[j]) and (path[j]!=':');j++);
lg = j - i;
strncpy(buf,path + i,(size_t)(lg));
if (lg == 0)
{
buf[lg++] = '.';
}
buf[lg++] = '/';
strcpy(buf + lg, gnuplotBin_.c_str());
if (access(buf, X_OK) == 0)
{
// found it!
break ;
}
buf[0] = 0;
i = j;
if (path[i] == ':') i++ ;
buf[lg++] = '.';
}
buf[lg++] = '/';
strcpy(buf + lg, gnuplotBin_.c_str());
if (access(buf, X_OK) == 0)
{
// found it!
break ;
}
buf[0] = 0;
i = j;
if (path[i] == ':') i++ ;
}
else
{
LATAN_ERROR(System, "PATH variable not set");
}
// if the buffer is still empty, the command was not found
if (buf[0] == 0)
{
LATAN_ERROR(System, "cannot find gnuplot in $PATH");
}
// otherwise truncate the command name to yield path only
}
// if the buffer is still empty, the command was not found
if (buf[0] != 0)
{
lg = (int)(strlen(buf) - 1);
while (buf[lg]!='/')
{
@@ -639,7 +637,19 @@ void Plot::getProgramPath(void)
}
buf[lg] = 0;
gnuplotPath_ = buf;
return gnuplotPath_;
}
// try in CWD, /usr/bin & /usr/local/bin
SEARCH_DIR(".");
SEARCH_DIR("/usr/bin");
SEARCH_DIR("/usr/local/bin");
// if this code is reached nothing was found
LATAN_ERROR(System, "cannot find gnuplot");
return "";
}
// plot parsing and output /////////////////////////////////////////////////////
@@ -653,10 +663,6 @@ void Plot::display(void)
string command;
ostringstream scriptBuf;
if (!getenv("DISPLAY"))
{
LATAN_ERROR(System, "cannot find DISPLAY variable: is it set ?");
}
getProgramPath();
command = gnuplotPath_ + "/" + gnuplotBin_ + " 2>/dev/null";
gnuplotPipe = popen(command.c_str(), "w");
@@ -683,7 +689,7 @@ void Plot::display(void)
}
}
void Plot::save(string dirName)
void Plot::save(string dirName, bool savePdf)
{
vector<string> commandBack;
string path, terminalBack, outputBack, gpCommand, scriptName;
@@ -691,11 +697,6 @@ void Plot::save(string dirName)
ofstream script;
mode755 = S_IRWXU|S_IRGRP|S_IXGRP|S_IROTH|S_IXOTH;
// backup I/O parameters
terminalBack = options_.terminal;
outputBack = options_.output;
commandBack = plotCommand_;
// generate directory
if (mkdir(dirName))
@@ -703,12 +704,20 @@ void Plot::save(string dirName)
LATAN_ERROR(Io, "impossible to create directory '" + dirName + "'");
}
// backup I/O parameters
terminalBack = options_.terminal;
outputBack = options_.output;
commandBack = plotCommand_;
// save PDF
options_.terminal = "pdf";
options_.output = dirName + "/plot.pdf";
display();
options_.terminal = terminalBack;
options_.output = outputBack;
if (savePdf)
{
options_.terminal = "pdf";
options_.output = dirName + "/plot.pdf";
display();
options_.terminal = terminalBack;
options_.output = outputBack;
}
// save script and datafiles
for (unsigned int i = 0; i < tmpFileName_.size(); ++i)
@@ -776,18 +785,25 @@ ostream & Latan::operator<<(ostream &out, const Plot &plot)
out << "yMin = " << plot.options_.scale[y].min << endl;
out << "yMax = " << plot.options_.scale[y].max << endl;
}
if (!plot.options_.title.empty())
{
out << "set title '" << plot.options_.title << "'" << endl;
}
out << "unset xrange" << endl;
if (plot.options_.scaleMode[x] & Plot::Scale::manual)
{
out << "set xrange [xMin:xMax]" << endl;
}
else
{
out << "set xrange [:]" << endl;
}
out << "unset yrange" << endl;
if (plot.options_.scaleMode[y] & Plot::Scale::manual)
{
out << "set yrange [yMin:yMax]" << endl;
}
else
{
out << "set yrange [:]" << endl;
}
out << "unset log" << endl;
if (plot.options_.scaleMode[x] & Plot::Scale::log)
{
out << "set log x" << endl;

8
lib/Core/Plot.hpp
View File

@@ -1,7 +1,7 @@
/*
* Plot.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -377,13 +377,13 @@ public:
Plot & operator<<(PlotModifier &&modifier);
// plot parsing and output
void display(void);
void save(std::string dirName);
void save(std::string dirName, bool savePdf = true);
friend std::ostream & operator<<(std::ostream &out, const Plot &plot);
// plot reset
void reset(void);
private:
// find gnuplot
void getProgramPath(void);
std::string getProgramPath(void);
private:
// default options
void initOptions(void);
private:

2
lib/Core/Utilities.cpp
View File

@@ -1,7 +1,7 @@
/*
* Utilities.cpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/Utilities.hpp
View File

@@ -1,7 +1,7 @@
/*
* Utilities.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Core/stdincludes.hpp
View File

@@ -1,7 +1,7 @@
/*
* stdincludes.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/CompiledFunction.cpp
View File

@@ -1,7 +1,7 @@
/*
* CompiledFunction.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/CompiledFunction.hpp
View File

@@ -1,7 +1,7 @@
/*
* CompiledFunction.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/CompiledModel.cpp
View File

@@ -1,7 +1,7 @@
/*
* CompiledModel.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/CompiledModel.hpp
View File

@@ -1,7 +1,7 @@
/*
* CompiledModel.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/Function.cpp
View File

@@ -1,7 +1,7 @@
/*
* Function.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/Function.hpp
View File

@@ -1,7 +1,7 @@
/*
* Function.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/Model.cpp
View File

@@ -1,7 +1,7 @@
/*
* Model.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/Model.hpp
View File

@@ -1,7 +1,7 @@
/*
* Model.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/TabFunction.cpp
View File

@@ -1,7 +1,7 @@
/*
* TabFunction.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Functional/TabFunction.hpp
View File

@@ -1,7 +1,7 @@
/*
* TabFunction.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Global.cpp
View File

@@ -1,7 +1,7 @@
/*
* Global.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Global.hpp
View File

@@ -1,7 +1,7 @@
/*
* Global.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/AsciiFile.cpp
View File

@@ -1,7 +1,7 @@
/*
* AsciiFile.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/AsciiFile.hpp
View File

@@ -1,7 +1,7 @@
/*
* AsciiFile.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/AsciiLexer.lpp
View File

@@ -1,7 +1,7 @@
/*
* AsciiLexer.lpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/AsciiParser.ypp
View File

@@ -1,7 +1,7 @@
/*
* AsciiParser.ypp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/File.cpp
View File

@@ -1,7 +1,7 @@
/*
* File.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/File.hpp
View File

@@ -1,7 +1,7 @@
/*
* File.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/Hdf5File.cpp
View File

@@ -1,7 +1,7 @@
/*
* Hdf5File.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
* Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/Hdf5File.hpp
View File

@@ -1,7 +1,7 @@
/*
* Hdf5File.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
* Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/Io.cpp
View File

@@ -1,7 +1,7 @@
/*
* Io.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/Io.hpp
View File

@@ -1,7 +1,7 @@
/*
* Io.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/IoObject.hpp
View File

@@ -1,7 +1,7 @@
/*
* IoObject.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/XmlReader.cpp
View File

@@ -1,7 +1,7 @@
/*
* XmlReader.cpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2014 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Io/XmlReader.hpp
View File

@@ -1,7 +1,7 @@
/*
* XmlReader.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2015 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

4
lib/Makefile.am
View File

@@ -54,6 +54,8 @@ libLatAnalyze_la_SOURCES = \
Numerical/Minimizer.cpp \
Numerical/RootFinder.cpp \
Numerical/Solver.cpp \
Physics/CorrelatorFitter.cpp \
Physics/EffectiveMass.cpp \
Statistics/FitInterface.cpp \
Statistics/Histogram.cpp \
Statistics/Random.cpp \
@@ -97,6 +99,8 @@ HPPFILES = \
Numerical/Minimizer.hpp \
Numerical/RootFinder.hpp \
Numerical/Solver.hpp \
Physics/CorrelatorFitter.hpp \
Physics/EffectiveMass.hpp \
Statistics/Dataset.hpp \
Statistics/FitInterface.hpp \
Statistics/Histogram.hpp \

2
lib/Numerical/Derivative.cpp
View File

@@ -1,7 +1,7 @@
/*
* Derivative.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/Derivative.hpp
View File

@@ -1,7 +1,7 @@
/*
* Derivative.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/FFT.hpp
View File

@@ -1,7 +1,7 @@
/*
* FFT.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2017 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

3
lib/Numerical/GslFFT.cpp
View File

@@ -1,7 +1,7 @@
/*
* GslFFT.cpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2017 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -43,6 +43,7 @@ void GslFFT::resize(const Index size)
{
if (size_ != size)
{
clear();
size_ = size;
wavetable_ = gsl_fft_complex_wavetable_alloc(size_);
workspace_ = gsl_fft_complex_workspace_alloc(size_);

2
lib/Numerical/GslFFT.hpp
View File

@@ -1,7 +1,7 @@
/*
* GslFFT.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2017 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/GslHybridRootFinder.cpp
View File

@@ -1,7 +1,7 @@
/*
* GslHybridRootFinder.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/GslHybridRootFinder.hpp
View File

@@ -1,7 +1,7 @@
/*
* GslHybridRootFinder.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/GslMinimizer.cpp
View File

@@ -1,7 +1,7 @@
/*
* GslMinimizer.cpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/GslMinimizer.hpp
View File

@@ -1,7 +1,7 @@
/*
* GslMinimizer.hpp, part of LatAnalyze
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/GslQagsIntegrator.cpp
View File

@@ -1,7 +1,7 @@
/*
* GslQagsIntegrator.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/GslQagsIntegrator.hpp
View File

@@ -1,7 +1,7 @@
/*
* GslQagsIntegrator.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/Integrator.hpp
View File

@@ -1,7 +1,7 @@
/*
* Integrator.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/Minimizer.cpp
View File

@@ -1,7 +1,7 @@
/*
* Minimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/Minimizer.hpp
View File

@@ -1,7 +1,7 @@
/*
* Minimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/MinuitMinimizer.cpp
View File

@@ -1,7 +1,7 @@
/*
* MinuitMinimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/MinuitMinimizer.hpp
View File

@@ -1,7 +1,7 @@
/*
* MinuitMinimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/NloptMinimizer.cpp
View File

@@ -1,7 +1,7 @@
/*
* NloptMinimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/NloptMinimizer.hpp
View File

@@ -1,7 +1,7 @@
/*
* NloptMinimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/RootFinder.cpp
View File

@@ -1,7 +1,7 @@
/*
* RootFinder.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/RootFinder.hpp
View File

@@ -1,7 +1,7 @@
/*
* RootFinder.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/Solver.cpp
View File

@@ -1,7 +1,7 @@
/*
* Solver.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Numerical/Solver.hpp
View File

@@ -1,7 +1,7 @@
/*
* Solver.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

417
lib/Physics/CorrelatorFitter.cpp Normal file
View File

@@ -0,0 +1,417 @@
/*
* CorrelatorFitter.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* Correlator models *
******************************************************************************/
DoubleModel CorrelatorModels::makeExpModel(const Index nState)
{
DoubleModel mod;
mod.setFunction([nState](const double *x, const double *p)
{
double res = 0.;
for (unsigned int i = 0; i < nState; ++i)
{
res += p[2*i + 1]*exp(-p[2*i]*x[0]);
}
return res;
}, 1, 2*nState);
for (unsigned int i = 0; i < nState; ++i)
{
mod.parName().setName(2*i, "E_" + strFrom(i));
mod.parName().setName(2*i + 1, "Z_" + strFrom(i));
}
return mod;
}
DoubleModel CorrelatorModels::makeCoshModel(const Index nState, const Index nt)
{
DoubleModel mod;
mod.setFunction([nState, nt](const double *x, const double *p)
{
double res = 0.;
for (unsigned int i = 0; i < nState; ++i)
{
res += p[2*i + 1]*(exp(-p[2*i]*x[0]) + exp(-p[2*i]*(nt - x[0])));
}
return res;
}, 1, 2*nState);
for (unsigned int i = 0; i < nState; ++i)
{
mod.parName().setName(2*i, "E_" + strFrom(i));
mod.parName().setName(2*i + 1, "Z_" + strFrom(i));
}
return mod;
}
DoubleModel CorrelatorModels::makeSinhModel(const Index nState, const Index nt)
{
DoubleModel mod;
mod.setFunction([nState, nt](const double *x, const double *p)
{
double res = 0.;
for (unsigned int i = 0; i < nState; ++i)
{
res += p[2*i + 1]*(exp(-p[2*i]*x[0]) - exp(-p[2*i]*(nt - x[0])));
}
return res;
}, 1, 2*nState);
for (unsigned int i = 0; i < nState; ++i)
{
mod.parName().setName(2*i, "E_" + strFrom(i));
mod.parName().setName(2*i + 1, "Z_" + strFrom(i));
}
return mod;
}
DoubleModel CorrelatorModels::makeConstModel(void)
{
DoubleModel mod;
mod.setFunction([](const double *x __dumb, const double *p __dumb)
{
return p[0];
}, 1, 1);
mod.parName().setName(0, "cst");
return mod;
}
DoubleModel CorrelatorModels::makeLinearModel(void)
{
DoubleModel mod;
mod.setFunction([](const double *x, const double *p)
{
return p[1] + p[0]*x[0];
}, 1, 2);
return mod;
}
CorrelatorModels::ModelPar CorrelatorModels::parseModel(const string s)
{
smatch sm;
ModelPar par;
if (regex_match(s, sm, regex("exp([0-9]+)")))
{
par.type = CorrelatorType::exp;
par.nState = strTo<Index>(sm[1].str());
}
else if (regex_match(s, sm, regex("cosh([0-9]+)")))
{
par.type = CorrelatorType::cosh;
par.nState = strTo<Index>(sm[1].str());
}
else if (regex_match(s, sm, regex("sinh([0-9]+)")))
{
par.type = CorrelatorType::sinh;
par.nState = strTo<Index>(sm[1].str());
}
else if (s == "linear")
{
par.type = CorrelatorType::linear;
par.nState = 1;
}
else if (s == "cst")
{
par.type = CorrelatorType::cst;
par.nState = 1;
}
else
{
par.type = CorrelatorType::undefined;
par.nState = 0;
}
return par;
}
DoubleModel CorrelatorModels::makeModel(const CorrelatorModels::ModelPar par,
const Index nt)
{
switch (par.type)
{
case CorrelatorType::undefined:
LATAN_ERROR(Argument, "correlator type is undefined");
break;
case CorrelatorType::exp:
return makeExpModel(par.nState);
break;
case CorrelatorType::cosh:
return makeCoshModel(par.nState, nt);
break;
case CorrelatorType::sinh:
return makeSinhModel(par.nState, nt);
break;
case CorrelatorType::linear:
return makeLinearModel();
break;
case CorrelatorType::cst:
return makeConstModel();
break;
}
}
DVec CorrelatorModels::parameterGuess(const DMatSample &corr,
const ModelPar par)
{
DVec init;
Index nt = corr[central].size();
switch (par.type)
{
case CorrelatorType::undefined:
LATAN_ERROR(Argument, "correlator type is undefined");
break;
case CorrelatorType::exp:
case CorrelatorType::cosh:
case CorrelatorType::sinh:
init.resize(2*par.nState);
init(0) = log(corr[central](nt/4)/corr[central](nt/4 + 1));
init(1) = corr[central](nt/4)/(exp(-init(0)*nt/4));
for (Index p = 2; p < init.size(); p += 2)
{
init(p) = 2*init(p - 2);
init(p + 1) = init(p - 1)/2.;
}
break;
case CorrelatorType::linear:
init.resize(2);
init(0) = corr[central](nt/4) - corr[central](nt/4 + 1, 0);
init(1) = corr[central](nt/4, 0) + nt/4*init(0);
break;
case CorrelatorType::cst:
init.resize(1);
init(0) = corr[central](nt/4);
break;
default:
break;
}
return init;
}
/******************************************************************************
* Correlator utilities *
******************************************************************************/
DMatSample CorrelatorUtils::shift(const DMatSample &c, const Index ts)
{
if (ts != 0)
{
const Index nt = c[central].rows();
DMatSample buf = c;
FOR_STAT_ARRAY(buf, s)
{
for (Index t = 0; t < nt; ++t)
{
buf[s]((t - ts + nt)%nt) = c[s](t);
}
}
return buf;
}
else
{
return c;
}
}
DMatSample CorrelatorUtils::fold(const DMatSample &c)
{
const Index nt = c[central].rows();
DMatSample buf = c;
FOR_STAT_ARRAY(buf, s)
{
for (Index t = 0; t < nt; ++t)
{
buf[s](t) = 0.5*(c[s](t) + c[s]((nt - t) % nt));
}
}
return buf;
}
DMatSample CorrelatorUtils::fourierTransform(const DMatSample &c, FFT &fft,
const unsigned int dir)
{
const Index nSample = c.size();
const Index nt = c[central].rows();
bool isComplex = (c[central].cols() > 1);
CMatSample buf(nSample, nt, 1);
DMatSample out(nSample, nt, 2);
fft.resize(nt);
FOR_STAT_ARRAY(buf, s)
{
buf[s].real() = c[s].col(0);
if (isComplex)
{
buf[s].imag() = c[s].col(1);
}
else
{
buf[s].imag() = DVec::Constant(nt, 0.);
}
fft(buf[s], dir);
out[s].col(0) = buf[s].real();
out[s].col(1) = buf[s].imag();
}
return out;
}
/******************************************************************************
* CorrelatorFitter implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
CorrelatorFitter::CorrelatorFitter(const DMatSample &corr)
{
setCorrelator(corr);
}
CorrelatorFitter::CorrelatorFitter(const std::vector<DMatSample> &corr)
{
setCorrelators(corr);
}
// access //////////////////////////////////////////////////////////////////////
XYSampleData & CorrelatorFitter::data(void)
{
return *data_;
}
void CorrelatorFitter::setCorrelator(const DMatSample &corr)
{
std::vector<DMatSample> vec;
vec.push_back(corr);
setCorrelators(vec);
}
void CorrelatorFitter::setCorrelators(const std::vector<DMatSample> &corr)
{
Index nSample = corr[0].size();
DMatSample tVec(nSample);
std::vector<const DMatSample *> ptVec;
nt_ = corr[0][central].rows();
tVec.fill(DVec::LinSpaced(nt_, 0, nt_ - 1));
for (auto &c: corr)
{
ptVec.push_back(&c);
}
data_.reset(new XYSampleData(corr[0].size()));
data_->addXDim(nt_, "t/a", true);
for (unsigned int i = 0; i < corr.size(); ++i)
{
data_->addYDim("C_" + strFrom(i) + "(t)");
}
data_->setUnidimData(tVec, ptVec);
model_.resize(corr.size());
range_.resize(corr.size(), make_pair(0, nt_ - 1));
thinning_.resize(corr.size(), 1);
}
void CorrelatorFitter::setModel(const DoubleModel &model, const Index i)
{
model_[i] = model;
}
const DoubleModel & CorrelatorFitter::getModel(const Index i) const
{
return model_.at(i);
}
void CorrelatorFitter::setFitRange(const Index tMin, const Index tMax,
const Index i)
{
range_[i] = make_pair(tMin, tMax);
refreshRanges();
}
void CorrelatorFitter::setCorrelation(const bool isCorrelated, const Index i,
const Index j)
{
data_->assumeYYCorrelated(isCorrelated, i, j);
}
DMat CorrelatorFitter::getVarianceMatrix(void) const
{
return data_->getFitVarMat();
}
void CorrelatorFitter::setThinning(const Index thinning, const Index i)
{
thinning_[i] = thinning;
refreshRanges();
}
// fit functions ///////////////////////////////////////////////////////////////
SampleFitResult CorrelatorFitter::fit(Minimizer &minimizer, const DVec &init)
{
vector<Minimizer *> vecPt = {&minimizer};
return fit(vecPt, init);
}
SampleFitResult CorrelatorFitter::fit(vector<Minimizer *> &minimizer,
const DVec &init)
{
vector<const DoubleModel *> vecPt(model_.size());
for (unsigned int i = 0; i < model_.size(); ++i)
{
vecPt[i] = &(model_[i]);
}
return data_->fit(minimizer, init, vecPt);
}
// internal function to refresh fit ranges /////////////////////////////////////
void CorrelatorFitter::refreshRanges(void)
{
for (unsigned int i = 0; i < range_.size(); ++i)
for (Index t = 0; t < nt_; ++t)
{
data_->fitPoint((t >= range_[i].first) and (t <= range_[i].second)
and ((t - range_[i].first) % thinning_[i] == 0), t);
}
}

104
lib/Physics/CorrelatorFitter.hpp Normal file
View File

@@ -0,0 +1,104 @@
/*
* CorrelatorFitter.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_CorrelatorFitter_hpp_
#define Latan_CorrelatorFitter_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Functional/Model.hpp>
#include <LatAnalyze/Numerical/FFT.hpp>
#include <LatAnalyze/Statistics/XYSampleData.hpp>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* Correlator types & models *
******************************************************************************/
enum class CorrelatorType {undefined, exp, cosh, sinh, linear, cst};
namespace CorrelatorModels
{
struct ModelPar
{
CorrelatorType type;
Index nState;
};
DoubleModel makeExpModel(const Index nState);
DoubleModel makeCoshModel(const Index nState, const Index nt);
DoubleModel makeSinhModel(const Index nState, const Index nt);
DoubleModel makeConstModel(void);
DoubleModel makeLinearModel(void);
ModelPar parseModel(const std::string s);
DoubleModel makeModel(const ModelPar par, const Index nt);
DVec parameterGuess(const DMatSample &corr, const ModelPar par);
};
/******************************************************************************
* Correlator utilities *
******************************************************************************/
namespace CorrelatorUtils
{
DMatSample shift(const DMatSample &c, const Index ts);
DMatSample fold(const DMatSample &c);
DMatSample fourierTransform(const DMatSample &c, FFT &fft,
const unsigned int dir = FFT::Forward);
};
/******************************************************************************
* Correlator fit utility class *
******************************************************************************/
class CorrelatorFitter
{
public:
// constructors
CorrelatorFitter(const DMatSample &corr);
CorrelatorFitter(const std::vector<DMatSample> &corr);
// destructor
virtual ~CorrelatorFitter(void) = default;
// access
XYSampleData & data(void);
void setCorrelator(const DMatSample &corr);
void setCorrelators(const std::vector<DMatSample> &corr);
const DMatSample & getCorrelator(const Index i = 0) const;
const std::vector<DMatSample> & getCorrelators(void) const;
void setModel(const DoubleModel &model, const Index i = 0);
const DoubleModel & getModel(const Index i = 0) const;
void setFitRange(const Index tMin, const Index tMax, const Index i = 0);
void setCorrelation(const bool isCorrelated, const Index i = 0,
const Index j = 0);
DMat getVarianceMatrix(void) const;
void setThinning(const Index thinning, const Index i = 0);
// fit functions
SampleFitResult fit(Minimizer &minimizer, const DVec &init);
SampleFitResult fit(std::vector<Minimizer *> &minimizer, const DVec &init);
private:
// internal function to refresh fit ranges
void refreshRanges(void);
private:
Index nt_{0};
std::unique_ptr<XYSampleData> data_;
std::vector<DoubleModel> model_;
std::vector<std::pair<Index, Index>> range_;
std::vector<Index> thinning_;
};
END_LATAN_NAMESPACE
#endif // Latan_CorrelatorFitter_hpp_

132
lib/Physics/EffectiveMass.cpp Normal file
View File

@@ -0,0 +1,132 @@
/*
* EffectiveMass.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#include <LatAnalyze/Physics/EffectiveMass.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* EffectiveMass implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
EffectiveMass::EffectiveMass(const CorrelatorType type)
{
setType(type);
}
// access //////////////////////////////////////////////////////////////////////
CorrelatorType EffectiveMass::getType(void) const
{
return type_;
}
void EffectiveMass::setType(const CorrelatorType type)
{
type_ = type;
}
DVec EffectiveMass::getTime(const Index nt) const
{
DVec tvec;
switch (type_)
{
case CorrelatorType::undefined:
LATAN_ERROR(Argument, "correlator type is undefined");
break;
case CorrelatorType::exp:
case CorrelatorType::linear:
tvec = DVec::LinSpaced(nt - 1, 0, nt - 2);
break;
case CorrelatorType::cosh:
case CorrelatorType::sinh:
tvec = DVec::LinSpaced(nt - 2, 1, nt - 2);
break;
case CorrelatorType::cst:
tvec = DVec::LinSpaced(nt, 0, nt - 1);
break;
}
return tvec;
}
// compute effective mass //////////////////////////////////////////////////////
DVec EffectiveMass::operator()(const DVec &corr) const
{
Index nt = corr.size();
DVec em;
if (nt < 2)
{
LATAN_ERROR(Size, "input vector has less than 2 elements");
}
switch (type_)
{
case CorrelatorType::undefined:
LATAN_ERROR(Argument, "correlator type is undefined");
break;
case CorrelatorType::exp:
em.resize(nt - 1);
for (Index t = 1; t < nt; ++t)
{
em(t - 1) = log(corr(t - 1)/corr(t));
}
break;
case CorrelatorType::cosh:
em.resize(nt - 2);
for (Index t = 1; t < nt - 1; ++t)
{
em(t - 1) = acosh((corr(t - 1) + corr(t + 1))/(2.*corr(t)));
}
break;
case CorrelatorType::sinh:
em.resize(nt - 2);
for (Index t = 1; t < nt - 1; ++t)
{
em(t - 1) = acosh((corr(t - 1) + corr(t + 1))/(2.*corr(t)));
}
break;
case CorrelatorType::linear:
em.resize(nt - 1);
for (Index t = 1; t < nt; ++t)
{
em(t - 1) = corr(t) - corr(t - 1);
}
break;
case CorrelatorType::cst:
em = corr;
break;
}
return em;
}
DMatSample EffectiveMass::operator()(const DMatSample &corr) const
{
DMatSample em(corr.size());
FOR_STAT_ARRAY(em, s)
{
em[s] = (*this)(corr[s]);
}
return em;
}

50
lib/Physics/EffectiveMass.hpp Normal file
View File

@@ -0,0 +1,50 @@
/*
* EffectiveMass.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* LatAnalyze 3 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LatAnalyze 3. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef Latan_EffectiveMass_hpp_
#define Latan_EffectiveMass_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Statistics/MatSample.hpp>
#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
BEGIN_LATAN_NAMESPACE
/******************************************************************************
* Effective mass class *
******************************************************************************/
class EffectiveMass
{
public:
// constructors
EffectiveMass(const CorrelatorType type = CorrelatorType::exp);
// access
CorrelatorType getType(void) const;
void setType(const CorrelatorType type);
DVec getTime(const Index nt) const;
// compute effective mass
DVec operator()(const DVec &corr) const;
DMatSample operator()(const DMatSample &corr) const;
private:
CorrelatorType type_;
};
END_LATAN_NAMESPACE
#endif // Latan_EffectiveMass_hpp_

2
lib/Statistics/Dataset.hpp
View File

@@ -1,7 +1,7 @@
/*
* Dataset.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/FitInterface.cpp
View File

@@ -1,7 +1,7 @@
/*
* FitInterface.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/FitInterface.hpp
View File

@@ -1,7 +1,7 @@
/*
* FitInterface.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/Histogram.cpp
View File

@@ -1,7 +1,7 @@
/*
* Histogram.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/Histogram.hpp
View File

@@ -1,7 +1,7 @@
/*
* Histogram.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/MatSample.hpp
View File

@@ -1,7 +1,7 @@
/*
* MatSample.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/Random.cpp
View File

@@ -1,7 +1,7 @@
/*
* Random.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/Random.hpp
View File

@@ -1,7 +1,7 @@
/*
* Random.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2019 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/StatArray.cpp
View File

@@ -1,7 +1,7 @@
/*
* StatArray.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/StatArray.hpp
View File

@@ -1,7 +1,7 @@
/*
* StatArray.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/XYSampleData.cpp
View File

@@ -1,7 +1,7 @@
/*
* XYSampleData.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/XYSampleData.hpp
View File

@@ -1,7 +1,7 @@
/*
* XYSampleData.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/XYStatData.cpp
View File

@@ -1,7 +1,7 @@
/*
* XYStatData.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/Statistics/XYStatData.hpp
View File

@@ -1,7 +1,7 @@
/*
* XYStatData.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
lib/includes.hpp
View File

@@ -1,7 +1,7 @@
/*
* includes.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

355
physics/2pt-fit.cpp
View File

@@ -1,11 +1,13 @@
#include <LatAnalyze/Core/Math.hpp>
#include <LatAnalyze/Core/OptParser.hpp>
#include <LatAnalyze/Core/Plot.hpp>
#include <LatAnalyze/Functional/CompiledModel.hpp>
#include <LatAnalyze/Io/Io.hpp>
#include <LatAnalyze/Statistics/MatSample.hpp>
#include <LatAnalyze/Core/Math.hpp>
#include <LatAnalyze/Numerical/MinuitMinimizer.hpp>
#include <LatAnalyze/Numerical/NloptMinimizer.hpp>
#include <LatAnalyze/Core/Plot.hpp>
#include <LatAnalyze/Physics/CorrelatorFitter.hpp>
#include <LatAnalyze/Physics/EffectiveMass.hpp>
#include <LatAnalyze/Statistics/MatSample.hpp>
#include <LatAnalyze/Statistics/XYSampleData.hpp>
using namespace std;
@@ -17,17 +19,6 @@ struct TwoPtFit
Index tMin, tMax;
};
void setFitRange(XYSampleData &data, const Index ti, const Index tf,
const Index thinning, const Index nt)
{
for (Index t = 0; t < nt; ++t)
{
data.fitPoint((t >= ti) and (t <= tf)
and ((t - ti) % thinning == 0), t);
}
}
int main(int argc, char *argv[])
{
// parse arguments /////////////////////////////////////////////////////////
@@ -47,7 +38,7 @@ int main(int argc, char *argv[])
opt.addOption("s", "shift" , OptParser::OptType::value , true,
"time variable shift", "0");
opt.addOption("m", "model" , OptParser::OptType::value , true,
"fit model (exp|exp2|exp3|sinh|cosh|cosh2|cosh3|explin|const|<interpreter code>)", "cosh");
"fit model (exp<n>|sinh<n>|cosh<n>|linear|cst|<interpreter code>)", "exp1");
opt.addOption("" , "nPar" , OptParser::OptType::value , true,
"number of model parameters for custom models "
"(-1 if irrelevant)", "-1");
@@ -114,115 +105,25 @@ int main(int argc, char *argv[])
DMatSample tmp, corr;
Index nSample, nt;
tmp = Io::load<DMatSample>(corrFileName);
nSample = tmp.size();
nt = tmp[central].rows();
tmp = tmp.block(0, 0, nt, 1);
corr = tmp;
FOR_STAT_ARRAY(corr, s)
{
for (Index t = 0; t < nt; ++t)
{
corr[s]((t - shift + nt)%nt) = tmp[s](t);
}
}
corr = Io::load<DMatSample>(corrFileName);
nSample = corr.size();
nt = corr[central].rows();
corr = corr.block(0, 0, nt, 1);
corr = CorrelatorUtils::shift(corr, shift);
if (fold)
{
tmp = corr;
FOR_STAT_ARRAY(corr, s)
{
for (Index t = 0; t < nt; ++t)
{
corr[s](t) = 0.5*(tmp[s](t) + tmp[s]((nt - t) % nt));
}
}
corr = CorrelatorUtils::fold(corr);
}
// make models /////////////////////////////////////////////////////////////
DoubleModel mod;
bool sinhModel = false, coshModel = false, linearModel = false, constModel = false;
// make model //////////////////////////////////////////////////////////////
CorrelatorFitter fitter(corr);
DoubleModel mod;
auto modelPar = CorrelatorModels::parseModel(model);
if ((model == "exp") or (model == "exp1"))
if (modelPar.type != CorrelatorType::undefined)
{
nPar = 2;
mod.setFunction([](const double *x, const double *p)
{
return p[1]*exp(-p[0]*x[0]);
}, 1, nPar);
}
else if (model == "exp2")
{
nPar = 4;
mod.setFunction([](const double *x, const double *p)
{
return p[1]*exp(-p[0]*x[0]) + p[3]*exp(-p[2]*x[0]);
}, 1, nPar);
}
else if (model == "exp3")
{
nPar = 6;
mod.setFunction([](const double *x, const double *p)
{
return p[1]*exp(-p[0]*x[0]) + p[3]*exp(-p[2]*x[0])
+ p[5]*exp(-p[4]*x[0]);
}, 1, nPar);
}
else if (model == "sinh")
{
sinhModel = true;
nPar = 2;
mod.setFunction([nt](const double *x, const double *p)
{
return p[1]*(exp(-p[0]*x[0])-exp(-p[0]*(nt-x[0])));
}, 1, nPar);
}
else if ((model == "cosh") or (model == "cosh1"))
{
coshModel = true;
nPar = 2;
mod.setFunction([nt](const double *x, const double *p)
{
return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0])));
}, 1, nPar);
}
else if (model == "cosh2")
{
coshModel = true;
nPar = 4;
mod.setFunction([nt](const double *x, const double *p)
{
return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0])))
+ p[3]*(exp(-p[2]*x[0])+exp(-p[2]*(nt-x[0])));
}, 1, nPar);
}
else if (model == "cosh3")
{
coshModel = true;
nPar = 6;
mod.setFunction([nt](const double *x, const double *p)
{
return p[1]*(exp(-p[0]*x[0])+exp(-p[0]*(nt-x[0])))
+ p[3]*(exp(-p[2]*x[0])+exp(-p[2]*(nt-x[0])))
+ p[5]*(exp(-p[2]*x[0])+exp(-p[4]*(nt-x[0])));
}, 1, nPar);
}
else if (model == "explin")
{
linearModel = true;
nPar = 2;
mod.setFunction([](const double *x, const double *p)
{
return p[1] - p[0]*x[0];
}, 1, nPar);
}
else if (model == "const")
{
constModel = true;
nPar = 1;
mod.setFunction([](const double *x __dumb, const double *p)
{
return p[0];
}, 1, nPar);
mod = CorrelatorModels::makeModel(modelPar, nt);
nPar = mod.getNPar();
}
else
{
@@ -240,81 +141,52 @@ int main(int argc, char *argv[])
}
// fit /////////////////////////////////////////////////////////////////////
DMatSample tvec(nSample);
XYSampleData data(nSample);
DVec init(nPar);
NloptMinimizer globMin(NloptMinimizer::Algorithm::GN_CRS2_LM);
MinuitMinimizer locMin;
vector<Minimizer *> unCorrMin{&globMin, &locMin};
FOR_STAT_ARRAY(tvec, s)
{
tvec[s] = DVec::LinSpaced(nt, 0, nt - 1);
}
data.addXDim(nt, "t/a", true);
data.addYDim("C(t)");
data.setUnidimData(tvec, corr);
// set parameter name ******************************************************
if(constModel)
{
mod.parName().setName(0, "const");
}
else
{
for (Index p = 0; p < nPar; p += 2)
{
mod.parName().setName(p, "E_" + strFrom(p/2));
mod.parName().setName(p + 1, "Z_" + strFrom(p/2));
}
}
// set initial values ******************************************************
if (linearModel)
{
init(0) = data.y(nt/4, 0)[central] - data.y(nt/4 + 1, 0)[central];
init(1) = data.y(nt/4, 0)[central] + nt/4*init(0);
}
else if(constModel)
{
init(0) = data.y(nt/4, 0)[central];
// set fitter **************************************************************
fitter.setModel(mod);
fitter.data().setSvdTolerance(svdTol);
fitter.setThinning(thinning);
// set initial values ******************************************************
if (modelPar.type != CorrelatorType::undefined)
{
init = CorrelatorModels::parameterGuess(corr, modelPar);
}
else
{
init(0) = log(data.y(nt/4, 0)[central]/data.y(nt/4 + 1, 0)[central]);
init(1) = data.y(nt/4, 0)[central]/(exp(-init(0)*nt/4));
}
for (Index p = 2; p < nPar; p += 2)
{
init(p) = 2*init(p - 2);
init(p + 1) = init(p - 1)/2.;
init.fill(0.1);
}
// set limits for minimisers ***********************************************
for (Index p = 0; p < nPar; p += 2)
{
if (linearModel)
{
globMin.setLowLimit(p, -10.*fabs(init(p)));
globMin.setHighLimit(p, 10.*fabs(init(p)));
}
else if(constModel)
{
globMin.setLowLimit(p, -10*fabs(init(0)));
locMin.setLowLimit(p, -10*fabs(init(0)));
globMin.setHighLimit(p, 10*fabs(init(0)));
locMin.setHighLimit(p, 10*fabs(init(0)));
}
else
if ((modelPar.type == CorrelatorType::exp) or
(modelPar.type == CorrelatorType::cosh) or
(modelPar.type == CorrelatorType::sinh))
{
globMin.setLowLimit(p, 0.);
locMin.setLowLimit(p, 0.);
globMin.setHighLimit(p, 10.*init(p));
}
if(!constModel)
{
globMin.setLowLimit(p + 1, -10.*fabs(init(p + 1)));
globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1)));
}
else if(modelPar.type == CorrelatorType::linear)
{
globMin.setLowLimit(p, -10.*fabs(init(p)));
locMin.setLowLimit(p, -10.*fabs(init(p)));
globMin.setHighLimit(p, 10.*init(p));
globMin.setLowLimit(p + 1, -10.*fabs(init(p + 1)));
globMin.setHighLimit(p + 1, 10.*fabs(init(p + 1)));
}
else
{
globMin.setLowLimit(p, -10*fabs(init(p)));
globMin.setHighLimit(p, 10*fabs(init(p)));
}
}
globMin.setPrecision(0.001);
globMin.setMaxIteration(100000);
@@ -322,28 +194,28 @@ int main(int argc, char *argv[])
locMin.setMaxIteration(1000000);
locMin.setVerbosity(verbosity);
// fit /////////////////////////////////////////////////////////////////////
// standard fit ////////////////////////////////////////////////////////////
if (!doScan)
{
// fit *****************************************************************
SampleFitResult fit;
setFitRange(data, ti, tf, thinning, nt);
fitter.setFitRange(ti, tf);
if (doCorr)
{
cout << "-- uncorrelated fit..." << endl;
}
cout << "using model '" << model << "'" << endl;
data.setSvdTolerance(svdTol);
data.assumeYYCorrelated(false, 0, 0);
fit = data.fit(unCorrMin, init, mod);
fitter.setCorrelation(false);
fit = fitter.fit(unCorrMin, init);
fit.print();
if (doCorr)
{
cout << "-- correlated fit..." << endl;
cout << "using model '" << model << "'" << endl;
init = fit[central];
data.assumeYYCorrelated(true, 0, 0);
fit = data.fit(locMin, init, mod);
fitter.setCorrelation(true);
fit = fitter.fit(locMin, init);
fit.print();
}
if (!outFileName.empty())
@@ -353,84 +225,50 @@ int main(int argc, char *argv[])
// plots ***************************************************************
if (doPlot)
{
if (!constModel)
DMatSample tvec(nSample);
tvec.fill(DVec::LinSpaced(nt, 0, nt - 1));
if (modelPar.type != CorrelatorType::cst)
{
Plot p;
p << PlotRange(Axis::x, 0, nt - 1);
if (!linearModel and !constModel)
if ((modelPar.type == CorrelatorType::exp) or
(modelPar.type == CorrelatorType::cosh) or
(modelPar.type == CorrelatorType::sinh))
{
p << LogScale(Axis::y);
}
p << Color("rgb 'blue'") << PlotPredBand(fit.getModel(_), 0, nt - 1);
p << Color("rgb 'blue'") << PlotFunction(fit.getModel(), 0, nt - 1);
p << Color("rgb 'red'") << PlotData(data.getData());
p << Color("rgb 'red'") << PlotData(fitter.data().getData());
p << Label("t/a", Axis::x) << Caption("Correlator");
p.display();
if(savePlot != "")
{
p.save(savePlot + "_corr");
}
}
if (modelPar.type != CorrelatorType::undefined)
{
Plot p;
DMatSample effMass(nSample);
DVec effMassT, fitErr;
Index maxT = (coshModel) ? (nt - 2) : (nt - 1);
double e0, e0Err;
Plot p;
EffectiveMass effMass(modelPar.type);
DMatSample em;
DVec fitErr, emtvec;
double e0, e0Err;
effMass.resizeMat(maxT, 1);
effMassT.setLinSpaced(maxT, 0, maxT-1);
emtvec = effMass.getTime(nt);
em = effMass(corr);
fitErr = fit.variance().cwiseSqrt();
e0 = fit[central](0);
e0Err = fitErr(0);
if (coshModel or sinhModel)
{
FOR_STAT_ARRAY(effMass, s)
{
for (Index t = 1; t < nt - 1; ++t)
{
effMass[s](t - 1) = acosh((corr[s](t-1) + corr[s](t+1))
/(2.*corr[s](t)));
}
}
}
else if (linearModel)
{
FOR_STAT_ARRAY(effMass, s)
{
for (Index t = 0; t < nt - 1; ++t)
{
effMass[s](t) = corr[s](t) - corr[s](t+1);
}
}
}
else if (constModel)
{
FOR_STAT_ARRAY(effMass, s)
{
for (Index t = 0; t < nt - 1; ++t)
{
effMass[s](t) = corr[s](t);
}
}
}
else
{
FOR_STAT_ARRAY(effMass, s)
{
for (Index t = 1; t < nt; ++t)
{
effMass[s](t - 1) = log(corr[s](t-1)/corr[s](t));
}
}
}
p.reset();
p << PlotRange(Axis::x, 0, maxT);
p << PlotRange(Axis::y, e0 - 20.*e0Err, e0 + 20.*e0Err);
p << Color("rgb 'blue'") << PlotBand(0, maxT, e0 - e0Err, e0 + e0Err);
p << PlotRange(Axis::x, 0, nt - 1);
p << PlotRange(Axis::y, e0 - 30.*e0Err, e0 + 30.*e0Err);
p << Color("rgb 'blue'") << PlotBand(0, nt - 1, e0 - e0Err, e0 + e0Err);
p << Color("rgb 'blue'") << PlotHLine(e0);
p << Color("rgb 'red'") << PlotData(effMassT, effMass);
p << Caption("Effective Mass");
p << Color("rgb 'red'") << PlotData(emtvec, em);
p << Label("t/a", Axis::x) << Caption("Effective Mass");
p.display();
if(savePlot != "")
{
@@ -440,16 +278,19 @@ int main(int argc, char *argv[])
if (doHeatmap)
{
Plot p;
Index n = data.getFitVarMat().rows();
DMat id = DMat::Identity(n, n);
Index n = fitter.data().getFitVarMat().rows();
DMat id = DMat::Identity(n, n),
var = fitter.data().getFitVarMat();
p << PlotMatrix(Math::varToCorr(data.getFitVarMat()));
p << PlotMatrix(Math::varToCorr(var));
p << Caption("correlation matrix");
p.display();
if (svdTol > 0.)
{
DMat proj = id - var*fitter.data().getFitVarMatPInv();
p.reset();
p << PlotMatrix(id - data.getFitVarMat()*data.getFitVarMatPInv());
p << PlotMatrix(proj);
p << Caption("singular space projector");
p.display();
}
@@ -460,8 +301,9 @@ int main(int argc, char *argv[])
// scan fits ///////////////////////////////////////////////////////////////
else
{
// fits ****************************************************************
Index nFit = 0, f = 0, ti0 = ti + (tf - ti)/4, tf0 = tf - (tf - ti)/4,
matSize = tf - ti - nPar + 1;
matSize = tf - ti + 1;
DMat err, pVal(matSize, matSize), relErr(matSize, matSize),
ccdf(matSize, matSize), val(matSize, matSize);
map<double, TwoPtFit> fit;
@@ -474,14 +316,13 @@ int main(int argc, char *argv[])
<< endl;
thinning = 1;
}
setFitRange(data, ti0, tf0, thinning, nt);
data.setSvdTolerance(svdTol);
data.assumeYYCorrelated(false, 0, 0);
tmpFit = data.fit(unCorrMin, init, mod);
fitter.setFitRange(ti0, tf0);
fitter.setCorrelation(false);
tmpFit = fitter.fit(unCorrMin, init);
tmpFit.print();
cout << "-- scanning all possible fit ranges..." << endl;
init = tmpFit[central];
data.assumeYYCorrelated(doCorr, 0, 0);
fitter.setCorrelation(doCorr);
pVal.fill(Math::nan);
relErr.fill(Math::nan);
val.fill(Math::nan);
@@ -496,8 +337,8 @@ int main(int argc, char *argv[])
{
Index i = ta - ti, j = tb - ti;
setFitRange(data, ta, tb, thinning, nt);
tmpFit = data.fit(locMin, init, mod);
fitter.setFitRange(ta, tb);
tmpFit = fitter.fit(locMin, init);
err = tmpFit.variance().cwiseSqrt();
pVal(i, j) = tmpFit.getPValue();
ccdf(i, j) = tmpFit.getCcdf();
@@ -531,8 +372,8 @@ int main(int argc, char *argv[])
p << PlotMatrix(pVal);
p << Caption("p-value matrix");
p << Label("tMin - " + strFrom(ti), Axis::x);
p << Label("tMax - " + strFrom(ti), Axis::y);
p << Label("tMax - " + strFrom(ti), Axis::x);
p << Label("tMin - " + strFrom(ti), Axis::y);
p.display();
if(savePlot != "")
{
@@ -541,8 +382,8 @@ int main(int argc, char *argv[])
p.reset();
p << PlotMatrix(relErr);
p << Caption("Relative error matrix");
p << Label("tMin - " + strFrom(ti), Axis::x);
p << Label("tMax - " + strFrom(ti), Axis::y);
p << Label("tMax - " + strFrom(ti), Axis::x);
p << Label("tMin - " + strFrom(ti), Axis::y);
p.display();
if(savePlot != "")
{
@@ -551,8 +392,8 @@ int main(int argc, char *argv[])
p.reset();
p << PlotMatrix(val);
p << Caption("Fit result matrix");
p << Label("tMin - " + strFrom(ti), Axis::x);
p << Label("tMax - " + strFrom(ti), Axis::y);
p << Label("tMax - " + strFrom(ti), Axis::x);
p << Label("tMin - " + strFrom(ti), Axis::y);
p.display();
if(savePlot != "")
{
@@ -561,8 +402,8 @@ int main(int argc, char *argv[])
p.reset();
p << PlotMatrix(ccdf);
p << Caption("chi^2 CCDF matrix");
p << Label("tMin - " + strFrom(ti), Axis::x);
p << Label("tMax - " + strFrom(ti), Axis::y);
p << Label("tMax - " + strFrom(ti), Axis::x);
p << Label("tMin - " + strFrom(ti), Axis::y);
p.display();
if(savePlot != "")
{

2
update_eigen.sh
View File

@@ -6,7 +6,7 @@ if (( $# != 1 )); then
fi
ARC=$1
INITDIR=`pwd`
INITDIR=$(pwd -P)
rm -rf lib/Eigen
ARCDIR=`tar -tf ${ARC} | head -n1 | sed -e 's@/.*@@'`
tar -xf ${ARC}

2
utils/resample.cpp
View File

@@ -1,7 +1,7 @@
/*
* resample.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
utils/sample-combine.cpp
View File

@@ -1,7 +1,7 @@
/*
* sample_combine.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
utils/sample-element.cpp
View File

@@ -1,7 +1,7 @@
/*
* sample-element.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli, Matt Spraggs
* Copyright (C) 2013 - 2020 Antonin Portelli, Matt Spraggs
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

2
utils/sample-fake.cpp
View File

@@ -1,7 +1,7 @@
/*
* make_fake_sample.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2016 Antonin Portelli
* Copyright (C) 2013 - 2020 Antonin Portelli
*
* LatAnalyze 3 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by

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