HyperPointSet.cpp

Go to the documentation of this file.

#include "Mint/HyperPointSet.h"
#include <math.h>

HyperPointSet::HyperPointSet(int dimension) : 
  _dimension(dimension)
{
  WELCOME_LOG << "Hello from the HyperPointSet() Constructor";
}

HyperPointSet::HyperPointSet(TString path) :
  _dimension(-1)
{

  load(path);
 
}

HyperPointSet::HyperPointSet(int npoints, const HyperPoint& point) :
  _dimension(point.getDimension())
{

  for (int i = 0; i < npoints; i++) push_back(point);

}


HyperPointSet::HyperPointSet(const HyperPoint& point) :
  _dimension(point.getDimension())
{
  push_back(point);
}

HyperPointSet::HyperPointSet(const HyperPoint& point1, const HyperPoint& point2) :
  _dimension(point1.getDimension())
{
  push_back(point1);
  push_back(point2);
}

HyperPointSet::HyperPointSet(const HyperPoint& point1, const HyperPoint& point2, const HyperPoint& point3) :
  _dimension(point1.getDimension())
{
  push_back(point1);
  push_back(point2);
  push_back(point3);
}

bool HyperPointSet::compatible(const HyperPoint& other, bool printError) const{

  if (_dimension == other.getDimension()) return true;

  if (printError) ERROR_LOG << "This HyperPoint is NOT compatible with this HyperPointSet";
  return false;
}

const HyperPoint& HyperPointSet::at(int i) const      { 
  return _points.at(i);
}

HyperPoint& HyperPointSet::at(int i)                  { 
  return _points.at(i);
}

unsigned int HyperPointSet::size() const              { 
  return _points.size();
}


void HyperPointSet::push_back(const HyperPoint& point){ 
  if (compatible(point) == false) return;
  _points.push_back(point); 
} 


double HyperPointSet::getSumW() const{
  
  double sumW = 0.0;

  for(unsigned i = 0; i < size(); i++){
    sumW += at(i).getWeight();
  }
  return sumW;
}

double HyperPointSet::getSumW2() const{

  double sumW2 = 0.0;
  double w = 0.0;
  for(unsigned i = 0; i < size(); i++){
    w = at(i).getWeight();
    sumW2 += w*w;
  }
  return sumW2;

}

HyperPointSet HyperPointSet::gramSchmidtProcess() const{

  int nPoints = size();
  if ( nPoints > getDimension() ) nPoints = getDimension();
  
  if (linearlyIndependant() == 0) {
    ERROR_LOG << "Trying to perform the Gram Schmidt Process on a set of L.D. vectors";
  }

  HyperPointSet orthogonalPoints(at(0));

  for (int i = 1; i < nPoints; i++){
    HyperPoint nextPoint = at(i);
    for (int j = 0; j < i; j++){
      nextPoint = nextPoint - orthogonalPoints.at(j).project(at(i));
    }
    orthogonalPoints.push_back(nextPoint);
  }
  
  orthogonalPoints.normalise();
  return orthogonalPoints;

}

bool HyperPointSet::linearlyIndependant() const{
  //need to implement this...
  return 1;
}

void HyperPointSet::normalise(){

  for(unsigned i = 0; i < size(); i++){
    at(i) = at(i)/at(i).norm();
  }

}

void HyperPointSet::sort(){

  std::sort(_points.begin(), _points.end());

}

void HyperPointSet::removeDuplicates(){

  //INFO_LOG << "Removing all duplicates from the HyperPointSet - did you remember to call sort() first?!?";

  HyperPointSet hyperPointSet(getDimension());
  
  HyperPoint currentHyperPoint = _points.at(0); 
  hyperPointSet.push_back(currentHyperPoint);

  for(unsigned int i = 1; i < size(); i++){
    if ( _points.at(i) != currentHyperPoint ) {
      currentHyperPoint = _points.at(i); 
      hyperPointSet.push_back(currentHyperPoint);
    }
  }
  
  _points.clear();
  *this = hyperPointSet;

}

void HyperPointSet::addHyperPointSet(const HyperPointSet& other){
  
  if (_dimension != other._dimension){
    ERROR_LOG << "Removing all duplicates from the HyperPointSet - did you remember to call sort() first?!?";

  }

  for (unsigned int i = 0; i < other.size(); i++){
    this->push_back(other.at(i));
  }

}

double HyperPointSet::getCorrelation(int i, int j) const{

  WidthFinder statsi;
  WidthFinder statsj;

  for (unsigned evt = 0; evt < size(); evt++){
    statsi.add(at(evt).at(i));
    statsj.add(at(evt).at(j));
  }
  
  double meani  = statsi.mean();
  double meanj  = statsj.mean();
  double widthi = statsi.width();
  double widthj = statsj.width();

  double correlation = 0.0;

  for (unsigned evt = 0; evt < size(); evt++){
    correlation += (at(evt).at(i) - meani)*(at(evt).at(j) - meanj);
  } 
  
  correlation = correlation / (widthi * widthj * size());
  
  return correlation;
}

double HyperPointSet::getCovarience(int i, int j) const{

  MeanFinder statsi;
  MeanFinder statsj;

  for (unsigned evt = 0; evt < size(); evt++){
    statsi.add(at(evt).at(i));
    statsj.add(at(evt).at(j));
  }
  
  double meani  = statsi.mean();
  double meanj  = statsj.mean();

  double correlation = 0.0;

  for (unsigned evt = 0; evt < size(); evt++){
    correlation += (at(evt).at(i) - meani)*(at(evt).at(j) - meanj);
  } 
  
  correlation = correlation / size();
  
  return correlation;
}

TMatrixD HyperPointSet::getCovarienceMatrix() const{

  TMatrixD matrix(getDimension(), getDimension());

  for (int i = 0; i < getDimension(); i++){
    for (int j = 0; j < getDimension(); j++){
      matrix(i,j) = getCovarience(i,j);
    }    
  }
  
  return matrix;
}

TMatrixD HyperPointSet::getCorrelationMatrix() const{

  TMatrixD matrix(getDimension(), getDimension());

  for (int i = 0; i < getDimension(); i++){
    for (int j = 0; j < getDimension(); j++){
      matrix(i,j) = getCorrelation(i,j);
    }    
  }
  
  return matrix;
}

void HyperPointSet::print(std::ostream& os) const{
  
  os << "---------- HyperPointSet ------------" << std::endl;
  for (unsigned i = 0; i < size(); i++){
    os << "HyperPoint " << i << ":   ";
    at(i).print(os, 0);
    os << std::endl;
  }

}

void HyperPointSet::save(TString path){

  TFile* file = new TFile(path, "RECREATE");
  
  if (file == 0) {
    ERROR_LOG << "Cannot open file at " << path;
    return;
  }

  save();

  file->Write();
  file->Close();  

}

void HyperPointSet::save(){
 
  TTree* tree = new TTree("HyperPointSet", "HyperPointSet");
  
  if (tree == 0) {
    ERROR_LOG << "Cannot create tree";
    return;
  }

  std::vector<double>  values;
  std::vector<double>  weights;

  tree->Branch("values" , &values );
  tree->Branch("weights", &weights);
  
  for(unsigned i = 0; i < size(); i++ ){
    weights = at(i).getWeightsVector();
    values  = at(i).getVector();
    tree->Fill();
  }
  
  tree->ResetBranchAddresses();

  tree->Write();

}

//HyperPointSet, so multiple files can be loaded
void HyperPointSet::load(TString path){
  
  TFile* file = new TFile(path, "READ");
  if (file == 0) {
    ERROR_LOG << "Cannot open file at " << path;
    return;
  }

  TTree* tree = (TTree*)file->Get("HyperPointSet");
  if (tree == 0) {
    ERROR_LOG << "Cannot open tree";
    return;
  }
  
  std::vector<double>*  values  = 0;
  std::vector<double>*  weights = 0;

  tree->SetBranchAddress("values" , &values );
  tree->SetBranchAddress("weights", &weights);
  
  int nEntries = tree->GetEntries();

  for(int i = 0; i < nEntries; i++ ){
    tree->GetEntry(i);
    
    if (_dimension == -1) _dimension = values->size();
    if ( int(values->size()) != _dimension) ERROR_LOG << "You are loading HyperPoints that do not have the correct dimensionality for this HyperPointSet";

    HyperPoint point(*values);
    for(unsigned w = 0; w < weights->size(); w++){
      point.addWeight(weights->at(w));
    }
    this->push_back(point);
  }
  
  tree->ResetBranchAddresses();
  file->Close();

}

//HyperCuboid HyperPointSet::getLimits() const{
//
//  HyperCuboid limits(getDimension());
//
//  for (int i = 0; i < getDimension(); i++){
//    MinMaxFinder minMax;
//    for (unsigned j = 0; j < size(); j++){
//      minMax.add( this->at(j).at(i) );
//    }
//    limits.getLowCorner() .at(i) = minMax.getMin();
//    limits.getHighCorner().at(i) = minMax.getMax();
//  }
//  
//  return limits;
//}

HyperPoint HyperPointSet::getMin() const{
  HyperPoint val(getDimension());

  for (int i = 0; i < getDimension(); i++){
    MinMaxFinder minMax;
    for (unsigned j = 0; j < size(); j++){
      minMax.add( this->at(j).at(i) );
    }
    val.at(i) = minMax.getMin();
  }
  return val;
}

HyperPoint HyperPointSet::getMax() const{
  HyperPoint val(getDimension());

  for (int i = 0; i < getDimension(); i++){
    MinMaxFinder minMax;
    for (unsigned j = 0; j < size(); j++){
      minMax.add( this->at(j).at(i) );
    }
    val.at(i) = minMax.getMax();
  }
  return val;
}


HyperPoint HyperPointSet::mean() const{
  HyperPoint temp(getDimension());
  for (unsigned i = 0; i < size(); i++){
    temp = temp + at(i);
  }
  return temp/double(size());
}

HyperPoint HyperPointSet::geometricMean() const{

  HyperPoint geoMean(getDimension());
  for (unsigned i = 0; i < size(); i++){
    HyperPoint temp(getDimension());
    for (int j= 0; j < getDimension(); j++){
      temp.at(j) = log( this->at(i).at(j) );
    }
    geoMean = geoMean + temp;
  }
  
  geoMean = geoMean/double(size());

  for (int j= 0; j < getDimension(); j++){
    geoMean.at(j) = exp( geoMean.at(j) );
  }
  
  return geoMean;

}

HyperPoint HyperPointSet::harmonicMean() const{

  HyperPoint geoMean(getDimension());
  for (unsigned i = 0; i < size(); i++){
    HyperPoint temp(getDimension());
    for (int j= 0; j < getDimension(); j++){
      temp.at(j) = 1.0 / this->at(i).at(j);
    }
    geoMean = geoMean + temp;
  }
  
  geoMean = geoMean/double(size());

  for (int j= 0; j < getDimension(); j++){
    geoMean.at(j) = 1.0 / geoMean.at(j);
  }
  
  return geoMean;

}

HyperPointSet::~HyperPointSet(){
  GOODBYE_LOG << "Goodbye from the HyperPointSet() Destructor";
}