Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
c2_function.hh
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
38 //
39 // $Id: c2_function.hh 105199 2017-07-14 12:18:09Z gcosmo $
40 
41 #ifndef __has_c2_function_hh
42 #define __has_c2_function_hh 1
43 
44 // MSVC does not automatically define numerical constants such as M_PI
45 // this came from the msdn website, so it should be right...
46 #ifdef _MSC_VER
47 #define _USE_MATH_DEFINES
48 #define c2_isnan _isnan
49 #define c2_isfinite _finite
50 #else
51 #define c2_isnan std::isnan
52 #define c2_isfinite std::isfinite
53 #endif
54 
55 #include <cmath>
56 #include <vector>
57 #include <utility>
58 #include <string>
59 #include <stdexcept>
60 #include <typeinfo>
61 #include <sstream>
62 #include <limits> // fails under gcc-4.3 without this here, was ok
63 
65 class c2_exception : public std::exception {
66 public:
69  c2_exception(const char msgcode[]) : info(msgcode) { }
70  virtual ~c2_exception() throw() { }
73  virtual const char* what() const throw() { return info.c_str(); }
74 private:
75  std::string info;
76 };
77 
78 // put these forward references here, and with a bogus typename to make swig
79 template <typename float_type> class c2_composed_function_p;
80 template <typename float_type> class c2_sum_p;
81 template <typename float_type> class c2_diff_p;
82 template <typename float_type> class c2_product_p;
83 template <typename float_type> class c2_ratio_p;
84 template <typename float_type> class c2_piecewise_function_p;
85 template <typename float_type> class c2_quadratic_p;
86 template <typename float_type> class c2_ptr;
101 template <typename float_type> class c2_fblock
105 {
106 public:
108  float_type x;
110  float_type y;
112  float_type yp;
114  float_type ypp;
117  bool ypbad;
119  bool yppbad;
120 };
121 
147 template <typename float_type=double> class c2_function {
148 public:
151  const std::string cvs_header_vers() const { return
152  "c2_function.hh 490 2012-04-10 19:05:40Z marcus ";
153  }
154 
157  const std::string cvs_file_vers() const ;
158 
159 public:
161  virtual ~c2_function() {
163  if(root_info) delete root_info;
164  if(owner_count) {
165  std::ostringstream outstr;
166  outstr << "attempt to delete an object with non-zero ownership in class";
167  outstr << typeid(*this).name() << std::endl;
168  //throw c2_exception(outstr.str().c_str());
169  }
170  }
171 
180  virtual float_type value_with_derivatives(float_type x, float_type *yprime,
181  float_type *yprime2) const /* throw(c2_exception) */ =0 ;
182 
186  inline float_type operator () (float_type x) const /* throw(c2_exception) */
187  { return value_with_derivatives(x, (float_type *)0, (float_type *)0); }
188 
195  inline float_type operator () (float_type x, float_type *yprime,
196  float_type *yprime2) const /* throw(c2_exception) */
197  { return value_with_derivatives(x, yprime, yprime2); }
198 
231  float_type find_root(float_type lower_bracket, float_type upper_bracket,
232  float_type start,
233  float_type value, int *error=0,
234  float_type *final_yprime=0,
235  float_type *final_yprime2=0 ) const /* throw(c2_exception) */;
267  float_type partial_integrals(std::vector<float_type> xgrid,
268  std::vector<float_type> *partials = 0,
269  float_type abs_tol=1e-12,
270  float_type rel_tol=1e-12,
271  int derivs=2, bool adapt=true,
272  bool extrapolate=true)
273  const /* throw(c2_exception) */;
274 
306  float_type integral(float_type amin, float_type amax,
307  std::vector<float_type> *partials = 0,
308  float_type abs_tol=1e-12, float_type rel_tol=1e-12,
309  int derivs=2, bool adapt=true, bool extrapolate=true)
310  const /* throw(c2_exception) */;
311 
384  float_type amin, float_type amax,
385  float_type abs_tol=1e-12, float_type rel_tol=1e-12,
386  int derivs=2, std::vector<float_type> *xvals=0,
387  std::vector<float_type> *yvals=0) const /* throw(c2_exception) */;
388 
389  inline float_type xmin() const { return fXMin; }
390  inline float_type xmax() const { return fXMax; }
391  void set_domain(float_type amin, float_type amax) { fXMin=amin; fXMax=amax; }
392 
395  size_t get_evaluations() const { return evaluations; }
397  void reset_evaluations() const { evaluations=0; }
399  inline void increment_evaluations() const { evaluations++; }
400 
408  bool check_monotonicity(const std::vector<float_type> &data,
409  const char message[]) const /* throw(c2_exception) */;
410 
427  virtual void set_sampling_grid(const std::vector<float_type> &grid)
428  /* throw(c2_exception) */;
429 
432  std::vector<float_type> *get_sampling_grid_pointer() const {
433  return sampling_grid; }
434 
435  virtual void get_sampling_grid(float_type amin, float_type amax,
436  std::vector<float_type> &grid) const ;
437 
439  void preen_sampling_grid(std::vector<float_type> *result) const;
440  void refine_sampling_grid(std::vector<float_type> &grid,
441  size_t refinement) const;
442 
450  float_type amin,
451  float_type amax,
452  float_type norm=1.0) const /* throw(c2_exception) */;
454  float_type amin, float_type amax,
455  float_type norm=1.0)
456  const /* throw(c2_exception) */;
465  float_type amin, float_type amax, const c2_function<float_type> &weight,
466  float_type norm=1.0)
467  const /* throw(c2_exception) */;
468 
473  c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs) const
474  { return *new c2_sum_p<float_type>(*this, rhs); }
480  { return *new c2_diff_p<float_type>(*this, rhs); }
485  c2_product_p<float_type> &operator *
486  (const c2_function<float_type> &rhs) const
487  { return *new c2_product_p<float_type>(*this, rhs); }
489  { return *new c2_ratio_p<float_type>(*this, rhs); }
495  (const c2_function<float_type> &inner) const
496  { return *new c2_composed_function_p<float_type>((*this), inner); }
497 
502  float_type get_trouble_point() const { return bad_x_point; }
503 
506  void claim_ownership() const { owner_count++; }
509  size_t release_ownership_for_return() const /* throw(c2_exception) */ {
510  if(!owner_count) {
511  std::ostringstream outstr;
512  outstr << "attempt to release ownership of an unowned function in class ";
513  outstr << typeid(*this).name() << std::endl;
514  throw c2_exception(outstr.str().c_str());
515  }
516  owner_count--;
517  return owner_count;
518  }
519  void release_ownership() const /* throw(c2_exception) */ {
520  if(!release_ownership_for_return()) delete this;
521  }
524  size_t count_owners() const { return owner_count; }
525 
526 protected:
528  : sampling_grid(0),
529  no_overwrite_grid(false),
530  fXMin(src.fXMin), fXMax(src.fXMax), root_info(0), owner_count(0)
531  {} // copy constructor only copies domain, and is only for internal use
534  fXMin(-std::numeric_limits<float_type>::max()),
535  fXMax(std::numeric_limits<float_type>::max()), root_info(0), owner_count(0)
536  {}
537 
538  virtual void set_sampling_grid_pointer(std::vector<float_type> &grid)
539  {
542  }
543 
544  std::vector<float_type> * sampling_grid;
546 
547  float_type fXMin, fXMax;
548  mutable size_t evaluations;
551  mutable float_type bad_x_point;
552 public:
556  inline void fill_fblock(c2_fblock<float_type> &fb) const /* throw(c2_exception) */
557  {
558  fb.y=value_with_derivatives(fb.x, &fb.yp, &fb.ypp);
559  fb.ypbad=c2_isnan(fb.yp) || !c2_isfinite(fb.yp);
560  fb.yppbad=c2_isnan(fb.ypp) || !c2_isfinite(fb.ypp);
562  }
563 
564 private:
567  struct recur_item {
568  c2_fblock<float_type> f1; size_t depth;
569  float_type previous_estimate, abs_tol, step_sum;
570  bool done;
571  size_t f0index, f2index;
572  };
573 
574 
579  struct c2_integrate_recur {
580  c2_fblock<float_type> *f0, *f1;
581  float_type abs_tol, rel_tol, eps_scale, extrap_coef, extrap2,
582  dx_tolerance, abs_tol_min;
583  std::vector< recur_item > *rb_stack;
584  int derivs;
585  bool adapt, extrapolate, inited;
586  };
587 
590  struct c2_sample_recur {
591  c2_fblock<float_type> *f0, *f1;
592  float_type abs_tol, rel_tol, dx_tolerance, abs_tol_min;
593  int derivs;
595  std::vector<float_type> *xvals, *yvals;
596  std::vector< recur_item > *rb_stack;
597  bool inited;
598  };
599 
602  struct c2_root_info {
603  c2_fblock<float_type> lower, upper;
604  bool inited;
605  };
606 
612  float_type integrate_step(struct c2_integrate_recur &rb)
613  const /* throw(c2_exception) */;
614 
620  void sample_step(struct c2_sample_recur &rb) const /* throw(c2_exception) */;
621 
628  mutable struct c2_root_info *root_info;
629 
630  mutable size_t owner_count;
631 };
632 
639 template <typename float_type=double> class c2_classic_function_p
640  : public c2_function<float_type> {
641 public:
644  c2_classic_function_p(const float_type (*c_func)(float_type))
645  : c2_function<float_type>(), func(c_func) {}
646 
649  virtual float_type value_with_derivatives(
650  float_type x, float_type *yprime,
651  float_type *yprime2) const /* throw(c2_exception) */
652  {
653  if(!func)
654  throw c2_exception("c2_classic_function called with null function");
655  if(yprime) *yprime=0;
656  if(yprime2) *yprime2=0;
657  return func(x);
658  }
660 
661 protected:
663  const float_type (*func)(float_type);
664 };
665 
681 template <typename float_type> class c2_const_ptr {
682 public:
684  c2_const_ptr() : func(0) {}
688  { this->set_function(&f); }
692  { this->set_function(src.get_ptr()); }
694  {
695  if(func) func->release_ownership();
696  func=f;
697  if(func) func->claim_ownership();
698  }
699 
702  const c2_const_ptr<float_type> & operator =
704  { this->set_function(f.get_ptr()); return f; }
707  const c2_function<float_type> & operator =
709  { this->set_function(&f); return f; }
720  void release_for_return() /* throw(c2_exception) */
721  {
722  if(func) func->release_ownership_for_return();
723  func=0;
724  }
730  void unset_function(void) { this->set_function(0); }
732  ~c2_const_ptr() { this->set_function(0); }
733 
735  inline const c2_function<float_type> &get() const /* throw(c2_exception) */
736  {
737  if(!func) throw c2_exception("c2_ptr accessed uninitialized");
738  return *func;
739  }
741  inline const c2_function<float_type> *get_ptr() const { return func; }
743  inline const c2_function<float_type> *operator -> () const
744  { return &get(); }
746  bool valid() const { return func != 0; }
747 
750  operator const c2_function<float_type>& () const { return this->get(); }
751 
755 
756  float_type operator()(float_type x) const /* throw(c2_exception) */
757  { return get()(x); }
767  float_type operator()(float_type x, float_type *yprime,
768  float_type *yprime2) const /* throw(c2_exception) */
769  { return get().value_with_derivatives(x, yprime, yprime2); }
775  const /* throw(c2_exception) */
776  { return *new c2_sum_p<float_type>(get(), rhs); }
778  const /* throw(c2_exception) */
779  { return *new c2_diff_p<float_type>(get(), rhs); }
781  const /* throw(c2_exception) */
782  { return *new c2_product_p<float_type>(get(), rhs); }
784  const /* throw(c2_exception) */
785  { return *new c2_ratio_p<float_type>(get(), rhs); }
791  const /* throw(c2_exception) */
792  { return *new c2_composed_function_p<float_type>(get(), inner); }
793 
794 protected:
796 };
797 
798 template <typename float_type> class c2_ptr : public c2_const_ptr<float_type >
799 {
800 public:
802  c2_ptr() : c2_const_ptr<float_type>() {}
806  c2_const_ptr<float_type>() { this->set_function(&f); }
809  c2_ptr(const c2_ptr<float_type> &src) :
810  c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
812  inline c2_function<float_type> &get() const /* throw(c2_exception) */
813  { return *const_cast<c2_function<float_type>*>(
817  { return const_cast<c2_function<float_type>*>(this->func); }
820  { return &get(); }
824  { this->set_function(f.get_ptr()); return f; }
828  { this->set_function(&f); return f; }
829 private:
831  void operator =(const c2_const_ptr<float_type> &) { }
833  void operator =(const c2_function<float_type> &) { }
834 };
835 
836 template <typename float_type, template <typename> class c2_class >
838  : public c2_const_ptr<float_type> {
839 public:
841  c2_typed_ptr() : c2_ptr<float_type>() {}
844  c2_typed_ptr(c2_class<float_type> &f)
845  : c2_const_ptr<float_type>() { this->set_function(&f); }
849  : c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
850 
852  inline c2_class<float_type> &get() const /* throw(c2_exception) */
853  {
854  return *static_cast<c2_class<float_type> *>
856  }
858  inline c2_class<float_type> *operator -> () const
859  { return &get(); }
861  inline c2_class<float_type> *get_ptr() const
862  { return static_cast<c2_class<float_type> *>(
863  const_cast<c2_function<float_type>*>(this->func)); }
864 
865  operator c2_class<float_type>&() const { return get(); }
869  { this->set_function(f.get_ptr()); }
872  void operator =(c2_class<float_type> &f)
873  { this->set_function(&f); }
874 private:
875  void operator =(const c2_const_ptr<float_type> &) { }
876  void operator =(const c2_function<float_type> &) { }
877 };
878 
879 template <typename float_type=double> class c2_plugin_function_p :
880  public c2_function<float_type> {
881 public:
883  c2_plugin_function_p() : c2_function<float_type>(), func() {}
886  c2_function<float_type>(),func(f) { }
887 
889  {
890  func.set_function(f);
891  if(f) this->set_domain(f->xmin(), f->xmax());
892  }
895  virtual float_type value_with_derivatives(
896  float_type x, float_type *yprime,
897  float_type *yprime2) const /* throw(c2_exception) */
898  {
899  if(!func.valid())
900  throw c2_exception("c2_plugin_function_p called uninitialized");
901  return func->value_with_derivatives(x, yprime, yprime2);
902  }
904  virtual ~c2_plugin_function_p() { }
905 
907  void unset_function() { func.unset_function(); }
908 
909  virtual void get_sampling_grid(float_type amin, float_type amax,
910  std::vector<float_type> &grid) const
911  {
912  if(!func.valid())
913  throw c2_exception("c2_plugin_function_p called uninitialized");
914  if(this->sampling_grid)
916  else func->get_sampling_grid(amin, amax, grid);
917  }
918 protected:
920 };
921 
922 template <typename float_type=double> class c2_const_plugin_function_p
923  : public c2_plugin_function_p<float_type> {
924 public:
929  c2_plugin_function_p<float_type>() { this->set_function(&f); }
932  const_cast<c2_function<float_type>*>(f)); }
935 
937  const c2_function<float_type> &get() const /* throw(c2_exception) */
938  { return this->func.get(); }
939 };
940 
941 template <typename float_type=double> class c2_binary_function
942  : public c2_function<float_type> {
943 public:
944  virtual float_type value_with_derivatives(
945  float_type x, float_type *yprime,
946  float_type *yprime2) const /* throw(c2_exception) */
947  {
948  if(stub)
949  throw c2_exception("attempt to evaluate a c2_binary_function stub");
950  return this->combine(*Left.get_ptr(), *Right.get_ptr(), x, yprime, yprime2);
951  }
952 
955  virtual ~c2_binary_function() { }
956 
957 protected:
958  c2_binary_function( float_type (*combiner)(
961  float_type x, float_type *yprime,
962  float_type *yprime2),
963  const c2_function<float_type> &left,
964  const c2_function<float_type> &right) :
965  c2_function<float_type>(), combine(combiner), Left(left),
966  Right(right), stub(false)
967  {
968  this->set_domain(
969  (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(),
970  (left.xmax() < right.xmax()) ? left.xmax() : right.xmax()
971  );
972  }
973  c2_binary_function(float_type (*combiner)(
976  float_type x, float_type *yprime, float_type *yprime2)
977  ) : c2_function<float_type>(), combine(combiner),
978  Left(), Right(), stub(true) { }
979 
980 public:
981  float_type (* const combine)(
984  float_type x, float_type *yprime, float_type *yprime2);
985 
986 protected:
988  bool stub;
989 };
990 
991 template <typename float_type=double> class c2_scaled_function_p
992  : public c2_function<float_type> {
993 public:
999  float_type scale) :
1000  c2_function<float_type>(), func(outer), yscale(scale) { }
1001 
1004  void reset(float_type scale) { yscale=scale; }
1005 
1006  virtual float_type value_with_derivatives(
1007  float_type x, float_type *yprime,
1008  float_type *yprime2) const /* throw (c2_exception) */
1009  {
1010  float_type y=this->func->value_with_derivatives(x, yprime, yprime2);
1011  if(yprime) (*yprime)*=yscale;
1012  if(yprime2) (*yprime2)*=yscale;
1013  return y*yscale;
1014  }
1015 
1016 protected:
1020  float_type yscale;
1021 };
1022 
1032 template <typename float_type=double> class c2_cached_function_p
1033  : public c2_function<float_type> {
1034 public:
1039  : c2_function<float_type>(),
1040  func(f), init(false) {}
1043  virtual float_type value_with_derivatives(
1044  float_type x, float_type *yprime,
1045  float_type *yprime2) const /* throw(c2_exception) */
1046  {
1047  if(!init || x != x0) {
1048  y=this->func->value_with_derivatives(x, &yp, &ypp);
1049  x0=x;
1050  init=true;
1051  }
1052  if(yprime) *yprime=yp;
1053  if(yprime2) *yprime2=ypp;
1054  return y;
1055  }
1056 
1057 protected:
1060  mutable bool init;
1061  mutable float_type x0, y, yp, ypp;
1062 
1063 };
1064 
1065 template <typename float_type=double> class c2_composed_function_p
1066  : public c2_binary_function<float_type> {
1067 public:
1069  const c2_function<float_type> &inner) :
1070  c2_binary_function<float_type>(combine, outer, inner) {
1071  this->set_domain(inner.xmin(), inner.xmax()); }
1074 
1076  static float_type combine(const c2_function<float_type> &left,
1078  float_type x, float_type *yprime,
1079  float_type *yprime2) /* throw(c2_exception) */
1080  {
1081  float_type y0, y1;
1082  if(yprime || yprime2) {
1083  float_type yp0, ypp0, yp1, ypp1;
1084  y0=right.value_with_derivatives(x, &yp0, &ypp0);
1085  y1=left.value_with_derivatives(y0, &yp1, &ypp1);
1086  if(yprime) *yprime=yp1*yp0;
1087  if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1;
1088  } else {
1089  y0=right(x);
1090  y1=left(y0);
1091  }
1092  return y1;
1093  }
1094 };
1095 
1096 template <typename float_type=double> class c2_sum_p
1097  : public c2_binary_function<float_type> {
1098 public:
1104  : c2_binary_function<float_type>(combine, left, right) {}
1106  c2_sum_p() : c2_binary_function<float_type>(combine) {} ;
1107 
1109  static float_type combine(const c2_function<float_type> &left,
1111  float_type x, float_type *yprime,
1112  float_type *yprime2) /* throw(c2_exception) */
1113  {
1114  float_type y0, y1;
1115  if(yprime || yprime2) {
1116  float_type yp0, ypp0, yp1, ypp1;
1117  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1118  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1119  if(yprime) *yprime=yp0+yp1;
1120  if(yprime2) *yprime2=ypp0+ypp1;
1121  } else {
1122  y0=left(x);
1123  y1=right(x);
1124  }
1125  return y0+y1;
1126  }
1127 };
1128 
1129 template <typename float_type=double> class c2_diff_p
1130  : public c2_binary_function<float_type> {
1131 public:
1137  : c2_binary_function<float_type>(combine, left, right) {}
1139  c2_diff_p() : c2_binary_function<float_type>(combine) {} ;
1140 
1142  static float_type combine(const c2_function<float_type> &left,
1144  float_type x, float_type *yprime,
1145  float_type *yprime2) /* throw(c2_exception) */
1146  {
1147  float_type y0, y1;
1148  if(yprime || yprime2) {
1149  float_type yp0, ypp0, yp1, ypp1;
1150  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1151  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1152  if(yprime) *yprime=yp0-yp1;
1153  if(yprime2) *yprime2=ypp0-ypp1;
1154  } else {
1155  y0=left(x);
1156  y1=right(x);
1157  }
1158  return y0-y1;
1159  }
1160 };
1161 
1162 
1166 template <typename float_type=double> class c2_product_p
1167  : public c2_binary_function<float_type> {
1168 public:
1174  : c2_binary_function<float_type>(combine, left, right) {}
1176  c2_product_p() : c2_binary_function<float_type>(combine) {} ;
1177 
1179  static float_type combine(const c2_function<float_type> &left,
1181  float_type x, float_type *yprime,
1182  float_type *yprime2) /* throw(c2_exception) */
1183  {
1184  float_type y0, y1;
1185  if(yprime || yprime2) {
1186  float_type yp0, ypp0, yp1, ypp1;
1187  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1188  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1189  if(yprime) *yprime=y1*yp0+y0*yp1;
1190  if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0;
1191  } else {
1192  y0=left(x);
1193  y1=right(x);
1194  }
1195  return y0*y1;
1196  }
1197 };
1198 
1199 
1203 template <typename float_type=double> class c2_ratio_p
1204  : public c2_binary_function<float_type> {
1205 public:
1211  : c2_binary_function<float_type>(combine, left, right) {}
1213  c2_ratio_p() : c2_binary_function<float_type>(combine) {} ;
1214 
1216  static float_type combine(const c2_function<float_type> &left,
1218  float_type x, float_type *yprime,
1219  float_type *yprime2) /* throw(c2_exception) */
1220  {
1221  float_type y0, y1;
1222  if(yprime || yprime2) {
1223  float_type yp0, ypp0, yp1, ypp1;
1224  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1225  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1226  if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio
1227  if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)
1228  /(y1*y1*y1);
1229  } else {
1230  y0=left(x);
1231  y1=right(x);
1232  }
1233  return y0/y1;
1234  }
1235 };
1236 
1241 template <typename float_type> class c2_constant_p
1242  : public c2_function<float_type> {
1243 public:
1244  c2_constant_p(float_type x) : c2_function<float_type>(), value(x) {}
1245  void reset(float_type val) { value=val; }
1246  virtual float_type value_with_derivatives(
1247  float_type, float_type *yprime,
1248  float_type *yprime2) const /* throw(c2_exception) */
1249  { if(yprime) *yprime=0; if(yprime2) *yprime2=0; return value; }
1250 
1251 private:
1252  float_type value;
1253 };
1254 
1257 template <typename float_type> class c2_transformation {
1258 public:
1265  c2_transformation(bool transformed,
1266  float_type (*xin)(float_type),
1267  float_type (*xinp)(float_type),
1268  float_type (*xinpp)(float_type),
1269  float_type (*xout)(float_type)
1270  ) :
1271  fTransformed(transformed), fHasStaticTransforms(true),
1272  pIn(xin), pInPrime(xinp), pInDPrime(xinpp), pOut(xout) { }
1273 
1277  c2_transformation(bool transformed) :
1278  fTransformed(transformed), fHasStaticTransforms(false),
1280  pOut(report_error) { }
1282  virtual ~c2_transformation() { }
1284  const bool fTransformed;
1288 
1294  float_type (* const pIn)(float_type);
1296  float_type (* const pInPrime)(float_type);
1298  float_type (* const pInDPrime)(float_type);
1300  float_type (* const pOut)(float_type);
1301 
1303  virtual float_type fIn(float_type x) const { return pIn(x); }
1305  virtual float_type fInPrime(float_type x) const { return pInPrime(x); }
1307  virtual float_type fInDPrime(float_type x) const { return pInDPrime(x); }
1309  virtual float_type fOut(float_type x) const { return pOut(x); }
1310 
1311 protected:
1313  static float_type report_error(float_type x) {
1314  throw c2_exception("use of improperly constructed axis transform");
1315  return x; }
1317  static float_type ident(float_type x) { return x; }
1319  static float_type one(float_type) { return 1; }
1321  static float_type zero(float_type) { return 0; }
1323  static float_type recip(float_type x) { return 1.0/x; }
1325  static float_type recip_prime(float_type x) { return -1/(x*x); }
1327  static float_type recip_prime2(float_type x) { return 2/(x*x*x); }
1328 };
1329 
1332 template <typename float_type> class c2_transformation_linear
1333  : public c2_transformation<float_type> {
1334 public:
1337  false, this->ident,
1338  this->one, this->zero,
1339  this->ident) { }
1342 };
1345 template <typename float_type> class c2_transformation_log
1346  : public c2_transformation<float_type> {
1347 public:
1350  true, std::log, this->recip,
1351  this->recip_prime, std::exp) { }
1354 };
1357 template <typename float_type> class c2_transformation_recip
1358  : public c2_transformation<float_type> {
1359 public:
1362  true, this->recip,
1363  this->recip_prime,
1364  this->recip_prime2, this->recip) { }
1367 };
1368 
1378 template <typename float_type>
1380 public:
1385  const c2_transformation<float_type> &xx,
1386  const c2_transformation<float_type> &yy) :
1387  isIdentity(!(xx.fTransformed || yy.fTransformed)), X(xx), Y(yy) { }
1389  virtual ~c2_function_transformation() { delete &X; delete &Y; }
1390  virtual float_type evaluate(
1391  float_type xraw,
1392  float_type y, float_type yp0, float_type ypp0,
1393  float_type *yprime, float_type *yprime2) const;
1394  const bool isIdentity;
1399 };
1400 
1404 template <typename float_type> class c2_lin_lin_function_transformation :
1405  public c2_function_transformation<float_type> {
1406 public:
1408  c2_function_transformation<float_type>(
1409  *new c2_transformation_linear<float_type>,
1410  *new c2_transformation_linear<float_type>
1411  ) { }
1413 };
1414 
1418 template <typename float_type> class c2_log_log_function_transformation :
1419  public c2_function_transformation<float_type> {
1420 public:
1422  c2_function_transformation<float_type>(
1423  *new c2_transformation_log<float_type>,
1424  *new c2_transformation_log<float_type>
1425  ) { }
1427 };
1428 
1432 template <typename float_type> class c2_lin_log_function_transformation :
1433  public c2_function_transformation<float_type> {
1434 public:
1436  c2_function_transformation<float_type>(
1437  *new c2_transformation_linear<float_type>,
1438  *new c2_transformation_log<float_type>
1439  ) { }
1441 };
1442 
1446 template <typename float_type> class c2_log_lin_function_transformation :
1447  public c2_function_transformation<float_type> {
1448 public:
1450  c2_function_transformation<float_type>(
1451  *new c2_transformation_log<float_type>,
1452  *new c2_transformation_linear<float_type>
1453  ) { }
1455 };
1456 
1461 template <typename float_type> class c2_arrhenius_function_transformation :
1462  public c2_function_transformation<float_type> {
1463 public:
1465  c2_function_transformation<float_type>(
1466  *new c2_transformation_recip<float_type>,
1467  *new c2_transformation_log<float_type>
1468  ) { }
1470 };
1471 
1513 template <typename float_type=double> class interpolating_function_p
1514  : public c2_function<float_type> {
1515 public:
1519  fTransform(*new c2_lin_lin_function_transformation<float_type>) { }
1520 
1525  transform)
1526  : c2_function<float_type>(),
1527  fTransform(transform) { }
1528 
1552  const std::vector<float_type> &x,
1553  const std::vector<float_type> &f,
1554  bool lowerSlopeNatural, float_type lowerSlope,
1555  bool upperSlopeNatural, float_type upperSlope, bool splined=true
1556  ) /* throw(c2_exception) */;
1557 
1574  std::vector<std::pair<float_type, float_type> > &data,
1575  bool lowerSlopeNatural, float_type lowerSlope,
1576  bool upperSlopeNatural, float_type upperSlope, bool splined=true
1577  ) /* throw(c2_exception) */;
1578 
1615  float_type amin, float_type amax,
1616  float_type abs_tol, float_type rel_tol,
1617  bool lowerSlopeNatural, float_type lowerSlope,
1618  bool upperSlopeNatural, float_type upperSlope
1619  ) /* throw(c2_exception) */;
1620 
1646  const std::vector<float_type> &bincenters,
1647  const c2_function<float_type> &binheights)
1648  /* throw(c2_exception) */;
1649 
1651  const std::vector<float_type> &bins,
1652  const std::vector<float_type> &binheights,
1653  bool splined=true)
1654  /* throw(c2_exception) */;
1655 
1656  virtual float_type value_with_derivatives(
1657  float_type x, float_type *yprime,
1658  float_type *yprime2) const /* throw(c2_exception) */;
1659 
1661  virtual ~interpolating_function_p() { delete &fTransform; }
1662 
1664  const /* throw(c2_exception) */
1665  { return *new interpolating_function_p<float_type>(); }
1666 
1667  void get_data(std::vector<float_type> &xvals,
1668  std::vector<float_type> &yvals) const /* throw() */ ;
1669 
1671  std::vector<float_type> &xvals,
1672  std::vector<float_type> &yvals,
1673  std::vector<float_type> &y2vals) const
1674  { xvals=X; yvals=F; y2vals=y2; }
1675 
1676  void set_lower_extrapolation(float_type bound);
1677  void set_upper_extrapolation(float_type bound);
1678 
1680  unary_operator(const c2_function<float_type> &source) const;
1681 
1684  const c2_binary_function<float_type> *combining_stub
1685  ) const;
1688  return binary_operator(rhs, new c2_sum_p<float_type>()); }
1691  return binary_operator(rhs, new c2_diff_p<float_type>()); }
1694  return binary_operator(rhs, new c2_product_p<float_type>()); }
1697  return binary_operator(rhs, new c2_ratio_p<float_type>()); }
1699  Xraw=rhs.Xraw; X=rhs.X; F=rhs.F; y2=rhs.y2;
1701  }
1702 
1704 
1705 protected:
1707  void spline(
1708  bool lowerSlopeNatural, float_type lowerSlope,
1709  bool upperSlopeNatural, float_type upperSlope
1710  ) /* throw(c2_exception) */;
1711 
1712  static bool comp_pair(std::pair<float_type,float_type> const &i,
1713  std::pair<float_type,float_type> const &j)
1714  {return i.first<j.first;}
1715 
1716  std::vector<float_type> Xraw, X, F, y2;
1719  mutable size_t lastKLow;
1720 };
1721 
1725 template <typename float_type=double> class log_lin_interpolating_function_p
1726  : public interpolating_function_p <float_type> {
1727 public:
1731  interpolating_function_p<float_type>(
1732  *new c2_log_lin_function_transformation<float_type>)
1733  { }
1735  const /* throw(c2_exception) */
1737 };
1738 
1739 
1744 template <typename float_type=double> class lin_log_interpolating_function_p
1745  : public interpolating_function_p <float_type> {
1746 public:
1750  : interpolating_function_p<float_type>(
1751  *new c2_lin_log_function_transformation<float_type>)
1752  { }
1754  const /* throw(c2_exception) */
1756 };
1757 
1758 
1764 template <typename float_type=double> class log_log_interpolating_function_p
1765  : public interpolating_function_p <float_type> {
1766 public:
1770  interpolating_function_p<float_type>(
1771  *new c2_log_log_function_transformation<float_type>)
1772  { }
1774  const /* throw(c2_exception) */
1776 };
1777 
1778 
1784 template <typename float_type=double> class arrhenius_interpolating_function_p
1785  : public interpolating_function_p <float_type> {
1786 public:
1790  : interpolating_function_p<float_type>(
1791  *new c2_arrhenius_function_transformation<float_type>)
1792  { }
1794  const /* throw(c2_exception) */
1796 };
1797 
1802 template <typename float_type=double> class c2_sin_p
1803  : public c2_function<float_type> {
1804 public:
1806  c2_sin_p() : c2_function<float_type>() {}
1807 
1808  virtual float_type value_with_derivatives(
1809  float_type x, float_type *yprime,
1810  float_type *yprime2) const /* throw(c2_exception) */
1811  { float_type q=std::sin(x);
1812  if(yprime) *yprime=std::cos(x);
1813  if(yprime2) *yprime2=-q;
1814  return q; }
1815 
1816  virtual void get_sampling_grid(float_type amin, float_type amax,
1817  std::vector<float_type> &grid) const;
1818 };
1819 
1824 template <typename float_type=double> class c2_cos_p
1825  : public c2_sin_p<float_type> {
1826 public:
1828  c2_cos_p() : c2_sin_p<float_type>() {}
1829 
1830  virtual float_type value_with_derivatives(
1831  float_type x, float_type *yprime,
1832  float_type *yprime2) const /* throw(c2_exception) */
1833  { float_type q=std::cos(x);
1834  if(yprime) *yprime=-std::sin(x);
1835  if(yprime2) *yprime2=-q;
1836  return q; }
1837 };
1838 
1843 template <typename float_type=double> class c2_tan_p
1844  : public c2_function<float_type> {
1845 public:
1847  c2_tan_p() : c2_function<float_type>() {}
1848 
1849  virtual float_type value_with_derivatives(
1850  float_type x, float_type *yprime,
1851  float_type *yprime2) const /* throw(c2_exception) */
1852  {
1853  float_type c=std::cos(x), ss=std::sin(x);
1854  float_type t=ss/c;
1855  float_type yp=1/(c*c);
1856  if(yprime) { *yprime=yp; }
1857  if(yprime2){*yprime2=2*t*yp; }
1858  return t;
1859  }
1860 };
1861 
1866 template <typename float_type=double> class c2_log_p
1867  : public c2_function<float_type> {
1868 public:
1870  c2_log_p() : c2_function<float_type>() {}
1871 
1872  virtual float_type value_with_derivatives(
1873  float_type x, float_type *yprime,
1874  float_type *yprime2) const /* throw(c2_exception) */
1875  { if(yprime) *yprime=1.0/x;
1876  if(yprime2) *yprime2=-1.0/(x*x);
1877  return std::log(x); }
1878 };
1879 
1884 template <typename float_type=double> class c2_exp_p
1885  : public c2_function<float_type> {
1886 public:
1888  c2_exp_p() : c2_function<float_type>() {}
1889 
1890  virtual float_type value_with_derivatives(
1891  float_type x, float_type *yprime,
1892  float_type *yprime2) const /* throw(c2_exception) */
1893  { float_type q=std::exp(x);
1894  if(yprime) *yprime=q;
1895  if(yprime2) *yprime2=q;
1896  return q; }
1897 };
1898 
1903 template <typename float_type=double> class c2_sqrt_p
1904  : public c2_function<float_type> {
1905 public:
1907  c2_sqrt_p() : c2_function<float_type>() {}
1908 
1909  virtual float_type value_with_derivatives(
1910  float_type x, float_type *yprime,
1911  float_type *yprime2) const /* throw(c2_exception) */
1912  { float_type q=std::sqrt(x);
1913  if(yprime) *yprime=0.5/q;
1914  if(yprime2) *yprime2=-0.25/(x*q);
1915  return q; }
1916 };
1917 
1922 template <typename float_type=double> class c2_recip_p
1923  : public c2_function<float_type> {
1924 public:
1926  c2_recip_p(float_type scale) : c2_function<float_type>(), rscale(scale) {}
1927 
1928  virtual float_type value_with_derivatives(
1929  float_type x, float_type *yprime,
1930  float_type *yprime2) const /* throw(c2_exception) */
1931  {
1932  float_type q=1.0/x;
1933  float_type y=rscale*q;
1934  if(yprime) *yprime=-y*q;
1935  if(yprime2) *yprime2=2*y*q*q;
1936  return y;
1937  }
1940  void reset(float_type scale) { rscale=scale; }
1941 private:
1942  float_type rscale;
1943 };
1944 
1949 template <typename float_type=double> class c2_identity_p
1950  : public c2_function<float_type> {
1951 public:
1953  c2_identity_p() : c2_function<float_type>() {}
1954 
1955  virtual float_type value_with_derivatives(
1956  float_type x, float_type *yprime,
1957  float_type *yprime2) const /* throw(c2_exception) */
1958  { if(yprime) *yprime=1.0; if(yprime2) *yprime2=0; return x; }
1959 };
1960 
1972 template <typename float_type=double> class c2_linear_p
1973  : public c2_function<float_type> {
1974 public:
1979  c2_linear_p(float_type x0, float_type y0, float_type slope) :
1980  c2_function<float_type>(), xint(x0), intercept(y0), m(slope) {}
1985  void reset(float_type x0, float_type y0, float_type slope)
1986  { xint=x0; intercept=y0; m=slope; }
1987  virtual float_type value_with_derivatives(
1988  float_type x, float_type *yprime,
1989  float_type *yprime2) const /* throw(c2_exception) */
1990  { if(yprime) *yprime=m;
1991  if(yprime2) *yprime2=0;
1992  return m*(x-xint)+intercept; }
1993 
1994 private:
1995  float_type xint, intercept, m;
1996 protected:
1998 };
1999 
2014 template <typename float_type=double> class c2_quadratic_p
2015  : public c2_function<float_type> {
2016 public:
2022  c2_quadratic_p(float_type x0, float_type y0,
2023  float_type xcoef, float_type x2coef) :
2024  c2_function<float_type>(), intercept(y0), center(x0), a(x2coef), b(xcoef) {}
2030  void reset(float_type x0, float_type y0, float_type xcoef,
2031  float_type x2coef) { intercept=y0; center=x0; a=x2coef; b=xcoef; }
2032  virtual float_type value_with_derivatives(
2033  float_type x, float_type *yprime,
2034  float_type *yprime2) const /* throw(c2_exception) */
2035  { float_type dx=x-center;
2036  if(yprime) *yprime=2*a*dx+b;
2037  if(yprime2) *yprime2=2*a;
2038  return a*dx*dx+b*dx+intercept; }
2039 
2040 private:
2041  float_type intercept, center, a, b;
2042 protected:
2044 };
2045 
2058 template <typename float_type=double> class c2_power_law_p
2059  : public c2_function<float_type> {
2060 public:
2064  c2_power_law_p(float_type scale, float_type power) :
2065  c2_function<float_type>(), a(scale), b(power) {}
2069  void reset(float_type scale, float_type power) { a=scale; b=power; }
2070  virtual float_type value_with_derivatives(
2071  float_type x, float_type *yprime,
2072  float_type *yprime2) const /* throw(c2_exception) */
2073  { float_type q=a*std::pow(x,b-2);
2074  if(yprime) *yprime=b*q*x;
2075  if(yprime2) *yprime2=b*(b-1)*q;
2076  return q*x*x; }
2077 
2078 private:
2079  float_type a, b;
2080 protected:
2082 };
2083 
2101 template <typename float_type=double> class c2_inverse_function_p
2102  : public c2_function<float_type> {
2103 public:
2107  virtual float_type value_with_derivatives(
2108  float_type x, float_type *yprime,
2109  float_type *yprime2) const /* throw(c2_exception) */;
2110 
2114  void set_start_hint(float_type hint) const { start_hint=hint; }
2115 
2116  virtual float_type get_start_hint(float_type x) const
2117  { return hinting_function.valid()? hinting_function(x) : start_hint; }
2118 
2145  { hinting_function.set_function(hint_func); }
2151  { hinting_function=hint_func; }
2152 
2153 protected:
2155  mutable float_type start_hint;
2158 };
2159 
2173 template <typename float_type=double> class accumulated_histogram
2174  : public interpolating_function_p <float_type> {
2175 public:
2185  accumulated_histogram(const std::vector<float_type>binedges,
2186  const std::vector<float_type> binheights,
2187  bool normalize=false,
2188  bool inverse_function=false, bool drop_zeros=true);
2189 
2190 };
2191 
2215 template <typename float_type=double> class c2_connector_function_p
2216  : public c2_function<float_type> {
2217 public:
2228  c2_connector_function_p(float_type x0, const c2_function<float_type> &f0,
2229  float_type x2, const c2_function<float_type> &f2,
2230  bool auto_center, float_type y1);
2247  float_type x0, float_type y0, float_type yp0, float_type ypp0,
2248  float_type x2, float_type y2, float_type yp2, float_type ypp2,
2249  bool auto_center, float_type y1);
2259  const c2_fblock<float_type> &fb0,
2260  const c2_fblock<float_type> &fb2,
2261  bool auto_center, float_type y1);
2262 
2264  virtual ~c2_connector_function_p();
2265  virtual float_type value_with_derivatives(
2266  float_type x, float_type *yprime,
2267  float_type *yprime2) const /* throw (c2_exception) */;
2268 protected:
2270  void init(
2271  const c2_fblock<float_type> &fb0,
2272  const c2_fblock<float_type> &fb2,
2273  bool auto_center, float_type y1);
2274 
2275  float_type fhinv, fy1, fa, fb, fc, fd, fe, ff;
2276 };
2277 
2303 template <typename float_type=double> class c2_piecewise_function_p
2304  : public c2_function<float_type> {
2305 public:
2309  virtual ~c2_piecewise_function_p();
2310  virtual float_type value_with_derivatives(
2311  float_type x, float_type *yprime,
2312  float_type *yprime2) const /* throw (c2_exception) */;
2325  void append_function(const c2_function<float_type> &func)
2326  /* throw (c2_exception) */;
2327 protected:
2328  std::vector<c2_const_ptr<float_type> > functions;
2329  mutable int lastKLow;
2330 };
2331 
2332 #include "c2_function.icc"
2333 
2334 #endif
G4double G4ParticleHPJENDLHEData::G4double result
const c2_function< float_type > * func
Definition: c2_function.hh:795
c2_recip_p(float_type scale)
constructor.
interpolating_function_p< float_type > & load_random_generator_bins(const std::vector< float_type > &bins, const std::vector< float_type > &binheights, bool splined=true)
static float_type zero(float_type)
utility function f(x)=0 useful in axis transforms
float_type operator()(float_type x) const
evaluate the function in the classic way, ignoring derivatives.
Definition: c2_function.hh:186
const std::string cvs_file_vers() const
get versioning information for the source file
void reset(float_type scale)
set a new scale factor
c2_typed_ptr(c2_class< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:844
const XML_Char XML_Encoding * info
Definition: expat.h:530
log_log_interpolating_function_p()
an empty log-log cubic-spline interpolating_function_p
c2_exp_p()
constructor.
virtual float_type fInDPrime(float_type x) const
virtual input X transform second derivative
c2_plugin_function_p(c2_function< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:885
virtual float_type fOut(float_type x) const
virtual output X transform
std::vector< c2_const_ptr< float_type > > functions
c2_linear_p(float_type x0, float_type y0, float_type slope)
Construct the operator f=y0 + slope * (x-x0)
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
#define inline
Definition: internal.h:71
void set_start_hint(float_type hint) const
give the function a hint as to where to look for its inverse
interpolating_function_p< float_type > & multiply_pointwise(const c2_function< float_type > &rhs) const
c2_function< float_type > * get_ptr() const
get an unchecked pointer to our owned function
Definition: c2_function.hh:816
c2_log_p()
constructor.
create a c2_function which is the ratio of two other c2_functions.This should always be constructed u...
Definition: c2_function.hh:83
virtual void get_sampling_grid(float_type amin, float_type amax, std::vector< float_type > &grid) const
void get_internal_data(std::vector< float_type > &xvals, std::vector< float_type > &yvals, std::vector< float_type > &y2vals) const
const bool fTransformed
flag to indicate if this transform is not the identity
float_type x
the abscissa
Definition: c2_function.hh:108
a transformation of a function in and out of Arrhenius (1/x vs. log(y)) space
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
Definition: c2_function.hh:944
a container into which any conventional c-style function can be dropped, to create a degenerate c2_fu...
Definition: c2_function.hh:639
const c2_const_ptr< float_type > Left
Definition: c2_function.hh:987
the identity transform
compute tan(x) with its derivatives.The factory function c2_factory::tan() creates *new c2_tan_p ...
interpolating_function_p< float_type > & add_pointwise(const c2_function< float_type > &rhs) const
float_type fXMax
Definition: c2_function.hh:547
float_type ypp
the second derivative at x
Definition: c2_function.hh:114
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
create the formal inverse function of another functionfor example, given a c2_function f ...
c2_ratio_p()
Create a stub just for the combiner to avoid statics.
c2_const_plugin_function_p(const c2_function< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:928
c2_const_ptr< float_type > sampler_function
log axis transform
void claim_ownership() const
increment our reference count. Destruction is only legal if the count is zero.
Definition: c2_function.hh:506
std::vector< float_type > y2
virtual ~c2_piecewise_function_p()
destructor
size_t get_evaluations() const
Definition: c2_function.hh:395
virtual float_type value_with_derivatives(float_type, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
create a linear mapping of another functionfor example, given a c2_function f
virtual ~c2_connector_function_p()
destructor
c2_ratio_p< float_type > & operator/(const c2_function< float_type > &rhs) const
Definition: c2_function.hh:488
const c2_const_ptr< float_type > Right
Definition: c2_function.hh:987
void reset(float_type x0, float_type y0, float_type xcoef, float_type x2coef)
Modify the coefficients after construction.
A spline with X in reciprocal space and Y transformed in log space.Most useful for thermodynamic type...
void preen_sampling_grid(std::vector< float_type > *result) const
The grid is modified in place.
interpolating_function_p< float_type > & divide_pointwise(const c2_function< float_type > &rhs) const
float_type(*const combine)(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2)
Definition: c2_function.hh:981
c2_function< float_type > & normalized_function(float_type amin, float_type amax, float_type norm=1.0) const
create a new c2_function from this one which is normalized on the interval
c2_diff_p< float_type > & operator-(const c2_function< float_type > &rhs) const
Definition: c2_function.hh:777
virtual const char * what() const
Definition: c2_function.hh:73
c2_function< float_type > & square_normalized_function(float_type amin, float_type amax, float_type norm=1.0) const
c2_cos_p()
constructor.
const c2_ptr< float_type > & operator=(const c2_ptr< float_type > &f)
fill the container from another container
Definition: c2_function.hh:823
interpolating_function_p< float_type > & subtract_pointwise(const c2_function< float_type > &rhs) const
const c2_transformation< float_type > & Y
the Y axis transform
c2_const_plugin_function_p()
construct the container with no function
Definition: c2_function.hh:926
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
static float_type recip(float_type x)
utility function f(x)=1/x useful in axis transforms
interpolating_function_p(const c2_function_transformation< float_type > &transform)
an empty cubic-spline interpolating_function_p with a specific transform
a transformation of a function in and out of log-log space
interpolating_function_p< float_type > & load(const std::vector< float_type > &x, const std::vector< float_type > &f, bool lowerSlopeNatural, float_type lowerSlope, bool upperSlopeNatural, float_type upperSlope, bool splined=true)
do the dirty work of constructing the spline from a function.
float_type get_trouble_point() const
Find out where a calculation ran into trouble, if it got a nan. If the most recent computation did no...
Definition: c2_function.hh:502
c2_diff_p< float_type > & operator-(const c2_function< float_type > &rhs) const
factory function to create a c2_diff_p from a regular algebraic expression.
Definition: c2_function.hh:479
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
virtual void get_sampling_grid(float_type amin, float_type amax, std::vector< float_type > &grid) const
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
Definition: c2_function.hh:895
c2_cached_function_p(const c2_function< float_type > &f)
construct the container
a transformation of a function in and out of lin-lin space
#define c2_isnan
Definition: c2_function.hh:51
virtual ~interpolating_function_p()
destructor
tuple x
Definition: test.py:50
void append_function(const c2_function< float_type > &func)
append a new function to the sequence
compute cos(x) with its derivatives.The factory function c2_factory::cos() creates *new c2_cos_p ...
c2_classic_function_p(const float_type(*c_func)(float_type))
construct the container
Definition: c2_function.hh:644
size_t evaluations
Definition: c2_function.hh:548
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
compute exp(x) with its derivatives.The factory function c2_factory::exp() creates *new c2_exp_p ...
create a c2_function which is the product of two other c2_functions.This should always be constructed...
Definition: c2_function.hh:82
virtual void get_sampling_grid(float_type amin, float_type amax, std::vector< float_type > &grid) const
Definition: c2_function.hh:909
const c2_function< float_type > * operator->() const
get a checked pointer to our owned function
Definition: c2_function.hh:743
static float_type ident(float_type x)
utility function f(x)=x useful in axis transforms
float_type(*const pOut)(float_type)
non-virtual pointer to output X transform
const G4ThreeVector const G4double const
void set_function(const c2_function< float_type > *f)
Definition: c2_function.hh:930
reciprocal axis transform
c2_class< float_type > * operator->() const
get a checked pointer to our owned function
Definition: c2_function.hh:858
c2_constant_p(float_type x)
float_type partial_integrals(std::vector< float_type > xgrid, std::vector< float_type > *partials=0, float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) const
const c2_transformation< float_type > & X
the X axis transform
virtual float_type evaluate(float_type xraw, float_type y, float_type yp0, float_type ypp0, float_type *yprime, float_type *yprime2) const
c2_sum_p()
Create a stub just for the combiner to avoid statics.
c2_piecewise_function_p()
construct the container
interpolating_function_p< float_type > & load_pairs(std::vector< std::pair< float_type, float_type > > &data, bool lowerSlopeNatural, float_type lowerSlope, bool upperSlopeNatural, float_type upperSlope, bool splined=true)
do the dirty work of constructing the spline from a function.
c2_product_p< float_type > & operator*(const c2_function< float_type > &rhs) const
Definition: c2_function.hh:780
std::vector< float_type > F
c2_plugin_function_p()
construct the container with no function
Definition: c2_function.hh:883
static float_type one(float_type)
utility function f(x)=1 useful in axis transforms
bool check_monotonicity(const std::vector< float_type > &data, const char message[]) const
check that a vector is monotonic, throw an exception if not, and return a flag if it is reversed ...
virtual ~c2_const_plugin_function_p()
destructor
Definition: c2_function.hh:934
virtual ~c2_transformation()
the destructor
c2_function(const c2_function< float_type > &src)
Definition: c2_function.hh:527
interpolating_function_p< float_type > & binary_operator(const c2_function< float_type > &rhs, const c2_binary_function< float_type > *combining_stub) const
std::vector< float_type > * get_sampling_grid_pointer() const
get the sampling grid, which may be a null pointer
Definition: c2_function.hh:432
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
virtual interpolating_function_p< float_type > & clone() const
void set_hinting_function(const c2_function< float_type > *hint_func)
set or unset the approximate function used to start the root finder
float_type(*const pInDPrime)(float_type)
non-virtual pointer to input X transform second derivative
virtual void set_sampling_grid(const std::vector< float_type > &grid)
establish a grid of &#39;interesting&#39; points on the function.
const XML_Char const XML_Char * data
Definition: expat.h:268
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
size_t release_ownership_for_return() const
decrement our reference count. Do not destroy at zero.
Definition: c2_function.hh:509
A spline with X and Y transformed into log space.Most useful for functions looking like y=x^n or any ...
c2_function_transformation(const c2_transformation< float_type > &xx, const c2_transformation< float_type > &yy)
construct this from two c2_transformation instances
compute scale/x with its derivatives.The factory function c2_factory::recip() creates *new c2_recip_p...
void spline(bool lowerSlopeNatural, float_type lowerSlope, bool upperSlopeNatural, float_type upperSlope)
create the spline coefficients
void reset(float_type scale)
reset the scale factor
static float_type combine(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2)
execute math necessary to do multiplication
const c2_function_transformation< float_type > & fTransform
create a cubic spline interpolation of a set of (x,y) pairsThis is one of the main reasons for c2_fun...
c2_composed_function_p()
Create a stub just for the combiner to avoid statics.
log_lin_interpolating_function_p()
an empty log-linear cubic-spline interpolating_function_p
bool valid() const
check if we have a valid function
Definition: c2_function.hh:746
A spline with X transformed into log space.
void set_hinting_function(const c2_const_ptr< float_type > hint_func)
set the hinting function from a pointer.
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
A container into which any other c2_function can be dropped.It allows a function to be pre-evaluated ...
c2_quadratic_p(float_type x0, float_type y0, float_type xcoef, float_type x2coef)
Construct the operator.
structure used to hold evaluated function data at a point.
Definition: c2_function.hh:104
void reset(float_type x0, float_type y0, float_type slope)
Change the slope and intercepts after construction.
c2_piecewise_function_p< float_type > * adaptively_sample(float_type amin, float_type amax, float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, std::vector< float_type > *xvals=0, std::vector< float_type > *yvals=0) const
create a c2_piecewise_function_p from c2_connector_function_p segments which is a representation of t...
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
c2_function< float_type > * operator->() const
get a checked pointer to our owned function
Definition: c2_function.hh:819
void set_function(c2_function< float_type > *f)
Definition: c2_function.hh:888
compute x with its derivatives.The factory function c2_factory::identity() creates *new c2_identity_p...
void set_upper_extrapolation(float_type bound)
c2_const_ptr(const c2_const_ptr< float_type > &src)
copy constructor
Definition: c2_function.hh:691
a transformation of a coordinate, including an inverse
float_type(*const pIn)(float_type)
non-virtual pointer to input X transform
c2_transformation(bool transformed, float_type(*xin)(float_type), float_type(*xinp)(float_type), float_type(*xinpp)(float_type), float_type(*xout)(float_type))
initialize all our function pointers
virtual ~c2_transformation_recip()
destructor
void set_domain(float_type amin, float_type amax)
Definition: c2_function.hh:391
virtual interpolating_function_p< float_type > & clone() const
c2_binary_function(float_type(*combiner)(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2), const c2_function< float_type > &left, const c2_function< float_type > &right)
Definition: c2_function.hh:958
virtual ~c2_function()
destructor
Definition: c2_function.hh:161
virtual ~c2_transformation_log()
destructor
a transformation of a function in and out of lin-log space
c2_tan_p()
constructor.
void release_for_return()
release the function without destroying it, so it can be returned from a function ...
Definition: c2_function.hh:720
std::vector< float_type > X
const c2_const_ptr< float_type > func
c2_diff_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left - right
const c2_function< float_type > & get() const
get a reference to our owned function
Definition: c2_function.hh:735
float_type yp
the derivative at x
Definition: c2_function.hh:112
static float_type combine(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2)
execute math necessary to do subtraction
create a container for a c2_function which handles the reference counting. It is useful as a smart co...
Definition: c2_function.hh:681
static float_type combine(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2)
execute math necessary to do division
c2_transformation(bool transformed)
initialize all our function pointers so that only the (overridden) virtual functions can be called wi...
c2_transformation_log()
constructor
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
c2_transformation_linear()
constructor
std::vector< float_type > Xraw
const c2_const_ptr< float_type > func
const float_type(* func)(float_type)
pointer to our function
Definition: c2_function.hh:663
interpolating_function_p< float_type > & sample_function(const c2_function< float_type > &func, float_type amin, float_type amax, float_type abs_tol, float_type rel_tol, bool lowerSlopeNatural, float_type lowerSlope, bool upperSlopeNatural, float_type upperSlope)
do the dirty work of constructing the spline from a function.
const std::string cvs_header_vers() const
get versioning information for the header file
Definition: c2_function.hh:151
the exception class for c2_function operations.
Definition: c2_function.hh:65
float_type integral(float_type amin, float_type amax, std::vector< float_type > *partials=0, float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true) const
a fully-automated integrator which uses the information provided by the get_sampling_grid() function ...
c2_ptr()
construct the container with no function
Definition: c2_function.hh:802
float_type find_root(float_type lower_bracket, float_type upper_bracket, float_type start, float_type value, int *error=0, float_type *final_yprime=0, float_type *final_yprime2=0) const
solve f(x)==value very efficiently, with explicit knowledge of derivatives of the function ...
const c2_function< float_type > * get_ptr() const
get an unchecked pointer to our owned function
Definition: c2_function.hh:741
create a c2_function which smoothly connects two other c2_functions.This takes two points and generat...
static float_type recip_prime(float_type x)
utility function f(x)=-1/x**2 useful in axis transforms
a transformation of a function in and out of a coordinate space, using 2 c2_transformations ...
c2_product_p()
Create a stub just for the combiner to avoid statics.
bool yppbad
is NaN of Inf
Definition: c2_function.hh:119
c2_transformation_recip()
constructor
the parent class for all c2_functions.c2_functions know their value, first, and second derivative at ...
Definition: c2_function.hh:147
arrhenius_interpolating_function_p()
an empty arrhenius cubic-spline interpolating_function_p
c2_sin_p()
constructor.
c2_product_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left * right
virtual ~c2_transformation_linear()
destructor
c2_typed_ptr(const c2_typed_ptr< float_type, c2_class > &src)
copy constructor
Definition: c2_function.hh:848
T max(const T t1, const T t2)
brief Return the largest of the two arguments
c2_ptr(c2_function< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:805
c2_const_ptr(const c2_function< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:687
c2_exception(const char msgcode[])
construct the exception with an error message
Definition: c2_function.hh:69
c2_ratio_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left / right
static float_type recip_prime2(float_type x)
utility function f(x)=2/x**3 useful in axis transforms
float_type fXMin
Definition: c2_function.hh:547
interpolating_function_p< float_type > & unary_operator(const c2_function< float_type > &source) const
c2_class< float_type > * get_ptr() const
get an unchecked pointer to our owned function
Definition: c2_function.hh:861
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
Definition: c2_function.hh:649
A spline with Y transformed into log space.Most useful for functions looking like y=exp(x) ...
c2_power_law_p(float_type scale, float_type power)
Construct the operator.
c2_composed_function_p(const c2_function< float_type > &outer, const c2_function< float_type > &inner)
float_type(*const pInPrime)(float_type)
non-virtual pointer to input X transform derivative
float_type operator()(float_type x, float_type *yprime, float_type *yprime2) const
convenience operator to make us look like a function
Definition: c2_function.hh:767
c2_binary_function(float_type(*combiner)(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2))
Definition: c2_function.hh:973
c2_scaled_function_p(const c2_function< float_type > &outer, float_type scale)
construct the function with its scale factor.
Definition: c2_function.hh:998
virtual interpolating_function_p< float_type > & clone() const
~c2_const_ptr()
destructor
Definition: c2_function.hh:732
c2_sum_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left + right
c2_ptr< float_type > func
Definition: c2_function.hh:919
virtual ~c2_function_transformation()
destructor
const c2_const_ptr< float_type > func
the scaling factor for the function
c2_sum_p< float_type > & operator+(const c2_function< float_type > &rhs) const
factory function to create a c2_sum_p from a regular algebraic expression.
Definition: c2_function.hh:774
c2_const_ptr()
construct the container with no function
Definition: c2_function.hh:684
An interpolating_function_p which is the cumulative integral of a histogram.Note than binedges should...
virtual float_type get_start_hint(float_type x) const
size_t count_owners() const
get the reference count, mostly for debugging
Definition: c2_function.hh:524
virtual interpolating_function_p< float_type > & clone() const
lin_log_interpolating_function_p()
an empty linear-log cubic-spline interpolating_function_p
float_type y
the value of the function at x
Definition: c2_function.hh:110
c2_connector_function_p(float_type x0, const c2_function< float_type > &f0, float_type x2, const c2_function< float_type > &f2, bool auto_center, float_type y1)
construct the container from two functions
virtual ~c2_exception()
Definition: c2_function.hh:70
void operator=(const c2_typed_ptr< float_type, c2_class > &f)
fill the container from another container
Definition: c2_function.hh:868
accumulated_histogram(const std::vector< float_type >binedges, const std::vector< float_type > binheights, bool normalize=false, bool inverse_function=false, bool drop_zeros=true)
Construct the integrated histogram.
static float_type combine(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2)
execute math necessary to do addition
void reset(float_type val)
virtual ~c2_plugin_function_p()
destructor
Definition: c2_function.hh:904
void refine_sampling_grid(std::vector< float_type > &grid, size_t refinement) const
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
static PROLOG_HANDLER error
Definition: xmlrole.cc:112
std::vector< float_type > * sampling_grid
Definition: c2_function.hh:544
compute sqrt(x) with its derivatives.The factory function c2_factory::sqrt() creates *new c2_sqrt_p()...
interpolating_function_p< float_type > & load_random_generator_function(const std::vector< float_type > &bincenters, const c2_function< float_type > &binheights)
initialize from a grid of points and a c2_function (un-normalized) to an interpolator which...
virtual ~c2_binary_function()
destructor releases ownership of member functions
Definition: c2_function.hh:955
c2_ratio_p< float_type > & operator/(const c2_function< float_type > &rhs) const
Definition: c2_function.hh:783
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
#define c2_isfinite
Definition: c2_function.hh:52
void increment_evaluations() const
count evaluations
Definition: c2_function.hh:399
float_type bad_x_point
this point may be used to record where a calculation ran into trouble
Definition: c2_function.hh:551
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const =0
get the value and derivatives.
bool no_overwrite_grid
Definition: c2_function.hh:545
tuple c
Definition: test.py:13
static float_type report_error(float_type x)
utility function for unimplemented conversion
void release_ownership() const
Definition: c2_function.hh:519
void clone_data(const interpolating_function_p< float_type > &rhs)
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
void fill_fblock(c2_fblock< float_type > &fb) const
fill in a c2_fblock&lt;float_type&gt;... a shortcut for the integrator &amp; sampler
Definition: c2_function.hh:556
c2_diff_p()
Create a stub just for the combiner to avoid statics.
float_type operator()(float_type x) const
convenience operator to make us look like a function
Definition: c2_function.hh:756
virtual interpolating_function_p< float_type > & clone() const
create a quadratic mapping of another functionfor example, given a c2_function f
Definition: c2_function.hh:85
float_type xmax() const
Definition: c2_function.hh:390
void reset_evaluations() const
reset the counter
Definition: c2_function.hh:397
void set_lower_extrapolation(float_type bound)
c2_identity_p()
constructor.
void reset(float_type scale, float_type power)
Modify the mapping after construction.
virtual float_type fIn(float_type x) const
virtual input X transform
compute sin(x) with its derivatives.The factory function c2_factory::sin() creates *new c2_sin_p ...
interpolating_function_p()
an empty linear-linear cubic-spline interpolating_function_p
void unset_function()
clear our function
Definition: c2_function.hh:907
c2_sqrt_p()
constructor.
c2_const_ptr< float_type > hinting_function
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
float_type xmin() const
Definition: c2_function.hh:389
virtual float_type fInPrime(float_type x) const
virtual input X transform derivative
void init(const c2_fblock< float_type > &fb0, const c2_fblock< float_type > &fb2, bool auto_center, float_type y1)
fill container numerically
virtual ~c2_classic_function_p()
Definition: c2_function.hh:659
static float_type combine(const c2_function< float_type > &left, const c2_function< float_type > &right, float_type x, float_type *yprime, float_type *yprime2)
execute math necessary to do composition
c2_typed_ptr()
construct the container with no function
Definition: c2_function.hh:841
c2_ptr(const c2_ptr< float_type > &src)
copy constructor
Definition: c2_function.hh:809
create a c2_function which is a piecewise assembly of other c2_functions.The functions must have incr...
Definition: c2_function.hh:84
create a power law mapping of another functionfor example, given a c2_function f
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
virtual void set_sampling_grid_pointer(std::vector< float_type > &grid)
Definition: c2_function.hh:538
void get_data(std::vector< float_type > &xvals, std::vector< float_type > &yvals) const
a transformation of a function in and out of log-lin space
a c2_function which is constantThe factory function c2_factory::constant() creates *new c2_constant_p...
static bool comp_pair(std::pair< float_type, float_type > const &i, std::pair< float_type, float_type > const &j)
const bool fHasStaticTransforms
flag to indicate if the static function pointers can be used for efficiency
void unset_function(void)
clear the function
Definition: c2_function.hh:730
compute log(x) with its derivatives.The factory function c2_factory::log() creates *new c2_log_p ...
void set_function(const c2_function< float_type > *f)
Definition: c2_function.hh:693