Commit 9e578a9f authored by Tobias WEBER's avatar Tobias WEBER
Browse files

forked jl module from private repo

parent e954d5b4
......@@ -29,7 +29,7 @@
#include <boost/property_tree/ini_parser.hpp>
#include <boost/optional.hpp>
#include "string.h"
#include "str.h"
#include "traits.h"
......
......@@ -23,7 +23,7 @@
#include <algorithm>
#include <limits>
#include "math.h"
#include "mat.h"
#include "log.h"
#include "traits.h"
......@@ -217,6 +217,9 @@ bool fit(t_func&& func,
{
try
{
for(auto x : vecYErr)
std::cout << x << ", ";
std::cout << std::endl;
if(!vecX.size() || !vecY.size() || !vecYErr.size())
{
log_err("No data given to fitter.");
......
......@@ -20,10 +20,10 @@
#include <iostream>
#include <iomanip>
#include "string.h"
#include "str.h"
#include "log.h"
#include "file.h"
#include "math.h"
#include "mat.h"
#include "phys.h"
......
/**
* tlibs2
* julia module
* @author Tobias Weber <tobias.weber@tum.de>, <tweber@ill.fr>
* @date 2017 -- 2018
* @license GPLv3, see 'LICENSE' file
* @desc Forked on 7-Nov-2018 from the privately and TUM-PhD-developed "tlibs" project (https://github.com/t-weber/tlibs).
*
* g++-8 -std=c++17 -shared -fPIC -O2 -march=native -I. -I/usr/local/include/julia -I/usr/include/julia -o tl2jl.so log.cpp jl.cpp -lboost_system -lboost_iostreams -lMinuit2 -lgomp -ljulia
*/
#include "jl.h"
#include "log.h"
#include "instr.h"
#include "fit.h"
using t_real = double;
bool g_bDebug = 0;
extern "C" void load_tl2(int bDebug)
{
g_bDebug = (bDebug != 0);
tl2::log_debug.SetEnabled(g_bDebug);
tl2::log_debug("Loaded tl2 jl module.");
}
// ----------------------------------------------------------------------------
/**
* loads an instrument data file
*/
extern "C" jl_array_t* load_instr(const char* pcFile)
{
tl2::FileInstrBase<t_real>* pInstr = tl2::FileInstrBase<t_real>::LoadInstr(pcFile);
if(!pInstr)
{
jl_array_t *pArrNull = jl_alloc_array_1d(
jl_apply_array_type(reinterpret_cast<jl_value_t*>(jl_any_type), 1), 0);
tl2::log_err("In ", __func__, ": Cannot load ", pcFile, ".");
return pArrNull;
}
// [ column names, data, keys, values ]
jl_array_t *pArr = jl_alloc_array_1d(
jl_apply_array_type(reinterpret_cast<jl_value_t*>(jl_any_type), 1), 4);
jl_array_t** pArrDat = reinterpret_cast<jl_array_t**>(jl_array_data(pArr));
// data column names
pArrDat[0] = tl2::make_jl_str_arr(pInstr->GetColNames());
// data matrix
pArrDat[1] = tl2::make_jl_2darr(pInstr->GetData());
// scan property map
std::tie(pArrDat[2], pArrDat[3]) = tl2::make_jl_strmap_arr(pInstr->GetAllParams());
return pArr;
}
// ----------------------------------------------------------------------------
// Fitting
/**
* internal helper function to call fitter using variable args
*/
template<std::size_t iNumArgs>
static inline bool _invoke_fit(void *_pFkt,
const std::vector<tl2::t_real_min>& vecX, const std::vector<tl2::t_real_min>& vecY,
const std::vector<tl2::t_real_min>& vecYerr,
const std::vector<std::string>& vecParamNames,
std::vector<tl2::t_real_min>& vecVals, std::vector<tl2::t_real_min>& vecErrs,
const std::vector<bool>& vecFixed)
{
auto *pFkt = reinterpret_cast<tl2::t_fkt_vararg<t_real, iNumArgs>>(_pFkt);
return tl2::fit<iNumArgs>(pFkt,
vecX, vecY, vecYerr, vecParamNames, vecVals, vecErrs, &vecFixed, g_bDebug);
}
/**
* function fitting
*/
extern "C" int fit(void *_pFkt, std::size_t iNumParams,
const t_real *pX, const t_real *pY, const t_real *pYerr, std::size_t iArrLen,
jl_array_t *parrParamNames, jl_array_t *parrFixed,
t_real *pValues, t_real *pErrors)
{
std::vector<tl2::t_real_min> vecX, vecY, vecYerr;
vecX.reserve(iArrLen);
vecY.reserve(iArrLen);
vecYerr.reserve(iArrLen);
// copy arrays to vectors
for(std::size_t i=0; i<iArrLen; ++i)
{
vecX.push_back(pX[i]);
vecY.push_back(pY[i]);
vecYerr.push_back(pYerr[i]);
}
// parameter names
std::vector<std::string> vecParamNames =
tl2::from_jl_arr<std::vector, std::string>(parrParamNames, 1);
// fixed parameters
std::vector<std::string> vecFixedParams =
tl2::from_jl_arr<std::vector, std::string>(parrFixed);
std::vector<bool> vecFixed;
for(const std::string& strParam : vecParamNames)
{
bool bFixed = std::find(vecFixedParams.begin(), vecFixedParams.end(), strParam)
!= vecFixedParams.end();
vecFixed.push_back(bFixed);
}
// values & errors
std::vector<tl2::t_real_min> vecVals, vecErrs;
for(std::size_t iIdx=0; iIdx<iNumParams; ++iIdx)
{
vecVals.push_back(pValues[iIdx]);
vecErrs.push_back(pErrors[iIdx]);
}
// ------------------------------------------------------------------------
// fill up missing parameters and hints
if(vecParamNames.size() < iNumParams)
{
for(std::size_t iArg=vecParamNames.size(); iArg<iNumParams; ++iArg)
{
std::ostringstream ostrArg;
ostrArg << "arg_" << iArg;
vecParamNames.push_back(ostrArg.str());
}
}
while(vecVals.size() < iNumParams) vecVals.push_back(0);
while(vecErrs.size() < iNumParams) vecErrs.push_back(0);
while(vecFixed.size() < iNumParams) vecFixed.push_back(0);
// ------------------------------------------------------------------------
#define __CALL_FIT(NUM) _invoke_fit<NUM+1>(_pFkt, vecX, vecY, vecYerr, vecParamNames, vecVals, vecErrs, vecFixed)
bool bOk = 0;
// stating all needed specialisations of the fit template function
if(iNumParams == 1) bOk = __CALL_FIT(1);
else if(iNumParams == 2) bOk = __CALL_FIT(2);
else if(iNumParams == 3) bOk = __CALL_FIT(3);
else if(iNumParams == 4) bOk = __CALL_FIT(4);
else if(iNumParams == 5) bOk = __CALL_FIT(5);
else if(iNumParams == 6) bOk = __CALL_FIT(6);
else if(iNumParams == 7) bOk = __CALL_FIT(7);
else if(iNumParams == 8) bOk = __CALL_FIT(8);
else if(iNumParams == 9) bOk = __CALL_FIT(9);
else if(iNumParams == 10) bOk = __CALL_FIT(10);
else if(iNumParams == 11) bOk = __CALL_FIT(11);
else if(iNumParams == 12) bOk = __CALL_FIT(12);
else if(iNumParams == 13) bOk = __CALL_FIT(13);
else if(iNumParams == 14) bOk = __CALL_FIT(14);
else if(iNumParams == 15) bOk = __CALL_FIT(15);
else tl2::log_err("In ", __func__, ": Invalid number of function arguments.");
// write back fitted values & errors
for(std::size_t iParam=0; iParam<iNumParams; ++iParam)
{
pValues[iParam] = vecVals[iParam];
pErrors[iParam] = vecErrs[iParam];
}
return int(bOk);
}
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
// Minimisation
/**
* internal helper function to call fitter using variable args
*/
template<std::size_t iNumArgs>
static inline bool _invoke_minimise(void *_pFkt,
const std::vector<std::string>& vecParamNames,
std::vector<tl2::t_real_min>& vecVals, std::vector<tl2::t_real_min>& vecErrs,
const std::vector<bool>& vecFixed)
{
auto *pFkt = reinterpret_cast<tl2::t_fkt_vararg<t_real, iNumArgs>>(_pFkt);
return tl2::minimise<iNumArgs>(pFkt, vecParamNames, vecVals, vecErrs, &vecFixed, g_bDebug);
}
/**
* function minimisation
*/
extern "C" int minimise(void *_pFkt, std::size_t iNumParams,
jl_array_t *parrParamNames, jl_array_t *parrFixed,
t_real *pValues, t_real *pErrors)
{
// parameter names
std::vector<std::string> vecParamNames =
tl2::from_jl_arr<std::vector, std::string>(parrParamNames);
// fixed parameters
std::vector<std::string> vecFixedParams =
tl2::from_jl_arr<std::vector, std::string>(parrFixed);
std::vector<bool> vecFixed;
for(const std::string& strParam : vecParamNames)
{
bool bFixed = std::find(vecFixedParams.begin(), vecFixedParams.end(), strParam)
!= vecFixedParams.end();
vecFixed.push_back(bFixed);
}
// values & errors
std::vector<tl2::t_real_min> vecVals, vecErrs;
for(std::size_t iIdx=0; iIdx<iNumParams; ++iIdx)
{
vecVals.push_back(pValues[iIdx]);
vecErrs.push_back(pErrors[iIdx]);
}
// ------------------------------------------------------------------------
// fill up missing parameters and hints
if(vecParamNames.size() < iNumParams)
{
for(std::size_t iArg=vecParamNames.size(); iArg<iNumParams; ++iArg)
{
std::ostringstream ostrArg;
ostrArg << "arg_" << iArg;
vecParamNames.push_back(ostrArg.str());
}
}
while(vecVals.size() < iNumParams) vecVals.push_back(0);
while(vecErrs.size() < iNumParams) vecErrs.push_back(0);
while(vecFixed.size() < iNumParams) vecFixed.push_back(0);
// ------------------------------------------------------------------------
#define __CALL_MINI(NUM) _invoke_minimise<NUM>(_pFkt, vecParamNames, vecVals, vecErrs, vecFixed)
bool bOk = 0;
// stating all needed specialisations of the fit template function
if(iNumParams == 1) bOk = __CALL_MINI(1);
else if(iNumParams == 2) bOk = __CALL_MINI(2);
else if(iNumParams == 3) bOk = __CALL_MINI(3);
else if(iNumParams == 4) bOk = __CALL_MINI(4);
else if(iNumParams == 5) bOk = __CALL_MINI(5);
else if(iNumParams == 6) bOk = __CALL_MINI(6);
else if(iNumParams == 7) bOk = __CALL_MINI(7);
else if(iNumParams == 8) bOk = __CALL_MINI(8);
else if(iNumParams == 9) bOk = __CALL_MINI(9);
else if(iNumParams == 10) bOk = __CALL_MINI(10);
else if(iNumParams == 11) bOk = __CALL_MINI(11);
else if(iNumParams == 12) bOk = __CALL_MINI(12);
else if(iNumParams == 13) bOk = __CALL_MINI(13);
else if(iNumParams == 14) bOk = __CALL_MINI(14);
else if(iNumParams == 15) bOk = __CALL_MINI(15);
else tl2::log_err("In ", __func__, ": Invalid number of function arguments.");
// write back fitted values & errors
for(std::size_t iParam=0; iParam<iNumParams; ++iParam)
{
pValues[iParam] = vecVals[iParam];
pErrors[iParam] = vecErrs[iParam];
}
return int(bOk);
}
/**
* tlibs2
* julia interface helpers
* @author Tobias Weber <tobias.weber@tum.de>, <tweber@ill.fr>
* @date 2017 -- 2018
* @license GPLv3, see 'LICENSE' file
* @desc Forked on 7-Nov-2018 from the privately and TUM-PhD-developed "tlibs" project (https://github.com/t-weber/tlibs).
*/
#ifndef __TLIBS2_JL_H__
#define __TLIBS2_JL_H__
#include <julia.h>
#include <tuple>
#include <limits>
#include <string>
namespace tl2
{
// ----------------------------------------------------------------------------
/**
* Julia data type traits
*/
template<typename T> struct jl_traits {};
template<> struct jl_traits<int64_t>
{
using value_type = int64_t;
static jl_datatype_t* get_type() { return jl_int64_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_int64(pVal); }
static jl_value_t* box(value_type val) { return jl_box_int64(val); }
};
template<> struct jl_traits<uint64_t>
{
using value_type = uint64_t;
static jl_datatype_t* get_type() { return jl_uint64_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_uint64(pVal); }
static jl_value_t* box(value_type val) { return jl_box_uint64(val); }
};
template<> struct jl_traits<int32_t>
{
using value_type = int32_t;
static jl_datatype_t* get_type() { return jl_int32_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_int32(pVal); }
static jl_value_t* box(value_type val) { return jl_box_int32(val); }
};
template<> struct jl_traits<uint32_t>
{
using value_type = uint32_t;
static jl_datatype_t* get_type() { return jl_uint32_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_uint32(pVal); }
static jl_value_t* box(value_type val) { return jl_box_uint32(val); }
};
template<> struct jl_traits<int16_t>
{
using value_type = int16_t;
static jl_datatype_t* get_type() { return jl_int16_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_int16(pVal); }
static jl_value_t* box(value_type val) { return jl_box_int16(val); }
};
template<> struct jl_traits<uint16_t>
{
using value_type = uint16_t;
static jl_datatype_t* get_type() { return jl_uint16_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_uint16(pVal); }
static jl_value_t* box(value_type val) { return jl_box_uint16(val); }
};
template<> struct jl_traits<int8_t>
{
using value_type = int8_t;
static jl_datatype_t* get_type() { return jl_int8_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_int8(pVal); }
static jl_value_t* box(value_type val) { return jl_box_int8(val); }
};
template<> struct jl_traits<uint8_t>
{
using value_type = uint8_t;
static jl_datatype_t* get_type() { return jl_uint8_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_uint8(pVal); }
static jl_value_t* box(value_type val) { return jl_box_uint8(val); }
};
template<> struct jl_traits<float>
{
using value_type = float;
static jl_datatype_t* get_type() { return jl_float32_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_float32(pVal); }
static jl_value_t* box(value_type val) { return jl_box_float32(val); }
};
template<> struct jl_traits<double>
{
using value_type = double;
static jl_datatype_t* get_type() { return jl_float64_type; }
static value_type unbox(jl_value_t *pVal) { return jl_unbox_float64(pVal); }
static jl_value_t* box(value_type val) { return jl_box_float64(val); }
};
template<> struct jl_traits<std::string>
{
using value_type = std::string;
static jl_datatype_t* get_type() { return jl_string_type; }
static value_type unbox(jl_value_t *pVal)
{
const char* pc = jl_string_ptr(pVal);
return std::string(pc);
}
static jl_value_t* box(value_type val) { return jl_cstr_to_string(val.c_str()); }
};
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
/**
* converts an stl container of containers into a julia 2d array
*/
template<template<class...> class t_cont/*=std::vector*/, class T/*=double*/>
jl_array_t* make_jl_2darr(const t_cont<t_cont<T>>& vecvec)
{
// number of columns and rows
std::size_t iCols = vecvec.size();
std::size_t iRows = std::numeric_limits<std::size_t>::max();
for(const auto& vec : vecvec)
iRows = std::min(iRows, vec.size());
if(!iCols) iRows = 0;
jl_array_t *pArr = jl_alloc_array_2d(
jl_apply_array_type(reinterpret_cast<jl_value_t*>(jl_traits<T>::get_type()), 2),
iRows, iCols);
T* pDat = reinterpret_cast<T*>(jl_array_data(pArr));
std::size_t iCurCol = 0;
for(const auto& vec : vecvec)
{
std::size_t iCurRow = 0;
for(const auto& d : vec)
{
pDat[iCurCol*iRows + iCurRow] = d;
++iCurRow;
}
++iCurCol;
}
return pArr;
}
/**
* converts an stl container of strings into a julia array of strings
*/
template<template<class...> class t_cont/*=std::vector*/, class t_str/*=std::string*/>
jl_array_t* make_jl_str_arr(const t_cont<t_str>& vecStr)
{
jl_array_t *pArr = jl_alloc_array_1d(jl_apply_array_type(
reinterpret_cast<jl_value_t*>(jl_string_type), 1), vecStr.size());
jl_value_t** pDat = reinterpret_cast<jl_value_t**>(jl_array_data(pArr));
std::size_t iIdx = 0;
for(const t_str& str : vecStr)
{
pDat[iIdx] = jl_cstr_to_string(str.c_str());
++iIdx;
}
return pArr;
}
/**
* converts a julia array into an stl container
*/
template<template<class...> class t_cont/*=std::vector*/, class t_val>
t_cont<t_val> from_jl_arr(jl_array_t *pArr, std::size_t iSkipFront = 0)
{
const std::size_t iSize = jl_array_len(pArr);
t_cont<t_val> vecRet;
vecRet.reserve(iSize);
for(std::size_t iElem=iSkipFront; iElem<iSize; ++iElem)
{
jl_value_t* pVal = jl_arrayref(pArr, iElem);
t_val val = jl_traits<t_val>::unbox(pVal);
vecRet.push_back(val);
}
return vecRet;
}
/**
* converts a map of strings into two julia arrays of strings (key & value)
*/
template<template<class...> class t_cont/*=std::map*/, class t_str/*=std::string*/>
std::tuple<jl_array_t*, jl_array_t*> make_jl_strmap_arr(const t_cont<t_str, t_str>& map)
{
jl_array_t *pArrKey = jl_alloc_array_1d(
jl_apply_array_type(reinterpret_cast<jl_value_t*>(jl_string_type), 1), map.size());
jl_array_t *pArrVal = jl_alloc_array_1d(
jl_apply_array_type(reinterpret_cast<jl_value_t*>(jl_string_type), 1), map.size());
jl_value_t** pDatKey = reinterpret_cast<jl_value_t**>(jl_array_data(pArrKey));
jl_value_t** pDatVal = reinterpret_cast<jl_value_t**>(jl_array_data(pArrVal));
std::size_t iIdx = 0;
for(const auto& pair : map)
{
pDatKey[iIdx] = jl_cstr_to_string(pair.first.c_str());
pDatVal[iIdx] = jl_cstr_to_string(pair.second.c_str());
++iIdx;
}
return std::make_tuple(pArrKey, pArrVal);
}
// ----------------------------------------------------------------------------
}
#endif
......@@ -9,8 +9,6 @@
#include "log.h"
#include <iomanip>
#include <cstring>
#include <chrono>
#include <boost/date_time/c_time.hpp>
......
......@@ -10,6 +10,9 @@
#ifndef __TLIBS2_LOGGER_H__
#define __TLIBS2_LOGGER_H__
#include <cstring>
#include <cmath>
#include <ctime>
#include <vector>
#include <unordered_map>
#include <iostream>
......@@ -17,15 +20,11 @@
#include <thread>
#include <mutex>
#include <utility>
#include <chrono>
#include <ctime>
#include <cmath>
#include <string>
#include <exception>
#include <boost/type_index.hpp>
#include "algos.h"
#include <boost/type_index.hpp>
namespace tl2 {
......