Isis Developer Reference
FunctionTools.h
Go to the documentation of this file.
1#ifndef FunctionTools_h
2#define FunctionTools_h
8/* SPDX-License-Identifier: CC0-1.0 */
9
10
11#include <math.h>
12#include <cfloat>
13
14#include <QList>
15#include <QString>
16
17#include "Constants.h"
18#include "IException.h"
19#include "IString.h"
20
21namespace Isis {
38 public:
39
40
48 static QList <double> realLinearRoots(double coeffLinearTerm, double coeffConstTerm) {
49 double m=coeffLinearTerm, b=coeffConstTerm;
50
51 QList<double> roots;
52
53 //if the slope is zero there are either 0 or infinite roots
54 // for the present I see no need to handle the infinite roots situation more elegantly...
55 if (m==0.0)
56 return roots;
57
58 roots << -b/m;
59
60 return roots;
61 }
62
71 static QList <double> realQuadraticRoots(double coeffQuadTerm, double coeffLinearTerm,
72 double coeffConstTerm) {
73 double A=coeffQuadTerm, B=coeffLinearTerm, C=coeffConstTerm;
74 double disc,q,temp; //helpers on the way to a solution
75
76 if (A==0.0)
77 return realLinearRoots(coeffLinearTerm, coeffConstTerm);
78
79 QList<double> roots;
80
81 disc = B*B-4*A*C;
82 if (disc < 0) return roots; //no solution, return empty set
83 q = -0.5*(B + (B < 0 ? -1.0 : 1.0)*sqrt(disc));
84
85 if (q == 0) roots << q/A;
86 else {
87 roots << q/A;
88 //after the first root make sure there are no duplicates
89 temp = C/q;
90 if (!roots.contains(temp)) roots << temp;
91 }
92 return roots;
93 }
94
95
96
106 static QList <double> realCubicRoots(double coeffCubicTerm, double coeffQuadTerm,
107 double coeffLinearTerm, double coeffConstTerm) {
108 double a=coeffQuadTerm, b=coeffLinearTerm, c=coeffConstTerm;
109 double Q,R; //helpers on the way to a solution
110
111 //first verify this is really a cubic
112 if (coeffCubicTerm == 0.0)
113 return realQuadraticRoots(coeffQuadTerm,coeffLinearTerm,coeffConstTerm);
114
115
116 //algorithm wants the leading coefficient to be 1.0
117 if ( coeffCubicTerm != 1.0) {
118 a /= coeffCubicTerm;
119 b /= coeffCubicTerm;
120 c /= coeffCubicTerm;
121 }
122
123 QList<double> roots;
124
125 Q = (a*a - 3.0*b) / 9.0;
126 R = (2*a*a*a - 9.0*a*b + 27.0*c) / 54.0;
127 //cout << "Q R: " << Q << " " << R << "\n";
128
129 if ( a == 0.0 && b == 0.0 ) { //one simple root
130 roots << (c < 0 ? 1.0 : -1.0)*pow(fabs(c),1.0/3.0);
131 }
132 else if (R*R <= Q*Q*Q) { //there are three roots (one of them can be a double root)
133 double theta = acos(R/sqrt(Q*Q*Q)),temp;
134 //cout << "three roots, theta: " << theta << "\n";
135 Q = -2.0*sqrt(Q); //just done to save some FLOPs
136 roots << Q*cos(theta/3.0) - a /3.0;
137 //after the first root make sure there are no duplicates
138 temp = Q*cos((theta + TWOPI)/3.0) - a / 3.0;
139 if (!roots.contains(temp)) roots << temp;
140 temp = Q*cos((theta - TWOPI)/3.0) - a / 3.0;
141 if (!roots.contains(temp)) roots << temp;
142 }
143 else { //there is a single real root
144 double A, B;
145 A = (R < 0 ? 1.0 : -1.0 ) * pow(fabs(R) + sqrt(R*R - Q*Q*Q),1.0/3.0);
146 B = A == 0.0 ? 0.0 : Q / A;
147
148 roots << (A + B) - a / 3.0;
149 }
150
151 return roots;
152 }
153
175 template <typename Functor>
176 static bool brentsRootFinder(Functor &func, const QList<double> pt1, const QList<double> pt2,
177 double tol, int maxIter, double &root) {
178 double a = pt1[0], b = pt2[0], c, d = 0;
179 double fa = pt1[1], fb = pt2[1], fc;
180 double p, q, r, s, t, tol1, bnew, fbnew, temp, deltaI, deltaF;
181
182 // offset used for improved numerical stability
183 double offset = (a + b) / 2.0;
184 a -= offset;
185 b -= offset;
186
187 // check to see if the points bracket a root(s), if the signs are equal they don't
188 if ( (fa > 0) - (fa < 0) == (fb > 0) - (fb < 0) ) {
189 QString msg = "The function evaluations of two bounding points passed to Brents Method "
190 "have the same sign. Therefore, they don't necessary bound a root. No "
191 "root finding will be attempted.\n";
193 return false;
194 }
195
196 bool mflag = true;
197
198 if (fabs(fa) < fabs(fb)) {
199 // if a is a better guess for the root than b, switch them
200 // b is always the current best guess for the root
201 temp = a;
202 a = b;
203 b = temp;
204 temp = fa;
205 fa = fb;
206 fb = temp;
207 }
208
209 c = a;
210 fc = fa;
211
212 for (int iter = 0; iter < maxIter; iter++) {
213 tol1 = DBL_EPSILON * 2.0 * fabs(b); // numerical tolerance
214 if (a != c && b != c) {
215 // inverse quadratic interpolation
216 r = fb / fc;
217 s = fb / fa;
218 t = fa / fc;
219 p = s * (t * (r - t) * (c - b) - (1 - r) * (b - a));
220 q = (t - 1) * (r - 1) * (s - 1);
221 bnew = b + p / q;
222 }
223 else {
224 // secant rule
225 bnew = b - fb * (b - a) / (fb - fa);
226 }
227
228 // five tests follow to determine if the iterpolation methods are working better than
229 // bisection. p and q are setup as the bounds we want the new root guess to fall within.
230 // this enforces that the new root guess be withing the 3/4 of the interval closest to
231 // b, the current best guess.
232 temp = (3 * a + b) / 4.0;
233 p = temp < b ? temp : b;
234 q = b > temp ? b : temp;
235 deltaI = fabs(b - bnew); // magnitude of the interpolated correctio
236 if ( // if the root isn't within the 3/4 of the interval closest to b (the current best guess)
237 (bnew < p || bnew > q ) ||
238 // or if the last iteration was a bisection
239 // and the new correction is greater magnitude than half the magnitude of last
240 // correction, ie it's doing less to narrow the root than a bisection would
241 (mflag && deltaI >= fabs(b - c) / 2.0) ||
242 // or if the last iteration was an interpolation
243 // and the new correction magnitude is greater than
244 // the half the magnitude of the correction from two iterations ago,
245 // i.e. it's not converging faster than bisection
246 (!mflag && deltaI >= fabs(c - d) / 2.0) ||
247 // or if the last iteration was a bisection
248 // and the last correction was less than the numerical tolerance
249 // ie we are reaching the limits of our numerical
250 // ability to find a better root so lets do bisection, it's numerical safer
251 (mflag && fabs(b - c) < tol1) ||
252 // or it the last iteration was an interpolation
253 // and the correction from two iteration ago was less than the current
254 // numerical tolerance, ie we are reaching the limits of our numerical
255 // ability to find a better root so lets do bisection, it's numerical safer
256 (!mflag && fabs(c - d) < tol1) ) {
257
258 // bisection method
259 bnew = (a + b) / 2.0;
260 mflag = true;
261 }
262 else {
263 mflag = false;
264 }
265 try {
266 fbnew = func(bnew + offset);
267 } catch (IException &e) {
268 QString msg = "Function evaluation failed at:" + toString(bnew) +
269 ". Function must be continuous and defined for the entire interval "
270 "inorder to gaurentee brentsRootFinder will work.";
272 }
273 d = c; // thus d always equals the best guess from two iterations ago
274 c = b; // thus c always equals the best guess from the previous iteration
275 fc = fb;
276
277 if ( (fa > 0) - (fa < 0) == (fbnew > 0) - (fbnew < 0) ) {
278 // if b and bnew bracket the root
279 deltaF = fabs(a - bnew); // the Final magnitude of the correction
280 a = bnew;
281 fa = fbnew;
282 }
283 else {
284 // a and bnew bracket the root
285 deltaF = fabs(b - bnew); // the Final magnitude of the correction
286 b = bnew;
287 fb = fbnew;
288 }
289
290 if (fabs(fa) < fabs(fb)) {
291 // if a is a better guess for the root than b, switch them
292 temp = a;
293 a = b;
294 b = temp;
295 temp = fa;
296 fa = fb;
297 fb = temp;
298 }
299 if (fabs(fb) < tol) {
300 // if the tolerance is meet
301 root = b + offset;
302 return true;
303 }
304 else if ( deltaF < tol1 && fabs(b) < 100.0 * tol) {
305 // we've reach the numerical limit to how well the root can be defined, and the function
306 // is at least approaching zero
307 // This was added specifically for the apollo pan camera, the camera classes (and maybe
308 // NAIF as well can not actually converge to zero for the extreme edges of some pan
309 // images (partial derivates WRT to line approach infinity). They can get close
310 // 'enough' however
311 root = b + offset;
312 return true;
313 }
314 else if (deltaF < tol1) {
315 // we've reached the limit of the numerical ability to refine the root and the function
316 // is not near zero. This is a clasically ill defined root (nearly vertical function)
317 // "This is not the [root] you're looking for"
318 return false;
319 }
320 }
321 // maximum number of iteration exceeded
322 return false;
323 } // end brentsRootFinder
324
325
326 private:
333
339
340 }; // End FuntionTools class definition
341
342}; // End namespace Isis
343
344#endif
#define _FILEINFO_
Macro for the filename and line number.
Definition IException.h:24
A collection of tools for mathmatical function root finding, maximization, etc (eventually) This clas...
Definition FunctionTools.h:37
static bool brentsRootFinder(Functor &func, const QList< double > pt1, const QList< double > pt2, double tol, int maxIter, double &root)
Van Wijngaarden-Dekker-Brent Method for root finding on a discreetly defined function meaning that we...
Definition FunctionTools.h:176
static QList< double > realLinearRoots(double coeffLinearTerm, double coeffConstTerm)
Find the real roots (0 or 1) of a linear equation Form: coeffLinearTerm*X + coeffConstTerm = 0....
Definition FunctionTools.h:48
static QList< double > realQuadraticRoots(double coeffQuadTerm, double coeffLinearTerm, double coeffConstTerm)
The correct way to find the real roots of a quadratic (0, 1, or 2) (according to numerical recipies 3...
Definition FunctionTools.h:71
static QList< double > realCubicRoots(double coeffCubicTerm, double coeffQuadTerm, double coeffLinearTerm, double coeffConstTerm)
Find the real roots of a cubic (1, 2, or 3) (see numerical recipies 3rd edtion page 227) Form: coeffC...
Definition FunctionTools.h:106
Isis exception class.
Definition IException.h:91
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition IException.h:146
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition IString.cpp:211
const double TWOPI
Two * PI, a complete revolution.
Definition Constants.h:42