#include "lsq.h"

/*
 *
 * Least squares fit:
 * 
 *   y = b0 + b1x
 *
 *        n*sum(xi*yi) - (sum(xi)*sum(yi))
 *   b1 = --------------------------------
 *        n*sum(xi^2) - (sum(xi))^2
 * 
 *
 *   b0 = sum(yi)/n - b1*(sum(xi)/n)
 *
 */
 
 
void 
least_squares(double *x, double *y, int n, double *m, double *b, struct lsq *p) 
{
   int i;
 
   p->sum_xi = p->sum_yi = p->sum_xi_2 = p->sum_xi_yi = 0.0;
   p->sum_n = n;
 
   for ( i = 0; i < n; i++ ) {
     p->sum_xi += x[i];
     p->sum_yi += y[i];
     p->sum_xi_2 += x[i] * x[i];
     p->sum_xi_yi += x[i] * y[i];
     }
 
   *m = ( (double)p->sum_n * p->sum_xi_yi - p->sum_xi * p->sum_yi ) / 
          ( (double)p->sum_n * p->sum_xi_2 - p->sum_xi * p->sum_xi );
   *b = (p->sum_yi / (double)p->sum_n) - (*m) * (p->sum_xi / (double)p->sum_n); 
}
 
 
 /*
  * incrementally update existing values with a new data point 
  */
 
void 
least_squares_update(double x, double y, double *m, double *b, struct lsq *p) 
{
   ++(p->sum_n);
 
   p->sum_xi += x;
   p->sum_yi += y;
   p->sum_xi_2 += x * x;
   p->sum_xi_yi += x * y;

   *m = ( (double)p->sum_n * p->sum_xi_yi - p->sum_xi * p->sum_yi ) 
          / ( (double)p->sum_n * p->sum_xi_2 - p->sum_xi * p->sum_xi );
   *b = (p->sum_yi / (double)p->sum_n) - (*m) * (p->sum_xi / (double)p->sum_n);
}
 
 
 /* 
  * return the least squares error:
  *
  *             (y[i] - y_hat[i])^2
  *             -------------------
  *                     n
y  */
 
double 
least_squares_error(double *x, double *y, int n, double m, double b) {
  int i;
  double error, sum;
 
  sum = 0.0;
 
  for ( i = 0; i < n; i++ ) {
    error = y[i] - (m * x[i] + b);
    sum += error * error;
    }
 
  return ( sum / (double)n );
}
 
 
 /* 
   return the maximum least squares error:
 
               (y[i] - y_hat[i])^2
  */
 double least_squares_max_error(double *x, double *y, int n, double m, double b){
     int i;
     double error, max_error;
 
     max_error = 0.0;
 
     for ( i = 0; i < n; i++ ) {
         error = y[i] - (m * x[i] + b);
         error = error * error;
         if ( error > max_error ) {
             max_error = error;
         }
     }
 
     return ( max_error );
 }