Isis 3 Programmer Reference
ZeroBufferFit.cpp
1 
7 /* SPDX-License-Identifier: CC0-1.0 */
8 
9 #include <cmath>
10 #include <string>
11 #include <vector>
12 #include <numeric>
13 #include <iostream>
14 #include <sstream>
15 
16 #include "ZeroBufferFit.h"
17 #include "IString.h"
18 #include "HiCalTypes.h"
19 #include "HiCalUtil.h"
20 #include "HiCalConf.h"
21 #include "LeastSquares.h"
22 #include "LowPassFilter.h"
23 #include "MultivariateStatistics.h"
24 #include "IException.h"
25 
26 using namespace std;
27 
28 namespace Isis {
29 
38  ZeroBufferFit::ZeroBufferFit(const HiCalConf &conf) :
39  NonLinearLSQ(), Module("ZeroBufferFit") {
40  DbProfile prof = conf.getMatrixProfile();
41  _history.add("Profile["+ prof.Name()+"]");
42  _timet.setBin(ToInteger(prof("Summing")));
43  _timet.setLineTime(ToDouble(prof("ScanExposureDuration")));
44 
45  _skipFit = IsEqual(ConfKey(prof, "ZeroBufferFitSkipFit", QString("TRUE")), "TRUE");
46  _useLinFit = IsTrueValue(prof, "ZeroBufferFitOnFailUseLinear");
47 
48 
49  _absErr = toDouble(ConfKey(prof, "AbsoluteError", QString("1.0E-4")));
50  _relErr = toDouble(ConfKey(prof, "RelativeError", QString("1.0E-4")));
51 
52  _sWidth = toInt(ConfKey(prof, "GuessFilterWidth", QString("17")));
53  _sIters = toInt(ConfKey(prof, "GuessFilterIterations", QString("1")));
54 
55  if ( prof.exists("MaximumIterations") ) {
56  setMaxIters(ToInteger(prof("MaximumIterations")));
57  }
58 
59  _maxLog = toDouble(ConfKey(prof, "MaximumLog", QString("709.0")));
60  _badLines = ToInteger(prof("TrimLines"))/ToInteger(prof("Summing"));
61  _minLines = toInt(ConfKey(prof,"ZeroBufferFitMinimumLines", QString("100")));
62 
63 
64  QString histstr = "ZeroBufferFit(AbsErr[" + ToString(_absErr) +
65  "],RelErr[" + ToString(_relErr) +
66  "],MaxIter[" + ToString(maxIters()) + "])";
67  _history.add(histstr);
68 
69  }
70 
91  HiVector ZeroBufferFit::Solve(const HiVector &d) {
92  ostringstream hist;
93  _data = d;
94  if ( _skipFit || (!gotGoodLines(d))) {
95  _b2 = _data;
96  _coefs = HiVector(4, 0.0);
97  _uncert = _coefs;
98  _cc = HiVector(2, 0.0);
99  _chisq = 0.0;
100  if ( !gotGoodLines(d) ) {
101  hist << "NotEnoughLines(GoodLines[" << goodLines(d)
102  << "],MinimumLines[" << _minLines << "]);";
103  }
104 
105  hist << "SkipFit(TRUE: Not using LMFit)";
106  _history.add(hist.str().c_str());
107  }
108  else {
109  hist << "Fit(";
110  _b2 = HiVector(goodLines(_data));
111  if ( success(curvefit()) ) {
112  _coefs = coefs();
113  _uncert = uncert();
114  hist << "Solved,#Iters[" << nIterations() << "],ChiSq[" << Chisq()
115  << "],DoF[" << DoF() << "])";
116  _history.add(hist.str().c_str());
117  _history.add("a0("+ToString(_coefs[0])+"+-"+ToString(_uncert[0])+")");
118  _history.add("a1("+ToString(_coefs[1])+"+-"+ToString(_uncert[1])+")");
119  _history.add("a2("+ToString(_coefs[2])+"+-"+ToString(_uncert[2])+")");
120  _history.add("a3("+ToString(_coefs[3])+"+-"+ToString(_uncert[3])+")");
121  }
122  else {
123  // Punt, fit a straight line to the data
124  _cc = poly_fit(d);
125  HiVector a(4);
126  a[0] = _cc[0];
127  a[1] = _cc[1];
128  a[2] = 0.0;
129  a[3] = 0.0;
130  _coefs = a;
131 
132  hist << "Failed::Reason("<< statusstr() << "),#Iters["
133  << nIterations() << "])";
134  _history.add(hist.str().c_str());
135  _history.add("a0("+ToString(_coefs[0])+")");
136  _history.add("a1("+ToString(_coefs[1])+")");
137  _history.add("a2("+ToString(_coefs[2])+")");
138  _history.add("a3("+ToString(_coefs[3])+")");
139  if ( _useLinFit ) {
140  _history.add("OnFailureUse(LinearFit(Zf))");
141  }
142  else {
143  _skipFit = true;
144  _history.add("OnFailureUse(ZfBuffer)");
145  }
146  }
147  }
148  return (Yfit());
149  }
150 
162  NonLinearLSQ::NLVector ZeroBufferFit::guess() {
163  int n = _data.dim();
164  int nb = n - _badLines;
165 
166  HiVector b1 = _data.subarray(0, nb-1);
167  LowPassFilter gfilter(b1, _history, _sWidth, _sIters);
168 
169  int nb2 = nb/2;
170  _b2 = gfilter.ref();
171  HiVector cc = poly_fit(_b2.subarray(nb2,_b2.dim()-1), nb2-1);
172 
173  // Compute the 3rd term guess by getting the average of the residual
174  // at both ends of the data set.
175  Statistics s;
176 
177  // Get the head of the data set
178  int n0 = MIN(nb, 20);
179  for ( int k = 0 ; k < n0 ; k++ ) {
180  double d = _b2[k] - (cc[0] + cc[1] * _timet(k));
181  s.AddData(&d, 1);
182  }
183  double head = s.Average();
184 
185  // Get the tail of the data set
186  s.Reset();
187  n0 = (int) (0.9 * nb);
188  for ( int l = n0 ; l < nb ; l++ ) {
189  double d = _b2[l] - (cc[0] + cc[1] * _timet(l));
190  s.AddData(&d, 1);
191  }
192  double tail = s.Average();
193 
194  // Populate the guess with the results
195  NLVector g(4, 0.0);
196  g[0] = cc[0];
197  g[1] = cc[1];
198  g[2] = head-tail;
199  g[3] = -5.0/_timet(nb-1);
200  _guess = g;
201  _history.add("Guess["+ToString(_guess[0])+ ","+
202  ToString(_guess[1])+ ","+
203  ToString(_guess[2])+ ","+
204  ToString(_guess[3])+ "]");
205  return (g);
206  }
207 
219  int ZeroBufferFit::checkIteration(const int Iter, const NLVector &fitcoefs,
220  const NLVector &uncerts, double cplxconj,
221  int Istatus) {
222  _chisq = pow(cplxconj, 2.0);
223  return (Istatus);
224  }
225 
227  NonLinearLSQ::NLVector ZeroBufferFit::f_x(const NLVector &a) {
228  double a0 = a[0];
229  double a1 = a[1];
230  double a2 = a[2];
231  double a3 = a[3];
232 
233  int n = _b2.dim();
234  NLVector f(n);
235  for (int i = 0 ; i < n ; i++) {
236  double lt = _timet(i);
237  double et = a3 * lt;
238  double Yi = a0 + (a1 * lt) + a2 * exp(MIN(et,_maxLog));
239  f[i] = (Yi - _b2[i]);
240  }
241  return (f);
242  }
243 
245  NonLinearLSQ::NLMatrix ZeroBufferFit::df_x(const NLVector &a) {
246  // double a0 = a[0];
247  // double a1 = a[1];
248  double a2 = a[2];
249  double a3 = a[3];
250 
251  int n = _b2.dim();
252  NLMatrix J(n, 4);
253  for (int i = 0; i < n; i++) {
254  double lt = _timet(i);
255  double et = a3 * lt;
256  double p0 = exp (MIN(et,_maxLog));
257  J[i][0] = 1.0;
258  J[i][1] = lt;
259  J[i][2] = p0;
260  J[i][3] = a2 * lt * p0;
261  }
262  return (J);
263  }
264 
266  HiVector ZeroBufferFit::Yfit() const {
267  if ( _skipFit || (!gotGoodLines(_data))) {
268  return (_data.copy());
269  }
270  else {
271  HiVector dcorr(_data.dim());
272  HiVector a = _coefs;
273  for ( int i = 0 ; i < dcorr.dim() ; i++ ) {
274  double lt = _timet(i);
275  dcorr[i] = a[0] + (a[1] * lt) + a[2] * exp(a[3] * lt);
276  }
277  return (dcorr);
278  }
279  }
280 
282  HiVector ZeroBufferFit::Normalize(const HiVector &v) {
283 
284  HiVector vNorm(v.dim());
285  double v0 = v[0];
286  for ( int i = 0 ; i < v.dim() ; i++ ) {
287  vNorm[i] = v[i] - v0;
288  }
289  _history.add("Normalize[" + ToString(v0) + "]");
290  return (vNorm);
291  }
292 
305  HiVector ZeroBufferFit::poly_fit(const HiVector &d, const double line0) const {
306  // Needs a linear fit to latter half of data
307  MultivariateStatistics fit;
308  int n = d.dim();
309  for ( int i = 0 ; i < n ; i++ ) {
310  double t = _timet(line0+i);
311  fit.AddData(&t, &d[i], 1);
312  }
313  NLVector cc(2);
314  fit.LinearRegression(cc[0], cc[1]);
315  return (cc);
316  }
317 
319  void ZeroBufferFit::printOn(std::ostream &o) const {
320  o << "# History = " << _history << endl;
321  // Write out the header
322  o << setw(_fmtWidth) << "Line"
323  << setw(_fmtWidth+1) << "Time"
324  << setw(_fmtWidth+1) << "Data"
325  << setw(_fmtWidth+1) << "Fit\n";
326 
327  HiVector fit = Yfit();
328  for (int i = 0 ; i < _data.dim() ; i++) {
329  o << formatDbl(i) << " "
330  << formatDbl(_timet(i)) << " "
331  << formatDbl(_data[i]) << " "
332  << formatDbl(fit[i]) << endl;
333  }
334  return;
335  }
336 
337 } // namespace Isis
Isis::Module::_history
HiHistory _history
Hierarchial component history.
Definition: Module.h:152
Isis::IsEqual
bool IsEqual(const QString &v1, const QString &v2="TRUE")
Shortened string equality test.
Definition: HiCalUtil.h:258
Isis::Module::_fmtWidth
int _fmtWidth
Default field with of double.
Definition: Module.h:153
Isis::Statistics
This class is used to accumulate statistics on double arrays.
Definition: Statistics.h:94
Isis::Statistics::AddData
void AddData(const double *data, const unsigned int count)
Add an array of doubles to the accumulators and counters.
Definition: Statistics.cpp:141
Isis::DbProfile::Name
QString Name() const
Returns the name of this property.
Definition: DbProfile.h:104
Isis::ZeroBufferFit::Solve
HiVector Solve(const HiVector &d)
Compute non-linear fit to (typically) ZeroBufferSmooth module.
Definition: ZeroBufferFit.cpp:91
Isis::ToString
QString ToString(const T &value)
Helper function to convert values to strings.
Definition: HiCalUtil.h:246
Isis::DbProfile::exists
bool exists(const QString &key) const
Checks for the existance of a keyword.
Definition: DbProfile.h:115
Isis::Module
Module manages HiRISE calibration vectors from various sources.
Definition: Module.h:39
Isis::Module::ref
const HiVector & ref() const
Return data via a const reference.
Definition: Module.h:111
Isis::ConfKey
T ConfKey(const DbProfile &conf, const QString &keyname, const T &defval, int index=0)
Find a keyword in a profile using default for non-existant keywords.
Definition: HiCalUtil.h:205
Isis::ZeroBufferFit::goodLines
int goodLines(const HiVector &d) const
Returns the number of good lines in the image.
Definition: ZeroBufferFit.h:144
Isis::Statistics::Reset
void Reset()
Reset all accumulators and counters to zero.
Definition: Statistics.cpp:113
Isis::NonLinearLSQ::uncert
NLVector uncert() const
Return uncertainties from last fit processing.
Definition: NonLinearLSQ.h:105
Isis::ZeroBufferFit::gotGoodLines
bool gotGoodLines(const HiVector &d) const
Determines if the vector contains any valid lines.
Definition: ZeroBufferFit.h:146
Isis::NonLinearLSQ::setMaxIters
void setMaxIters(int m)
Sets the maximum number of iterations.
Definition: NonLinearLSQ.h:62
Isis::LowPassFilter
Compute a low pass filter from a Module class content.
Definition: LowPassFilter.h:30
Isis::MultivariateStatistics::AddData
void AddData(const double *x, const double *y, const unsigned int count)
Add two arrays of doubles to the accumulators and counters.
Definition: MultivariateStatistics.cpp:66
Isis::NonLinearLSQ::coefs
NLVector coefs() const
Return coefficients from last fit processing.
Definition: NonLinearLSQ.h:103
Isis::ZeroBufferFit::DoF
int DoF() const
Returns the Degrees of Freedom.
Definition: ZeroBufferFit.h:117
Isis::Module::formatDbl
QString formatDbl(const double &value) const
Properly format values that could be special pixels.
Definition: Module.h:169
Isis::ZeroBufferFit::Normalize
HiVector Normalize(const HiVector &v)
Compute normalized solution vector from result.
Definition: ZeroBufferFit.cpp:282
Isis::NonLinearLSQ::success
bool success() const
Determine success from last fit processing.
Definition: NonLinearLSQ.h:83
Isis::NonLinearLSQ
NonLinearLSQ Computes a fit using a Levenberg-Marquardt algorithm.
Definition: NonLinearLSQ.h:42
Isis::MultivariateStatistics
Container of multivariate statistics.
Definition: MultivariateStatistics.h:54
Isis::toInt
int toInt(const QString &string)
Global function to convert from a string to an integer.
Definition: IString.cpp:93
Isis::NonLinearLSQ::maxIters
int maxIters() const
Maximum number iterations for valid solution.
Definition: NonLinearLSQ.h:70
Isis::ToInteger
int ToInteger(const T &value)
Helper function to convert values to Integers.
Definition: HiCalUtil.h:222
Isis::ZeroBufferFit::printOn
virtual void printOn(std::ostream &o) const
Provides virtualized dump of data from this module.
Definition: ZeroBufferFit.cpp:319
Isis::ZeroBufferFit::df_x
NLMatrix df_x(const NLVector &a)
Computes the first derivative of the function at the current iteration.
Definition: ZeroBufferFit.cpp:245
Isis::IsTrueValue
bool IsTrueValue(const DbProfile &prof, const QString &key, const QString &value="TRUE")
Determines if the keyword value is the expected value.
Definition: HiCalUtil.h:276
Isis::MIN
T MIN(const T &A, const T &B)
Implement templatized MIN fumnction.
Definition: PhotometricFunction.h:34
Isis::DbProfile
A DbProfile is a container for access parameters to a database.
Definition: DbProfile.h:51
Isis::HiVector
TNT::Array1D< double > HiVector
1-D Buffer
Definition: HiCalTypes.h:27
Isis::NonLinearLSQ::statusstr
std::string statusstr() const
Return error message pertaining to last fit procesing.
Definition: NonLinearLSQ.h:87
Isis::Statistics::Average
double Average() const
Computes and returns the average.
Definition: Statistics.cpp:300
Isis::ToDouble
double ToDouble(const T &value)
Helper function to convert values to doubles.
Definition: HiCalUtil.h:234
Isis::MultivariateStatistics::LinearRegression
void LinearRegression(double &a, double &b) const
Fits a line.
Definition: MultivariateStatistics.cpp:222
Isis::ZeroBufferFit::Yfit
HiVector Yfit() const
Computes the solution vector using current coefficents.
Definition: ZeroBufferFit.cpp:266
Isis::ZeroBufferFit::poly_fit
HiVector poly_fit(const HiVector &d, const double line0=0.0) const
Compute a polyonomial fit using multivariate statistics.
Definition: ZeroBufferFit.cpp:305
Isis::toDouble
double toDouble(const QString &string)
Global function to convert from a string to a double.
Definition: IString.cpp:149
std
Namespace for the standard library.
Isis::ZeroBufferFit::checkIteration
int checkIteration(const int Iter, const NLVector &fitcoefs, const NLVector &uncerts, double cplxconj, int Istatus)
Computes the interation check for convergence.
Definition: ZeroBufferFit.cpp:219
Isis::ZeroBufferFit::f_x
NLVector f_x(const NLVector &a)
Computes the function value at the current iteration.
Definition: ZeroBufferFit.cpp:227
Isis::ZeroBufferFit::guess
NLVector guess()
Compute the initial guess of the fit.
Definition: ZeroBufferFit.cpp:162
Isis::NonLinearLSQ::nIterations
int nIterations() const
Return number of iterations from last fit processing.
Definition: NonLinearLSQ.h:107
Isis::ZeroBufferFit::Chisq
double Chisq() const
Returns the Chi-Square value of the fit solution.
Definition: ZeroBufferFit.h:115
Isis
This is free and unencumbered software released into the public domain.
Definition: Apollo.h:16