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
387 double sum = va * va * Daa + 2 * (va * vb * Dab + va * vc * Dac + vb * vc * Dbc) + vc * vc * Dcc;
447// Calculates the second directional derivative vector for the parabola equation f(x) = a + b*(x-c)^2
473double gaufxy(double x, double y, double a, double x0, double y0, double A, double B, double C, double b)
475 return b + a * exp(-(A * (pow(x - x0, 2.0)) + 2.0 * B * (x - x0) * (y - y0) + C * (pow(y - y0, 2.0))));
549// The matrix is symmetric, as ddf/dParm1dParm2 = ddf/dParm2dParm1 so we only need to differentiate
646 double sum = 2 * va * ((vx0 * Dax0) + (vy0 * Day0) + (vA * DaA) + (vB * DaB) + (vC * DaC)) + // a diffs
647 vx0 * ((vx0 * Dx0x0) + 2 * ((vy0 * Dx0y0) + (vA * Dx0A) + (vB * Dx0B) + (vC * Dx0C))) + // x0 diffs
714// The matrix is symmetric, as ddf/dParm1dParm2 = ddf/dParm2dParm1 so we only need to differentiate
769void CurveFitting::fitCurve(const FittingGoal goal, const QVector<int> &x_, const QVector<double> &y_,
773 if ((x_.size() != y_.size()) || (x_.size() != weight_.size()) || (x_.size() != outliers_.size()))
774 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve inconsistent parameters. x=%1, y=%2, weight=%3, outliers=%4")
804 // Something went wrong, log an error and reset state so solver starts from scratch if called again
805 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve called with curveFit=%1").arg(curveFit);
818 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D called with useWeights = true. Not currently supported");
833 // Something went wrong, log an error and reset state so solver starts from scratch if called again
834 qCDebug(KSTARS_EKOS_FOCUS) << QString("CurveFitting::fitCurve3D called with curveFit=%1").arg(curveFit);
863 y = hypfx(x, m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX], m_coefficients[D_IDX]);
867 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error: CurveFitting::f called with m_CurveType = %1 m_coefficients.size = %2")
877 z = gaufxy(x, y, m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX], m_coefficients[D_IDX],
882 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error: CurveFitting::f3D called with m_CurveType = %1 m_coefficients.size = %2")
888QVector<double> CurveFitting::polynomial_fit(const double *const data_x, const double *const data_y, const int n,
936QVector<double> CurveFitting::hyperbola_fit(FittingGoal goal, const QVector<double> data_x, const QVector<double> data_y,
937 const QVector<double> data_weights, const QVector<bool> outliers, const bool useWeights, const OptimisationDirection optDir)
955 gsl_multifit_nlinear_workspace *w = gsl_multifit_nlinear_alloc(gsl_multifit_nlinear_trust, ¶ms, dataPoints.dps.size(),
972 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
985 qCDebug(KSTARS_EKOS_FOCUS) << QString("iter %1: A=%2, B=%3, C=%4, D=%5 rcond(J)=%6, avratio=%7, |f(x)|=%8")
1003 // fails on its first step where we will always retry after adjusting parameters. It helps with
1023 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=hyperbola, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1047 QString("LM solver (Hyperbola): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7], retry=%8")
1048 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1061 QString("LM Solver (Hyperbola): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6, D=%7")
1062 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(getLMReasonCode(info)).arg(vc[A_IDX]).arg(vc[B_IDX])
1094void CurveFitting::hypSetupParams(FittingGoal goal, gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol,
1106 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1107 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1108 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1112 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1116 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1152// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1153// If we found a solution before and we're just adding more datapoints use the last solution as the guess.
1154// If we don't have a solution use the datapoints to approximate. Work out the min and max points and use these
1156void CurveFitting::hypMakeGuess(const int attempt, const DataPointT &dataPoints, gsl_vector * guess)
1158 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1162 if (!m_FirstSolverRun && (m_LastCurveType == FOCUS_HYPERBOLA) && (m_LastCoefficients.size() == NUM_HYPERBOLA_PARAMS))
1242 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM solver (hyp) initial params: perturbation=%1, A=%2, B=%3, C=%4, D=%5")
1252QVector<double> CurveFitting::parabola_fit(FittingGoal goal, const QVector<double> data_x, const QVector<double> data_y,
1253 const QVector<double> data_weights, const QVector<bool> outliers, bool useWeights, const OptimisationDirection optDir)
1271 gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc (gsl_multifit_nlinear_trust, ¶ms, dataPoints.dps.size(),
1288 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
1301 qCDebug(KSTARS_EKOS_FOCUS) << QString("iter %1: A=%2, B=%3, C=%4, rcond(J)=%5, avratio=%6 |f(x)|=%7")
1317 // fails on its first step where we will always retry after adjusting parameters. It helps with
1337 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=parabola, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1361 QString("LM solver (Parabola): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7], retry=%8")
1362 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1374 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Parabola): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6")
1375 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(getLMReasonCode(info)).arg(vc[A_IDX]).arg(vc[B_IDX])
1393void CurveFitting::parSetupParams(FittingGoal goal, gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol,
1405 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1406 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1407 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1411 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1415 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1451// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1452// If we found a solution before and we're just adding more datapoints use the last solution as the guess.
1453// If we don't have a solution use the datapoints to approximate. Work out the min and max points and use these
1455void CurveFitting::parMakeGuess(const int attempt, const DataPointT &dataPoints, gsl_vector * guess)
1457 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1461 if (!m_FirstSolverRun && (m_LastCurveType == FOCUS_PARABOLA) && (m_LastCoefficients.size() == NUM_PARABOLA_PARAMS))
1509 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM solver (par) initial params: perturbation=%1, A=%2, B=%3, C=%4")
1527 gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc (gsl_multifit_nlinear_trust, ¶ms, data.dps.size(),
1550 // Setup for multiple solve attempts. We won't worry too much if the solver fails as there should be
1551 // plenty of stars, but if the solver fails on its first step then adjust parameters and retry as this
1573 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
1574 /*auto callback = [](const size_t iter, void* _params, const gsl_multifit_nlinear_workspace * w)
1605 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=gaussian, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1616 QString("LM solver (Gaussian): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7]")
1617 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1629 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Gaussian): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6, "
1630 "D=%7, E=%8, F=%9, G=%10").arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(getLMReasonCode(info))
1648// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1649// Since we have already run some HFR calcs on the star, use these values to calculate the guess
1650void CurveFitting::gauMakeGuess(const int attempt, const StarParams &starParams, gsl_vector * guess)
1652 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1676 QString("LM Solver (Gaussian): Guess perturbation=%1, A=%2, B=%3, C=%4, D=%5, E=%6, F=%7, G=%8")
1689void CurveFitting::gauSetupParams(gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol, double *gtol,
1700 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1705 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1706 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1707 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1711 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1732 gsl_multifit_nlinear_workspace* w = gsl_multifit_nlinear_alloc (gsl_multifit_nlinear_trust, ¶ms, data.dps.size(),
1776 // To activate, uncomment the callback lambda and change the call to gsl_multifit_nlinear_driver
1777 /*auto callback = [](const size_t iter, void* _params, const gsl_multifit_nlinear_workspace * w)
1804 qCDebug(KSTARS_EKOS_FOCUS) << QString("Starting LM solver, fit=plane, solver=%1, scale=%2, trs=%3, iters=%4, xtol=%5,"
1815 QString("LM solver (Plane): Failed after %1ms iters=%2 [attempt=%3] with status=%4 [%5] and info=%6 [%7]")
1816 .arg(timer.elapsed()).arg(gsl_multifit_nlinear_niter(w)).arg(attempt + 1).arg(status).arg(gsl_strerror(status))
1828 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Plane): Solution found after %1ms %2 iters (%3). A=%4, B=%5, C=%6")
1846// Initialise parameters before starting the solver. Its important to start with a guess as near to the solution as possible
1850 // If we are retrying then perturb the initial conditions. The hope is that by doing this the solver
1854 if (!m_FirstSolverRun && (m_LastCurveType == FOCUS_PLANE) && (m_LastCoefficients.size() == NUM_PLANE_PARAMS) && attempt < 3)
1868 qCDebug(KSTARS_EKOS_FOCUS) << QString("LM Solver (Plane): Guess perturbation=%1, A=%2, B=%3, C=%4")
1877void CurveFitting::plaSetupParams(gsl_multifit_nlinear_parameters *params, int *numIters, double *xtol, double *gtol,
1888 // avmax is the max allowed ratio of 2nd order term (acceleration, a) to the 1st order term (velocity, v)
1893 // - gsl_multifit_nlinear_scale_more uses. More strategy. Good for problems with parameters with widely different scales
1894 // - gsl_multifit_nlinear_scale_levenberg. Levensberg strategy. Can be good but not scale invariant.
1895 // - gsl_multifit_nlinear_scale_marquardt. Marquardt strategy. Considered inferior to More and Levensberg
1899 // - gsl_multifit_nlinear_solver_qr produces reliable results but needs more iterations than Cholesky
1910bool CurveFitting::findMinMax(double expected, double minPosition, double maxPosition, double *position, double *value,
1930 // 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
1935bool CurveFitting::minimumQuadratic(double expected, double minPosition, double maxPosition, double *position,
1984bool CurveFitting::minMaxHyperbola(double expected, double minPosition, double maxPosition, double *position,
1991 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::minimumHyperbola coefficients.size()=%1").arg(
2002 // Check 2: At the minimum position (=c), the value of f(x) (which is the HFR) given by b+d is > 0
2006 ((optDir == OPTIMISATION_MINIMISE && b > 0.0) || (optDir == OPTIMISATION_MAXIMISE && b < 0.0)))
2016bool CurveFitting::minMaxParabola(double expected, double minPosition, double maxPosition, double *position,
2023 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::minimumParabola coefficients.size()=%1").arg(
2034 // Check 2: At the minimum position (=c), the value of f(x) (which is the HFR) given by a is > 0
2038 ((optDir == OPTIMISATION_MINIMISE && b > 0.0) || (optDir == OPTIMISATION_MAXIMISE && b < 0.0)))
2061 // 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
2072 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::getGaussianParams coefficients.size()=%1").arg(
2088 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::getGaussianParams a=%1 b=%2 x0=%3 y0=%4").arg(a).arg(b)
2110 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::getGaussianParams FWHM=%1").arg(FWHM);
2126// R2 (R squared) is the coefficient of determination gives a measure of how well the curve fits the datapoints.
2127// It lies between 0 and 1. 1 means that all datapoints will lie on the curve which therefore exactly describes the
2140 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 called for Quadratic");
2147 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 hyperbola coefficients.size()=").arg(
2157 curvePoints.push_back(hypfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX],
2170 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 parabola coefficients.size()=").arg(
2180 curvePoints.push_back(parfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX], m_coefficients[C_IDX]));
2192 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 gaussian coefficients.size()=").arg(
2200 curvePoints.push_back(gaufxy(m_dataPoints.dps[i].x, m_dataPoints.dps[i].y, m_coefficients[A_IDX],
2214 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 plane coefficients.size()=").arg(
2222 curvePoints.push_back(plafxy(m_dataPoints.dps[i].x, m_dataPoints.dps[i].y, m_coefficients[A_IDX],
2232 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateR2 curveFit=%1").arg(curveFit);
2241double CurveFitting::calcR2(const QVector<double> dataPoints, const QVector<double> curvePoints,
2265 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calcR2 tss=%1").arg(totalSumSquares);
2273void CurveFitting::calculateCurveDeltas(CurveFit curveFit, std::vector<std::pair<int, double>> &curveDeltas)
2281 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas called for Quadratic");
2287 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas hyperbola coefficients.size()=").arg(
2294 curveDeltas.push_back(std::make_pair(i, abs(m_y[i] - hypfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX],
2301 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas parabola coefficients.size()=").arg(
2308 curveDeltas.push_back(std::make_pair(i, abs(m_y[i] - parfx(m_x[i], m_coefficients[A_IDX], m_coefficients[B_IDX],
2314 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas called for Gaussian");
2318 qCDebug(KSTARS_EKOS_FOCUS) << QString("Error in CurveFitting::calculateCurveDeltas curveFit=%1").arg(curveFit);
KGuiItem ok()
NETWORKMANAGERQT_EXPORT NetworkManager::Status status()
void start()
qsizetype size() const const
QString & append(QChar ch)
QString arg(Args &&... args) const const
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 Tue Mar 26 2024 11:19:02 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006
Documentation copyright © 1996-2024 The KDE developers.
Generated on Tue Mar 26 2024 11:19:02 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006
KDE's Doxygen guidelines are available online.