/*
 * Sres.c
 *
 *  Created on: 2009/03/04
 *      Author: H.Shouno
 *  This is a wrapper program for S-cell, using Sres() 
 *  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 US_OUT		lhs[0] // output US layer cell planes. 

#define USSIZE_IN	rhs[0] // US plane size
#define UC_IN		rhs[1] // previous Uc layers.
#define A_IN		rhs[2] // excitatory connections (cell matrix (size nKs x nKc), and each element is a squared matrix)
#define B_IN		rhs[3] // inhibitory connections (double valued array ( size nKs ) )
#define C_IN		rhs[4] // fixed connections to inhibitory plane called USV
#define TH_IN		rhs[5] // Threshold value (optional)

//
// available expressions...
// Us = Sres( USSIZE, UC, A, B, C, [Threshold] );
//

/*			layer:	   U0 UG S1 C1 S2 C2 S3 C3 S4 C4
	int		Nc[]	= {   71,   37,   21,   13,    1};	// size of Uc(or U0) cell-plane, Nc[L+1]
	int		Ns[]	= {      68,   38,   22,    5   };	// size of Us cell-planes,  Ns[L]
	int		Ks[]	= {      16,   80,   80,   30   };	// no. of S-cell-planes, init.val., changeable,  Ks[L]
	int		Ksmax[]	= {      16,  200,  500,  400   };	// max. possible no. of S-cell-planes,  Ksmax[L]

	double	Ag		=     3.3;							// radius of connect. ag (off-surround)
	double	Ago		=     1.5;							// radius of on-center of ag
	double	gamcg	=  0.05;								// make_cg(cg): spatial decay at x=Ag
*/

#define KSMAX 512
#define NMAX 512
#define DEFAULT_THETA 0.5
#define NDIMS 3
//
// Gateway function
//
void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] )
{
	int i, j;
	double ***us, ***uc;
	int  nks, nkc, ns, nc;
	double ****as, bs[KSMAX], **cs;
	int na = -1;
	double theta;
	double *coefp;
	int ndims;
	int dims[NDIMS];
	const int *ucdims;

	if( nlhs != 1 )		mexErrMsgTxt( "Sres() requires 1 output: Us = Sres(...)" );
	if( nrhs < 5 ) 		mexErrMsgTxt( "Sres() requires 5 or more inputs: Sres( UsPlaneSize, UC, A, B, C, [Thteta] )" );
	

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

	// 2nd Uc planes.
	ndims = mxGetNumberOfDimensions( UC_IN );
	if( ndims != 3 )	mexErrMsgTxt( "UC layer must be a 3 dimensional array sized as NC x NC x NKC" );
	ucdims = mxGetDimensions( UC_IN ); // Do not rewrite ucdims
	nc = ucdims[0];	nkc = ucdims[2];
	if( nc != ucdims[1] || !mxIsDouble( UC_IN ) ){
		mexErrMsgTxt( "UC should be a 3 dimensional array like a piled double valued square matrix." );
	}
	if( nc <= 0 || nkc <= 0 ){
		mexErrMsgTxt( "Uc size must be a positive value." );
	}
	uc = geta3D( nkc, nc, nc, nkc );
	coefp = mxGetPr( UC_IN );
	
	memmove( &uc[0][0][0], coefp, sizeof(double)*nc*nc*nkc );

	// 3rd Coefficient A
	if( !mxIsCell(A_IN) ){
		mexErrMsgTxt( "Coef A must be a Cell Matrix sized (nKs x nKc)." );
	}
	nks = mxGetM( A_IN );
	if( nkc != mxGetN( A_IN ) ){
		mexPrintf( "Uc is (%d x %d x '%d'), but A is (%d x '%d'). Oops, size mismatch!", 
			   nc, nc, nkc, nks, mxGetM( A_IN ) );
		mexErrMsgTxt( "Coef size and Uc size are not identical." );
	}
	as = geta4D1( nks );
	add4D2( as, nks, nkc, nks );
	for( i = 0 ; i < nks ; i++ ){
		for( j = 0 ; j < nkc ; j++ ){
			mxArray *p = mxGetCell( A_IN, i + j * nks );
			if( p != NULL ){
				if( !mxIsDouble( p ) ){
					mexErrMsgTxt( "Coefficeint Matrix A should be double valued array." );
				}
				if( na < 0 )	na = mxGetN( p );
				else if( na != mxGetM( p ) ){
					mexErrMsgTxt( "some Coefficient matrix in A is not square." );
				}
				coefp = mxGetPr( p );
				as[i][j] = get2D( na, na );
				memmove( &as[i][j][0][0], coefp, na*na*sizeof(double) );
			}
			else{
				as[i][j] = NULL;
			}
		}
	}


	
	// 4th Coefficient B
	if( mxGetM( B_IN ) != nks || mxGetN( B_IN ) != 1 ){
		mexErrMsgTxt( "Coefficients B: size mismatched" );
	}
	if( !mxIsDouble( B_IN ) ){
		mexErrMsgTxt( "Coefficients B should be a double valued array" );
	}	
	coefp = mxGetPr( B_IN );
	memmove( bs, coefp, nks*sizeof(double) );

	// 5th Coef C  Uc --> Usv
	if( na != mxGetN( C_IN ) || na != mxGetM( C_IN ) ){
		mexErrMsgTxt( "Coefficeints C size is mismatched. it should be identical to the element of A" );
	}
	if( !mxIsDouble( C_IN ) ){
		mexErrMsgTxt( "Coefficients C must be a double valued matrix." );
	}
	coefp = mxGetPr( C_IN );
	cs = get2D( na, na );
	memmove( &cs[0][0], coefp, na*na*sizeof(double) );

	// 6th Threshold value (optional)
	theta = DEFAULT_THETA;
	if( nrhs >= 6 ){
		if( !mxIsDouble( TH_IN ) ){
			mexErrMsgTxt( "Threshold value should be a double value" );
		}
		theta = mxGetScalar( TH_IN );
	}

	ndims = NDIMS;
	dims[0] = dims[1] = ns;
	dims[2] = nks;
	US_OUT = mxCreateNumericArray( ndims, dims, mxDOUBLE_CLASS, mxREAL );
	coefp = mxGetPr( US_OUT );

	us = geta3D( nks, ns, ns, nks ); 
	Sres( us, uc, as, cs, bs, nc, nkc, ns, nks, na, nc, ns, theta ); // The last nc and ns parameter is for working memory

	memmove( coefp, &us[0][0][0], ns*ns*nks*sizeof(double) );
	return;
}
