portelli/LatAnalyze
1
0
Fork 0
mirror of https://github.com/aportelli/LatAnalyze.git synced 2025-06-20 00:06:55 +01:00
Code Releases Activity

reorganizing sources...

Browse Source
This commit is contained in:
Antonin Portelli
2014-03-13 18:51:01 +00:00
parent 0109379dc3
commit 5b65dc8fc9
68 changed files with 245 additions and 245 deletions
Download Patch File Download Diff File Expand all files Collapse all files

187
lib/AsciiFile.cpp Normal file
View File

@ -0,0 +1,187 @@
/*
* AsciiFile.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/AsciiFile.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* AsciiFile implementation *
******************************************************************************/
// AsciiParserState constructor ////////////////////////////////////////////////
AsciiFile::AsciiParserState::AsciiParserState(istream *stream, string *name,
IoDataTable *data)
: ParserState<IoDataTable>(stream, name, data)
{
initScanner();
}
// AsciiParserState destructor /////////////////////////////////////////////////
AsciiFile::AsciiParserState::~AsciiParserState(void)
{
destroyScanner();
}
// constructor /////////////////////////////////////////////////////////////////
AsciiFile::AsciiFile(const string &name, const unsigned int mode)
{
open(name, mode);
}
// destructor //////////////////////////////////////////////////////////////////
AsciiFile::~AsciiFile(void)
{
close();
}
// access //////////////////////////////////////////////////////////////////////
void AsciiFile::save(const DMat &m, const std::string &name)
{
checkWritability();
isParsed_ = false;
fileStream_ << "#L latan_begin mat " << name << endl;
fileStream_ << m.cols() << endl;
fileStream_ << scientific << m << endl;
fileStream_ << "#L latan_end mat " << endl;
}
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)
{
save(s[i], name + "_S_" + strFrom(i));
}
fileStream_ << "#L latan_end rs_sample " << endl;
}
void AsciiFile::save(const RandGenState &state, const std::string &name)
{
checkWritability();
isParsed_ = false;
fileStream_ << "#L latan_begin rg_state " << name << endl;
fileStream_ << state << endl;
fileStream_ << "#L latan_end rg_state " << endl;
}
// tests ///////////////////////////////////////////////////////////////////////
bool AsciiFile::isOpen() const
{
return fileStream_.is_open();
}
// IO //////////////////////////////////////////////////////////////////////////
void AsciiFile::close(void)
{
state_.reset();
if (isOpen())
{
fileStream_.close();
}
name_ = "";
mode_ = Mode::null;
isParsed_ = false;
deleteData();
}
void AsciiFile::open(const string &name, const unsigned int mode)
{
if (isOpen())
{
LATAN_ERROR(Io, "file already opened with name '" + name_ + "'");
}
else
{
ios_base::openmode stdMode = static_cast<ios_base::openmode>(0);
if (mode & Mode::write)
{
stdMode |= ios::out|ios::trunc;
}
if (mode & Mode::read)
{
stdMode |= ios::in;
}
if (mode & Mode::append)
{
stdMode |= ios::out|ios::app;
}
name_ = name;
mode_ = mode;
isParsed_ = false;
fileStream_.open(name_.c_str(), stdMode);
if (mode_ & Mode::read)
{
state_.reset(new AsciiParserState(&fileStream_, &name_, &data_));
}
else
{
state_.reset();
}
}
}
std::string AsciiFile::load(const string &name)
{
if ((mode_ & Mode::read)&&(isOpen()))
{
if (!isParsed_)
{
state_->isFirst = true;
parse();
}
if (name.empty())
{
return state_->first;
}
else
{
return name;
}
}
else
{
if (isOpen())
{
LATAN_ERROR(Io, "file '" + name_ + "' is not opened in read mode");
}
else
{
LATAN_ERROR(Io, "file not opened");
}
}
}
// parser //////////////////////////////////////////////////////////////////////
//// Bison/Flex parser declaration
int _Ascii_parse(AsciiFile::AsciiParserState *state);
void AsciiFile::parse()
{
fileStream_.seekg(0);
_Ascii_parse(state_.get());
isParsed_ = true;
}

88
lib/AsciiFile.hpp Normal file
View File

@ -0,0 +1,88 @@
/*
* AsciiFile.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_AsciiFile_hpp_
#define Latan_AsciiFile_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/File.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <fstream>
BEGIN_NAMESPACE
/******************************************************************************
* ASCII datafile class *
******************************************************************************/
class AsciiFile: public File
{
public:
class AsciiParserState: public ParserState<IoDataTable>
{
public:
// constructor
AsciiParserState(std::istream *stream, std::string *name,
IoDataTable *data);
// destructor
virtual ~AsciiParserState(void);
// first element reference
bool isFirst;
std::string first;
// parsing buffers
RandGenState stateBuf;
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);
virtual void destroyScanner(void);
};
public:
// constructors
AsciiFile(void) = default;
AsciiFile(const std::string &name, const unsigned int mode);
// destructor
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);
// tests
virtual bool isOpen(void) const;
// IO
virtual void close(void);
virtual void open(const std::string &name, const unsigned int mode);
private:
// IO
virtual std::string load(const std::string &name = "");
// parser
void parse(void);
private:
std::fstream fileStream_;
bool isParsed_{false};
std::unique_ptr<AsciiParserState> state_{nullptr};
};
END_NAMESPACE
#endif // Latan_AsciiFile_hpp_

95
lib/AsciiLexer.lpp Normal file
View File

@ -0,0 +1,95 @@
/*
* AsciiLexer.lpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/>.
*/
%option reentrant
%option prefix="_Ascii_"
%option bison-bridge
%option bison-locations
%option noyywrap
%option yylineno
%{
#include <iostream>
#include <LatAnalyze/AsciiFile.hpp>
#include <LatAnalyze/AsciiParser.hpp>
#pragma GCC diagnostic ignored "-Wsign-compare"
#pragma GCC diagnostic ignored "-Wunused-function"
#pragma GCC diagnostic ignored "-Wunused-parameter"
using namespace std;
using namespace Latan;
#define YY_EXTRA_TYPE AsciiFile::AsciiParserState*
#define YY_USER_ACTION \
yylloc->first_line = yylloc->last_line = yylineno;\
yylloc->first_column = yylloc->last_column + 1;\
yylloc->last_column = yylloc->first_column + yyleng - 1;
#define YY_INPUT(buf, result, max_size) \
{ \
(*yyextra->stream).read(buf,max_size);\
result = (*yyextra->stream).gcount();\
}
#define YY_DEBUG 0
#if (YY_DEBUG == 1)
#define RETTOK(tok) cout << #tok << "(" << yytext << ")" << endl; return tok
#else
#define RETTOK(tok) return tok
#endif
%}
DIGIT [0-9]
ALPHA [a-zA-Z_+.-]
SIGN \+|-
EXP e|E
INT {SIGN}?{DIGIT}+
FLOAT {SIGN}?(({DIGIT}+\.{DIGIT}*)|({DIGIT}*\.{DIGIT}+))({EXP}{SIGN}?{INT}+)?
LMARK #L
BLANK [ \t]
%x MARK TYPE
%%
{LMARK} {BEGIN(MARK);}
{INT} {yylval->val_int = strTo<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);}
<MARK>latan_end {BEGIN(TYPE); RETTOK(CLOSE);}
<TYPE>mat {BEGIN(INITIAL); RETTOK(MAT);}
<TYPE>rs_sample {BEGIN(INITIAL); RETTOK(SAMPLE);}
<TYPE>rg_state {BEGIN(INITIAL); RETTOK(RG_STATE);}
<*>\n {yylloc->last_column = 0;}
<*>[ \t]
<*>. {yylval->val_char = yytext[0]; RETTOK(ERR);}
%%
void AsciiFile::AsciiParserState::initScanner()
{
yylex_init(&scanner);
yyset_extra(this, scanner);
}
void AsciiFile::AsciiParserState::destroyScanner()
{
yylex_destroy(scanner);
}

187
lib/AsciiParser.ypp Normal file
View File

@ -0,0 +1,187 @@
/*
* AsciiParser.ypp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Global.hpp>
#include <LatAnalyze/AsciiFile.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <iostream>
#include <sstream>
#include <utility>
#include <cstring>
using namespace std;
using namespace Latan;
#define TEST_FIRST(name) \
if (state->isFirst)\
{\
state->first = (name);\
state->isFirst = false;\
}
%}
%pure-parser
%name-prefix="_Ascii_"
%locations
%defines
%error-verbose
%parse-param { Latan::AsciiFile::AsciiParserState* state }
%initial-action {yylloc.last_column = 0;}
%lex-param { void* scanner }
%union
{
int val_int;
double val_double;
char val_char;
char val_str[256];
}
%token <val_char> ERR
%token <val_double> FLOAT
%token <val_int> INT
%token <val_str> ID
%token OPEN CLOSE MAT SAMPLE RG_STATE
%type <val_str> mat sample rg_state
%{
int _Ascii_lex(YYSTYPE* lvalp, YYLTYPE* llocp, void* scanner);
void _Ascii_error(YYLTYPE* locp, AsciiFile::AsciiParserState* state,
const char* err)
{
stringstream buf;
buf << *state->streamName << ":" << locp->first_line << ":"\
<< locp->first_column << ": " << err;
LATAN_ERROR(Parsing, buf.str());
}
#define scanner state->scanner
%}
%%
datas:
/* empty string */
| datas data
;
data:
mat
{
TEST_FIRST($1);
(*state->data)[$1].reset(new DMat(state->dMatQueue.front()));
state->dMatQueue.pop();
}
| sample
{
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:
OPEN MAT ID INT floats CLOSE MAT
{
const unsigned int nRow = state->doubleQueue.size()/$INT, nCol = $INT;
Index i, j, r = 0;
if (state->doubleQueue.size() != nRow*nCol)
{
LATAN_ERROR(Size, "matrix '" + *state->streamName + ":" + $ID +
"' has a wrong size");
}
state->dMatQueue.push(DMat(nRow, nCol));
while (!state->doubleQueue.empty())
{
j = r % nCol;
i = (r - j)/nCol;
state->dMatQueue.back()(i, j) = state->doubleQueue.front();
state->doubleQueue.pop();
r++;
}
strcpy($$, $ID);
}
;
sample:
OPEN SAMPLE ID INT mats CLOSE SAMPLE
{
const unsigned int nSample = $INT;
if (state->dMatQueue.size() != nSample + DMatSample::offset)
{
LATAN_ERROR(Size, "sample '" + *state->streamName + ":" + $ID +
"' has a wrong size");
}
state->dMatSampleBuf.resize(nSample);
state->dMatSampleBuf[central] = state->dMatQueue.front();
state->dMatQueue.pop();
for (int i = 0; i < $INT; ++i)
{
state->dMatSampleBuf[i] = state->dMatQueue.front();
state->dMatQueue.pop();
}
strcpy($$, $ID);
}
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
+ ":" + $ID + "' has a wrong size");
}
for (Index i = 0; i < RLXG_STATE_SIZE; ++i)
{
state->stateBuf[i] = state->intQueue.front();
state->intQueue.pop();
}
strcpy($$, $ID);
}
;
mats:
mats mat
| mat
;
floats:
floats FLOAT {state->doubleQueue.push($FLOAT);}
| floats INT {state->doubleQueue.push(static_cast<double>($INT));}
| FLOAT {state->doubleQueue.push($FLOAT);}
| INT {state->doubleQueue.push(static_cast<double>($INT));}
;
ints:
ints INT {state->intQueue.push($INT);}
| INT {state->intQueue.push($INT);}
;

301
lib/Chi2Function.cpp Normal file
View File

@ -0,0 +1,301 @@
/*
* Chi2Function.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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)
: DoubleFunction(0)
, 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)
{
if (static_cast<Index>(model_.size()) != data_.getYDim())
{
model_.resize(static_cast<unsigned int>(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<unsigned int>(j)] = &model;
nPar_ = model.getNPar();
}
void Chi2Function::setModel(const vector<const DoubleModel *> &modelVector)
{
if (static_cast<Index>(model_.size()) != data_.getYDim())
{
model_.resize(static_cast<unsigned int>(data_.getYDim()));
}
if (modelVector.size() != static_cast<unsigned int>(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<unsigned int>(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<DMatBase> 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.inverse().eval();
buffer_->isInit = true;
}
// function call ///////////////////////////////////////////////////////////////
double Chi2Function::operator()(const double *arg) const
{
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> p(arg, nPar_, 1);
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
is = 0;
for (Index j = 0; j < yDim; ++j)
FOR_VEC(buffer_->dInd, k)
{
const DoubleModel *f = model_[static_cast<unsigned int>(j)];
double f_jk, y_jk = data_.y(j, buffer_->dInd(k));
if (!f)
{
LATAN_ERROR(Memory, "null model");
}
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++;
}
}
f_jk = (*f)(buffer_->x, p);
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(i*nPoint + k) = xi_ik - x_ik;
}
// compute result
res = buffer_->v.dot(buffer_->invVar*buffer_->v);
return res;
}

81
lib/Chi2Function.hpp Normal file
View File

@ -0,0 +1,81 @@
/*
* Chi2Function.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_NAMESPACE
/******************************************************************************
* <Class> *
******************************************************************************/
// forward declaration of XYStatData
class XYStatData;
class Chi2Function: public DoubleFunction
{
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
using DoubleFunction::operator();
protected:
// function call
virtual double operator()(const double *arg) const;
private:
// access
void resizeBuffer(void) const;
// compute variance matrix inverse
void setVarianceBlock(const Index l1, const Index l2,
ConstBlock<DMatBase> m) const;
void initBuffer(void) const;
private:
const XYStatData &data_;
std::shared_ptr<Chi2FunctionBuffer> buffer_;
std::vector<const DoubleModel *> model_;
Index nPar_{-1};
};
END_NAMESPACE
#endif // Latan_Chi2Function_hpp_

81
lib/CompiledFunction.cpp Normal file
View File

@ -0,0 +1,81 @@
/*
* CompiledFunction.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/CompiledFunction.hpp>
#include <LatAnalyze/Math.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* Compiled double function implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
CompiledDoubleFunction::CompiledDoubleFunction(const unsigned nArg)
: DoubleFunction(nArg, nullptr)
{}
CompiledDoubleFunction::CompiledDoubleFunction(const unsigned nArg,
const string &code)
: CompiledDoubleFunction(nArg)
{
setCode(code);
}
// access //////////////////////////////////////////////////////////////////////
void CompiledDoubleFunction::setCode(const string &code)
{
interpreter_.reset(new MathInterpreter(code));
context_.reset(new RunContext);
StdMath::addStdMathFunc(context_->fTable);
}
// function call ///////////////////////////////////////////////////////////////
double CompiledDoubleFunction::operator()(const double *arg) const
{
double result;
for (Index i = 0; i < getNArg(); ++i)
{
context_->vTable["x_" + strFrom(i)] = arg[i];
}
(*interpreter_)(*context_);
if (!context_->dStack.empty())
{
result = context_->dStack.top();
context_->dStack.pop();
}
else
{
result = 0.0;
LATAN_ERROR(Program, "program execution resulted in an empty stack");
}
return result;
}
// IO //////////////////////////////////////////////////////////////////////////
ostream & Latan::operator<<(ostream &out, CompiledDoubleFunction &f)
{
f.interpreter_->compile();
out << *(f.interpreter_);
return out;
}

57
lib/CompiledFunction.hpp Normal file
View File

@ -0,0 +1,57 @@
/*
* CompiledFunction.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_CompiledFunction_hpp_
#define Latan_CompiledFunction_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/MathInterpreter.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* compiled double function class *
******************************************************************************/
class CompiledDoubleFunction: public DoubleFunction
{
public:
// constructors
explicit CompiledDoubleFunction(const unsigned nArg);
CompiledDoubleFunction(const unsigned nArg, const std::string &code);
// destructor
virtual ~CompiledDoubleFunction(void) = default;
// access
void setCode(const std::string &code);
// function call
using DoubleFunction::operator();
virtual double operator()(const double *arg) const;
// IO
friend std::ostream & operator<<(std::ostream &out,
CompiledDoubleFunction &f);
private:
std::shared_ptr<MathInterpreter> interpreter_;
std::shared_ptr<RunContext> context_;
};
std::ostream & operator<<(std::ostream &out, CompiledDoubleFunction &f);
END_NAMESPACE
#endif // Latan_CompiledFunction_hpp_

88
lib/CompiledModel.cpp Normal file
View File

@ -0,0 +1,88 @@
/*
* CompiledModel.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/CompiledModel.hpp>
#include <LatAnalyze/Math.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* CompiledDoubleModel implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
CompiledDoubleModel::CompiledDoubleModel(const unsigned nArg,
const unsigned nPar)
: DoubleModel(nArg, nPar, nullptr)
{}
CompiledDoubleModel::CompiledDoubleModel(const unsigned nArg,
const unsigned nPar,
const string &code)
: CompiledDoubleModel(nArg, nPar)
{
setCode(code);
}
// access //////////////////////////////////////////////////////////////////////
void CompiledDoubleModel::setCode(const std::string &code)
{
interpreter_.reset(new MathInterpreter(code));
context_.reset(new RunContext);
StdMath::addStdMathFunc(context_->fTable);
}
// function call ///////////////////////////////////////////////////////////////
double CompiledDoubleModel::operator()(const double *arg,
const double *par) const
{
double result;
for (Index i = 0; i < getNArg(); ++i)
{
context_->vTable["x_" + strFrom(i)] = arg[i];
}
for (Index j = 0; j < getNPar(); ++j)
{
context_->vTable["p_" + strFrom(j)] = par[j];
}
(*interpreter_)(*context_);
if (!context_->dStack.empty())
{
result = context_->dStack.top();
context_->dStack.pop();
}
else
{
result = 0.0;
LATAN_ERROR(Program, "program execution resulted in an empty stack");
}
return result;
}
// IO //////////////////////////////////////////////////////////////////////////
ostream & Latan::operator<<(std::ostream &out, CompiledDoubleModel &m)
{
m.interpreter_->compile();
out << *(m.interpreter_);
return out;
}

59
lib/CompiledModel.hpp Normal file
View File

@ -0,0 +1,59 @@
/*
* CompiledModel.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_CompiledModel_hpp_
#define Latan_CompiledModel_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Model.hpp>
#include <LatAnalyze/MathInterpreter.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* compiled double model class *
******************************************************************************/
class CompiledDoubleModel: public DoubleModel
{
public:
// constructor
explicit CompiledDoubleModel(const unsigned nArg, const unsigned nPar);
CompiledDoubleModel(const unsigned nArg, const unsigned nPar,
const std::string &code);
// destructor
virtual ~CompiledDoubleModel(void) = default;
// access
void setCode(const std::string &code);
// function call
using DoubleModel::operator();
virtual double operator()(const double *arg, const double *par) const;
// IO
friend std::ostream & operator<<(std::ostream &out,
CompiledDoubleModel &f);
private:
std::shared_ptr<MathInterpreter> interpreter_;
std::shared_ptr<RunContext> context_;
};
std::ostream & operator<<(std::ostream &out, CompiledDoubleModel &f);
END_NAMESPACE
#endif // Latan_CompiledModel_hpp_

141
lib/Dataset.hpp Normal file
View File

@ -0,0 +1,141 @@
/*
* Dataset.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Dataset_hpp_
#define Latan_Dataset_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/File.hpp>
#include <LatAnalyze/StatArray.hpp>
#include <LatAnalyze/RandGen.hpp>
#include <fstream>
#include <vector>
BEGIN_NAMESPACE
/******************************************************************************
* Dataset class *
******************************************************************************/
template <typename T>
class Dataset: public StatArray<T>
{
private:
typedef StatArray<T> Base;
public:
// constructors
Dataset(void) = default;
template <typename FileType>
Dataset(const std::string &listFileName, const std::string &dataName);
using Base::Base;
// destructor
virtual ~Dataset(void) = default;
// IO
template <typename FileType>
void load(const std::string &listFileName, const std::string &dataName);
// resampling
Sample<T> bootstrapMean(const Index nSample, RandGen& generator);
private:
// mean from pointer vector for resampling
void ptVectorMean(T &m, const std::vector<const T *> &v);
};
/******************************************************************************
* Dataset template implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename T>
template <typename FileType>
Dataset<T>::Dataset(const std::string &listFileName,
const std::string &dataName)
{
load<FileType>(listFileName, dataName);
}
// IO //////////////////////////////////////////////////////////////////////////
template <typename T>
template <typename FileType>
void Dataset<T>::load(const std::string &listFileName,
const std::string &dataName)
{
FileType file;
std::ifstream listFile;
char dataFileNameBuf[MAX_PATH_LENGTH];
std::vector<std::string> dataFileName;
listFile.open(listFileName, std::ios::in);
while (!listFile.eof())
{
listFile.getline(dataFileNameBuf, MAX_PATH_LENGTH);
if (!std::string(dataFileNameBuf).empty())
{
dataFileName.push_back(dataFileNameBuf);
}
}
listFile.close();
this->resize(dataFileName.size());
for (Index i = 0; i < static_cast<Index>(dataFileName.size()); ++i)
{
file.open(dataFileName[i], File::Mode::read);
(*this)[i] = file.template read<T>(dataName);
file.close();
}
}
// resampling //////////////////////////////////////////////////////////////////
template <typename T>
Sample<T> Dataset<T>::bootstrapMean(const Index nSample, RandGen& generator)
{
Index nData = this->size();
std::vector<const T *> data(nData);
Sample<T> s(nSample);
for (Index j = 0; j < nData; ++j)
{
data[j] = &((*this)[j]);
}
ptVectorMean(s[central], data);
for (Index i = 0; i < nSample; ++i)
{
for (Index j = 0; j < nData; ++j)
{
data[j] = &((*this)[generator.discreteUniform(nData)]);
}
ptVectorMean(s[i], data);
}
return s;
}
template <typename T>
void Dataset<T>::ptVectorMean(T &m, const std::vector<const T *> &v)
{
if (v.size())
{
m = *(v[0]);
for (unsigned int i = 1; i < v.size(); ++i)
{
m += *(v[i]);
}
m /= static_cast<double>(v.size());
}
}
END_NAMESPACE
#endif // Latan_Dataset_hpp_

50
lib/Exceptions.cpp Normal file
View File

@ -0,0 +1,50 @@
/*
* Exceptions.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Exceptions.hpp>
#include <LatAnalyze/includes.hpp>
#ifndef ERR_SUFF
#define ERR_SUFF " (" + loc + ")"
#endif
#define CONST_EXC(name, init) \
name::name(string msg, string loc)\
:init\
{}
using namespace std;
using namespace Latan;
using namespace Exceptions;
// logic errors
CONST_EXC(Logic, logic_error(Env::msgPrefix + msg + ERR_SUFF))
CONST_EXC(Definition, Logic("definition error: " + msg, loc))
CONST_EXC(Implementation, Logic("implementation error: " + msg, loc))
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(Compilation, Runtime("compilation error: " + msg, loc))
CONST_EXC(Io, Runtime("IO error: " + msg, loc))
CONST_EXC(Memory, Runtime("memory error: " + msg, loc))
CONST_EXC(Parsing, Runtime(msg, loc))
CONST_EXC(Program, Runtime(msg, loc))
CONST_EXC(Syntax, Runtime("syntax error: " + msg, loc))
CONST_EXC(System, Runtime("system error: " + msg, loc))

65
lib/Exceptions.hpp Normal file
View File

@ -0,0 +1,65 @@
/*
* Exceptions.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Exceptions_hpp_
#define Latan_Exceptions_hpp_
#include <stdexcept>
#ifndef LATAN_GLOBAL_HPP_
#include <LatAnalyze/Global.hpp>
#endif
#define SRC_LOC strFrom(__FUNCTION__) + " at " + strFrom(__FILE__) + ":"\
+ strFrom(__LINE__)
#define LATAN_ERROR(exc,msg) throw(Exceptions::exc(msg, SRC_LOC))
#define LATAN_WARNING(msg) \
std::cerr << Env::msgPrefix << "warning: " << msg\
<< " (" << SRC_LOC << ")" << std::endl
#define DECL_EXC(name, base) \
class name: public base\
{\
public:\
name(std::string msg, std::string loc);\
}
BEGIN_NAMESPACE
namespace Exceptions
{
// logic errors
DECL_EXC(Logic, std::logic_error);
DECL_EXC(Definition, Logic);
DECL_EXC(Implementation, Logic);
DECL_EXC(Range, Logic);
DECL_EXC(Size, Logic);
// runtime errors
DECL_EXC(Runtime, std::runtime_error);
DECL_EXC(Compilation, Runtime);
DECL_EXC(Io, Runtime);
DECL_EXC(Memory, Runtime);
DECL_EXC(Parsing, Runtime);
DECL_EXC(Program, Runtime);
DECL_EXC(Syntax, Runtime);
DECL_EXC(System, Runtime);
}
END_NAMESPACE
#endif // Latan_Exceptions_hpp_

68
lib/File.cpp Normal file
View File

@ -0,0 +1,68 @@
/*
* File.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/File.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* File implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
File::File(const string &name, const unsigned int mode)
: name_(name)
, mode_(mode)
{}
// destructor //////////////////////////////////////////////////////////////////
File::~File(void)
{
deleteData();
}
// access //////////////////////////////////////////////////////////////////////
const string & File::getName(void) const
{
return name_;
}
unsigned int File::getMode(void) const
{
return mode_;
}
// internal functions //////////////////////////////////////////////////////////
void File::deleteData(void)
{
for (auto &i : data_)
{
i.second.reset();
}
data_.clear();
}
void File::checkWritability(void)
{
if (!((mode_ & Mode::write)||(mode_ & Mode::append))||!isOpen())
{
LATAN_ERROR(Io, "file '" + name_ + "' is not writable");
}
}

157
lib/File.hpp Normal file
View File

@ -0,0 +1,157 @@
/*
* File.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Io_hpp_
#define Latan_Io_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/IoObject.hpp>
#include <LatAnalyze/ParserState.hpp>
#include <queue>
#include <unordered_map>
BEGIN_NAMESPACE
// forward declaration of IO types
class DMat;
class DMatSample;
class RandGen;
class RandGenState;
/******************************************************************************
* Abstract datafile class *
******************************************************************************/
typedef std::unordered_map<std::string, std::unique_ptr<IoObject>> IoDataTable;
class File
{
public:
class Mode
{
public:
enum
{
null = 0,
write = 1 << 0,
read = 1 << 1,
append = 1 << 2
};
};
public:
// constructors
File(void) = default;
File(const std::string &name, const unsigned int mode);
// destructor
virtual ~File(void);
// access
const std::string & getName(void) const;
unsigned int getMode(void) const;
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;
// tests
virtual bool isOpen(void) const = 0;
// IO
virtual void close(void) = 0;
virtual void open(const std::string &name, const unsigned int mode) = 0;
// static IO functions
protected:
// access
void setName(const std::string &name);
void setMode(const unsigned int mode);
// data access
void deleteData(void);
// error checking
void checkWritability(void);
private:
// data access
template <typename IoT>
const IoT& getData(const std::string &name = "") const;
// IO
virtual std::string load(const std::string &name = "") = 0;
protected:
std::string name_{""};
unsigned int mode_{Mode::null};
IoDataTable data_;
};
// Template implementations
template <typename IoT>
const IoT& File::read(const std::string &name)
{
std::string dataName;
dataName = load(name);
return getData<IoT>(dataName);
}
template <typename IoT>
const IoT& File::getData(const std::string &name) const
{
try
{
return dynamic_cast<const IoT &>(*(data_.at(name)));
}
catch(std::out_of_range)
{
LATAN_ERROR(Definition, "no data with name '" + name + "' in file "
+ name_);
}
}
/******************************************************************************
* Static IO functions *
******************************************************************************/
class Io
{
public:
template <typename IoT, typename FileType>
static IoT load(const std::string &fileName, const std::string &name = "");
template <typename IoT, typename FileType>
static void save(const IoT &data, const std::string &fileName,
const unsigned int mode = File::Mode::append,
const std::string &name = "");
};
// template implementation /////////////////////////////////////////////////////
template <typename IoT, typename FileType>
IoT Io::load(const std::string &fileName, const std::string &name)
{
FileType file(fileName, File::Mode::read);
return file.template read<IoT>(name);
}
template <typename IoT, typename FileType>
void Io::save(const IoT &data, const std::string &fileName,
const unsigned int mode, const std::string &name)
{
FileType file(fileName, mode);
std::string realName = (name.empty()) ? fileName : name;
file.template save(data, realName);
}
END_NAMESPACE
#endif

178
lib/FitInterface.cpp Normal file
View File

@ -0,0 +1,178 @@
/*
* FitInterface.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/FitInterface.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* FitInterface implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
FitInterface::FitInterface(const Index nData, const Index xDim,
const Index yDim)
{
resize(nData, xDim, yDim);
}
// access //////////////////////////////////////////////////////////////////////
void FitInterface::assumeXExact(const Index i, const bool isExact)
{
isXExact_(i) = (isExact) ? 1 : 0;
}
void FitInterface::assumeXXCorrelated(const Index i1, const Index i2,
const bool isCorrelated)
{
isXXCorr_(i1, i2) = (isCorrelated) ? 1 : 0;
}
void FitInterface::assumeYYCorrelated(const Index j1, const Index j2,
const bool isCorrelated)
{
isYYCorr_(j1, j2) = (isCorrelated) ? 1 : 0;
}
void FitInterface::assumeYXCorrelated(const Index j, const Index i,
const bool isCorrelated)
{
isYXCorr_(j, i) = (isCorrelated) ? 1 : 0;
}
void FitInterface::assumeDataCorrelated(const Index k1, const Index k2,
const bool isCorrelated)
{
isDataCorr_(k1, k2) = (isCorrelated) ? 1 : 0;
}
void FitInterface::fitPoint(const Index i, const bool isFitPoint)
{
isFitPoint_(i) = (isFitPoint) ? 1 : 0;
}
void FitInterface::fitPointRange(const Index k1, const Index k2,
const bool isFitPoint)
{
int size = static_cast<int>(k2-k1+1);
isFitPoint_.segment(k1, size) = IVec::Constant(size, (isFitPoint) ? 1 : 0);
}
void FitInterface::fitAllPoints(const bool isFitPoint)
{
fitPointRange(0, getNData()-1, isFitPoint);
}
Index FitInterface::getNData(void) const
{
return isFitPoint_.size();
}
Index FitInterface::getNFitPoint(void) const
{
return isFitPoint_.sum();
}
Index FitInterface::getXDim(void) const
{
return isXXCorr_.rows();
}
Index FitInterface::getYDim(void) const
{
return isYYCorr_.rows();
}
Index FitInterface::getStatXDim(void) const
{
return isXExact_.size() - isXExact_.sum();
}
void FitInterface::setFitInterface(const FitInterface &fitInterface)
{
if (&fitInterface != this)
{
resize(fitInterface.getNData(), fitInterface.getXDim(),
fitInterface.getYDim());
isXExact_ = fitInterface.isXExact_;
isFitPoint_ = fitInterface.isFitPoint_;
isXXCorr_ = fitInterface.isXXCorr_;
isYYCorr_ = fitInterface.isYYCorr_;
isYXCorr_ = fitInterface.isYXCorr_;
isDataCorr_ = fitInterface.isDataCorr_;
}
}
void FitInterface::setNData(const Index nData)
{
resize(nData, getXDim(), getYDim());
}
void FitInterface::setXDim(const Index xDim)
{
resize(getNData(), xDim, getYDim());
}
void FitInterface::setYDim(const Index yDim)
{
resize(getNData(), getXDim(), yDim);
}
void FitInterface::resize(const Index nData, const Index xDim, const Index yDim)
{
isXExact_.setConstant(xDim, 0);
isFitPoint_.setConstant(nData, 0);
isXXCorr_.setIdentity(xDim, xDim);
isYYCorr_.setIdentity(xDim, xDim);
isYXCorr_.setConstant(yDim, xDim, 0);
isDataCorr_.setIdentity(nData, nData);
}
// test ////////////////////////////////////////////////////////////////////////
bool FitInterface::isFitPoint(const Index k) const
{
return (isFitPoint_[k] == 1);
}
bool FitInterface::isXExact(const Index i) const
{
return (isXExact_[i] == 1);
}
bool FitInterface::isXXCorrelated(const Index i1, const Index i2) const
{
return (isXXCorr_(i1, i2) == 1);
}
bool FitInterface::isYYCorrelated(const Index j1, const Index j2) const
{
return (isYYCorr_(j1, j2) == 1);
}
bool FitInterface::isYXCorrelated(const Index j, const Index i) const
{
return (isYXCorr_(j, i) == 1);
}
bool FitInterface::isDataCorrelated(const Index k1, const Index k2) const
{
return (isDataCorr_(k1, k2) == 1);
}

79
lib/FitInterface.hpp Normal file
View File

@ -0,0 +1,79 @@
/*
* FitInterface.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_FitInterface_hpp_
#define Latan_FitInterface_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Minimizer.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* base class for data fit *
******************************************************************************/
class FitInterface
{
public:
typedef Minimizer::Verbosity FitVerbosity;
public:
// constructors
FitInterface(void) = default;
FitInterface(const Index nData, const Index xDim, const Index yDim);
// destructor
virtual ~FitInterface(void) = default;
// 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 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;
void setFitInterface(const FitInterface &fitInterface);
void setNData(const Index nData);
void setXDim(const Index xDim);
void setYDim(const Index yDim);
void resize(const Index nData, const Index xDim, const Index yDim);
// test
bool isFitPoint(const Index k) 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;
private:
IVec isXExact_, isFitPoint_;
IMat isXXCorr_, isYYCorr_, isYXCorr_, isDataCorr_;
};
END_NAMESPACE
#endif // Latan_FitInterface_hpp_

118
lib/Function.cpp Normal file
View File

@ -0,0 +1,118 @@
/*
* Function.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Function.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* DoubleFunction implementation *
******************************************************************************/
const DoubleFunction::vecFunc DoubleFunction::nullFunction_ = nullptr;
// constructor /////////////////////////////////////////////////////////////////
DoubleFunction::DoubleFunction(const Index nArg, const vecFunc &f)
: buffer_(new DVec)
{
setFunction(f, nArg);
}
// access //////////////////////////////////////////////////////////////////////
Index DoubleFunction::getNArg(void) const
{
return buffer_->size();
}
void DoubleFunction::setFunction(const vecFunc &f, const Index nArg)
{
buffer_->resize(nArg);
f_ = f;
}
// error checking //////////////////////////////////////////////////////////////
void DoubleFunction::checkSize(const Index nPar) const
{
if (nPar != getNArg())
{
LATAN_ERROR(Size, "function argument vector has a wrong size (expected "
+ strFrom(getNArg()) + ", got " + strFrom(nPar)
+ ")");
}
}
// function call ///////////////////////////////////////////////////////////////
double DoubleFunction::operator()(const double *arg) const
{
return f_(arg);
}
double DoubleFunction::operator()(const DVec &arg) const
{
checkSize(arg.size());
return (*this)(arg.data());
}
double DoubleFunction::operator()(const std::vector<double> &arg) const
{
checkSize(static_cast<Index>(arg.size()));
return (*this)(arg.data());
}
double DoubleFunction::operator()(std::stack<double> &arg) const
{
for (Index i = 0; i < getNArg(); ++i)
{
if (arg.empty())
{
LATAN_ERROR(Size, "function argument stack is empty (expected "
+ strFrom(getNArg()) + "arguments, got " + strFrom(i)
+ ")");
}
(*buffer_)(getNArg() - i - 1) = arg.top();
arg.pop();
}
return (*this)(*buffer_);
}
double DoubleFunction::operator()(void) const
{
checkSize(0);
return (*this)(nullptr);
}
/******************************************************************************
* DoubleFunctionSample implementation *
******************************************************************************/
DSample DoubleFunctionSample::operator()(const DMatSample &arg) const
{
DSample result(arg.size());
FOR_STAT_ARRAY(arg, s)
{
result[s] = (*this)[s](arg[s]);
}
return result;
}

95
lib/Function.hpp Normal file
View File

@ -0,0 +1,95 @@
/*
* Function.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Function_hpp_
#define Latan_Function_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <functional>
#include <stack>
#include <vector>
BEGIN_NAMESPACE
/******************************************************************************
* Double function class *
******************************************************************************/
class DoubleFunction
{
private:
typedef std::function<double(const double *)> vecFunc;
public:
// constructor
explicit DoubleFunction(const Index nArg = 0,
const vecFunc &f = nullFunction_);
// destructor
virtual ~DoubleFunction(void) = default;
// access
virtual Index getNArg(void) const;
void setFunction(const vecFunc &f, const Index nArg);
// function call
virtual double operator()(const double *arg) const;
double operator()(const DVec &arg) const;
double operator()(const std::vector<double> &arg) const;
double operator()(std::stack<double> &arg) const;
double operator()(void) const;
template <typename... Ts>
double operator()(const double arg0, const Ts... args) const;
private:
// error checking
void checkSize(const Index nPar) const;
private:
std::shared_ptr<DVec> buffer_{nullptr};
vecFunc f_;
static const vecFunc nullFunction_;
};
template <typename... Ts>
double DoubleFunction::operator()(const double arg0, const Ts... args) const
{
static_assert(static_or<std::is_trivially_assignable<double, Ts>::value...>::value,
"DoubleFunction arguments are not compatible with double");
const double arg[] = {arg0, args...};
checkSize(sizeof...(args) + 1);
return (*this)(arg);
}
/******************************************************************************
* DoubleFunctionSample class *
******************************************************************************/
class DoubleFunctionSample: public Sample<DoubleFunction>
{
public:
// constructors
using Sample<DoubleFunction>::Sample;
// destructor
virtual ~DoubleFunctionSample(void) = default;
// function call
DSample operator()(const DMatSample &arg) const;
};
END_NAMESPACE
#endif // Latan_Function_hpp_

30
lib/Global.cpp Normal file
View File

@ -0,0 +1,30 @@
/*
* Global.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Global.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
const string Env::fullName = PACKAGE_STRING;
const string Env::name = PACKAGE_NAME;
const string Env::version = PACKAGE_VERSION;
const string Env::msgPrefix = "[" + strFrom(PACKAGE_NAME) + " v"
+ strFrom(PACKAGE_VERSION) + "] ";

189
lib/Global.hpp Normal file
View File

@ -0,0 +1,189 @@
/*
* Global.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Global_hpp_
#define Latan_Global_hpp_
// supress warning for the osbolete use of 'register' keyword in Eigen
#pragma GCC diagnostic ignored "-Wdeprecated-register"
#include <LatAnalyze/Eigen/Dense>
#include <memory>
#include <string>
#include <sstream>
#include <type_traits>
#include <cstdlib>
#define BEGIN_NAMESPACE namespace Latan {
#define END_NAMESPACE }
// macro utilities
#define unique_arg(...) __VA_ARGS__
#define DEBUG_VAR(x) std::cout << #x << "= " << x << std::endl
#define DEBUG_MAT(m) std::cout << #m << "=\n" << m << std::endl
// attribute to switch off unused warnings with gcc
#ifndef __GNUC__
#define __unused
#endif
// max length for paths
#define MAX_PATH_LENGTH 512u
// copy/assignement from Eigen expression
#define EIGEN_EXPR_CTOR(ctorName, Class, Base, EigenBase) \
template <typename Derived>\
ctorName(const Eigen::EigenBase<Derived> &m): Base(m) {}\
template<typename Derived>\
Class & operator=(const Eigen::EigenBase<Derived> &m)\
{\
this->Base::operator=(m);\
return *this;\
}
BEGIN_NAMESPACE
// Eigen type aliases //////////////////////////////////////////////////////////
const int dynamic = -1;
// array types
template <typename T, int nRow = dynamic, int nCol = dynamic>
using Array = Eigen::Array<T, nRow, nCol>;
// matrix types
template <typename T, int nRow = dynamic, int nCol = dynamic>
using Mat = Eigen::Matrix<T, nRow, nCol>;
typedef Mat<int> IMat;
typedef Mat<double> DMatBase;
// vector types
template <typename T>
using Vec = Mat<T, dynamic, 1>;
typedef Vec<int> IVec;
typedef Vec<double> DVec;
#define FOR_VEC(vec, i) for (Latan::Index i = 0; i < (vec).size(); ++i)
#define FOR_ARRAY(ar, i) FOR_VEC(ar, i)
// block types
template <typename Derived>
using Block = Eigen::Block<Derived>;
template <typename Derived>
using ConstBlock = const Eigen::Block<const Derived>;
template <typename Derived>
using Row = typename Derived::RowXpr;
template <typename Derived>
using ConstRow = typename Derived::ConstRowXpr;
template <typename Derived>
using Col = typename Derived::ColXpr;
template <typename Derived>
using ConstCol = typename Derived::ConstColXpr;
// map type
template <typename Derived>
using Map = Eigen::Map<Derived>;
template <typename Derived>
using ConstMap = Eigen::Map<const Derived>;
// Index type //////////////////////////////////////////////////////////////////
typedef DMatBase::Index Index;
// Placeholder type ////////////////////////////////////////////////////////////
struct PlaceHolder {};
extern PlaceHolder _;
// Type utilities //////////////////////////////////////////////////////////////
// pointer type test
template <typename Derived, typename Base>
inline bool isDerivedFrom(const Base *pt)
{
return (dynamic_cast<const Derived *>(pt) != nullptr);
}
// static logical or
template <bool... b>
struct static_or;
template <bool... tail>
struct static_or<true, tail...> : static_or<tail...> {};
template <bool... tail>
struct static_or<false, tail...> : std::false_type {};
template <>
struct static_or<> : std::true_type {};
// Environment /////////////////////////////////////////////////////////////////
namespace Env
{
extern const std::string fullName;
extern const std::string name;
extern const std::string version;
extern const std::string msgPrefix;
}
// String conversions //////////////////////////////////////////////////////////
template <typename T>
inline T strTo(const std::string &str)
{
T buf;
std::istringstream stream(str);
stream >> buf;
return buf;
}
// optimized specializations
template <>
inline float strTo<float>(const std::string &str)
{
return strtof(str.c_str(), (char **)NULL);
}
template <>
inline double strTo<double>(const std::string &str)
{
return strtod(str.c_str(), (char **)NULL);
}
template <>
inline int strTo<int>(const std::string &str)
{
return (int)(strtol(str.c_str(), (char **)NULL, 10));
}
template <typename T>
inline std::string strFrom(const T x)
{
std::ostringstream stream;
stream << x;
return stream.str();
}
END_NAMESPACE
#include <LatAnalyze/Exceptions.hpp>
#endif // Latan_Global_hpp_

49
lib/IoObject.hpp Normal file
View File

@ -0,0 +1,49 @@
/*
* IoObject.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_IoObject_hpp_
#define Latan_IoObject_hpp_
#include <LatAnalyze/Global.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* Abstract class for IO objects *
******************************************************************************/
class IoObject
{
public:
enum class IoType
{
noType = 0,
dMat = 1,
dMatSample = 2,
rgState = 3
};
public:
// destructor
virtual ~IoObject(void) = default;
// access
virtual IoType getType(void) const = 0;
};
END_NAMESPACE
#endif // Latan_IoObject_hpp_

1
lib/LatAnalyze Symbolic link
View File

@ -0,0 +1 @@
.

85
lib/Makefile.am Normal file
View File

@ -0,0 +1,85 @@
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
else
if CXX_INTEL
COM_CXXFLAGS = -Wall
endif
endif
AM_LFLAGS = -olex.yy.c
AM_YFLAGS = -d
include eigen_files.mk
nobase_dist_pkginclude_HEADERS = $(eigen_files)
BUILT_SOURCES = AsciiParser.hpp MathParser.hpp
lib_LTLIBRARIES = libLatAnalyze.la
libLatAnalyze_la_SOURCES =\
AsciiFile.cpp \
AsciiParser.ypp \
AsciiLexer.lpp \
Chi2Function.cpp \
CompiledFunction.cpp \
CompiledModel.cpp \
Exceptions.cpp \
Function.cpp \
Global.cpp \
includes.hpp \
File.cpp \
FitInterface.cpp \
Mat.cpp \
Math.cpp \
MathInterpreter.cpp \
MathParser.ypp \
MathLexer.lpp \
MatSample.cpp \
Minimizer.cpp \
Model.cpp \
Plot.cpp \
RandGen.cpp \
XYSampleData.cpp \
XYStatData.cpp \
../config.h
libLatAnalyze_ladir = $(pkgincludedir)
libLatAnalyze_la_HEADERS =\
AsciiFile.hpp \
Chi2Function.hpp \
CompiledFunction.hpp \
CompiledModel.hpp \
Dataset.hpp \
FitInterface.hpp \
Function.hpp \
Global.hpp \
File.hpp \
IoObject.hpp \
Mat.hpp \
Math.hpp \
MathInterpreter.hpp \
MatSample.hpp \
Minimizer.hpp \
Model.hpp \
ParserState.hpp \
Plot.hpp \
RandGen.hpp \
StatArray.hpp \
XYSampleData.hpp \
XYStatData.hpp
if HAVE_MINUIT
libLatAnalyze_la_SOURCES += MinuitMinimizer.cpp
libLatAnalyze_la_HEADERS += MinuitMinimizer.hpp
endif
libLatAnalyze_la_CFLAGS = $(COM_CFLAGS)
libLatAnalyze_la_CXXFLAGS = $(COM_CXXFLAGS)
ACLOCAL_AMFLAGS = -I .buildutils/m4

37
lib/Mat.cpp Normal file
View File

@ -0,0 +1,37 @@
/*
* Mat.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Mat.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* DMat implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
DMat::DMat(const Index nRow, const Index nCol)
: Base(nRow, nCol)
{}
IoObject::IoType DMat::getType(void) const
{
return IoType::dMat;
}

53
lib/Mat.hpp Normal file
View File

@ -0,0 +1,53 @@
/*
* Mat.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Mat_hpp_
#define Latan_Mat_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/IOObject.hpp>
#define FOR_MAT(mat, i, j) \
for (Latan::Index j = 0; j < mat.cols(); ++j)\
for (Latan::Index i = 0; i < mat.rows(); ++i)
BEGIN_NAMESPACE
/******************************************************************************
* double matrix type *
******************************************************************************/
class DMat: public DMatBase, public IoObject
{
private:
typedef DMatBase Base;
public:
// constructors
DMat(void) = default;
DMat(const Index nRow, const Index nCol);
EIGEN_EXPR_CTOR(DMat, DMat, Base, MatrixBase)
// destructor
virtual ~DMat(void) = default;
// IO
virtual IoType getType(void) const;
};
END_NAMESPACE
#endif // Latan_Mat_hpp_

117
lib/MatSample.cpp Normal file
View File

@ -0,0 +1,117 @@
/*
* MatSample.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/MatSample.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* DMatSample implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
DMatSample::DMatSample(const Index nSample, const Index nRow,
const Index nCol)
: Sample<DMat>(nSample)
{
resizeMat(nRow, nCol);
}
DMatSample::DMatSample(ConstBlock &sampleBlock)
: DMatSample(sampleBlock.getSample().size(), sampleBlock.getNRow(),
sampleBlock.getNCol())
{
const DMatSample &sample = sampleBlock.getSample();
this->resize(sample.size());
FOR_STAT_ARRAY(*this, s)
{
(*this)[s] = sample[s].block(sampleBlock.getStartRow(),
sampleBlock.getStartCol(),
sampleBlock.getNRow(),
sampleBlock.getNCol());
}
}
DMatSample::DMatSample(ConstBlock &&sampleBlock)
: DMatSample(sampleBlock)
{}
// assignement operator ////////////////////////////////////////////////////////
DMatSample & DMatSample::operator=(Block &sampleBlock)
{
DMatSample tmp(sampleBlock);
this->swap(tmp);
return *this;
}
DMatSample & DMatSample::operator=(Block &&sampleBlock)
{
*this = sampleBlock;
return *this;
}
DMatSample & DMatSample::operator=(ConstBlock &sampleBlock)
{
DMatSample tmp(sampleBlock);
this->swap(tmp);
return *this;
}
DMatSample & DMatSample::operator=(ConstBlock &&sampleBlock)
{
*this = sampleBlock;
return *this;
}
// block access ////////////////////////////////////////////////////////////////
DMatSample::ConstBlock DMatSample::block(const Index i, const Index j,
const Index nRow,
const Index nCol) const
{
return ConstBlock(*this, i, j, nRow, nCol);
}
DMatSample::Block DMatSample::block(const Index i, const Index j,
const Index nRow, const Index nCol)
{
return Block(*this, i, j, nRow, nCol);
}
// resize all matrices /////////////////////////////////////////////////////////
void DMatSample::resizeMat(const Index nRow, const Index nCol)
{
FOR_STAT_ARRAY(*this, s)
{
(*this)[s].resize(nRow, nCol);
}
}
// IO type /////////////////////////////////////////////////////////////////////
IoObject::IoType DMatSample::getType(void) const
{
return IoType::dMatSample;
}

181
lib/MatSample.hpp Normal file
View File

@ -0,0 +1,181 @@
/*
* MatSample.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_MatSample_hpp_
#define Latan_MatSample_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/StatArray.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* DMat sample class *
******************************************************************************/
class DMatSample: public Sample<DMat>, public IoObject
{
public:
// block type template
template <class S>
class BlockTemplate
{
private:
typedef typename std::remove_const<S>::type NonConstType;
public:
// constructors
BlockTemplate(S &sample, const Index i, const Index j, const Index nRow,
const Index nCol);
BlockTemplate(BlockTemplate<NonConstType> &b);
BlockTemplate(BlockTemplate<NonConstType> &&b);
// destructor
~BlockTemplate(void) = default;
// access
S & getSample(void);
const S & getSample(void) const;
Index getStartRow(void) const;
Index getStartCol(void) const;
Index getNRow(void) const;
Index getNCol(void) const;
// assignement operators
BlockTemplate<S> & operator=(const S &sample);
BlockTemplate<S> & operator=(const S &&sample);
private:
S &sample_;
const Index i_, j_, nRow_, nCol_;
};
// block types
typedef BlockTemplate<DMatSample> Block;
typedef const BlockTemplate<const DMatSample> ConstBlock;
public:
// constructors
using Sample<DMat>::Sample;
DMatSample(void) = default;
DMatSample(const Index nSample, const Index nRow, const Index nCol);
DMatSample(ConstBlock &sampleBlock);
DMatSample(ConstBlock &&sampleBlock);
// destructor
virtual ~DMatSample(void) = default;
// assignement operator
DMatSample & operator=(Block &sampleBlock);
DMatSample & operator=(Block &&sampleBlock);
DMatSample & operator=(ConstBlock &sampleBlock);
DMatSample & operator=(ConstBlock &&sampleBlock);
// block access
ConstBlock block(const Index i, const Index j, const Index nRow,
const Index nCol) const;
Block block(const Index i, const Index j, const Index nRow,
const Index nCol);
// resize all matrices
void resizeMat(const Index nRow, const Index nCol);
// IO type
virtual IoType getType(void) const;
};
/******************************************************************************
* Block template implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
template <class S>
DMatSample::BlockTemplate<S>::BlockTemplate(S &sample, const Index i,
const Index j, const Index nRow,
const Index nCol)
: sample_(sample)
, i_(i)
, j_(j)
, nRow_(nRow)
, nCol_(nCol)
{}
template <class S>
DMatSample::BlockTemplate<S>::BlockTemplate(BlockTemplate<NonConstType> &b)
: sample_(b.getSample())
, i_(b.getStartRow())
, j_(b.getStartCol())
, nRow_(b.getNRow())
, nCol_(b.getNCol())
{}
template <class S>
DMatSample::BlockTemplate<S>::BlockTemplate(BlockTemplate<NonConstType> &&b)
: BlockTemplate(b)
{}
// access //////////////////////////////////////////////////////////////////////
template <class S>
S & DMatSample::BlockTemplate<S>::getSample(void)
{
return sample_;
}
template <class S>
const S & DMatSample::BlockTemplate<S>::getSample(void) const
{
return sample_;
}
template <class S>
Index DMatSample::BlockTemplate<S>::getStartRow(void) const
{
return i_;
}
template <class S>
Index DMatSample::BlockTemplate<S>::getStartCol(void) const
{
return j_;
}
template <class S>
Index DMatSample::BlockTemplate<S>::getNRow(void) const
{
return nRow_;
}
template <class S>
Index DMatSample::BlockTemplate<S>::getNCol(void) const
{
return nCol_;
}
// assignement operators ///////////////////////////////////////////////////////
template <class S>
DMatSample::BlockTemplate<S> &
DMatSample::BlockTemplate<S>::operator=(const S &sample)
{
FOR_STAT_ARRAY(sample_, s)
{
sample_[s].block(i_, j_, nRow_, nCol_) = sample[s];
}
return *this;
}
template <class S>
DMatSample::BlockTemplate<S> &
DMatSample::BlockTemplate<S>::operator=(const S &&sample)
{
*this = sample;
return *this;
}
END_NAMESPACE
#endif // Latan_MatSample_hpp_

152
lib/Math.cpp Normal file
View File

@ -0,0 +1,152 @@
/*
* Math.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Math.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* Standard C functions *
******************************************************************************/
#define DEF_STD_FUNC_1ARG(name) \
auto name##VecFunc = [](const double arg[1]){return (name)(arg[0]);};\
DoubleFunction STDMATH_NAMESPACE::name(1, name##VecFunc);
#define DEF_STD_FUNC_2ARG(name) \
auto name##VecFunc = [](const double arg[2]){return (name)(arg[0], arg[1]);};\
DoubleFunction STDMATH_NAMESPACE::name(2, name##VecFunc);
// Trigonometric functions
DEF_STD_FUNC_1ARG(cos)
DEF_STD_FUNC_1ARG(sin)
DEF_STD_FUNC_1ARG(tan)
DEF_STD_FUNC_1ARG(acos)
DEF_STD_FUNC_1ARG(asin)
DEF_STD_FUNC_1ARG(atan)
DEF_STD_FUNC_2ARG(atan2)
// Hyperbolic functions
DEF_STD_FUNC_1ARG(cosh)
DEF_STD_FUNC_1ARG(sinh)
DEF_STD_FUNC_1ARG(tanh)
DEF_STD_FUNC_1ARG(acosh)
DEF_STD_FUNC_1ARG(asinh)
DEF_STD_FUNC_1ARG(atanh)
// Exponential and logarithmic functions
DEF_STD_FUNC_1ARG(exp)
DEF_STD_FUNC_1ARG(log)
DEF_STD_FUNC_1ARG(log10)
DEF_STD_FUNC_1ARG(exp2)
DEF_STD_FUNC_1ARG(expm1)
DEF_STD_FUNC_1ARG(log1p)
DEF_STD_FUNC_1ARG(log2)
// Power functions
DEF_STD_FUNC_2ARG(pow)
DEF_STD_FUNC_1ARG(sqrt)
DEF_STD_FUNC_1ARG(cbrt)
DEF_STD_FUNC_2ARG(hypot)
// Error and gamma functions
DEF_STD_FUNC_1ARG(erf)
DEF_STD_FUNC_1ARG(erfc)
DEF_STD_FUNC_1ARG(tgamma)
DEF_STD_FUNC_1ARG(lgamma)
// Rounding and remainder functions
DEF_STD_FUNC_1ARG(ceil)
DEF_STD_FUNC_1ARG(floor)
DEF_STD_FUNC_2ARG(fmod)
DEF_STD_FUNC_1ARG(trunc)
DEF_STD_FUNC_1ARG(round)
DEF_STD_FUNC_1ARG(rint)
DEF_STD_FUNC_1ARG(nearbyint)
DEF_STD_FUNC_2ARG(remainder)
// Minimum, maximum, difference functions
DEF_STD_FUNC_2ARG(fdim)
DEF_STD_FUNC_2ARG(fmax)
DEF_STD_FUNC_2ARG(fmin)
// Absolute value
DEF_STD_FUNC_1ARG(fabs)
#define ADD_FUNC(func) fTable[#func] = &STDMATH_NAMESPACE::func
void STDMATH_NAMESPACE::addStdMathFunc(FunctionTable &fTable)
{
// Trigonometric functions
ADD_FUNC(cos);
ADD_FUNC(sin);
ADD_FUNC(tan);
ADD_FUNC(acos);
ADD_FUNC(asin);
ADD_FUNC(atan);
ADD_FUNC(atan2);
// Hyperbolic functions
ADD_FUNC(cosh);
ADD_FUNC(sinh);
ADD_FUNC(tanh);
ADD_FUNC(acosh);
ADD_FUNC(asinh);
ADD_FUNC(atanh);
// Exponential and logarithmic functions
ADD_FUNC(exp);
ADD_FUNC(log);
ADD_FUNC(log10);
ADD_FUNC(exp2);
ADD_FUNC(expm1);
ADD_FUNC(log1p);
ADD_FUNC(log2);
// Power functions
ADD_FUNC(pow);
ADD_FUNC(sqrt);
ADD_FUNC(cbrt);
ADD_FUNC(hypot);
// Error and gamma functions
ADD_FUNC(erf);
ADD_FUNC(erfc);
ADD_FUNC(tgamma);
ADD_FUNC(lgamma);
// Rounding and remainder functions
ADD_FUNC(ceil);
ADD_FUNC(floor);
ADD_FUNC(fmod);
ADD_FUNC(trunc);
ADD_FUNC(round);
ADD_FUNC(rint);
ADD_FUNC(nearbyint);
ADD_FUNC(remainder);
// Minimum, maximum, difference functions
ADD_FUNC(fdim);
ADD_FUNC(fmax);
ADD_FUNC(fmin);
// Absolute value
ADD_FUNC(fabs);
}

104
lib/Math.hpp Normal file
View File

@ -0,0 +1,104 @@
/*
* Math.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Math_hpp_
#define Latan_Math_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/MathInterpreter.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* Standard C functions *
******************************************************************************/
#define STDMATH_NAMESPACE StdMath
#define DECL_STD_FUNC(name) \
namespace STDMATH_NAMESPACE\
{\
extern DoubleFunction name;\
}
// Trigonometric functions
DECL_STD_FUNC(cos)
DECL_STD_FUNC(sin)
DECL_STD_FUNC(tan)
DECL_STD_FUNC(acos)
DECL_STD_FUNC(asin)
DECL_STD_FUNC(atan)
DECL_STD_FUNC(atan2)
// Hyperbolic functions
DECL_STD_FUNC(cosh)
DECL_STD_FUNC(sinh)
DECL_STD_FUNC(tanh)
DECL_STD_FUNC(acosh)
DECL_STD_FUNC(asinh)
DECL_STD_FUNC(atanh)
// Exponential and logarithmic functions
DECL_STD_FUNC(exp)
DECL_STD_FUNC(log)
DECL_STD_FUNC(log10)
DECL_STD_FUNC(exp2)
DECL_STD_FUNC(expm1)
DECL_STD_FUNC(log1p)
DECL_STD_FUNC(log2)
// Power functions
DECL_STD_FUNC(pow)
DECL_STD_FUNC(sqrt)
DECL_STD_FUNC(cbrt)
DECL_STD_FUNC(hypot)
// Error and gamma functions
DECL_STD_FUNC(erf)
DECL_STD_FUNC(erfc)
DECL_STD_FUNC(tgamma)
DECL_STD_FUNC(lgamma)
// Rounding and remainder functions
DECL_STD_FUNC(ceil)
DECL_STD_FUNC(floor)
DECL_STD_FUNC(fmod)
DECL_STD_FUNC(trunc)
DECL_STD_FUNC(round)
DECL_STD_FUNC(rint)
DECL_STD_FUNC(nearbyint)
DECL_STD_FUNC(remainder)
// Minimum, maximum, difference functions
DECL_STD_FUNC(fdim)
DECL_STD_FUNC(fmax)
DECL_STD_FUNC(fmin)
// Absolute value
DECL_STD_FUNC(fabs)
// Add standard math functions to a table for the math compiler
namespace STDMATH_NAMESPACE
{
void addStdMathFunc(FunctionTable &fTable);
}
END_NAMESPACE
#endif // Latan_Math_hpp_

502
lib/MathInterpreter.cpp Normal file
View File

@ -0,0 +1,502 @@
/*
* MathInterpreter.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/MathInterpreter.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* Instruction set *
******************************************************************************/
#define CODE_WIDTH 6
#define CODE_MOD setw(CODE_WIDTH) << left
// Instruction operator ////////////////////////////////////////////////////////
ostream &Latan::operator<<(ostream& out, const Instruction& ins)
{
ins.print(out);
return out;
}
// Push constructors ///////////////////////////////////////////////////////////
Push::Push(const double val)
: type_(ArgType::Constant)
, val_(val)
, name_("")
{}
Push::Push(const string &name)
: type_(ArgType::Variable)
, val_(0.0)
, name_(name)
{}
// Push execution //////////////////////////////////////////////////////////////
void Push::operator()(RunContext &context) const
{
if (type_ == ArgType::Constant)
{
context.dStack.push(val_);
}
else
{
try
{
context.dStack.push(context.vTable.at(name_));
}
catch (out_of_range)
{
LATAN_ERROR(Range, "unknown variable '" + name_ + "'");
}
}
context.insIndex++;
}
// Push print //////////////////////////////////////////////////////////////////
void Push::print(ostream &out) const
{
out << CODE_MOD << "push";
if (type_ == ArgType::Constant)
{
out << CODE_MOD << val_;
}
else
{
out << CODE_MOD << name_;
}
}
// Pop constructor /////////////////////////////////////////////////////////////
Pop::Pop(const string &name)
: name_(name)
{}
// Pop execution ///////////////////////////////////////////////////////////////
void Pop::operator()(RunContext &context) const
{
if (!name_.empty())
{
context.vTable[name_] = context.dStack.top();
}
context.dStack.pop();
context.insIndex++;
}
// Pop print ///////////////////////////////////////////////////////////////////
void Pop::print(ostream &out) const
{
out << CODE_MOD << "pop" << CODE_MOD << name_;
}
// Store constructor ///////////////////////////////////////////////////////////
Store::Store(const string &name)
: name_(name)
{}
// Store execution /////////////////////////////////////////////////////////////
void Store::operator()(RunContext &context) const
{
if (!name_.empty())
{
context.vTable[name_] = context.dStack.top();
}
context.insIndex++;
}
// Store print /////////////////////////////////////////////////////////////////
void Store::print(ostream &out) const
{
out << CODE_MOD << "store" << CODE_MOD << name_;
}
// Call constructor ////////////////////////////////////////////////////////////
Call::Call(const string &name)
: name_(name)
{}
// Call execution //////////////////////////////////////////////////////////////
void Call::operator()(RunContext &context) const
{
try
{
context.dStack.push((*context.fTable.at(name_))(context.dStack));
}
catch (out_of_range)
{
LATAN_ERROR(Range, "unknown function '" + name_ + "'");
}
context.insIndex++;
}
// Call print //////////////////////////////////////////////////////////////////
void Call::print(ostream &out) const
{
out << CODE_MOD << "call" << CODE_MOD << name_;
}
// Math operations /////////////////////////////////////////////////////////////
#define DEF_OP(name, nArg, exp, insName)\
void name::operator()(RunContext &context) const\
{\
double x[nArg];\
for (int i = 0; i < nArg; ++i)\
{\
x[nArg-1-i] = context.dStack.top();\
context.dStack.pop();\
}\
context.dStack.push(exp);\
context.insIndex++;\
}\
void name::print(ostream &out) const\
{\
out << CODE_MOD << insName;\
}
DEF_OP(Neg, 1, -x[0], "neg")
DEF_OP(Add, 2, x[0] + x[1], "add")
DEF_OP(Sub, 2, x[0] - x[1], "sub")
DEF_OP(Mul, 2, x[0]*x[1], "mul")
DEF_OP(Div, 2, x[0]/x[1], "div")
DEF_OP(Pow, 2, pow(x[0],x[1]), "pow")
/******************************************************************************
* ExprNode implementation *
******************************************************************************/
// ExprNode constructors ///////////////////////////////////////////////////////
ExprNode::ExprNode(const string &name)
: name_(name)
, parent_(nullptr)
{}
// ExprNode access /////////////////////////////////////////////////////////////
const string &ExprNode::getName(void) const
{
return name_;
}
Index ExprNode::getNArg(void) const
{
return static_cast<Index>(arg_.size());
}
const ExprNode * ExprNode::getParent(void) const
{
return parent_;
}
Index ExprNode::getLevel(void) const
{
if (getParent())
{
return getParent()->getLevel() + 1;
}
else
{
return 0;
}
}
void ExprNode::setName(const std::string &name)
{
name_ = name;
}
void ExprNode::pushArg(ExprNode *node)
{
if (node)
{
node->parent_ = this;
arg_.push_back(unique_ptr<ExprNode>(node));
}
}
// ExprNode operators //////////////////////////////////////////////////////////
const ExprNode &ExprNode::operator[](const Index i) const
{
return *arg_[static_cast<unsigned int>(i)];
}
ostream &Latan::operator<<(ostream &out, const ExprNode &n)
{
Index level = n.getLevel();
for (Index i = 0; i <= level; ++i)
{
if (i == level)
{
out << "_";
}
else if (i == level - 1)
{
out << "|";
}
else
{
out << " ";
}
}
out << " " << n.getName() << endl;
for (Index i = 0; i < n.getNArg(); ++i)
{
out << n[i];
}
return out;
}
#define PUSH_INS(program, type, ...) \
program.push_back(unique_ptr<type>(new type(__VA_ARGS__)))
// VarNode compile /////////////////////////////////////////////////////////////
void VarNode::compile(Program &program) const
{
PUSH_INS(program, Push, getName());
}
// CstNode compile /////////////////////////////////////////////////////////////
void CstNode::compile(Program &program) const
{
PUSH_INS(program, Push, strTo<double>(getName()));
}
// SemicolonNode compile ///////////////////////////////////////////////////////
void SemicolonNode::compile(Program &program) const
{
auto &n = *this;
for (Index i = 0; i < getNArg(); ++i)
{
bool isAssign = isDerivedFrom<AssignNode>(&n[i]);
bool isSemiColumn = isDerivedFrom<SemicolonNode>(&n[i]);
bool isKeyword = isDerivedFrom<KeywordNode>(&n[i]);
if (isAssign||isSemiColumn||isKeyword)
{
n[i].compile(program);
}
}
}
// AssignNode compile //////////////////////////////////////////////////////////
void AssignNode::compile(Program &program) const
{
auto &n = *this;
if (isDerivedFrom<VarNode>(&n[0]))
{
bool hasSemicolonParent = isDerivedFrom<SemicolonNode>(getParent());
n[1].compile(program);
if (hasSemicolonParent)
{
PUSH_INS(program, Pop, n[0].getName());
}
else
{
PUSH_INS(program, Store, n[0].getName());
}
}
else
{
LATAN_ERROR(Compilation, "invalid LHS for '='");
}
}
// MathOpNode compile //////////////////////////////////////////////////////////
#define IFNODE(name, nArg) if ((n.getName() == (name))&&(n.getNArg() == nArg))
#define ELIFNODE(name, nArg) else IFNODE(name, nArg)
#define ELSE else
void MathOpNode::compile(Program &program) const
{
auto &n = *this;
for (Index i = 0; i < n.getNArg(); ++i)
{
n[i].compile(program);
}
IFNODE("-", 1) PUSH_INS(program, Neg,);
ELIFNODE("+", 2) PUSH_INS(program, Add,);
ELIFNODE("-", 2) PUSH_INS(program, Sub,);
ELIFNODE("*", 2) PUSH_INS(program, Mul,);
ELIFNODE("/", 2) PUSH_INS(program, Div,);
ELIFNODE("^", 2) PUSH_INS(program, Pow,);
ELSE LATAN_ERROR(Compilation, "unknown operator '" + getName() + "'");
}
// FuncNode compile ////////////////////////////////////////////////////////////
void FuncNode::compile(Program &program) const
{
auto &n = *this;
for (Index i = 0; i < n.getNArg(); ++i)
{
n[i].compile(program);
}
PUSH_INS(program, Call, getName());
}
// ReturnNode compile ////////////////////////////////////////////////////////////
void ReturnNode::compile(Program &program) const
{
auto &n = *this;
n[0].compile(program);
program.push_back(nullptr);
}
/******************************************************************************
* MathInterpreter implementation *
******************************************************************************/
// MathParserState constructor /////////////////////////////////////////////////
MathInterpreter::MathParserState::MathParserState
(istream *stream, string *name, std::unique_ptr<ExprNode> *data)
: ParserState<std::unique_ptr<ExprNode>>(stream, name, data)
{
initScanner();
}
// MathParserState destructor //////////////////////////////////////////////////
MathInterpreter::MathParserState::~MathParserState(void)
{
destroyScanner();
}
// constructors ////////////////////////////////////////////////////////////////
MathInterpreter::MathInterpreter(const std::string &code)
: codeName_("<string>")
{
setCode(code);
}
// access //////////////////////////////////////////////////////////////////////
const Instruction * MathInterpreter::operator[](const Index i) const
{
return program_[static_cast<unsigned int>(i)].get();
}
const ExprNode * MathInterpreter::getAST(void) const
{
return root_.get();
}
void MathInterpreter::push(const Instruction *i)
{
program_.push_back(unique_ptr<const Instruction>(i));
}
// initialization //////////////////////////////////////////////////////////////
void MathInterpreter::setCode(const std::string &code)
{
if (status_)
{
reset();
}
code_.reset(new stringstream(code));
codeName_ = "<string>";
state_.reset(new MathParserState(code_.get(), &codeName_, &root_));
program_.clear();
status_ = Status::initialised;
}
void MathInterpreter::reset(void)
{
code_.reset();
codeName_ = "<no_code>";
state_.reset();
root_.reset();
program_.clear();
status_ = 0;
}
// parser //////////////////////////////////////////////////////////////////////
// Bison/Flex parser declaration
int _math_parse(MathInterpreter::MathParserState *state);
void MathInterpreter::parse(void)
{
_math_parse(state_.get());
}
// interpreter /////////////////////////////////////////////////////////////////
void MathInterpreter::compile(void)
{
bool gotReturn = false;
if (!(status_ & Status::parsed))
{
parse();
status_ |= Status::parsed;
status_ -= status_ & Status::compiled;
}
if (root_)
{
root_->compile(program_);
for (unsigned int i = 0; i < program_.size(); ++i)
{
if (!program_[i])
{
gotReturn = true;
program_.resize(i);
program_.shrink_to_fit();
break;
}
}
}
if (!root_||!gotReturn)
{
LATAN_ERROR(Syntax, "expected 'return' in program '" + codeName_ + "'");
}
status_ |= Status::compiled;
}
// execution ///////////////////////////////////////////////////////////////////
void MathInterpreter::operator()(RunContext &context)
{
if (!(status_ & Status::compiled))
{
compile();
}
execute(context);
}
void MathInterpreter::execute(RunContext &context) const
{
context.insIndex = 0;
while (context.insIndex != program_.size())
{
(*(program_[context.insIndex]))(context);
}
}
// IO //////////////////////////////////////////////////////////////////////////
ostream &Latan::operator<<(ostream &out, const MathInterpreter &program)
{
for (unsigned int i = 0; i < program.program_.size(); ++i)
{
out << *(program.program_[i]) << endl;
}
return out;
}

278
lib/MathInterpreter.hpp Normal file
View File

@ -0,0 +1,278 @@
/*
* MathInterpreter.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_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>
#define MAXIDLENGTH 256
BEGIN_NAMESPACE
/******************************************************************************
* Instruction classes *
******************************************************************************/
typedef std::map<std::string, double> VarTable;
typedef std::map<std::string, DoubleFunction *> FunctionTable;
struct RunContext
{
unsigned int insIndex;
std::stack<double> dStack;
VarTable vTable;
FunctionTable fTable;
};
// Abstract base
class Instruction
{
public:
// destructor
virtual ~Instruction(void) = default;
// instruction execution
virtual void operator()(RunContext &context) const = 0;
friend std::ostream & operator<<(std::ostream &out, const Instruction &ins);
private:
virtual void print(std::ostream &out) const = 0;
};
std::ostream & operator<<(std::ostream &out, const Instruction &ins);
// Instruction container
typedef std::vector<std::unique_ptr<const Instruction>> Program;
// Push
class Push: public Instruction
{
private:
enum class ArgType
{
Constant = 0,
Variable = 1
};
public:
//constructors
explicit Push(const double val);
explicit Push(const std::string &name);
// instruction execution
virtual void operator()(RunContext &context) const;
private:
virtual void print(std::ostream& out) const;
private:
ArgType type_;
double val_;
std::string name_;
};
// Pop
class Pop: public Instruction
{
public:
//constructor
explicit Pop(const std::string &name);
// instruction execution
virtual void operator()(RunContext &context) const;
private:
virtual void print(std::ostream& out) const;
private:
std::string name_;
};
// Store
class Store: public Instruction
{
public:
//constructor
explicit Store(const std::string &name);
// instruction execution
virtual void operator()(RunContext &context) const;
private:
virtual void print(std::ostream& out) const;
private:
std::string name_;
};
// Call function
class Call: public Instruction
{
public:
//constructor
explicit Call(const std::string &name);
// instruction execution
virtual void operator()(RunContext &context) const;
private:
virtual void print(std::ostream& out) const;
private:
std::string name_;
};
// Floating point operations
#define DECL_OP(name)\
class name: public Instruction\
{\
public:\
virtual void operator()(RunContext &context) const;\
private:\
virtual void print(std::ostream &out) const;\
}
DECL_OP(Neg);
DECL_OP(Add);
DECL_OP(Sub);
DECL_OP(Mul);
DECL_OP(Div);
DECL_OP(Pow);
/******************************************************************************
* Expression node classes *
******************************************************************************/
class ExprNode
{
public:
// constructors
explicit ExprNode(const std::string &name);
// destructor
virtual ~ExprNode() = default;
// access
const std::string& getName(void) const;
Index getNArg(void) const;
const ExprNode * getParent(void) const;
Index getLevel(void) const;
void setName(const std::string &name);
void pushArg(ExprNode *node);
// operator
const ExprNode &operator[](const Index i) const;
// compile
virtual void compile(Program &program) const = 0;
private:
std::string name_;
std::vector<std::unique_ptr<ExprNode>> arg_;
const ExprNode * parent_;
};
std::ostream &operator<<(std::ostream &out, const ExprNode &n);
#define DECL_NODE(base, name) \
class name: public base\
{\
public:\
using base::base;\
virtual void compile(Program &program) const;\
}
DECL_NODE(ExprNode, VarNode);
DECL_NODE(ExprNode, CstNode);
DECL_NODE(ExprNode, SemicolonNode);
DECL_NODE(ExprNode, AssignNode);
DECL_NODE(ExprNode, MathOpNode);
DECL_NODE(ExprNode, FuncNode);
class KeywordNode: public ExprNode
{
public:
using ExprNode::ExprNode;
virtual void compile(Program &program) const = 0;
};
DECL_NODE(KeywordNode, ReturnNode);
/******************************************************************************
* Interpreter class *
******************************************************************************/
class MathInterpreter
{
public:
// parser state
class MathParserState: public ParserState<std::unique_ptr<ExprNode>>
{
public:
// constructor
explicit MathParserState(std::istream *stream, std::string *name,
std::unique_ptr<ExprNode> *data);
// destructor
virtual ~MathParserState(void);
private:
// allocation/deallocation functions defined in MathLexer.lpp
virtual void initScanner(void);
virtual void destroyScanner(void);
};
private:
// status flags
class Status
{
public:
enum
{
none = 0,
initialised = 1 << 0,
parsed = 1 << 1,
compiled = 1 << 2
};
};
public:
// constructors
MathInterpreter(void) = default;
MathInterpreter(const std::string &code);
// destructor
~MathInterpreter(void) = default;
// access
const Instruction * operator[](const Index i) const;
const ExprNode * getAST(void) const;
// initialization
void setCode(const std::string &code);
// interpreter
void compile(void);
// execution
void operator()(RunContext &context);
// IO
friend std::ostream & operator<<(std::ostream &out,
const MathInterpreter &program);
private:
// initialization
void reset(void);
// access
void push(const Instruction *i);
// parser
void parse(void);
// interpreter
void compileNode(const ExprNode &node);
// execution
void execute(RunContext &context) const;
private:
std::unique_ptr<std::istream> code_{nullptr};
std::string codeName_{"<no_code>"};
std::unique_ptr<MathParserState> state_{nullptr};
std::unique_ptr<ExprNode> root_{nullptr};
Program program_;
unsigned int status_{Status::none};
};
std::ostream & operator<<(std::ostream &out, const MathInterpreter &program);
END_NAMESPACE
#endif // Latan_MathInterpreter_hpp_

93
lib/MathLexer.lpp Normal file
View File

@ -0,0 +1,93 @@
/*
* MathLexer.lpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/>.
*/
%option reentrant
%option prefix="_math_"
%option bison-bridge
%option bison-locations
%option noyywrap
%option yylineno
%{
#include <iostream>
#include <LatAnalyze/MathInterpreter.hpp>
#include <LatAnalyze/MathParser.hpp>
#pragma GCC diagnostic ignored "-Wsign-compare"
#pragma GCC diagnostic ignored "-Wunused-function"
#pragma GCC diagnostic ignored "-Wunused-parameter"
using namespace std;
using namespace Latan;
#define YY_EXTRA_TYPE MathInterpreter::MathParserState *
#define YY_USER_ACTION \
yylloc->first_line = yylloc->last_line = yylineno;\
yylloc->first_column = yylloc->last_column + 1;\
yylloc->last_column = yylloc->first_column + yyleng - 1;
#define YY_INPUT(buf, result, max_size) \
{ \
(*yyextra->stream).read(buf, max_size);\
result = (*yyextra->stream).gcount();\
}
#define YY_DEBUG 0
#if (YY_DEBUG == 1)
#define RET(var) cout << (var) << "(" << yytext << ")" << endl; return (var)
#define RETTOK(tok) cout << #tok << "(" << yytext << ")" << endl; return tok
#else
#define RET(var) return (var)
#define RETTOK(tok) return tok
#endif
%}
DIGIT [0-9]
ALPHA [a-zA-Z_]
FLOAT (({DIGIT}+(\.{DIGIT}*)?)|({DIGIT}*\.{DIGIT}+))([eE][+-]?{DIGIT}+)?
SPECIAL [;,()+\-*/^={}]
BLANK [ \t]
%%
{FLOAT} {
strncpy(yylval->val_str,yytext,MAXIDLENGTH);
RETTOK(FLOAT);
}
{SPECIAL} {RET(*yytext);}
return {RETTOK(RETURN);}
{ALPHA}({ALPHA}|{DIGIT})* {
strncpy(yylval->val_str,yytext,MAXIDLENGTH);
RETTOK(ID);
}
<*>\n {yylloc->last_column = 0;}
<*>{BLANK}
<*>. {yylval->val_char = yytext[0]; RETTOK(ERR);}
%%
void MathInterpreter::MathParserState::initScanner()
{
yylex_init(&scanner);
yyset_extra(this, scanner);
}
void MathInterpreter::MathParserState::destroyScanner()
{
yylex_destroy(scanner);
}

145
lib/MathParser.ypp Normal file
View File

@ -0,0 +1,145 @@
/*
* MathParser.ypp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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 <iostream>
#include <sstream>
#include <cstring>
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/MathInterpreter.hpp>
using namespace std;
using namespace Latan;
%}
%pure-parser
%name-prefix="_math_"
%locations
%defines
%error-verbose
%parse-param { Latan::MathInterpreter::MathParserState *state }
%initial-action {yylloc.last_column = 0;}
%lex-param { void* scanner }
%union
{
double val_double;
char val_char;
char val_str[MAXIDLENGTH];
Latan::ExprNode *val_node;
}
%token <val_char> ERR
%token <val_str> FLOAT
%token <val_str> ID
%token RETURN
%left '='
%left '+' '-'
%left '*' '/'
%nonassoc UMINUS
%left '^'
%type <val_node> stmt stmt_list expr func_args
%{
int _math_lex(YYSTYPE *lvalp, YYLTYPE *llocp, void *scanner);
void _math_error(YYLTYPE *locp, MathInterpreter::MathParserState *state,
const char *err)
{
stringstream buf;
buf << *(state->streamName) << ":" << locp->first_line << ":"\
<< locp->first_column << ": " << err;
LATAN_ERROR(Parsing, buf.str());
}
void _math_warning(YYLTYPE *locp, MathInterpreter::MathParserState *state,
const char *err)
{
stringstream buf;
buf << *(state->streamName) << ":" << locp->first_line << ":"\
<< locp->first_column << ": " << err;
LATAN_WARNING(buf.str());
}
#define scanner state->scanner
%}
%%
program:
/* empty string */
| stmt_list {(*(state->data)).reset($1);}
;
stmt:
';'
{$$ = new SemicolonNode(";");}
| expr ';'
{$$ = nullptr; _math_warning(&yylloc, state, "useless statement removed");}
| ID '=' expr ';'
{$$ = new AssignNode("="); $$->pushArg(new VarNode($ID)); $$->pushArg($3);}
| RETURN expr ';'
{$$ = new ReturnNode("return"); $$->pushArg($2);}
| '{' stmt_list '}'
{$$ = $2;}
;
stmt_list:
stmt
{$$ = $1;}
| stmt_list stmt
{$$ = new SemicolonNode(";"); $$->pushArg($1); $$->pushArg($2);}
;
expr:
FLOAT
{$$ = new CstNode($FLOAT);}
| ID
{$$ = new VarNode($ID);}
| '-' expr %prec UMINUS
{$$ = new MathOpNode("-"); $$->pushArg($2);}
| expr '+' expr
{$$ = new MathOpNode("+"); $$->pushArg($1); $$->pushArg($3);}
| expr '-' expr
{$$ = new MathOpNode("-"); $$->pushArg($1); $$->pushArg($3);}
| expr '*' expr
{$$ = new MathOpNode("*"); $$->pushArg($1); $$->pushArg($3);}
| expr '/' expr
{$$ = new MathOpNode("/"); $$->pushArg($1); $$->pushArg($3);}
| expr '^' expr
{$$ = new MathOpNode("^"); $$->pushArg($1); $$->pushArg($3);}
| '(' expr ')'
{$$ = $2;}
| ID '(' func_args ')'
{$$ = $3; $$->setName($ID);}
;
func_args:
/* empty string */
{$$ = new FuncNode("");}
| expr
{$$ = new FuncNode(""); $$->pushArg($1);}
| func_args ',' expr
{$$ = $1; $$->pushArg($3);}
;

58
lib/Minimizer.cpp Normal file
View File

@ -0,0 +1,58 @@
/*
* Minimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Minimizer.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* Minimizer implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
Minimizer::Minimizer(Verbosity verbosity)
: verbosity_(verbosity)
{}
// access //////////////////////////////////////////////////////////////////////
Index Minimizer::getDim(void) const
{
return x_.size();
}
DVec & Minimizer::getState(void)
{
return x_;
}
Minimizer::Verbosity Minimizer::getVerbosity(void) const
{
return verbosity_;
}
void Minimizer::setInit(const DVec &x0)
{
x_ = x0;
}
void Minimizer::setVerbosity(const Verbosity verbosity)
{
verbosity_ = verbosity;
}

64
lib/Minimizer.hpp Normal file
View File

@ -0,0 +1,64 @@
/*
* Minimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Minimizer_hpp_
#define Latan_Minimizer_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Mat.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* Abstract minimizer class *
******************************************************************************/
class Minimizer
{
public:
enum class Verbosity
{
Silent = 0,
Normal = 1,
Debug = 2
};
public:
// constructor
explicit Minimizer(Verbosity verbosity = Verbosity::Silent);
// destructor
virtual ~Minimizer(void) = default;
// access
Index getDim(void) const;
Verbosity getVerbosity(void) const;
void setInit(const DVec &x0);
void setVerbosity(const Verbosity verbosity);
// minimization
virtual const DVec & operator()(const DoubleFunction &f) = 0;
protected:
// access
DVec & getState(void);
private:
DVec x_;
Verbosity verbosity_;
};
END_NAMESPACE
#endif // Latan_Minimizer_hpp_

129
lib/MinuitMinimizer.cpp Normal file
View File

@ -0,0 +1,129 @@
/*
* MinuitMinimizer.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/MinuitMinimizer.hpp>
#include <LatAnalyze/includes.hpp>
#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>
using namespace std;
using namespace ROOT;
using namespace Minuit2;
using namespace Latan;
/******************************************************************************
* 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.;
}
// minimization ////////////////////////////////////////////////////////////////
const DVec & MinuitMinimizer::operator()(const DoubleFunction &f)
{
DVec &x = getState();
Verbosity verbosity = getVerbosity();
if (f.getNArg() != x.size())
{
x.conservativeResize(f.getNArg());
}
// set parameters
MinuitFunction minuitF(f);
MnUserParameters parameters;
for (Index i = 0; i < x.size(); ++i)
{
parameters.Add("x_" + strFrom(i), x(i), 0.1*fabs(x(i)));
}
// pre-minimization
MnMigrad migrad1(minuitF, parameters, 1);
FunctionMinimum min = migrad1();
if (verbosity >= Verbosity::Debug)
{
cout << "-- MINUIT pre-minimization -----------------------------";
cout << min;
cout << "--------------------------------------------------------";
cout << endl;
}
for (unsigned int i = 0; i < x.size(); ++i)
{
parameters.SetValue(i, min.UserParameters().Value(i));
parameters.SetError(i, min.UserParameters().Error(i));
}
// minimization and output
MnMigrad migrad2(minuitF, parameters, 2);
min = migrad2();
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, min.UserParameters(), 2);
cout << "-- MINUIT scan -----------------------------------------";
cout << endl;
for (unsigned int i = 0; i < x.size(); i++)
{
cout << "Parameter x_" << i << endl;
scanRes = scanner.Scan(i);
plot(scanRes);
}
cout << "--------------------------------------------------------";
cout << endl;
}
return x;
}

61
lib/MinuitMinimizer.hpp Normal file
View File

@ -0,0 +1,61 @@
/*
* MinuitMinimizer.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_MinuitMinimizer_hpp_
#define Latan_MinuitMinimizer_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <Minuit2/FCNBase.h>
BEGIN_NAMESPACE
/******************************************************************************
* Minuit minimizer *
******************************************************************************/
class MinuitMinimizer: public Minimizer
{
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_;
};
public:
// constructors
using Minimizer::Minimizer;
// destructor
virtual ~MinuitMinimizer(void) = default;
// minimization
virtual const DVec & operator()(const DoubleFunction &f);
};
END_NAMESPACE
#endif // Latan_MinuitMinimizer_hpp_

103
lib/Model.cpp Normal file
View File

@ -0,0 +1,103 @@
/*
* Model.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Model.hpp>
#include <LatAnalyze/includes.hpp>
#include <functional>
using namespace std;
using namespace std::placeholders;
using namespace Latan;
/******************************************************************************
* Model implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
DoubleModel::DoubleModel(const Index nArg, const Index nPar, const vecFunc &f)
: size_(new ModelSize)
{
setFunction(nArg, nPar, f);
}
// access //////////////////////////////////////////////////////////////////////
Index DoubleModel::getNArg(void) const
{
return size_->nArg;
}
Index DoubleModel::getNPar(void) const
{
return size_->nPar;
}
void DoubleModel::setFunction(const Index nArg, const Index nPar,
const vecFunc &f)
{
size_->nArg = nArg;
size_->nPar = nPar;
f_ = f;
}
// error checking //////////////////////////////////////////////////////////////
void DoubleModel::checkSize(const Index nArg, const Index nPar) const
{
if (nArg != getNArg())
{
LATAN_ERROR(Size, "model argument vector has a wrong size (expected "
+ strFrom(getNArg()) + ", got " + strFrom(nArg)
+ ")");
}
if (nPar != getNPar())
{
LATAN_ERROR(Size, "model parameter vector has a wrong size (expected "
+ strFrom(getNPar()) + ", got " + strFrom(nPar)
+ ")");
}
}
// function call ///////////////////////////////////////////////////////////////
double DoubleModel::operator()(const DVec &arg, const DVec &par) const
{
checkSize(arg.size(), par.size());
return (*this)(arg.data(), par.data());
}
double DoubleModel::operator()(const vector<double> &arg,
const vector<double> &par) const
{
checkSize(static_cast<Index>(arg.size()), static_cast<Index>(par.size()));
return (*this)(arg.data(), par.data());
}
double DoubleModel::operator()(const double *data, const double *par) const
{
return (*this)(data, par);
}
// model bind //////////////////////////////////////////////////////////////////
DoubleFunction DoubleModel::getBind(const DVec &par) const
{
auto modelWithVec = [this](const double *_arg, const DVec &_par)
{return (*this)(_arg, _par.data());};
auto modelBind = bind(modelWithVec, _1, par);
return DoubleFunction(getNArg(), modelBind);
}

67
lib/Model.hpp Normal file
View File

@ -0,0 +1,67 @@
/*
* Model.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Model_hpp_
#define Latan_Model_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/Mat.hpp>
#include <vector>
BEGIN_NAMESPACE
/******************************************************************************
* Double model class *
******************************************************************************/
class DoubleModel
{
private:
typedef std::function<double(const double *, const double *)> vecFunc;
struct ModelSize{Index nArg, nPar;};
public:
// constructor
DoubleModel(const Index nArg = 0, const Index nPar = 0,
const vecFunc &f = nullFunction_);
// destructor
virtual ~DoubleModel(void) = default;
// access
virtual Index getNArg(void) const;
virtual Index getNPar(void) const;
void setFunction(const Index nArg = 0, const Index nPar = 0,
const vecFunc &f = nullFunction_);
// function call
double operator()(const DVec &data, const DVec &par) const;
double operator()(const std::vector<double> &data,
const std::vector<double> &par) const;
virtual double operator()(const double *data, const double *par) const;
// model bind
DoubleFunction getBind(const DVec &par) const;
private:
// error checking
void checkSize(const Index nArg, const Index nPar) const;
private:
std::shared_ptr<ModelSize> size_;
vecFunc f_;
static const vecFunc nullFunction_;
};
END_NAMESPACE
#endif // Latan_Model_hpp_

59
lib/ParserState.hpp Normal file
View File

@ -0,0 +1,59 @@
/*
* ParserState.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_ParserState_hpp_
#define Latan_ParserState_hpp_
#include <LatAnalyze/Global.hpp>
#include <iostream>
BEGIN_NAMESPACE
template <typename DataObj>
class ParserState
{
public:
// constructor
ParserState(std::istream *streamPt, std::string *namePt, DataObj *dataPt);
// destructor
virtual ~ParserState(void) = default;
private:
// scanner allocation/deallocation
virtual void initScanner(void) = 0;
virtual void destroyScanner(void) = 0;
public:
DataObj *data;
void *scanner;
std::istream *stream;
std::string *streamName;
};
template <typename DataObj>
ParserState<DataObj>::ParserState(std::istream *streamPt, std::string *namePt,
DataObj *dataPt)
: data(dataPt)
, scanner(nullptr)
, stream(streamPt)
, streamName(namePt)
{}
END_NAMESPACE
#endif // Latan_ParserState_hpp_

380
lib/Plot.cpp Normal file
View File

@ -0,0 +1,380 @@
/*
* Plot.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/Plot.hpp>
#include <LatAnalyze/includes.hpp>
#include <LatAnalyze/Mat.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* Plot objects *
******************************************************************************/
// PlotObject access ///////////////////////////////////////////////////////////
const string & PlotObject::getCommand(void) const
{
return command_;
}
void PlotObject::setCommand(const string &command)
{
command_ = command;
}
string PlotObject::popTmpFile(void)
{
string res = tmpFileName_.top();
tmpFileName_.pop();
return res;
}
void PlotObject::pushTmpFile(const std::string &fileName)
{
tmpFileName_.push(fileName);
}
// PlotObject dump a matrix to a temporary file ////////////////////////////////
string PlotObject::dumpToTmpFile(const DMat &m)
{
char tmpFileName[MAX_PATH_LENGTH];
int fd;
FILE *tmpFile;
for (Index j = 0; j < m.cols(); ++j)
{
}
strcpy(tmpFileName, "./latan_plot_tmp.XXXXXX.dat");
fd = mkstemps(tmpFileName, 4);
if (fd == -1)
{
LATAN_ERROR(System, "impossible to create a temporary file from template "
+ strFrom(tmpFileName));
}
tmpFile = fdopen(fd, "w");
for (Index i = 0; i < m.rows(); ++i)
{
for (Index j = 0; j < m.cols(); ++j)
{
fprintf(tmpFile, "%e ", m(i, j));
}
fprintf(tmpFile, "\n");
}
fclose(tmpFile);
return string(tmpFileName);
}
// PlotObject test /////////////////////////////////////////////////////////////
bool PlotObject::gotTmpFile(void) const
{
return !tmpFileName_.empty();
}
// PlotCommand constructor /////////////////////////////////////////////////////
PlotCommand::PlotCommand(const string &command)
{
setCommand(command);
}
// PlotData constructor ////////////////////////////////////////////////////////
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(i, i).diagonal().array().sqrt();
usingCmd = (data.isXExact(i)) ? "u 1:3:4 w yerr" : "u 1:2:3:4 w xyerr";
tmpFileName = dumpToTmpFile(d);
pushTmpFile(tmpFileName);
setCommand("'" + tmpFileName + "' " + usingCmd);
}
// PlotFunction constructor ////////////////////////////////////////////////////
PlotFunction::PlotFunction(const DoubleFunction &function, const double xMin,
const double xMax, const unsigned int nSample)
{
DMat d(nSample, 2);
string tmpFileName;
double dx = (xMax - xMin)/static_cast<double>(nSample - 1);
for (Index i = 0; i < nSample; ++i)
{
d(i, 0) = xMin + i*dx;
d(i, 1) = function(d(i, 0));
}
tmpFileName = dumpToTmpFile(d);
pushTmpFile(tmpFileName);
setCommand("'" + tmpFileName + "' u 1:2 w lines");
}
/******************************************************************************
* Plot modifiers *
******************************************************************************/
// Color constructor ///////////////////////////////////////////////////////////
Color::Color(const string &color)
: color_(color)
{}
// Color modifier //////////////////////////////////////////////////////////////
void Color::operator()(PlotOptions &option) const
{
option.lineColor = color_;
}
// LogScale constructor ////////////////////////////////////////////////////////
LogScale::LogScale(const Axis axis)
: axis_(axis)
{}
// Logscale modifier ///////////////////////////////////////////////////////////
void LogScale::operator()(PlotOptions &option) const
{
option.scaleMode[static_cast<int>(axis_)] |= Plot::Scale::log;
}
// PlotRange constructor ////////////////////////////////////////////////////////
PlotRange::PlotRange(const Axis axis, const double min, const double max)
: axis_(axis)
, min_(min)
, max_(max)
{}
// PlotRange modifier ///////////////////////////////////////////////////////////
void PlotRange::operator()(PlotOptions &option) const
{
int a = static_cast<int>(axis_);
option.scaleMode[a] |= Plot::Scale::manual;
option.scale[a].min = min_;
option.scale[a].max = max_;
}
/******************************************************************************
* Plot implementation *
******************************************************************************/
// destructor //////////////////////////////////////////////////////////////////
Plot::~Plot(void)
{
while (!tmpFileName_.empty())
{
if (remove(tmpFileName_.top().c_str()))
{
LATAN_ERROR(System, "impossible to remove temporary file '" +
tmpFileName_.top() + "'");
}
tmpFileName_.pop();
}
}
// plot objects ////////////////////////////////////////////////////////////////
Plot & Plot::operator<<(PlotObject &&command)
{
string commandStr;
while (command.gotTmpFile())
{
tmpFileName_.push(command.popTmpFile());
}
commandStr = command.getCommand();
if (!options_.lineColor.empty())
{
commandStr += " lc " + options_.lineColor;
options_.lineColor = "";
}
if (options_.title.empty())
{
commandStr += " notitle";
}
else
{
commandStr += " t '" + options_.title + "'";
options_.title = "";
}
plotCommand_.push_back(commandStr);
return *this;
}
Plot & Plot::operator<<(PlotModifier &&modifier)
{
modifier(options_);
return *this;
}
// find gnuplot ////////////////////////////////////////////////////////////////
void 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)
{
for (i=0;path[i];)
{
for (j=i;(path[j])&&(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++ ;
}
}
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 */
lg = (int)(strlen(buf) - 1);
while (buf[lg]!='/')
{
buf[lg] = 0;
lg--;
}
buf[lg] = 0;
gnuplotPath_ = buf;
}
}
// plot parsing and output /////////////////////////////////////////////////////
void Plot::display(void)
{
std::string command;
FILE *gnuplotPipe;
if (!getenv("DISPLAY"))
{
LATAN_ERROR(System, "cannot find DISPLAY variable: is it set ?");
}
getProgramPath();
command = gnuplotPath_ + "/" + gnuplotBin_ + " " + gnuplotArgs_;
gnuplotPipe = popen(command.c_str(), "w");
if (!gnuplotPipe)
{
LATAN_ERROR(System, "error starting gnuplot (command was '" + command
+ "')");
}
commandBuffer_.str("");
commandBuffer_ << *this;
fprintf(gnuplotPipe, "%s", commandBuffer_.str().c_str());
if (pclose(gnuplotPipe) == -1)
{
LATAN_ERROR(System, "problem closing communication to gnuplot");
}
}
ostream & Latan::operator<<(ostream &out, const Plot &plot)
{
std::string begin, end;
int x = static_cast<int>(Axis::x), y = static_cast<int>(Axis::y);
if (!plot.options_.terminal.empty())
{
out << "set term " << plot.options_.terminal << endl;
}
if (!plot.options_.output.empty())
{
out << "set output '" << plot.options_.terminal << "'" << endl;
}
if (plot.options_.scaleMode[x] & Plot::Scale::manual)
{
out << "xMin = " << plot.options_.scale[x].min << endl;
out << "xMax = " << plot.options_.scale[x].max << endl;
}
if (plot.options_.scaleMode[y] & Plot::Scale::manual)
{
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;
}
if (plot.options_.scaleMode[x] & Plot::Scale::manual)
{
out << "set xrange [xMin:xMax]" << endl;
}
if (plot.options_.scaleMode[y] & Plot::Scale::manual)
{
out << "set yrange [yMin:yMax]" << endl;
}
if (plot.options_.scaleMode[x] & Plot::Scale::log)
{
out << "set log x" << endl;
}
if (plot.options_.scaleMode[y] & Plot::Scale::log)
{
out << "set log y" << endl;
}
if (!plot.options_.label[x].empty())
{
out << "set xlabel '" << plot.options_.label[x] << "'" << endl;
}
if (!plot.options_.label[y].empty())
{
out << "set ylabel '" << plot.options_.label[y] << "'" << endl;
}
for (unsigned int i = 0; i < plot.headCommand_.size(); ++i)
{
out << plot.headCommand_[i] << endl;
}
for (unsigned int i = 0; i < plot.plotCommand_.size(); ++i)
{
begin = (i == 0) ? "plot " : " ";
end = (i == plot.plotCommand_.size() - 1) ? "" : ",\\";
out << begin << plot.plotCommand_[i] << end << endl;
}
return out;
}

213
lib/Plot.hpp Normal file
View File

@ -0,0 +1,213 @@
/*
* Plot.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_Plot_hpp_
#define Latan_Plot_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/Mat.hpp>
#include <LatAnalyze/XYStatData.hpp>
#include <sstream>
#include <stack>
#include <vector>
// gnuplot default parameters
#ifndef GNUPLOT_BIN
#define GNUPLOT_BIN "gnuplot"
#endif
#ifndef GNUPLOT_ARGS
#define GNUPLOT_ARGS "-p"
#endif
BEGIN_NAMESPACE
/******************************************************************************
* Plot objects *
******************************************************************************/
class PlotObject
{
public:
// destructor
virtual ~PlotObject(void) = default;
// access
std::string popTmpFile(void);
const std::string & getCommand(void) const;
// test
bool gotTmpFile(void) const;
protected:
// access
void pushTmpFile(const std::string &fileName);
void setCommand(const std::string &command);
// dump a matrix to a temporary file
std::string dumpToTmpFile(const DMat &m);
private:
// plot command
std::string command_;
// stack of created temporary files
std::stack<std::string> tmpFileName_;
};
class PlotCommand: public PlotObject
{
public:
// constructor
explicit PlotCommand(const std::string &command);
// destructor
virtual ~PlotCommand(void) = default;
};
class PlotData: public PlotObject
{
public:
// constructor
PlotData(const XYStatData &data, const Index i = 0, const Index j = 0);
// destructor
virtual ~PlotData(void) = default;
};
class PlotFunction: public PlotObject
{
public:
// constructor
PlotFunction(const DoubleFunction &function, const double xMin,
const double xMax, const unsigned int nSample = 1000);
// destructor
virtual ~PlotFunction(void) = default;
};
/******************************************************************************
* Plot modifiers *
******************************************************************************/
enum class Axis {x = 0, y = 1};
struct Range
{
double min, max;
};
struct PlotOptions
{
std::string terminal {""};
std::string output {""};
std::string title {""};
unsigned int scaleMode[2] {0,0};
Range scale[2] {{0.0,0.0},{0.0,0.0}};
std::string label[2] {"",""};
std::string lineColor {""};
};
class PlotModifier
{
public:
// destructor
virtual ~PlotModifier(void) = default;
// modifier
virtual void operator()(PlotOptions &option) const = 0;
};
class Color: public PlotModifier
{
public:
// constructor
explicit Color(const std::string &color);
// destructor
virtual ~Color(void) = default;
// modifier
virtual void operator()(PlotOptions &option) const;
private:
const std::string color_;
};
class LogScale: public PlotModifier
{
public:
// constructor
explicit LogScale(const Axis axis);
// destructor
virtual ~LogScale(void) = default;
// modifier
virtual void operator()(PlotOptions &option) const;
private:
const Axis axis_;
};
class PlotRange: public PlotModifier
{
public:
// constructor
PlotRange(const Axis axis, const double min, const double max);
// destructor
virtual ~PlotRange(void) = default;
// modifier
virtual void operator()(PlotOptions &option) const;
private:
const Axis axis_;
const double min_, max_;
};
/******************************************************************************
* Plot class *
******************************************************************************/
class Plot
{
public:
class Scale
{
public:
enum
{
reset = 0,
manual = 1 << 0,
log = 1 << 1
};
};
public:
// constructor/destructor
Plot(void) = default;
virtual ~Plot(void);
// plot operations
Plot & operator<<(PlotObject &&command);
Plot & operator<<(PlotModifier &&modifier);
// plot parsing and output
void display(void);
friend std::ostream & operator<<(std::ostream &out, const Plot &plot);
private:
// find gnuplot
void getProgramPath(void);
private:
// gnuplot execution parameters
std::string gnuplotBin_ {GNUPLOT_BIN};
std::string gnuplotArgs_ {GNUPLOT_ARGS};
std::string gnuplotPath_ {""};
// string buffer for commands
std::ostringstream commandBuffer_;
// stack of created temporary files
std::stack<std::string> tmpFileName_;
// plot content
PlotOptions options_;
std::vector<std::string> headCommand_;
std::vector<std::string> plotCommand_;
};
std::ostream & operator<<(std::ostream &out, const Plot &plot);
END_NAMESPACE
#endif // Latan_Plot_hpp_

688
lib/RandGen.cpp Normal file
View File

@ -0,0 +1,688 @@
/*
* RandGen.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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>
#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;
}
double 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 Normal file
View File

@ -0,0 +1,103 @@
/*
* RandGen.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_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);
double discreteUniform(const unsigned int n);
double gaussian(const double mean = 0.0, const double sigma = 1.0);
private:
RanLxd generator_;
};
END_NAMESPACE
#endif // Latan_RandGen_hpp_

266
lib/StatArray.hpp Normal file
View File

@ -0,0 +1,266 @@
/*
* StatArray.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_StatArray_hpp_
#define Latan_StatArray_hpp_
#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)
BEGIN_NAMESPACE
/******************************************************************************
* Array class with statistics *
******************************************************************************/
template <typename T, Index os = 0>
class StatArray: public Array<T, dynamic, 1>
{
private:
typedef Array<T, dynamic, 1> Base;
public:
// constructors
StatArray(void);
explicit StatArray(const Index size);
EIGEN_EXPR_CTOR(StatArray, unique_arg(StatArray<T, os>), Base, ArrayBase)
// destructor
virtual ~StatArray(void) = default;
// access
Index size(void) const;
void resize(const Index size);
// operators
T & operator[](const Index s);
const T & operator[](const Index s) const;
// statistics
void bin(Index binSize);
T mean(const Index pos = 0, const Index n = -1) const;
T mean(void) const;
T covariance(const StatArray<T, os> &array, const Index pos = 0,
const Index n = -1) const;
T covarianceMatrix(const StatArray<T, os> &array, const Index pos = 0,
const Index n = -1) const;
T variance(const Index pos = 0, const Index n = -1) const;
T varianceMatrix(const Index pos = 0, const Index n = -1) const;
public:
static const Index offset = os;
};
// reduction operations
namespace ReducOp
{
// general templates
template <typename T>
inline T prod(const T &a, const T &b);
template <typename T>
inline T tensProd(const T &v1, const T &v2);
template <typename T>
inline T sum(const T &a, const T &b);
// matrix specializations
template <>
inline DMat prod(const DMat &a, const DMat &b);
template <>
inline DMat tensProd(const DMat &v1, const DMat &v2);
}
// Sample types
#define SAMPLE_OFFSET 1
const int central = -SAMPLE_OFFSET;
template <typename T>
using Sample = StatArray<T, SAMPLE_OFFSET>;
typedef Sample<double> DSample;
/******************************************************************************
* StatArray class template implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
template <typename T, Index os>
StatArray<T, os>::StatArray(void)
: Base(static_cast<typename Base::Index>(os))
{}
template <typename T, Index os>
StatArray<T, os>::StatArray(const Index size)
: Base(static_cast<typename Base::Index>(size + os))
{}
// access //////////////////////////////////////////////////////////////////////
template <typename T, Index os>
Index StatArray<T, os>::size(void) const
{
return Base::size() - os;
}
template <typename T, Index os>
void StatArray<T, os>::resize(const Index size)
{
Base::resize(size + os);
}
// operators ///////////////////////////////////////////////////////////////////
template <typename T, Index os>
T & StatArray<T, os>::operator[](const Index s)
{
return Base::operator[](s + os);
}
template <typename T, Index os>
const T & StatArray<T, os>::operator[](const Index s) const
{
return Base::operator[](s + os);
}
// statistics //////////////////////////////////////////////////////////////////
template <typename T, Index os>
void StatArray<T, os>::bin(Index binSize)
{
Index q = size()/binSize, r = size()%binSize;
for (Index i = 0; i < q; ++i)
{
(*this)[i] = mean(i*binSize, binSize);
}
if (r != 0)
{
(*this)[q] = mean(q*binSize, r);
this->conservativeResize(os + q + 1);
}
else
{
this->conservativeResize(os + q);
}
}
template <typename T, Index os>
T StatArray<T, os>::mean(const Index pos, const Index n) const
{
T result;
const Index m = (n >= 0) ? n : size();
if (m)
{
result = this->segment(pos+os, m).redux(&ReducOp::sum<T>);
}
return result/static_cast<double>(n);
}
template <typename T, Index os>
T StatArray<T, os>::covariance(const StatArray<T, os> &array, const Index pos,
const Index n) const
{
T s1, s2, prs, res;
const Index m = (n >= 0) ? n : size();
if (m)
{
auto arraySeg = array.segment(pos+os, m);
auto thisSeg = this->segment(pos+os, m);
s1 = thisSeg.redux(&ReducOp::sum<T>);
s2 = arraySeg.redux(&ReducOp::sum<T>);
prs = thisSeg.binaryExpr(arraySeg, &ReducOp::prod<T>)
.redux(&ReducOp::sum<T>);
res = prs - ReducOp::prod(s1, s2)/static_cast<double>(m);
}
return res/static_cast<double>(m - 1);
}
template <typename T, Index os>
T StatArray<T, os>::covarianceMatrix(const StatArray<T, os> &array,
const Index pos, const Index n) const
{
T s1, s2, prs, res;
const Index m = (n >= 0) ? n : size();
if (m)
{
auto arraySeg = array.segment(pos+os, m);
auto thisSeg = this->segment(pos+os, m);
s1 = thisSeg.redux(&ReducOp::sum<T>);
s2 = arraySeg.redux(&ReducOp::sum<T>);
prs = thisSeg.binaryExpr(arraySeg, &ReducOp::tensProd<T>)
.redux(&ReducOp::sum<T>);
res = prs - ReducOp::tensProd(s1, s2)/static_cast<double>(m);
}
return res/static_cast<double>(m - 1);
}
template <typename T, Index os>
T StatArray<T, os>::variance(const Index pos, const Index n) const
{
return covariance(*this, pos, n);
}
template <typename T, Index os>
T StatArray<T, os>::varianceMatrix(const Index pos, const Index n) const
{
return covarianceMatrix(*this, pos, n);
}
// reduction operations ////////////////////////////////////////////////////////
template <typename T>
inline T ReducOp::sum(const T &a, const T &b)
{
return a + b;
}
template <typename T>
inline T ReducOp::prod(const T &a, const T &b)
{
return a*b;
}
template <typename T>
inline T ReducOp::tensProd(const T &v1 __unused, const T &v2 __unused)
{
LATAN_ERROR(Implementation,
"tensorial product not implemented for this type");
}
template <>
inline DMat ReducOp::prod(const DMat &a, const DMat &b)
{
return a.cwiseProduct(b);
}
template <>
inline DMat ReducOp::tensProd(const DMat &v1, const DMat &v2)
{
if ((v1.cols() != 1)||(v2.cols() != 1))
{
LATAN_ERROR(Size,
"tensorial product is only valid with column vectors");
}
return v1*v2.transpose();
}
END_NAMESPACE
#endif // Latan_StatArray_hpp_

294
lib/XYSampleData.cpp Normal file
View File

@ -0,0 +1,294 @@
/*
* XYSampleData.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/XYSampleData.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* SampleFitResult implementation *
******************************************************************************/
double SampleFitResult::getChi2(const Index s) const
{
return chi2_[s];
}
const DSample & SampleFitResult::getChi2(const PlaceHolder ph __unused) const
{
return chi2_;
}
double SampleFitResult::getChi2PerDof(const Index s) const
{
return chi2_[s]/static_cast<double>(nDof_);
}
DSample SampleFitResult::getChi2PerDof(const PlaceHolder ph __unused) const
{
return chi2_/static_cast<double>(nDof_);
}
const DoubleFunction & SampleFitResult::getModel(const Index s,
const Index j) const
{
return model_[static_cast<unsigned int>(j)][s];
}
const DoubleFunctionSample & SampleFitResult::getModel(
const PlaceHolder ph __unused,
const Index j) const
{
return model_[static_cast<unsigned int>(j)];
}
/******************************************************************************
* XYSampleData implementation *
******************************************************************************/
// constructors ////////////////////////////////////////////////////////////////
XYSampleData::XYSampleData(const Index nData, const Index xDim,
const Index yDim, const Index nSample)
{
resize(nData, xDim, yDim, nSample);
}
// access //////////////////////////////////////////////////////////////////////
const XYStatData & XYSampleData::getData(const Index s)
{
setDataToSample(s);
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 __unused,
const PlaceHolder ph2 __unused)
{
isCovarianceInit_ = false;
return x_.block(0, 0, getNData(), getXDim());
}
XYSampleData::ConstSampleBlock XYSampleData::x(const PlaceHolder ph1 __unused,
const PlaceHolder ph2 __unused)
const
{
return x_.block(0, 0, getNData(), getXDim());
}
XYSampleData::SampleBlock XYSampleData::x(const Index i,
const PlaceHolder ph2 __unused)
{
isCovarianceInit_ = false;
return x_.block(0, i, getNData(), 1);
}
XYSampleData::ConstSampleBlock XYSampleData::x(const Index i,
const PlaceHolder ph2 __unused)
const
{
return x_.block(0, i, getNData(), 1);
}
XYSampleData::SampleBlock XYSampleData::x(const PlaceHolder ph1 __unused,
const Index k)
{
isCovarianceInit_ = false;
return x_.block(k, 0, 1, getXDim());
}
XYSampleData::ConstSampleBlock XYSampleData::x(const PlaceHolder ph1 __unused,
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 __unused,
const PlaceHolder ph2 __unused)
{
isCovarianceInit_ = false;
return y_.block(0, 0, getNData(), getXDim());
}
XYSampleData::ConstSampleBlock XYSampleData::y(const PlaceHolder ph1 __unused,
const PlaceHolder ph2 __unused)
const
{
return y_.block(0, 0, getNData(), getXDim());
}
XYSampleData::SampleBlock XYSampleData::y(const Index j,
const PlaceHolder ph2 __unused)
{
isCovarianceInit_ = false;
return y_.block(0, j, getNData(), 1);
}
XYSampleData::ConstSampleBlock XYSampleData::y(const Index j,
const PlaceHolder ph2 __unused)
const
{
return y_.block(0, j, getNData(), 1);
}
XYSampleData::SampleBlock XYSampleData::y(const PlaceHolder ph1 __unused,
const Index k)
{
isCovarianceInit_ = false;
return y_.block(k, 0, 1, getXDim());
}
XYSampleData::ConstSampleBlock XYSampleData::y(const PlaceHolder ph1 __unused,
const Index k) const
{
return y_.block(k, 0, 1, getXDim());
}
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(
const std::vector<const DoubleModel *> &modelVector,
Minimizer &minimizer, const DVec &init,
const FitVerbosity verbosity)
{
const Index nSample = x_.size();
FitResult sampleResult;
SampleFitResult result;
bool initChi2;
// copy interface to sample data
data_.setFitInterface(*this);
// sample loop
result.resize(nSample);
result.chi2_.resize(nSample);
FOR_STAT_ARRAY(x_, s)
{
// set data
setDataToSample(s);
// initialize chi^2 only once
initChi2 = (s == central);
// fit
sampleResult = data_.fit(modelVector, minimizer, init, initChi2,
verbosity);
// store result
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)
{
result.model_[j].resize(nSample);
result.model_[j][s] = sampleResult.getModel(j);
}
}
return result;
}
SampleFitResult XYSampleData::fit(const DoubleModel &model,
Minimizer &minimizer,
const DVec &init,
const FitVerbosity verbosity)
{
vector<const DoubleModel *> modelVector(1);
modelVector[0] = &model;
return fit(modelVector, minimizer, init, verbosity);
}
void XYSampleData::setDataToSample(const Index s)
{
// 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;
}
// set data
data_.x() = x_[s];
data_.y() = y_[s];
}

113
lib/XYSampleData.hpp Normal file
View File

@ -0,0 +1,113 @@
/*
* XYSampleData.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_XYSampleData_hpp_
#define Latan_XYSampleData_hpp_
#include <LatAnalyze/Global.hpp>
#include <LatAnalyze/FitInterface.hpp>
#include <LatAnalyze/Function.hpp>
#include <LatAnalyze/MatSample.hpp>
#include <LatAnalyze/Minimizer.hpp>
#include <LatAnalyze/XYStatData.hpp>
BEGIN_NAMESPACE
/******************************************************************************
* object for fit result *
******************************************************************************/
class SampleFitResult: public DMatSample
{
friend class XYSampleData;
public:
// constructors
using DMatSample::DMatSample;
SampleFitResult(void) = default;
// destructor
virtual ~SampleFitResult(void) = default;
// access
double getChi2(const Index s = central) const;
const DSample & getChi2(const PlaceHolder ph) const;
double getChi2PerDof(const Index s = central) const;
DSample getChi2PerDof(const PlaceHolder ph) const;
const DoubleFunction & getModel(const Index s = central,
const Index j = 0) const;
const DoubleFunctionSample & getModel(const PlaceHolder ph,
const Index j = 0) const;
private:
DSample chi2_;
Index nDof_{0};
std::vector<DoubleFunctionSample> model_;
};
/******************************************************************************
* 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);
// 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 i, const PlaceHolder ph2 = _);
ConstSampleBlock y(const Index i, 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 i, const Index k);
ConstSampleBlock y(const Index i, const Index k) const;
// fit
SampleFitResult fit(const std::vector<const DoubleModel *> &modelVector,
Minimizer &minimizer, const DVec &init,
const FitVerbosity verbosity = FitVerbosity::Silent);
SampleFitResult fit(const DoubleModel &model, Minimizer &minimizer,
const DVec &init,
const FitVerbosity verbosity = FitVerbosity::Silent);
private:
void setDataToSample(const Index s);
private:
bool isCovarianceInit_;
DMatSample x_, y_;
XYStatData data_;
};
END_NAMESPACE
#endif // Latan_XYSampleData_hpp_

264
lib/XYStatData.cpp Normal file
View File

@ -0,0 +1,264 @@
/*
* XYStatData.cpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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/XYStatData.hpp>
#include <LatAnalyze/includes.hpp>
using namespace std;
using namespace Latan;
/******************************************************************************
* FitResult implementation *
******************************************************************************/
// access //////////////////////////////////////////////////////////////////////
double FitResult::getChi2(void) const
{
return chi2_;
}
Index FitResult::getNDof(void) const
{
return nDof_;
}
double FitResult::getChi2PerDof(void) const
{
return chi2_/static_cast<double>(nDof_);
}
const DoubleFunction & FitResult::getModel(const Index j) const
{
return model_[static_cast<unsigned int>(j)];
}
/******************************************************************************
* XYStatData implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
XYStatData::XYStatData(void)
: chi2_(*this)
{}
XYStatData::XYStatData(const Index nData, const Index xDim, const Index yDim)
: XYStatData()
{
resize(nData, xDim, yDim);
}
// access //////////////////////////////////////////////////////////////////////
void XYStatData::resize(const Index nData, const Index xDim, const Index yDim)
{
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)
{
var_[xx](i1, i2).resize(nData, nData);
}
FOR_MAT(var_[yy], j1, j2)
{
var_[yy](j1, j2).resize(nData, nData);
}
FOR_MAT(var_[yx], j, i)
{
var_[yx](j, i).resize(nData, nData);
}
}
Block<DMatBase> XYStatData::x(const PlaceHolder ph1 __unused,
const PlaceHolder ph2 __unused)
{
return x_.block(0, 0, getNData(), getXDim());
}
ConstBlock<DMatBase> XYStatData::x(const PlaceHolder ph1 __unused,
const PlaceHolder ph2 __unused) const
{
return x_.block(0, 0, getNData(), getXDim());
}
Block<DMatBase> XYStatData::x(const Index i, const PlaceHolder ph2 __unused)
{
return x_.block(0, i, getNData(), 1);
}
ConstBlock<DMatBase> XYStatData::x(const Index i,
const PlaceHolder ph2 __unused) const
{
return x_.block(0, i, getNData(), 1);
}
Block<DMatBase> XYStatData::x(const PlaceHolder ph1 __unused, const Index k)
{
return x_.block(k, 0, 1, getXDim());
}
ConstBlock<DMatBase> XYStatData::x(const PlaceHolder ph1 __unused,
const Index k) const
{
return x_.block(k, 0, 1, getXDim());
}
double & XYStatData::x(const Index i, const Index k)
{
return x_(k, i);
}
const double & XYStatData::x(const Index i, const Index k) const
{
return x_(k, i);
}
Block<DMatBase> XYStatData::y(const PlaceHolder ph1 __unused,
const PlaceHolder ph2 __unused)
{
return y_.block(0, 0, getNData(), getYDim());
}
ConstBlock<DMatBase> XYStatData::y(const PlaceHolder ph1 __unused,
const PlaceHolder ph2 __unused) const
{
return y_.block(0, 0, getNData(), getYDim());
}
Block<DMatBase> XYStatData::y(const Index j, const PlaceHolder ph2 __unused)
{
return y_.block(0, j, getNData(), 1);
}
ConstBlock<DMatBase> XYStatData::y(const Index j,
const PlaceHolder ph2 __unused) const
{
return y_.block(0, j, getNData(), 1);
}
Block<DMatBase> XYStatData::y(const PlaceHolder ph1 __unused, const Index k)
{
return y_.block(k, 0, 1, getYDim());
}
ConstBlock<DMatBase> XYStatData::y(const PlaceHolder ph1 __unused,
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<DMatBase> XYStatData::xxVar(const Index i1, const Index i2)
{
return FULL_BLOCK(var_[xx](i1, i2));
}
ConstBlock<DMatBase> XYStatData::xxVar(const Index i1, const Index i2) const
{
return FULL_BLOCK(var_[xx](i1, i2));
}
Block<DMatBase> XYStatData::yyVar(const Index j1, const Index j2)
{
return FULL_BLOCK(var_[yy](j1, j2));
}
ConstBlock<DMatBase> XYStatData::yyVar(const Index j1, const Index j2) const
{
return FULL_BLOCK(var_[yy](j1, j2));
}
Block<DMatBase> XYStatData::yxVar(const Index j, const Index i)
{
return FULL_BLOCK(var_[yx](j, i));
}
ConstBlock<DMatBase> XYStatData::yxVar(const Index j, const Index i) const
{
return FULL_BLOCK(var_[yx](j, i));
}
// fit /////////////////////////////////////////////////////////////////////////
FitResult XYStatData::fit(const vector<const DoubleModel *> &modelVector,
Minimizer &minimizer, const DVec &init,
const bool reinitChi2,
const FitVerbosity verbosity)
{
// initialization
chi2_.setModel(modelVector);
if (reinitChi2)
{
chi2_.requestInit();
}
minimizer.setVerbosity(verbosity);
// 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)
if (!isXExact(i))
{
for (Index k = 0; k < getNData(); ++k)
if (isFitPoint(k))
{
fullInit(chi2_.getNPar() + nPoint*is + kf) = x(i, k);
kf++;
}
is++;
}
minimizer.setInit(fullInit);
// fit
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.model_[j] = modelVector[j]->getBind(result);
}
return result;
}
FitResult XYStatData::fit(const DoubleModel &model, Minimizer &minimizer,
const DVec &init, const bool reinitChi2,
const FitVerbosity verbosity)
{
vector<const DoubleModel *> modelVector(1);
modelVector[0] = &model;
return fit(modelVector, minimizer, init, reinitChi2, verbosity);
}

122
lib/XYStatData.hpp Normal file
View File

@ -0,0 +1,122 @@
/*
* XYData.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_XYData_hpp_
#define Latan_XYData_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_NAMESPACE
/******************************************************************************
* object for fit result *
******************************************************************************/
class FitResult: public DVec
{
friend class XYStatData;
public:
// constructors
FitResult(void) = default;
EIGEN_EXPR_CTOR(FitResult, FitResult, Base, MatrixBase)
// destructor
virtual ~FitResult(void) = default;
// access
double getChi2(void) const;
Index getNDof(void) const;
double getChi2PerDof(void) const;
const DoubleFunction & getModel(const Index j = 0) const;
private:
double chi2_{0.0};
Index nDof_{0};
std::vector<DoubleFunction> model_;
};
/******************************************************************************
* object 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);
// destructor
virtual ~XYStatData(void) = default;
// access
void resize(const Index nData, const Index xDim,
const Index yDim);
Block<DMatBase> x(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _);
ConstBlock<DMatBase> x(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _) const;
Block<DMatBase> x(const Index i, const PlaceHolder ph2 = _);
ConstBlock<DMatBase> x(const Index i, const PlaceHolder ph2 = _) const;
Block<DMatBase> x(const PlaceHolder ph1, const Index k);
ConstBlock<DMatBase> 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<DMatBase> y(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _);
ConstBlock<DMatBase> y(const PlaceHolder ph1 = _,
const PlaceHolder ph2 = _) const;
Block<DMatBase> y(const Index i, const PlaceHolder ph2 = _);
ConstBlock<DMatBase> y(const Index i, const PlaceHolder ph2 = _) const;
Block<DMatBase> y(const PlaceHolder ph1, const Index k);
ConstBlock<DMatBase> 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<DMatBase> xxVar(const Index i1, const Index i2);
ConstBlock<DMatBase> xxVar(const Index i1, const Index i2) const;
Block<DMatBase> yyVar(const Index j1, const Index j2);
ConstBlock<DMatBase> yyVar(const Index j1, const Index j2) const;
Block<DMatBase> yxVar(const Index j, const Index i);
ConstBlock<DMatBase> yxVar(const Index j, const Index i) const;
// fit
FitResult fit(const std::vector<const DoubleModel *> &modelVector,
Minimizer &minimizer, const DVec &init,
const bool reinitChi2 = true,
const FitVerbosity verbosity = FitVerbosity::Silent);
FitResult fit(const DoubleModel &model, Minimizer &minimizer,
const DVec &init, const bool reinitChi2 = true,
const FitVerbosity verbosity = FitVerbosity::Silent);
private:
DMat x_, y_;
Mat<DMat> var_[3];
IVec isXExact_, isFitPoint_;
Chi2Function chi2_;
};
END_NAMESPACE
#endif // Latan_XYData_hpp_

39
lib/includes.hpp Normal file
View File

@ -0,0 +1,39 @@
/*
* includes.hpp, part of LatAnalyze 3
*
* Copyright (C) 2013 - 2014 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_includes_hpp_
#define Latan_includes_hpp_
#include <fstream>
#include <iostream>
#include <iomanip>
#include <iostream>
#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_