cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TestIterTrack.cpp
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 #include "cdstd.h"
4 #include <UnitTest++.h>
5 #include "cddefines.h"
6 #include "iter_track.h"
7 
8 namespace {
9  // test a function that would have gone into a limit cycle
10  TEST(IterTrackBasicFloat)
11  {
13  sys_float x = 1.f;
14  sys_float xnew = 2.f;
15  while( abs(x-xnew) > 2.f*FLT_EPSILON*abs(x) )
16  {
17  x = xnew;
18  // derivative at root is -1
19  xnew = track.next_val( x, 2.f/x );
20  }
21  CHECK( fp_equal( xnew, (sys_float)sqrt(2.) ) );
22  // now clear the tracker and try again
23  track.clear();
24  x = 1.f;
25  xnew = 2.f;
26  while( abs(x-xnew) > 2.f*FLT_EPSILON*abs(x) )
27  {
28  x = xnew;
29  // derivative at root is -1
30  xnew = track.next_val( x, 3.f/x );
31  }
32  CHECK( fp_equal( xnew, (sys_float)sqrt(3.) ) );
33  }
34 
35  // now test the same for doubles
36  TEST(IterTrackBasicDouble)
37  {
39  double x = 1.;
40  double xnew = 2.;
41  while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
42  {
43  x = xnew;
44  // derivative at root is -1
45  xnew = track.next_val( x, 2./x );
46  }
47  CHECK( fp_equal( xnew, sqrt(2.) ) );
48  }
49 
50  // test a function that would have diverged
51  TEST(IterTrackBasicUnstable)
52  {
54  double x = 1.;
55  double xnew = 2.;
56  while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
57  {
58  x = xnew;
59  // derivative at root is -5
60  xnew = track.next_val( x, 1./x - 2.*x );
61  }
62  // possible solutions are +/-sqrt(1/3), either one may be found
63  CHECK( fp_equal( abs(xnew), 1./sqrt(3.) ) );
64  }
65 
66  // test a function that would have converged anyway
67  TEST(IterTrackBasicStableNeg)
68  {
70  double x = 1.;
71  double xnew = 2.;
72  // this would have diverged without iter_tracking
73  while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
74  {
75  x = xnew;
76  // derivative at root is -1/3
77  xnew = track.next_val( x, 1./x + x/3. );
78  }
79  CHECK( fp_equal( xnew, sqrt(1.5) ) );
80  }
81 
82  TEST(IterTrackBasicStablePos)
83  {
85  double x = 1.;
86  double xnew = 2.;
87  // this would have diverged without iter_tracking
88  while( abs(x-xnew) > 2.*DBL_EPSILON*abs(x) )
89  {
90  x = xnew;
91  // derivative at root is 1/3
92  xnew = track.next_val( x, 1./x + 2.*x/3. );
93  }
94  CHECK( fp_equal( xnew, sqrt(3.) ) );
95  }
96 
97  double testfun(double x)
98  {
99  return sin(x)-0.5;
100  }
101 
102  TEST(IterTrack)
103  {
104  double x1, fx1, x2, fx2, x3, fx3;
105  iter_track track;
106  x1 = 0.;
107  fx1 = testfun(x1);
108  x3 = 1.5;
109  fx3 = testfun(x3);
110  double tol = 1.e-12;
111  track.set_tol(tol);
112  CHECK_EQUAL( 0, track.init_bracket(x1,fx1,x3,fx3) );
113  CHECK( fp_equal( track.bracket_width(), abs(x3-x1) ) );
114  CHECK( !track.lgConverged() );
115  x2 = 0.5*(x1+x3);
116  fx2 = testfun(x2);
117  track.add( x2, fx2 );
118  double xnew = track.next_val(0.01);
119  CHECK( fp_equal_tol( abs(xnew/x2-1.), 0.01, 1.e-12 ) );
120  const int navg = 5;
121  vector<double> xvals( navg ); // keep track of the last navg x-values
122  for( int i=0; i < 100 && !track.lgConverged(); ++i )
123  {
124  x2 = track.next_val();
125  fx2 = testfun(x2);
126  track.add( x2, fx2 );
127  // use xvals as circular buffer
128  xvals[i%navg] = x2;
129  }
130  CHECK( track.lgConverged() );
131  double exact_root = asin(0.5);
132  CHECK( fp_equal_tol( track.root(), exact_root, tol ) );
133  double sigma;
134  double val = track.deriv( navg, sigma );
135  double delta_lo = *min_element( xvals.begin(), xvals.end() ) - exact_root;
136  double delta_hi = *max_element( xvals.begin(), xvals.end() ) - exact_root;
137  CHECK( delta_lo < 0. );
138  CHECK( delta_hi > 0. );
139  // the exact derivative at the root is sqrt(3)/2 = 0.8660254...
140  // the exact 2nd derivative at the root is -1/2
141  double err_lo = -0.5*delta_lo;
142  double err_hi = -0.5*delta_hi;
143  CHECK( fp_bound( sqrt(3.)/2.+err_hi, val, sqrt(3.)/2.+err_lo ) );
144  // the tangent at the exact root is given by asin(0.5) + sqrt(3)/2*(x-x0)
145  // if we subtract that from the Taylor expansion of testfun we get:
146  // residual = -1/4*(x-x0)^2 + O((x-x0)^3), hence sigma should be less
147  // than the maximum absolute value of -1/4*(x-x0)^2 (the actual fit
148  // should run slightly closer to the maximum deviant value than the tangent).
149  CHECK( sigma < max( pow2(err_lo), pow2(err_hi) ) );
150  // ask for more points than are available to see if that is handled correctly
151  val = track.deriv( 200 );
152  double val2 = track.deriv();
153  CHECK( fp_equal( val, val2 ) );
154  val = track.deriv( 200, sigma );
155  double sigma2;
156  val = track.deriv( sigma2 );
157  CHECK( fp_equal( sigma, sigma2 ) );
158 
159  // now do the same thing for the zero_fit() methods...
160  val = track.zero_fit( navg, sigma );
161  // the exact root is asin(0.5) = 0.52359877...
162  CHECK( fp_equal_tol( exact_root, val, 2.*sigma ) );
163  val = track.zero_fit( 200 );
164  val2 = track.zero_fit();
165  CHECK( fp_equal( val, val2 ) );
166  val = track.zero_fit( 200, sigma );
167  val = track.zero_fit( sigma2 );
168  CHECK( fp_equal( sigma, sigma2 ) );
169  }
170 
171  // this is the short version of the above...
172  TEST(AmsterdamMethod)
173  {
174  double x1, fx1, x2, fx2;
175  x1 = 0.;
176  fx1 = testfun(x1);
177  x2 = 1.5;
178  fx2 = testfun(x2);
179  double tol = 1.e-12;
180  int err = -1;
181  double x = Amsterdam_Method( testfun, x1, fx1, x2, fx2, tol, 1000, err );
182  CHECK_EQUAL( 0, err );
183  CHECK( fp_equal_tol( x, asin(0.5), tol ) );
184  }
185 
186  // now test some unstable functions
187  double testfun2(double x)
188  {
189  return exp(x)-3.;
190  }
191 
192  // the derivative at the root is 3
193  TEST(AmsterdamMethod2)
194  {
195  double x1, fx1, x2, fx2;
196  x1 = 0.;
197  fx1 = testfun2(x1);
198  x2 = 3.;
199  fx2 = testfun2(x2);
200  double tol = 1.e-12;
201  int err = -1;
202  double x = Amsterdam_Method( testfun2, x1, fx1, x2, fx2, tol, 1000, err );
203  CHECK_EQUAL( 0, err );
204  CHECK( fp_equal_tol( x, log(3.), tol ) );
205  }
206 
207  double testfun3(double x)
208  {
209  return 1./x - 2.*x;
210  }
211 
212  // the derivative at the root is -4
213  TEST(AmsterdamMethod3)
214  {
215  double x1, fx1, x2, fx2;
216  x1 = 0.1;
217  fx1 = testfun3(x1);
218  x2 = 3.;
219  fx2 = testfun3(x2);
220  double tol = 1.e-12;
221  int err = -1;
222  double x = Amsterdam_Method( testfun3, x1, fx1, x2, fx2, tol, 1000, err );
223  CHECK_EQUAL( 0, err );
224  CHECK( fp_equal_tol( x, sqrt(0.5), tol ) );
225  }
226 }
static double x2[63]
double bracket_width() const
Definition: iter_track.h:85
static double x1[83]
double next_val()
Definition: iter_track.cpp:57
void add(double x, double fx)
Definition: iter_track.h:120
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:904
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
T next_val(T current, T next_est)
Definition: iter_track.h:250
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
int init_bracket(double x1, double fx1, double x2, double fx2)
Definition: iter_track.h:104
bool lgConverged()
Definition: iter_track.h:89
float sys_float
Definition: cddefines.h:127
long max(int a, long b)
Definition: cddefines.h:817
bool fp_bound(sys_float lo, sys_float x, sys_float hi, int n=3)
Definition: cddefines.h:927
T pow2(T a)
Definition: cddefines.h:981
double zero_fit(int n, double &sigma) const
Definition: iter_track.cpp:203
double deriv(int n, double &sigma) const
Definition: iter_track.cpp:184
void set_tol(double tol)
Definition: iter_track.h:81
double root() const
Definition: iter_track.h:100