00001
00002
00003
00004 #ifndef ITER_TRACK_H_
00005 #define ITER_TRACK_H_
00006
00008
00009
00010
00011
00013
00014 double Amsterdam_Method( double (*f)(double), double a, double fa, double c, double fc,
00015 double tol, int max_iter, int& err );
00016
00017 class iter_track
00018 {
00019 vector< pair<double,double> > p_history;
00020 double p_result;
00021 double p_tol;
00022 int p_a;
00023 int p_b;
00024 int p_c;
00025 bool p_lgRootFound;
00026
00027 void p_clear0()
00028 {
00029 p_history.clear();
00030 }
00031 void p_clear1()
00032 {
00033 p_history.reserve( 10 );
00034 set_NaN( p_result );
00035 p_tol = numeric_limits<double>::max();
00036 p_a = -1;
00037 p_b = -1;
00038 p_c = -1;
00039 p_lgRootFound = false;
00040 }
00041 void p_set_root(double x)
00042 {
00043 p_result = x;
00044 p_lgRootFound = true;
00045 }
00046 double p_x(int ip) const
00047 {
00048 return p_history[ip].first;
00049 }
00050 double p_y(int ip) const
00051 {
00052 return p_history[ip].second;
00053 }
00054 double p_midpoint() const
00055 {
00056 return 0.5*(p_x(p_a) + p_x(p_c));
00057 }
00058 double p_numerator(double dab, double dcb, double fa, double fb, double fc)
00059 {
00060 return fb*(dab*fc*(fc-fb)-fa*dcb*(fa-fb));
00061 }
00062 double p_denominator(double fa, double fb, double fc)
00063 {
00064 return (fc-fb)*(fa-fb)*(fa-fc);
00065 }
00066
00067 public:
00068 iter_track()
00069 {
00070 p_clear1();
00071 }
00072 ~iter_track()
00073 {
00074 p_clear0();
00075 }
00076 void clear()
00077 {
00078 p_clear0();
00079 p_clear1();
00080 }
00081 void set_tol(double tol)
00082 {
00083 p_tol = tol;
00084 }
00085 double bracket_width() const
00086 {
00087 return p_x(p_c) - p_x(p_a);
00088 }
00089 bool lgConverged()
00090 {
00091 if( p_lgRootFound )
00092 return true;
00093 if( bracket_width() < p_tol )
00094 {
00095 p_result = p_midpoint();
00096 return true;
00097 }
00098 return false;
00099 }
00100 double root() const
00101 {
00102 return p_result;
00103 }
00104 int init_bracket( double x1, double fx1, double x2, double fx2 )
00105 {
00106
00107 int s1 = sign3(fx1);
00108 int test = s1*sign3(fx2);
00109 if( test > 0 )
00110 return -1;
00111 if( test == 0 )
00112 p_set_root( ( s1 == 0 ) ? x1 : x2 );
00113
00114 p_history.push_back( pair<double,double>(x1,fx1) );
00115 p_history.push_back( pair<double,double>(x2,fx2) );
00116 p_a = ( x1 < x2 ) ? 0 : 1;
00117 p_c = ( x1 < x2 ) ? 1 : 0;
00118 return 0;
00119 }
00120 void add( double x, double fx )
00121 {
00122 p_history.push_back( pair<double,double>(x,fx) );
00123 p_b = p_history.size()-1;
00124 if( fx == 0. )
00125 p_set_root( x );
00126 }
00127 double next_val();
00128 double next_val(double max_rel_step)
00129 {
00130 double next = next_val();
00131 double last = p_history.back().first;
00132 double rel_step = safe_div( next, last ) - 1.;
00133 rel_step = sign( min(abs(rel_step),abs(max_rel_step)), rel_step );
00134 return (1.+rel_step)*last;
00135 }
00136
00137
00138 double deriv(int n, double& sigma) const;
00139 double deriv(double& sigma) const
00140 {
00141 return deriv( p_history.size(), sigma );
00142 }
00143 double deriv(int n) const
00144 {
00145 double sigma;
00146 return deriv( n, sigma );
00147 }
00148 double deriv() const
00149 {
00150 double sigma;
00151 return deriv( p_history.size(), sigma );
00152 }
00153
00154
00155 double zero_fit(int n, double& sigma) const;
00156 double zero_fit(double& sigma) const
00157 {
00158 return zero_fit( p_history.size(), sigma );
00159 }
00160 double zero_fit(int n) const
00161 {
00162 double sigma;
00163 return zero_fit( n, sigma );
00164 }
00165 double zero_fit() const
00166 {
00167 double sigma;
00168 return zero_fit( p_history.size(), sigma );
00169 }
00170
00171 void print_status() const
00172 {
00173 dprintf( ioQQQ, "a %i %.15e %.15e\n", p_a, p_x(p_a), p_y(p_a) );
00174 dprintf( ioQQQ, "b %i %.15e %.15e\n", p_b, p_x(p_b), p_y(p_b) );
00175 dprintf( ioQQQ, "c %i %.15e %.15e\n", p_c, p_x(p_c), p_y(p_c) );
00176 }
00177 void print_history() const
00178 {
00179 fprintf(ioQQQ," Temperature C-H history\n");
00180 for( unsigned int i=0; i < p_history.size(); ++i )
00181 fprintf( ioQQQ, "%.15e %.15e\n", p_x(i), p_y(i) );
00182 }
00183 };
00184
00185 template<class T>
00186 class iter_track_basic
00187 {
00188 T p_lo_bound;
00189 T p_hi_bound;
00190 void p_clear1()
00191 {
00192 p_lo_bound = numeric_limits<T>::max();
00193 p_hi_bound = numeric_limits<T>::min();
00194 }
00195 public:
00196 iter_track_basic()
00197 {
00198 p_clear1();
00199 }
00200 void clear()
00201 {
00202 p_clear1();
00203 }
00204 T next_val( T current, T next_est )
00205 {
00206 if( next_est < current )
00207 p_hi_bound = current;
00208 else
00209 p_lo_bound = current;
00210 if( p_lo_bound < p_hi_bound )
00211 return (T)0.5*(p_lo_bound + p_hi_bound);
00212 else
00213 return next_est;
00214 }
00215 };
00216
00217 #endif