ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
c2_function.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file c2_function.hh
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 
40 #ifndef __has_c2_function_hh
41 #define __has_c2_function_hh 1
42 
43 // MSVC does not automatically define numerical constants such as M_PI
44 // this came from the msdn website, so it should be right...
45 #ifdef _MSC_VER
46 #define _USE_MATH_DEFINES
47 #define c2_isnan _isnan
48 #define c2_isfinite _finite
49 #else
50 #define c2_isnan std::isnan
51 #define c2_isfinite std::isfinite
52 #endif
53 
54 #include <cmath>
55 #include <vector>
56 #include <utility>
57 #include <string>
58 #include <stdexcept>
59 #include <typeinfo>
60 #include <sstream>
61 #include <limits> // fails under gcc-4.3 without this here, was ok
62 
64 class c2_exception : public std::exception {
65 public:
68  c2_exception(const char msgcode[]) : info(msgcode) { }
69  virtual ~c2_exception() throw() { }
72  virtual const char* what() const throw() { return info.c_str(); }
73 private:
74  std::string info;
75 };
76 
77 // put these forward references here, and with a bogus typename to make swig
78 template <typename float_type> class c2_composed_function_p;
79 template <typename float_type> class c2_sum_p;
80 template <typename float_type> class c2_diff_p;
81 template <typename float_type> class c2_product_p;
82 template <typename float_type> class c2_ratio_p;
83 template <typename float_type> class c2_piecewise_function_p;
84 template <typename float_type> class c2_quadratic_p;
85 template <typename float_type> class c2_ptr;
100 
101 
102 
103 template <typename float_type> class c2_fblock
104 {
105 public:
107  float_type x;
109  float_type y;
111  float_type yp;
113  float_type ypp;
116  bool ypbad;
118  bool yppbad;
119 };
120 
146 template <typename float_type=double> class c2_function {
147 public:
150  const std::string cvs_header_vers() const { return
151  "c2_function.hh 490 2012-04-10 19:05:40Z marcus ";
152  }
153 
156  const std::string cvs_file_vers() const ;
157 
158 public:
160  virtual ~c2_function() {
161  if(sampling_grid && !no_overwrite_grid) delete sampling_grid;
162  if(root_info) delete root_info;
163  if(owner_count) {
164  std::ostringstream outstr;
165  outstr << "attempt to delete an object with non-zero ownership in class";
166  outstr << typeid(*this).name() << std::endl;
167  //throw c2_exception(outstr.str().c_str());
168  }
169  }
170 
179  virtual float_type value_with_derivatives(float_type x, float_type *yprime,
180  float_type *yprime2) const /* throw(c2_exception) */ =0 ;
181 
185  inline float_type operator () (float_type x) const /* throw(c2_exception) */
186  { return value_with_derivatives(x, (float_type *)0, (float_type *)0); }
187 
194  inline float_type operator () (float_type x, float_type *yprime,
195  float_type *yprime2) const /* throw(c2_exception) */
196  { return value_with_derivatives(x, yprime, yprime2); }
197 
230  float_type find_root(float_type lower_bracket, float_type upper_bracket,
231  float_type start,
232  float_type value, int *error=0,
233  float_type *final_yprime=0,
234  float_type *final_yprime2=0 ) const /* throw(c2_exception) */;
266  float_type partial_integrals(std::vector<float_type> xgrid,
267  std::vector<float_type> *partials = 0,
268  float_type abs_tol=1e-12,
269  float_type rel_tol=1e-12,
270  int derivs=2, bool adapt=true,
271  bool extrapolate=true)
272  const /* throw(c2_exception) */;
273 
305  float_type integral(float_type amin, float_type amax,
306  std::vector<float_type> *partials = 0,
307  float_type abs_tol=1e-12, float_type rel_tol=1e-12,
308  int derivs=2, bool adapt=true, bool extrapolate=true)
309  const /* throw(c2_exception) */;
310 
382  c2_piecewise_function_p<float_type> *adaptively_sample(
383  float_type amin, float_type amax,
384  float_type abs_tol=1e-12, float_type rel_tol=1e-12,
385  int derivs=2, std::vector<float_type> *xvals=0,
386  std::vector<float_type> *yvals=0) const /* throw(c2_exception) */;
387 
388  inline float_type xmin() const { return fXMin; }
389  inline float_type xmax() const { return fXMax; }
390  void set_domain(float_type amin, float_type amax) { fXMin=amin; fXMax=amax; }
391 
394  size_t get_evaluations() const { return evaluations; }
396  void reset_evaluations() const { evaluations=0; }
398  inline void increment_evaluations() const { evaluations++; }
399 
407  bool check_monotonicity(const std::vector<float_type> &data,
408  const char message[]) const /* throw(c2_exception) */;
409 
426  virtual void set_sampling_grid(const std::vector<float_type> &grid)
427  /* throw(c2_exception) */;
428 
431  std::vector<float_type> *get_sampling_grid_pointer() const {
432  return sampling_grid; }
433 
434  virtual void get_sampling_grid(float_type amin, float_type amax,
435  std::vector<float_type> &grid) const ;
436 
438  void preen_sampling_grid(std::vector<float_type> *result) const;
439  void refine_sampling_grid(std::vector<float_type> &grid,
440  size_t refinement) const;
441 
448  c2_function<float_type> &normalized_function(
449  float_type amin,
450  float_type amax,
451  float_type norm=1.0) const /* throw(c2_exception) */;
452  c2_function<float_type> &square_normalized_function(
453  float_type amin, float_type amax,
454  float_type norm=1.0)
455  const /* throw(c2_exception) */;
463  c2_function<float_type> &square_normalized_function(
464  float_type amin, float_type amax, const c2_function<float_type> &weight,
465  float_type norm=1.0)
466  const /* throw(c2_exception) */;
467 
472  c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs) const
473  { return *new c2_sum_p<float_type>(*this, rhs); }
478  c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs) const
479  { return *new c2_diff_p<float_type>(*this, rhs); }
484  c2_product_p<float_type> &operator *
485  (const c2_function<float_type> &rhs) const
486  { return *new c2_product_p<float_type>(*this, rhs); }
487  c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs) const
488  { return *new c2_ratio_p<float_type>(*this, rhs); }
493  c2_composed_function_p<float_type> & operator ()
494  (const c2_function<float_type> &inner) const
495  { return *new c2_composed_function_p<float_type>((*this), inner); }
496 
501  float_type get_trouble_point() const { return bad_x_point; }
502 
505  void claim_ownership() const { owner_count++; }
508  size_t release_ownership_for_return() const /* throw(c2_exception) */ {
509  if(!owner_count) {
510  std::ostringstream outstr;
511  outstr << "attempt to release ownership of an unowned function in class ";
512  outstr << typeid(*this).name() << std::endl;
513  throw c2_exception(outstr.str().c_str());
514  }
515  owner_count--;
516  return owner_count;
517  }
518  void release_ownership() const /* throw(c2_exception) */ {
519  if(!release_ownership_for_return()) delete this;
520  }
523  size_t count_owners() const { return owner_count; }
524 
525 protected:
526  c2_function(const c2_function<float_type> &src)
527  : sampling_grid(0),
528  no_overwrite_grid(false),
529  fXMin(src.fXMin), fXMax(src.fXMax), root_info(0), owner_count(0)
530  {} // copy constructor only copies domain, and is only for internal use
531  c2_function() :
532  sampling_grid(0), no_overwrite_grid(0),
533  fXMin(-std::numeric_limits<float_type>::max()),
534  fXMax(std::numeric_limits<float_type>::max()), root_info(0), owner_count(0)
535  {}
536 
537  virtual void set_sampling_grid_pointer(std::vector<float_type> &grid)
538  {
539  if (sampling_grid && !no_overwrite_grid) delete sampling_grid;
540  sampling_grid=&grid; no_overwrite_grid=1;
541  }
542 
543  std::vector<float_type> * sampling_grid;
544  bool no_overwrite_grid;
545 
546  float_type fXMin, fXMax;
547  mutable size_t evaluations;
550  mutable float_type bad_x_point;
551 public:
555  inline void fill_fblock(c2_fblock<float_type> &fb) const /* throw(c2_exception) */
556  {
557  fb.y=value_with_derivatives(fb.x, &fb.yp, &fb.ypp);
558  fb.ypbad=c2_isnan(fb.yp) || !c2_isfinite(fb.yp);
559  fb.yppbad=c2_isnan(fb.ypp) || !c2_isfinite(fb.ypp);
560  increment_evaluations();
561  }
562 
563 private:
566  struct recur_item {
567  c2_fblock<float_type> f1; size_t depth;
568  float_type previous_estimate, abs_tol, step_sum;
569  bool done;
570  size_t f0index, f2index;
571  };
572 
573 
578  struct c2_integrate_recur {
579  c2_fblock<float_type> *f0, *f1;
580  float_type abs_tol, rel_tol, eps_scale, extrap_coef, extrap2,
581  dx_tolerance, abs_tol_min;
582  std::vector< recur_item > *rb_stack;
583  int derivs;
584  bool adapt, extrapolate, inited;
585  };
586 
589  struct c2_sample_recur {
590  c2_fblock<float_type> *f0, *f1;
591  float_type abs_tol, rel_tol, dx_tolerance, abs_tol_min;
592  int derivs;
593  c2_piecewise_function_p<float_type> *out;
594  std::vector<float_type> *xvals, *yvals;
595  std::vector< recur_item > *rb_stack;
596  bool inited;
597  };
598 
601  struct c2_root_info {
602  c2_fblock<float_type> lower, upper;
603  bool inited;
604  };
605 
611  float_type integrate_step(struct c2_integrate_recur &rb)
612  const /* throw(c2_exception) */;
613 
619  void sample_step(struct c2_sample_recur &rb) const /* throw(c2_exception) */;
620 
627  mutable struct c2_root_info *root_info;
628 
629  mutable size_t owner_count;
630 };
631 
638 template <typename float_type=double> class c2_classic_function_p
639  : public c2_function<float_type> {
640 public:
643  c2_classic_function_p(const float_type (*c_func)(float_type))
644  : c2_function<float_type>(), func(c_func) {}
645 
648  virtual float_type value_with_derivatives(
649  float_type x, float_type *yprime,
650  float_type *yprime2) const /* throw(c2_exception) */
651  {
652  if(!func)
653  throw c2_exception("c2_classic_function called with null function");
654  if(yprime) *yprime=0;
655  if(yprime2) *yprime2=0;
656  return func(x);
657  }
658  virtual ~c2_classic_function_p() { }
659 
660 protected:
662  const float_type (*func)(float_type);
663 };
664 
680 template <typename float_type> class c2_const_ptr {
681 public:
683  c2_const_ptr() : func(0) {}
686  c2_const_ptr(const c2_function<float_type> &f) : func(0)
687  { this->set_function(&f); }
690  c2_const_ptr(const c2_const_ptr<float_type> &src) : func(0)
691  { this->set_function(src.get_ptr()); }
692  void set_function(const c2_function<float_type> *f)
693  {
694  if(func) func->release_ownership();
695  func=f;
696  if(func) func->claim_ownership();
697  }
698 
701  const c2_const_ptr<float_type> & operator =
702  (const c2_const_ptr<float_type> &f)
703  { this->set_function(f.get_ptr()); return f; }
706  const c2_function<float_type> & operator =
707  (const c2_function<float_type> &f)
708  { this->set_function(&f); return f; }
719  void release_for_return() /* throw(c2_exception) */
720  {
721  if(func) func->release_ownership_for_return();
722  func=0;
723  }
729  void unset_function(void) { this->set_function(0); }
731  ~c2_const_ptr() { this->set_function(0); }
732 
734  inline const c2_function<float_type> &get() const /* throw(c2_exception) */
735  {
736  if(!func) throw c2_exception("c2_ptr accessed uninitialized");
737  return *func;
738  }
740  inline const c2_function<float_type> *get_ptr() const { return func; }
742  inline const c2_function<float_type> *operator -> () const
743  { return &get(); }
745  bool valid() const { return func != 0; }
746 
749  operator const c2_function<float_type>& () const { return this->get(); }
750 
754 
755  float_type operator()(float_type x) const /* throw(c2_exception) */
756  { return get()(x); }
766  float_type operator()(float_type x, float_type *yprime,
767  float_type *yprime2) const /* throw(c2_exception) */
768  { return get().value_with_derivatives(x, yprime, yprime2); }
773  c2_sum_p<float_type> &operator + (const c2_function<float_type> &rhs)
774  const /* throw(c2_exception) */
775  { return *new c2_sum_p<float_type>(get(), rhs); }
776  c2_diff_p<float_type> &operator - (const c2_function<float_type> &rhs)
777  const /* throw(c2_exception) */
778  { return *new c2_diff_p<float_type>(get(), rhs); }
779  c2_product_p<float_type> &operator * (const c2_function<float_type> &rhs)
780  const /* throw(c2_exception) */
781  { return *new c2_product_p<float_type>(get(), rhs); }
782  c2_ratio_p<float_type> &operator / (const c2_function<float_type> &rhs)
783  const /* throw(c2_exception) */
784  { return *new c2_ratio_p<float_type>(get(), rhs); }
788  c2_composed_function_p<float_type> &
789  operator ()(const c2_function<float_type> &inner)
790  const /* throw(c2_exception) */
791  { return *new c2_composed_function_p<float_type>(get(), inner); }
792 
793 protected:
794  const c2_function<float_type> * func;
795 };
796 
797 template <typename float_type> class c2_ptr : public c2_const_ptr<float_type >
798 {
799 public:
801  c2_ptr() : c2_const_ptr<float_type>() {}
804  c2_ptr(c2_function<float_type> &f) :
805  c2_const_ptr<float_type>() { this->set_function(&f); }
808  c2_ptr(const c2_ptr<float_type> &src) :
809  c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
811  inline c2_function<float_type> &get() const /* throw(c2_exception) */
812  { return *const_cast<c2_function<float_type>*>(
813  &c2_const_ptr<float_type>::get()); }
815  inline c2_function<float_type> *get_ptr() const
816  { return const_cast<c2_function<float_type>*>(this->func); }
818  inline c2_function<float_type> *operator -> () const
819  { return &get(); }
822  const c2_ptr<float_type> & operator =(const c2_ptr<float_type> &f)
823  { this->set_function(f.get_ptr()); return f; }
826  c2_function<float_type> & operator =(c2_function<float_type> &f)
827  { this->set_function(&f); return f; }
828 private:
830  void operator =(const c2_const_ptr<float_type> &) { }
832  void operator =(const c2_function<float_type> &) { }
833 };
834 
835 template <typename float_type, template <typename> class c2_class >
836 class c2_typed_ptr
837  : public c2_const_ptr<float_type> {
838 public:
840  c2_typed_ptr() : c2_ptr<float_type>() {}
843  c2_typed_ptr(c2_class<float_type> &f)
844  : c2_const_ptr<float_type>() { this->set_function(&f); }
847  c2_typed_ptr(const c2_typed_ptr<float_type, c2_class> &src)
848  : c2_const_ptr<float_type>() { this->set_function(src.get_ptr()); }
849 
851  inline c2_class<float_type> &get() const /* throw(c2_exception) */
852  {
853  return *static_cast<c2_class<float_type> *>
854  (const_cast<c2_function<float_type>*>(&c2_const_ptr<float_type>::get()));
855  }
857  inline c2_class<float_type> *operator -> () const
858  { return &get(); }
860  inline c2_class<float_type> *get_ptr() const
861  { return static_cast<c2_class<float_type> *>(
862  const_cast<c2_function<float_type>*>(this->func)); }
863 
864  operator c2_class<float_type>&() const { return get(); }
867  void operator =(const c2_typed_ptr<float_type, c2_class> &f)
868  { this->set_function(f.get_ptr()); }
871  void operator =(c2_class<float_type> &f)
872  { this->set_function(&f); }
873 private:
874  void operator =(const c2_const_ptr<float_type> &) { }
875  void operator =(const c2_function<float_type> &) { }
876 };
877 
878 template <typename float_type=double> class c2_plugin_function_p :
879  public c2_function<float_type> {
880 public:
882  c2_plugin_function_p() : c2_function<float_type>(), func() {}
884  c2_plugin_function_p(c2_function<float_type> &f) :
885  c2_function<float_type>(),func(f) { }
886 
887  void set_function(c2_function<float_type> *f)
888  {
889  func.set_function(f);
890  if(f) this->set_domain(f->xmin(), f->xmax());
891  }
894  virtual float_type value_with_derivatives(
895  float_type x, float_type *yprime,
896  float_type *yprime2) const /* throw(c2_exception) */
897  {
898  if(!func.valid())
899  throw c2_exception("c2_plugin_function_p called uninitialized");
900  return func->value_with_derivatives(x, yprime, yprime2);
901  }
903  virtual ~c2_plugin_function_p() { }
904 
906  void unset_function() { func.unset_function(); }
907 
908  virtual void get_sampling_grid(float_type amin, float_type amax,
909  std::vector<float_type> &grid) const
910  {
911  if(!func.valid())
912  throw c2_exception("c2_plugin_function_p called uninitialized");
913  if(this->sampling_grid)
914  c2_function<float_type>::get_sampling_grid(amin, amax, grid);
915  else func->get_sampling_grid(amin, amax, grid);
916  }
917 protected:
918  c2_ptr<float_type> func;
919 };
920 
921 template <typename float_type=double> class c2_const_plugin_function_p
922  : public c2_plugin_function_p<float_type> {
923 public:
925  c2_const_plugin_function_p() : c2_plugin_function_p<float_type>() {}
927  c2_const_plugin_function_p(const c2_function<float_type> &f) :
928  c2_plugin_function_p<float_type>() { this->set_function(&f); }
929  void set_function(const c2_function<float_type> *f)
930  { c2_plugin_function_p<float_type>::set_function(
931  const_cast<c2_function<float_type>*>(f)); }
933  virtual ~c2_const_plugin_function_p() { }
934 
936  const c2_function<float_type> &get() const /* throw(c2_exception) */
937  { return this->func.get(); }
938 };
939 
940 template <typename float_type=double> class c2_binary_function
941  : public c2_function<float_type> {
942 public:
943  virtual float_type value_with_derivatives(
944  float_type x, float_type *yprime,
945  float_type *yprime2) const /* throw(c2_exception) */
946  {
947  if(stub)
948  throw c2_exception("attempt to evaluate a c2_binary_function stub");
949  return this->combine(*Left.get_ptr(), *Right.get_ptr(), x, yprime, yprime2);
950  }
951 
954  virtual ~c2_binary_function() { }
955 
956 protected:
957  c2_binary_function( float_type (*combiner)(
958  const c2_function<float_type> &left,
959  const c2_function<float_type> &right,
960  float_type x, float_type *yprime,
961  float_type *yprime2),
962  const c2_function<float_type> &left,
963  const c2_function<float_type> &right) :
964  c2_function<float_type>(), combine(combiner), Left(left),
965  Right(right), stub(false)
966  {
967  this->set_domain(
968  (left.xmin() > right.xmin()) ? left.xmin() : right.xmin(),
969  (left.xmax() < right.xmax()) ? left.xmax() : right.xmax()
970  );
971  }
972  c2_binary_function(float_type (*combiner)(
973  const c2_function<float_type> &left,
974  const c2_function<float_type> &right,
975  float_type x, float_type *yprime, float_type *yprime2)
976  ) : c2_function<float_type>(), combine(combiner),
977  Left(), Right(), stub(true) { }
978 
979 public:
980  float_type (* const combine)(
981  const c2_function<float_type> &left,
982  const c2_function<float_type> &right,
983  float_type x, float_type *yprime, float_type *yprime2);
984 
985 protected:
986  const c2_const_ptr<float_type> Left, Right;
987  bool stub;
988 };
989 
990 template <typename float_type=double> class c2_scaled_function_p
991  : public c2_function<float_type> {
992 public:
997  c2_scaled_function_p(const c2_function<float_type> &outer,
998  float_type scale) :
999  c2_function<float_type>(), func(outer), yscale(scale) { }
1000 
1003  void reset(float_type scale) { yscale=scale; }
1004 
1005  virtual float_type value_with_derivatives(
1006  float_type x, float_type *yprime,
1007  float_type *yprime2) const /* throw (c2_exception) */
1008  {
1009  float_type y=this->func->value_with_derivatives(x, yprime, yprime2);
1010  if(yprime) (*yprime)*=yscale;
1011  if(yprime2) (*yprime2)*=yscale;
1012  return y*yscale;
1013  }
1014 
1015 protected:
1016  c2_scaled_function_p<float_type>() : func() {}
1018  const c2_const_ptr<float_type> func;
1019  float_type yscale;
1020 };
1021 
1031 template <typename float_type=double> class c2_cached_function_p
1032  : public c2_function<float_type> {
1033 public:
1037  c2_cached_function_p(const c2_function<float_type> &f)
1038  : c2_function<float_type>(),
1039  func(f), init(false) {}
1042  virtual float_type value_with_derivatives(
1043  float_type x, float_type *yprime,
1044  float_type *yprime2) const /* throw(c2_exception) */
1045  {
1046  if(!init || x != x0) {
1047  y=this->func->value_with_derivatives(x, &yp, &ypp);
1048  x0=x;
1049  init=true;
1050  }
1051  if(yprime) *yprime=yp;
1052  if(yprime2) *yprime2=ypp;
1053  return y;
1054  }
1055 
1056 protected:
1057  c2_cached_function_p() : func() {}
1058  const c2_const_ptr<float_type> func;
1059  mutable bool init;
1060  mutable float_type x0, y, yp, ypp;
1061 
1062 };
1063 
1064 template <typename float_type=double> class c2_composed_function_p
1065  : public c2_binary_function<float_type> {
1066 public:
1067  c2_composed_function_p(const c2_function<float_type> &outer,
1068  const c2_function<float_type> &inner) :
1069  c2_binary_function<float_type>(combine, outer, inner) {
1070  this->set_domain(inner.xmin(), inner.xmax()); }
1072  c2_composed_function_p() : c2_binary_function<float_type>(combine) {}
1073 
1075  static float_type combine(const c2_function<float_type> &left,
1076  const c2_function<float_type> &right,
1077  float_type x, float_type *yprime,
1078  float_type *yprime2) /* throw(c2_exception) */
1079  {
1080  float_type y0, y1;
1081  if(yprime || yprime2) {
1082  float_type yp0, ypp0, yp1, ypp1;
1083  y0=right.value_with_derivatives(x, &yp0, &ypp0);
1084  y1=left.value_with_derivatives(y0, &yp1, &ypp1);
1085  if(yprime) *yprime=yp1*yp0;
1086  if(yprime2) *yprime2=ypp0*yp1+yp0*yp0*ypp1;
1087  } else {
1088  y0=right(x);
1089  y1=left(y0);
1090  }
1091  return y1;
1092  }
1093 };
1094 
1095 template <typename float_type=double> class c2_sum_p
1096  : public c2_binary_function<float_type> {
1097 public:
1101  c2_sum_p(const c2_function<float_type> &left,
1102  const c2_function<float_type> &right)
1103  : c2_binary_function<float_type>(combine, left, right) {}
1105  c2_sum_p() : c2_binary_function<float_type>(combine) {} ;
1106 
1108  static float_type combine(const c2_function<float_type> &left,
1109  const c2_function<float_type> &right,
1110  float_type x, float_type *yprime,
1111  float_type *yprime2) /* throw(c2_exception) */
1112  {
1113  float_type y0, y1;
1114  if(yprime || yprime2) {
1115  float_type yp0, ypp0, yp1, ypp1;
1116  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1117  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1118  if(yprime) *yprime=yp0+yp1;
1119  if(yprime2) *yprime2=ypp0+ypp1;
1120  } else {
1121  y0=left(x);
1122  y1=right(x);
1123  }
1124  return y0+y1;
1125  }
1126 };
1127 
1128 template <typename float_type=double> class c2_diff_p
1129  : public c2_binary_function<float_type> {
1130 public:
1134  c2_diff_p(const c2_function<float_type> &left,
1135  const c2_function<float_type> &right)
1136  : c2_binary_function<float_type>(combine, left, right) {}
1138  c2_diff_p() : c2_binary_function<float_type>(combine) {} ;
1139 
1141  static float_type combine(const c2_function<float_type> &left,
1142  const c2_function<float_type> &right,
1143  float_type x, float_type *yprime,
1144  float_type *yprime2) /* throw(c2_exception) */
1145  {
1146  float_type y0, y1;
1147  if(yprime || yprime2) {
1148  float_type yp0, ypp0, yp1, ypp1;
1149  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1150  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1151  if(yprime) *yprime=yp0-yp1;
1152  if(yprime2) *yprime2=ypp0-ypp1;
1153  } else {
1154  y0=left(x);
1155  y1=right(x);
1156  }
1157  return y0-y1;
1158  }
1159 };
1160 
1161 
1165 template <typename float_type=double> class c2_product_p
1166  : public c2_binary_function<float_type> {
1167 public:
1171  c2_product_p(const c2_function<float_type> &left,
1172  const c2_function<float_type> &right)
1173  : c2_binary_function<float_type>(combine, left, right) {}
1175  c2_product_p() : c2_binary_function<float_type>(combine) {} ;
1176 
1178  static float_type combine(const c2_function<float_type> &left,
1179  const c2_function<float_type> &right,
1180  float_type x, float_type *yprime,
1181  float_type *yprime2) /* throw(c2_exception) */
1182  {
1183  float_type y0, y1;
1184  if(yprime || yprime2) {
1185  float_type yp0, ypp0, yp1, ypp1;
1186  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1187  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1188  if(yprime) *yprime=y1*yp0+y0*yp1;
1189  if(yprime2) *yprime2=ypp0*y1+2.0*yp0*yp1+ypp1*y0;
1190  } else {
1191  y0=left(x);
1192  y1=right(x);
1193  }
1194  return y0*y1;
1195  }
1196 };
1197 
1198 
1202 template <typename float_type=double> class c2_ratio_p
1203  : public c2_binary_function<float_type> {
1204 public:
1208  c2_ratio_p(const c2_function<float_type> &left,
1209  const c2_function<float_type> &right)
1210  : c2_binary_function<float_type>(combine, left, right) {}
1212  c2_ratio_p() : c2_binary_function<float_type>(combine) {} ;
1213 
1215  static float_type combine(const c2_function<float_type> &left,
1216  const c2_function<float_type> &right,
1217  float_type x, float_type *yprime,
1218  float_type *yprime2) /* throw(c2_exception) */
1219  {
1220  float_type y0, y1;
1221  if(yprime || yprime2) {
1222  float_type yp0, ypp0, yp1, ypp1;
1223  y0=left.value_with_derivatives(x, &yp0, &ypp0);
1224  y1=right.value_with_derivatives(x, &yp1, &ypp1);
1225  if(yprime) *yprime=(yp0*y1-y0*yp1)/(y1*y1); // first deriv of ratio
1226  if(yprime2) *yprime2=(y1*y1*ypp0+y0*(2*yp1*yp1-y1*ypp1)-2*y1*yp0*yp1)
1227  /(y1*y1*y1);
1228  } else {
1229  y0=left(x);
1230  y1=right(x);
1231  }
1232  return y0/y1;
1233  }
1234 };
1235 
1240 template <typename float_type> class c2_constant_p
1241  : public c2_function<float_type> {
1242 public:
1243  c2_constant_p(float_type x) : c2_function<float_type>(), value(x) {}
1244  void reset(float_type val) { value=val; }
1245  virtual float_type value_with_derivatives(
1246  float_type, float_type *yprime,
1247  float_type *yprime2) const /* throw(c2_exception) */
1248  { if(yprime) *yprime=0; if(yprime2) *yprime2=0; return value; }
1249 
1250 private:
1251  float_type value;
1252 };
1253 
1256 template <typename float_type> class c2_transformation {
1257 public:
1264  c2_transformation(bool transformed,
1265  float_type (*xin)(float_type),
1266  float_type (*xinp)(float_type),
1267  float_type (*xinpp)(float_type),
1268  float_type (*xout)(float_type)
1269  ) :
1270  fTransformed(transformed), fHasStaticTransforms(true),
1271  pIn(xin), pInPrime(xinp), pInDPrime(xinpp), pOut(xout) { }
1272 
1276  c2_transformation(bool transformed) :
1277  fTransformed(transformed), fHasStaticTransforms(false),
1278  pIn(report_error), pInPrime(report_error), pInDPrime(report_error),
1279  pOut(report_error) { }
1281  virtual ~c2_transformation() { }
1283  const bool fTransformed;
1286  const bool fHasStaticTransforms;
1287 
1293  float_type (* const pIn)(float_type);
1295  float_type (* const pInPrime)(float_type);
1297  float_type (* const pInDPrime)(float_type);
1299  float_type (* const pOut)(float_type);
1300 
1302  virtual float_type fIn(float_type x) const { return pIn(x); }
1304  virtual float_type fInPrime(float_type x) const { return pInPrime(x); }
1306  virtual float_type fInDPrime(float_type x) const { return pInDPrime(x); }
1308  virtual float_type fOut(float_type x) const { return pOut(x); }
1309 
1310 protected:
1312  static float_type report_error(float_type x) {
1313  throw c2_exception("use of improperly constructed axis transform");
1314  return x; }
1316  static float_type ident(float_type x) { return x; }
1318  static float_type one(float_type) { return 1; }
1320  static float_type zero(float_type) { return 0; }
1322  static float_type recip(float_type x) { return 1.0/x; }
1324  static float_type recip_prime(float_type x) { return -1/(x*x); }
1326  static float_type recip_prime2(float_type x) { return 2/(x*x*x); }
1327 };
1328 
1331 template <typename float_type> class c2_transformation_linear
1332  : public c2_transformation<float_type> {
1333 public:
1335  c2_transformation_linear() : c2_transformation<float_type>(
1336  false, this->ident,
1337  this->one, this->zero,
1338  this->ident) { }
1340  virtual ~c2_transformation_linear() { }
1341 };
1344 template <typename float_type> class c2_transformation_log
1345  : public c2_transformation<float_type> {
1346 public:
1348  c2_transformation_log() : c2_transformation<float_type>(
1349  true, std::log, this->recip,
1350  this->recip_prime, std::exp) { }
1352  virtual ~c2_transformation_log() { }
1353 };
1356 template <typename float_type> class c2_transformation_recip
1357  : public c2_transformation<float_type> {
1358 public:
1360  c2_transformation_recip() : c2_transformation<float_type>(
1361  true, this->recip,
1362  this->recip_prime,
1363  this->recip_prime2, this->recip) { }
1365  virtual ~c2_transformation_recip() { }
1366 };
1367 
1377 template <typename float_type>
1378 class c2_function_transformation {
1379 public:
1383  c2_function_transformation(
1384  const c2_transformation<float_type> &xx,
1385  const c2_transformation<float_type> &yy) :
1386  isIdentity(!(xx.fTransformed || yy.fTransformed)), X(xx), Y(yy) { }
1388  virtual ~c2_function_transformation() { delete &X; delete &Y; }
1389  virtual float_type evaluate(
1390  float_type xraw,
1391  float_type y, float_type yp0, float_type ypp0,
1392  float_type *yprime, float_type *yprime2) const;
1393  const bool isIdentity;
1395  const c2_transformation<float_type> &X;
1397  const c2_transformation<float_type> &Y;
1398 };
1399 
1403 template <typename float_type> class c2_lin_lin_function_transformation :
1404  public c2_function_transformation<float_type> {
1405 public:
1406  c2_lin_lin_function_transformation() :
1407  c2_function_transformation<float_type>(
1408  *new c2_transformation_linear<float_type>,
1409  *new c2_transformation_linear<float_type>
1410  ) { }
1411  virtual ~c2_lin_lin_function_transformation() { }
1412 };
1413 
1417 template <typename float_type> class c2_log_log_function_transformation :
1418  public c2_function_transformation<float_type> {
1419 public:
1420  c2_log_log_function_transformation() :
1421  c2_function_transformation<float_type>(
1422  *new c2_transformation_log<float_type>,
1423  *new c2_transformation_log<float_type>
1424  ) { }
1425  virtual ~c2_log_log_function_transformation() { }
1426 };
1427 
1431 template <typename float_type> class c2_lin_log_function_transformation :
1432  public c2_function_transformation<float_type> {
1433 public:
1434  c2_lin_log_function_transformation() :
1435  c2_function_transformation<float_type>(
1436  *new c2_transformation_linear<float_type>,
1437  *new c2_transformation_log<float_type>
1438  ) { }
1439  virtual ~c2_lin_log_function_transformation() { }
1440 };
1441 
1445 template <typename float_type> class c2_log_lin_function_transformation :
1446  public c2_function_transformation<float_type> {
1447 public:
1448  c2_log_lin_function_transformation() :
1449  c2_function_transformation<float_type>(
1450  *new c2_transformation_log<float_type>,
1451  *new c2_transformation_linear<float_type>
1452  ) { }
1453  virtual ~c2_log_lin_function_transformation() { }
1454 };
1455 
1460 template <typename float_type> class c2_arrhenius_function_transformation :
1461  public c2_function_transformation<float_type> {
1462 public:
1463  c2_arrhenius_function_transformation() :
1464  c2_function_transformation<float_type>(
1465  *new c2_transformation_recip<float_type>,
1466  *new c2_transformation_log<float_type>
1467  ) { }
1468  virtual ~c2_arrhenius_function_transformation() { }
1469 };
1470 
1512 template <typename float_type=double> class interpolating_function_p
1513  : public c2_function<float_type> {
1514 public:
1517  interpolating_function_p() : c2_function<float_type>(),
1518  fTransform(*new c2_lin_lin_function_transformation<float_type>) { }
1519 
1523  interpolating_function_p(const c2_function_transformation<float_type> &
1524  transform)
1525  : c2_function<float_type>(),
1526  fTransform(transform) { }
1527 
1550  interpolating_function_p<float_type> & load(
1551  const std::vector<float_type> &x,
1552  const std::vector<float_type> &f,
1553  bool lowerSlopeNatural, float_type lowerSlope,
1554  bool upperSlopeNatural, float_type upperSlope, bool splined=true
1555  ) /* throw(c2_exception) */;
1556 
1572  interpolating_function_p<float_type> & load_pairs(
1573  std::vector<std::pair<float_type, float_type> > &data,
1574  bool lowerSlopeNatural, float_type lowerSlope,
1575  bool upperSlopeNatural, float_type upperSlope, bool splined=true
1576  ) /* throw(c2_exception) */;
1577 
1612  interpolating_function_p<float_type> &
1613  sample_function(const c2_function<float_type> &func,
1614  float_type amin, float_type amax,
1615  float_type abs_tol, float_type rel_tol,
1616  bool lowerSlopeNatural, float_type lowerSlope,
1617  bool upperSlopeNatural, float_type upperSlope
1618  ) /* throw(c2_exception) */;
1619 
1644  interpolating_function_p<float_type> & load_random_generator_function(
1645  const std::vector<float_type> &bincenters,
1646  const c2_function<float_type> &binheights)
1647  /* throw(c2_exception) */;
1648 
1649  interpolating_function_p<float_type> & load_random_generator_bins(
1650  const std::vector<float_type> &bins,
1651  const std::vector<float_type> &binheights,
1652  bool splined=true)
1653  /* throw(c2_exception) */;
1654 
1655  virtual float_type value_with_derivatives(
1656  float_type x, float_type *yprime,
1657  float_type *yprime2) const /* throw(c2_exception) */;
1658 
1660  virtual ~interpolating_function_p() { delete &fTransform; }
1661 
1662  virtual interpolating_function_p<float_type> &clone()
1663  const /* throw(c2_exception) */
1664  { return *new interpolating_function_p<float_type>(); }
1665 
1666  void get_data(std::vector<float_type> &xvals,
1667  std::vector<float_type> &yvals) const /* throw() */ ;
1668 
1669  void get_internal_data(
1670  std::vector<float_type> &xvals,
1671  std::vector<float_type> &yvals,
1672  std::vector<float_type> &y2vals) const
1673  { xvals=X; yvals=F; y2vals=y2; }
1674 
1675  void set_lower_extrapolation(float_type bound);
1676  void set_upper_extrapolation(float_type bound);
1677 
1678  interpolating_function_p <float_type> &
1679  unary_operator(const c2_function<float_type> &source) const;
1680 
1681  interpolating_function_p <float_type> &
1682  binary_operator(const c2_function<float_type> &rhs,
1683  const c2_binary_function<float_type> *combining_stub
1684  ) const;
1685  interpolating_function_p <float_type> &
1686  add_pointwise (const c2_function<float_type> &rhs) const {
1687  return binary_operator(rhs, new c2_sum_p<float_type>()); }
1688  interpolating_function_p <float_type> &
1689  subtract_pointwise (const c2_function<float_type> &rhs) const {
1690  return binary_operator(rhs, new c2_diff_p<float_type>()); }
1691  interpolating_function_p <float_type> &
1692  multiply_pointwise (const c2_function<float_type> &rhs) const {
1693  return binary_operator(rhs, new c2_product_p<float_type>()); }
1694  interpolating_function_p <float_type> &
1695  divide_pointwise (const c2_function<float_type> &rhs) const {
1696  return binary_operator(rhs, new c2_ratio_p<float_type>()); }
1697  void clone_data(const interpolating_function_p <float_type> &rhs) {
1698  Xraw=rhs.Xraw; X=rhs.X; F=rhs.F; y2=rhs.y2;
1699  set_sampling_grid_pointer(Xraw);
1700  }
1701 
1702  const c2_function_transformation<float_type> &fTransform;
1703 
1704 protected:
1706  void spline(
1707  bool lowerSlopeNatural, float_type lowerSlope,
1708  bool upperSlopeNatural, float_type upperSlope
1709  ) /* throw(c2_exception) */;
1710 
1711  static bool comp_pair(std::pair<float_type,float_type> const &i,
1712  std::pair<float_type,float_type> const &j)
1713  {return i.first<j.first;}
1714 
1715  std::vector<float_type> Xraw, X, F, y2;
1716  c2_const_ptr<float_type> sampler_function;
1717  bool xInverted;
1718  mutable size_t lastKLow;
1719 };
1720 
1724 template <typename float_type=double> class log_lin_interpolating_function_p
1725  : public interpolating_function_p <float_type> {
1726 public:
1729  log_lin_interpolating_function_p() :
1730  interpolating_function_p<float_type>(
1731  *new c2_log_lin_function_transformation<float_type>)
1732  { }
1733  virtual interpolating_function_p<float_type> &clone()
1734  const /* throw(c2_exception) */
1735  { return *new log_lin_interpolating_function_p<float_type>(); }
1736 };
1737 
1738 
1743 template <typename float_type=double> class lin_log_interpolating_function_p
1744  : public interpolating_function_p <float_type> {
1745 public:
1748  lin_log_interpolating_function_p()
1749  : interpolating_function_p<float_type>(
1750  *new c2_lin_log_function_transformation<float_type>)
1751  { }
1752  virtual interpolating_function_p<float_type> &clone()
1753  const /* throw(c2_exception) */
1754  { return *new lin_log_interpolating_function_p<float_type>(); }
1755 };
1756 
1757 
1763 template <typename float_type=double> class log_log_interpolating_function_p
1764  : public interpolating_function_p <float_type> {
1765 public:
1768  log_log_interpolating_function_p() :
1769  interpolating_function_p<float_type>(
1770  *new c2_log_log_function_transformation<float_type>)
1771  { }
1772  virtual interpolating_function_p<float_type> &clone()
1773  const /* throw(c2_exception) */
1774  { return *new log_log_interpolating_function_p<float_type>(); }
1775 };
1776 
1777 
1783 template <typename float_type=double> class arrhenius_interpolating_function_p
1784  : public interpolating_function_p <float_type> {
1785 public:
1788  arrhenius_interpolating_function_p()
1789  : interpolating_function_p<float_type>(
1790  *new c2_arrhenius_function_transformation<float_type>)
1791  { }
1792  virtual interpolating_function_p<float_type> &clone()
1793  const /* throw(c2_exception) */
1794  { return *new arrhenius_interpolating_function_p<float_type>(); }
1795 };
1796 
1801 template <typename float_type=double> class c2_sin_p
1802  : public c2_function<float_type> {
1803 public:
1805  c2_sin_p() : c2_function<float_type>() {}
1806 
1807  virtual float_type value_with_derivatives(
1808  float_type x, float_type *yprime,
1809  float_type *yprime2) const /* throw(c2_exception) */
1810  { float_type q=std::sin(x);
1811  if(yprime) *yprime=std::cos(x);
1812  if(yprime2) *yprime2=-q;
1813  return q; }
1814 
1815  virtual void get_sampling_grid(float_type amin, float_type amax,
1816  std::vector<float_type> &grid) const;
1817 };
1818 
1823 template <typename float_type=double> class c2_cos_p
1824  : public c2_sin_p<float_type> {
1825 public:
1827  c2_cos_p() : c2_sin_p<float_type>() {}
1828 
1829  virtual float_type value_with_derivatives(
1830  float_type x, float_type *yprime,
1831  float_type *yprime2) const /* throw(c2_exception) */
1832  { float_type q=std::cos(x);
1833  if(yprime) *yprime=-std::sin(x);
1834  if(yprime2) *yprime2=-q;
1835  return q; }
1836 };
1837 
1842 template <typename float_type=double> class c2_tan_p
1843  : public c2_function<float_type> {
1844 public:
1846  c2_tan_p() : c2_function<float_type>() {}
1847 
1848  virtual float_type value_with_derivatives(
1849  float_type x, float_type *yprime,
1850  float_type *yprime2) const /* throw(c2_exception) */
1851  {
1852  float_type c=std::cos(x), ss=std::sin(x);
1853  float_type t=ss/c;
1854  float_type yp=1/(c*c);
1855  if(yprime) { *yprime=yp; }
1856  if(yprime2){*yprime2=2*t*yp; }
1857  return t;
1858  }
1859 };
1860 
1865 template <typename float_type=double> class c2_log_p
1866  : public c2_function<float_type> {
1867 public:
1869  c2_log_p() : c2_function<float_type>() {}
1870 
1871  virtual float_type value_with_derivatives(
1872  float_type x, float_type *yprime,
1873  float_type *yprime2) const /* throw(c2_exception) */
1874  { if(yprime) *yprime=1.0/x;
1875  if(yprime2) *yprime2=-1.0/(x*x);
1876  return std::log(x); }
1877 };
1878 
1883 template <typename float_type=double> class c2_exp_p
1884  : public c2_function<float_type> {
1885 public:
1887  c2_exp_p() : c2_function<float_type>() {}
1888 
1889  virtual float_type value_with_derivatives(
1890  float_type x, float_type *yprime,
1891  float_type *yprime2) const /* throw(c2_exception) */
1892  { float_type q=std::exp(x);
1893  if(yprime) *yprime=q;
1894  if(yprime2) *yprime2=q;
1895  return q; }
1896 };
1897 
1902 template <typename float_type=double> class c2_sqrt_p
1903  : public c2_function<float_type> {
1904 public:
1906  c2_sqrt_p() : c2_function<float_type>() {}
1907 
1908  virtual float_type value_with_derivatives(
1909  float_type x, float_type *yprime,
1910  float_type *yprime2) const /* throw(c2_exception) */
1911  { float_type q=std::sqrt(x);
1912  if(yprime) *yprime=0.5/q;
1913  if(yprime2) *yprime2=-0.25/(x*q);
1914  return q; }
1915 };
1916 
1921 template <typename float_type=double> class c2_recip_p
1922  : public c2_function<float_type> {
1923 public:
1925  c2_recip_p(float_type scale) : c2_function<float_type>(), rscale(scale) {}
1926 
1927  virtual float_type value_with_derivatives(
1928  float_type x, float_type *yprime,
1929  float_type *yprime2) const /* throw(c2_exception) */
1930  {
1931  float_type q=1.0/x;
1932  float_type y=rscale*q;
1933  if(yprime) *yprime=-y*q;
1934  if(yprime2) *yprime2=2*y*q*q;
1935  return y;
1936  }
1939  void reset(float_type scale) { rscale=scale; }
1940 private:
1941  float_type rscale;
1942 };
1943 
1948 template <typename float_type=double> class c2_identity_p
1949  : public c2_function<float_type> {
1950 public:
1952  c2_identity_p() : c2_function<float_type>() {}
1953 
1954  virtual float_type value_with_derivatives(
1955  float_type x, float_type *yprime,
1956  float_type *yprime2) const /* throw(c2_exception) */
1957  { if(yprime) *yprime=1.0; if(yprime2) *yprime2=0; return x; }
1958 };
1959 
1971 template <typename float_type=double> class c2_linear_p
1972  : public c2_function<float_type> {
1973 public:
1978  c2_linear_p(float_type x0, float_type y0, float_type slope) :
1979  c2_function<float_type>(), xint(x0), intercept(y0), m(slope) {}
1984  void reset(float_type x0, float_type y0, float_type slope)
1985  { xint=x0; intercept=y0; m=slope; }
1986  virtual float_type value_with_derivatives(
1987  float_type x, float_type *yprime,
1988  float_type *yprime2) const /* throw(c2_exception) */
1989  { if(yprime) *yprime=m;
1990  if(yprime2) *yprime2=0;
1991  return m*(x-xint)+intercept; }
1992 
1993 private:
1994  float_type xint, intercept, m;
1995 protected:
1996  c2_linear_p() {}
1997 };
1998 
2013 template <typename float_type=double> class c2_quadratic_p
2014  : public c2_function<float_type> {
2015 public:
2021  c2_quadratic_p(float_type x0, float_type y0,
2022  float_type xcoef, float_type x2coef) :
2023  c2_function<float_type>(), intercept(y0), center(x0), a(x2coef), b(xcoef) {}
2029  void reset(float_type x0, float_type y0, float_type xcoef,
2030  float_type x2coef) { intercept=y0; center=x0; a=x2coef; b=xcoef; }
2031  virtual float_type value_with_derivatives(
2032  float_type x, float_type *yprime,
2033  float_type *yprime2) const /* throw(c2_exception) */
2034  { float_type dx=x-center;
2035  if(yprime) *yprime=2*a*dx+b;
2036  if(yprime2) *yprime2=2*a;
2037  return a*dx*dx+b*dx+intercept; }
2038 
2039 private:
2040  float_type intercept, center, a, b;
2041 protected:
2042  c2_quadratic_p() {}
2043 };
2044 
2057 template <typename float_type=double> class c2_power_law_p
2058  : public c2_function<float_type> {
2059 public:
2063  c2_power_law_p(float_type scale, float_type power) :
2064  c2_function<float_type>(), a(scale), b(power) {}
2068  void reset(float_type scale, float_type power) { a=scale; b=power; }
2069  virtual float_type value_with_derivatives(
2070  float_type x, float_type *yprime,
2071  float_type *yprime2) const /* throw(c2_exception) */
2072  { float_type q=a*std::pow(x,b-2);
2073  if(yprime) *yprime=b*q*x;
2074  if(yprime2) *yprime2=b*(b-1)*q;
2075  return q*x*x; }
2076 
2077 private:
2078  float_type a, b;
2079 protected:
2080  c2_power_law_p() {}
2081 };
2082 
2100 template <typename float_type=double> class c2_inverse_function_p
2101  : public c2_function<float_type> {
2102 public:
2105  c2_inverse_function_p(const c2_function<float_type> &source);
2106  virtual float_type value_with_derivatives(
2107  float_type x, float_type *yprime,
2108  float_type *yprime2) const /* throw(c2_exception) */;
2109 
2113  void set_start_hint(float_type hint) const { start_hint=hint; }
2114 
2115  virtual float_type get_start_hint(float_type x) const
2116  { return hinting_function.valid()? hinting_function(x) : start_hint; }
2117 
2143  void set_hinting_function(const c2_function<float_type> *hint_func)
2144  { hinting_function.set_function(hint_func); }
2149  void set_hinting_function(const c2_const_ptr<float_type> hint_func)
2150  { hinting_function=hint_func; }
2151 
2152 protected:
2153  c2_inverse_function_p() {}
2154  mutable float_type start_hint;
2155  const c2_const_ptr<float_type> func;
2156  c2_const_ptr<float_type> hinting_function;
2157 };
2158 
2172 template <typename float_type=double> class accumulated_histogram
2173  : public interpolating_function_p <float_type> {
2174 public:
2184  accumulated_histogram(const std::vector<float_type>binedges,
2185  const std::vector<float_type> binheights,
2186  bool normalize=false,
2187  bool inverse_function=false, bool drop_zeros=true);
2188 
2189 };
2190 
2214 template <typename float_type=double> class c2_connector_function_p
2215  : public c2_function<float_type> {
2216 public:
2227  c2_connector_function_p(float_type x0, const c2_function<float_type> &f0,
2228  float_type x2, const c2_function<float_type> &f2,
2229  bool auto_center, float_type y1);
2245  c2_connector_function_p(
2246  float_type x0, float_type y0, float_type yp0, float_type ypp0,
2247  float_type x2, float_type y2, float_type yp2, float_type ypp2,
2248  bool auto_center, float_type y1);
2257  c2_connector_function_p(
2258  const c2_fblock<float_type> &fb0,
2259  const c2_fblock<float_type> &fb2,
2260  bool auto_center, float_type y1);
2261 
2263  virtual ~c2_connector_function_p();
2264  virtual float_type value_with_derivatives(
2265  float_type x, float_type *yprime,
2266  float_type *yprime2) const /* throw (c2_exception) */;
2267 protected:
2269  void init(
2270  const c2_fblock<float_type> &fb0,
2271  const c2_fblock<float_type> &fb2,
2272  bool auto_center, float_type y1);
2273 
2274  float_type fhinv, fy1, fa, fb, fc, fd, fe, ff;
2275 };
2276 
2302 template <typename float_type=double> class c2_piecewise_function_p
2303  : public c2_function<float_type> {
2304 public:
2306  c2_piecewise_function_p();
2308  virtual ~c2_piecewise_function_p();
2309  virtual float_type value_with_derivatives(
2310  float_type x, float_type *yprime,
2311  float_type *yprime2) const /* throw (c2_exception) */;
2324  void append_function(const c2_function<float_type> &func)
2325  /* throw (c2_exception) */;
2326 protected:
2327  std::vector<c2_const_ptr<float_type> > functions;
2328  mutable int lastKLow;
2329 };
2330 
2331 #include "c2_function.icc"
2332 
2333 #endif