/**
* @file LinRegression.cpp
* @brief Algorithm for linear regression
* @author Frank Prüfer
* @date 08/13/2013

*/  

#include "vislearning/regression/linregression/LinRegression.h"
#include "core/vector/Algorithms.h"
#include "LinRegression.h"

using namespace OBJREC;

using namespace std;

using namespace NICE;

LinRegression::LinRegression()
{
  dim = 0;
}

LinRegression::LinRegression(uint dimension)
{
  dim = dimension;
}

LinRegression::LinRegression ( const LinRegression & src ) : 
RegressionAlgorithm ( src )
{
  dim = src.dim;
}

LinRegression::~LinRegression()
{
}

LinRegression* LinRegression::clone ( void ) const
{
  return new LinRegression(*this);
}


void LinRegression::teach ( const NICE::VVector & x, const NICE::Vector & y ){
  
  if (dim == 0)           //dimension not specified via constructor
  {
    dim = x[0].size()+1;  //use full dimension of data
  }
  
  cerr<<"dim: "<<dim<<endl;
  cerr<<"examples: "<<x.size()<<endl;
  
  for ( uint i = 0;i < dim;i++ )  //initialize vector of model parameters
  {
    modelParams.push_back(0.0);
  }
  
  if ( dim == 2 )         //two-dimensional least squares
  {  
    double meanX;
    double meanY = y.Mean();
    double sumX = 0.0;
    
    for ( uint i = 0;i < x.size();i++ )
      sumX += x[i][0];

    meanX = sumX / (double)x.size();
    
    for ( uint i = 0; i < x.size(); i++ )
      modelParams[1] += x[i][0] * y[i];

    modelParams[1] -= x.size() * meanX * meanY;
    
    double tmp = 0.0;
    for ( uint i = 0; i < x.size(); i++ )
      tmp += x[i][0] * x[i][0];

    tmp -= x.size() * meanX * meanX;
    
    modelParams[1] /= tmp;
    
    modelParams[0] = meanY - modelParams[1] * meanX;
  }
  else {  //N-dimensional least squares
    NICE::Matrix X, tmp, G;
    NICE::Vector params;
    
    x.toMatrix(X);
    
    NICE::Matrix Xtmp(X.rows(),X.cols()+1,1.0);
    
    // attach front column with ones

    for(uint row = 0;row<X.rows();row++)
    {
      for(uint col = 0;col<X.cols();col++)
      {
        Xtmp(row,col+1) = X(row,col);
      }
    }

    // modelParams =(X'X)^-1 * X'y
    NICE::Matrix tmpInv;
    NICE::Vector rhs;
      
    rhs.multiply(Xtmp,y,true);
    tmp.multiply(Xtmp,Xtmp,true);
      
    choleskyDecomp(tmp,G);
    choleskyInvert(G,tmpInv);

    params.multiply(tmpInv,rhs);

    modelParams = params.std_vector();
  }
}

std::vector<double> LinRegression::getModelParams(){
  return modelParams;
}

double LinRegression::predict ( const NICE::Vector & x ){
  double y;
  if ( dim == 2 )     //two-dimensional least squares
  {  
    y = modelParams[0] + modelParams[1] * x[0];
  }
  else {
    // y = x * modelParams
    NICE::Vector nModel(modelParams);
    NICE:: Vector xTmp(1,1.0);
    xTmp.append(x);
    y = xTmp.scalarProduct(nModel);
  }
  
  return y;
}