/***************************************************************** * Copyright (c) 2020 by Douglas Wilhelm Harder. * All rights reserved. * * To use, please contact dwharder@uwaterloo.ca *****************************************************************/ #ifndef CA_UWATERLOO_DWHARDER_LEAST_SQUARES_ESTIMATORS_H #define CA_UWATERLOO_DWHARDER_LEAST_SQUARES_ESTIMATORS_H /* You have the option of changing this default value, * either here or in your source code. */ #ifndef CA_UWATERLOO_DWHARDER_LSE_N #define CA_UWATERLOO_DWH_LSQ_N 10 #endif /**************************/ /* Structure declarations */ /**************************/ struct linear_estimator_t; struct quadratic_estimator_t; /***************************************************************** * These classes store data that can than be queried for * information * * Data points are assumed to be sampled periodically * and this class stores information that allows it to * calculate the least-squares best-fitting linear and * quadratic polynomials, respectively. * * Each time the add(...) member function is called, it * is assumed that this is the most recent sample, after * which any of the calls to the queries will use that * assumption: * * current() * Evaluate the least-squares polynomial * at the current sample time, used to * best estimate the value of the signal * at the current time. * * next() * Evaluate the least-squares polynomial * at the next time sample, used to best * estimate what the value of the signal * is likely to be one period in the * future. * * change_per_period() * Evaluate the instanteous change per * period of the least-squares * polynomial at the current time sample. * * - This can be used to estimate the * rate of change per unit time by * dividing this value by the period. * * average_last_period() * Evaluate the integral of the least- * squares polynomial over the last * time step, which will thus equal * the average value of the signal over * that period. * * - This can be used to estimate the * integral over the last period by * multiplying this average by the * period. * * average_next_period() * Evaluate the integral of the least- * squares polynomial over the next * time step, which will thus estimate * the average value of the signal over * the next period. * * - This can be used to estimate the * integral over the next period by * multiplying this average by the * period. * * zero() * Find the zero of the least-squares * polynomial. * * extremum() * Find the estimator of the extrema of * the sampled data by finding the extremum * of the least-squares quadratic polynomial. * * Mathematical Model * ================== * * As long as the samples are made periodically, * there is no need to consider the period, and * thus we assume that the most recent sample is * at k = 0, and previous readings are made at * k = -1, -2, -3, ..., 1 - N. Thus, if the current * reading is the nth reading, we are dealing * with the points: * * (1 - N, y[n - N + 1]), (2 - N, y[n - N + 2]), * ..., (-2, y[n - 2]), (-1, y[n - 1]), (0, y[n]) * * / 0 0 \ / 0 \ * | ----- ----- | | ----- | * | \ 2 \ | | \ | * | ) k ) k | | ) k y[n + k] | * | / / | | / | * | ----- ----- | / \ | ----- | * | k = 1 - N k = 1 - N || a | | k = 1 - N | * | || | = | | * | 0 0 || b | | 0 | * | ----- ----- | \ / | ----- | * | \ \ | | \ | * | ) k ) 1 | | ) y[n + k] | * | / / | | / | * | ----- ----- | | ----- | * \k = 1 - N k = 1 - N / \k = 1 - N / * * The left-hand matrix simplifies to * * / N (2 N - 1) (N - 1) N (1 - N) \ * | ------------------- --------- | * | 6 2 | * | | * | N (1 - N) | * | --------- N | * \ 2 / * * Thus, if we represent the right-hand vector as * * / SP \ * | xy | * | | * | S | * \ y / * * we may solve this system of linear equations to get * (N - 1) S + 2 SP * y xy * a = 6 ------------------- * (N - 1) N (N + 1) * * (2 N - 1) S + 3 SP * y xy * b = 2 --------------------- * N (N + 1) * * Thus, all that is necessary is to calculate S and SP . * y xy * * The problem here, is the calculation of the two sums are * an O(n) operation. Thus sum of the ys can, however, be * easily calculated as given the sum of the most recent * N values from -N to -1, we can calulate the sum from * 1 - N to 0 by subtracting y[n - N] and adding y[n] to * the previously calculated sum: * * Previous value: * S = y[n-N] + y[n-N+1] + ... + y[n-2] + y[n-1] * y * * Desired updated value: * S = y[n-N+1] + ... + y[n-2] + y[n-1] + y[n] * y * * Updating operation: * S <- y[n] + S - y[n-N] * y y * * To calculate the sum of the products, this becomes slightly * more complex, as we have * * Previous value: * SP = - (N - 1)y[n-N] - (N - 2)y[n-N+1] - (N - 3)y[n-N+2] * ... - 2y[n-3] - y[n-2] + 0 * xy * * Desired updated value: * SP = - (N - 1)y[n-N+1] - (N - 2)y[n-N+2] - * ... - 3y[n-3] - 2y[n-2] - y[n-1] + 0 * xy * * This can, hoewver, be calculated by the following three steps: * * Partially update the sum of the ys: * S <- S - y[n-N] * y y * * Remove the last term from the sum of the products and subtract * the partially updated sum of ys: * SP <- SP + (N - 1)y[n-N] - S * xy xy y * * Conclude by fully updating the sum of the ys: * S <- S - y[n-N] * y y * * Thus, S , SP , a and b can all be calculated in O(1) time with each period. * y xy * * Similar formulas can be used to caclculate S , SP , SP , a, b, and c * y xy xxy * when finding the interpolating quadratic polynomial, all of which can be * done in O(1) time. * * This class therefore allows the user to query properties of these * least-squares polynomials also in O(1) time. * * JITTER * ====== * Suppose that the current sample is not taken exactly at the current * time, but rather at the current time plus epsilon times the period. * Such an error in the sampling is referred to as jitter. If we were * to keep all values as is, at the times they were taken, then it * could not be possible to use these simplifed formulas as the * solutions to the least-squares polynomials would no longer have * a constant matrix. Thus, we must adjust for the jitter to have * the best estimator of what the value of the sample would be if * the sample was taken at the correct time. * * We will have one constrant, and that will allow us to find the * formula for this new value: specifically, we must find a y value * y* such that the estimator of the current value with the error * in the sample equals the estimator using y* at the correct time. * * The add(...) member function allows you to specify the jitter, * in which case, the return value will be the estimator of what the * sample would have been at the correct time to get the same estimator * at the current time. * * WEAKNESS * ======== * * Each time S , SP or SP is updated, it involves a number of * y xy xxy * FLOPs, each of which potentially introduces numerical error into * the next value. If all the y-values are within an order of * magnitude of each other, then after a thousand steps, the * actual value of these sums may drift from the calculated values. * * Thus, it is useful every thousand periods or so to recalculate * these sums by recalculating these sums, which will be an O(N) * operation. The recalculate_sums() member function does this * and returns the 2-norm error between the previous calculated * coefficients and the newly calculated coefficients of the * least-squares polynomial. *****************************************************************/ /*************************/ /* Structure definitions */ /*************************/ struct linear_estimator_t { #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CURRENT double current; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_NEXT double next; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_PER_PERIOD double change_per_period; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_AVERAGE_LAST_PERIOD double average_last_period; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_AVERAGE_NEXT_PERIOD double average_next_period; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_ZERO double zero; #endif double y_[CA_UWATERLOO_DWH_N]; double s_y_; double sp_xy_; std::size_t index_; /* The least-squares polynomial is ax + b */ double a_; double b_; }; class quadratic_estimator_t { // Estimators #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_CURRENT double current; #endif #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_NEXT double next; #endif #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_CHANGE_PER_PERIOD double change_per_period; #endif #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_CHANGE_PER_PERIOD_PER_PERIOD double change_per_period_per_period; #endif #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_AVERGE_LAST_PERIOD double average_last_period; #endif #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_AVERGE_NEXT_PERIOD double average_next_period; #endif #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_ZERO double zero1; double zero2; #endif #ifndef CA_UWATERLOO_DWH_NO_QUADRATIC_EXTREMUM double extremum; #endif double y_[CA_UWATERLOO_DWH_N]; double s_y_; double sp_xy_; double sp_xxy_; std::size_t index_; // 2 // The least-squares polynomial is ax + bx + c double a_; double b_; double c_; }; /*************************/ /* Function declarations */ /*************************/ void linear_estimator_init( struct least_squares_estimator_t *p_e, double y ); double linear_estimator_add( struct least_squares_estimator_t *p_e, double y, double epsilon ); double linear_estimator_recalculate_sums( struct least_squares_estimator_t *p_e ); void quadratic_estimator_init( struct least_squares_estimator_t *p_e, double y ); double quadratic_estimator_add( struct least_squares_estimator_t *p_e, double y, double epsilon ); double quadratic_estimator_recalculate_sums( struct least_squares_estimator_t *p_e ); #endif