Geant4  10.02.p01
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 66587 2012-12-21 11:06:44Z ihrivnac $
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 without this.
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 in c2_function.cc before
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 happy.
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;
98 template <typename float_type> class c2_fblock
102 {
103 public:
105  float_type x;
107  float_type y;
109  float_type yp;
111  float_type ypp;
113  bool ypbad;
115  bool yppbad;
116 };
117 
138 template <typename float_type=double> class c2_function {
139 public:
142  const std::string cvs_header_vers() const { return
143  "c2_function.hh 490 2012-04-10 19:05:40Z marcus ";
144  }
145 
148  const std::string cvs_file_vers() const ;
149 
150 public:
152  virtual ~c2_function() {
154  if(root_info) delete root_info;
155  if(owner_count) {
156  std::ostringstream outstr;
157  outstr << "attempt to delete an object with non-zero ownership in class ";
158  outstr << typeid(*this).name() << std::endl;
159  throw c2_exception(outstr.str().c_str());
160  }
161  }
162 
171  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception) =0 ; // { return 0; };
172 
176  inline float_type operator () (float_type x) const throw(c2_exception)
177  { return value_with_derivatives(x, (float_type *)0, (float_type *)0); }
178 
185  inline float_type operator () (float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
186  { return value_with_derivatives(x, yprime, yprime2); }
187 
206  float_type find_root(float_type lower_bracket, float_type upper_bracket, float_type start,
207  float_type value, int *error=0,
208  float_type *final_yprime=0, float_type *final_yprime2=0 ) const throw(c2_exception) ; // solve f(x)=value
209 
228  float_type partial_integrals(std::vector<float_type> xgrid, std::vector<float_type> *partials = 0,
229  float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true)
230  const throw(c2_exception);
231 
251  float_type integral(float_type amin, float_type amax, std::vector<float_type> *partials = 0,
252  float_type abs_tol=1e-12, float_type rel_tol=1e-12, int derivs=2, bool adapt=true, bool extrapolate=true)
253  const throw(c2_exception);
254 
292  c2_piecewise_function_p<float_type> *adaptively_sample(float_type amin, float_type amax,
293  float_type abs_tol=1e-12, float_type rel_tol=1e-12,
294  int derivs=2, std::vector<float_type> *xvals=0, std::vector<float_type> *yvals=0) const throw(c2_exception);
295 
297  inline float_type xmin() const { return fXMin; }
299  inline float_type xmax() const { return fXMax; }
301  void set_domain(float_type amin, float_type amax) { fXMin=amin; fXMax=amax; }
302 
307  size_t get_evaluations() const { return evaluations; }
309  void reset_evaluations() const { evaluations=0; } // evaluations are 'invisible' to constant
311  inline void increment_evaluations() const { evaluations++; }
312 
318  bool check_monotonicity(const std::vector<float_type> &data, const char message[]) const throw(c2_exception);
319 
328  virtual void set_sampling_grid(const std::vector<float_type> &grid) throw(c2_exception);
329 
332  std::vector<float_type> *get_sampling_grid_pointer() const { return sampling_grid; }
333 
340  virtual void get_sampling_grid(float_type amin, float_type amax, std::vector<float_type> &grid) const ;
341 
345  void preen_sampling_grid(std::vector<float_type> *result) const;
349  void refine_sampling_grid(std::vector<float_type> &grid, size_t refinement) const;
350 
356  c2_function<float_type> &normalized_function(float_type amin, float_type amax, float_type norm=1.0) const throw(c2_exception);
362  c2_function<float_type> &square_normalized_function(float_type amin, float_type amax, float_type norm=1.0) const throw(c2_exception);
370  float_type amin, float_type amax, const c2_function<float_type> &weight, float_type norm=1.0)
371  const throw(c2_exception);
372 
376  c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs) const
377  { return *new c2_sum_p<float_type>(*this, rhs); }
382  { return *new c2_diff_p<float_type>(*this, rhs); }
387  { return *new c2_product_p<float_type>(*this, rhs); }
392  { return *new c2_ratio_p<float_type>(*this, rhs); }
398  { return *new c2_composed_function_p<float_type>((*this), inner); }
399 
403  float_type get_trouble_point() const { return bad_x_point; }
404 
406  void claim_ownership() const { owner_count++; }
410  if(!owner_count) {
411  std::ostringstream outstr;
412  outstr << "attempt to release ownership of an unowned function in class ";
413  outstr << typeid(*this).name() << std::endl;
414  throw c2_exception(outstr.str().c_str());
415  }
416  owner_count--;
417  return owner_count;
418  }
420  void release_ownership() const throw(c2_exception) {
421  if(!release_ownership_for_return()) delete this;
422  }
425  size_t count_owners() const { return owner_count; }
426 
427 protected:
429  no_overwrite_grid(false),
430  fXMin(src.fXMin), fXMax(src.fXMax), root_info(0), owner_count(0)
431  {} // copy constructor only copies domain, and is only for internal use
434  fXMin(-std::numeric_limits<float_type>::max()),
435  fXMax(std::numeric_limits<float_type>::max()), root_info(0), owner_count(0)
436  {} // prevent accidental naked construction (impossible any since this has pure virtual methods)
437 
438  // this should only be called very early on, by a constructor, before anyone else
439  // sets a sampling grid, or it will leak memory
440  virtual void set_sampling_grid_pointer(std::vector<float_type> &grid)
441  {
442  if (sampling_grid && !no_overwrite_grid) delete sampling_grid; // grid was ours, lose it.
444  }
445 
446  std::vector<float_type> * sampling_grid;
448 
449  float_type fXMin, fXMax;
450  mutable size_t evaluations;
452  mutable float_type bad_x_point;
453 public:
456  inline void fill_fblock(c2_fblock<float_type> &fb) const throw(c2_exception)
457  {
458  fb.y=value_with_derivatives(fb.x, &fb.yp, &fb.ypp);
459  fb.ypbad=c2_isnan(fb.yp) || !c2_isfinite(fb.yp);
460  fb.yppbad=c2_isnan(fb.ypp) || !c2_isfinite(fb.ypp);
462  }
463 
464 private:
466  struct recur_item {
469  bool done;
470  size_t f0index, f2index;
471  };
472 
473 
481  std::vector< recur_item > *rb_stack;
482  int derivs;
484  };
485 
488  struct c2_sample_recur {
491  int derivs;
493  std::vector<float_type> *xvals, *yvals;
494  std::vector< recur_item > *rb_stack;
495  bool inited;
496  };
497 
500  struct c2_root_info {
502  bool inited;
503  };
504 
510  float_type integrate_step(struct c2_integrate_recur &rb) const throw(c2_exception);
511 
517  void sample_step(struct c2_sample_recur &rb) const throw(c2_exception);
518 
523  mutable struct c2_root_info *root_info;
524 
525  mutable size_t owner_count;
526 };
527 
534 template <typename float_type=double> class c2_classic_function_p : public c2_function<float_type> {
535 public:
538  c2_classic_function_p(const float_type (*c_func)(float_type)) : c2_function<float_type>(), func(c_func) {}
539 
542  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
543  {
544  if(!func) throw c2_exception("c2_classic_function called with null function");
545  if(yprime) *yprime=0;
546  if(yprime2) *yprime2=0;
547  return func(x);
548  }
550 
551 protected:
553  const float_type (*func)(float_type);
554 };
555 
566 template <typename float_type> class c2_const_ptr {
567 public:
569  c2_const_ptr() : func(0) {}
573  { this->set_function(&f); }
577  { this->set_function(src.get_ptr()); }
581  {
582  if(func) func->release_ownership();
583  func=f;
584  if(func) func->claim_ownership();
585  }
586 
590  { this->set_function(f.get_ptr()); return f; }
594  { this->set_function(&f); return f; }
602  {
603  if(func) func->release_ownership_for_return();
604  func=0;
605  }
609  void unset_function(void) { this->set_function(0); }
611  ~c2_const_ptr() { this->set_function(0); }
612 
614  inline const c2_function<float_type> &get() const throw(c2_exception)
615  {
616  if(!func) throw c2_exception("c2_ptr accessed uninitialized");
617  return *func;
618  }
620  inline const c2_function<float_type> *get_ptr() const { return func; }
622  inline const c2_function<float_type> *operator -> () const
623  { return &get(); }
625  bool valid() const { return func != 0; }
626 
628  operator const c2_function<float_type>& () const { return this->get(); }
629 
636  float_type operator()(float_type x) const throw(c2_exception) { return get()(x); }
645  float_type operator()(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
646  { return get().value_with_derivatives(x, yprime, yprime2); }
651  { return *new c2_sum_p<float_type>(get(), rhs); }
656  { return *new c2_diff_p<float_type>(get(), rhs); }
661  { return *new c2_product_p<float_type>(get(), rhs); }
666  { return *new c2_ratio_p<float_type>(get(), rhs); }
671  { return *new c2_composed_function_p<float_type>(get(), inner); }
672 
673 protected:
675 };
676 
681 template <typename float_type> class c2_ptr : public c2_const_ptr<float_type >
682 {
683 public:
685  c2_ptr() : c2_const_ptr<float_type>() {}
689  c2_const_ptr<float_type>() { this->set_function(&f); }
692  c2_ptr(const c2_ptr<float_type> &src) :
693  c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
695  inline c2_function<float_type> &get() const throw(c2_exception)
696  { return *const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get()); }
699  { return const_cast<c2_function<float_type>*>(this->func); }
702  { return &get(); }
706  { this->set_function(f.get_ptr()); return f; }
710  { this->set_function(&f); return f; }
711 private:
716 };
717 
725 template <typename float_type, template <typename> class c2_class > class c2_typed_ptr : public c2_const_ptr<float_type> {
726 public:
728  c2_typed_ptr() : c2_ptr<float_type>() {}
731  c2_typed_ptr(c2_class<float_type> &f)
732  : c2_const_ptr<float_type>() { this->set_function(&f); }
736  : c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
737 
739  inline c2_class<float_type> &get() const throw(c2_exception)
740  {
741  return *static_cast<c2_class<float_type> *>(const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get()));
742  }
744  inline c2_class<float_type> *operator -> () const
745  { return &get(); }
747  inline c2_class<float_type> *get_ptr() const
748  { return static_cast<c2_class<float_type> *>(const_cast<c2_function<float_type>*>(this->func)); }
750  operator c2_class<float_type>&() const { return get(); }
754  { this->set_function(f.get_ptr()); }
757  void operator =(c2_class<float_type> &f)
758  { this->set_function(&f); }
759 private:
764 };
765 
782 template <typename float_type=double> class c2_plugin_function_p :
783  public c2_function<float_type> {
784 public:
786  c2_plugin_function_p() : c2_function<float_type>(), func() {}
789  c2_function<float_type>(),func(f) { }
793  {
794  func.set_function(f);
795  if(f) this->set_domain(f->xmin(), f->xmax());
796  }
799  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
800  {
801  if(!func.valid()) throw c2_exception("c2_plugin_function_p called uninitialized");
802  return func->value_with_derivatives(x, yprime, yprime2);
803  }
805  virtual ~c2_plugin_function_p() { }
806 
808  void unset_function() { func.unset_function(); }
809 
810  virtual void get_sampling_grid(float_type amin, float_type amax, std::vector<float_type> &grid) const
811  {
812  if(!func.valid()) throw c2_exception("c2_plugin_function_p called uninitialized");
813  if(this->sampling_grid) c2_function<float_type>::get_sampling_grid(amin, amax, grid);
814  else func->get_sampling_grid(amin, amax, grid);
815  }
816 protected:
818 };
819 
824 template <typename float_type=double> class c2_const_plugin_function_p : public c2_plugin_function_p<float_type> {
825 public:
830  c2_plugin_function_p<float_type>() { this->set_function(&f); }
836 
838  const c2_function<float_type> &get() const throw(c2_exception)
839  { return this->func.get(); }
840 };
841 
845 template <typename float_type=double> class c2_binary_function : public c2_function<float_type> {
846 public:
849 
850  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception)
851  {
852  if(stub) throw c2_exception("attempt to evaluate a c2_binary_function stub");
853  return this->combine(*Left.get_ptr(), *Right.get_ptr(), x, yprime, yprime2);
854  }
855 
858  virtual ~c2_binary_function() { }
859 
860 protected:
866  float_type (*combiner)(const c2_function<float_type> &left, const c2_function<float_type> &right,
867  float_type x, float_type *yprime, float_type *yprime2),
868  const c2_function<float_type> &left, const c2_function<float_type> &right) :
869  c2_function<float_type>(), combine(combiner), Left(left), Right(right), stub(false)
870  {
871  this->set_domain(
872  (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(),
873  (left.xmax() < right.xmax()) ? left.xmax() : right.xmax()
874  );
875  }
876 
880  float_type (*combiner)(const c2_function<float_type> &left, const c2_function<float_type> &right,
881  float_type x, float_type *yprime, float_type *yprime2)
882  ) : c2_function<float_type>(), combine(combiner), Left(), Right(), stub(true) { }
883 
884 public:
885  float_type (* const combine)(const c2_function<float_type> &left, const c2_function<float_type> &right,
886  float_type x, float_type *yprime, float_type *yprime2);
887 
888 protected:
891  bool stub;
892 
893 };
894 
899 template <typename float_type=double> class c2_scaled_function_p : public c2_function<float_type> {
900 public:
905  c2_scaled_function_p(const c2_function<float_type> &outer, float_type scale) :
906  c2_function<float_type>(), func(outer), yscale(scale) { }
907 
910  void reset(float_type scale) { yscale=scale; }
911 
915  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception)
916  {
917  float_type y=this->func->value_with_derivatives(x, yprime, yprime2);
918  if(yprime) (*yprime)*=yscale;
919  if(yprime2) (*yprime2)*=yscale;
920  return y*yscale;
921  }
922 
923 protected:
924  c2_scaled_function_p<float_type>() : func() {} // hide default constructor, since its use is almost always an error.
927  float_type yscale;
928 };
929 
937 template <typename float_type=double> class c2_cached_function_p : public c2_function<float_type> {
938 public:
943  func(f), init(false) {}
948  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
949  {
950  if(!init || x != x0) {
951  y=this->func->value_with_derivatives(x, &yp, &ypp);
952  x0=x;
953  init=true;
954  }
955  if(yprime) *yprime=yp;
956  if(yprime2) *yprime2=ypp;
957  return y;
958  }
959 
960 protected:
961  c2_cached_function_p() : func() {} // hide default constructor, since its use is almost always an error.
963  mutable bool init;
964  mutable float_type x0, y, yp, ypp;
965 
966 };
967 
973 template <typename float_type=double> class c2_composed_function_p : public c2_binary_function<float_type> {
974 public:
975 
981  c2_binary_function<float_type>(combine, outer, inner) { this->set_domain(inner.xmin(), inner.xmax()); }
984 
987  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
988  {
989  float_type y0, y1;
990  if(yprime || yprime2) {
991  float_type yp0, ypp0, yp1, ypp1;
992  y0=right.value_with_derivatives(x, &yp0, &ypp0);
993  y1=left.value_with_derivatives(y0, &yp1, &ypp1);
994  if(yprime) *yprime=yp1*yp0;
995  if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1;
996  } else {
997  y0=right(x);
998  y1=left(y0);
999  }
1000  return y1;
1001  }
1002 };
1003 
1007 template <typename float_type=double> class c2_sum_p : public c2_binary_function<float_type> {
1008 public:
1014  c2_sum_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
1015 
1018  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
1019  {
1020  float_type y0, y1;
1021  if(yprime || yprime2) {
1022  float_type yp0, ypp0, yp1, ypp1;
1023  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1024  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1025  if(yprime) *yprime=yp0+yp1;
1026  if(yprime2) *yprime2=ypp0+ypp1;
1027  } else {
1028  y0=left(x);
1029  y1=right(x);
1030  }
1031  return y0+y1;
1032  }
1033 };
1034 
1035 
1039 template <typename float_type=double> class c2_diff_p : public c2_binary_function<float_type> {
1040 public:
1046  c2_diff_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
1047 
1050  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
1051  {
1052  float_type y0, y1;
1053  if(yprime || yprime2) {
1054  float_type yp0, ypp0, yp1, ypp1;
1055  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1056  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1057  if(yprime) *yprime=yp0-yp1;
1058  if(yprime2) *yprime2=ypp0-ypp1;
1059  } else {
1060  y0=left(x);
1061  y1=right(x);
1062  }
1063  return y0-y1;
1064  }
1065 };
1066 
1067 
1071 template <typename float_type=double> class c2_product_p : public c2_binary_function<float_type> {
1072 public:
1078  c2_product_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
1079 
1082  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
1083  {
1084  float_type y0, y1;
1085  if(yprime || yprime2) {
1086  float_type yp0, ypp0, yp1, ypp1;
1087  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1088  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1089  if(yprime) *yprime=y1*yp0+y0*yp1;
1090  if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0;
1091  } else {
1092  y0=left(x);
1093  y1=right(x);
1094  }
1095  return y0*y1;
1096  }
1097 };
1098 
1099 
1103 template <typename float_type=double> class c2_ratio_p : public c2_binary_function<float_type> {
1104 public:
1110  c2_ratio_p() : c2_binary_function<float_type>(combine) {} ; // create a stub just for the combiner to avoid statics
1111 
1114  float_type x, float_type *yprime, float_type *yprime2) throw(c2_exception)
1115  {
1116  float_type y0, y1;
1117  if(yprime || yprime2) {
1118  float_type yp0, ypp0, yp1, ypp1;
1119  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1120  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1121  if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio
1122  if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)/(y1*y1*y1);
1123  } else {
1124  y0=left(x);
1125  y1=right(x);
1126  }
1127  return y0/y1;
1128  }
1129 
1130 };
1131 
1136 template <typename float_type> class c2_constant_p : public c2_function<float_type> {
1137 public:
1138  c2_constant_p(float_type x) : c2_function<float_type>(), value(x) {}
1139  void reset(float_type val) { value=val; }
1140  virtual float_type value_with_derivatives(float_type, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1141  { if(yprime) *yprime=0; if(yprime2) *yprime2=0; return value; }
1142 
1143 private:
1144  float_type value;
1145 };
1146 
1149 template <typename float_type> class c2_transformation {
1150 public:
1157  c2_transformation(bool transformed,
1158  float_type (*xin)(float_type), float_type (*xinp)(float_type), float_type (*xinpp)(float_type), float_type (*xout)(float_type)
1159  ) :
1160  fTransformed(transformed), fHasStaticTransforms(true),
1161  pIn(xin), pInPrime(xinp), pInDPrime(xinpp), pOut(xout) { }
1162 
1165  c2_transformation(bool transformed) :
1166  fTransformed(transformed), fHasStaticTransforms(false),
1169  virtual ~c2_transformation() { }
1171  const bool fTransformed;
1174 
1178  float_type (* const pIn)(float_type);
1180  float_type (* const pInPrime)(float_type);
1182  float_type (* const pInDPrime)(float_type);
1184  float_type (* const pOut)(float_type);
1185 
1187  virtual float_type fIn(float_type x) const { return pIn(x); }
1189  virtual float_type fInPrime(float_type x) const { return pInPrime(x); }
1191  virtual float_type fInDPrime(float_type x) const { return pInDPrime(x); }
1193  virtual float_type fOut(float_type x) const { return pOut(x); }
1194 
1195 protected:
1197  static float_type report_error(float_type x) { throw c2_exception("use of improperly constructed axis transform"); return x; }
1199  static float_type ident(float_type x) { return x; }
1201  static float_type one(float_type) { return 1; }
1203  static float_type zero(float_type) { return 0; }
1205  static float_type recip(float_type x) { return 1.0/x; }
1207  static float_type recip_prime(float_type x) { return -1/(x*x); }
1209  static float_type recip_prime2(float_type x) { return 2/(x*x*x); }
1210 
1211 };
1212 
1215 template <typename float_type> class c2_transformation_linear : public c2_transformation<float_type> {
1216 public:
1218  c2_transformation_linear() : c2_transformation<float_type>(false, this->ident, this->one, this->zero, this->ident) { }
1221 };
1224 template <typename float_type> class c2_transformation_log : public c2_transformation<float_type> {
1225 public:
1227  c2_transformation_log() : c2_transformation<float_type>(true, std::log, this->recip, this->recip_prime, std::exp) { }
1230 };
1233 template <typename float_type> class c2_transformation_recip : public c2_transformation<float_type> {
1234 public:
1236  c2_transformation_recip() : c2_transformation<float_type>(true, this->recip, this->recip_prime, this->recip_prime2, this->recip) { }
1239 };
1240 
1246 template <typename float_type>
1248 public:
1254  isIdentity(!(xx.fTransformed || yy.fTransformed)), X(xx), Y(yy) { }
1256  virtual ~c2_function_transformation() { delete &X; delete &Y; }
1265  virtual float_type evaluate(float_type xraw,
1266  float_type y, float_type yp0, float_type ypp0,
1267  float_type *yprime, float_type *yprime2) const;
1269  const bool isIdentity;
1274 };
1275 
1279 template <typename float_type> class c2_lin_lin_function_transformation :
1280  public c2_function_transformation<float_type> {
1281 public:
1283  c2_function_transformation<float_type>(
1284  *new c2_transformation_linear<float_type>,
1285  *new c2_transformation_linear<float_type>
1286  ) { }
1288 };
1289 
1293 template <typename float_type> class c2_log_log_function_transformation :
1294  public c2_function_transformation<float_type> {
1295 public:
1297  c2_function_transformation<float_type>(
1298  *new c2_transformation_log<float_type>,
1299  *new c2_transformation_log<float_type>
1300  ) { }
1302 };
1303 
1307 template <typename float_type> class c2_lin_log_function_transformation :
1308  public c2_function_transformation<float_type> {
1309 public:
1311  c2_function_transformation<float_type>(
1312  *new c2_transformation_linear<float_type>,
1313  *new c2_transformation_log<float_type>
1314  ) { }
1316 };
1317 
1321 template <typename float_type> class c2_log_lin_function_transformation :
1322  public c2_function_transformation<float_type> {
1323 public:
1325  c2_function_transformation<float_type>(
1326  *new c2_transformation_log<float_type>,
1327  *new c2_transformation_linear<float_type>
1328  ) { }
1330 };
1331 
1335 template <typename float_type> class c2_arrhenius_function_transformation :
1336  public c2_function_transformation<float_type> {
1337 public:
1339  c2_function_transformation<float_type>(
1340  *new c2_transformation_recip<float_type>,
1341  *new c2_transformation_log<float_type>
1342  ) { }
1344 };
1345 
1377 template <typename float_type=double> class interpolating_function_p : public c2_function<float_type> {
1378 public:
1384  fTransform(*new c2_lin_lin_function_transformation<float_type>) { }
1385 
1389  fTransform(transform) { }
1390 
1401  interpolating_function_p<float_type> & load(const std::vector<float_type> &x, const std::vector<float_type> &f,
1402  bool lowerSlopeNatural, float_type lowerSlope,
1403  bool upperSlopeNatural, float_type upperSlope, bool splined=true
1404  ) throw(c2_exception);
1405 
1414  interpolating_function_p<float_type> & load_pairs(
1415  std::vector<std::pair<float_type, float_type> > &data,
1416  bool lowerSlopeNatural, float_type lowerSlope,
1417  bool upperSlopeNatural, float_type upperSlope, bool splined=true
1418  ) throw(c2_exception);
1419 
1437  interpolating_function_p<float_type> & sample_function(const c2_function<float_type> &func,
1438  float_type amin, float_type amax, float_type abs_tol, float_type rel_tol,
1439  bool lowerSlopeNatural, float_type lowerSlope,
1440  bool upperSlopeNatural, float_type upperSlope
1441  ) throw(c2_exception);
1442 
1443 
1458  const std::vector<float_type> &bincenters, const c2_function<float_type> &binheights)
1459  throw(c2_exception);
1460 
1479  const std::vector<float_type> &bins, const std::vector<float_type> &binheights, bool splined=true)
1480  throw(c2_exception);
1481 
1482  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception);
1483 
1485  virtual ~interpolating_function_p() { delete &fTransform; }
1486 
1489  { return *new interpolating_function_p<float_type>(); }
1490 
1497  void get_data(std::vector<float_type> &xvals, std::vector<float_type> &yvals) const throw() ;
1498 
1506  std::vector<float_type> &xvals,
1507  std::vector<float_type> &yvals,
1508  std::vector<float_type> &y2vals) const
1509  { xvals=X; yvals=F; y2vals=y2; }
1510 
1516  void set_lower_extrapolation(float_type bound);
1522  void set_upper_extrapolation(float_type bound);
1523 
1524  // these functions correctly combine the interpolating function with another interpolating function
1525  // preserving the X bounds and mapping functions of the host (left hand) function.
1526 
1534 
1543  const c2_binary_function<float_type> *combining_stub
1544  ) const;
1549  return binary_operator(rhs, new c2_sum_p<float_type>()); }
1554  return binary_operator(rhs, new c2_diff_p<float_type>()); }
1559  return binary_operator(rhs, new c2_product_p<float_type>()); }
1564  return binary_operator(rhs, new c2_ratio_p<float_type>()); }
1569  Xraw=rhs.Xraw; X=rhs.X; F=rhs.F; y2=rhs.y2;
1571  }
1572 
1574 
1575 protected:
1577  void spline(
1578  bool lowerSlopeNatural, float_type lowerSlope,
1579  bool upperSlopeNatural, float_type upperSlope
1580  ) throw(c2_exception);
1581 
1582  // This is for sorting the data. It must be static if it's going to be a class member.
1583  static bool comp_pair(std::pair<float_type,float_type> const &i, std::pair<float_type,float_type> const &j) {return i.first<j.first;}
1584 
1585  std::vector<float_type> Xraw, X, F, y2;
1588  mutable size_t lastKLow;
1589 };
1590 
1597 template <typename float_type=double> class log_lin_interpolating_function_p : public interpolating_function_p <float_type> {
1598 public:
1602  { }
1606 };
1607 
1608 
1614 template <typename float_type=double> class lin_log_interpolating_function_p : public interpolating_function_p <float_type> {
1615 public:
1619  { }
1623 };
1624 
1625 
1631 template <typename float_type=double> class log_log_interpolating_function_p : public interpolating_function_p <float_type> {
1632 public:
1636  { }
1640 };
1641 
1642 
1649 template <typename float_type=double> class arrhenius_interpolating_function_p : public interpolating_function_p <float_type> {
1650 public:
1654  { }
1658 };
1659 
1664 template <typename float_type=double> class c2_sin_p : public c2_function<float_type> {
1665 public:
1667  c2_sin_p() : c2_function<float_type>() {}
1668 
1669  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1670  { float_type q=std::sin(x); if(yprime) *yprime=std::cos(x); if(yprime2) *yprime2=-q; return q; }
1671 
1676  virtual void get_sampling_grid(float_type amin, float_type amax, std::vector<float_type> &grid) const;
1677 };
1678 
1683 template <typename float_type=double> class c2_cos_p : public c2_sin_p<float_type> {
1684 public:
1686  c2_cos_p() : c2_sin_p<float_type>() {}
1687 
1688  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1689  { float_type q=std::cos(x); if(yprime) *yprime=-std::sin(x); if(yprime2) *yprime2=-q; return q; }
1690 };
1691 
1696 template <typename float_type=double> class c2_tan_p : public c2_function<float_type> {
1697 public:
1699  c2_tan_p() : c2_function<float_type>() {}
1700 
1701  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1702  {
1703  float_type c=std::cos(x), ss=std::sin(x);
1704  float_type t=ss/c;
1705  float_type yp=1/(c*c);
1706  if(yprime) *yprime=yp; if(yprime2) *yprime2=2*t*yp;
1707  return t;
1708  }
1709 };
1710 
1715 template <typename float_type=double> class c2_log_p : public c2_function<float_type> {
1716 public:
1718  c2_log_p() : c2_function<float_type>() {}
1719 
1720  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1721  { if(yprime) *yprime=1.0/x; if(yprime2) *yprime2=-1.0/(x*x); return std::log(x); }
1722 };
1723 
1728 template <typename float_type=double> class c2_exp_p : public c2_function<float_type> {
1729 public:
1731  c2_exp_p() : c2_function<float_type>() {}
1732 
1733  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1734  { float_type q=std::exp(x); if(yprime) *yprime=q; if(yprime2) *yprime2=q; return q; }
1735 };
1736 
1741 template <typename float_type=double> class c2_sqrt_p : public c2_function<float_type> {
1742 public:
1744  c2_sqrt_p() : c2_function<float_type>() {}
1745 
1746  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1747  { float_type q=std::sqrt(x); if(yprime) *yprime=0.5/q; if(yprime2) *yprime2=-0.25/(x*q); return q; }
1748 };
1749 
1754 template <typename float_type=double> class c2_recip_p : public c2_function<float_type> {
1755 public:
1757  c2_recip_p(float_type scale) : c2_function<float_type>(), rscale(scale) {}
1758 
1759  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1760  {
1761  float_type q=1.0/x;
1762  float_type y=rscale*q;
1763  if(yprime) *yprime=-y*q;
1764  if(yprime2) *yprime2=2*y*q*q;
1765  return y;
1766  }
1769  void reset(float_type scale) { rscale=scale; }
1770 private:
1771  float_type rscale;
1772 };
1773 
1778 template <typename float_type=double> class c2_identity_p : public c2_function<float_type> {
1779 public:
1781  c2_identity_p() : c2_function<float_type>() {}
1782 
1783  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1784  { if(yprime) *yprime=1.0; if(yprime2) *yprime2=0; return x; }
1785 };
1786 
1798 template <typename float_type=double> class c2_linear_p : public c2_function<float_type> {
1799 public:
1804  c2_linear_p(float_type x0, float_type y0, float_type slope) :
1805  c2_function<float_type>(), xint(x0), intercept(y0), m(slope) {}
1810  void reset(float_type x0, float_type y0, float_type slope) { xint=x0; intercept=y0; m=slope; }
1811  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1812  { if(yprime) *yprime=m; if(yprime2) *yprime2=0; return m*(x-xint)+intercept; }
1813 
1814 private:
1815  float_type xint, intercept, m;
1816 protected:
1817  c2_linear_p() {} // do not allow naked construction... it is usually an accident.
1818 };
1819 
1833 template <typename float_type=double> class c2_quadratic_p : public c2_function<float_type> {
1834 public:
1840  c2_quadratic_p(float_type x0, float_type y0, float_type xcoef, float_type x2coef) :
1841  c2_function<float_type>(), intercept(y0), center(x0), a(x2coef), b(xcoef) {}
1847  void reset(float_type x0, float_type y0, float_type xcoef, float_type x2coef) { intercept=y0; center=x0; a=x2coef; b=xcoef; }
1848  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1849  { float_type dx=x-center; if(yprime) *yprime=2*a*dx+b; if(yprime2) *yprime2=2*a; return a*dx*dx+b*dx+intercept; }
1850 
1851 private:
1852  float_type intercept, center, a, b;
1853 protected:
1854  c2_quadratic_p() {} // do not allow naked construction... it is usually an accident.
1855 };
1856 
1869 template <typename float_type=double> class c2_power_law_p : public c2_function<float_type> {
1870 public:
1874  c2_power_law_p(float_type scale, float_type power) :
1875  c2_function<float_type>(), a(scale), b(power) {}
1879  void reset(float_type scale, float_type power) { a=scale; b=power; }
1880  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception)
1881  { float_type q=a*std::pow(x,b-2); if(yprime) *yprime=b*q*x; if(yprime2) *yprime2=b*(b-1)*q; return q*x*x; }
1882 
1883 private:
1884  float_type a, b;
1885 protected:
1886  c2_power_law_p() {} // do not allow naked construction... it is usually an accident.
1887 };
1888 
1907 template <typename float_type=double> class c2_inverse_function_p : public c2_function<float_type> {
1908 public:
1912  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw(c2_exception);
1913 
1916  void set_start_hint(float_type hint) const { start_hint=hint; }
1917 
1923  virtual float_type get_start_hint(float_type x) const
1924  { return hinting_function.valid()? hinting_function(x) : start_hint; }
1925 
1943  { hinting_function.set_function(hint_func); }
1949  { hinting_function=hint_func; }
1950 
1951 protected:
1952  c2_inverse_function_p() {} // do not allow naked construction... it is usually an accident.
1953  mutable float_type start_hint;
1956 };
1957 
1968 template <typename float_type=double> class accumulated_histogram : public interpolating_function_p <float_type> {
1969 public:
1976  accumulated_histogram(const std::vector<float_type>binedges, const std::vector<float_type> binheights,
1977  bool normalize=false, bool inverse_function=false, bool drop_zeros=true);
1978 
1979 };
1980 
1996 template <typename float_type=double> class c2_connector_function_p : public c2_function<float_type> {
1997 public:
2006  c2_connector_function_p(float_type x0, const c2_function<float_type> &f0, float_type x2, const c2_function<float_type> &f2,
2007  bool auto_center, float_type y1);
2022  float_type x0, float_type y0, float_type yp0, float_type ypp0,
2023  float_type x2, float_type y2, float_type yp2, float_type ypp2,
2024  bool auto_center, float_type y1);
2032  const c2_fblock<float_type> &fb0,
2033  const c2_fblock<float_type> &fb2,
2034  bool auto_center, float_type y1);
2035 
2037  virtual ~c2_connector_function_p();
2038  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception);
2039 protected:
2041  void init(
2042  const c2_fblock<float_type> &fb0,
2043  const c2_fblock<float_type> &fb2,
2044  bool auto_center, float_type y1);
2045 
2046  float_type fhinv, fy1, fa, fb, fc, fd, fe, ff;
2047 };
2048 
2068 template <typename float_type=double> class c2_piecewise_function_p : public c2_function<float_type> {
2069 public:
2073  virtual ~c2_piecewise_function_p();
2074  virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const throw (c2_exception);
2082  void append_function(const c2_function<float_type> &func) throw (c2_exception);
2083 protected:
2084  std::vector<c2_const_ptr<float_type> > functions;
2085  mutable int lastKLow;
2086 };
2087 
2088 #include "c2_function.icc"
2089 
2090 #endif
const c2_function< float_type > * func
Definition: c2_function.hh:674
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)
initialize from a grid of points and an std::vector of probability densities (un-normalized) to an in...
create a non-generic container for a c2_function which handles the reference counting.
Definition: c2_function.hh:725
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:176
const std::string cvs_file_vers() const
get versioning information for the source file
void reset(float_type scale)
set a new scale factor
Definition: c2_function.hh:910
c2_typed_ptr(c2_class< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:731
bool ypbad
flag, filled in by c2_function::fill_fblock(), indicating the derivative is NaN of Inf ...
Definition: c2_function.hh:113
bool stub
if true, we don't own any functions, we are just a source of a combining function.
Definition: c2_function.hh:891
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:788
virtual float_type fOut(float_type x) const
virtual output X transform
c2_linear_p(float_type x0, float_type y0, float_type slope)
Construct the operator f=y0 + slope * (x-x0)
void reset(float_type x0, float_type y0, float_type xcoef, float_type x2coef)
Modify the coefficients after construction.
void set_start_hint(float_type hint) const
give the function a hint as to where to look for its inverse
float_type xint
interpolating_function_p< float_type > & multiply_pointwise(const c2_function< float_type > &rhs) const
produce a newly resampled interpolating_function_p which is the specified product.
c2_piecewise_function_p()
construct the container
c2_function< float_type > * operator->() const
get a checked pointer to our owned function
Definition: c2_function.hh:701
c2_log_p()
constructor.
virtual void get_sampling_grid(float_type amin, float_type amax, std::vector< float_type > &grid) const
return a grid dynamically, suitable for use with trig functions with period 2*pi
void get_internal_data(std::vector< float_type > &xvals, std::vector< float_type > &yvals, std::vector< float_type > &y2vals) const
retrieve copies of the transformed x, y and y2 tables from which this was built
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
const bool fTransformed
flag to indicate if this transform is not the identity
float_type x
the abscissa
Definition: c2_function.hh:105
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
function to manage the binary operation, used by c2_binary_function::value_with_derivatives() ...
Definition: c2_function.hh:850
create a c2_function which is the ratio of two other c2_functions.This should always be constructed u...
Definition: c2_function.hh:83
Create a very lightweight method to return a scalar multiple of another function. \ \The factory func...
Definition: c2_function.hh:899
a container into which any conventional c-style function can be dropped, to create a degenerate c2_fu...
Definition: c2_function.hh:534
const c2_const_ptr< float_type > Left
Definition: c2_function.hh:889
the identity transform
c2_ptr(const c2_ptr< float_type > &src)
copy constructor
Definition: c2_function.hh:692
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
produce a newly resampled interpolating_function_p which is the specified sum.
float_type fXMax
Definition: c2_function.hh:449
static const G4double f2
float_type ypp
the second derivative at x
Definition: c2_function.hh:111
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_const_plugin_function_p(const c2_function< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:829
c2_const_ptr< float_type > sampler_function
log axis transform
c2_quadratic_p(float_type x0, float_type y0, float_type xcoef, float_type x2coef)
Construct the operator.
void claim_ownership() const
increment our reference count. Destruction is only legal if the count is zero.
Definition: c2_function.hh:406
std::vector< float_type > * xvals
Definition: c2_function.hh:493
c2_diff_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left - right
std::vector< float_type > y2
create a c2_function which is the difference of two other c2_functions.This should always be construc...
Definition: c2_function.hh:81
size_t get_evaluations() const
this is a counter owned by the function but which can be used to monitor efficiency of algorithms...
Definition: c2_function.hh:307
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
factory function to create a c2_ratio_p from a regular algebraic expression.
Definition: c2_function.hh:391
const c2_const_ptr< float_type > Right
Definition: c2_function.hh:889
create a container for a c2_function which handles the reference counting.
Definition: c2_function.hh:86
float_type rscale
A spline with X in reciprocal space and Y transformed in log space.Most useful for thermodynamic type...
c2_fblock< float_type > * f1
Definition: c2_function.hh:489
void preen_sampling_grid(std::vector< float_type > *result) const
clean up endpoints on a grid of points
interpolating_function_p< float_type > & divide_pointwise(const c2_function< float_type > &rhs) const
produce a newly resampled interpolating_function_p which is the specified ratio.
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:885
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
factory function to create a c2_diff_p from a regular algebraic expression.
Definition: c2_function.hh:655
virtual const char * what() const
Returns a C-style character string describing the general cause of the current error.
Definition: c2_function.hh:73
c2_function< float_type > & square_normalized_function(float_type amin, float_type amax, float_type norm=1.0) const
create a new c2_function from this one which is square-normalized on the interval ...
c2_cos_p()
constructor.
interpolating_function_p< float_type > & subtract_pointwise(const c2_function< float_type > &rhs) const
produce a newly resampled interpolating_function_p which is the specified difference.
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:827
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
c2_fblock< float_type > upper
Definition: c2_function.hh:501
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
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
c2_ratio_p()
Create a stub just for the combiner to avoid statics.
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:403
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:381
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
return the grid of 'interesting' points along this function which lie in the region requested ...
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:799
c2_fblock< float_type > f1
Definition: c2_function.hh:467
c2_cached_function_p(const c2_function< float_type > &f)
construct the container
Definition: c2_function.hh:942
a transformation of a function in and out of lin-lin space
virtual ~c2_piecewise_function_p()
destructor
#define c2_isnan
Definition: c2_function.hh:51
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
Definition: c2_function.hh:986
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:538
size_t evaluations
Definition: c2_function.hh:450
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 ...
c2_piecewise_function_p< float_type > * out
Definition: c2_function.hh:492
virtual void get_sampling_grid(float_type amin, float_type amax, std::vector< float_type > &grid) const
return the grid of 'interesting' points along this function which lie in the region requested ...
Definition: c2_function.hh:810
const c2_function< float_type > * operator->() const
get a checked pointer to our owned function
Definition: c2_function.hh:622
float_type integrate_step(struct c2_integrate_recur &rb) const
Carry out the recursive subdivision and integration.
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
void set_function(const c2_function< float_type > *f)
fill the container with a new function, or clear it with a null pointer
Definition: c2_function.hh:832
reciprocal axis transform
c2_class< float_type > * operator->() const
get a checked pointer to our owned function
Definition: c2_function.hh:744
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
for points in xgrid, adaptively return Integral[f(x),{x,xgrid[i],xgrid[i+1]}] and return in vector...
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
evaluate the transformation from internal coordinates to external coordinates
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.
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
float_type value
c2_product_p< float_type > & operator*(const c2_function< float_type > &rhs) const
factory function to create a c2_product_p from a regular algebraic expression.
Definition: c2_function.hh:660
std::vector< float_type > F
c2_plugin_function_p()
construct the container with no function
Definition: c2_function.hh:786
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:835
virtual ~c2_transformation()
the destructor
c2_function(const c2_function< float_type > &src)
Definition: c2_function.hh:428
float_type b
interpolating_function_p< float_type > & binary_operator(const c2_function< float_type > &rhs, const c2_binary_function< float_type > *combining_stub) const
create a new interpolating_function_p which is the parent interpolating_function_p combined with rhs ...
std::vector< float_type > * get_sampling_grid_pointer() const
get the sampling grid, which may be a null pointer
Definition: c2_function.hh:332
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
create a new, empty interpolating function of this type (virtual constructor)
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 'interesting' points on the function.
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:409
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
float_type a
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
c2_sum_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left + right
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...
float_type b
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:625
A spline with X transformed into log space.Most useful for functions looking like y=log(x) or any oth...
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.
Definition: c2_function.hh:915
A container into which any other c2_function can be dropped.It allows a function to be pre-evaluated ...
Definition: c2_function.hh:937
c2_ptr()
construct the container with no function
Definition: c2_function.hh:685
structure used to hold evaluated function data at a point.
Definition: c2_function.hh:101
create a quadratic mapping of another functionfor example, given a c2_function f
Definition: c2_function.hh:85
void reset(float_type x0, float_type y0, float_type slope)
Change the slope and intercepts after construction.
a c2_plugin_function_p which promises not to fiddle with the plugged function.The factory function c2...
Definition: c2_function.hh:824
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.
void append_function(const c2_function< float_type > &func)
append a new function to the sequence
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
void set_function(c2_function< float_type > *f)
fill the container with a new function, or clear it with a null pointer and copy our domain ...
Definition: c2_function.hh:792
compute x with its derivatives.The factory function c2_factory::identity() creates *new c2_identity_p...
size_t owner_count
Definition: c2_function.hh:525
void set_upper_extrapolation(float_type bound)
enable extrapolation of the function above the tabulated data.
c2_const_ptr(const c2_const_ptr< float_type > &src)
copy constructor
Definition: c2_function.hh:576
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)
set the domain for this function.
Definition: c2_function.hh:301
create a c2_function which is the product of two other c2_functions.This should always be constructed...
Definition: c2_function.hh:82
virtual interpolating_function_p< float_type > & clone() const
create a new, empty interpolating function of this type (virtual constructor)
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)
construct the binary function
Definition: c2_function.hh:865
virtual ~c2_function()
destructor
Definition: c2_function.hh:152
virtual ~c2_transformation_log()
destructor
structure used to hold root bracketing information
Definition: c2_function.hh:500
a transformation of a function in and out of lin-log space
c2_tan_p()
constructor.
float_type center
void release_for_return()
release the function without destroying it, so it can be returned from a function ...
Definition: c2_function.hh:601
std::vector< float_type > X
const c2_const_ptr< float_type > func
Definition: c2_function.hh:962
const c2_const_ptr< float_type > & operator=(const c2_const_ptr< float_type > &f)
fill the container from another container
Definition: c2_function.hh:589
std::vector< recur_item > * rb_stack
Definition: c2_function.hh:481
a container into which any other c2_function can be dropped, to allow expressions with replacable com...
Definition: c2_function.hh:782
const c2_function< float_type > & get() const
get a reference to our owned function
Definition: c2_function.hh:614
c2_product_p()
Create a stub just for the combiner to avoid statics.
float_type yp
the derivative at x
Definition: c2_function.hh:109
structure used to pass information recursively in integrator.
Definition: c2_function.hh:478
create a container for a c2_function which handles the reference counting.It is useful as a smart con...
Definition: c2_function.hh:566
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:553
virtual float_type value_with_derivatives(float_type x, float_type *yprime, float_type *yprime2) const
get the value and derivatives.
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:142
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 ...
void sample_step(struct c2_sample_recur &rb) const
Carry out the recursive subdivision for sampling.
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:620
c2_sum_p()
Create a stub just for the combiner to avoid statics.
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 ...
bool yppbad
flag, filled in by c2_function::fill_fblock(), indicating the second derivative is NaN of Inf ...
Definition: c2_function.hh:115
std::vector< c2_const_ptr< float_type > > functions
c2_transformation_recip()
constructor
the parent class for all c2_functions.
Definition: c2_function.hh:138
arrhenius_interpolating_function_p()
an empty arrhenius cubic-spline interpolating_function_p
create a c2_function which is a piecewise assembly of other c2_functions.The functions must have incr...
Definition: c2_function.hh:84
c2_sin_p()
constructor.
virtual ~c2_transformation_linear()
destructor
c2_typed_ptr(const c2_typed_ptr< float_type, c2_class > &src)
copy constructor
Definition: c2_function.hh:735
T max(const T t1, const T t2)
brief Return the largest of the two arguments
c2_const_ptr(const c2_function< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:572
c2_exception(const char msgcode[])
construct the exception with an error message
Definition: c2_function.hh:69
float_type intercept
the data element for the internal recursion stack for the sampler and integrator
Definition: c2_function.hh:466
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:449
c2_diff_p()
Create a stub just for the combiner to avoid statics.
interpolating_function_p< float_type > & unary_operator(const c2_function< float_type > &source) const
create a new interpolating_function_p which is the source function applied to every point in the inte...
c2_ratio_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left / right
structure used to pass information recursively in sampler.
Definition: c2_function.hh:488
c2_class< float_type > * get_ptr() const
get an unchecked pointer to our owned function
Definition: c2_function.hh:747
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:542
A spline with Y transformed into log space.Most useful for functions looking like y=exp(x) ...
c2_product_p< float_type > & operator*(const c2_function< float_type > &rhs) const
factory function to create a c2_product_p from a regular algebraic expression.
Definition: c2_function.hh:386
c2_power_law_p(float_type scale, float_type power)
Construct the operator.
float_type intercept
const G4double x[NPOINTSGL]
float_type(*const pInPrime)(float_type)
non-virtual pointer to input X transform derivative
const bool isIdentity
flag indicating of the transform is the identity, and can be skipped for efficiency ...
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:645
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))
construct a 'stub' c2_binary_function, which provides access to the combine() function ...
Definition: c2_function.hh:879
c2_scaled_function_p(const c2_function< float_type > &outer, float_type scale)
construct the function with its scale factor.
Definition: c2_function.hh:905
virtual interpolating_function_p< float_type > & clone() const
create a new, empty interpolating function of this type (virtual constructor)
~c2_const_ptr()
destructor
Definition: c2_function.hh:611
c2_ptr< float_type > func
Definition: c2_function.hh:817
virtual ~c2_function_transformation()
destructor
const c2_const_ptr< float_type > func
the scaling factor for the function
Definition: c2_function.hh:926
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:650
c2_product_p(const c2_function< float_type > &left, const c2_function< float_type > &right)
construct left * right
c2_const_ptr()
construct the container with no function
Definition: c2_function.hh:569
An interpolating_function_p which is the cumulative integral of a histogram.
std::vector< recur_item > * rb_stack
Definition: c2_function.hh:494
virtual float_type get_start_hint(float_type x) const
get the starting hint.
size_t count_owners() const
get the reference count, mostly for debugging
Definition: c2_function.hh:425
virtual interpolating_function_p< float_type > & clone() const
create a new, empty interpolating function of this type (virtual constructor)
lin_log_interpolating_function_p()
an empty linear-log cubic-spline interpolating_function_p
c2_fblock< float_type > * f0
Definition: c2_function.hh:489
c2_fblock< float_type > * f0
Definition: c2_function.hh:479
float_type y
the value of the function at x
Definition: c2_function.hh:107
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:753
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.
void reset(float_type val)
virtual ~c2_plugin_function_p()
destructor
Definition: c2_function.hh:805
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 refine_sampling_grid(std::vector< float_type > &grid, size_t refinement) const
refine a grid by splitting each interval into more intervals
c2_function< float_type > * get_ptr() const
get an unchecked pointer to our owned function
Definition: c2_function.hh:698
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:948
static PROLOG_HANDLER error
Definition: xmlrole.cc:112
std::vector< float_type > * sampling_grid
Definition: c2_function.hh:446
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:858
float_type a
c2_ratio_p< float_type > & operator/(const c2_function< float_type > &rhs) const
factory function to create a c2_ratio_p from a regular algebraic expression.
Definition: c2_function.hh:665
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:311
Provides function composition (nesting)This allows evaluation of f(g(x)) where f and g are c2_functio...
Definition: c2_function.hh:79
float_type bad_x_point
this point may be used to record where a calculation ran into trouble
Definition: c2_function.hh:452
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:447
create a c2_function which is the sum of two other c2_function objects.This should always be construc...
Definition: c2_function.hh:80
static float_type report_error(float_type x)
utility function for unimplemented conversion
struct c2_root_info * root_info
this carry a memory of the last root bracketing, to avoid the necessity of evaluating the function on...
Definition: c2_function.hh:523
void release_ownership() const
decrement our reference count. If the count reaches zero, destroy ourself.
Definition: c2_function.hh:420
void clone_data(const interpolating_function_p< float_type > &rhs)
copy data from another interpolating function. This only makes sense if the source function has the s...
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... a shortcut for the integrator & sampler
Definition: c2_function.hh:456
std::vector< float_type > * yvals
Definition: c2_function.hh:493
float_type operator()(float_type x) const
convenience operator to make us look like a function
Definition: c2_function.hh:636
virtual interpolating_function_p< float_type > & clone() const
create a new, empty interpolating function of this type (virtual constructor)
float_type xmax() const
return the upper bound of the domain for this function as set by set_domain()
Definition: c2_function.hh:299
c2_composed_function_p(const c2_function< float_type > &outer, const c2_function< float_type > &inner)
construct outer( inner (x))
Definition: c2_function.hh:980
c2_composed_function_p()
Create a stub just for the combiner to avoid statics.
Definition: c2_function.hh:983
void reset_evaluations() const
reset the counter
Definition: c2_function.hh:309
void set_lower_extrapolation(float_type bound)
enable extrapolation of the function below the tabulated data.
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
std::string info
Definition: c2_function.hh:75
void unset_function()
clear our function
Definition: c2_function.hh:808
float_type m
c2_sqrt_p()
constructor.
const c2_ptr< float_type > & operator=(const c2_ptr< float_type > &f)
fill the container from another container
Definition: c2_function.hh:705
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
return the lower bound of the domain for this function as set by set_domain()
Definition: c2_function.hh:297
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:549
c2_typed_ptr()
construct the container with no function
Definition: c2_function.hh:728
c2_fblock< float_type > lower
Definition: c2_function.hh:501
c2_ptr(c2_function< float_type > &f)
construct the container with a pre-defined function
Definition: c2_function.hh:688
create a power law mapping of another functionfor example, given a c2_function f
Provides support for c2_function objects which are constructed from two other c2_function objects...
Definition: c2_function.hh:845
virtual void set_sampling_grid_pointer(std::vector< float_type > &grid)
Definition: c2_function.hh:440
void get_data(std::vector< float_type > &xvals, std::vector< float_type > &yvals) const
retrieve copies of the x & y tables from which this was built
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)
float_type previous_estimate
Definition: c2_function.hh:468
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:609
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
compute log(x) with its derivatives.The factory function c2_factory::log() creates *new c2_log_p ...
c2_fblock< float_type > * f1
Definition: c2_function.hh:479
void set_function(const c2_function< float_type > *f)
fill the container with a new function, or clear it with a null pointer
Definition: c2_function.hh:580