/**
* @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"

using namespace OBJREC;

using namespace std;

using namespace NICE;

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

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

LinRegression::~LinRegression()
{
}

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 alpha-vector
    alpha.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++ ){
	alpha[1] += x[i][0] * y[i];
    }
    
    alpha[1] -= x.size() * meanX * meanY;
    
    double tmpAlpha = 0.0;
    for ( uint i = 0; i < x.size(); i++ ){
      tmpAlpha += x[i][0] * x[i][0];
    }
    tmpAlpha -= x.size() * meanX * meanX;
    
    alpha[1] /= tmpAlpha;
    
    alpha[0] = meanY - alpha[1] * meanX;
    
//     cerr<<"Alpha-Vektor: ";
//     for(uint i=0;i<alpha.size();i++){
//       cerr<<alpha[i]<<" ";
//     }
//     cerr<<endl;
  }
  else {  //N-dimensional least squares
    NICE::Matrix X, tmpAlpha, 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);
      }
    }

    // alpha =(X'X)^-1 * X'y
    NICE::Matrix tmpAlphaInv;
    NICE::Vector rhs;
      
    rhs.multiply(Xtmp,y,true);
      
    tmpAlpha.multiply(Xtmp,Xtmp,true);
      
    choleskyDecomp(tmpAlpha,G);
    choleskyInvert(G,tmpAlphaInv);
      
    params.multiply(tmpAlphaInv,rhs);

    alpha = params.std_vector();
  }
}

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

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