/***************************************************************************
 * 
 * Augustus Lidaka, Shitanshu Aggarwal, and Jerod Weinman
 *
 *
 * CUDA_VALUE.CU:
 *
 *    This program coordinates the CUDA calls for evaluating value and
 *    gradient for the objective function. It is the interface between
 *    MATLAB and the GPU device.
 *
 *
 * Preconditions:
 *
 *    This function should recieve pointers to input variables and output
 *    variables from MATLAB. The input variables should remain constant
 *    while this function is being invoked. The input variable G should be
 *    initialized with the observed energy Gobs. eps should be the machine
 *    epsilon of single point representation calculated by MATLAB.
 *
 *
 * Postconditions:
 *
 *     The program calculates the value and the gradient of training
 *     instances given a weight vector on the GPU. This program will handle
 *     errors and set the value accordingly if necessary. The gradient is
 *     calculated iteratively by subtracting a partial calculated gradient
 *     from the input variable G. 
 * 
 *
 * Legend:
 *   
 *    N  -  Number of instances
 *    D  -  Number of features
 *    L  -  Number of labels
 *
 *
 * Input Parameters:
 * 
 *    const float *X      Feature vector, contains response of each instance
 *                        to each feature. Size: NxD
 *    const unsigned int *Y  Label vector, contains the true label of each
 *                           instance. Size: Nx1
 *    const float *W      Weights matrix, contains the weight of each label
 *                        for every feature . Size: (D+1)xL
 *    float *G            Gradient matrix, initialized with the observed
 *                        "energy". Size: DxL
 *    const float eps     Matlab's internal epsilon value. Should be machine
 *                        epsilon for single precision floats.
 *
 *    Optional:
 *    const float *I      Instance weight matrix (optional), contains the
 *                        weights for each instance. Size: Nx1
 *    bool debug          This is a flag to toggle debugging bhavior.
 *
 * Output:
 *
 *    float V             Value, we are trying to maximize this
 *                        variable. Size: 1x1
 *    float* G            Gradient matrix, contains the difference between
 *                        the observed and the calculated energy. Size: DxL 
 *    float *gradient_bias  Vector containing the gradient on the 
 *                          bias weights.
 *
 **************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include "cublas.h"
#include <unistd.h>

#ifndef __TEST_DRIVER__
#include "mex.h"
#endif

#include "cuda_value.h"
#include "parmarginals.cu"



typedef struct {
  float *X;
  unsigned int *Y;
  float *I;
  unsigned int N; //Number of instances in this chunk
  int L; //Number of labels
  int D; //Number of features
} slice_t;



/* sliceInstances increments the "breadloaf"; it modifies its first
 * parameter, breadloaf, to point to the next "slice" we will use, based on
 * information in the slice_t struct containing pointers to the current set
 * of data, slice. */
void sliceInstances(slice_t *breadloaf, slice_t *slice, int instances)
{

  // Setting up the slice
  *slice = *breadloaf;


  if (instances >= breadloaf->N)
  {
    //This is the last chunk of bread
    breadloaf->X = NULL;
    breadloaf->Y = NULL;
    breadloaf->I = NULL;
    breadloaf->N = 0;
  }

  else
  {
    slice->N = instances;

    //Incrementing the head of the breadloaf
    breadloaf->X += slice->D * slice->N;
    breadloaf->Y += slice->N;
    if (NULL != breadloaf->I)
      breadloaf->I += slice->N;
    breadloaf->N -= slice->N;
  }

}



/* Function recieves the slice and three pointers which hold device memory
   addresses. Allocates the slice struct to the pointers to the memory on
   the device. */
bool copySliceToDevice(slice_t *slice, float *Xd, unsigned int *Yd,
                       float *Id) 
{    

  cudaError status;
  status = cudaMemcpy(Xd, slice->X, slice->D * slice->N * sizeof(float),
                      cudaMemcpyHostToDevice);
  if (status != cudaSuccess)
  {
    printf("Could not copy X chunk to device: %s\n",
           cudaGetErrorString(status));
    return 0;
  }

  status = cudaMemcpy(Yd, slice->Y,
                      slice->N * sizeof(unsigned int),
                      cudaMemcpyHostToDevice);
  if (status != cudaSuccess)
  {
    printf("Could not copy Y chunk to device: %s!\n",
           cudaGetErrorString(status));
    return 0;
  }

  if (NULL != slice->I)
  {
    status = cudaMemcpy(Id, slice->I, slice->N * sizeof(float),
                        cudaMemcpyHostToDevice);
    if (status != cudaSuccess)
    {
      printf("Could not copy I chunk to device: %s\n",
             cudaGetErrorString(status));
      return 0;
    }
  }
  else
    Id = NULL;

  return 1;

}



/* Recieves pointers pointing to the pointers which will hold adresses to
 * the memory allocated on the device. Then mallocs stuff on the device and
 * populates the pointers.*/
bool mallocDeviceWorkspace(unsigned int N, unsigned int D, unsigned int L, 
     	                   float **Xd, unsigned int **Yd, 
			   float **Id, float **Wd, float **Bd, float **Gd, 
			   float **vv, float **Vd, float **weight_bias_d,
			   float **ones_d, float **gradient_bias_d)
{
    cudaError status;

    /* NOTE: CUDA-Malloc fails to report device memory limits. Caution must 
     * be taken	taken while using this function. The total memory should not
     * exceed 3.8GB */	
    status = cudaMalloc((void **)Xd,D*N*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("\ncudaMalloc failed for Xd: %s\n", cudaGetErrorString(status));
      return 0;
    }
       
    status = cudaMalloc((void **)Wd,L*D*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for Wd: %s\n", cudaGetErrorString(status));
      return 0;
    }
    
    status = cudaMalloc((void **)weight_bias_d,L*1*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for weight_bias_d: %s\n",
             cudaGetErrorString(status));
      return 0;
    }
 
    status = cudaMalloc((void **)ones_d,1*N*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for ones_d: %s\n",
             cudaGetErrorString(status));
      return 0;
    }
    
    status = cudaMalloc((void **)Bd,L*N*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for Bd: %s\n", cudaGetErrorString(status));
      return 0;
    }
    
    status = cudaMalloc((void **)Yd,1*N*sizeof(unsigned int));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for Yd: %s\n", cudaGetErrorString(status));
      return 0;
    }

    if (NULL != *Id)
    {
      status = cudaMalloc((void **)Id,1*N*sizeof(float));
      if (status != cudaSuccess)
      {
        printf("cudaMalloc failed for Id: %s\n", cudaGetErrorString(status));
        return 0;
      }
    }

    status = cudaMalloc((void **)vv,1*N*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for vv: %s\n", cudaGetErrorString(status));
      return 0;
    }

    status = cudaMalloc((void **)Vd,1*1*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for Vd: %s\n", cudaGetErrorString(status));
      return 0;
    }

    status = cudaMalloc((void **)Gd,L*D*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for Gd: %s\n", cudaGetErrorString(status));
      return 0;
    }

    status = cudaMalloc((void **)gradient_bias_d, L*1*sizeof(float));
    if (status != cudaSuccess)
    {
      printf("cudaMalloc failed for gradient_bias_d: %s\n",
             cudaGetErrorString(status));
      return 0;
    }


    return 1;

}



/* Free all the memory referenced by the input pointers */
void freeDeviceWorkspace(float *Xd, unsigned int *Yd, float *Id, float *Wd, 
                         float *Bd, float *Gd, float *vv, float *Vd,
                         float *weight_bias_d, float *ones_d, 
                         float *gradient_bias_d)
{

  cudaError_t err;
  
	err = cudaFree(Xd);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
	err = cudaFree(Yd);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
  if (NULL != Id)
  {
    err = cudaFree(Id);
    if (err != cudaSuccess)
        printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
  }
 	err = cudaFree(Wd);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
 	err = cudaFree(Bd);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
 	err = cudaFree(Gd);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
 	err = cudaFree(vv);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
 	err = cudaFree(Vd);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
 	err = cudaFree(weight_bias_d);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
 	err = cudaFree(ones_d);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));
	err = cudaFree(gradient_bias_d);
	if (err != cudaSuccess)
	    printf("Error freeing workspace: %s.\n", cudaGetErrorString(err));

	return;

}



/* Executes the calculation of the raw energy on the device */ 

/* NOTE: if the debug flag is raised, the function returns a pointer to the 
 * raw energy copied to the host */
float* rawEnergy(float *Bd, float *Wd, float *Xd, const float *Wh, 
                 float *weight_bias_d, float *ones_d, int N, int D, int L,
                 bool debug)
{

  /* Fetch weight_bias_d and generate ones_d */
  cublasStatus status;
  const float *weight_bias_h = Wh + (L*D);
    
  /* Copy weights bias viector to device */
  status = cublasSetVector(L, sizeof(float), weight_bias_h, 1, 
			                     weight_bias_d, 1);
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("Failed to copy weights bias vector to device: code %d\n",status);
    return NULL;
  }
  
  /* Clear last error */
  cublasGetError();
  
  /* Initializing Bd with bias vector */
  cublasSgemm('n','n', L, N, 1, 1, weight_bias_d, L, ones_d, 1, 0, Bd, L);
  cudaThreadSynchronize();
  status = cublasGetError();
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("cublasSgemm failed initializing B: code %d\n", status);
    return NULL;
  }


  
  /* Calculate [raw] energy */
  cublasSgemm('n','n', L, N, D, 1, Wd, L, Xd, D, 1, Bd, L);
  cudaThreadSynchronize();
  status = cublasGetError();
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("CUBLAS kernel execution failed while calculating energy!\n");
    return NULL;
  }

  if (debug) {
  
    /* Copy energy matrix from device to host*/
    
    float *Bh = (float *)malloc(sizeof(float)*L*N);
    if (!Bh)
    {
      printf("Could not allocate host memory to store result B!\n");
      return NULL;
    }
    
    status = cublasGetMatrix(L, N, sizeof(float), Bd, L, Bh, L);
    if (status != CUBLAS_STATUS_SUCCESS)
    {
      printf("Could not copy result B from device to host!\n");
      return NULL;
    }
    
    return Bh;
  }
  
  return NULL;

}
 


/* Executes the calculation of the gradient. */ 
void calculateGradient(float *Gd, float *Bd, float *Xd, float *ones_d, 
		       float *gradient_bias_d, unsigned int N, int D, int L)
 {

   cublasStatus status;

   /* Clear last error */
   cublasGetError();

  /* Calculate gradient on device */
  cublasSgemm('n','t', D, L, N, -1, Xd, D, Bd, L, 1, Gd, D);
  cudaThreadSynchronize();
  status = cublasGetError();
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("CUBLAS kernel execution failed while calculating gradient!\n");
    return;
  }



  /* Calculating bias feature on device */
  cublasSgemv('n', L, N, -1, Bd, L, ones_d, 1, 1, gradient_bias_d, 1);
  cudaThreadSynchronize();
  status = cublasGetError();
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("CUBLAS kernel execution failed while calculating bias vector!\n");
    return;
  }

 return;
}



/* cuda_value does basically all the work, including setting the output
 * (presumably it is passed pointers to the mexFunction output variables).
 *
 * cuda_value(X, Y, W, Gobs, gradient_bias_in, N, D, L, eps, I, G, V,
 *            gradient_bias)
 */
void cuda_value(float *X, unsigned int *Y, float *W, float *Gobs, float
                *gradient_bias_in, unsigned int N, unsigned int D, unsigned
                int L, float eps, float *I, float *G, float *V, float
                *gradient_bias)
{

  // Misc variables
  float *ones_h;
  cudaError_t err;
  int i;
  cublasStatus status;



  // Initialize CUBLAS
  status = cublasInit();
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("CUBLAS initialization error!\n");
    return;
  }



  /* Get device properties (total memory and multiprocessor count) */
  int dev;
  cudaDeviceProp *prop;

  err = cudaGetDevice(&dev);
  if (err != CUDA_SUCCESS)
  {
    printf("Could not get CUDA device ID: %s\n", cudaGetErrorString(err));
    return;
  }

  prop = (cudaDeviceProp *)malloc(sizeof(cudaDeviceProp));
  if (NULL == prop)
  {
    printf("Error malloc'ing device property struct prop.\n");
    return;
  }

  err = cudaGetDeviceProperties(prop, dev);
  if (err != CUDA_SUCCESS)
  {
    printf("Error getting device properties: %s\n", cudaGetErrorString(err));
    return;
  }

  size_t device_total = prop->totalGlobalMem;
  /* We don't want to use too close to 100% memory; could cause unforseeable
   * problems. */
  device_total *= .90;

  free(prop);
  prop = NULL;



  /* Calculate size of slice:
   * The dimensions of all the variables to be allocated on the device are
   * added up and set equal to the available memory to calculate the largest 
   * N for a sigle slice of instances.
   *
   * The equation to calculate Nslices is a rearrangement of the equation:
   * 
   * memory = sizeof(float)* dimensions( X + Y + W + I + B + G + vv +
   *                                     V + ones_d + gradient_bias_d + 
   *                                     weights_bias_d )
   * i.e.,
   * device_total = sizeof(float) * ( D*N + N*1 + L*D + N + 
   *                                  L*N + L*D + N + 1 + N + L + L )
  */
  unsigned int Nslice;
  Nslice = (device_total/sizeof(float) - 2*L*(D+1) - 1) / (D + L + 4);
  if (Nslice <= 0)
  {
    printf("Not enough memory for static-sized variables!\n");
    return;
  }

  // Hardlimit on size of slice (each instance gets one block)
  // Potential bottleneck, since we could be excessively doing memory transfers
  if (Nslice > MAX_NUM_BLOCKS)
    Nslice = MAX_NUM_BLOCKS;



  /* Allocate Memory on device */
  float *Xd, *Id, *Wd, *Bd, *Gd, *vv, *Vd;
  float *weight_bias_d, *ones_d, *gradient_bias_d;
  unsigned int *Yd;

  if (NULL == I)
    Id = NULL;

  if(1 != mallocDeviceWorkspace(Nslice, D, L, &Xd, &Yd, &Id, &Wd, &Bd, &Gd,
                                &vv, &Vd, &weight_bias_d, &ones_d,
                                &gradient_bias_d))
  {
    printf("Could not allocate workspace device memory!\n");
    return;
  }



  /* Setting up breadloaf */
  slice_t breadloaf = {X,Y,I,N,L,D};
  slice_t slice;




  /* Copy constant, persistent variables to device */
  // Weights W
  status = cublasSetMatrix(L, D, sizeof(float), W, L, Wd, L);
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("Could not copy W from host to device memory!\n");
    return;
  }

  // Initialize G with Gobs
  status = cublasSetMatrix(D, L, sizeof(float), Gobs, D, Gd, D);
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("Could not copy Gobs from host to device memory!\n");
    return;
  }

  // Construct and copy a ones vector
  ones_h = (float *)malloc(sizeof(float) * N);
  for (i = 0; i < N; i++)
    ones_h[i] = 1;

  status = cublasSetVector(N, sizeof(float), ones_h, 1, ones_d, 1);
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("Could not copy ones_h from host to device memory!\n");
    return;
  }

  // Gradient bias
  status = cublasSetVector(L, sizeof(float), gradient_bias_in, 1,
                           gradient_bias_d, 1);
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("Could not copy ones_h from host to device memory!\n");
    return;
  }



  /* Initialize several more variables before looping over slices */
  unsigned int num_slices = N/Nslice + (N % Nslice ? 1 : 0);
  unsigned int block_size;
  int *error_d;
  int error_h = 0;
  float *vv_d = NULL;
  float Vh;
  *V = 0;



  /* Loop over slices */
  for (i=0; i < num_slices; i++)
  {

    // Slice a chunk
    sliceInstances(&breadloaf, &slice, Nslice);


        
    // Copy to device
    status = copySliceToDevice(&slice, Xd, Yd, Id);
    if (1 != status)
    {
      printf("copySliceToDevice failed!\n");
      *V = NAN;
      return;
    }



    // rawEnergy
    rawEnergy(Bd, Wd, Xd, W, weight_bias_d, ones_d, slice.N, D, L, 0);



    // parmarginals
    // calculate block size (thread per block)
    if (L < MAX_THREADS_PER_BLOCK)
      block_size = L;
    else
      block_size = MAX_THREADS_PER_BLOCK;

    // make space for vv_d, one of the side effects of parmarginals
    err = cudaMalloc((void **)&vv_d, sizeof(float)*slice.N);
    if (err != cudaSuccess)
    {
      printf("Error cudaMalloc'ing vv_d: %s\n", cudaGetErrorString(err));
      return;
    }

    // make space for the error flag on the device
    err = cudaMalloc((void **)&error_d, sizeof(int));
    if (err != cudaSuccess)
    {
      printf("Error cudaMalloc'ing error_d: %s\n", cudaGetErrorString(err));
      return;
    }

    // copy over the error flag as it is now, i.e., set it to 0
    err = cudaMemcpy(error_d, &error_h, sizeof(int), cudaMemcpyHostToDevice);
    if (err != cudaSuccess)
    {
      printf("Error cudaMemcpy'ing error_d: %s\n", cudaGetErrorString(err));
      return;
    }

    // call parmarginals
    parmarginals <<< slice.N, block_size >>>
                  (Bd, Yd, Id, slice.N, L, eps, vv_d, error_d);
    cudaThreadSynchronize();
    err = cudaGetLastError();
    if (err != cudaSuccess)
    {
      printf("Error in parmarginals: %s\n", cudaGetErrorString(err));
      return;
    }

    /* Check for errors (detected NaN) values in the marginals */
    err = cudaMemcpy(&error_h, error_d, sizeof(int), cudaMemcpyDeviceToHost);
    if (err != cudaSuccess)
    {
      printf("Error cudaMemcpy'ing error_h: %s\n", cudaGetErrorString(err));
      return;
    }
    else if (error_h)
    {
      /* If error_h is set. */
      printf("Error: ");

      switch (error_h)
      {
      case UINT_MAX_IS_NAN:          printf("UINT_MAX_IS_NAN\n"); break;
      case FLOAT_MAX_IS_NAN:         printf("FLOAT_MAX_IS_NAN\n"); break;
      case THREAD_Z_SANS_MAX_IS_NEG: printf("THREAD_Z_SANS_MAX_IS_NEG\n"); break;
      case THREAD_Z_SANS_MAX_IS_INF: printf("THREAD_Z_SANS_MAX_IS_INF\n"); break;
      case EXP_Z_SANS_MAX_IS_NAN:    printf("EXP_Z_SANS_MAX_IS_NAN\n"); break;
      case EXP_Z_SANS_MAX_IS_ZERO:   printf("EXP_Z_SANS_MAX_IS_ZERO\n"); break;
      case EXP_Z_SANS_MAX_IS_NEG:    printf("EXP_Z_SANS_MAX_IS_NEG\n"); break;
      case EXP_Z_SANS_MAX_IS_INF:    printf("EXP_Z_SANS_MAX_IS_INF\n"); break;
      case LOG_Z_SANS_MAX_IS_NAN:    printf("LOG_Z_SANS_MAX_IS_NAN\n"); break;
      case LOG_Z_IS_NAN:             printf("LOG_Z_IS_NAN\n"); break;
      case LOG_B_IS_NAN:             printf("LOG_B_IS_NAN\n"); break;
      default:
        /* Unknown error */
        printf("Unknown error in parmarginals: %d.\n", error_h);
      }
      
      /*  In that case, return NaN value. The gradient will be
       * meaningless, so use the regularization gradient (done in
       * parvalue). */
      
      cudaFree(error_d);
      cudaFree(vv_d);
      freeDeviceWorkspace(Xd, Yd, Id, Wd, Bd, Gd, vv, Vd, weight_bias_d,
                          ones_d, gradient_bias_d);
      free(ones_h);

      *V = NAN;
      return;
    }
    
    // free error_d
    cudaFree(error_d);
    


    /** calculate total value for slice **/

    if (1 || NULL==I)
      // Since all values are log(x), x\in[0,1], the values should be negative.
      // Thus, the absolute value is always of a negative number, the abs sum
      // has a sign that simply needs to be flipped.
      Vh =  - cublasSasum( slice.N, vv_d, 1);
    else
      Vh = cublasSdot( slice.N, vv_d, 1, Id, 1);


    // free vv_d
    cudaFree(vv_d);
    if (err != cudaSuccess)
    {
      printf("Error cudaFree'ing vv_d: %s.\n", cudaGetErrorString(err));
      return;
    }
    vv_d = NULL;



    // calcGradient
    calculateGradient(Gd, Bd, Xd, ones_d, gradient_bias_d, slice.N, D, L);





    // accumulate V 
    // (NOTE: G & gradient_bias are already being updated on the GPU)
    *V += Vh;

  }



  // Copy final G back to host
  err = cudaMemcpy(G,Gd,sizeof(float)*D*L,cudaMemcpyDeviceToHost);
  if (err != cudaSuccess)
  {
    printf("Error cudaMemcpy'ing G: %s.\n", cudaGetErrorString(err));
    return;
  }



  // Copy final gradient bias vector to host
  cudaMemcpy(gradient_bias,gradient_bias_d, sizeof(float)*L*1, 
	           cudaMemcpyDeviceToHost);
  if (err != cudaSuccess)
  {
    printf("Error cudaMemcpy'ing gradient_bias: %s.\n",
           cudaGetErrorString(err));
    return;
  }
  


  /* Shutdown CUBLAS */
  status = cublasShutdown();
  if (status != CUBLAS_STATUS_SUCCESS)
  {
    printf("Error calling cublasShutdown!\n");
    return;
  }
  


  // Free memory
  freeDeviceWorkspace(Xd, Yd, Id, Wd, Bd, Gd, vv, Vd, weight_bias_d, ones_d,
                      gradient_bias_d);
  free(ones_h);



  return;
   
}



#ifndef __TEST_DRIVER__
/* This function contains the 4 kernel calls */
/* [V G gradient_bias] = cuda_value(X, Y, W, Gobs, gradient_bias_in, eps, I,
 *                                  debug) */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

  // Dimensions
  unsigned int N, D, L;

  // Input variables
  float *eps;
  float *X, *I, *W, *Gobs, *gradient_bias_in;
  unsigned int *Y;

  // Get input data from MATLAB
  X = (float *)mxGetData(prhs[0]);
  Y = (unsigned int*)mxGetData(prhs[1]);
  W = (float *)mxGetData(prhs[2]); //W: L*(D+1)
  Gobs = (float *)mxGetData(prhs[3]); //Gobs: D*L
  gradient_bias_in = (float *)mxGetData(prhs[4]); //gradient_bias_in:  L
  eps = (float *)mxGetData(prhs[5]);

  if (mxIsEmpty(prhs[6]))
    I = NULL;
  else
    I = (float *)mxGetData(prhs[6]);

  // Get matrix dimensions
  D = mxGetM(prhs[0]); //Number of features, height of X
  N = mxGetN(prhs[0]); //Number of instances, width of X
  L = mxGetM(prhs[2]); //Number of labels, height of W

  // Output variables
  float *G, *V, *gradient_bias;

  //Put V in the output
  mwSize Vdims[] = {1,1};
  plhs[0] = mxCreateNumericArray(2, Vdims, mxSINGLE_CLASS, mxREAL);
  V = (float *)mxGetData(plhs[0]);
  
  // Copy final G back to host
  mwSize Gdims[] = {D,L};
  plhs[1] = mxCreateNumericArray(2, Gdims, mxSINGLE_CLASS, mxREAL);
  G = (float *)mxGetData(plhs[1]);

  // Copy final gradient bias vector to host
  // Technically, gradient_bias will be Lx1 following the CUBLAS call, but
  // we want to return it to Matlab as a 1xL vector to agree with the (DxL)
  // gradient matrix.
  mwSize Bdims[] = {1,L};
  plhs[2] = mxCreateNumericArray(2, Bdims, mxSINGLE_CLASS, mxREAL);
  gradient_bias = (float *)mxGetData(plhs[2]);



  cuda_value(X, Y, W, Gobs, gradient_bias_in, N, D, L, *eps, I, G, V,
             gradient_bias);

}
#endif



/*

  Copyright 2009,2010 Augustus Lidaka, Shitanshu Aggarwal, Jerod Weinman

  This file is part of CUDAMAXENT.

    CUDAMAXENT is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    CUDAMAXENT is distributed in the hope that it will be useful, but
    WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with CUDAMAXENT.  If not, see <http://www.gnu.org/licenses/>.

*/
