NICE_VisLearning/regression/splineregression/CRSplineReg.cpp

/**
* @file CRSplineReg.cpp
* @brief Implementation of Catmull-Rom-Splines for regression purposes
* @author Frank Prüfer
* @date 09/03/2013

*/  
#ifdef NICE_USELIB_OPENMP
#include <omp.h>
#endif

#include <iostream>

#include "vislearning/regression/splineregression/CRSplineReg.h"
#include "vislearning/regression/linregression/LinRegression.h"

#include "vislearning/math/mathbase/FullVector.h"

using namespace OBJREC;

using namespace std;
using namespace NICE;

CRSplineReg::CRSplineReg ( )
{
  tau = 0.5;
}

CRSplineReg::CRSplineReg ( const CRSplineReg & src ) : RegressionAlgorithm ( src )
{
  tau = src.tau;
  dataSet = src.dataSet;
  labelSet = src.labelSet;
}

CRSplineReg::~CRSplineReg()
{
}

void CRSplineReg::teach ( const NICE::VVector & _dataSet, const NICE::Vector & _labelSet)
{
    fprintf (stderr, "teach using all !\n");
    //NOTE this is crucial if we clear _teachSet afterwards!
    //therefore, take care NOT to call _techSet.clear() somewhere out of this method
    this->dataSet = _dataSet;
    this->labelSet = _labelSet.std_vector();
    
    std::cerr << "number of known training samples: " << this->dataSet.size() << std::endl;   
    
}

void CRSplineReg::teach ( const NICE::Vector & x, const double & y )
{
    std::cerr << "CRSplineReg::teach one new example" << std::endl;
    
    for ( size_t i = 0 ; i < x.size() ; i++ )
      if ( isnan(x[i]) ) 
      {
          fprintf (stderr, "There is a NAN value in within this vector: x[%d] = %f\n", (int)i, x[i]);
          cerr << x << endl;
          exit(-1);
      }

    dataSet.push_back ( x );
    
    labelSet.push_back ( y );
    
    std::cerr << "number of known training samples: " << dataSet.size()<< std::endl;
}

double CRSplineReg::predict ( const NICE::Vector & x )
{
  
  if ( dataSet.size() <= 0 ) {
    fprintf (stderr, "CRSplineReg: please use the train method first\n");
    exit(-1);
  }
  int dimension = dataSet[0].size();
  int sortDim = 0;

  FullVector data ( dataSet.size()+1 );
  
#pragma omp parallel for  
  for ( uint i = 0; i < dataSet.size(); i++ ){
    data[i] = dataSet[i][sortDim];
  }
  data[dataSet.size()] = x[sortDim];
    
  std::vector<int> sortedInd;
  data.getSortedIndices(sortedInd);
    
  int index;
   
  for ( uint i = 0; i < sortedInd.size(); i++ ){
    if ( sortedInd[i] == dataSet.size() ){
      index = i;
      break;
    }
  }

  NICE::Matrix points (4,dimension+1,0.0);
  if ( index >= 2 && index < (sortedInd.size() - 2) ){	//everything is okay
    points.setRow(0,dataSet[sortedInd[index-2]]);
    points(0,dimension) = labelSet[sortedInd[index-2]];
    points.setRow(1,dataSet[sortedInd[index-1]]);
    points(1,dimension) = labelSet[sortedInd[index-1]];      
    points.setRow(2,dataSet[sortedInd[index+1]]);
    points(2,dimension) = labelSet[sortedInd[index+1]];      
    points.setRow(3,dataSet[sortedInd[index+2]]);
    points(3,dimension) = labelSet[sortedInd[index+2]];           
  }
  else if ( index == 1 ){	//just one point left from x
    points.setRow(0,dataSet[sortedInd[index-1]]);
    points(0,dimension) = labelSet[sortedInd[index-1]];
    points.setRow(1,dataSet[sortedInd[index-1]]);
    points(1,dimension) = labelSet[sortedInd[index-1]];      
    points.setRow(2,dataSet[sortedInd[index+1]]);
    points(2,dimension) = labelSet[sortedInd[index+1]];      
    points.setRow(3,dataSet[sortedInd[index+2]]);
    points(3,dimension) = labelSet[sortedInd[index+2]];      
  }
  else if ( index == 0 ){	//x is the farthest left point
    points.setRow(0,dataSet[sortedInd[index+1]]);
    points(0,dimension) = labelSet[sortedInd[index+1]];
    points.setRow(1,dataSet[sortedInd[index+1]]);
    points(1,dimension) = labelSet[sortedInd[index+1]];      
    points.setRow(2,dataSet[sortedInd[index+1]]);
    points(2,dimension) = labelSet[sortedInd[index+1]];      
    points.setRow(3,dataSet[sortedInd[index+2]]);
    points(3,dimension) = labelSet[sortedInd[index+2]]; 
  }
  else if ( index == (sortedInd.size() - 2) ){	//just one point right from x
    points.setRow(0,dataSet[sortedInd[index-2]]);
    points(0,dimension) = labelSet[sortedInd[index-2]];
    points.setRow(1,dataSet[sortedInd[index-1]]);
    points(1,dimension) = labelSet[sortedInd[index-1]];      
    points.setRow(2,dataSet[sortedInd[index+1]]);
    points(2,dimension) = labelSet[sortedInd[index+1]];      
    points.setRow(3,dataSet[sortedInd[index+1]]);
    points(3,dimension) = labelSet[sortedInd[index+1]];   
  }
  else if ( index == (sortedInd.size() - 1) ){	//x is the farthest right point
    points.setRow(0,dataSet[sortedInd[index-2]]);
    points(0,dimension) = labelSet[sortedInd[index-2]];
    points.setRow(1,dataSet[sortedInd[index-1]]);
    points(1,dimension) = labelSet[sortedInd[index-1]];      
    points.setRow(2,dataSet[sortedInd[index-1]]);
    points(2,dimension) = labelSet[sortedInd[index-1]];      
    points.setRow(3,dataSet[sortedInd[index-1]]);
    points(3,dimension) = labelSet[sortedInd[index-1]];     
  }

//   NICE::Vector P(dimension);
//   double y;
//   double diff,bestT;
//   double bestDiff = std::numeric_limits<double>::max();
//   double b0,b1,b2,b3;
// 
//   for ( double t = 0.0; t <= 1.0; t += 0.02 ){
//     b0 = tau * (-(t*t*t) + 2*t*t - t);
//     b1 = tau * (3*t*t*t - 5*t*t + 2);
//     b2 = tau * (-3*t*t*t + 4*t*t + t);
//     b3 = tau * (t*t*t - t*t);    
// 
// #pragma omp parallel for  
//     for ( uint i = 0; i < dimension; i++ ){
//       P[i] = b0*points(0,i) + b1*points(1,i) + b2*points(2,i) + b3*points(3,i);
//     }
//     
//     diff = (P-x).normL2();
//     if ( diff < bestDiff ){
//       bestDiff = diff;
//       bestT = t;
//     }
//   }
//   b0 = tau * (-(bestT*bestT*bestT) + 2*bestT*bestT - bestT);
//   b1 = tau * (3*bestT*bestT*bestT - 5*bestT*bestT + 2);
//   b2 = tau * (-3*bestT*bestT*bestT + 4*bestT*bestT + bestT);
//   b3 = tau * (bestT*bestT*bestT - bestT*bestT); 

  double t = (x[sortDim]-points(1,sortDim)) / (points(2,sortDim)-points(1,sortDim));
//   cerr<<"t: "<<t<<endl;

  //P(t) = b0*P0 + b1*P1 + b2*P2 + b3*P3    
  NICE::Vector P(dimension);
  double y;
  double b0,b1,b2,b3;
    
  b0 = tau * (-(t*t*t) + 2*t*t - t);
  b1 = tau * (3*t*t*t - 5*t*t + 2);
  b2 = tau * (-3*t*t*t + 4*t*t + t);
  b3 = tau * (t*t*t - t*t);

#pragma omp parallel for  
  for ( uint i = 0; i < dimension; i++ ){
    P[i] = b0*points(0,i) + b1*points(1,i) + b2*points(2,i) + b3*points(3,i);
  }
  
  double diff1 = (P-x).normL2();
  uint counter = 1;
  while ( diff1 > 1e-5 && counter <= 21){
    double tmp = t;;
    if (tmp > 0.5)
      tmp = 1 - tmp;
    t += tmp/counter;
     
    b0 = tau * (-(t*t*t) + 2*t*t - t);
    b1 = tau * (3*t*t*t - 5*t*t + 2);
    b2 = tau * (-3*t*t*t + 4*t*t + t);
    b3 = tau * (t*t*t - t*t);
      
    for ( uint i = 0; i < dimension; i++ ){
      P[i] = b0*points(0,i) + b1*points(1,i) + b2*points(2,i) + b3*points(3,i);
    }
    
    double diff2 = (P-x).normL2();
    if ( diff2 > diff1 && t > 0) {
      t -= 2*tmp/counter;    
	
      b0 = tau * (-(t*t*t) + 2*t*t - t);
      b1 = tau * (3*t*t*t - 5*t*t + 2);
      b2 = tau * (-3*t*t*t + 4*t*t + t);
      b3 = tau * (t*t*t - t*t);

#pragma omp parallel for      
      for ( uint i = 0; i < dimension; i++ ){
	P[i] = b0*points(0,i) + b1*points(1,i) + b2*points(2,i) + b3*points(3,i);
      }
      diff1 = (P-x).normL2();
    }
    counter++;
  }
  
  y = b0*points(0,dimension) + b1*points(1,dimension) + b2*points(2,dimension) + b3*points(3,dimension);
  return y;
  
}