/*
 * Cres.c
 *
 *  Created on: 2009/03/07
 *      Author: H.Shouno
 *  This is a wrapper program for C-cell, using Cres() 
 *  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 UC_OUT		lhs[0] // output UC layer cell planes. 

#define UCSIZE_IN	rhs[0] // UC plane size
#define US_IN		rhs[1] // previous Us layer.
#define D_IN		rhs[2] // excitatory connections (double valued squared matrix)

//
// available expressions...
// Uc = Cres( UCSIZE, US, D, [SubSampling rate] );
//

#define NDIMS 3
#define MC_DEFAULT 2	// sub-sampling rate (thinning-out)
//
// Gateway function
//
void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] )
{
	double ***us, ***uc;
	int  nks, ns, nc;
	double **d;
	int nd;
	int mc = MC_DEFAULT;
	double *coefp;
	int ndims;
	int dims[NDIMS];
	const int *usdims;
	
	if( nlhs != 1 )		mexErrMsgTxt( "Cres() requires 1 output: Uc = Cres(...)" );
	if( nrhs < 3 ) 		mexErrMsgTxt( "Cres() requires 3 inputs: Cres( UcPlaneSize, US, D )" );
	

	// 1st Us plane size
	if( !mxIsDouble( UCSIZE_IN ) ){
		mexErrMsgTxt( "1st argumtent is a Us plane size, which should be a double value " );
	}
	nc = (int)mxGetScalar( UCSIZE_IN );
	if( nc < 0 ){
		mexErrMsgTxt( "Us plane size should be a positive integer." );
	}

	// 2nd 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 );

	// 3rd Coefficient D   US --> UC
	nd = mxGetN( D_IN );
	if( nd != mxGetM( D_IN ) ){
		mexErrMsgTxt( "Coefficeints D size is mismatched. it should be a double valued swuare matrix" );
	}
	if( !mxIsDouble( D_IN ) ){
		mexErrMsgTxt( "Coefficients D must be a double valued matrix." );
	}
	coefp = mxGetPr( D_IN );
	d = get2D( nd, nd );
	memmove( &d[0][0], coefp, nd*nd*sizeof(double) );

	ndims = NDIMS;
	dims[0] = dims[1] = nc;
	dims[2] = nks;		// # of Uc planes 
	UC_OUT = mxCreateNumericArray( ndims, dims, mxDOUBLE_CLASS, mxREAL );
	coefp = mxGetPr( UC_OUT );

	uc = geta3D( nks, nc, nc, nks ); 
	Cres( uc, us, d, ns, nc, mc, nks, nd, ns );  // The last arguments is max( ns, nd ) for work memory
	memmove( coefp, &uc[0][0][0], nc*nc*nks*sizeof(double) );
	return;
}
