gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
daFunk.hpp
Go to the documentation of this file.
1 /*
2  * _______
3  * \ ___ `'. _..._ .
4  * ' |--.\ \ _.._ .' '. .'|
5  * | | \ ' .' .._| . .-. . .' |
6  * | | | ' __ | ' | ' ' |< |
7  * | | | | .:--.'. __| |__ _ _ | | | | | | ____
8  * | | ' .'/ | \ ||__ __|| ' / | | | | | | | \ .'
9  * | |___.' /' `" __ | | | | .' | .' | | | | | | |/ .
10  * /_______.'/ .'.''| | | | / | / | | | | | | /\ \
11  * \_______|/ / / | |_ | | | `'. | | | | | | | \ \
12  * \ \._,\ '/ | | ' .'| '/| | | | ' \ \ \
13  * `--' `" |_| `-' `--' '--' '--''------' '---'
14  *
15  * daFunk - dynamisch allokierbare Funktionen
16  *
17  * v0.1 Dec 2014
18  * v0.2 Mar 2015 - Completely rewritten internal structure
19  * v0.3 May 2016 - Extensions
20  *
21  * Christoph Weniger, created Dec 2014, edited until May 2016
22  * <c.weniger@uva.nl>
23  *
24  * with contributions related to thread-safety from
25  * Lars A. Dal, Apr, Jun 2015
26  * <l.a.dal@fys.uio.no>
27  *
28  *
29  *
30  * The MIT License (MIT)
31  *
32  * Copyright (c) 2016 Christoph Weniger
33  *
34  * Permission is hereby granted, free of charge, to any person obtaining a copy of
35  * this software and associated documentation files (the "Software"), to deal in
36  * the Software without restriction, including without limitation the rights to
37  * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
38  * of the Software, and to permit persons to whom the Software is furnished to do
39  * so, subject to the following conditions:
40  *
41  * The above copyright notice and this permission notice shall be included in all
42  * copies or substantial portions of the Software.
43  *
44  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
45  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
46  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
47  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
48  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
49  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
50  * SOFTWARE.
51  */
52 
53 #ifndef __DAFUNK_HPP__
54 #define __DAFUNK_HPP__
55 
56 #include <iostream>
57 #include <vector>
58 #include <algorithm>
59 #include <map>
60 #include <set>
61 #include <cmath>
62 #include <stdexcept>
63 
64 //#define NDEBUG
65 #include <assert.h>
66 
67 //#define GAMBIT_DIR
68 
69 #ifdef GAMBIT_DIR
70 #include <omp.h>
73 #endif
74 
75 #include "boost/shared_ptr.hpp"
76 #include "boost/enable_shared_from_this.hpp"
77 
78 //#include <memory>
79 //using std::shared_ptr;
80 //using std::enable_shared_from_this;
81 
82 namespace daFunk {class FunkPlain;}
83 
84 #define DEF_FUNKTRAIT(C) \
85 class C { \
86  public: \
87  static daFunk::FunkPlain* ptr; \
88  static void set(daFunk::FunkPlain* new_ptr) \
89  { \
90  delete ptr; \
91  ptr = new_ptr; \
92  } \
93 }; \
94 daFunk::FunkPlain* C::ptr = NULL;
95 
96 // Extensions
97 #include <gsl/gsl_integration.h>
98 #include <gsl/gsl_errno.h>
99 
100 namespace daFunk
101 {
102  //
103  // Type declarations etc.
104  //
105 
106  using boost::shared_ptr;
107  using boost::enable_shared_from_this;
108  using boost::dynamic_pointer_cast;
109  using boost::static_pointer_cast;
110 
111  class FunkBase;
112  class FunkBound;
114 
115  typedef shared_ptr<FunkBase> Funk;
116  typedef shared_ptr<FunkBound> BoundFunk;
117  typedef std::vector<std::string> ArgsType;
118  typedef std::map<std::string, std::vector<std::pair<Funk, Funk>>> Singularities;
119 
120  //template <typename... Args>
121  //using PlainPtr = double(*)(Args&...);
122  typedef double(*PlainPtr1)(double&);
123  typedef double(*PlainPtr2)(double&,double&);
124  typedef double(*PlainPtr3)(double&,double&,double&);
125  typedef double(*PlainPtr4)(double&,double&,double&,double&);
126 
127  //template <typename... Args>
128  //using PlainPtrs = std::pair<double(*)(Args...,void*), void*>;
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;
133 
134 
135  //
136  // Vector initialization from argument list
137  //
138  // Usage: std::vector<T> v = vec<T>(v1, v2, v3, ...);
139  //
140 
141  template <typename T>
142  inline std::vector<T> vec(std::vector<T> vector)
143  {
144  return vector;
145  }
146  template <typename T, typename... Args>
147  inline std::vector<T> vec(std::vector<T> vector, T value, Args... args)
148  {
149  vector.push_back(value);
150  return vec<T>(vector, args...);
151  }
152  template <typename T, typename... Args>
153  inline std::vector<T> vec(T value, Args... args)
154  {
155  std::vector<T> vector;
156  vector.push_back(value);
157  vector = vec<T>(vector, args...);
158  return vector;
159  }
160  template <typename T>
161  inline std::vector<T> vec()
162  {
163  std::vector<T> vector;
164  return vector;
165  }
166 
167  //
168  // Math helper functions
169  //
170 
171  // Generate linear 1-dim grid
172  inline std::vector<double> linspace(double x0, double x1, unsigned int n)
173  {
174  std::vector<double> ret;
175  double dx = 0;
176  if (n > 1)
177  dx = (x1-x0)/(n-1);
178  for (unsigned int i = 0; i<n; i++)
179  {
180  ret.push_back(x0 + i * dx);
181  }
182  return ret;
183  }
184 
185  // Generate logarithmic 1-dim grid
186  inline std::vector<double> logspace(double x0, double x1, unsigned int n)
187  {
188  std::vector<double> ret;
189  double dx = 0;
190  if (n > 1)
191  dx = (x1-x0)/(n-1);
192  for (unsigned int i = 0; i<n; i++)
193  {
194  ret.push_back(pow(10, x0 + i * dx));
195  }
196  return ret;
197  }
198 
199 
200  //
201  // Helper functions for internal calculations
202  //
203 
204  inline bool args_match(ArgsType args1, ArgsType args2)
205  {
206  bool m = true;
207  for ( auto it = args1.begin(); it!=args1.end(); it++ )
208  {
209  if ( std::find(args2.begin(), args2.end(), *it) == args2.end() )
210  m = false;
211  }
212  if ( args1.size() != args2.size() ) m = false;
213  return m;
214  }
215 
216  inline std::string args_string(ArgsType args)
217  {
218  std::string msg = "";
219  for ( auto it = args.begin(); it!=args.end(); it++ )
220  {
221  msg += *it;
222  if ( it != args.end() - 1)
223  msg += ", ";
224  }
225  return msg;
226  }
227 
228  inline ArgsType joinArgs(ArgsType args1, ArgsType args2)
229  {
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());
233  return args1;
234  }
235 
236  inline ArgsType eraseArg(ArgsType args, std::string arg)
237  {
238  auto it = std::find(args.begin(), args.end(), arg);
239  if (it!=args.end()) args.erase(it);
240  return args;
241  }
242 
243  inline Singularities joinSingl(Singularities s1, Singularities s2)
244  {
245  for ( auto it = s2.begin(); it != s2.end(); ++it )
246  {
247  if ( s1.find(it->first) == s1.end() )
248  s1[it->first] = it->second;
249  else
250  {
251  for ( auto it2 = it->second.begin(); it2 != it->second.end(); ++it2 )
252  {
253  if ( std::find(s1[it->first].begin(), s1[it->first].end(), *it2) == s1[it->first].end() )
254  s1[it->first].push_back(*it2);
255  }
256  }
257  }
258  return s1;
259  }
260 
261 
262  //
263  // Index lists (taken from stackoverflow)
264  //
265 
266  // The structure that encapsulates index lists
267  template <size_t... Is>
268  struct index_list
269  {
270  };
271 
272  // Collects internal details for generating index ranges [MIN, MAX)
273  namespace detail
274  {
275  // Declare primary template for index range builder
276  template <size_t MIN, size_t N, size_t... Is>
278 
279  // Base step
280  template <size_t MIN, size_t... Is>
281  struct range_builder<MIN, MIN, Is...>
282  {
283  typedef index_list<Is...> type;
284  };
285 
286  // Induction step
287  template <size_t MIN, size_t N, size_t... Is>
288  struct range_builder : public range_builder<MIN, N - 1, N - 1, Is...>
289  {
290  };
291  }
292 
293 
294  //
295  // **** The central virtual base class ****
296  //
297 
298  class FunkBase: public enable_shared_from_this<FunkBase>
299  {
300  public:
301  FunkBase() {};
302  virtual ~FunkBase() {};
303 
304  // Standard setting handles
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();
309 
310  // Standard binding handles
311  template <typename... Args> shared_ptr<FunkBound> bind(Args... args);
312 
313  // Convenience functions
314  const std::vector<std::string> & getArgs() { return this->arguments; };
315  std::size_t getNArgs() {return this->arguments.size();};
316  bool hasArg(std::string);
317  bool hasArgs();
318  Funk help();
319  template <typename... Args> bool assert_args(Args... args);
320 
321  // Return value & standard resolve
322  virtual double value(const std::vector<double> &, size_t bindID) = 0;
323 
324  // datamap maps the required function arguments onto the specific
325  // entries in the double-valued data array generated by eval().
326  // datalen acts like a pointer on the last entry of that array, and
327  // can be increased if more workspace is required. bindID ensures
328  // that resolution for various different binds can happen in
329  // parallel with the same Funk objects.
330  virtual void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap);
331 
332 
333  // Singularities handling
334  Singularities getSingl() { return singularities; }
335  private: // Note: This works currently only for constant functions
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);
339  public:
340  Funk set_singularity(std::string arg, double pos, double width);
341 
342  // Print message
343  Funk print(std::string arg);
344 
345 
346  //
347  // "External" functions
348  //
349 
350  // Integration
351  template <typename... Args> shared_ptr<FunkIntegrate_gsl1d> gsl_integration(Args... args);
352 
353  // Plain function generators
354  PlainPtrs1 plain(std::string);
355  PlainPtrs2 plain(std::string, std::string);
356  PlainPtrs3 plain(std::string, std::string, std::string);
357  PlainPtrs4 plain(std::string, std::string, std::string, std::string);
358  template <typename T>
359  PlainPtr1 plain(std::string);
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);
366 
367  protected:
368  std::vector<Funk> functions; // Dependent functions
369  ArgsType arguments; // Argument names
370  std::vector<std::vector<size_t> > indices; // Indices for data object
371  size_t datalen;
373  };
374 
375  // A vector class with global knowledge about its health status.
376  // (BoundFunk objects are occasionally destructed *after* livingVector has
377  // been destructed, causing segfaults if not catched properly.)
378  class livingVector: public std::vector<size_t>
379  {
380  public:
382  {
383  livingVector::is_dead() = true;
384  }
385  static bool & is_dead()
386  {
387  static bool dead = false;
388  return dead;
389  }
390  };
391 
392  class FunkBound
393  {
394  public:
395  FunkBound(Funk f, size_t datalen, size_t bindID) : f(f), datalen(datalen), bindID(bindID) {};
396  ~FunkBound() {bindID_manager(bindID,false);};
397  double value(std::vector<double> & map, size_t bindID) {(void)bindID; (void)map; return 0;};
398 
399  template <typename... Args> inline double eval(Args... argss)
400  {
401  auto data = vec<double>(argss...);
402  data.resize(datalen);
403  return f->value(data, bindID);
404  }
405 
406  template <typename... Args> inline std::vector<double> vect(Args... argss)
407  {
408  std::vector<std::vector<double>> coll;
409  return this->vect2(coll, argss...);
410  }
411 
412  private:
413  template <typename... Args>
414  friend shared_ptr<FunkBound> FunkBase::bind(Args... argss);
415  // Function for generating unique bindIDs.
416  // IDs are sequential, starting from 0.
417  static void bindID_manager(size_t &bindID, bool bind)
418  {
419  static size_t n_idx = 0;
420  static livingVector free;
421  if(bind)
422  {
423  #pragma omp critical (bindID_allocation)
424  {
425  if(free.empty())
426  {
427  // Increment after assignment (n_idx++ instead of ++n_idx)
428  bindID = n_idx++;
429  }
430  else
431  {
432  bindID = free.back();
433  free.pop_back();
434  }
435  }
436  }
437  else
438  {
439  if ( not livingVector::is_dead() )
440  {
441  #pragma omp critical (bindID_allocation)
442  free.push_back(bindID);
443  }
444  }
445  }
446  template <typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll)
447  {
448  size_t size = 1;
449  std::vector<bool> vec_flag;
450  for ( auto it = coll.begin(); it != coll.end(); ++it )
451  {
452  if ( it->size() == 1 )
453  {
454  vec_flag.push_back(false);
455  continue;
456  }
457  vec_flag.push_back(true);
458  if ( size == 1 ) size = it->size();
459  if ( size != it->size() )
460  {
461  std::cerr << "daFunk::FunkBase WARNING: Inconsistent vector lengths." << std::endl;
462  return vec<double>();
463  }
464  }
465  auto r = vec<double>();
466  auto data = vec<double>();
467  data.resize(datalen);
468  for ( size_t i = 0; i != size; ++i )
469  {
470  for ( size_t j = 0; j != coll.size(); ++j )
471  {
472  if ( vec_flag[j] )
473  data[j] = coll[j][i];
474  else
475  data[j] = coll[j][0];
476  }
477  r.push_back(f->value(data, bindID));
478  }
479  return r;
480  }
481 
482  template <typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll, double x, Args... argss)
483  {
484  coll.push_back(vec<double>(x));
485  return this->vect2(coll, argss...);
486  }
487 
488  template <typename... Args> inline std::vector<double> vect2(std::vector<std::vector<double>> & coll, std::vector<double> v, Args... argss)
489  {
490  coll.push_back(v);
491  return this->vect2(coll, argss...);
492  }
493 
494  Funk f; // bound function
495 
496  // datalen is the length of the double-valued data array that is
497  // needed as workspace for function evaluation, and that is created
498  // on the heap for each eval separately to ensure thread-safety.
499  size_t datalen;
500 
501  // bindID has the purpose of allowing bound functions (instances of
502  // FunkBase and daughter classes) to be bound by various binding
503  // functions simultaneously.
504  size_t bindID;
505  };
506 
507 
508  //
509  // Derived class with (templated) static member functions as plain function
510  // prototypes.
511  //
512 
513  class FunkPlain: public FunkBase
514  {
515  public:
516  FunkPlain(Funk fin, std::string arg1) : f(fin->bind(arg1)) {}
517  FunkPlain(Funk fin, std::string arg1, std::string arg2) : f(fin->bind(arg1, arg2)) {}
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)) {}
520 
521  static double plain1p(double x1, void* ptr)
522  {
523  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
524  return funkPtrPtr->f->eval(x1);
525  }
526  static double plain2p(double x1, double x2, void* ptr)
527  {
528  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
529  return funkPtrPtr->f->eval(x1, x2);
530  }
531  static double plain3p(double x1, double x2, double x3, void* ptr)
532  {
533  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
534  return funkPtrPtr->f->eval(x1, x2, x3);
535  }
536  static double plain4p(double x1, double x2, double x3, double x4, void* ptr)
537  {
538  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(ptr);
539  return funkPtrPtr->f->eval(x1, x2, x3, x4);
540  }
541 
542  template <typename T>
543  static double plain1(double& x1)
544  {
545  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
546  return funkPtrPtr->f->eval(x1);
547  }
548  template <typename T>
549  static double plain2(double& x1, double& x2)
550  {
551  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
552  return funkPtrPtr->f->eval(x1, x2);
553  }
554  template <typename T>
555  static double plain3(double& x1, double& x2, double& x3)
556  {
557  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
558  return funkPtrPtr->f->eval(x1, x2, x3);
559  }
560  template <typename T>
561  static double plain4(double& x1, double& x2, double& x3, double& x4)
562  {
563  FunkPlain * funkPtrPtr = static_cast<FunkPlain*>(T::ptr);
564  return funkPtrPtr->f->eval(x1, x2, x3, x4);
565  }
566 
567  double value(const std::vector<double> & args, size_t bindID)
568  {
569  (void)args;
570  (void)bindID;
571  assert ( 0 == 1 ); // This function should never be called
572  return 0;
573  }
574 
575  private:
576  shared_ptr<FunkBound> f; // bound function
577  std::string arg1, arg2, arg3, arg4;
578  };
579 
580  //
581  // Derived class that implements constant
582  //
583 
584  class FunkConst: public FunkBase
585  {
586  public:
587  template <typename... Args>
588  FunkConst(double c, Args ...argss) : c(c) { arguments = vec<std::string>(argss...); }
589  FunkConst(double c) : c(c) { arguments.resize(0); }
590 
591  double value(const std::vector<double> & data, size_t bindID)
592  {
593  (void)data;
594  (void)bindID;
595  return c;
596  }
597 
598  private:
599  double c;
600  };
601  template <typename... Args>
602  inline Funk one(Args... argss) { return Funk(new FunkConst(1., argss...)); }
603  template <typename... Args>
604  inline Funk zero(Args... argss) { return Funk(new FunkConst(0., argss...)); }
605  template <typename... Args>
606  inline Funk cnst(double x, Args... argss) { return Funk(new FunkConst(x, argss...)); }
607 
608 
609  //
610  // Derived class that implements setting of parameters
611  //
612 
613  class FunkDerived: public FunkBase
614  {
615  public:
616  FunkDerived(Funk f, std::string arg, Funk g) : my_arg(arg)
617  {
618  setup(f, arg, g);
619  };
620 
621  FunkDerived(Funk f, std::string arg, double x) : my_arg(arg)
622  {
623  setup(f, arg, cnst(x));
624  }
625 
626  // We need to sneak in an additional parameter
627  void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap)
628  {
629  functions[1]->resolve(datamap, datalen, bindID, argmap); // resolve g
630  // add new slot for result from of g
631  if(my_index.size() <= bindID)
632  {
633  my_index.resize(bindID+1);
634  // TODO: Introduce informative error message if bind fails
635  // because of inconsistencies, e.g. along the lines:
636  //
637  // std::cout << "FATAL ERROR: bind() attempts to resolve dependencies for FunkDerived" << std::endl;
638  // std::cout << "object inconsistently." << std::endl;
639  // std::cout << "Encountered while resolving " << my_arg << " and:" << std::endl;
640  // for (auto it = datamap.begin(); it != datamap.end(); it++)
641  // {
642  // std::cout << " " << it->first << std::endl;
643  // }
644  // exit(1);
645  }
646  if ( argmap.find(my_arg) == argmap.end() )
647  {
648  if(datamap.find(my_arg) == datamap.end())
649  {
650  my_index[bindID] = datalen;
651  argmap[my_arg] = datalen;
652  ++datalen;
653  }
654  else
655  {
656  my_index[bindID] = datamap[my_arg];
657  }
658  }
659  else
660  {
661  my_index[bindID] = argmap[my_arg];
662  }
663  datamap[my_arg] = my_index[bindID]; // add or overwrite entry in datamap
664  functions[0]->resolve(datamap, datalen, bindID, argmap); // resolve f
665  }
666 
667  double value(const std::vector<double> & data, size_t bindID)
668  {
669  std::vector<double> data2(data);
670  data2[my_index[bindID]] = functions[1]->value(data, bindID);
671  return functions[0]->value(data2, bindID);
672  }
673 
674  private:
675  std::string my_arg;
676 
677 
678  void setup(Funk f, std::string arg, Funk g)
679  {
680  functions = vec(f, g);
681  Singularities tmp_singl = f->getSingl();
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);
685  arguments = joinArgs(eraseArg(f->getArgs(), arg), g->getArgs());
686  };
687  std::vector<size_t> my_index;
688  };
689 
690 
691  //
692  // Derived class for the import of plain functions
693  //
694 
695  template <bool threadsafe, typename... funcargs>
696  class FunkFunc : public FunkBase
697  {
698  public:
699  template <typename... Args>
700  FunkFunc(double (*f)(funcargs...), Args... argss)
701  {
702  ptr = f;
703  digest_input(argss...);
704  }
705 
706  double value(const std::vector<double> & data, size_t bindID)
707  {
708  std::tuple<typename std::remove_reference<funcargs>::type...> my_input;
709  double result;
710  size_t i = 0;
711  #pragma omp critical (FunkFunc_setInput)
712  {
713  for ( auto f = functions.begin(); f != functions.end(); ++f, ++i)
714  {
715  *map[i] = (*f)->value(data, bindID);
716  }
717  my_input = input;
718  }
719  if(threadsafe)
720  {
721  result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
722  }
723  else
724  {
725  #pragma omp critical (FunkFunc_externalFunctionCall)
726  {
727  result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
728  }
729  }
730  return result;
731  }
732 
733  template <size_t... Args>
734  double ppp(index_list<Args...>, std::tuple<typename std::remove_reference<funcargs>::type...> & my_input)
735  {
736  return (*ptr)(std::get<Args>(my_input)...);
737  }
738 
739  private:
740  std::tuple<typename std::remove_reference<funcargs>::type...> input;
741  std::vector<double*> map;
742  double (*ptr)(funcargs...);
743 
744  // Digest input parameters
745  // (forwarding everything except daFunk::Funk types, which is mapped onto
746  // funktion parameters)
747  template<typename T, typename... Args>
748  void digest_input(T x, Args... argss)
749  {
750  const int i = sizeof...(funcargs) - sizeof...(argss) - 1;
751  std::get<i>(input) = x;
752  digest_input(argss...);
753  }
754  template<typename... Args>
755  void digest_input(Funk f, Args... argss)
756  {
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...);
763  }
764  void digest_input() {};
765  };
766 
767  template <typename... funcargs, typename... Args>
768  Funk func(double (*f)(funcargs...), Args... args) {
769  return Funk(new FunkFunc<false, funcargs...>(f, args...));
770  }
771  // Version that assumes the function to be threadsafe
772  template <typename... funcargs, typename... Args>
773  Funk func_fromThreadsafe(double (*f)(funcargs...), Args... args) {
774  return Funk(new FunkFunc<true, funcargs...>(f, args...));
775  }
776 
777 
778  template <bool threadsafe, typename O, typename... funcargs>
779  class FunkFuncM : public FunkBase
780  {
781  public:
782  template <typename... Args>
783  FunkFuncM(O* obj, double (O::* f)(funcargs...), Args... argss) : obj(obj)
784  {
785  ptr = f;
786  digest_input(argss...);
787  }
788 
789  template <typename... Args>
790  FunkFuncM(shared_ptr<O> obj, double (O::* f)(funcargs...), Args... argss) : shared_obj(obj), obj(&*obj)
791  {
792  ptr = f;
793  digest_input(argss...);
794  }
795 
796  double value(const std::vector<double> & data, size_t bindID)
797  {
798  std::tuple<typename std::remove_reference<funcargs>::type...> my_input;
799  double result;
800  size_t i = 0;
801  #pragma omp critical(FunkFuncM_value)
802  {
803  for ( auto f = functions.begin(); f != functions.end(); ++f, ++i)
804  {
805  *map[i] = (*f)->value(data, bindID);
806  }
807  my_input = input;
808  }
809  if(threadsafe)
810  {
811  result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
812  }
813  else
814  {
815  #pragma omp critical (FunkFuncM_objectFunctionCall)
816  {
817  result = ppp(typename detail::range_builder<0, sizeof...(funcargs)>::type(), my_input);
818  }
819  }
820  return result;
821  }
822 
823  template <size_t... Args>
824  double ppp(index_list<Args...>, std::tuple<typename std::remove_reference<funcargs>::type...> & my_input)
825  {
826  return (*obj.*ptr)(std::get<Args>(my_input)...);
827  }
828 
829  private:
830  std::tuple<typename std::remove_reference<funcargs>::type...> input;
831  std::vector<double*> map;
832  double (O::* ptr)(funcargs...);
833  shared_ptr<O> shared_obj;
834  O* obj;
835 
836  // Digest input parameters
837  // (forwarding everything except daFunk::Funk types, which is mapped onto
838  // funktion parameters)
839  template<typename T, typename... Args>
840  void digest_input(T x, Args... argss)
841  {
842  const int i = sizeof...(funcargs) - sizeof...(argss) - 1;
843  std::get<i>(input) = x;
844  digest_input(argss...);
845  }
846  template<typename... Args>
847  void digest_input(Funk f, Args... argss)
848  {
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...);
855  }
856  void digest_input() {};
857  };
858 
859 
860  template <typename O, typename... funcargs, typename... Args>
861  Funk funcM(O* obj, double (O::* f)(funcargs...), Args... args) {
862  return Funk(new FunkFuncM<false, O, funcargs...>(obj, f, args...));
863  }
864  template <typename O, typename... funcargs, typename... Args>
865  Funk funcM(shared_ptr<O> obj, double (O::* f)(funcargs...), Args... args) {
866  return Funk(new FunkFuncM<false, O, funcargs...>(obj, f, args...));
867  }
868 
869  // Versions that assume the object to be threadsafe
870  template <typename O, typename... funcargs, typename... Args>
871  Funk funcM_fromThreadsafe(O* obj, double (O::* f)(funcargs...), Args... args) {
872  return Funk(new FunkFuncM<true, O, funcargs...>(obj, f, args...));
873  }
874  template <typename O, typename... funcargs, typename... Args>
875  Funk funcM_fromThreadsafe(shared_ptr<O> obj, double (O::* f)(funcargs...), Args... args) {
876  return Funk(new FunkFuncM<true, O, funcargs...>(obj, f, args...));
877  }
878 
879  //
880  // Derived class that implements delta function
881  //
882 
883  class FunkDelta: public FunkBase
884  {
885  public:
886  FunkDelta(std::string arg, double pos, double width) : pos(pos), width(width)
887  {
888  arguments = vec(arg);
889  this->set_singularity("v", pos, width);
890  singularities[arg].push_back(std::pair<Funk, Funk>(cnst(pos), cnst(width)));
891  }
892 
893  double value(const std::vector<double> & data, size_t bindID)
894  {
895  double x = data[indices[bindID][0]];
896  return exp(-pow(x-pos,2)/pow(width,2)/2)/sqrt(2*M_PI)/width;
897  }
898 
899  private:
900  double pos, width;
901  };
902  inline Funk delta(std::string arg, double pos, double width) { return Funk(new FunkDelta(arg, pos, width)); }
903 
904  //
905  // Derived class that implements simple linear variable
906  //
907 
908  class FunkVar: public FunkBase
909  {
910  public:
911  FunkVar(std::string arg)
912  {
913  arguments = vec(arg);
914  }
915 
916  double value(const std::vector<double> & data, size_t bindID)
917  {
918  return data[indices[bindID][0]];
919  }
920  };
921  inline Funk var(std::string arg) { return Funk(new FunkVar(arg)); }
922 
923 
924  //
925  // Definition of FunkBase member functions
926  //
927 
928  inline Funk FunkBase::set_singularity(std::string arg, Funk pos, Funk width)
929  {
930  singularities[arg].push_back(std::pair<Funk, Funk>(pos, width));
931  return shared_from_this();
932  };
933  inline Funk FunkBase::set_singularity(std::string arg, double pos, Funk width)
934  { return shared_from_this()->set_singularity(arg, cnst(pos), width); }
935  inline Funk FunkBase::set_singularity(std::string arg, double pos, double width)
936  { return shared_from_this()->set_singularity(arg, cnst(pos), cnst(width)); }
937  inline Funk FunkBase::set_singularity(std::string arg, Funk pos, double width)
938  { return shared_from_this()->set_singularity(arg, pos, cnst(width)); }
939 
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)
941  {
942  // Resolve my dependencies
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 )
948  {
949  try
950  {
951  *it2 = datamap.at(*it1);
952  }
953  catch (std::out_of_range & e)
954  {
955  std::string msg = "FunkBase::resolve() encountered internal problem when resolving " + *it1 + ".\n";
956  msg+= " --- Actual arguments of object: " + args_string(arguments);
957 #ifdef GAMBIT_DIR
958  Gambit::utils_error().raise(LOCAL_INFO, msg);
959 #else
960  throw std::invalid_argument(msg);
961 #endif
962  }
963  }
964 
965  // Set indices
966  if(indices.size() <= bindID)
967  {
968  indices.resize(bindID+1);
969  // TODO: Throw error if problems are encountered
970  //std::cout << "FATAL ERROR: bind() attempts to resolve dependencies of FunkBase" << std::endl;
971  //std::cout << "object inconsistently." << std::endl;
972  //std::cout << "Encountered while resolving: " << std::endl;
973  //for (auto it = datamap.begin(); it != datamap.end(); it++)
974  //{
975  // std::cout << " " << it->first << std::endl;
976  //}
977  //exit(1);
978  }
979  indices[bindID] = tmp_indices;
980 
981  // Resolve other dependencies
982  for (auto it = functions.begin(); it != functions.end(); ++it)
983  {
984  (*it)->resolve(datamap, datalen, bindID, argmap);
985  }
986 
987  }
988 
989  template <typename... Args> inline bool FunkBase::assert_args(Args... args)
990  {
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 )
994  {
995  std::set<std::string> theirargs(it->begin(), it->end());
996  if ( myargs == theirargs )
997  return true;
998  }
999  return false;
1000 
1001  }
1002 
1003  template <typename... Args> inline Funk FunkBase::set(std::string arg, Funk g, Args... args)
1004  {
1005  Funk f = shared_from_this();
1006  if ( std::find(arguments.begin(), arguments.end(), arg) != arguments.end() )
1007  {
1008  f = Funk(new FunkDerived(f, arg, g));
1009  }
1010  else
1011  {
1012  std::cout << "daFunk::FunkBase WARNING: Ignoring \"" << arg << "\" = function." << std::endl;
1013  }
1014  return f->set(args...);
1015  }
1016 
1017  template <typename... Args> inline Funk FunkBase::set(std::string arg, std::string arg1, Args... args)
1018  { return shared_from_this()->set(arg, var(arg1))->set(args...); }
1019 
1020  template <typename... Args> inline Funk FunkBase::set(std::string arg, double x, Args... args)
1021  { return shared_from_this()->set(arg, cnst(x))->set(args...); }
1022 
1023  template <> inline Funk FunkBase::set()
1024  { return shared_from_this(); }
1025 
1026  template <typename... Args> inline shared_ptr<FunkBound> FunkBase::bind(Args... argss)
1027  {
1028  size_t bindID;
1029  FunkBound::bindID_manager(bindID,true);
1030  std::map<std::string, size_t> datamap;
1031  size_t i;
1032  auto bound_arguments = vec<std::string>(argss...);
1033  datalen = bound_arguments.size();
1034  i = 0;
1035  for ( auto it = bound_arguments.begin(); it != bound_arguments.end(); ++it, ++i )
1036  {
1037  datamap[*it] = i;
1038  }
1039  std::map<std::string,size_t> argmap;
1040  if ( not args_match(arguments, bound_arguments) )
1041  {
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);
1045 #ifdef GAMBIT_DIR
1046  Gambit::utils_error().raise(LOCAL_INFO, msg);
1047 #else
1048  throw std::invalid_argument(msg);
1049 #endif
1050  }
1051  this->resolve(datamap, datalen, bindID, argmap);
1052  return shared_ptr<FunkBound>(new FunkBound(shared_from_this(), datalen, bindID));
1053  }
1054 
1055  inline bool FunkBase::hasArg(std::string arg)
1056  {
1057  return ( std::find(arguments.begin(), arguments.end(), arg) != arguments.end() );
1058  }
1059 
1060  inline bool FunkBase::hasArgs()
1061  {
1062  return ( this->arguments.size() != 0 );
1063  }
1064 
1065  template <typename... Args> inline shared_ptr<FunkIntegrate_gsl1d> FunkBase::gsl_integration(Args... args)
1066  {
1067  return getIntegrate_gsl1d(shared_from_this(), args...);
1068  }
1069 
1071  {
1072  std::cout << "Arguments:";
1073  for ( auto it = arguments.begin(); it != arguments.end(); it++ )
1074  {
1075  std::cout << " \"" << *it << "\"";
1076  }
1077  if ( arguments.size() == 0 )
1078  std::cout << " none";
1079  std::cout << std::endl;
1080  for ( auto it = singularities.begin() ; it != singularities.end(); ++it )
1081  {
1082  std::cout << "Singularities in " << it->first << ": " << it->second.size() << std::endl;
1083  }
1084  return shared_from_this();
1085  }
1086 
1087  inline PlainPtrs1 FunkBase::plain(std::string arg1)
1088  {
1089  void* ptr = new FunkPlain(shared_from_this(), arg1);
1090  return PlainPtrs1(&FunkPlain::plain1p, ptr);
1091  }
1092  inline PlainPtrs2 FunkBase::plain(std::string arg1, std::string arg2)
1093  {
1094  void* ptr = new FunkPlain(shared_from_this(), arg1, arg2);
1095  return PlainPtrs2(&FunkPlain::plain2p, ptr);
1096  }
1097  inline PlainPtrs3 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3)
1098  {
1099  void* ptr = new FunkPlain(shared_from_this(), arg1, arg2, arg3);
1100  return PlainPtrs3(&FunkPlain::plain3p, ptr);
1101  }
1102  inline PlainPtrs4 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3, std::string arg4)
1103  {
1104  void* ptr = new FunkPlain(shared_from_this(), arg1, arg2, arg3, arg4);
1105  return PlainPtrs4(&FunkPlain::plain4p, ptr);
1106  }
1107 
1108  template <typename T>
1109  inline PlainPtr1 FunkBase::plain(std::string arg1)
1110  {
1111  T::set(new FunkPlain(shared_from_this(), arg1));
1112  return &FunkPlain::plain1<T>;
1113  }
1114  template <typename T>
1115  inline PlainPtr2 FunkBase::plain(std::string arg1, std::string arg2)
1116  {
1117  T::set(new FunkPlain(shared_from_this(), arg1, arg2));
1118  return &FunkPlain::plain2<T>;
1119  }
1120  template <typename T>
1121  inline PlainPtr3 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3)
1122  {
1123  T::set(new FunkPlain(shared_from_this(), arg1, arg2, arg3));
1124  return &FunkPlain::plain3<T>;
1125  }
1126  template <typename T>
1127  inline PlainPtr4 FunkBase::plain(std::string arg1, std::string arg2, std::string arg3, std::string arg4)
1128  {
1129  T::set(new FunkPlain(shared_from_this(), arg1, arg2, arg3, arg4));
1130  return &FunkPlain::plain4<T>;
1131  }
1132 
1133 
1134  //
1135  // Mathematical functions from cmath
1136  //
1137 
1138  // Unary minus sign
1139  class FunkMath_umin: public FunkBase
1140  {
1141  public:
1143  {
1144  functions = vec(f);
1145  singularities = f->getSingl();
1146  arguments = f->getArgs();
1147  }
1148  double value(const std::vector<double> & data, size_t bindID)
1149  {
1150  return -(functions[0]->value(data, bindID));
1151  }
1152  };
1153  inline Funk operator - (Funk f) { return Funk(new FunkMath_umin(f)); }
1154 
1155  // Unary operations
1156 #define MATH_OPERATION(OPERATION) \
1157  class FunkMath_##OPERATION: public FunkBase \
1158  { \
1159  public: \
1160  FunkMath_##OPERATION(Funk f) \
1161  { \
1162  functions = vec(f); \
1163  arguments = f->getArgs(); \
1164  singularities = f->getSingl(); \
1165  } \
1166  double value(const std::vector<double> & data, size_t bindID) \
1167  { \
1168  return OPERATION(functions[0]->value(data, bindID)); \
1169  } \
1170  }; \
1171  inline Funk OPERATION (Funk f) { return Funk(new FunkMath_##OPERATION(f)); }
1172  MATH_OPERATION(cos)
1173  MATH_OPERATION(sin)
1174  MATH_OPERATION(tan)
1175  MATH_OPERATION(acos)
1176  MATH_OPERATION(asin)
1177  MATH_OPERATION(atan)
1178  MATH_OPERATION(cosh)
1179  MATH_OPERATION(sinh)
1180  MATH_OPERATION(tanh)
1181  MATH_OPERATION(acosh)
1182  MATH_OPERATION(asinh)
1183  MATH_OPERATION(atanh)
1184  MATH_OPERATION(exp)
1185  MATH_OPERATION(log)
1186  MATH_OPERATION(log10)
1187  MATH_OPERATION(sqrt)
1188  MATH_OPERATION(fabs)
1189 #undef MATH_OPERATION
1190 
1191  // Standard binary operations
1192 #define MATH_OPERATION(OPERATION, SYMBOL) \
1193  class FunkMath_##OPERATION: public FunkBase \
1194  { \
1195  public: \
1196  FunkMath_##OPERATION(Funk f1, Funk f2) \
1197  { \
1198  functions = vec(f1, f2); \
1199  arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
1200  singularities = joinSingl(f1->getSingl(), f2->getSingl());\
1201  } \
1202  FunkMath_##OPERATION(double x, Funk f2) \
1203  { \
1204  auto f1 = cnst(x); \
1205  functions = vec(f1, f2); \
1206  arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
1207  singularities = joinSingl(f1->getSingl(), f2->getSingl());\
1208  } \
1209  FunkMath_##OPERATION(Funk f1, double x) \
1210  { \
1211  auto f2 = cnst(x); \
1212  functions = vec(f1, f2); \
1213  arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
1214  singularities = joinSingl(f1->getSingl(), f2->getSingl());\
1215  } \
1216  double value(const std::vector<double> & data, size_t bindID) \
1217  { \
1218  return functions[0]->value(data, bindID) SYMBOL functions[1]->value(data, bindID); \
1219  } \
1220  }; \
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)); }
1224  MATH_OPERATION(Sum,+)
1225  MATH_OPERATION(Mul,*)
1226  MATH_OPERATION(Div,/)
1227  MATH_OPERATION(Dif,-)
1228 #undef MATH_OPERATION
1229 
1230  // More binary operations
1231 #define MATH_OPERATION(OPERATION) \
1232  class FunkMath_##OPERATION: public FunkBase \
1233  { \
1234  public: \
1235  FunkMath_##OPERATION(Funk f1, Funk f2) \
1236  { \
1237  functions = vec(f1, f2); \
1238  arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
1239  singularities = joinSingl(f1->getSingl(), f2->getSingl()); \
1240  } \
1241  FunkMath_##OPERATION(double x, Funk f2) \
1242  { \
1243  auto f1 = cnst(x); \
1244  functions = vec(f1, f2); \
1245  arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
1246  singularities = joinSingl(f1->getSingl(), f2->getSingl()); \
1247  } \
1248  FunkMath_##OPERATION(Funk f1, double x) \
1249  { \
1250  auto f2 = cnst(x); \
1251  functions = vec(f1, f2); \
1252  arguments = joinArgs(f1->getArgs(), f2->getArgs()); \
1253  singularities = joinSingl(f1->getSingl(), f2->getSingl()); \
1254  } \
1255  double value(const std::vector<double> & data, size_t bindID) \
1256  { \
1257  return OPERATION(functions[0]->value(data, bindID), functions[1]->value(data, bindID)); \
1258  } \
1259  }; \
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)); }
1264  MATH_OPERATION(fmin)
1265  MATH_OPERATION(fmax)
1266 #undef MATH_OPERATION
1267 
1269  // *** End of core implementation ***
1271 
1272 
1273 
1275  // *** Extensions ***
1277 
1278 
1279  //
1280  // Derived class: 1dim linear or logarithmic interpolation
1281  //
1282 
1283  class FunkInterp : public FunkBase
1284  {
1285  public:
1286  FunkInterp(Funk f, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin")
1287  {
1288  setup(f, Xgrid, Ygrid, mode);
1289  }
1290  FunkInterp(std::string arg, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin")
1291  {
1292  setup(var(arg), Xgrid, Ygrid, mode);
1293  }
1294  FunkInterp(double x, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode = "lin")
1295  {
1296  setup(cnst(x), Xgrid, Ygrid, mode);
1297  }
1298 
1299  double value(const std::vector<double> & data, size_t bindID)
1300  {
1301  functions[0]->value(data, bindID);
1302  return (this->*ptr)(data[indices[bindID][0]]);
1303  }
1304 
1305  private:
1306  void setup(Funk f, std::vector<double> & Xgrid, std::vector<double> & Ygrid, std::string mode)
1307  {
1308  // TODO: Catch invalid setup
1309  functions = vec(f);
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;
1316  }
1317 
1318  double logInterp(double x)
1319  {
1320  // Linear interpolation in log-log space
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));
1329  }
1330 
1331  double linearInterp(double x)
1332  {
1333  // Linear interpolation in lin-lin space
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);
1342  }
1343 
1344  double(FunkInterp::*ptr)(double);
1345  std::vector<double> Xgrid;
1346  std::vector<double> Ygrid;
1347  std::string mode;
1348  };
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)); }
1350 
1351 
1352  //
1353  // Basic if..else clause
1354  //
1355 
1356  class FunkIfElse: public FunkBase
1357  {
1358  public:
1360  {
1361  functions = vec(f, g, h);
1362  singularities = joinSingl(joinSingl(f->getSingl(), g->getSingl()), h->getSingl());
1363  arguments = joinArgs(joinArgs(f->getArgs(), g->getArgs()), h->getArgs());
1364  }
1365  double value(const std::vector<double> & data, size_t bindID)
1366  {
1367  if ( functions[0]->value(data,bindID) >= 0. )
1368  return functions[1]->value(data,bindID);
1369  else
1370  return functions[2]->value(data,bindID);
1371  }
1372  };
1373  inline Funk ifelse(Funk f, Funk g, Funk h) { return Funk(new FunkIfElse(f, g, h)); }
1374  inline Funk ifelse(Funk f, double g, Funk h) { return Funk(new FunkIfElse(f, cnst(g), h)); }
1375  inline Funk ifelse(Funk f, double g, double h) { return Funk(new FunkIfElse(f, cnst(g), cnst(h))); }
1376  inline Funk ifelse(Funk f, Funk g, double h) { return Funk(new FunkIfElse(f, g, cnst(h))); }
1377 
1378 
1379  //
1380  // Throw errors when called
1381  //
1382 
1383  class ThrowError: public FunkBase
1384  {
1385  public:
1386  ThrowError(std::string msg) : msg(msg)
1387  {
1388  }
1389  double value(const std::vector<double> & data, size_t bindID)
1390  {
1391  (void)bindID;
1392  (void)data;
1393 #ifdef GAMBIT_DIR
1394  if ( omp_get_level() == 0 ) // Outside of OMP blocks
1395  {
1396  Gambit::utils_error().raise(LOCAL_INFO, "daFunk::ThrowError says: " + msg);
1397  }
1398  else // Inside of OMP blocks
1399  {
1400  Gambit::piped_errors.request(LOCAL_INFO, "daFunk::ThrowError says: " + msg);
1401  }
1402 #else
1403  throw std::invalid_argument("daFunk::ThrowError says: " + msg);
1404 #endif
1405  return 0;
1406  }
1407 
1408  private:
1409  std::string msg; // Error message to throw when function is called
1410  };
1411  inline Funk throwError(std::string msg) { return Funk(new ThrowError(msg)); }
1412 
1413 #ifdef GAMBIT_DIR
1414  class RaiseInvalidPoint: public FunkBase
1415  {
1416  public:
1417  RaiseInvalidPoint(std::string msg) : msg(msg)
1418  {
1419  }
1420  double value(const std::vector<double> & data, size_t bindID)
1421  {
1422  (void)bindID;
1423  (void)data;
1424  if ( omp_get_level() == 0 ) // Outside of OMP blocks
1425  {
1426  Gambit::utils_warning().raise(LOCAL_INFO, "daFunk::RaiseInvalidPoint says: " + msg);
1427  Gambit::invalid_point().raise("daFunk::RaiseInvalidPoint says: " + msg);
1428  }
1429  else // Inside OMP blocks
1430  {
1431  Gambit::utils_warning().raise(LOCAL_INFO, "daFunk::RaisePipedInvalidPoint says: " + msg);
1432  Gambit::piped_invalid_point.request("daFunk::RaisePipedInvalidPoint says: " + msg);
1433  }
1434  return 0;
1435  }
1436 
1437  private:
1438  std::string msg; // Error message to throw when function is called
1439  };
1440  inline Funk raiseInvalidPoint(std::string msg) { return Funk(new RaiseInvalidPoint(msg)); }
1441 #endif
1442 
1443 
1444  //
1445  // Prints message when called
1446  //
1447 
1448  class Bottle: public FunkBase
1449  {
1450  public:
1451  Bottle(Funk f, std::string msg) : msg(msg)
1452  {
1453  functions = vec(f);
1454  singularities = f->getSingl();
1455  arguments = f->getArgs();
1456  }
1457  double value(const std::vector<double> & data, size_t bindID)
1458  {
1459  std::cout << "daFunk::Message says:\n" << msg << std::endl;
1460  return functions[0]->value(data, bindID);
1461  }
1462 
1463  private:
1464  std::string msg; // Message in a bottle. Ha!
1465  };
1466  //inline Funk print(std::string msg) { return Funk(new Bottle(msg)); }
1467  inline Funk FunkBase::print(std::string msg)
1468  { return Funk(new Bottle(shared_from_this(), msg)); };
1469 
1470 
1471  //
1472  // GSL integration
1473  //
1474 
1475  class FunkIntegrate_gsl1d: public FunkBase, public gsl_function
1476  {
1477  public:
1478  FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f1, Funk f2)
1479  {
1480  setup(f0, arg, f1, f2);
1481  }
1482  FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, Funk f)
1483  {
1484  setup(f0, arg, cnst(x), f);
1485  }
1486  FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, double y)
1487  {
1488  setup(f0, arg, cnst(x), cnst(y));
1489  }
1490  FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, double x)
1491  {
1492  setup(f0, arg, f, cnst(x));
1493  }
1494  FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, Funk f)
1495  {
1496  setup(f0, arg, var(x), f);
1497  }
1498  FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, std::string y)
1499  {
1500  setup(f0, arg, var(x), var(y));
1501  }
1502  FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, std::string x)
1503  {
1504  setup(f0, arg, f, var(x));
1505  }
1506  FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, double y)
1507  {
1508  setup(f0, arg, var(x), cnst(y));
1509  }
1510  FunkIntegrate_gsl1d(Funk f0, std::string arg, double y, std::string x)
1511  {
1512  setup(f0, arg, cnst(y), var(x));
1513  }
1514 
1515  void resolve(std::map<std::string, size_t> datamap, size_t & datalen, size_t bindID, std::map<std::string,size_t> &argmap)
1516  {
1517  functions[1]->resolve(datamap, datalen, bindID, argmap); // Resolve boundary 0
1518  functions[2]->resolve(datamap, datalen, bindID, argmap); // Resolve boundary 1
1519 
1520  // Set indices
1521  if(index.size() <= bindID)
1522  {
1523  index.resize(bindID+1);
1524  // TODO: Throw error if problems are encountered
1525  //std::cout << "FATAL ERROR: bind() attempts to resolve dependencies of FunkIntegrate_gsl1d" << std::endl;
1526  //std::cout << "object inconsistently." << std::endl;
1527  //std::cout << "Encountered while resolving " << arg << " and:" << std::endl;
1528  //for (auto it = datamap.begin(); it != datamap.end(); it++)
1529  //{
1530  // std::cout << " " << it->first << std::endl;
1531  //}
1532  //exit(1);
1533  }
1534  if ( argmap.find(arg) == argmap.end() )
1535  {
1536  if(datamap.find(arg) == datamap.end())
1537  {
1538  index[bindID] = datalen;
1539  argmap[arg] = datalen;
1540  ++datalen;
1541  }
1542  else
1543  {
1544  index[bindID] = datamap[arg];
1545  }
1546  }
1547  else
1548  {
1549  index[bindID] = argmap[arg];
1550  }
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 )
1554  {
1555  it->first->resolve(datamap, datalen, bindID, argmap);
1556  it->second->resolve(datamap, datalen, bindID, argmap);
1557  }
1558  }
1559 
1561  {
1562  gsl_integration_workspace_free(gsl_workspace);
1563  }
1564 
1565  shared_ptr<FunkIntegrate_gsl1d> set_epsrel(double epsrel)
1566  { this->epsrel = epsrel; return static_pointer_cast<FunkIntegrate_gsl1d>(this->FunkIntegrate_gsl1d::shared_from_this()); }
1567  shared_ptr<FunkIntegrate_gsl1d> set_epsabs(double epsabs)
1568  { this->epsabs = epsabs; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
1569  shared_ptr<FunkIntegrate_gsl1d> set_limit(size_t limit)
1570  { this->limit = limit; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
1571  shared_ptr<FunkIntegrate_gsl1d> set_singularity_factor(double f)
1572  { this->singl_factor = f; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
1573  shared_ptr<FunkIntegrate_gsl1d> set_use_log_fallback(bool flag)
1574  { this->use_log_fallback = flag; return static_pointer_cast<FunkIntegrate_gsl1d>(this->shared_from_this()); }
1575 
1576  double value(const std::vector<double> & data, size_t bindID)
1577  {
1578  double result;
1579  #pragma omp critical(FunkIntegrate_gsl1d_integration)
1580  {
1581  local_data = data;
1582  local_bindID = bindID;
1583  double error;
1584  function=&FunkIntegrate_gsl1d::invoke;
1585  params=this;
1586  double x0 = functions[1]->value(data, bindID);
1587  double x1 = functions[2]->value(data, bindID);
1588  gsl_set_error_handler_off();
1589  int status = 0;
1590  if ( my_singularities.size() == 0 )
1591  {
1592  status = gsl_integration_qags(this, x0, x1, epsabs, epsrel, limit, gsl_workspace, &result, &error);
1593  }
1594  else
1595  {
1596  double s = 0;
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 )
1601  {
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;
1606  if ( z0 == z1 )
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);
1610  }
1611  std::sort(ranges.begin(), ranges.end());
1612  for ( auto it = ranges.begin(); it != ranges.end()-1; ++it )
1613  {
1614  status = gsl_integration_qags(this, *it, *(it+1), epsabs, epsrel, limit, gsl_workspace, &result, &error);
1615  s += result;
1616  if (status) break;
1617  }
1618  result = s;
1619  }
1620  if (status and this->use_log_fallback)
1621  {
1622  // The last resort: A cheap integration on log grid, linear interpolation
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++)
1630  {
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;
1635  y0 = y1;
1636  }
1637  result = sum;
1638  }
1639  // TODO: Implement flags to optionally throw an error
1640  if (status and not this->use_log_fallback)
1641  {
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();
1647 // std::cout << "Dumping integrand:" << std::endl;
1648 // for ( double x = x0; x <= x1; x = (x0>0) ? x*1.01 : x+(x1-x0)/1000)
1649 // std::cerr << " " << x << " " << invoke(x, this) << std::endl;
1650  std::cerr << "Returning zero." << std::endl;
1651  result = 0.;
1652  }
1653  }
1654  return result;
1655  }
1656 
1657  private:
1658  void setup(Funk f0, std::string arg, Funk f1, Funk f2)
1659  {
1660  this->functions = vec(f0, f1, f2);
1661 
1662  singularities = joinSingl(f1->getSingl(), f2->getSingl());
1663  if ( f0->getSingl().find(arg) != f0->getSingl().end() )
1664  my_singularities = f0->getSingl()[arg];
1665  Singularities tmp_singl = f0->getSingl();
1666  tmp_singl.erase(arg);
1667  singularities = joinSingl(singularities, tmp_singl);
1668 
1669  arguments = joinArgs(eraseArg(f0->getArgs(), arg), joinArgs(f1->getArgs(), f2->getArgs()));
1670  gsl_workspace = gsl_integration_workspace_alloc (100000);
1671 
1672  this->arg = arg;
1673  limit = 100;
1674  epsrel = 1e-2;
1675  epsabs = 1e-2;
1676  use_log_fallback = false;
1677  singl_factor = 4;
1678  }
1679 
1680  // Static member function that invokes integrand
1681  static double invoke(double x, void *params) {
1682  FunkIntegrate_gsl1d * ptr = static_cast<FunkIntegrate_gsl1d*>(params);
1683  ptr->local_data[ptr->index[ptr->local_bindID]] = x;
1684  return ptr->functions[0]->value(ptr->local_data, ptr->local_bindID);
1685  }
1686 
1687  // Required for rewiring input parameters
1688  std::vector<double> local_data;
1690  std::vector<std::pair<Funk, Funk>> my_singularities;
1691 
1692  // Integration range and function pointer
1693  std::string arg;
1694 
1695  // GSL workspace and parameters
1696  gsl_integration_workspace * gsl_workspace;
1697  size_t limit;
1698  std::vector<size_t> index;
1699  double epsrel;
1700  double epsabs;
1702 
1704  };
1705 
1706  // Standard behaviour
1707  template <typename T1, typename T2>
1708  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)); }
1709 
1710  //
1711  // Helper function for producing singularity-augmented x-grids
1712  //
1713 
1714  inline std::vector<double> augmentSingl(const std::vector<double> & xgrid, Funk f, int N = 100, double sigma = 3)
1715  {
1716  std::vector<double> result = xgrid;
1717  double xmin = result.front();
1718  double xmax = result.back();
1719 
1720  if ( f->getNArgs() != 1 )
1721  {
1722  std::string msg = "augment_with_singularities(): takes only functions with one argument.\n";
1723  msg+= " --- Actual arguments are: " + args_string(f->getArgs());
1724 #ifdef GAMBIT_DIR
1725  Gambit::utils_error().raise(LOCAL_INFO, msg);
1726 #else
1727  throw std::invalid_argument(msg);
1728 #endif
1729  }
1730 
1731  std::string arg = f->getArgs()[0];
1732  Singularities singlsMap = f->getSingl();
1733 
1734  if ( singlsMap.find(arg) != singlsMap.end() )
1735  {
1736  auto singls = singlsMap.at(arg);
1737  for (auto it = singls.begin(); it != singls.end(); it ++)
1738  {
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());
1743  }
1744  std::sort(result.begin(), result.end());
1745  }
1746 
1747  return result;
1748  }
1749 }
1750 
1751 
1752 #endif // __DAFUNK_HPP__
FunkFuncM(shared_ptr< O > obj, double(O::*f)(funcargs...), Args... argss)
Definition: daFunk.hpp:790
shared_ptr< FunkIntegrate_gsl1d > set_use_log_fallback(bool flag)
Definition: daFunk.hpp:1573
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)
Definition: daFunk.hpp:488
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.
greatScanData data
Definition: great.cpp:38
FunkPlain(Funk fin, std::string arg1, std::string arg2, std::string arg3, std::string arg4)
Definition: daFunk.hpp:519
double(* PlainPtr3)(double &, double &, double &)
Definition: daFunk.hpp:124
double(* fptr)(int &)
Pointer to a function that takes an integer by reference and returns a double.
Funk throwError(std::string msg)
Definition: daFunk.hpp:1411
std::pair< double(*)(double, void *), void * > PlainPtrs1
Definition: daFunk.hpp:129
FunkIfElse(Funk f, Funk g, Funk h)
Definition: daFunk.hpp:1359
std::vector< double > vect(Args... argss)
Definition: daFunk.hpp:406
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)
Definition: daFunk.hpp:518
General small utility macros.
double ppp(index_list< Args... >, std::tuple< typename std::remove_reference< funcargs >::type... > &my_input)
Definition: daFunk.hpp:734
Bottle(Funk f, std::string msg)
Definition: daFunk.hpp:1451
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)
Definition: daFunk.hpp:627
FunkIntegrate_gsl1d(Funk f0, std::string arg, double y, std::string x)
Definition: daFunk.hpp:1510
virtual void resolve(std::map< std::string, size_t > datamap, size_t &datalen, size_t bindID, std::map< std::string, size_t > &argmap)
Definition: daFunk.hpp:940
shared_ptr< FunkIntegrate_gsl1d > getIntegrate_gsl1d(Funk fptr, std::string arg, T1 x1, T2 x2)
Definition: daFunk.hpp:1708
FunkDerived(Funk f, std::string arg, double x)
Definition: daFunk.hpp:621
ArgsType arguments
Definition: daFunk.hpp:369
static double plain3p(double x1, double x2, double x3, void *ptr)
Definition: daFunk.hpp:531
FunkConst(double c, Args ...argss)
Definition: daFunk.hpp:588
std::string args_string(ArgsType args)
Definition: daFunk.hpp:216
Funk print(std::string arg)
Definition: daFunk.hpp:1467
FunkFuncM(O *obj, double(O::*f)(funcargs...), Args... argss)
Definition: daFunk.hpp:783
const std::vector< std::string > & getArgs()
Definition: daFunk.hpp:314
Funk delta(std::string arg, double pos, double width)
Definition: daFunk.hpp:902
ArgsType eraseArg(ArgsType args, std::string arg)
Definition: daFunk.hpp:236
bool args_match(ArgsType args1, ArgsType args2)
Definition: daFunk.hpp:204
Funk funcM(O *obj, double(O::*f)(funcargs...), Args... args)
Definition: daFunk.hpp:861
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:667
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:893
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, double x)
Definition: daFunk.hpp:1490
void digest_input(T x, Args... argss)
Definition: daFunk.hpp:748
void request(std::string origin, std::string message)
Request an exception.
Definition: exceptions.cpp:547
#define LOCAL_INFO
Definition: local_info.hpp:34
Singularities getSingl()
Definition: daFunk.hpp:334
FunkBound(Funk f, size_t datalen, size_t bindID)
Definition: daFunk.hpp:395
static bool & is_dead()
Definition: daFunk.hpp:385
std::pair< double(*)(double, double, double, double, void *), void * > PlainPtrs4
Definition: daFunk.hpp:132
std::string msg
Definition: daFunk.hpp:1464
void digest_input(Funk f, Args... argss)
Definition: daFunk.hpp:847
std::vector< double > linspace(double x0, double x1, unsigned int n)
Definition: daFunk.hpp:172
FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, Funk f)
Definition: daFunk.hpp:1482
std::vector< double > local_data
Definition: daFunk.hpp:1688
FunkFunc(double(*f)(funcargs...), Args... argss)
Definition: daFunk.hpp:700
std::string msg
Definition: daFunk.hpp:1409
gsl_integration_workspace * gsl_workspace
Definition: daFunk.hpp:1696
Funk var(std::string arg)
Definition: daFunk.hpp:921
void digest_input()
Definition: daFunk.hpp:856
START_MODEL x4
Definition: demo.hpp:52
static double plain2p(double x1, double x2, void *ptr)
Definition: daFunk.hpp:526
bool hasArg(std::string)
Definition: daFunk.hpp:1055
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, Funk f)
Definition: daFunk.hpp:1494
shared_ptr< FunkIntegrate_gsl1d > set_limit(size_t limit)
Definition: daFunk.hpp:1569
Piped_invalid_point piped_invalid_point
Global instance of piped invalid point class.
Definition: exceptions.cpp:544
shared_ptr< FunkIntegrate_gsl1d > set_epsrel(double epsrel)
Definition: daFunk.hpp:1565
std::tuple< typename std::remove_reference< funcargs >::type... > input
Definition: daFunk.hpp:830
Funk set_singularity(std::string arg, Funk pos, Funk width)
Definition: daFunk.hpp:928
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:1457
Funk funcM_fromThreadsafe(O *obj, double(O::*f)(funcargs...), Args... args)
Definition: daFunk.hpp:871
double value(const std::vector< double > &args, size_t bindID)
Definition: daFunk.hpp:567
void digest_input(Funk f, Args... argss)
Definition: daFunk.hpp:755
std::map< std::string, std::vector< std::pair< Funk, Funk > > > Singularities
Definition: daFunk.hpp:118
std::vector< double > vect2(std::vector< std::vector< double >> &coll)
Definition: daFunk.hpp:446
std::tuple< typename std::remove_reference< funcargs >::type... > input
Definition: daFunk.hpp:740
Funk func_fromThreadsafe(double(*f)(funcargs...), Args... args)
Definition: daFunk.hpp:773
std::vector< double * > map
Definition: daFunk.hpp:741
static double plain2(double &x1, double &x2)
Definition: daFunk.hpp:549
void digest_input(T x, Args... argss)
Definition: daFunk.hpp:840
std::size_t getNArgs()
Definition: daFunk.hpp:315
shared_ptr< O > shared_obj
Definition: daFunk.hpp:833
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:1148
shared_ptr< FunkBound > bind(Args... args)
Definition: daFunk.hpp:1026
START_MODEL x2
Definition: demo.hpp:42
ThrowError(std::string msg)
Definition: daFunk.hpp:1386
const double sigma
Definition: SM_Z.hpp:43
std::vector< Funk > functions
Definition: daFunk.hpp:368
double eval(Args... argss)
Definition: daFunk.hpp:399
virtual void raise(const std::string &)
Raise the exception, i.e. throw it. Exact override of base method.
Definition: exceptions.cpp:422
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)
Definition: daFunk.hpp:1349
shared_ptr< FunkBound > BoundFunk
Definition: daFunk.hpp:116
void resolve(std::map< std::string, size_t > datamap, size_t &datalen, size_t bindID, std::map< std::string, size_t > &argmap)
Definition: daFunk.hpp:1515
shared_ptr< FunkIntegrate_gsl1d > set_singularity_factor(double f)
Definition: daFunk.hpp:1571
ArgsType joinArgs(ArgsType args1, ArgsType args2)
Definition: daFunk.hpp:228
EXPORT_SYMBOLS warning & utils_warning()
Utility warnings.
dictionary fin
PlainPtrs1 plain(std::string)
Definition: daFunk.hpp:1087
static double plain4p(double x1, double x2, double x3, double x4, void *ptr)
Definition: daFunk.hpp:536
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f, std::string x)
Definition: daFunk.hpp:1502
std::pair< double(*)(double, double, double, void *), void * > PlainPtrs3
Definition: daFunk.hpp:131
void setup(Funk f0, std::string arg, Funk f1, Funk f2)
Definition: daFunk.hpp:1658
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:1365
FunkIntegrate_gsl1d(Funk f0, std::string arg, double x, double y)
Definition: daFunk.hpp:1486
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:1389
std::string my_arg
Definition: daFunk.hpp:675
std::vector< double > logspace(double x0, double x1, unsigned int n)
Definition: daFunk.hpp:186
std::vector< std::pair< Funk, Funk > > my_singularities
Definition: daFunk.hpp:1690
Funk func(double(*f)(funcargs...), Args... args)
Definition: daFunk.hpp:768
FunkIntegrate_gsl1d(Funk f0, std::string arg, Funk f1, Funk f2)
Definition: daFunk.hpp:1478
void setup(Model &mssm)
shared_ptr< FunkBound > f
Definition: daFunk.hpp:576
std::vector< double * > map
Definition: daFunk.hpp:831
virtual ~FunkBase()
Definition: daFunk.hpp:302
hb_ModelParameters void
std::vector< size_t > index
Definition: daFunk.hpp:1698
Funk operator-(Funk f)
Definition: daFunk.hpp:1153
START_MODEL dNur_CMB r
FunkIntegrate_gsl1d(Funk f0, std::string arg, std::string x, double y)
Definition: daFunk.hpp:1506
FunkDelta(std::string arg, double pos, double width)
Definition: daFunk.hpp:886
START_MODEL x3
Definition: demo.hpp:42
void setup(Funk f, std::string arg, Funk g)
Definition: daFunk.hpp:678
bool assert_args(Args... args)
Definition: daFunk.hpp:989
double(* PlainPtr4)(double &, double &, double &, double &)
Definition: daFunk.hpp:125
std::vector< std::vector< size_t > > indices
Definition: daFunk.hpp:370
Singularities singularities
Definition: daFunk.hpp:372
FunkConst(double c)
Definition: daFunk.hpp:589
static void bindID_manager(size_t &bindID, bool bind)
Definition: daFunk.hpp:417
std::vector< double > vect2(std::vector< std::vector< double >> &coll, double x, Args... argss)
Definition: daFunk.hpp:482
FunkDerived(Funk f, std::string arg, Funk g)
Definition: daFunk.hpp:616
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:706
double(* PlainPtr1)(double &)
Definition: daFunk.hpp:122
Funk cnst(double x, Args... argss)
Definition: daFunk.hpp:606
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)
Definition: daFunk.hpp:1498
void request(std::string message)
Request an exception.
Definition: exceptions.cpp:471
std::vector< double > augmentSingl(const std::vector< double > &xgrid, Funk f, int N=100, double sigma=3)
Definition: daFunk.hpp:1714
FunkPlain(Funk fin, std::string arg1, std::string arg2)
Definition: daFunk.hpp:517
double value(std::vector< double > &map, size_t bindID)
Definition: daFunk.hpp:397
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:1576
invalid_point_exception & invalid_point()
Invalid point exceptions.
void digest_input()
Definition: daFunk.hpp:764
Funk zero(Args... argss)
Definition: daFunk.hpp:604
static double plain1p(double x1, void *ptr)
Definition: daFunk.hpp:521
std::pair< double(*)(double, double, void *), void * > PlainPtrs2
Definition: daFunk.hpp:130
static double plain3(double &x1, double &x2, double &x3)
Definition: daFunk.hpp:555
std::vector< T > vec(std::vector< T > vector)
Definition: daFunk.hpp:142
double(* PlainPtr2)(double &, double &)
Definition: daFunk.hpp:123
FunkVar(std::string arg)
Definition: daFunk.hpp:911
size_t datalen
Definition: daFunk.hpp:371
Funk ifelse(Funk f, Funk g, Funk h)
Definition: daFunk.hpp:1373
static double plain4(double &x1, double &x2, double &x3, double &x4)
Definition: daFunk.hpp:561
shared_ptr< FunkIntegrate_gsl1d > set_epsabs(double epsabs)
Definition: daFunk.hpp:1567
shared_ptr< FunkBase > Funk
Definition: daFunk.hpp:113
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:796
std::vector< size_t > my_index
Definition: daFunk.hpp:686
std::vector< std::string > ArgsType
Definition: daFunk.hpp:117
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:591
FunkPlain(Funk fin, std::string arg1)
Definition: daFunk.hpp:516
shared_ptr< FunkIntegrate_gsl1d > gsl_integration(Args... args)
Definition: daFunk.hpp:1065
std::string arg4
Definition: daFunk.hpp:577
static double plain1(double &x1)
Definition: daFunk.hpp:543
double value(const std::vector< double > &data, size_t bindID)
Definition: daFunk.hpp:916
Singularities joinSingl(Singularities s1, Singularities s2)
Definition: daFunk.hpp:243
Funk one(Args... argss)
Definition: daFunk.hpp:602
static double invoke(double x, void *params)
Definition: daFunk.hpp:1681
double ppp(index_list< Args... >, std::tuple< typename std::remove_reference< funcargs >::type... > &my_input)
Definition: daFunk.hpp:824