///////////////////////////////////////////////////////////////////////////
// ASVR: Active Set Support Vector Regression                            //
// Copyright 2002 David R. Musicant and Alexander Feinberg               //
// This software is free for academic use only.                          //
// For commercial use, contact dmusican@carleton.edu.                    //
///////////////////////////////////////////////////////////////////////////
// asvrmex.cpp: Mex file for linking MATLAB with ASVR method.           //
///////////////////////////////////////////////////////////////////////////

#include <iostream.h>
#include <math.h>
#include "mex.h"
#include "matrix.h"
#include "asvrTrain.cpp"
#include "Version.H"

#include <iomanip>

#define PRINT mexPrintf
#define FLUSH
#define TRUE 1

// Wrapper class for MATLAB matrices
class MatlabInputMatrix {
 private:
  double *matrix;
 public:
  int m,n;

  MatlabInputMatrix(const mxArray *pointer) {
    matrix = (double *)mxGetPr(pointer);
    m = mxGetM(pointer);
    n = mxGetN(pointer);
  }

  inline double operator()(int i, int j) {
    return matrix[i+m*j];
  }
};

// Wrapper class for MATLAB matrices
class MatlabOutputMatrix {
private:
  double *matrix;
 public:
  int m,n;

  MatlabOutputMatrix(mxArray **pointer, int rows, int cols) {
    m = rows;
    n = cols;
    *pointer = mxCreateDoubleMatrix(m,n,mxREAL);
    matrix = (double *)mxGetPr(*pointer);
  }

  inline double operator()(int i, int j) {
    return matrix[i+m*j];
  }

  inline void Assign(int i,int j, double value) {
    matrix[i+m*j] = value;
  }
};

// Exceptions do not work as expected in a mexfile, so instead just
// spit out an error message
void ThrowIt(const char *text) {
  mexErrMsgTxt(text);
}


template <class T>
void matlabToMatrix(MatlabInputMatrix &in, MatrixTemplate<T> *out)
{
	
	DenseMatrixType<T> *outD;
	
	if ((in.m != out->rows) || (in.n != out->cols))
	{	throw (const char *)"matlabToMatrix: matrix sizes don't match!";	}
	
	if (outD = dynamic_cast<DenseMatrixType<T>*>(out))
	{
		int i, j, m = in.m, n = in.n;
		
		double **dataOut = outD->dataMatrix;
		
		for (i = 0; i < m; i++)
		{
			
			for (j = 0; j < n; j++)
			{
				dataOut[i][j] = in(i,j);
			}
		}
	}
	else
	{	throw (const char *)"matlabToMatrix: bad matrix type!";	}
}	



template <class T>
void matrixToMatlab(MatrixTemplate<T> *in, MatlabOutputMatrix &out)
{
	DenseMatrixType<T> *inD;
	
	if (inD = dynamic_cast<DenseMatrixType<T>*>(in))
	{
		int i, j, m = inD->rows, n = inD->cols;
		
		double **dataIn = inD->dataMatrix, *dataInRow;
		
		for (i=0; i<m; i++)
		{
			dataInRow = dataIn[i];
			
			for (j=0; j<n; j++)
			{	out.Assign(i, j, dataInRow[j]);	}
		}
	}
}



void matlabToVector(MatlabInputMatrix &in, Vector *out)
{
	DenseVector *outD;
	
	if (outD = dynamic_cast<DenseVector*>(out))
	{
		double *dataMatrix = outD->dataMatrix;
		
		int i, max;
		
		if (in.m > in.n)
		{
			max = in.m;
			
			for (int i = 0; i < max; i++)
			{	dataMatrix[i] = in(i,0);	}
		}
		else
		{
			max = in.n;
			
			for (int i = 0; i < max; i++)
			{	dataMatrix[i] = in(0,i);	}
		}
	}
	else
	{	throw (const char *)"matlabToVector: bad vector type!";	}
}



void vectorToMatlab(Vector *in, MatlabOutputMatrix &out)
{
	DenseVector *inD;
	
	if (inD = dynamic_cast<DenseVector*>(in))
	{
		if (in->cols > 0)
		{
			double *dataMatrix = inD->dataMatrix;
			
			int i, max;
			
			if (out.m > out.n)
			{
				max = out.m;
				
				for (int i = 0; i < max; i++)
				{	out.Assign(i,0, dataMatrix[i]);	}
			}
			else
			{
				max = out.n;
				
				for (int i = 0; i < max; i++)
				{	out.Assign(0,i, dataMatrix[i]);	}
			}
		}
	}
	else
	{	throw (const char *)"matlabToVector: bad vector type!";	}
}


int longSide(MatlabInputMatrix &m)
{
	if (m.m > m.n)
	{	return(m.m);	}
	else
	{	return(m.n);	}
}
	
		
void mexFunction(int nlhs,       mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
  MatlabInputMatrix dataIn(prhs[0]), yIn(prhs[1]);
  double eps = 0.1, C = 0.1, D = 1, printOn = 1.0;
  double tol = 1e-3, converge = 1,
    normalize = 0, prompt = 1;
  int maxiter = 100;
	
  switch(nrhs)
    {
    case 10:
      prompt = mxGetScalar(prhs[9]);
    case 9:
      normalize = mxGetScalar(prhs[8]);
    case 8:
      converge = mxGetScalar(prhs[7]);
    case 7:
      maxiter = (int)mxGetScalar(prhs[6]);
    case 6:
      tol = mxGetScalar(prhs[5]);
    case 5:
      D = (double)mxGetScalar(prhs[4]);
    case 4:
      C = (double)mxGetScalar(prhs[3]);
    case 3:
      eps = (double)mxGetScalar(prhs[2]);
      break;
    }
	
  DenseMatrixType<double> *A = new DenseMatrixType<double>(dataIn.m, dataIn.n);
  DenseVector *y = new DenseVector(longSide(yIn)),
    *w = new DenseVector(dataIn.n);
  double gamma, objective, time;
  int iters;
	
  matlabToMatrix(dataIn, A);
  matlabToVector(yIn, y);

  // Normalize if needed
  if (normalize) {
    double *mean = new double[A->cols];
    double *std = new double[A->cols];
    A->GetMeanStd(mean,std);
    A->Normalize(mean,std);
  }
	
  try
    {
      computeSolution(A, y, eps, C, D, w, gamma, objective,
		      false, (printOn == 0.0), time, maxiter,tol,iters);
    }
  catch (const char *err)
    {	mexErrMsgTxt(err);	}
	
  cout << "ASVR finished successfully.\n" << flush;
	
  MatlabOutputMatrix *wOut = NULL, *gammaOut = NULL, *alphaOut = NULL,
    *alphaHatOut = NULL, *objOut = NULL, *stepsOut = NULL,
    *timeOut = NULL;
	
  switch(nlhs)
    {
    case 7:
      timeOut = new MatlabOutputMatrix(&plhs[6],1,1);
      timeOut->Assign(0,0,time);

    case 6:
      stepsOut = new MatlabOutputMatrix(&plhs[5],1,1);
      stepsOut->Assign(0,0,iters);
	  
    case 5:
      objOut = new MatlabOutputMatrix(&plhs[4],1,1);
      objOut->Assign(0,0,objective);
	  
    case 4:
      alphaHatOut = new MatlabOutputMatrix(&plhs[3],0,0);
	  
    case 3:
      alphaOut = new MatlabOutputMatrix(&plhs[2],0,0);
	  
    case 2:
      gammaOut = new MatlabOutputMatrix(&plhs[1],1,1);
      gammaOut->Assign(0,0,gamma);
	  
    case 1:
    case 0:
      wOut = new MatlabOutputMatrix(&plhs[0],w->cols,1);
      vectorToMatlab(w, *wOut);
	  
      break;
    }
  delete(A);
  delete(y);
  delete(w);
	
  if (wOut != NULL)
    {	delete(wOut);	}
  if (gammaOut != NULL)
    {	delete(gammaOut);	}
  if (alphaOut != NULL)
    {	delete(alphaOut);	}
  if (alphaHatOut != NULL)
    {	delete(alphaHatOut);	}
  if (objOut != NULL)
    {	delete(objOut); 	}
  if (stepsOut != NULL)
    {	delete(stepsOut);	}
  if (timeOut != NULL)
    {	delete(timeOut);	}
}
