|
GAMBIT
v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
|
Go to the documentation of this file. 53 #ifndef __DAFUNK_HPP__ 54 #define __DAFUNK_HPP__ 75 #include "boost/shared_ptr.hpp" 76 #include "boost/enable_shared_from_this.hpp" 84 #define DEF_FUNKTRAIT(C) \ 87 static daFunk::FunkPlain* ptr; \ 88 static void set(daFunk::FunkPlain* new_ptr) \ 94 daFunk::FunkPlain* C::ptr = NULL; 97 #include <gsl/gsl_integration.h> 98 #include <gsl/gsl_errno.h> 106 using boost::shared_ptr; 107 using boost::enable_shared_from_this; 108 using boost::dynamic_pointer_cast; 109 using boost::static_pointer_cast; 115 typedef shared_ptr<FunkBase> Funk; 118 typedef std::map<std::string, std::vector<std::pair<Funk, Funk>>> Singularities; 129 typedef std::pair<double(*)(double,void*), void*> PlainPtrs1; 130 typedef std::pair<double(*)(double,double,void*), void*> PlainPtrs2; 131 typedef std::pair<double(*)(double,double,double,void*), void*> PlainPtrs3; 132 typedef std::pair<double(*)(double,double,double,double,void*), void*> PlainPtrs4; 141 template < typename T> 142 inline std::vector<T> vec(std::vector<T> vector) 146 template < typename T, typename... Args> 147 inline std::vector<T> vec(std::vector<T> vector, T value, Args... args) 149 vector.push_back(value); 150 return vec<T>(vector, args...); 152 template < typename T, typename... Args> 153 inline std::vector<T> vec(T value, Args... args) 155 std::vector<T> vector; 156 vector.push_back(value); 157 vector = vec<T>(vector, args...); 160 template < typename T> 161 inline std::vector<T> vec() 163 std::vector<T> vector; 172 inline std::vector<double> linspace( double x0, double x1, unsigned int n) 174 std::vector<double> ret; 178 for ( unsigned int i = 0; i< n; i++) 180 ret.push_back(x0 + i * dx); 186 inline std::vector<double> logspace( double x0, double x1, unsigned int n) 188 std::vector<double> ret; 192 for ( unsigned int i = 0; i< n; i++) 194 ret.push_back( pow(10, x0 + i * dx)); 207 for ( auto it = args1.begin(); it!=args1.end(); it++ ) 209 if ( std::find(args2.begin(), args2.end(), *it) == args2.end() ) 212 if ( args1.size() != args2.size() ) m = false; 218 std::string msg = ""; 219 for ( auto it = args.begin(); it!=args.end(); it++ ) 222 if ( it != args.end() - 1) 228 inline ArgsType joinArgs(ArgsType args1, ArgsType args2) 230 args1.insert(args1.end(), args2.begin(), args2.end()); 231 std::set<std::string> argsset(args1.begin(), args1.end()); 232 args1.assign(argsset.begin(), argsset.end()); 238 auto it = std::find(args.begin(), args.end(), arg); 239 if (it!=args.end()) args.erase(it); 243 inline Singularities joinSingl(Singularities s1, Singularities s2) 245 for ( auto it = s2.begin(); it != s2.end(); ++it ) 247 if ( s1.find(it->first) == s1.end() ) 248 s1[it->first] = it->second; 251 for ( auto it2 = it->second.begin(); it2 != it->second.end(); ++it2 ) 253 if ( std::find(s1[it->first].begin(), s1[it->first].end(), *it2) == s1[it->first].end() ) 254 s1[it->first].push_back(*it2); 267 template < size_t... Is> 276 template < size_t MIN, size_t N, size_t... Is> 280 template < size_t MIN, size_t... Is> 287 template < size_t MIN, size_t N, size_t... Is> 298 class FunkBase: public enable_shared_from_this<FunkBase> 305 template < typename... Args> Funk set(std::string arg, Funk g, Args... args); 306 template < typename... Args> Funk set(std::string arg, double x, Args... args); 307 template < typename... Args> Funk set(std::string arg, std::string arg1, Args... args); 308 template < typename... Args> Funk set(); 311 template < typename... Args> shared_ptr<FunkBound> bind(Args... args); 314 const std::vector<std::string> & getArgs() { return this->arguments; }; 315 std::size_t getNArgs() { return this->arguments.size();}; 316 bool hasArg(std::string); 319 template < typename... Args> bool assert_args(Args... args); 322 virtual double value( const std::vector<double> &, size_t bindID) = 0; 330 virtual void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap); 336 Funk set_singularity(std::string arg, Funk pos, Funk width); 337 Funk set_singularity(std::string arg, double pos, Funk width); 338 Funk set_singularity(std::string arg, Funk pos, double width); 340 Funk set_singularity(std::string arg, double pos, double width); 351 template < typename... Args> shared_ptr<FunkIntegrate_gsl1d> gsl_integration(Args... args); 356 PlainPtrs3 plain(std::string, std::string, std::string); 357 PlainPtrs4 plain(std::string, std::string, std::string, std::string); 358 template < typename T> 360 template < typename T> 361 PlainPtr2 plain(std::string, std::string); 362 template < typename T> 363 PlainPtr3 plain(std::string, std::string, std::string); 364 template < typename T> 365 PlainPtr4 plain(std::string, std::string, std::string, std::string); 387 static bool dead = false; 395 FunkBound( Funk f, size_t datalen, size_t bindID) : f(f), datalen(datalen), bindID(bindID) {}; 397 double value(std::vector<double> & map, size_t bindID) {( void)bindID; ( void)map; return 0;}; 399 template < typename... Args> inline double eval(Args... argss) 401 auto data = vec<double>(argss...); 402 data.resize(datalen); 403 return f->value( data, bindID); 406 template < typename... Args> inline std::vector<double> vect(Args... argss) 408 std::vector<std::vector<double>> coll; 409 return this->vect2(coll, argss...); 413 template < typename... Args> 419 static size_t n_idx = 0; 423 #pragma omp critical (bindID_allocation) 432 bindID = free.back(); 441 #pragma omp critical (bindID_allocation) 442 free.push_back(bindID); 446 template < typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll) 449 std::vector<bool> vec_flag; 450 for ( auto it = coll.begin(); it != coll.end(); ++it ) 452 if ( it->size() == 1 ) 454 vec_flag.push_back( false); 457 vec_flag.push_back( true); 458 if ( size == 1 ) size = it->size(); 459 if ( size != it->size() ) 461 std::cerr << "daFunk::FunkBase WARNING: Inconsistent vector lengths." << std::endl; 462 return vec<double>(); 465 auto r = vec<double>(); 466 auto data = vec<double>(); 467 data.resize(datalen); 468 for ( size_t i = 0; i != size; ++i ) 470 for ( size_t j = 0; j != coll.size(); ++j ) 473 data[j] = coll[j][i]; 475 data[j] = coll[j][0]; 477 r.push_back( f->value(data, bindID)); 482 template < typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll, double x, Args... argss) 484 coll.push_back(vec<double>(x)); 485 return this->vect2(coll, argss...); 488 template < typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll, std::vector<double> v, Args... argss) 491 return this->vect2(coll, argss...); 518 FunkPlain( Funk fin, std::string arg1, std::string arg2, std::string arg3) : f(fin->bind(arg1, arg2, arg3)) {} 519 FunkPlain( Funk fin, std::string arg1, std::string arg2, std::string arg3, std::string arg4) : f(fin->bind(arg1, arg2, arg3, arg4)) {} 524 return funkPtrPtr-> f->eval(x1); 529 return funkPtrPtr-> f->eval(x1, x2); 534 return funkPtrPtr-> f->eval(x1, x2, x3); 536 static double plain4p( double x1, double x2, double x3, double x4, void* ptr) 539 return funkPtrPtr-> f->eval(x1, x2, x3, x4); 542 template < typename T> 546 return funkPtrPtr-> f->eval(x1); 548 template < typename T> 552 return funkPtrPtr-> f->eval(x1, x2); 554 template < typename T> 558 return funkPtrPtr-> f->eval(x1, x2, x3); 560 template < typename T> 564 return funkPtrPtr-> f->eval(x1, x2, x3, x4); 567 double value( const std::vector<double> & args, size_t bindID) 576 shared_ptr<FunkBound> f; 577 std::string arg1, arg2, arg3, arg4; 587 template < typename... Args> 588 FunkConst( double c, Args ...argss) : c(c) { arguments = vec<std::string>(argss...); } 591 double value( const std::vector<double> & data, size_t bindID) 601 template < typename... Args> 603 template < typename... Args> 605 template < typename... Args> 627 void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap) 629 functions[1]->resolve(datamap, datalen, bindID, argmap); 631 if(my_index.size() <= bindID) 633 my_index.resize(bindID+1); 646 if ( argmap.find(my_arg) == argmap.end() ) 648 if(datamap.find(my_arg) == datamap.end()) 650 my_index[bindID] = datalen; 651 argmap[my_arg] = datalen; 656 my_index[bindID] = datamap[my_arg]; 661 my_index[bindID] = argmap[my_arg]; 663 datamap[my_arg] = my_index[bindID]; 664 functions[0]->resolve(datamap, datalen, bindID, argmap); 667 double value( const std::vector<double> & data, size_t bindID) 669 std::vector<double> data2(data); 670 data2[my_index[bindID]] = functions[1]->value(data, bindID); 671 return functions[0]->value(data2, bindID); 680 functions = vec(f, g); 682 if ( tmp_singl.erase(arg) > 0 ) 683 std::cout << "daFunk::FunkBase WARNING: Loosing singularity information while setting " << arg << std::endl; 684 singularities = joinSingl(g->getSingl(), tmp_singl); 695 template < bool threadsafe, typename... funcargs> 699 template < typename... Args> 703 digest_input(argss...); 706 double value( const std::vector<double> & data, size_t bindID) 708 std::tuple<typename std::remove_reference<funcargs>::type...> my_input; 711 #pragma omp critical (FunkFunc_setInput) 713 for ( auto f = functions.begin(); f != functions.end(); ++ f, ++i) 715 *map[i] = (*f)->value(data, bindID); 725 #pragma omp critical (FunkFunc_externalFunctionCall) 733 template < size_t... Args> 736 return (*ptr)(std::get<Args>(my_input)...); 740 std::tuple<typename std::remove_reference<funcargs>::type...> input; 747 template< typename T, typename... Args> 750 const int i = sizeof...(funcargs) - sizeof...(argss) - 1; 751 std::get<i>(input) = x; 752 digest_input(argss...); 754 template< typename... Args> 757 const int i = sizeof...(funcargs) - sizeof...(argss) - 1; 758 map.push_back(&std::get<i>(input)); 759 arguments = joinArgs(arguments, f->getArgs()); 760 functions.push_back(f); 761 singularities = joinSingl(singularities, f->getSingl()); 762 digest_input(argss...); 767 template < typename... funcargs, typename... Args> 772 template < typename... funcargs, typename... Args> 778 template < bool threadsafe, typename O, typename... funcargs> 782 template < typename... Args> 783 FunkFuncM(O* obj, double (O::* f)(funcargs...), Args... argss) : obj(obj) 786 digest_input(argss...); 789 template < typename... Args> 790 FunkFuncM(shared_ptr<O> obj, double (O::* f)(funcargs...), Args... argss) : shared_obj(obj), obj(&*obj) 793 digest_input(argss...); 796 double value( const std::vector<double> & data, size_t bindID) 798 std::tuple<typename std::remove_reference<funcargs>::type...> my_input; 801 #pragma omp critical(FunkFuncM_value) 803 for ( auto f = functions.begin(); f != functions.end(); ++ f, ++i) 805 *map[i] = (*f)->value(data, bindID); 815 #pragma omp critical (FunkFuncM_objectFunctionCall) 823 template < size_t... Args> 826 return (*obj.*ptr)(std::get<Args>(my_input)...); 830 std::tuple<typename std::remove_reference<funcargs>::type...> input; 839 template< typename T, typename... Args> 842 const int i = sizeof...(funcargs) - sizeof...(argss) - 1; 843 std::get<i>(input) = x; 844 digest_input(argss...); 846 template< typename... Args> 849 const int i = sizeof...(funcargs) - sizeof...(argss) - 1; 850 map.push_back(&std::get<i>(input)); 851 arguments = joinArgs(arguments, f->getArgs()); 852 functions.push_back(f); 853 singularities = joinSingl(singularities, f->getSingl()); 854 digest_input(argss...); 860 template < typename O, typename... funcargs, typename... Args> 864 template < typename O, typename... funcargs, typename... Args> 870 template < typename O, typename... funcargs, typename... Args> 874 template < typename O, typename... funcargs, typename... Args> 886 FunkDelta(std::string arg, double pos, double width) : pos(pos), width(width) 888 arguments = vec(arg); 889 this->set_singularity( "v", pos, width); 890 singularities[arg].push_back(std::pair<Funk, Funk>( cnst(pos), cnst(width))); 893 double value( const std::vector<double> & data, size_t bindID) 895 double x = data[indices[bindID][0]]; 896 return exp(- pow(x-pos,2)/ pow(width,2)/2)/sqrt(2*M_PI)/width; 913 arguments = vec(arg); 916 double value( const std::vector<double> & data, size_t bindID) 918 return data[indices[bindID][0]]; 930 singularities[arg].push_back(std::pair<Funk, Funk>(pos, width)); 931 return shared_from_this(); 934 { return shared_from_this()->set_singularity(arg, cnst(pos), width); } 936 { return shared_from_this()->set_singularity(arg, cnst(pos), cnst(width)); } 938 { return shared_from_this()->set_singularity(arg, pos, cnst(width)); } 940 inline void FunkBase::resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap) 943 auto it1 = arguments.begin(); 944 std::vector<size_t> tmp_indices; 945 tmp_indices.resize(arguments.size()); 946 auto it2 = tmp_indices.begin(); 947 for (; it1 != arguments.end() && it2 != tmp_indices.end(); ++it1, ++it2 ) 951 *it2 = datamap.at(*it1); 953 catch (std::out_of_range & e) 955 std::string msg = "FunkBase::resolve() encountered internal problem when resolving " + *it1 + ".\n"; 956 msg+= " --- Actual arguments of object: " + args_string(arguments); 960 throw std::invalid_argument(msg); 966 if(indices.size() <= bindID) 968 indices.resize(bindID+1); 979 indices[bindID] = tmp_indices; 982 for ( auto it = functions.begin(); it != functions.end(); ++it) 984 (*it)->resolve(datamap, datalen, bindID, argmap); 991 std::vector<std::vector<std::string>> list = vec<std::vector<std::string>>( args...); 992 std::set<std::string> myargs(arguments.begin(), arguments.end()); 993 for ( auto it = list.begin(); it != list.end(); ++it ) 995 std::set<std::string> theirargs(it->begin(), it->end()); 996 if ( myargs == theirargs ) 1005 Funk f = shared_from_this(); 1006 if ( std::find(arguments.begin(), arguments.end(), arg) != arguments.end() ) 1012 std::cout << "daFunk::FunkBase WARNING: Ignoring \"" << arg << "\" = function." << std::endl; 1014 return f->set( args...); 1018 { return shared_from_this()->set(arg, var(arg1))->set( args...); } 1021 { return shared_from_this()->set(arg, cnst(x))->set( args...); } 1024 { return shared_from_this(); } 1026 template < typename... Args> inline shared_ptr<FunkBound> FunkBase::bind(Args... argss) 1030 std::map<std::string, size_t> datamap; 1032 auto bound_arguments = vec<std::string>(argss...); 1033 datalen = bound_arguments.size(); 1035 for ( auto it = bound_arguments.begin(); it != bound_arguments.end(); ++it, ++i ) 1039 std::map<std::string,size_t> argmap; 1040 if ( not args_match(arguments, bound_arguments) ) 1042 std::string msg = "FunkBase::bind() tries to resolve wrong arguments.\n"; 1043 msg+= " --- Arguments that are supposed to be bound: " + args_string(bound_arguments) + "\n"; 1044 msg+= " --- Actual arguments of object: " + args_string(arguments); 1048 throw std::invalid_argument(msg); 1051 this->resolve(datamap, datalen, bindID, argmap); 1052 return shared_ptr<FunkBound>( new FunkBound(shared_from_this(), datalen, bindID)); 1057 return ( std::find(arguments.begin(), arguments.end(), arg) != arguments.end() ); 1062 return ( this->arguments.size() != 0 ); 1072 std::cout << "Arguments:"; 1073 for ( auto it = arguments.begin(); it != arguments.end(); it++ ) 1075 std::cout << " \"" << *it << "\""; 1077 if ( arguments.size() == 0 ) 1078 std::cout << " none"; 1079 std::cout << std::endl; 1080 for ( auto it = singularities.begin() ; it != singularities.end(); ++it ) 1082 std::cout << "Singularities in " << it->first << ": " << it->second.size() << std::endl; 1084 return shared_from_this(); 1089 void* ptr = new FunkPlain(shared_from_this(), arg1); 1094 void* ptr = new FunkPlain(shared_from_this(), arg1, arg2); 1099 void* ptr = new FunkPlain(shared_from_this(), arg1, arg2, arg3); 1104 void* ptr = new FunkPlain(shared_from_this(), arg1, arg2, arg3, arg4); 1108 template < typename T> 1111 T::set( new FunkPlain(shared_from_this(), arg1)); 1112 return &FunkPlain::plain1<T>; 1114 template < typename T> 1117 T::set( new FunkPlain(shared_from_this(), arg1, arg2)); 1118 return &FunkPlain::plain2<T>; 1120 template < typename T> 1123 T::set( new FunkPlain(shared_from_this(), arg1, arg2, arg3)); 1124 return &FunkPlain::plain3<T>; 1126 template < typename T> 1129 T::set( new FunkPlain(shared_from_this(), arg1, arg2, arg3, arg4)); 1130 return &FunkPlain::plain4<T>; 1145 singularities = f->getSingl(); 1146 arguments = f->getArgs(); 1148 double value( const std::vector<double> & data, size_t bindID) 1150 return -(functions[0]->value(data, bindID)); 1156 #define MATH_OPERATION(OPERATION) \ 1157 class FunkMath_##OPERATION: public FunkBase \ 1160 FunkMath_##OPERATION(Funk f) \ 1162 functions = vec(f); \ 1163 arguments = f->getArgs(); \ 1164 singularities = f->getSingl(); \ 1166 double value(const std::vector<double> & data, size_t bindID) \ 1168 return OPERATION(functions[0]->value(data, bindID)); \ 1171 inline Funk OPERATION (Funk f) { return Funk(new FunkMath_##OPERATION(f)); } 1189 #undef MATH_OPERATION 1192 #define MATH_OPERATION(OPERATION, SYMBOL) \ 1193 class FunkMath_##OPERATION: public FunkBase \ 1196 FunkMath_##OPERATION(Funk f1, Funk f2) \ 1198 functions = vec(f1, f2); \ 1199 arguments = joinArgs(f1->getArgs(), f2->getArgs()); \ 1200 singularities = joinSingl(f1->getSingl(), f2->getSingl());\ 1202 FunkMath_##OPERATION(double x, Funk f2) \ 1204 auto f1 = cnst(x); \ 1205 functions = vec(f1, f2); \ 1206 arguments = joinArgs(f1->getArgs(), f2->getArgs()); \ 1207 singularities = joinSingl(f1->getSingl(), f2->getSingl());\ 1209 FunkMath_##OPERATION(Funk f1, double x) \ 1211 auto f2 = cnst(x); \ 1212 functions = vec(f1, f2); \ 1213 arguments = joinArgs(f1->getArgs(), f2->getArgs()); \ 1214 singularities = joinSingl(f1->getSingl(), f2->getSingl());\ 1216 double value(const std::vector<double> & data, size_t bindID) \ 1218 return functions[0]->value(data, bindID) SYMBOL functions[1]->value(data, bindID); \ 1221 inline Funk operator SYMBOL (Funk f1, Funk f2) { return Funk(new FunkMath_##OPERATION(f1, f2)); } \ 1222 inline Funk operator SYMBOL (double x, Funk f) { return Funk(new FunkMath_##OPERATION(x, f)); } \ 1223 inline Funk operator SYMBOL (Funk f, double x) { return Funk(new FunkMath_##OPERATION(f, x)); } 1228 #undef MATH_OPERATION 1231 #define MATH_OPERATION(OPERATION) \ 1232 class FunkMath_##OPERATION: public FunkBase \ 1235 FunkMath_##OPERATION(Funk f1, Funk f2) \ 1237 functions = vec(f1, f2); \ 1238 arguments = joinArgs(f1->getArgs(), f2->getArgs()); \ 1239 singularities = joinSingl(f1->getSingl(), f2->getSingl()); \ 1241 FunkMath_##OPERATION(double x, Funk f2) \ 1243 auto f1 = cnst(x); \ 1244 functions = vec(f1, f2); \ 1245 arguments = joinArgs(f1->getArgs(), f2->getArgs()); \ 1246 singularities = joinSingl(f1->getSingl(), f2->getSingl()); \ 1248 FunkMath_##OPERATION(Funk f1, double x) \ 1250 auto f2 = cnst(x); \ 1251 functions = vec(f1, f2); \ 1252 arguments = joinArgs(f1->getArgs(), f2->getArgs()); \ 1253 singularities = joinSingl(f1->getSingl(), f2->getSingl()); \ 1255 double value(const std::vector<double> & data, size_t bindID) \ 1257 return OPERATION(functions[0]->value(data, bindID), functions[1]->value(data, bindID)); \ 1260 inline Funk OPERATION (Funk f1, Funk f2) { return Funk(new FunkMath_##OPERATION(f1, f2)); } \ 1261 inline Funk OPERATION (double x, Funk f) { return Funk(new FunkMath_##OPERATION(x, f)); } \ 1262 inline Funk OPERATION (Funk f, double x) { return Funk(new FunkMath_##OPERATION(f, x)); } 1266 #undef MATH_OPERATION 1283 class FunkInterp : public FunkBase 1286 FunkInterp( Funk f, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin") 1288 setup(f, Xgrid, Ygrid, mode); 1290 FunkInterp(std::string arg, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin") 1292 setup( var(arg), Xgrid, Ygrid, mode); 1294 FunkInterp( double x, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin") 1299 double value( const std::vector<double> & data, size_t bindID) 1301 functions[0]->value(data, bindID); 1302 return (this->*ptr)(data[indices[bindID][0]]); 1306 void setup( Funk f, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode) 1310 singularities = f->getSingl(); 1311 arguments = f->getArgs(); 1312 this->Xgrid = Xgrid; 1313 this->Ygrid = Ygrid; 1314 if ( mode == "lin" ) this->ptr = &FunkInterp::linearInterp; 1315 else if ( mode == "log" ) this->ptr = &FunkInterp::logInterp; 1318 double logInterp( double x) 1321 int i = 0; int imax = Xgrid.size() - 1; 1322 if (x<Xgrid[0] or x>Xgrid[imax]) return 0; 1323 for (; i < imax; i++) { if (Xgrid[i] > x) break;}; 1324 double x0 = Xgrid[i-1]; 1325 double x1 = Xgrid[i]; 1326 double y0 = Ygrid[i-1]; 1327 double y1 = Ygrid[i]; 1328 return y0 * std::exp(std::log(y1/y0) * std::log(x/x0) / std::log(x1/x0)); 1331 double linearInterp( double x) 1334 int i = 0; int imax = Xgrid.size() - 1; 1335 if (x<Xgrid[0] or x>Xgrid[imax]) return 0; 1336 for (; i < imax; i++) { if (Xgrid[i] > x) break;}; 1337 double x0 = Xgrid[i-1]; 1338 double x1 = Xgrid[i]; 1339 double y0 = Ygrid[i-1]; 1340 double y1 = Ygrid[i]; 1341 return y0 + (x-x0)/(x1-x0)*(y1-y0); 1345 std::vector<double> Xgrid; 1346 std::vector<double> Ygrid; 1349 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)); } 1361 functions = vec(f, g, h); 1365 double value( const std::vector<double> & data, size_t bindID) 1367 if ( functions[0]->value(data,bindID) >= 0. ) 1368 return functions[1]->value(data,bindID); 1370 return functions[2]->value(data,bindID); 1389 double value( const std::vector<double> & data, size_t bindID) 1394 if ( omp_get_level() == 0 ) 1403 throw std::invalid_argument( "daFunk::ThrowError says: " + msg); 1414 class RaiseInvalidPoint: public FunkBase 1417 RaiseInvalidPoint(std::string msg) : msg(msg) 1420 double value( const std::vector<double> & data, size_t bindID) 1424 if ( omp_get_level() == 0 ) 1440 inline Funk raiseInvalidPoint(std::string msg) { return Funk( new RaiseInvalidPoint(msg)); } 1454 singularities = f->getSingl(); 1455 arguments = f->getArgs(); 1457 double value( const std::vector<double> & data, size_t bindID) 1459 std::cout << "daFunk::Message says:\n" << msg << std::endl; 1460 return functions[0]->value(data, bindID); 1468 { return Funk( new Bottle(shared_from_this(), msg)); }; 1480 setup(f0, arg, f1, f2); 1515 void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap) 1517 functions[1]->resolve(datamap, datalen, bindID, argmap); 1518 functions[2]->resolve(datamap, datalen, bindID, argmap); 1521 if( index.size() <= bindID) 1523 index.resize(bindID+1); 1534 if ( argmap.find(arg) == argmap.end() ) 1536 if(datamap.find(arg) == datamap.end()) 1538 index[bindID] = datalen; 1539 argmap[arg] = datalen; 1544 index[bindID] = datamap[arg]; 1549 index[bindID] = argmap[arg]; 1551 datamap[arg] = index[bindID]; 1552 functions[0]->resolve(datamap, datalen, bindID, argmap); 1553 for ( auto it = my_singularities.begin(); it != my_singularities.end(); ++it ) 1555 it->first->resolve(datamap, datalen, bindID, argmap); 1556 it->second->resolve(datamap, datalen, bindID, argmap); 1562 gsl_integration_workspace_free(gsl_workspace); 1566 { this->epsrel = epsrel; return static_pointer_cast< FunkIntegrate_gsl1d>(this->FunkIntegrate_gsl1d::shared_from_this()); } 1568 { this->epsabs = epsabs; return static_pointer_cast< FunkIntegrate_gsl1d>(this->shared_from_this()); } 1570 { this->limit = limit; return static_pointer_cast< FunkIntegrate_gsl1d>(this->shared_from_this()); } 1572 { this->singl_factor = f; return static_pointer_cast< FunkIntegrate_gsl1d>(this->shared_from_this()); } 1574 { this->use_log_fallback = flag; return static_pointer_cast< FunkIntegrate_gsl1d>(this->shared_from_this()); } 1576 double value( const std::vector<double> & data, size_t bindID) 1579 #pragma omp critical(FunkIntegrate_gsl1d_integration) 1582 local_bindID = bindID; 1586 double x0 = functions[1]->value(data, bindID); 1587 double x1 = functions[2]->value(data, bindID); 1588 gsl_set_error_handler_off(); 1590 if ( my_singularities.size() == 0 ) 1592 status = gsl_integration_qags( this, x0, x1, epsabs, epsrel, limit, gsl_workspace, &result, &error); 1597 std::vector<double> ranges; 1598 ranges.push_back(x0); 1599 ranges.push_back(x1); 1600 for ( auto it = my_singularities.begin(); it != my_singularities.end(); ++it ) 1602 double mean = it->first->value(data, bindID); 1603 double sigma = it->second->value(data, bindID); 1604 double z0 = mean - singl_factor* sigma; 1605 double z1 = mean + singl_factor* sigma; 1607 std::cout << "daFunk::FunkBase WARNING: Singularity width is beyond machine precision." << std::endl; 1608 if ( z0 > x0 and z0 < x1 ) ranges.push_back(z0); 1609 if ( z1 > x0 and z1 < x1 ) ranges.push_back(z1); 1611 std::sort(ranges.begin(), ranges.end()); 1612 for ( auto it = ranges.begin(); it != ranges.end()-1; ++it ) 1614 status = gsl_integration_qags( this, *it, *(it+1), epsabs, epsrel, limit, gsl_workspace, &result, &error); 1620 if (status and this->use_log_fallback) 1623 const double N = 300; 1624 std::vector<double> Xgrid = 1625 logspace(std::log10(x0), std::log10(x1), N); 1626 double sum = 0, y0, y1, dx; 1627 local_data[ index[bindID]] = Xgrid[0]; 1628 y0 = functions[0]->value(local_data, bindID); 1629 for ( size_t i = 0; i<N-1; i++) 1631 local_data[index[bindID]] = Xgrid[i+1]; 1632 y1 = functions[0]->value(local_data, bindID); 1633 dx = Xgrid[i+1]-Xgrid[i]; 1634 sum += dx*(y0+y1)/2; 1640 if (status and not this->use_log_fallback) 1642 std::cerr << "daFunk::FunkIntegrate_gsl1d WARNING: " << gsl_strerror(status) << std::endl; 1643 std::cerr << "Attempt to integrate from " << x0 << " to " << x1 << std::endl; 1644 std::cerr << "Attempt to integrate from " << x0 << " to " << x1 << std::endl; 1645 std::cerr << "Details about the integrand:" << std::endl; 1646 functions[0]->help(); 1650 std::cerr << "Returning zero." << std::endl; 1660 this->functions = vec(f0, f1, f2); 1662 singularities = joinSingl(f1->getSingl(), f2->getSingl()); 1663 if ( f0->getSingl().find(arg) != f0->getSingl().end() ) 1664 my_singularities = f0->getSingl()[arg]; 1666 tmp_singl.erase(arg); 1667 singularities = joinSingl(singularities, tmp_singl); 1670 gsl_workspace = gsl_integration_workspace_alloc (100000); 1676 use_log_fallback = false; 1681 static double invoke( double x, void *params) { 1707 template < typename T1, typename T2> 1716 std::vector<double> result = xgrid; 1717 double xmin = result.front(); 1718 double xmax = result.back(); 1720 if ( f->getNArgs() != 1 ) 1722 std::string msg = "augment_with_singularities(): takes only functions with one argument.\n"; 1723 msg+= " --- Actual arguments are: " + args_string(f->getArgs()); 1727 throw std::invalid_argument(msg); 1731 std::string arg = f->getArgs()[0]; 1734 if ( singlsMap.find(arg) != singlsMap.end() ) 1736 auto singls = singlsMap.at(arg); 1737 for ( auto it = singls.begin(); it != singls.end(); it ++) 1739 double position= it->first->bind()->eval(); 1740 double width = it->second->bind()->eval(); 1741 std::vector<double> singlgrid = linspace(std::max(position-width* sigma, xmin), std::min(position+width*sigma, xmax), N); 1742 result.insert(result.end(), singlgrid.begin(), singlgrid.end()); 1744 std::sort(result.begin(), result.end()); 1752 #endif // __DAFUNK_HPP__ FunkFuncM(shared_ptr< O > obj, double(O::*f)(funcargs...), Args... argss)
shared_ptr< FunkIntegrate_gsl1d > set_use_log_fallback(bool flag)
void print(MixMatrix A) Helper function to print a matrix.
std::vector< double > vect2(std::vector< std::vector< double >> &coll, std::vector< double > v, Args... argss)
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
Piped_exceptions piped_errors Global instance of Piped_exceptions class for errors.
FunkPlain(Funk fin, std::string arg1, std::string arg2, std::string arg3, std::string arg4)
double(* PlainPtr3)(double &, double &, double &)
double(* fptr)(int &) Pointer to a function that takes an integer by reference and returns a double.
Funk throwError(std::string msg)
std::pair< double(*)(double, void *), void * > PlainPtrs1
FunkIfElse(Funk f, Funk g, Funk h)
std::vector< double > vect(Args... argss)
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) MATH_OPERATION(Sum
FunkPlain(Funk fin, std::string arg1, std::string arg2, std::string arg3)
General small utility macros.
double ppp(index_list< Args... >, std::tuple< typename std::remove_reference< funcargs >::type... > &my_input)
Bottle(Funk f, std::string msg)
EXPORT_SYMBOLS error & utils_error() Utility errors.
void resolve(std::map< std::string, size_t > datamap, size_t &datalen, size_t bindID, std::map< std::string, size_t > &argmap)
FunkIntegrate_gsl1d(Funk f0, std::string arg, double y, std::string x)
virtual void resolve(std::map< std::string, size_t > datamap, size_t &datalen, size_t bindID, std::map< std::string, size_t > &argmap)
shared_ptr< FunkIntegrate_gsl1d > getIntegrate_gsl1d(Funk fptr, std::string arg, T1 x1, T2 x2)
FunkDerived(Funk f, std::string arg, double x)
static double plain3p(double x1, double x2, double x3, void *ptr)
FunkConst(double c, Args ...argss)
std::string args_string(ArgsType args)
Funk print(std::string arg)
FunkFuncM(O *obj, double(O::*f)(funcargs...), Args... argss)
const std::vector< std::string > & getArgs()
Funk delta(std::string arg, double pos, double width)
ArgsType eraseArg(ArgsType args, std::string arg)
bool args_match(ArgsType args1, ArgsType args2)
Funk funcM(O *obj, double(O::*f)(funcargs...), Args... args)
double value(const std::vector< double > &data, size_t bindID)
double value(const std::vector< double > &data, size_t bindID)
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, double x)
void digest_input(T x, Args... argss)
void request(std::string origin, std::string message) Request an exception.
FunkBound(Funk f, size_t datalen, size_t bindID)
std::pair< double(*)(double, double, double, double, void *), void * > PlainPtrs4
void digest_input(Funk f, Args... argss)
std::vector< double > linspace(double x0, double x1, unsigned int n)
FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, Funk f)
std::vector< double > local_data
FunkFunc(double(*f)(funcargs...), Args... argss)
gsl_integration_workspace * gsl_workspace
Funk var(std::string arg)
static double plain2p(double x1, double x2, void *ptr)
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, Funk f)
shared_ptr< FunkIntegrate_gsl1d > set_limit(size_t limit)
Piped_invalid_point piped_invalid_point Global instance of piped invalid point class.
shared_ptr< FunkIntegrate_gsl1d > set_epsrel(double epsrel)
std::tuple< typename std::remove_reference< funcargs >::type... > input
Funk set_singularity(std::string arg, Funk pos, Funk width)
double value(const std::vector< double > &data, size_t bindID)
Funk funcM_fromThreadsafe(O *obj, double(O::*f)(funcargs...), Args... args)
double value(const std::vector< double > &args, size_t bindID)
void digest_input(Funk f, Args... argss)
std::map< std::string, std::vector< std::pair< Funk, Funk > > > Singularities
std::vector< double > vect2(std::vector< std::vector< double >> &coll)
std::tuple< typename std::remove_reference< funcargs >::type... > input
Funk func_fromThreadsafe(double(*f)(funcargs...), Args... args)
std::vector< double * > map
static double plain2(double &x1, double &x2)
void digest_input(T x, Args... argss)
shared_ptr< O > shared_obj
double value(const std::vector< double > &data, size_t bindID)
shared_ptr< FunkBound > bind(Args... args)
ThrowError(std::string msg)
std::vector< Funk > functions
double eval(Args... argss)
virtual void raise(const std::string &) Raise the exception, i.e. throw it. Exact override of base method.
MATH_OPERATION(Dif,-) MATH_OPERATION(pow) MATH_OPERATION(fmin) MATH_OPERATION(fmax) class FunkInterp shared_ptr< FunkInterp > interp(T f, std::vector< double > x, std::vector< double > y)
shared_ptr< FunkBound > BoundFunk
void resolve(std::map< std::string, size_t > datamap, size_t &datalen, size_t bindID, std::map< std::string, size_t > &argmap)
shared_ptr< FunkIntegrate_gsl1d > set_singularity_factor(double f)
ArgsType joinArgs(ArgsType args1, ArgsType args2)
EXPORT_SYMBOLS warning & utils_warning() Utility warnings.
PlainPtrs1 plain(std::string)
static double plain4p(double x1, double x2, double x3, double x4, void *ptr)
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, std::string x)
std::pair< double(*)(double, double, double, void *), void * > PlainPtrs3
void setup(Funk f0, std::string arg, Funk f1, Funk f2)
double value(const std::vector< double > &data, size_t bindID)
FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, double y)
double value(const std::vector< double > &data, size_t bindID)
std::vector< double > logspace(double x0, double x1, unsigned int n)
std::vector< std::pair< Funk, Funk > > my_singularities
Funk func(double(*f)(funcargs...), Args... args)
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f1, Funk f2)
shared_ptr< FunkBound > f
std::vector< double * > map
std::vector< size_t > index
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, double y)
FunkDelta(std::string arg, double pos, double width)
void setup(Funk f, std::string arg, Funk g)
bool assert_args(Args... args)
double(* PlainPtr4)(double &, double &, double &, double &)
std::vector< std::vector< size_t > > indices
Singularities singularities
static void bindID_manager(size_t &bindID, bool bind)
std::vector< double > vect2(std::vector< std::vector< double >> &coll, double x, Args... argss)
FunkDerived(Funk f, std::string arg, Funk g)
double value(const std::vector< double > &data, size_t bindID)
double(* PlainPtr1)(double &)
Funk cnst(double x, Args... argss)
double pow(const double &a) Outputs a^i.
Exception objects required for standalone compilation.
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, std::string y)
void request(std::string message) Request an exception.
std::vector< double > augmentSingl(const std::vector< double > &xgrid, Funk f, int N=100, double sigma=3)
FunkPlain(Funk fin, std::string arg1, std::string arg2)
double value(std::vector< double > &map, size_t bindID)
double value(const std::vector< double > &data, size_t bindID)
invalid_point_exception & invalid_point() Invalid point exceptions.
static double plain1p(double x1, void *ptr)
std::pair< double(*)(double, double, void *), void * > PlainPtrs2
static double plain3(double &x1, double &x2, double &x3)
std::vector< T > vec(std::vector< T > vector)
double(* PlainPtr2)(double &, double &)
Funk ifelse(Funk f, Funk g, Funk h)
static double plain4(double &x1, double &x2, double &x3, double &x4)
shared_ptr< FunkIntegrate_gsl1d > set_epsabs(double epsabs)
shared_ptr< FunkBase > Funk
double value(const std::vector< double > &data, size_t bindID)
std::vector< size_t > my_index
std::vector< std::string > ArgsType
double value(const std::vector< double > &data, size_t bindID)
FunkPlain(Funk fin, std::string arg1)
shared_ptr< FunkIntegrate_gsl1d > gsl_integration(Args... args)
static double plain1(double &x1)
double value(const std::vector< double > &data, size_t bindID)
Singularities joinSingl(Singularities s1, Singularities s2)
static double invoke(double x, void *params)
double ppp(index_list< Args... >, std::tuple< typename std::remove_reference< funcargs >::type... > &my_input)
|