file Elements/daFunk.hpp #
[No description available]
Namespaces #
Classes #
Defines #
Macros Documentation #
define DEF_FUNKTRAIT #
#define DEF_FUNKTRAIT(
C
)
class C { \
public: \
static daFunk::FunkPlain* ptr; \
static void set(daFunk::FunkPlain* new_ptr) \
{ \
delete ptr; \
ptr = new_ptr; \
} \
}; \
daFunk::FunkPlain* C::ptr = NULL;
define MATH_OPERATION #
#define MATH_OPERATION(
OPERATION
)
class FunkMath_##OPERATION: public FunkBase \
{ \
public: \
FunkMath_##OPERATION(Funk f) \
{ \
functions = vec(f); \
arguments = f->getArgs(); \
singularities = f->getSingl(); \
} \
double value(const std::vector<double> & data, size_t bindID) \
{ \
return OPERATION(functions[0]->value(data, bindID)); \
} \
}; \
inline Funk OPERATION (Funk f) { return Funk(new FunkMath_##OPERATION(f)); }
define MATH_OPERATION #
#define MATH_OPERATION(
OPERATION,
SYMBOL
)
define MATH_OPERATION #
#define MATH_OPERATION(
OPERATION
)
class FunkMath_##OPERATION: public FunkBase \
{ \
public: \
FunkMath_##OPERATION(Funk f) \
{ \
functions = vec(f); \
arguments = f->getArgs(); \
singularities = f->getSingl(); \
} \
double value(const std::vector<double> & data, size_t bindID) \
{ \
return OPERATION(functions[0]->value(data, bindID)); \
} \
}; \
inline Funk OPERATION (Funk f) { return Funk(new FunkMath_##OPERATION(f)); }
Source code #
/*
* _______
* \ ___ `'. _..._ .
* ' |--.\ \ _.._ .' '. .'|
* | | \ ' .' .._| . .-. . .' |
* | | | ' __ | ' | ' ' |< |
* | | | | .:--.'. __| |__ _ _ | | | | | | ____
* | | ' .'/ | \ ||__ __|| ' / | | | | | | | \ .'
* | |___.' /' `" __ | | | | .' | .' | | | | | | |/ .
* /_______.'/ .'.''| | | | / | / | | | | | | /\ \
* \_______|/ / / | |_ | | | `'. | | | | | | | \ \
* \ \._,\ '/ | | ' .'| '/| | | | ' \ \ \
* `--' `" |_| `-' `--' '--' '--''------' '---'
*
* daFunk - dynamisch allokierbare Funktionen
*
* v0.1 Dec 2014
* v0.2 Mar 2015 - Completely rewritten internal structure
* v0.3 May 2016 - Extensions
*
* Christoph Weniger, created Dec 2014, edited until May 2016
* <c.weniger@uva.nl>
*
* with contributions related to thread-safety from
* Lars A. Dal, Apr, Jun 2015
* <l.a.dal@fys.uio.no>
*
*
*
* The MIT License (MIT)
*
* Copyright (c) 2016 Christoph Weniger
*
* Permission is hereby granted, free of charge, to any person obtaining a copy of
* this software and associated documentation files (the "Software"), to deal in
* the Software without restriction, including without limitation the rights to
* use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
* of the Software, and to permit persons to whom the Software is furnished to do
* so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in all
* copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
#ifndef __DAFUNK_HPP__
#define __DAFUNK_HPP__
#include <iostream>
#include <vector>
#include <algorithm>
#include <map>
#include <set>
#include <cmath>
#include <stdexcept>
//#define NDEBUG
#include <assert.h>
//#define GAMBIT_DIR
#ifdef GAMBIT_DIR
#include <omp.h>
#include "gambit/Utils/standalone_error_handlers.hpp"
#include "gambit/Utils/util_macros.hpp"
#endif
#include "boost/shared_ptr.hpp"
#include "boost/enable_shared_from_this.hpp"
//#include <memory>
//using std::shared_ptr;
//using std::enable_shared_from_this;
namespace daFunk {class FunkPlain;}
#define DEF_FUNKTRAIT(C) \
class C { \
public: \
static daFunk::FunkPlain* ptr; \
static void set(daFunk::FunkPlain* new_ptr) \
{ \
delete ptr; \
ptr = new_ptr; \
} \
}; \
daFunk::FunkPlain* C::ptr = NULL;
// Extensions
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
namespace daFunk
{
//
// Type declarations etc.
//
using boost::shared_ptr;
using boost::enable_shared_from_this;
using boost::dynamic_pointer_cast;
using boost::static_pointer_cast;
class FunkBase;
class FunkBound;
class FunkIntegrate_gsl1d;
typedef shared_ptr<FunkBase> Funk;
typedef shared_ptr<FunkBound> BoundFunk;
typedef std::vector<std::string> ArgsType;
typedef std::map<std::string, std::vector<std::pair<Funk, Funk>>> Singularities;
//template <typename... Args>
//using PlainPtr = double(*)(Args&...);
typedef double(*PlainPtr1)(double&);
typedef double(*PlainPtr2)(double&,double&);
typedef double(*PlainPtr3)(double&,double&,double&);
typedef double(*PlainPtr4)(double&,double&,double&,double&);
//template <typename... Args>
//using PlainPtrs = std::pair<double(*)(Args...,void*), void*>;
typedef std::pair<double(*)(double,void*), void*> PlainPtrs1;
typedef std::pair<double(*)(double,double,void*), void*> PlainPtrs2;
typedef std::pair<double(*)(double,double,double,void*), void*> PlainPtrs3;
typedef std::pair<double(*)(double,double,double,double,void*), void*> PlainPtrs4;
//
// Vector initialization from argument list
//
// Usage: std::vector<T> v = vec<T>(v1, v2, v3, ...);
//
template <typename T>
inline std::vector<T> vec(std::vector<T> vector)
{
return vector;
}
template <typename T, typename... Args>
inline std::vector<T> vec(std::vector<T> vector, T value, Args... args)
{
vector.push_back(value);
return vec<T>(vector, args...);
}
template <typename T, typename... Args>
inline std::vector<T> vec(T value, Args... args)
{
std::vector<T> vector;
vector.push_back(value);
vector = vec<T>(vector, args...);
return vector;
}
template <typename T>
inline std::vector<T> vec()
{
std::vector<T> vector;
return vector;
}
//
// Math helper functions
//
// Generate linear 1-dim grid
inline std::vector<double> linspace(double x0, double x1, unsigned int n)
{
std::vector<double> ret;
double dx = 0;
if (n > 1)
dx = (x1-x0)/(n-1);
for (unsigned int i = 0; i<n; i++)
{
ret.push_back(x0 + i * dx);
}
return ret;
}
// Generate logarithmic 1-dim grid
inline std::vector<double> logspace(double x0, double x1, unsigned int n)
{
std::vector<double> ret;
double dx = 0;
if (n > 1)
dx = (x1-x0)/(n-1);
for (unsigned int i = 0; i<n; i++)
{
ret.push_back(pow(10, x0 + i * dx));
}
return ret;
}
//
// Helper functions for internal calculations
//
inline bool args_match(ArgsType args1, ArgsType args2)
{
bool m = true;
for ( auto it = args1.begin(); it!=args1.end(); it++ )
{
if ( std::find(args2.begin(), args2.end(), *it) == args2.end() )
m = false;
}
if ( args1.size() != args2.size() ) m = false;
return m;
}
inline std::string args_string(ArgsType args)
{
std::string msg = "";
for ( auto it = args.begin(); it!=args.end(); it++ )
{
msg += *it;
if ( it != args.end() - 1)
msg += ", ";
}
return msg;
}
inline ArgsType joinArgs(ArgsType args1, ArgsType args2)
{
args1.insert(args1.end(), args2.begin(), args2.end());
std::set<std::string> argsset(args1.begin(), args1.end());
args1.assign(argsset.begin(), argsset.end());
return args1;
}
inline ArgsType eraseArg(ArgsType args, std::string arg)
{
auto it = std::find(args.begin(), args.end(), arg);
if (it!=args.end()) args.erase(it);
return args;
}
inline Singularities joinSingl(Singularities s1, Singularities s2)
{
for ( auto it = s2.begin(); it != s2.end(); ++it )
{
if ( s1.find(it->first) == s1.end() )
s1[it->first] = it->second;
else
{
for ( auto it2 = it->second.begin(); it2 != it->second.end(); ++it2 )
{
if ( std::find(s1[it->first].begin(), s1[it->first].end(), *it2) == s1[it->first].end() )
s1[it->first].push_back(*it2);
}
}
}
return s1;
}
//
// Index lists (taken from stackoverflow)
//
// The structure that encapsulates index lists
template <size_t... Is>
struct index_list
{
};
// Collects internal details for generating index ranges [MIN, MAX)
namespace detail
{
// Declare primary template for index range builder
template <size_t MIN, size_t N, size_t... Is>
struct range_builder;
// Base step
template <size_t MIN, size_t... Is>
struct range_builder<MIN, MIN, Is...>
{
typedef index_list<Is...> type;
};
// Induction step
template <size_t MIN, size_t N, size_t... Is>
struct range_builder : public range_builder<MIN, N - 1, N - 1, Is...>
{
};
}
//
// **** The central virtual base class ****
//
class FunkBase: public enable_shared_from_this<FunkBase>
{
public:
FunkBase() {};
virtual ~FunkBase() {};
// Standard setting handles
template <typename... Args> Funk set(std::string arg, Funk g, Args... args);
template <typename... Args> Funk set(std::string arg, double x, Args... args);
template <typename... Args> Funk set(std::string arg, std::string arg1, Args... args);
template <typename... Args> Funk set();
// Standard binding handles
template <typename... Args> shared_ptr<FunkBound> bind(Args... args);
// Convenience functions
const std::vector<std::string> & getArgs() { return this->arguments; };
std::size_t getNArgs() {return this->arguments.size();};
bool hasArg(std::string);
bool hasArgs();
Funk help();
template <typename... Args> bool assert_args(Args... args);
// Return value & standard resolve
virtual double value(const std::vector<double> &, size_t bindID) = 0;
// datamap maps the required function arguments onto the specific
// entries in the double-valued data array generated by eval().
// datalen acts like a pointer on the last entry of that array, and
// can be increased if more workspace is required. bindID ensures
// that resolution for various different binds can happen in
// parallel with the same Funk objects.
virtual void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap);
// Singularities handling
Singularities getSingl() { return singularities; }
private: // Note: This works currently only for constant functions
Funk set_singularity(std::string arg, Funk pos, Funk width);
Funk set_singularity(std::string arg, double pos, Funk width);
Funk set_singularity(std::string arg, Funk pos, double width);
public:
Funk set_singularity(std::string arg, double pos, double width);
// Print message
Funk print(std::string arg);
//
// "External" functions
//
// Integration
template <typename... Args> shared_ptr<FunkIntegrate_gsl1d> gsl_integration(Args... args);
// Plain function generators
PlainPtrs1 plain(std::string);
PlainPtrs2 plain(std::string, std::string);
PlainPtrs3 plain(std::string, std::string, std::string);
PlainPtrs4 plain(std::string, std::string, std::string, std::string);
template <typename T>
PlainPtr1 plain(std::string);
template <typename T>
PlainPtr2 plain(std::string, std::string);
template <typename T>
PlainPtr3 plain(std::string, std::string, std::string);
template <typename T>
PlainPtr4 plain(std::string, std::string, std::string, std::string);
protected:
std::vector<Funk> functions; // Dependent functions
ArgsType arguments; // Argument names
std::vector<std::vector<size_t> > indices; // Indices for data object
size_t datalen;
Singularities singularities;
};
// A vector class with global knowledge about its health status.
// (BoundFunk objects are occasionally destructed *after* livingVector has
// been destructed, causing segfaults if not catched properly.)
class livingVector: public std::vector<size_t>
{
public:
~livingVector()
{
livingVector::is_dead() = true;
}
static bool & is_dead()
{
static bool dead = false;
return dead;
}
};
class FunkBound
{
public:
FunkBound(Funk f, size_t datalen, size_t bindID) : f(f), datalen(datalen), bindID(bindID) {};
~FunkBound() {bindID_manager(bindID,false);};
double value(std::vector<double> & map, size_t bindID) {(void)bindID; (void)map; return 0;};
template <typename... Args> inline double eval(Args... argss)
{
auto data = vec<double>(argss...);
data.resize(datalen);
return f->value(data, bindID);
}
template <typename... Args> inline std::vector<double> vect(Args... argss)
{
std::vector<std::vector<double>> coll;
return this->vect2(coll, argss...);
}
private:
template <typename... Args>
friend shared_ptr<FunkBound> FunkBase::bind(Args... argss);
// Function for generating unique bindIDs.
// IDs are sequential, starting from 0.
static void bindID_manager(size_t &bindID, bool bind)
{
static size_t n_idx = 0;
static livingVector free;
if(bind)
{
#pragma omp critical (bindID_allocation)
{
if(free.empty())
{
// Increment after assignment (n_idx++ instead of ++n_idx)
bindID = n_idx++;
}
else
{
bindID = free.back();
free.pop_back();
}
}
}
else
{
if ( not livingVector::is_dead() )
{
#pragma omp critical (bindID_allocation)
free.push_back(bindID);
}
}
}
template <typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll)
{
size_t size = 1;
std::vector<bool> vec_flag;
for ( auto it = coll.begin(); it != coll.end(); ++it )
{
if ( it->size() == 1 )
{
vec_flag.push_back(false);
continue;
}
vec_flag.push_back(true);
if ( size == 1 ) size = it->size();
if ( size != it->size() )
{
std::cerr << "daFunk::FunkBase WARNING: Inconsistent vector lengths." << std::endl;
return vec<double>();
}
}
auto r = vec<double>();
auto data = vec<double>();
data.resize(datalen);
for ( size_t i = 0; i != size; ++i )
{
for ( size_t j = 0; j != coll.size(); ++j )
{
if ( vec_flag[j] )
data[j] = coll[j][i];
else
data[j] = coll[j][0];
}
r.push_back(f->value(data, bindID));
}
return r;
}
template <typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll, double x, Args... argss)
{
coll.push_back(vec<double>(x));
return this->vect2(coll, argss...);
}
template <typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll, std::vector<double> v, Args... argss)
{
coll.push_back(v);
return this->vect2(coll, argss...);
}
Funk f; // bound function
// datalen is the length of the double-valued data array that is
// needed as workspace for function evaluation, and that is created
// on the heap for each eval separately to ensure thread-safety.
size_t datalen;
// bindID has the purpose of allowing bound functions (instances of
// FunkBase and daughter classes) to be bound by various binding
// functions simultaneously.
size_t bindID;
};
//
// Derived class with (templated) static member functions as plain function
// prototypes.
//
class FunkPlain: public FunkBase
{
public:
FunkPlain(Funk fin, std::string arg1) : f(fin->bind(arg1)) {}
FunkPlain(Funk fin, std::string arg1, std::string arg2) : f(fin->bind(arg1, arg2)) {}
FunkPlain(Funk fin, std::string arg1, std::string arg2, std::string arg3) : f(fin->bind(arg1, arg2, arg3)) {}
FunkPlain(Funk fin, std::string arg1, std::string arg2, std::string arg3, std::string arg4) : f(fin->bind(arg1, arg2, arg3, arg4)) {}
static double plain1p(double x1, void* ptr)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
return funkPtrPtr->f->eval(x1);
}
static double plain2p(double x1, double x2, void* ptr)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
return funkPtrPtr->f->eval(x1, x2);
}
static double plain3p(double x1, double x2, double x3, void* ptr)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
return funkPtrPtr->f->eval(x1, x2, x3);
}
static double plain4p(double x1, double x2, double x3, double x4, void* ptr)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
return funkPtrPtr->f->eval(x1, x2, x3, x4);
}
template <typename T>
static double plain1(double& x1)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
return funkPtrPtr->f->eval(x1);
}
template <typename T>
static double plain2(double& x1, double& x2)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
return funkPtrPtr->f->eval(x1, x2);
}
template <typename T>
static double plain3(double& x1, double& x2, double& x3)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
return funkPtrPtr->f->eval(x1, x2, x3);
}
template <typename T>
static double plain4(double& x1, double& x2, double& x3, double& x4)
{
FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
return funkPtrPtr->f->eval(x1, x2, x3, x4);
}
double value(const std::vector<double> & args, size_t bindID)
{
(void)args;
(void)bindID;
assert ( 0 == 1 ); // This function should never be called
return 0;
}
private:
shared_ptr<FunkBound> f; // bound function
std::string arg1, arg2, arg3, arg4;
};
//
// Derived class that implements constant
//
class FunkConst: public FunkBase
{
public:
template <typename... Args>
FunkConst(double c, Args ...argss) : c(c) { arguments = vec<std::string>(argss...); }
FunkConst(double c) : c(c) { arguments.resize(0); }
double value(const std::vector<double> & data, size_t bindID)
{
(void)data;
(void)bindID;
return c;
}
private:
double c;
};
template <typename... Args>
inline Funk one(Args... argss) { return Funk(new FunkConst(1., argss...)); }
template <typename... Args>
inline Funk zero(Args... argss) { return Funk(new FunkConst(0., argss...)); }
template <typename... Args>
inline Funk cnst(double x, Args... argss) { return Funk(new FunkConst(x, argss...)); }
//
// Derived class that implements setting of parameters
//
class FunkDerived: public FunkBase
{
public:
FunkDerived(Funk f, std::string arg, Funk g) : my_arg(arg)
{
setup(f, arg, g);
};
FunkDerived(Funk f, std::string arg, double x) : my_arg(arg)
{
setup(f, arg, cnst(x));
}
// We need to sneak in an additional parameter
void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap)
{
functions[1]->resolve(datamap, datalen, bindID, argmap); // resolve g
// add new slot for result from of g
if(my_index.size() <= bindID)
{
my_index.resize(bindID+1);
// TODO: Introduce informative error message if bind fails
// because of inconsistencies, e.g. along the lines:
//
// std::cout << "FATAL ERROR: bind() attempts to resolve dependencies for FunkDerived" << std::endl;
// std::cout << "object inconsistently." << std::endl;
// std::cout << "Encountered while resolving " << my_arg << " and:" << std::endl;
// for (auto it = datamap.begin(); it != datamap.end(); it++)
// {
// std::cout << " " << it->first << std::endl;
// }
// exit(1);
}
if ( argmap.find(my_arg) == argmap.end() )
{
if(datamap.find(my_arg) == datamap.end())
{
my_index[bindID] = datalen;
argmap[my_arg] = datalen;
++datalen;
}
else
{
my_index[bindID] = datamap[my_arg];
}
}
else
{
my_index[bindID] = argmap[my_arg];
}
datamap[my_arg] = my_index[bindID]; // add or overwrite entry in datamap
functions[0]->resolve(datamap, datalen, bindID, argmap); // resolve f
}
double value(const std::vector<double> & data, size_t bindID)
{
std::vector<double> data2(data);
data2[my_index[bindID]] = functions[1]->value(data, bindID);
return functions[0]->value(data2, bindID);
}
private:
std::string my_arg;
void setup(Funk f, std::string arg, Funk g)
{
functions = vec(f, g);
Singularities tmp_singl = f->getSingl();
if ( tmp_singl.erase(arg) > 0 )
std::cout << "daFunk::FunkBase WARNING: Loosing singularity information while setting " << arg << std::endl;
singularities = joinSingl(g->getSingl(), tmp_singl);
arguments = joinArgs(eraseArg(f->getArgs(), arg), g->getArgs());
};
std::vector<size_t> my_index;
};
//
// Derived class for the import of plain functions
//
template <bool threadsafe, typename... funcargs>
class FunkFunc : public FunkBase
{
public:
template <typename... Args>
FunkFunc(double (*f)(funcargs...), Args... argss)
{
ptr = f;
digest_input(argss...);
}
double value(const std::vector<double> & data, size_t bindID)
{
std::tuple<typename std::remove_reference<funcargs>::type...> my_input;
double result;
size_t i = 0;
#pragma omp critical (FunkFunc_setInput)
{
for ( auto f = functions.begin(); f != functions.end(); ++f, ++i)
{
*map[i] = (*f)->value(data, bindID);
}
my_input = input;
}
if(threadsafe)
{
result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
}
else
{
#pragma omp critical (FunkFunc_externalFunctionCall)
{
result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
}
}
return result;
}
template <size_t... Args>
double ppp(index_list<Args...>, std::tuple<typename std::remove_reference<funcargs>::type...> & my_input)
{
return (*ptr)(std::get<Args>(my_input)...);
}
private:
std::tuple<typename std::remove_reference<funcargs>::type...> input;
std::vector<double*> map;
double (*ptr)(funcargs...);
// Digest input parameters
// (forwarding everything except daFunk::Funk types, which is mapped onto
// funktion parameters)
template<typename T, typename... Args>
void digest_input(T x, Args... argss)
{
const int i = sizeof...(funcargs) - sizeof...(argss) - 1;
std::get<i>(input) = x;
digest_input(argss...);
}
template<typename... Args>
void digest_input(Funk f, Args... argss)
{
const int i = sizeof...(funcargs) - sizeof...(argss) - 1;
map.push_back(&std::get<i>(input));
arguments = joinArgs(arguments, f->getArgs());
functions.push_back(f);
singularities = joinSingl(singularities, f->getSingl());
digest_input(argss...);
}
void digest_input() {};
};
template <typename... funcargs, typename... Args>
Funk func(double (*f)(funcargs...), Args... args) {
return Funk(new FunkFunc<false, funcargs...>(f, args...));
}
// Version that assumes the function to be threadsafe
template <typename... funcargs, typename... Args>
Funk func_fromThreadsafe(double (*f)(funcargs...), Args... args) {
return Funk(new FunkFunc<true, funcargs...>(f, args...));
}
template <bool threadsafe, typename O, typename... funcargs>
class FunkFuncM : public FunkBase
{
public:
template <typename... Args>
FunkFuncM(O* obj, double (O::* f)(funcargs...), Args... argss) : obj(obj)
{
ptr = f;
digest_input(argss...);
}
template <typename... Args>
FunkFuncM(shared_ptr<O> obj, double (O::* f)(funcargs...), Args... argss) : shared_obj(obj), obj(&*obj)
{
ptr = f;
digest_input(argss...);
}
double value(const std::vector<double> & data, size_t bindID)
{
std::tuple<typename std::remove_reference<funcargs>::type...> my_input;
double result;
size_t i = 0;
#pragma omp critical(FunkFuncM_value)
{
for ( auto f = functions.begin(); f != functions.end(); ++f, ++i)
{
*map[i] = (*f)->value(data, bindID);
}
my_input = input;
}
if(threadsafe)
{
result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
}
else
{
#pragma omp critical (FunkFuncM_objectFunctionCall)
{
result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
}
}
return result;
}
template <size_t... Args>
double ppp(index_list<Args...>, std::tuple<typename std::remove_reference<funcargs>::type...> & my_input)
{
return (*obj.*ptr)(std::get<Args>(my_input)...);
}
private:
std::tuple<typename std::remove_reference<funcargs>::type...> input;
std::vector<double*> map;
double (O::* ptr)(funcargs...);
shared_ptr<O> shared_obj;
O* obj;
// Digest input parameters
// (forwarding everything except daFunk::Funk types, which is mapped onto
// funktion parameters)
template<typename T, typename... Args>
void digest_input(T x, Args... argss)
{
const int i = sizeof...(funcargs) - sizeof...(argss) - 1;
std::get<i>(input) = x;
digest_input(argss...);
}
template<typename... Args>
void digest_input(Funk f, Args... argss)
{
const int i = sizeof...(funcargs) - sizeof...(argss) - 1;
map.push_back(&std::get<i>(input));
arguments = joinArgs(arguments, f->getArgs());
functions.push_back(f);
singularities = joinSingl(singularities, f->getSingl());
digest_input(argss...);
}
void digest_input() {};
};
template <typename O, typename... funcargs, typename... Args>
Funk funcM(O* obj, double (O::* f)(funcargs...), Args... args) {
return Funk(new FunkFuncM<false, O, funcargs...>(obj, f, args...));
}
template <typename O, typename... funcargs, typename... Args>
Funk funcM(shared_ptr<O> obj, double (O::* f)(funcargs...), Args... args) {
return Funk(new FunkFuncM<false, O, funcargs...>(obj, f, args...));
}
// Versions that assume the object to be threadsafe
template <typename O, typename... funcargs, typename... Args>
Funk funcM_fromThreadsafe(O* obj, double (O::* f)(funcargs...), Args... args) {
return Funk(new FunkFuncM<true, O, funcargs...>(obj, f, args...));
}
template <typename O, typename... funcargs, typename... Args>
Funk funcM_fromThreadsafe(shared_ptr<O> obj, double (O::* f)(funcargs...), Args... args) {
return Funk(new FunkFuncM<true, O, funcargs...>(obj, f, args...));
}
//
// Derived class that implements delta function
//
class FunkDelta: public FunkBase
{
public:
FunkDelta(std::string arg, double pos, double width) : pos(pos), width(width)
{
arguments = vec(arg);
singularities[arg].push_back(std::pair<Funk, Funk>(cnst(pos), cnst(width)));
}
double value(const std::vector<double> & data, size_t bindID)
{
double x = data[indices[bindID][0]];
return exp(-pow(x-pos,2)/pow(width,2)/2)/sqrt(2*M_PI)/width;
}
private:
double pos, width;
};
inline Funk delta(std::string arg, double pos, double width) { return Funk(new FunkDelta(arg, pos, width)); }
//
// Derived class that implements simple linear variable
//
class FunkVar: public FunkBase
{
public:
FunkVar(std::string arg)
{
arguments = vec(arg);
}
double value(const std::vector<double> & data, size_t bindID)
{
return data[indices[bindID][0]];
}
};
inline Funk var(std::string arg) { return Funk(new FunkVar(arg)); }
//
// Definition of FunkBase member functions
//
inline Funk FunkBase::set_singularity(std::string arg, Funk pos, Funk width)
{
singularities[arg].push_back(std::pair<Funk, Funk>(pos, width));
return shared_from_this();
};
inline Funk FunkBase::set_singularity(std::string arg, double pos, Funk width)
{ return shared_from_this()->set_singularity(arg, cnst(pos), width); }
inline Funk FunkBase::set_singularity(std::string arg, double pos, double width)
{ return shared_from_this()->set_singularity(arg, cnst(pos), cnst(width)); }
inline Funk FunkBase::set_singularity(std::string arg, Funk pos, double width)
{ return shared_from_this()->set_singularity(arg, pos, cnst(width)); }
inline void FunkBase::resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap)
{
// Resolve my dependencies
auto it1 = arguments.begin();
std::vector<size_t> tmp_indices;
tmp_indices.resize(arguments.size());
auto it2 = tmp_indices.begin();
for (; it1 != arguments.end() && it2 != tmp_indices.end(); ++it1, ++it2 )
{
try
{
*it2 = datamap.at(*it1);
}
catch (std::out_of_range & e)
{
std::string msg = "FunkBase::resolve() encountered internal problem when resolving " + *it1 + ".\n";
msg+= " --- Actual arguments of object: " + args_string(arguments);
#ifdef GAMBIT_DIR
Gambit::utils_error().raise(LOCAL_INFO, msg);
#else
throw std::invalid_argument(msg);
#endif
}
}
// Set indices
if(indices.size() <= bindID)
{
indices.resize(bindID+1);
// TODO: Throw error if problems are encountered
//std::cout << "FATAL ERROR: bind() attempts to resolve dependencies of FunkBase" << std::endl;
//std::cout << "object inconsistently." << std::endl;
//std::cout << "Encountered while resolving: " << std::endl;
//for (auto it = datamap.begin(); it != datamap.end(); it++)
//{
// std::cout << " " << it->first << std::endl;
//}
//exit(1);
}
indices[bindID] = tmp_indices;
// Resolve other dependencies
for (auto it = functions.begin(); it != functions.end(); ++it)
{
(*it)->resolve(datamap, datalen, bindID, argmap);
}
}
template <typename... Args> inline bool FunkBase::assert_args(Args... args)
{
std::vector<std::vector<std::string>> list = vec<std::vector<std::string>>(args...);
std::set<std::string> myargs(arguments.begin(), arguments.end());
for ( auto it = list.begin(); it != list.end(); ++it )
{
std::set<std::string> theirargs(it->begin(), it->end());
if ( myargs == theirargs )
return true;
}
return false;
}
template <typename... Args> inline Funk FunkBase::set(std::string arg, Funk g, Args... args)
{
Funk f = shared_from_this();
if ( std::find(arguments.begin(), arguments.end(), arg) != arguments.end() )
{
f = Funk(new FunkDerived(f, arg, g));
}
else
{
std::cout << "daFunk::FunkBase WARNING: Ignoring \"" << arg << "\" = function." << std::endl;
}
return f->set(args...);
}
template <typename... Args> inline Funk FunkBase::set(std::string arg, std::string arg1, Args... args)
{ return shared_from_this()->set(arg, var(arg1))->set(args...); }
template <typename... Args> inline Funk FunkBase::set(std::string arg, double x, Args... args)
{ return shared_from_this()->set(arg, cnst(x))->set(args...); }
template <> inline Funk FunkBase::set()
{ return shared_from_this(); }
template <typename... Args> inline shared_ptr<FunkBound> FunkBase::bind(Args... argss)
{
size_t bindID;
FunkBound::bindID_manager(bindID,true);
std::map<std::string, size_t> datamap;
size_t i;
auto bound_arguments = vec<std::string>(argss...);
datalen = bound_arguments.size();
i = 0;
for ( auto it = bound_arguments.begin(); it != bound_arguments.end(); ++it, ++i )
{
datamap[*it] = i;
}
std::map<std::string,size_t> argmap;
if ( not args_match(arguments, bound_arguments) )
{
std::string msg = "FunkBase::bind() tries to resolve wrong arguments.\n";
msg+= " --- Arguments that are supposed to be bound: " + args_string(bound_arguments) + "\n";
msg+= " --- Actual arguments of object: " + args_string(arguments);
#ifdef GAMBIT_DIR
Gambit::utils_error().raise(LOCAL_INFO, msg);
#else
throw std::invalid_argument(msg);
#endif
}
this->resolve(datamap, datalen, bindID, argmap);
return shared_ptr<FunkBound>(new FunkBound(shared_from_this(), datalen, bindID));
}
inline bool FunkBase::hasArg(std::string arg)
{
return ( std::find(arguments.begin(), arguments.end(), arg) != arguments.end() );
}
inline bool FunkBase::hasArgs()
{
return ( this->arguments.size() != 0 );
}
template <typename... Args> inline shared_ptr<FunkIntegrate_gsl1d> FunkBase::gsl_integration(Args... args)
{
return getIntegrate_gsl1d(shared_from_this(), args...);
}
inline Funk FunkBase::help()
{
std::cout << "Arguments:";
for ( auto it = arguments.begin(); it != arguments.end(); it++ )
{
std::cout << " \"" << *it << "\"";
}
if ( arguments.size() == 0 )
std::cout << " none";
std::cout << std::endl;
for ( auto it = singularities.begin() ; it != singularities.end(); ++it )
{
std::cout << "Singularities in " << it->first << ": " << it->second.size() << std::endl;
}
return shared_from_this();
}
inline PlainPtrs1 FunkBase::plain(std::string arg1)
{
void* ptr = new FunkPlain(shared_from_this(), arg1);
return PlainPtrs1(&FunkPlain::plain1p, ptr);
}
inline PlainPtrs2 FunkBase::plain(std::string arg1, std::string arg2)
{
void* ptr = new FunkPlain(shared_from_this(), arg1, arg2);
return PlainPtrs2(&FunkPlain::plain2p, ptr);
}
inline PlainPtrs3 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3)
{
void* ptr = new FunkPlain(shared_from_this(), arg1, arg2, arg3);
return PlainPtrs3(&FunkPlain::plain3p, ptr);
}
inline PlainPtrs4 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3, std::string arg4)
{
void* ptr = new FunkPlain(shared_from_this(), arg1, arg2, arg3, arg4);
return PlainPtrs4(&FunkPlain::plain4p, ptr);
}
template <typename T>
inline PlainPtr1 FunkBase::plain(std::string arg1)
{
T::set(new FunkPlain(shared_from_this(), arg1));
return &FunkPlain::plain1<T>;
}
template <typename T>
inline PlainPtr2 FunkBase::plain(std::string arg1, std::string arg2)
{
T::set(new FunkPlain(shared_from_this(), arg1, arg2));
return &FunkPlain::plain2<T>;
}
template <typename T>
inline PlainPtr3 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3)
{
T::set(new FunkPlain(shared_from_this(), arg1, arg2, arg3));
return &FunkPlain::plain3<T>;
}
template <typename T>
inline PlainPtr4 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3, std::string arg4)
{
T::set(new FunkPlain(shared_from_this(), arg1, arg2, arg3, arg4));
return &FunkPlain::plain4<T>;
}
//
// Mathematical functions from cmath
//
// Unary minus sign
class FunkMath_umin: public FunkBase
{
public:
FunkMath_umin(Funk f)
{
functions = vec(f);
singularities = f->getSingl();
arguments = f->getArgs();
}
double value(const std::vector<double> & data, size_t bindID)
{
return -(functions[0]->value(data, bindID));
}
};
inline Funk operator - (Funk f) { return Funk(new FunkMath_umin(f)); }
// Unary operations
#define MATH_OPERATION(OPERATION) \
class FunkMath_##OPERATION: public FunkBase \
{ \
public: \
FunkMath_##OPERATION(Funk f) \
{ \
functions = vec(f); \
arguments = f->getArgs(); \
singularities = f->getSingl(); \
} \
double value(const std::vector<double> & data, size_t bindID) \
{ \
return OPERATION(functions[0]->value(data, bindID)); \
} \
}; \
inline Funk OPERATION (Funk f) { return Funk(new FunkMath_##OPERATION(f)); }
MATH_OPERATION(cos)
MATH_OPERATION(sin)
MATH_OPERATION(tan)
MATH_OPERATION(acos)
MATH_OPERATION(asin)
MATH_OPERATION(atan)
MATH_OPERATION(cosh)
MATH_OPERATION(sinh)
MATH_OPERATION(tanh)
MATH_OPERATION(acosh)
MATH_OPERATION(asinh)
MATH_OPERATION(atanh)
MATH_OPERATION(exp)
MATH_OPERATION(log)
MATH_OPERATION(log10)
MATH_OPERATION(sqrt)
MATH_OPERATION(fabs)
#undef MATH_OPERATION
// Standard binary operations
#define MATH_OPERATION(OPERATION, SYMBOL) \
class FunkMath_##OPERATION: public FunkBase \
{ \
public: \
FunkMath_##OPERATION(Funk f1, Funk f2) \
{ \
functions = vec(f1, f2); \
arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
singularities = joinSingl(f1->getSingl(), f2->getSingl());\
} \
FunkMath_##OPERATION(double x, Funk f2) \
{ \
auto f1 = cnst(x); \
functions = vec(f1, f2); \
arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
singularities = joinSingl(f1->getSingl(), f2->getSingl());\
} \
FunkMath_##OPERATION(Funk f1, double x) \
{ \
auto f2 = cnst(x); \
functions = vec(f1, f2); \
arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
singularities = joinSingl(f1->getSingl(), f2->getSingl());\
} \
double value(const std::vector<double> & data, size_t bindID) \
{ \
return functions[0]->value(data, bindID) SYMBOL functions[1]->value(data, bindID); \
} \
}; \
inline Funk operator SYMBOL (Funk f1, Funk f2) { return Funk(new FunkMath_##OPERATION(f1, f2)); } \
inline Funk operator SYMBOL (double x, Funk f) { return Funk(new FunkMath_##OPERATION(x, f)); } \
inline Funk operator SYMBOL (Funk f, double x) { return Funk(new FunkMath_##OPERATION(f, x)); }
MATH_OPERATION(Sum,+)
MATH_OPERATION(Mul,*)
MATH_OPERATION(Div,/)
MATH_OPERATION(Dif,-)
#undef MATH_OPERATION
// More binary operations
#define MATH_OPERATION(OPERATION) \
class FunkMath_##OPERATION: public FunkBase \
{ \
public: \
FunkMath_##OPERATION(Funk f1, Funk f2) \
{ \
functions = vec(f1, f2); \
arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
singularities = joinSingl(f1->getSingl(), f2->getSingl()); \
} \
FunkMath_##OPERATION(double x, Funk f2) \
{ \
auto f1 = cnst(x); \
functions = vec(f1, f2); \
arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
singularities = joinSingl(f1->getSingl(), f2->getSingl()); \
} \
FunkMath_##OPERATION(Funk f1, double x) \
{ \
auto f2 = cnst(x); \
functions = vec(f1, f2); \
arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
singularities = joinSingl(f1->getSingl(), f2->getSingl()); \
} \
double value(const std::vector<double> & data, size_t bindID) \
{ \
return OPERATION(functions[0]->value(data, bindID), functions[1]->value(data, bindID)); \
} \
}; \
inline Funk OPERATION (Funk f1, Funk f2) { return Funk(new FunkMath_##OPERATION(f1, f2)); } \
inline Funk OPERATION (double x, Funk f) { return Funk(new FunkMath_##OPERATION(x, f)); } \
inline Funk OPERATION (Funk f, double x) { return Funk(new FunkMath_##OPERATION(f, x)); }
MATH_OPERATION(pow)
MATH_OPERATION(fmin)
MATH_OPERATION(fmax)
#undef MATH_OPERATION
////////////////////////////////////////
// *** End of core implementation ***
////////////////////////////////////////
////////////////////////////////////////
// *** Extensions ***
////////////////////////////////////////
//
// Derived class: 1dim linear or logarithmic interpolation
//
class FunkInterp : public FunkBase
{
public:
FunkInterp(Funk f, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin")
{
setup(f, Xgrid, Ygrid, mode);
}
FunkInterp(std::string arg, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin")
{
setup(var(arg), Xgrid, Ygrid, mode);
}
FunkInterp(double x, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin")
{
setup(cnst(x), Xgrid, Ygrid, mode);
}
double value(const std::vector<double> & data, size_t bindID)
{
functions[0]->value(data, bindID);
return (this->*ptr)(data[indices[bindID][0]]);
}
private:
void setup(Funk f, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode)
{
// TODO: Catch invalid setup
functions = vec(f);
singularities = f->getSingl();
arguments = f->getArgs();
this->Xgrid = Xgrid;
this->Ygrid = Ygrid;
if ( mode == "lin" ) this->ptr = &FunkInterp::linearInterp;
else if ( mode == "log" ) this->ptr = &FunkInterp::logInterp;
}
double logInterp(double x)
{
// Linear interpolation in log-log space
int i = 0; int imax = Xgrid.size() - 1;
if (x<Xgrid[0] or x>Xgrid[imax]) return 0;
for (; i < imax; i++) {if (Xgrid[i] > x) break;};
double x0 = Xgrid[i-1];
double x1 = Xgrid[i];
double y0 = Ygrid[i-1];
double y1 = Ygrid[i];
return y0 * std::exp(std::log(y1/y0) * std::log(x/x0) / std::log(x1/x0));
}
double linearInterp(double x)
{
// Linear interpolation in lin-lin space
int i = 0; int imax = Xgrid.size() - 1;
if (x<Xgrid[0] or x>Xgrid[imax]) return 0;
for (; i < imax; i++) {if (Xgrid[i] > x) break;};
double x0 = Xgrid[i-1];
double x1 = Xgrid[i];
double y0 = Ygrid[i-1];
double y1 = Ygrid[i];
return y0 + (x-x0)/(x1-x0)*(y1-y0);
}
double(FunkInterp::*ptr)(double);
std::vector<double> Xgrid;
std::vector<double> Ygrid;
std::string mode;
};
template <typename T> inline shared_ptr<FunkInterp> interp(T f, std::vector<double> x, std::vector<double> y) { return shared_ptr<FunkInterp>(new FunkInterp(f, x, y)); }
//
// Basic if..else clause
//
class FunkIfElse: public FunkBase
{
public:
FunkIfElse(Funk f, Funk g, Funk h)
{
functions = vec(f, g, h);
singularities = joinSingl(joinSingl(f->getSingl(), g->getSingl()), h->getSingl());
arguments = joinArgs(joinArgs(f->getArgs(), g->getArgs()), h->getArgs());
}
double value(const std::vector<double> & data, size_t bindID)
{
if ( functions[0]->value(data,bindID) >= 0. )
return functions[1]->value(data,bindID);
else
return functions[2]->value(data,bindID);
}
};
inline Funk ifelse(Funk f, Funk g, Funk h) { return Funk(new FunkIfElse(f, g, h)); }
inline Funk ifelse(Funk f, double g, Funk h) { return Funk(new FunkIfElse(f, cnst(g), h)); }
inline Funk ifelse(Funk f, double g, double h) { return Funk(new FunkIfElse(f, cnst(g), cnst(h))); }
inline Funk ifelse(Funk f, Funk g, double h) { return Funk(new FunkIfElse(f, g, cnst(h))); }
//
// Throw errors when called
//
class ThrowError: public FunkBase
{
public:
ThrowError(std::string msg) : msg(msg)
{
}
double value(const std::vector<double> & data, size_t bindID)
{
(void)bindID;
(void)data;
#ifdef GAMBIT_DIR
if ( omp_get_level() == 0 ) // Outside of OMP blocks
{
Gambit::utils_error().raise(LOCAL_INFO, "daFunk::ThrowError says: " + msg);
}
else // Inside of OMP blocks
{
Gambit::piped_errors.request(LOCAL_INFO, "daFunk::ThrowError says: " + msg);
}
#else
throw std::invalid_argument("daFunk::ThrowError says: " + msg);
#endif
return 0;
}
private:
std::string msg; // Error message to throw when function is called
};
inline Funk throwError(std::string msg) { return Funk(new ThrowError(msg)); }
#ifdef GAMBIT_DIR
class RaiseInvalidPoint: public FunkBase
{
public:
RaiseInvalidPoint(std::string msg) : msg(msg)
{
}
double value(const std::vector<double> & data, size_t bindID)
{
(void)bindID;
(void)data;
if ( omp_get_level() == 0 ) // Outside of OMP blocks
{
Gambit::utils_warning().raise(LOCAL_INFO, "daFunk::RaiseInvalidPoint says: " + msg);
Gambit::invalid_point().raise("daFunk::RaiseInvalidPoint says: " + msg);
}
else // Inside OMP blocks
{
Gambit::utils_warning().raise(LOCAL_INFO, "daFunk::RaisePipedInvalidPoint says: " + msg);
Gambit::piped_invalid_point.request("daFunk::RaisePipedInvalidPoint says: " + msg);
}
return 0;
}
private:
std::string msg; // Error message to throw when function is called
};
inline Funk raiseInvalidPoint(std::string msg) { return Funk(new RaiseInvalidPoint(msg)); }
#endif
//
// Prints message when called
//
class Bottle: public FunkBase
{
public:
Bottle(Funk f, std::string msg) : msg(msg)
{
functions = vec(f);
singularities = f->getSingl();
arguments = f->getArgs();
}
double value(const std::vector<double> & data, size_t bindID)
{
std::cout << "daFunk::Message says:\n" << msg << std::endl;
return functions[0]->value(data, bindID);
}
private:
std::string msg; // Message in a bottle. Ha!
};
//inline Funk print(std::string msg) { return Funk(new Bottle(msg)); }
inline Funk FunkBase::print(std::string msg)
{ return Funk(new Bottle(shared_from_this(), msg)); };
//
// GSL integration
//
class FunkIntegrate_gsl1d: public FunkBase, public gsl_function
{
public:
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f1, Funk f2)
{
setup(f0, arg, f1, f2);
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, Funk f)
{
setup(f0, arg, cnst(x), f);
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, double y)
{
setup(f0, arg, cnst(x), cnst(y));
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, double x)
{
setup(f0, arg, f, cnst(x));
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, Funk f)
{
setup(f0, arg, var(x), f);
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, std::string y)
{
setup(f0, arg, var(x), var(y));
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, std::string x)
{
setup(f0, arg, f, var(x));
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, double y)
{
setup(f0, arg, var(x), cnst(y));
}
FunkIntegrate_gsl1d(Funk f0, std::string arg, double y, std::string x)
{
setup(f0, arg, cnst(y), var(x));
}
void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap)
{
functions[1]->resolve(datamap, datalen, bindID, argmap); // Resolve boundary 0
functions[2]->resolve(datamap, datalen, bindID, argmap); // Resolve boundary 1
// Set indices
if(index.size() <= bindID)
{
index.resize(bindID+1);
// TODO: Throw error if problems are encountered
//std::cout << "FATAL ERROR: bind() attempts to resolve dependencies of FunkIntegrate_gsl1d" << std::endl;
//std::cout << "object inconsistently." << std::endl;
//std::cout << "Encountered while resolving " << arg << " and:" << std::endl;
//for (auto it = datamap.begin(); it != datamap.end(); it++)
//{
// std::cout << " " << it->first << std::endl;
//}
//exit(1);
}
if ( argmap.find(arg) == argmap.end() )
{
if(datamap.find(arg) == datamap.end())
{
index[bindID] = datalen;
argmap[arg] = datalen;
++datalen;
}
else
{
index[bindID] = datamap[arg];
}
}
else
{
index[bindID] = argmap[arg];
}
datamap[arg] = index[bindID];
functions[0]->resolve(datamap, datalen, bindID, argmap);
for ( auto it = my_singularities.begin(); it != my_singularities.end(); ++it )
{
it->first->resolve(datamap, datalen, bindID, argmap);
it->second->resolve(datamap, datalen, bindID, argmap);
}
}
~FunkIntegrate_gsl1d()
{
gsl_integration_workspace_free(gsl_workspace);
}
shared_ptr<FunkIntegrate_gsl1d> set_epsrel(double epsrel)
{ this->epsrel = epsrel; return static_pointer_cast<FunkIntegrate_gsl1d>(this->FunkIntegrate_gsl1d::shared_from_this()); }
shared_ptr<FunkIntegrate_gsl1d> set_epsabs(double epsabs)
{ this->epsabs = epsabs; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
shared_ptr<FunkIntegrate_gsl1d> set_limit(size_t limit)
{ this->limit = limit; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
shared_ptr<FunkIntegrate_gsl1d> set_singularity_factor(double f)
{ this->singl_factor = f; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
shared_ptr<FunkIntegrate_gsl1d> set_use_log_fallback(bool flag)
{ this->use_log_fallback = flag; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
double value(const std::vector<double> & data, size_t bindID)
{
double result;
#pragma omp critical(FunkIntegrate_gsl1d_integration)
{
local_data = data;
local_bindID = bindID;
double error;
function=&FunkIntegrate_gsl1d::invoke;
params=this;
double x0 = functions[1]->value(data, bindID);
double x1 = functions[2]->value(data, bindID);
gsl_set_error_handler_off();
int status = 0;
if ( my_singularities.size() == 0 )
{
status = gsl_integration_qags(this, x0, x1, epsabs, epsrel, limit, gsl_workspace, &result, &error);
}
else
{
double s = 0;
std::vector<double> ranges;
ranges.push_back(x0);
ranges.push_back(x1);
for ( auto it = my_singularities.begin(); it != my_singularities.end(); ++it )
{
double mean = it->first->value(data, bindID);
double sigma = it->second->value(data, bindID);
double z0 = mean - singl_factor*sigma;
double z1 = mean + singl_factor*sigma;
if ( z0 == z1 )
std::cout << "daFunk::FunkBase WARNING: Singularity width is beyond machine precision." << std::endl;
if ( z0 > x0 and z0 < x1 ) ranges.push_back(z0);
if ( z1 > x0 and z1 < x1 ) ranges.push_back(z1);
}
std::sort(ranges.begin(), ranges.end());
for ( auto it = ranges.begin(); it != ranges.end()-1; ++it )
{
status = gsl_integration_qags(this, *it, *(it+1), epsabs, epsrel, limit, gsl_workspace, &result, &error);
s += result;
if (status) break;
}
result = s;
}
if (status and this->use_log_fallback)
{
// The last resort: A cheap integration on log grid, linear interpolation
const double N = 300;
std::vector<double> Xgrid =
logspace(std::log10(x0), std::log10(x1), N);
double sum = 0, y0, y1, dx;
local_data[index[bindID]] = Xgrid[0];
y0 = functions[0]->value(local_data, bindID);
for (size_t i = 0; i<N-1; i++)
{
local_data[index[bindID]] = Xgrid[i+1];
y1 = functions[0]->value(local_data, bindID);
dx = Xgrid[i+1]-Xgrid[i];
sum += dx*(y0+y1)/2;
y0 = y1;
}
result = sum;
}
// TODO: Implement flags to optionally throw an error
if (status and not this->use_log_fallback)
{
std::cerr << "daFunk::FunkIntegrate_gsl1d WARNING: " << gsl_strerror(status) << std::endl;
std::cerr << "Attempt to integrate from " << x0 << " to " << x1 << std::endl;
std::cerr << "Attempt to integrate from " << x0 << " to " << x1 << std::endl;
std::cerr << "Details about the integrand:" << std::endl;
functions[0]->help();
// std::cout << "Dumping integrand:" << std::endl;
// for ( double x = x0; x <= x1; x = (x0>0) ? x*1.01 : x+(x1-x0)/1000)
// std::cerr << " " << x << " " << invoke(x, this) << std::endl;
std::cerr << "Returning zero." << std::endl;
result = 0.;
}
}
return result;
}
private:
void setup(Funk f0, std::string arg, Funk f1, Funk f2)
{
this->functions = vec(f0, f1, f2);
singularities = joinSingl(f1->getSingl(), f2->getSingl());
if ( f0->getSingl().find(arg) != f0->getSingl().end() )
my_singularities = f0->getSingl()[arg];
Singularities tmp_singl = f0->getSingl();
tmp_singl.erase(arg);
singularities = joinSingl(singularities, tmp_singl);
arguments = joinArgs(eraseArg(f0->getArgs(), arg), joinArgs(f1->getArgs(), f2->getArgs()));
gsl_workspace = gsl_integration_workspace_alloc (100000);
this->arg = arg;
limit = 100;
epsrel = 1e-2;
epsabs = 1e-2;
use_log_fallback = false;
singl_factor = 4;
}
// Static member function that invokes integrand
static double invoke(double x, void *params) {
FunkIntegrate_gsl1d * ptr = static_cast<FunkIntegrate_gsl1d*>(params);
ptr->local_data[ptr->index[ptr->local_bindID]] = x;
return ptr->functions[0]->value(ptr->local_data, ptr->local_bindID);
}
// Required for rewiring input parameters
std::vector<double> local_data;
size_t local_bindID;
std::vector<std::pair<Funk, Funk>> my_singularities;
// Integration range and function pointer
std::string arg;
// GSL workspace and parameters
gsl_integration_workspace * gsl_workspace;
size_t limit;
std::vector<size_t> index;
double epsrel;
double epsabs;
bool use_log_fallback;
double singl_factor;
};
// Standard behaviour
template <typename T1, typename T2>
inline shared_ptr<FunkIntegrate_gsl1d> getIntegrate_gsl1d(Funk fptr, std::string arg, T1 x1, T2 x2) { return shared_ptr<FunkIntegrate_gsl1d>(new FunkIntegrate_gsl1d(fptr, arg, x1, x2)); }
//
// Helper function for producing singularity-augmented x-grids
//
inline std::vector<double> augmentSingl(const std::vector<double> & xgrid, Funk f, int N = 100, double sigma = 3)
{
std::vector<double> result = xgrid;
double xmin = result.front();
double xmax = result.back();
if ( f->getNArgs() != 1 )
{
std::string msg = "augment_with_singularities(): takes only functions with one argument.\n";
msg+= " --- Actual arguments are: " + args_string(f->getArgs());
#ifdef GAMBIT_DIR
Gambit::utils_error().raise(LOCAL_INFO, msg);
#else
throw std::invalid_argument(msg);
#endif
}
std::string arg = f->getArgs()[0];
Singularities singlsMap = f->getSingl();
if ( singlsMap.find(arg) != singlsMap.end() )
{
auto singls = singlsMap.at(arg);
for (auto it = singls.begin(); it != singls.end(); it ++)
{
double position= it->first->bind()->eval();
double width = it->second->bind()->eval();
std::vector<double> singlgrid = linspace(std::max(position-width*sigma, xmin), std::min(position+width*sigma, xmax), N);
result.insert(result.end(), singlgrid.begin(), singlgrid.end());
}
}
// Sort the values only once
std::sort(result.begin(), result.end());
// Only use unique values, i.e. do not resolve the same singularity twice.
result.erase(std::unique(result.begin(), result.end()), result.end());
return result;
}
}
#endif // __DAFUNK_HPP__
Updated on 2023-06-26 at 21:36:55 +0000