/*
 * CresL.c
 *
 *  Created on: 2009/03/09
 *      Author: H.Shouno
 *  This is a wrapper program for C-cell, using CresL() 
 *  in Neocog.c by Prof.Fukushima.
 *
 */

#include <stdio.h>
#include <string.h>
#include <sys/stat.h>
#include "mex.h"
#include "neocog_proto.h"

void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] );

/*
 * Input & Output Arguments
 */

#define CAT_OUT		lhs[0] // output of category
#define UC_OUT		lhs[1] // outputu of layer

#define US_IN		rhs[0] // previous Us layer.
#define CAT_MAP_IN	rhs[1] // category map obtained by ReadCategory.c

//
// available expressions...
// Cat = CresL( US, CAT_MAP );
// [Cat Uc] = CresL( US, CAT_MAP );
//

#define NDIMS 3
//
// Gateway function
//
void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] )
{
	int i;
	double ***us, ***uc;
	int  nks, ns, nkc;
	int *cat_map, *catp;
	int nd, cat, catmin, catmax;
	double *coefp;
	int ndims;
	int dims[NDIMS];
	const int *usdims;
	
	if( nlhs < 1 )		mexErrMsgTxt( "CresL() requires 1 or 2 output(s): [Cat Uc] = Cres(...)" );
	if( nrhs < 2 ) 		mexErrMsgTxt( "Cres() requires 2 inputs: CresL( US, CAT_MAP )" );
	
	// 1st Us planes.
	ndims = mxGetNumberOfDimensions( US_IN );
	if( ndims != 3 )	mexErrMsgTxt( "US layer must be a 3 dimensional array sized as NS x NS x NKS" );
	usdims = mxGetDimensions( US_IN ); // Do not rewrite usdims
	ns = usdims[0];	nks = usdims[2];
	if( ns != usdims[1] || !mxIsDouble( US_IN ) ){
		mexErrMsgTxt( "US should be a 3 dimensional array like a piled double valued square matrix." );
	}
	if( ns <= 0 || nks <= 0 ){
		mexErrMsgTxt( "Us size must be a positive value." );
	}
	us = geta3D( nks, ns, ns, nks );
	coefp = mxGetPr( US_IN );
	memmove( &us[0][0][0], coefp, sizeof(double)*ns*ns*nks );

	// 2nd CATEGORY map US[LH-1] --> UC[LH-1]
	nd = mxGetM( CAT_MAP_IN );
	if( nd != nks ){
		mexErrMsgTxt( "Category size is mismatched. it should be a int32 valued matrix" );
	}
	if( !mxIsInt32( CAT_MAP_IN ) ){
		mexErrMsgTxt( "Category map must be a int32 valued matrix." );
	}
	catp = (int*)mxGetPr( CAT_MAP_IN );
	cat_map = get1Di( nd );
	memmove( cat_map, catp, nd*sizeof(int) );
	catmin = catmax = catp[0];
	for( i = 1 ; i < nd ; i++ ){
		if( catmin > catp[i] )	catmin = catp[i];
		if( catmax < catp[i] )	catmax = catp[i];
	}
	if( catmin != 0 ){
		mexErrMsgTxt( "Minimum of category number is not zero." );
	}
	if( catmax < 0 ){
		mexErrMsgTxt( "Maximum of category number is negative." );
	}
	nkc = catmax + 1;

	uc = geta3D( nkc, 1, 1, nkc ); 	
	cat = CresL(uc, us, cat_map, ns, nks, nkc);

	// 1st output: category
	CAT_OUT = mxCreateNumericMatrix( 1, 1, mxINT32_CLASS, mxREAL );
	catp = (int*)mxGetPr( CAT_OUT );
	*catp = cat;

	// 2nd output: Uc array
	if( nlhs > 1 ){
		ndims = NDIMS;
		dims[0] = dims[1] = 1;
		dims[2] = nkc;		// # of Uc planes 
		UC_OUT = mxCreateNumericArray( ndims, dims, mxDOUBLE_CLASS, mxREAL );
		coefp = mxGetPr( UC_OUT );

		memmove( coefp, &uc[0][0][0], 1*1*nkc*sizeof(double) );
	}
	return;
}
