CRSplineRegression is now able to handle multi-dimensional data.

regression/splineregression/CRSplineReg.cpp

@@ -5,10 +5,14 @@
 * @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"
 
@@ -71,109 +75,163 @@ double CRSplineReg::predict ( const NICE::Vector & x )
     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 );
   
-  if ( dataSet[0].size() == 1 ){	//one-dimensional case
-    FullVector data ( dataSet.size()+1 );
-    for ( uint i = 0; i < dataSet.size(); i++ ){
-      data[i] = dataSet[i][0];
-    }
-    cerr<<"data x: "<<x[0]<<endl;
-    data[dataSet.size()] = x[0];
-    
-    std::vector<int> ind;
-    data.getSortedIndices(ind);
-    
-    int index;
-    for ( uint i = 0; i < ind.size(); i++ ){
-      if ( ind[i] == dataSet.size() ){
-	index = i;
-	break;
-      }
-    }
+#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);
     
-    NICE::Matrix points (4,2,0.0);
-    if ( index >= 2 && index < (ind.size() - 2) ){	//everything is okay
-      points(0,0) = data[ind[index-2]];
-      points(0,1) = labelSet[ind[index-2]];
-      points(1,0) = data[ind[index-1]];
-      points(1,1) = labelSet[ind[index-1]];
-      points(2,0) = data[ind[index+1]];
-      points(2,1) = labelSet[ind[index+1]];
-      points(3,0) = data[ind[index+2]];
-      points(3,1) = labelSet[ind[index+2]];      
-    }
-    else if ( index == 1 ){	//just one point left from x
-      points(0,0) = data[ind[index-1]];
-      points(0,1) = labelSet[ind[index-1]];
-      points(1,0) = data[ind[index-1]];
-      points(1,1) = labelSet[ind[index-1]];
-      points(2,0) = data[ind[index+1]];
-      points(2,1) = labelSet[ind[index+1]];
-      points(3,0) = data[ind[index+2]];
-      points(3,1) = labelSet[ind[index+2]];       
-    }
-    else if ( index == 0 ){	//x is the farthest left point
-      //do linear approximation
+  int index;
+   
+  for ( uint i = 0; i < sortedInd.size(); i++ ){
+    if ( sortedInd[i] == dataSet.size() ){
+      index = i;
+      break;
     }
-    else if ( index == (ind.size() - 2) ){	//just one point right from x
-      points(0,0) = data[ind[index-2]];
-      points(0,1) = labelSet[ind[index-2]];
-      points(1,0) = data[ind[index-1]];
-      points(1,1) = labelSet[ind[index-1]];
-      points(2,0) = data[ind[index+1]];
-      points(2,1) = labelSet[ind[index+1]];
-      points(3,0) = data[ind[index+1]];
-      points(3,1) = labelSet[ind[index+1]];  
-    }
-    else if ( index == (ind.size() - 1) ){	//x is the farthest right point
-      //do linear approximation
-    }
-    
-    double t = (x[0] - points(1,0)) / (points(2,0) - points(1,0));
-    cerr<<"t: "<<t<<endl;
-    
-//     NICE::Vector vecT(4,1.0);
-//     
-//     vecT[1] = t;
-//     vecT[2] = t*t;
-//     vecT[3] = t*t*t;
-//     
-//     Matrix coeffMatrix (4,4,0.0);	// M = (0 2 0 0
-//     coeffMatrix(0,1) = 2.0;		//	-1 0 1 0
-//     coeffMatrix(1,0) = -1.0;		//	2 -5 4 -1
-//     coeffMatrix(1,2) = 1.0;		//	-1 3 -3 1)
-//     coeffMatrix(2,0) = 2.0;
-//     coeffMatrix(2,1) = -5.0;
-//     coeffMatrix(2,2) = 4.0;
-//     coeffMatrix(2,3) = -1.0;
-//     coeffMatrix(3,0) = -1.0;
-//     coeffMatrix(3,1) = 3.0;
-//     coeffMatrix(3,2) = -3.0;
-//     coeffMatrix(3,3) = 1.0;
+  }
+
+  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);
+//     }
 //     
-//     // P(t) = tau * vecT * coeffMatrix * points;
-//     NICE::Vector P;
-//     NICE::Matrix tmp;
-//     tmp.multiply(coeffMatrix,points);
-//     P.multiply(vecT,tmp);
-//     P *= tau;
+//     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;
     
-    //P(t) = b0*P0 + b1*P1 + b2*P2 + b3*P3
-    NICE::Vector P(2);
-    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);
+    }
     
-    P[0] = b0*points(0,0) + b1*points(1,0) + b2*points(2,0) + b3*points(3,0);
-    P[1] = b0*points(0,1) + b1*points(1,1) + b2*points(2,1) + b3*points(3,1);
-    
-    cerr<<"Response x: "<<P[0]<<endl;
-    cerr<<"Response y: "<<P[1]<<endl;
-    
-    return P[1];  
+    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;
+  
 }