MultivariateStatistics.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include <float.h>
#include <string>
#include <iostream>
 
#include "IException.h"
#include "MultivariateStatistics.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
#include "PvlObject.h"
 
 
namespace Isis {
 
  MultivariateStatistics::MultivariateStatistics() {
    Reset();
  }
 
 
  MultivariateStatistics::MultivariateStatistics(const PvlObject &inStats) {
    Reset();
    fromPvl(inStats);
  }
 
 
  void MultivariateStatistics::Reset() {
    p_x.Reset();
    p_y.Reset();
    p_sumxy = 0.0;
 
    p_validPixels = 0;
    p_invalidPixels = 0;
    p_totalPixels = 0;
  }
 
 
  MultivariateStatistics::~MultivariateStatistics() {};
 
 
  void MultivariateStatistics::AddData(const double *x, const double *y,
                                       const unsigned int count) {
    for(unsigned int i = 0; i < count; i++) {
      double yVal = y[i];
      double xVal = x[i];
      p_totalPixels++;
 
      if(Isis::IsValidPixel(xVal) && Isis::IsValidPixel(yVal)) {
        p_x.AddData(xVal);
        p_y.AddData(yVal);
        p_sumxy += xVal * yVal;
        p_validPixels++;
      }
      else {
        p_invalidPixels++;
      }
    }
  }
 
 
  void MultivariateStatistics::AddData(double x, double y, unsigned int count) {
    p_totalPixels += count;
 
    if(IsValidPixel(x) && IsValidPixel(y)) {
      p_sumxy += x * y * count;
      p_validPixels += count;
 
      for (unsigned int i = 0; i < count; i++) {
        p_x.AddData(x);
        p_y.AddData(y);
      }
    }
    else {
      p_invalidPixels += count;
    }
  }
 
 
  void MultivariateStatistics::RemoveData(const double *x, const double *y,
                                          const unsigned int count) {
    for(unsigned int i = 0; i < count; i++) {
      p_totalPixels--;
 
      if(Isis::IsValidPixel(x[i]) && Isis::IsValidPixel(y[i])) {
        p_x.RemoveData(&x[i], 1);
        p_y.RemoveData(&y[i], 1);
        p_sumxy -= x[i] * y[i];
        p_validPixels--;
      }
      else {
        p_invalidPixels--;
      }
    }
 
    if(p_totalPixels < 0) {
      std::string m = "You are removing non-existant data in [MultivariateStatistics::RemoveData]";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }
  }
 
 
  double MultivariateStatistics::Covariance() const {
    if(p_validPixels <= 1) return Isis::NULL8;
    double covar = p_sumxy - p_y.Average() * p_x.Sum() - p_x.Average() * p_y.Sum() +
                   p_x.Average() * p_y.Average() * (double)p_validPixels;
    return covar / (double)(p_validPixels - 1);
  }
 
 
  double MultivariateStatistics::Correlation() const {
 
    if(p_validPixels <= 1) return Isis::NULL8;
    double covar = Covariance();
    double stdX = p_x.StandardDeviation();
    double stdY = p_y.StandardDeviation();
    if(stdX == 0.0 || stdX == Isis::NULL8) return Isis::NULL8;
    if(stdY == 0.0 || stdY == Isis::NULL8) return Isis::NULL8;
    if(covar == Isis::NULL8) return Isis::NULL8;
    return covar / (stdX * stdY);
  }
 
 
  BigInt MultivariateStatistics::TotalPixels() const {
    return p_totalPixels;
  }
 
 
  BigInt MultivariateStatistics::ValidPixels() const {
    return p_validPixels;
  }
 
 
  BigInt MultivariateStatistics::InvalidPixels() const {
    return p_invalidPixels;
  }
 
 
  void MultivariateStatistics::LinearRegression(double &a, double &b) const {
    // From Modern Elementary Statistics - 5th edition, Freund, pp 367
    double denom = (double)p_validPixels * p_x.SumSquare() - p_x.Sum() * p_x.Sum();
    if(denom == 0.0) {
      std::string msg = "Unable to compute linear regression in Multivariate Statistics";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    a = p_y.Sum() * p_x.SumSquare() - p_x.Sum() * p_sumxy;
    a = a / denom;
 
    b = (double)p_validPixels * p_sumxy - p_x.Sum() * p_y.Sum();
    b = b / denom;
  }
 
 
  double MultivariateStatistics::SumXY() const {
    return p_sumxy;
  }
 
 
  Isis::Statistics MultivariateStatistics::X() const {
    return p_x;
  };
 
 
  Isis::Statistics MultivariateStatistics::Y() const {
    return p_y;
  };
 
 
  void MultivariateStatistics::fromPvl(const PvlObject &inStats) {
    p_sumxy = inStats["SumXY"];
    p_validPixels = inStats["ValidPixels"];
    p_invalidPixels = inStats["InvalidPixels"];
    p_totalPixels = inStats["TotalPixels"];
 
    // unserialize the X and Y Statistics as well
    PvlGroup xStats = inStats.findGroup("XStatistics"); //
    p_x = Statistics(xStats);
    PvlGroup yStats = inStats.findGroup("YStatistics");
    p_y = Statistics(yStats);
  }
 
 
  PvlObject MultivariateStatistics::toPvl(QString name) const {
    if (name.isEmpty()) {
      name = "MultivariateStatistics";
    }
    PvlObject mStats(name); 
    mStats += PvlKeyword("Covariance" , toString(Covariance()));
    mStats += PvlKeyword("Correlation", toString(Correlation()));
    mStats += PvlKeyword("SumXY", toString(SumXY()));
    mStats += PvlKeyword("ValidPixels", toString(ValidPixels()));
    mStats += PvlKeyword("InvalidPixels", toString(InvalidPixels()));
    mStats += PvlKeyword("TotalPixels", toString(TotalPixels()));
 
    PvlKeyword linReg("LinearRegression");
    double a, b;
    try {
      LinearRegression(a, b);
      linReg += toString(a);
      linReg += toString(b);
    } catch (IException &e) {
      // It is possible one of the overlaps was constant and therefore
      // the regression would be a vertical line (x=c instead of y=ax+b)
    }
    mStats += linReg; 
 
    PvlGroup xStats = X().toPvl("XStatistics");
    PvlGroup yStats = Y().toPvl("YStatistics");
    mStats.addGroup(xStats);
    mStats.addGroup(yStats);
 
    return mStats;
  }
 
}