cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iter_track.h
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 
4 #ifndef ITER_TRACK_H_
5 #define ITER_TRACK_H_
6 
8 // The class iter_track was derived from this original source: //
9 // http://www.mymathlib.com/roots/amsterdam.html //
10 // //
11 // The original code was heavily modified by: Peter van Hoof, ROB //
13 
14 double Amsterdam_Method( double (*f)(double), double a, double fa, double c, double fc,
15  double tol, int max_iter, int& err );
16 
18 {
19  vector< pair<double,double> > p_history;
20  double p_result;
21  double p_tol;
22  int p_a;
23  int p_b;
24  int p_c;
26 
27  void p_clear0()
28  {
29  p_history.clear();
30  }
31  void p_clear1()
32  {
33  p_history.reserve( 10 );
34  set_NaN( p_result );
36  p_a = -1;
37  p_b = -1;
38  p_c = -1;
39  p_lgRootFound = false;
40  }
41  void p_set_root(double x)
42  {
43  p_result = x;
44  p_lgRootFound = true;
45  }
46  double p_x(int ip) const
47  {
48  return p_history[ip].first;
49  }
50  double p_y(int ip) const
51  {
52  return p_history[ip].second;
53  }
54  double p_midpoint() const
55  {
56  return 0.5*(p_x(p_a) + p_x(p_c));
57  }
58  double p_numerator(double dab, double dcb, double fa, double fb, double fc)
59  {
60  return fb*(dab*fc*(fc-fb)-fa*dcb*(fa-fb));
61  }
62  double p_denominator(double fa, double fb, double fc)
63  {
64  return (fc-fb)*(fa-fb)*(fa-fc);
65  }
66 
67 public:
69  {
70  p_clear1();
71  }
73  {
74  p_clear0();
75  }
76  void clear()
77  {
78  p_clear0();
79  p_clear1();
80  }
81  void set_tol(double tol)
82  {
83  p_tol = tol;
84  }
85  double bracket_width() const
86  {
87  return p_x(p_c) - p_x(p_a);
88  }
89  bool lgConverged()
90  {
91  if( p_lgRootFound )
92  return true;
93  if( bracket_width() < p_tol )
94  {
95  p_result = p_midpoint();
96  return true;
97  }
98  return false;
99  }
100  double root() const
101  {
102  return p_result;
103  }
104  int init_bracket( double x1, double fx1, double x2, double fx2 )
105  {
106  // fx1 and fx2 must have opposite sign, or be zero
107  int s1 = sign3(fx1);
108  int test = s1*sign3(fx2);
109  if( test > 0 )
110  return -1;
111  if( test == 0 )
112  p_set_root( ( s1 == 0 ) ? x1 : x2 );
113 
114  p_history.push_back( pair<double,double>(x1,fx1) );
115  p_history.push_back( pair<double,double>(x2,fx2) );
116  p_a = ( x1 < x2 ) ? 0 : 1;
117  p_c = ( x1 < x2 ) ? 1 : 0;
118  return 0;
119  }
120  void add( double x, double fx )
121  {
122  p_history.push_back( pair<double,double>(x,fx) );
123  p_b = p_history.size()-1;
124  if( fx == 0. )
125  p_set_root( x );
126  }
127  double next_val();
128  double next_val(double max_rel_step)
129  {
130  double next = next_val();
131  double last = p_history.back().first;
132  double rel_step = safe_div( next, last ) - 1.;
133  rel_step = sign( min(abs(rel_step),abs(max_rel_step)), rel_step );
134  return (1.+rel_step)*last;
135  }
136  // these routines return a numerical estimate of the derivative
137  // by making a linear least squares fit of y(x) to the last n steps
138  double deriv(int n, double& sigma) const;
139  double deriv(double& sigma) const
140  {
141  return deriv( p_history.size(), sigma );
142  }
143  double deriv(int n) const
144  {
145  double sigma;
146  return deriv( n, sigma );
147  }
148  double deriv() const
149  {
150  double sigma;
151  return deriv( p_history.size(), sigma );
152  }
153  // these routines return a numerical estimate of the root
154  // by making a linear least squares fit of x(y) to the last n steps
155  double zero_fit(int n, double& sigma) const;
156  double zero_fit(double& sigma) const
157  {
158  return zero_fit( p_history.size(), sigma );
159  }
160  double zero_fit(int n) const
161  {
162  double sigma;
163  return zero_fit( n, sigma );
164  }
165  double zero_fit() const
166  {
167  double sigma;
168  return zero_fit( p_history.size(), sigma );
169  }
170  int in_bounds(double x) const
171  {
172  if( x < p_x(p_a) )
173  return -1;
174  else if( x > p_x(p_c) )
175  return 1;
176  else
177  return 0;
178  }
179  // the following two methods are debugging aids
180  void print_status() const
181  {
182  dprintf( ioQQQ, "a %i %.15e %.15e\n", p_a, p_x(p_a), p_y(p_a) );
183  dprintf( ioQQQ, "b %i %.15e %.15e\n", p_b, p_x(p_b), p_y(p_b) );
184  dprintf( ioQQQ, "c %i %.15e %.15e\n", p_c, p_x(p_c), p_y(p_c) );
185  }
186  void print_history() const
187  {
188  fprintf( ioQQQ, " x(i) y(i) iter_track history\n" );
189  for( unsigned int i=0; i < p_history.size(); ++i )
190  fprintf( ioQQQ, "%.15e %.15e\n", p_x(i), p_y(i) );
191  }
192 };
193 
194 //
229 //
230 template<class T>
232 {
235  void p_clear1()
236  {
237  // invalidate the bracket
240  }
241 public:
243  {
244  p_clear1();
245  }
246  void clear()
247  {
248  p_clear1();
249  }
250  T next_val( T current, T next_est )
251  {
252  // update the bounds of the bracket
253  if( next_est < current )
254  p_hi_bound = current;
255  else
256  p_lo_bound = current;
257  // if the bracket has been established, do a bisection
258  // otherwise simply return next_est.
259  if( p_lo_bound < p_hi_bound )
260  return (T)0.5*(p_lo_bound + p_hi_bound);
261  else
262  return next_est;
263  }
264  static const int PREV_ITER=2;
265 };
266 
268 {
272 
273  void p_clear1()
274  {
275  m_confidence = 1.0;
276  m_lastnext = 0.0;
277  m_lastcurr = 0.0;
278  }
279 public:
281  {
282  p_clear1();
283  }
284  void clear()
285  {
286  p_clear1();
287  }
288  realnum next_val( realnum current, realnum next_est );
289  static const int PREV_ITER=1;
290 };
291 
292 #endif
double p_midpoint() const
Definition: iter_track.h:54
double deriv(double &sigma) const
Definition: iter_track.h:139
realnum m_lastnext
Definition: iter_track.h:270
static double x2[63]
realnum m_lastcurr
Definition: iter_track.h:271
double bracket_width() const
Definition: iter_track.h:85
static double x1[83]
void set_NaN(sys_float &x)
Definition: cpu.cpp:906
double next_val()
Definition: iter_track.cpp:57
void add(double x, double fx)
Definition: iter_track.h:120
void print_history() const
Definition: iter_track.h:186
double p_x(int ip) const
Definition: iter_track.h:46
double p_denominator(double fa, double fb, double fc)
Definition: iter_track.h:62
double deriv(int n) const
Definition: iter_track.h:143
void p_clear0()
Definition: iter_track.h:27
T sign(T x, T y)
Definition: cddefines.h:842
double zero_fit(int n) const
Definition: iter_track.h:160
realnum m_confidence
Definition: iter_track.h:269
static const int PREV_ITER
Definition: iter_track.h:264
bool p_lgRootFound
Definition: iter_track.h:25
double Amsterdam_Method(double(*f)(double), double a, double fa, double c, double fc, double tol, int max_iter, int &err)
Definition: iter_track.cpp:277
FILE * ioQQQ
Definition: cddefines.cpp:7
T next_val(T current, T next_est)
Definition: iter_track.h:250
double p_y(int ip) const
Definition: iter_track.h:50
realnum next_val(realnum current, realnum next_est)
Definition: iter_track.cpp:19
double p_result
Definition: iter_track.h:20
vector< pair< double, double > > p_history
Definition: iter_track.h:19
int init_bracket(double x1, double fx1, double x2, double fx2)
Definition: iter_track.h:104
static const int PREV_ITER
Definition: iter_track.h:289
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1102
void clear()
Definition: iter_track.h:284
void p_clear1()
Definition: iter_track.h:31
double zero_fit() const
Definition: iter_track.h:165
bool lgConverged()
Definition: iter_track.h:89
float realnum
Definition: cddefines.h:124
long max(int a, long b)
Definition: cddefines.h:817
int sign3(T x)
Definition: cddefines.h:850
long min(int a, long b)
Definition: cddefines.h:762
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1011
double zero_fit(double &sigma) const
Definition: iter_track.h:156
void print_status() const
Definition: iter_track.h:180
double p_tol
Definition: iter_track.h:21
void clear()
Definition: iter_track.h:76
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
void p_clear1()
Definition: iter_track.h:273
int in_bounds(double x) const
Definition: iter_track.h:170
void p_set_root(double x)
Definition: iter_track.h:41
double deriv() const
Definition: iter_track.h:148
void set_tol(double tol)
Definition: iter_track.h:81
double p_numerator(double dab, double dcb, double fa, double fb, double fc)
Definition: iter_track.h:58
double root() const
Definition: iter_track.h:100
double next_val(double max_rel_step)
Definition: iter_track.h:128