/***************************************************************** * Copyright (c) 2020 by Douglas Wilhelm Harder. * All rights reserved. * * To use, please contact dwharder@uwaterloo.ca *****************************************************************/ #include static void linear_calculate_coefficients( linear_estimator_t *p_e ) { static void quadratic_calculate_coefficients( quadratic_estimator_t *p_e ) { static void linear_calculate_estimators( linear_estimator_t *p_e ); static void quadratic_calculate_estimators( quadratic_estimator_t *p_e ); // The constructor // // We assume that the initial state of the system is some constant // that is passed to the constructor, a parameter that has a default // value of 0.0. // // Assuming that all previous values were the passed value, then // // S = y N // y // N (1 - N) // SP = y --------- // xy 2 // // and the least-squares linear polynomial passing through these poitns // is p(k) = 0k + y void linear_estimator_init( linear_estimator_t *p_e, double y ) { p_e->s_y_ = y*CA_UWATERLOO_DWH_N; p_e->sp_xy_ = y*CA_UWATERLOO_DWH_N*(1.0 - CA_UWATERLOO_DWH_N)/2.0}, p_e->index_ = 0; p_e->a_ = 0.0; p_e->b_ = y; size_t k; for ( k = 0; k < CA_UWATERLOO_DWH_N; ++k ) { y_[k] = y; } linear_calculate_estimators( p_e ); } static void linear_calculate_coefficients( linear_estimator_t *p_e ) { // From these, derive the coefficients of the // least-squares best-fitting polynomial ax + b // - 6 FLOPs + 1 ASSIGN double const MA{6.0/((CA_UWATERLOO_DWH_N - 1.0)*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N + 1.0))}; p_e->a_ = ( CA_UWATERLOO_DWH_N - 1.0)*MA*p_e->s_y_ + 2.0*MA*p_e->sp_xy_; double const MB{2.0/(CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N + 1.0))}; p_e->b_ = (2.0*CA_UWATERLOO_DWH_N - 1.0)*MB*p_e->s_y_ + 3.0*MB*p_e->sp_xy_; } double linear_estimator_add( linear_estimator_t *p_e, double y, double epsilon ) { p_e->sp_xy_ += CA_UWATERLOO_DWH_N*p_e->y_[p_e->index_] - p_e->s_y_; p_e->s_y_ += -p_e->y_[p_e->index_] + y; p_e->y_[p_e->index_] = y; if ( epsilon != 0.0 ) { y = (CA_UWATERLOO_DWH_N*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(2.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N + 1.0)*y + ( (6.0*(3.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*y - 6.0*(3.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*p_e->s_y_ + 36.0*(1.0 - CA_UWATERLOO_DWH_N)*p_e->sp_xy_)*epsilon + (3.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*((CA_UWATERLOO_DWH_N + 9.0)*CA_UWATERLOO_DWH_N - 4.0)*y - 12.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(2.0*CA_UWATERLOO_DWH_N - 1.0)*p_e->s_y_ - 6.0*CA_UWATERLOO_DWH_N*(7.0*CA_UWATERLOO_DWH_N - 5.0)*p_e->sp_xy_) )*epsilon)/( (2.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 1.0)*CA_UWATERLOO_DWH_N*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N + 1.0) + 12.0*(2.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 1.0)*epsilon*(CA_UWATERLOO_DWH_N + epsilon) ); // Update the sum of the ys // - no need to upsdate the sum of the products xys p_e->s_y_ += -p_e->s_[p_e->sndex_] + y; p_e->s_[index_] = y; } // If after incrementation, index_ == N, reset N to 0 if ( ++p_e->sndex_ == CA_UWATERLOO_DWH_N ) { p_e->sndex_ = 0; } linear_calculate_coefficients( p_e ); linear_calculate_estimators( p_e ); return y; } double linear_estimator_recalculate_sums( linear_estimator_t *p_e ) { double old_a{ p_e->a_ }; double old_b{ p_e->b_ }; p_e->s_y_ = 0.0; p_e->sp_xy_ = 0.0; size_t k; for ( k = 0; k < CA_UWATERLOO_DWH_N; ++k ) { p_e->s_y_ += p_e->y_[k]; p_e->sp_xy_ += (-1.0*k)*p_e->y_[(CA_UWATERLOO_DWH_N + p_e->index_ - k) % CA_UWATERLOO_DWH_N]; } linear_estimator_calculate_estimators( p_e ); return std::sqrt( (p_e->a_ - old_a)*(p_e->a_ - old_a) + (p_e->b_ - old_b)*(p_e->b_ - old_b) ); } static void linear_calculate_estimators( linear_estimator_t *p_e ) { #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CURRENT p_e-current = p_e->b_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_NEXT p_e-next = p_e->b_ + p_e->a_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_PER_PERIOD p_e->change_per_period = p_e->a_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_AVERAGE_LAST_PERIOD p_e->average_last_period = p_e->b_ - p_e->a_/2.0; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_AVERAGE_NEXT_PERIOD p_e->average_next_period = p_e->b_ + p_e->a_/2.0; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_ZERO p_e->zero = -p_e->b_/p_e->a_; #endif } void quadratic_estimator_init( quadratic_estimator_t *p_e, double y ) { p_e->s_y_{y*CA_UWATERLOO_DWH_N}; p_e->sp_xy_{y*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)/2.0}; p_e->sp_xxy_{y*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(2.0*CA_UWATERLOO_DWH_N - 1.0)/6.0}; p_e->index_{0}; p_e->a_{0.0}; p_e->b_{0.0}; p_e->c_{y}; size_t k; for ( k = 0; k < CA_UWATERLOO_DWH_N; ++k ) { y_[k] = y; } quadratic_calculate_estimators( p_e ); } void quadratic_calculate_coefficients() { // 15 FLOPs and 3 ASSIGNs double const ma{30.0/((CA_UWATERLOO_DWH_N - 2.0)*(CA_UWATERLOO_DWH_N - 1.0)*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N + 1.0)*(CA_UWATERLOO_DWH_N + 2.0))}; p_e->a_ = (CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*ma*p_e->s_y_ + 6.0*(CA_UWATERLOO_DWH_N - 1.0)*ma*p_e->sp_xy_ + 6.0*ma*p_e->sp_xxy_; double const mb{6.0/((CA_UWATERLOO_DWH_N - 2.0)*(CA_UWATERLOO_DWH_N - 1.0)*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N + 1.0)*(CA_UWATERLOO_DWH_N + 2.0))}; p_e->b_ = (((6.0*CA_UWATERLOO_DWH_N - 21.0)*CA_UWATERLOO_DWH_N + 21.0)*CA_UWATERLOO_DWH_N - 6.0)*mb*p_e->s_y_ + 2.0*((16.0*CA_UWATERLOO_DWH_N - 30.0)*CA_UWATERLOO_DWH_N + 11.0)*mb*p_e->sp_xy_ + 30.0*(CA_UWATERLOO_DWH_N - 1.0)*mb*p_e->sp_xxy_; double const mc{3.0/(CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N + 1.0)*(CA_UWATERLOO_DWH_N + 2.0))}; p_e->c_ = ((3.0*CA_UWATERLOO_DWH_N - 3.0)*CA_UWATERLOO_DWH_N + 2.0)*mc*p_e->s_y_ + 6.0*(2.0*CA_UWATERLOO_DWH_N - 1.0)*mc*p_e->sp_xy_ + 10.0*mc*p_e->sp_xxy_; } double quadratic_estimator_add( quadratic_estimator_t *p_e, double y, double epsilon ) { p_e->sp_xxy_ += -CA_UWATERLOO_DWH_N*CA_UWATERLOO_DWH_N*p_e->y_[p_e->index_] - 2.0*p_e->sp_xy_ + p_e->s_y_; p_e->sp_xy_ += CA_UWATERLOO_DWH_N*p_e->y_[p_e->index_] - p_e->s_y_; p_e->s_y_ += -p_e->y_[p_e->index_] + y; p_e->y_[p_e->index_] = y; if ( epsilon != 0.0 ) { y = ( ( ( ( ( 60.0*((((5.0*CA_UWATERLOO_DWH_N - 46.0)*CA_UWATERLOO_DWH_N + 49.0)*CA_UWATERLOO_DWH_N - 32.0)*CA_UWATERLOO_DWH_N + 12.0)*y + 60.0*((((-5.0*CA_UWATERLOO_DWH_N + 46.0)*CA_UWATERLOO_DWH_N - 49.0)*CA_UWATERLOO_DWH_N + 32.0)*CA_UWATERLOO_DWH_N - 12.0)*p_e->s_y_ + 360.0*(((24.0 - 5.0*CA_UWATERLOO_DWH_N)*CA_UWATERLOO_DWH_N - 19.0)*CA_UWATERLOO_DWH_N + 6.0)*p_e->sp_xy_ - 1800.0*(CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*p_e->sp_xxy_ )*epsilon + 120.0*( 6.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(((CA_UWATERLOO_DWH_N - 6.0)*CA_UWATERLOO_DWH_N + 3.0)*CA_UWATERLOO_DWH_N - 2.0)*y - 6.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(((CA_UWATERLOO_DWH_N - 6.0)*CA_UWATERLOO_DWH_N + 3.0)*CA_UWATERLOO_DWH_N - 2.0)*p_e->s_y_ - CA_UWATERLOO_DWH_N*(2.0*CA_UWATERLOO_DWH_N - 1.0)*((17.0*CA_UWATERLOO_DWH_N - 57.0)*CA_UWATERLOO_DWH_N + 34.0)*p_e->sp_xy_ - 3.0*CA_UWATERLOO_DWH_N*((11.0*CA_UWATERLOO_DWH_N - 27.0)*CA_UWATERLOO_DWH_N + 22.0)*p_e->sp_xxy_ ) )*epsilon + ( 2.0*(((((((5.0*CA_UWATERLOO_DWH_N + 311.0)*CA_UWATERLOO_DWH_N - 1591.0)*CA_UWATERLOO_DWH_N + 2015.0)*CA_UWATERLOO_DWH_N - 802.0)*CA_UWATERLOO_DWH_N - 58.0)*CA_UWATERLOO_DWH_N + 264.0)*CA_UWATERLOO_DWH_N - 72.0)*y + 12.0*(((((3.0*(87.0 - 17.0*CA_UWATERLOO_DWH_N)*CA_UWATERLOO_DWH_N - 340.0)*CA_UWATERLOO_DWH_N + 137.0)*CA_UWATERLOO_DWH_N + 13.0)*CA_UWATERLOO_DWH_N - 44.0)*CA_UWATERLOO_DWH_N + 12.0)*p_e->s_y_ - 36.0*(CA_UWATERLOO_DWH_N - 1.0)*((((89.0*CA_UWATERLOO_DWH_N - 183.0)*CA_UWATERLOO_DWH_N + 52.0)*CA_UWATERLOO_DWH_N + 42.0)*CA_UWATERLOO_DWH_N - 12.0)*p_e->sp_xy_ + 120.0*(((25.0*(2.0 - CA_UWATERLOO_DWH_N)*CA_UWATERLOO_DWH_N - 32.0)*CA_UWATERLOO_DWH_N - 11.0)*CA_UWATERLOO_DWH_N + 6.0)*p_e->sp_xxy_ ) )*epsilon + ( 6.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*(2.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N + 1.0)*(((CA_UWATERLOO_DWH_N + 21.0)*CA_UWATERLOO_DWH_N - 16.0)*CA_UWATERLOO_DWH_N + 12.0)*y - 36.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*(2.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N + 1.0)*(3.0*(CA_UWATERLOO_DWH_N - 1.0)*CA_UWATERLOO_DWH_N + 2.0)*p_e->s_y_ - 24.0*CA_UWATERLOO_DWH_N*(2.0*CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N + 1.0)*(((21.0*CA_UWATERLOO_DWH_N - 60.0)*CA_UWATERLOO_DWH_N + 56.0)*CA_UWATERLOO_DWH_N - 20.0)*p_e->sp_xy_ - 180.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N + 1.0)*((5.0*CA_UWATERLOO_DWH_N - 8.0)*CA_UWATERLOO_DWH_N + 4.0)*p_e->sp_xxy_ ) )*epsilon + (CA_UWATERLOO_DWH_N*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*(2.0 + CA_UWATERLOO_DWH_N)*(3.0*(CA_UWATERLOO_DWH_N - 1.0)*CA_UWATERLOO_DWH_N + 2.0)*(CA_UWATERLOO_DWH_N + 1.0)*(CA_UWATERLOO_DWH_N + 1.0)*y) ) / ( ( ( ( 180.0*(CA_UWATERLOO_DWH_N - 2.0)*(CA_UWATERLOO_DWH_N - 1.0) *((3.0*CA_UWATERLOO_DWH_N - 3.0)*CA_UWATERLOO_DWH_N + 2.0)*epsilon + 360.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 2.0)*(CA_UWATERLOO_DWH_N - 1.0)*((3.0*CA_UWATERLOO_DWH_N - 3.0)*CA_UWATERLOO_DWH_N + 2.0) )*epsilon + 36.0*(CA_UWATERLOO_DWH_N - 2.0)*((7.0*CA_UWATERLOO_DWH_N + 1.0)*CA_UWATERLOO_DWH_N - 1.0) *(CA_UWATERLOO_DWH_N - 1.0)*((3.0*CA_UWATERLOO_DWH_N - 3.0)*CA_UWATERLOO_DWH_N + 2.0) )*epsilon + 36.0*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*(2.0*CA_UWATERLOO_DWH_N - 1.0) *(CA_UWATERLOO_DWH_N + 1.0)*((3.0*CA_UWATERLOO_DWH_N - 3.0)*CA_UWATERLOO_DWH_N + 2.0) )*epsilon + CA_UWATERLOO_DWH_N*CA_UWATERLOO_DWH_N*(CA_UWATERLOO_DWH_N - 1.0)*(CA_UWATERLOO_DWH_N - 2.0)*(2.0 + CA_UWATERLOO_DWH_N) *((3.0*CA_UWATERLOO_DWH_N - 3.0)*CA_UWATERLOO_DWH_N + 2.0)*(CA_UWATERLOO_DWH_N + 1.0)*(CA_UWATERLOO_DWH_N + 1.0) ); // Update the sum of the ys // - no need to upsdate the sum of the products xys p_e->s_y_ += -p_e->y_[p_e->index_] + y; p_e->y_[p_e->index_] = y; } // If after incrementation, index_ == N, reset N to 0 if ( ++p_e->index_ == CA_UWATERLOO_DWH_N ) { p_e->index_ = 0; } quadratic_calculate_coefficients( p_e ); quadratic_calculate_estimators( p_e ); return y; } double quadratic_estimator_recalculate_sums( quadratic_estimator_t *p_e ) { double old_a{ p_e->a_ }; double old_b{ p_e->b_ }; double old_c{ p_e->c_ }; p_e-> s_y_ = 0.0; p_e->sp_xy_ = 0.0; p_e->sp_xxy_ = 0.0; size_t k; for ( k = 0; k < CA_UWATERLOO_DWH_N; ++k ) { p_e-> s_y_ += p_e->y_[k]; p_e->sp_xy_ += (-1.0*k)*p_e->y_[(CA_UWATERLOO_DWH_N + p_e->index_ - k) % CA_UWATERLOO_DWH_N]; p_e->sp_xxy_ += k*k*p_e->y_[(CA_UWATERLOO_DWH_N + p_e->index_ - k) % CA_UWATERLOO_DWH_N]; } calculate_coefficients(); return std::sqrt( (p_e->a_ - old_a)*(p_e->a_ - old_a) + (p_e->b_ - old_b)*(p_e->b_ - old_b) + (p_e->c_ - old_c)*(p_e->c_ - old_c) ); } static void quadratic_calculate_estimators( linear_estimator_t *p_e ) { #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CURRENT p_e->current = p_e->c_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_NEXT p_e->next p_e->c_ + p_e->b_ + p_e->a_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_PER_PERIOD p_e->change_per_period = p_e->b_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_PER_PERIOD_PER_PERIOD p_e->change_per_period_per_period = 2.0*p_e->a_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_AVERAGE_LAST_PERIOD p_e->average_last_period = p_e->a_/3.0 - p_e->b_/2.0 + p_e->c_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_CHANGE_AVERAGE_NEXT_PERIOD p_e->average_next_period = p_e->a_/3.0 + p_e->b_/2.0 + p_e->c_; #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_ZERO p_e->zero1 = std::nan( "" ); p_e->zero2 = std::nan( "" ); #endif #ifndef CA_UWATERLOO_DWH_NO_LINEAR_ZERO p_e->extremum = -p_e->b_/(2.0*p_e->a_); #endif }