Kstars
curvefit.cpp
24// The functions here fit a number of different curves to the incoming data points using the Lehvensberg-Marquart
25// solver with geodesic acceleration as provided the Gnu Science Library (GSL). The following sources of information are useful:
26// www.wikipedia.org/wiki/Levenberg–Marquardt_algorithm - overview of the mathematics behind the algo
27// www.gnu.org/software/gsl/doc/html/nls.html - GSL manual describing Nonlinear Least-Squares fitting.
31// The Levensberg-Marquart algorithm is a non-linear least-squares solver and thus suitable for mnay
32// different equations. The basic is idea is to adjust the equation y = f(x,P) so that the computed
33// y values are as close as possible to the y values of the datapoints provided, so that the resultant
34// curve fits the data as best as it can. P is a set of parameters that are varied by the solver in order
35// to find the best fit. The solver measures how far away the curve is at each data point, squares the
38// The solver is supplied with an initial guess for the parameters, P. It calculates S, makes an adjustment
39// P and calculates a new S1. Provided S1 < S the solver has moved in the right direction and reduced S.
45// The solver is capable of solving either an unweighted or weighted set of datapoints. In essence, an
46// unweighted set of data gives equal weight to each datapoint when trying to fit a curve. An alternative is to
47// weight each datapoint with a measure that corresponds to how accurate the measurement of the datapoint
48// actually is. Given we are calculating HFRs of stars we can calculate the variance (standard deviation
49// squared) in the measurement of HFR and use 1 / SD^2 as weighting. What this means is that if we have a
50// datapoint where the SD in HFR measurements is small we would give this datapoint a relatively high
51// weighting when fitting the curve, and if there was a datapoint with a higher SD in HFR measurements
54// There are several optimisations to LM that speed up convergence for certain types of equations. This
55// code uses the LM solver with geodesic acceleration; a technique to accelerate convergence by using the
56// second derivative of the fitted function. An optimisation that has been implemented is that since, in normal
57// operation, the solver is run multiple times with the same or similar parameters, the initial guess for
60// The documents referenced provide the maths of the LM algorithm. What is required is the function to be
61// used f(x) and the derivative of f(x) with respect to the parameters of f(x). This partial derivative forms
62// a matrix called the Jacobian, J(x). LM then finds a minimum. Mathematically it will find a minimum local
63// to the initial guess, but not necessarily the global minimum but for the equations used this is not a
64// problem under normal operation as there will only be 1 minimum. If the data is poor, however, the
65// solver may find an "n" solution rather than a "u" especially at the start of the process where there
66// are a small number of points and the error in measurement is high. We'll ignore and continue and this
69// The original algorithm used in Ekos is a quadratic equation, which represents a parabolic curve. In order
70// not to disturb historic code the original solution will be called Quadratic. It uses a linear least-squares
115// 2. b must be > 0 for a V shaped curve (b < 0 for an inverted V and b = 0 for a horizontal line).
128// 2. b must be > 0 for a V shaped curve (b < 0 for an inverted V and b = 0 for a horizontal line).
164// The matrix is symmetric, as ddf/dParm1dParm2 = ddf/dParm2dParm1 so we only need to differentiate
214// The matrix is symmetric, as ddf/dParm1dParm2 = ddf/dParm2dParm1 so we only need to differentiate
229// There following are key parameters that drive the LM solver. Currently these settings are coded into
230// the program. It would be better to store them in a configuration file somewhere, but they aren't the
397 double sum = va * va * Daa + 2 * (va * vb * Dab + va * vc * Dac + vb * vc * Dbc) + vc * vc * Dcc;
466// Calculates the second directional derivative vector for the parabola equation f(x) = a + b*(x-c)^2
492double gaufxy(double x, double y, double a, double x0, double y0, double A, double B, double C, double b)
494 return b + a * exp(-(A * (pow(x - x0, 2.0)) + 2.0 * B * (x - x0) * (y - y0) + C * (pow(y - y0, 2.0))));
568// The matrix is symmetric, as ddf/dParm1dParm2 = ddf/dParm2dParm1 so we only need to differentiate
665 double sum = 2 * va * ((vx0 * Dax0) + (vy0 * Day0) + (vA * DaA) + (vB * DaB) + (vC * DaC)) + // a diffs
666 vx0 * ((vx0 * Dx0x0) + 2 * ((vy0 * Dx0y0) + (vA * Dx0A) + (vB * Dx0B) + (vC * Dx0C))) + // x0 diffs
733// The matrix is symmetric, as ddf/dParm1dParm2 = ddf/dParm2dParm1 so we only need to differentiate
788void CurveFitting::fitCurve(const FittingGoal goal, const QVector<int> &x_, const QVector<double> &y_,
792 if ((x_.size() != y_.size()) || (x_.size() != weight_.size()) || (x_.size() != outliers_.size()))
793 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve inconsistent parameters. x=%1, y=%2, weight=%3, outliers=%4")
823 // Something went wrong, log an error and reset state so solver starts from scratch if called again
824 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve called with curveFit=%1").arg(curveFit);
837 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D called with useWeights = true. Not currently supported");
852 // Something went wrong, log an error and reset state so solver starts from scratch if called again
853 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D called with curveFit=%1").arg(curveFit);
882 y = hypfx(x, m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX], m_coefficients[D_IDX]);
886 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error: CurveFitting::f called with m_CurveType = %1 m_coefficients.size = %2")
896 x = hypfy(y, m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX], m_coefficients[D_IDX]);
900 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error: CurveFitting::fy called with m_CurveType = %1 m_coefficients.size = %2")
910 z = gaufxy(x, y, m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX], m_coefficients[D_IDX],
915 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error: CurveFitting::f3D called with m_CurveType = %1 m_coefficients.size = %2")
921QVector<double> CurveFitting::polynomial_fit(const double *const data_x, const double *const data_y, const int n,
969QVector<double> CurveFitting::hyperbola_fit(FittingGoal goal, const QVector<double> data_x, const QVector<double> data_y,
970 const QVector<double> data_weights, const QVector<bool> outliers, const bool useWeights, const OptimisationDirection optDir)
988 gsl_multifit_nlinear_workspace *w = gsl_multifit_nlinear_alloc(gsl_multifit_nlinear_trust, ¶ms, dataPoints.dps.size(),
1005 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
1018 qCDebug(KSTARS_EKOS_FOCUS) << QString("iter %1: A=%2, B=%3, C=%4, D=%5 rcond(J)=%6, avratio=%7, |f(x)|=%8")
1036 // fails on its first step where we will always retry after adjusting parameters. It helps with
1056 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=hyperbola, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1080 QString("LM solver (Hyperbola): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7], retry=%8")
1081 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1094 QString("LM Solver (Hyperbola): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6, D=%7")
1095 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(getLMReasonCode(info)).arg(vc[A_IDX]).arg(vc[B_IDX])
1127void CurveFitting::hypSetupParams(FittingGoal goal, gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol,
1139 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1140 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1141 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1145 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1149 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1185// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1186// If we found a solution before and we're just adding more datapoints use the last solution as the guess.
1187// If we don't have a solution use the datapoints to approximate. Work out the min and max points and use these
1189void CurveFitting::hypMakeGuess(const int attempt, const DataPointT &dataPoints, gsl_vector * guess)
1191 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1195 if (!m_FirstSolverRun && (m_LastCurveType == FOCUS_HYPERBOLA) && (m_LastCoefficients.size() == NUM_HYPERBOLA_PARAMS))
1275 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM solver (hyp) initial params: perturbation=%1, A=%2, B=%3, C=%4, D=%5")
1285QVector<double> CurveFitting::parabola_fit(FittingGoal goal, const QVector<double> data_x, const QVector<double> data_y,
1286 const QVector<double> data_weights, const QVector<bool> outliers, bool useWeights, const OptimisationDirection optDir)
1304 gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc (gsl_multifit_nlinear_trust, ¶ms, dataPoints.dps.size(),
1321 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
1334 qCDebug(KSTARS_EKOS_FOCUS) << QString("iter %1: A=%2, B=%3, C=%4, rcond(J)=%5, avratio=%6 |f(x)|=%7")
1350 // fails on its first step where we will always retry after adjusting parameters. It helps with
1370 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=parabola, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1394 QString("LM solver (Parabola): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7], retry=%8")
1395 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1407 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Parabola): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6")
1408 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(getLMReasonCode(info)).arg(vc[A_IDX]).arg(vc[B_IDX])
1426void CurveFitting::parSetupParams(FittingGoal goal, gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol,
1438 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1439 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1440 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1444 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1448 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1484// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1485// If we found a solution before and we're just adding more datapoints use the last solution as the guess.
1486// If we don't have a solution use the datapoints to approximate. Work out the min and max points and use these
1488void CurveFitting::parMakeGuess(const int attempt, const DataPointT &dataPoints, gsl_vector * guess)
1490 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1494 if (!m_FirstSolverRun && (m_LastCurveType == FOCUS_PARABOLA) && (m_LastCoefficients.size() == NUM_PARABOLA_PARAMS))
1542 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM solver (par) initial params: perturbation=%1, A=%2, B=%3, C=%4")
1560 gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc (gsl_multifit_nlinear_trust, ¶ms, data.dps.size(),
1583 // Setup for multiple solve attempts. We won't worry too much if the solver fails as there should be
1584 // plenty of stars, but if the solver fails on its first step then adjust parameters and retry as this
1606 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
1607 /*auto callback = [](const size_t iter, void* _params, const gsl_multifit_nlinear_workspace * w)
1638 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=gaussian, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1649 QString("LM solver (Gaussian): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7]")
1650 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1662 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Gaussian): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6, "
1663 "D=%7, E=%8, F=%9, G=%10").arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(getLMReasonCode(info))
1681// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1682// Since we have already run some HFR calcs on the star, use these values to calculate the guess
1683void CurveFitting::gauMakeGuess(const int attempt, const StarParams &starParams, gsl_vector * guess)
1685 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1709 QString("LM Solver (Gaussian): Guess perturbation=%1, A=%2, B=%3, C=%4, D=%5, E=%6, F=%7, G=%8")
1722void CurveFitting::gauSetupParams(gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol, double *gtol,
1733 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1738 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1739 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1740 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1744 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1765 gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc (gsl_multifit_nlinear_trust, ¶ms, data.dps.size(),
1809 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
1810 /*auto callback = [](const size_t iter, void* _params, const gsl_multifit_nlinear_workspace * w)
1837 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=plane, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1848 QString("LM solver (Plane): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7]")
1849 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1861 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Plane): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6")
1879// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1883 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1887 if (!m_FirstSolverRun && (m_LastCurveType == FOCUS_PLANE) && (m_LastCoefficients.size() == NUM_PLANE_PARAMS) && attempt < 3)
1901 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Plane): Guess perturbation=%1, A=%2, B=%3, C=%4")
1910void CurveFitting::plaSetupParams(gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol, double *gtol,
1921 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1926 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1927 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1928 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1932 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1943bool CurveFitting::findMinMax(double expected, double minPosition, double maxPosition, double *position, double *value,
1963 // If we couldn't fit a curve then something's wrong so reset coefficients which will force the next run of the LM solver to start from scratch
1968bool CurveFitting::minimumQuadratic(double expected, double minPosition, double maxPosition, double *position,
2017bool CurveFitting::minMaxHyperbola(double expected, double minPosition, double maxPosition, double *position,
2024 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::minimumHyperbola coefficients.size()=%1").arg(
2035 // Check 2: At the minimum position (=c), the value of f(x) (which is the HFR) given by b+d is > 0
2039 ((optDir == OPTIMISATION_MINIMISE && b > 0.0) || (optDir == OPTIMISATION_MAXIMISE && b < 0.0)))
2049bool CurveFitting::minMaxParabola(double expected, double minPosition, double maxPosition, double *position,
2056 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::minimumParabola coefficients.size()=%1").arg(
2067 // Check 2: At the minimum position (=c), the value of f(x) (which is the HFR) given by a is > 0
2071 ((optDir == OPTIMISATION_MINIMISE && b > 0.0) || (optDir == OPTIMISATION_MAXIMISE && b < 0.0)))
2094 // If we couldn't fit a curve then something's wrong so reset coefficients which will force the next run of the LM solver to start from scratch
2105 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::getGaussianParams coefficients.size()=%1").arg(
2121 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::getGaussianParams a=%1 b=%2 x0=%3 y0=%4").arg(a).arg(b)
2143 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::getGaussianParams FWHM=%1").arg(FWHM);
2159// R2 (R squared) is the coefficient of determination gives a measure of how well the curve fits the datapoints.
2160// It lies between 0 and 1. 1 means that all datapoints will lie on the curve which therefore exactly describes the
2173 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 called for Quadratic");
2180 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 hyperbola coefficients.size()=").arg(
2190 curvePoints.push_back(hypfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX],
2203 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 parabola coefficients.size()=").arg(
2213 curvePoints.push_back(parfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX]));
2225 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 gaussian coefficients.size()=").arg(
2233 curvePoints.push_back(gaufxy(m_dataPoints.dps[i].x, m_dataPoints.dps[i].y, m_coefficients[A_IDX],
2247 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 plane coefficients.size()=").arg(
2255 curvePoints.push_back(plafxy(m_dataPoints.dps[i].x, m_dataPoints.dps[i].y, m_coefficients[A_IDX],
2265 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 curveFit=%1").arg(curveFit);
2274double CurveFitting::calcR2(const QVector<double> dataPoints, const QVector<double> curvePoints,
2298 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calcR2 tss=%1").arg(totalSumSquares);
2306void CurveFitting::calculateCurveDeltas(CurveFit curveFit, std::vector<std::pair<int, double>> &curveDeltas)
2314 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas called for Quadratic");
2320 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas hyperbola coefficients.size()=").arg(
2327 curveDeltas.push_back(std::make_pair(i, abs(m_y[i] - hypfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX],
2334 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas parabola coefficients.size()=").arg(
2341 curveDeltas.push_back(std::make_pair(i, abs(m_y[i] - parfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX],
2347 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas called for Gaussian");
2351 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas curveFit=%1").arg(curveFit);
KGuiItem ok()
NETWORKMANAGERQT_EXPORT NetworkManager::Status status()
void start()
void append(QList< T > &&value)
void clear()
void push_back(parameter_type value)
qsizetype size() const const
QString & append(QChar ch)
QString arg(Args &&... args) const const
void push_back(QChar ch)
QStringList split(QChar sep, Qt::SplitBehavior behavior, Qt::CaseSensitivity cs) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Jul 26 2024 11:59:51 by doxygen 1.11.0 written by Dimitri van Heesch, © 1997-2006
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri Jul 26 2024 11:59:51 by doxygen 1.11.0 written by Dimitri van Heesch, © 1997-2006
KDE's Doxygen guidelines are available online.