/**
* @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;    
  }
  else {  //N-dimensional least squares
    NICE::Matrix X, Xtmp, tmpAlpha, tmpAlphaInv, G;
    NICE::Vector params;
    uint X_rows,X_cols;
    
    x.toMatrix(X);
    X_rows = X.rows();
    X_cols = X.cols();
    cout<<X_rows<<" "<<X_cols<<endl;
    
    // attach front column with ones (UGLY)
    NICE::Vector oneVector(x.size(),1.0);
    Xtmp = X.transpose();
    cout<<Xtmp.rows()<<" "<<Xtmp.cols()<<endl;
    Xtmp.resize(X_cols+1,X_rows);
    cout<<Xtmp.rows()<<" "<<Xtmp.cols()<<endl;
    cout<<"Size of vector; "<<oneVector.size()<<endl;
    Xtmp.setRow(Xtmp.rows()-1,oneVector);
    Xtmp.exchangeRows(0,Xtmp.rows()-1);
    cout<<Xtmp.rows()<<" "<<Xtmp.cols()<<endl;
    X = Xtmp.transpose();
    
    cout<<X.rows()<<" "<<X.cols()<<endl;

    // alpha =(X'X)^-1 * X'y
    tmpAlpha.multiply(X,X,true);
    choleskyDecomp(tmpAlpha,G);
    choleskyInvert(G,tmpAlphaInv);
    params.multiply(X,y,true);
    params.multiply(tmpAlphaInv,params);

    alpha = params.std_vector();
    
    cout<<"Alpha-Vektor: ";
    for(uint i=0;i<alpha.size();i++){
      cout<<alpha[i]<<" ";
    }
    cout<<endl;
  }
}

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);
    y = x.scalarProduct(nAlpha);
  }
  
  return y;
}