View of /trunk/gadget/regressionline.cc

#include "regressionline.h"
#include "mathfunc.h"
#include "errorhandler.h"
#include "gadget.h"
#include "global.h"

// ********************************************************
// Functions for Regression
// ********************************************************
Regression::Regression() {
  linetype = FREE;
  error = 1;
  useweights = 0;
  sse = a = b = 0.0;
}

Regression::Regression(LineType ltype) {
  linetype = ltype;
  error = 0;
  useweights = 0;
  sse = a = b = 0.0;
}

void Regression::setWeights(const DoubleVector& weights) {
  if (!useweights) {
    handle.logMessage(LOGWARN, "Warning in regression - unexpected use of weights");
    error = 1;
    return;
  }
  w = weights;
}

double Regression::getSSE() {
  if (error)
    return verybig;
  return sse;
}

void Regression::calcFit() {
  if (error)
    return;

  switch (linetype) {
    case FREE:
      //need to calculate both the slope and the intercept
      this->calcSlopeIntercept();
      break;
    case FIXEDSLOPE:
      //only need to calculate the intercept
      this->calcIntercept();
      break;
    case FIXEDINTERCEPT:
      //only need to calculate the slope
      this->calcSlope();
      break;
    case FIXED:
      //nothing to be done here
      break;
    default:
      handle.logMessage(LOGWARN, "Warning in regression - unrecognised linetype", linetype);
      break;
  }

  //finally calculate the SSE value
  if (useweights)
    this->calcSSEWeights();
  else
    this->calcSSE();
}

void Regression::calcSSE() {
  if (error)
    return;

  int i;
  double tmp;
  sse = 0.0;
  for (i = 0; i < x.Size(); i++) {
    tmp = y[i] - (a + b * x[i]);
    sse += tmp * tmp;
  }
}

void Regression::calcSSEWeights() {
  if (error)
    return;

  int i;
  double tmp;
  sse = 0.0;
  for (i = 0; i < x.Size(); i++) {
    tmp = y[i] - (a + b * x[i]);
    sse += w[i] * tmp * tmp;
  }
}

void Regression::calcSlope() {
  if (error)
    return;

  int i;
  double sumX, sumY;
  sumX = sumY = 0.0;
  for (i = 0; i < x.Size(); i++) {
    sumX += x[i];
    sumY += y[i];
  }
  if (isZero(sumX))
    b = 0.0;
  else
    b = (sumY - (a * x.Size())) / sumX;

  //JMB - if there is a negative slope for the regression then things are going wrong
  if (b < 0.0) {
    handle.logMessage(LOGWARN, "Warning in regression - negative slope for regression line", b);
    error = 1;
  }
}

void Regression::calcIntercept() {
  if (error)
    return;

  int i;
  double sumX, sumY;
  sumX = sumY = 0.0;
  for (i = 0; i < x.Size(); i++) {
    sumX += x[i];
    sumY += y[i];
  }
  a = (sumY - (b * sumX)) / x.Size();
}

void Regression::calcSlopeIntercept() {
  if (error)
    return;

  int i;
  double sumX, sumY, nom, denom;
  sumX = sumY = nom = denom = 0.0;
  for (i = 0; i < x.Size(); i++) {
    sumX += x[i];
    sumY += y[i];
  }
  sumX /= x.Size();
  sumY /= y.Size();

  for (i = 0; i < x.Size(); i++) {
    nom += (x[i] - sumX) * (y[i] - sumY);
    denom += (x[i] - sumX) * (x[i] - sumX);
  }

  if (isZero(denom)) {
    b = 0.0;
    a = sumY;
  } else {
    b = nom / denom;
    a = sumY - (b * sumX);
  }

  //JMB - if there is a negative slope for the regression then things are going wrong
  if (b < 0.0) {
    handle.logMessage(LOGWARN, "Warning in regression - negative slope for regression line", b);
    error = 1;
  }
}

// ********************************************************
// Functions for LinearRegression
// ********************************************************
LinearRegression::LinearRegression(LineType ltype) : Regression(ltype) {
}

void LinearRegression::storeVectors(const DoubleVector& modData, const DoubleVector& obsData) {
  error = 0;  //begin by cleaning up error status
  if ((modData.Size() != obsData.Size()) || (modData.Size() < 2)) {
    handle.logMessage(LOGWARN, "Warning in linear regression - invalid vector sizes");
    error = 1;
    return;
  }
  x = modData;
  y = obsData;
}

// ********************************************************
// Functions for LogLinearRegression
// ********************************************************
LogLinearRegression::LogLinearRegression(LineType ltype) : Regression(ltype) {
}

void LogLinearRegression::storeVectors(const DoubleVector& modData, const DoubleVector& obsData) {
  error = 0;  //begin by cleaning up error status
  if ((modData.Size() != obsData.Size()) || (modData.Size() < 2)) {
    handle.logMessage(LOGWARN, "Warning in log linear regression - invalid vector sizes");
    error = 1;
    return;
  }

  x.Reset();
  x.resize(modData.Size(), 0.0);
  y.Reset();
  y.resize(obsData.Size(), 0.0);

  int i, l = 0;
  for (i = 0; i < x.Size(); i++, l++) {
    if (isZero(modData[i]) && isZero(obsData[i])) {
      //omit the point (0.0, 0.0)
      x.Delete(l);
      y.Delete(l);
      l--;
    } else if ((modData[i] < verysmall) || (obsData[i] < verysmall)) {
      handle.logMessage(LOGWARN, "Warning in log linear regession - received invalid values");
      error = 1;
      return;
    } else {
      x[l] = log(modData[i]);
      y[l] = log(obsData[i]);
    }
  }
}

// ********************************************************
// Functions for WeightRegression
// ********************************************************
WeightRegression::WeightRegression(LineType ltype) : LinearRegression(ltype) {
  useweights = 1;
}

void WeightRegression::storeVectors(const DoubleVector& modData, const DoubleVector& obsData) {
  LinearRegression::storeVectors(modData, obsData);
  if (x.Size() != w.Size()) {
    handle.logMessage(LOGWARN, "Warning in weight regression - invalid vector sizes");
    error = 1;
  }
}

// ********************************************************
// Functions for LogWeightRegression
// ********************************************************
LogWeightRegression::LogWeightRegression(LineType ltype) : LogLinearRegression(ltype) {
  useweights = 1;
}

void LogWeightRegression::storeVectors(const DoubleVector& modData, const DoubleVector& obsData) {
  LogLinearRegression::storeVectors(modData, obsData);
  if (x.Size() != w.Size()) {
    handle.logMessage(LOGWARN, "Warning in log weight regression - invalid vector sizes");
    error = 1;
  }
}