/*
 * Delta1stStage.c
 *
 *  Created on: 2009/02/20
 *      Author: H.Shouno
 *  This is a wrapper program for 1st layer training, using delta() 
 *  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 A_OUT	lhs[0] // output connection A (cell matrix)
#define B_OUT	lhs[1] // output connection B (array)

#define CS_IN		rhs[0] // fixed V-cell connection from Uc[0].
#define US_NUMS_IN	rhs[1] // extracting features in directions (should be a even number)
#define CG_IN		rhs[2] // on-center off-surround connection(optional)
#define Q_IN		rhs[3] // Deteminig of the Learing speed (optional)


//
// available expressions...
// [A B] = Delta1stStage( CS, US_NUMS, [CG], [Q] );
//
//

/*			layer:	   U0 UG S1 C1 S2 C2 S3 C3 S4 C4
	int		N0      =  65;								// size of input layer U0
	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]
	int		Kc[]	= {    2,   16,    0,    0,   10};	// no. of C-cell-planes, init.val., changeable,  Kc[L+1]

	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 MAXDIMS 16
#define KCMAX 512
#define QS_DEFAULT 100.0
// I/O PLANE parameters
#define NUC_DEFAULT 71
#define KUC_DEFAULT 2    // UC[0] should be a 71 x 71 x 2 planes
#define N0_DEFAULT 65
#define K0_DEFAULT 1    // U_INPUT is 65x65x1
#define NCG_DEFAULT 

// COEFFICIENTS parameters
#define CG_AG_DEFAULT	3.3
#define CG_AGO_DEFAULT	1.5
#define CG_GAMMA_DEFAULT 0.05

//
// Gateway function
//
void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] )
{
	int i, j, k;
	int na, nks;
	double ***u0, ***uc[1];
	int  k0, nu0, nkc0, nuc0;
	double **cg, **cs;
	int ng;
	double q;
	double ****as[1], *bs[1];
	double *coefp;

	if( nlhs != 2 )		mexErrMsgTxt( "Delta1stStage() requires 2 outputs: [A B] = Delta1stStage(...)" );
	if( nrhs < 2 ) 		mexErrMsgTxt( "Delta1stStage() requires 2 or 3 inputs: Delta( CS, US_NUMS, [CG] [Q] )" );
	

	// 1st fixed connection from Uc to Usv ---> size infomation is embedded
	na = (int)mxGetM( CS_IN );
	if( na != (int)mxGetN( CS_IN ) || !mxIsDouble( CS_IN ) ){
		mexErrMsgTxt( "1st arguments must be a double valued square matrix" );
	}
	if( na == 0 )		mexErrMsgTxt( "size of Coefficient matrix must be a positive." );
	cs = get2D( na, na );
	coefp = mxGetPr( CS_IN );
	memmove( cs[0], coefp, na*na*sizeof(double) );

	// 2nd input US
	nks = (int)mxGetScalar( US_NUMS_IN );
	if( nks < 0 )		mexErrMsgTxt( "# of Us[1] layers must be a positive." );	

	// 3rd On-Off conection
	if( nrhs >= 3 ){
		ng = (int)mxGetM( CG_IN );
		if( ng != (int)mxGetN( CG_IN ) || mxIsDouble(CG_IN) ){
			mexErrMsgTxt( "3rd arguments must be a double valued square matrix" );
		}
		if( ng == 0 )		mexErrMsgTxt( "size of Coefficient matrix must be a positive." );
		coefp = mxGetPr( CG_IN );
		cg = get2D( ng, ng );
		memmove( cg[0], coefp, ng*ng*sizeof(double) );
	}
	else{
		// fix coefficents
		ng = ((int) CG_AG_DEFAULT) * 2 + 1;	// size of on-center connect. cg[Ng][Ng]
		cg = get2D( ng, ng );
		make_cg( cg, ng, CG_AG_DEFAULT, CG_AGO_DEFAULT, CG_GAMMA_DEFAULT );
	}

	// 4th 
	if( nrhs >= 4 ){
		q = mxGetScalar( Q_IN );
		if( q <= 0.0 )	mexErrMsgTxt( "accelaration factor q must be a postive." );
	}
	else{
		q = QS_DEFAULT;
	}
	
	// calculating planes
	nu0 = N0_DEFAULT;	// # of Input Plane
	k0 = K0_DEFAULT;	
	u0 = geta3D( k0, nu0, nu0, k0 );

	nkc0 = KUC_DEFAULT;	// # of Uc(on/off) PLANES
	nuc0 = NUC_DEFAULT;	// Uc PLANE size
	uc[0] = geta3D( nkc0, nuc0, nuc0, nkc0 );

	as[0] = geta4D1( nks );
	bs[0] = get1D( nks );
	add4D2( as[0], nks, nkc0, nks );
	for( k = 0 ; k < nks ; k++ ){
		angle( u0[0], nu0, nks, k, 180.0, 0.5, -0.5 );
		onoff( uc[0], u0[0], cg, nu0, nuc0, ng );
		delta(as[0], bs[0], uc[0], cs, na, nuc0, nkc0, nks, k, nks/2, nks/2, KSMAX, q);
	}

	A_OUT = mxCreateCellMatrix( nks, nkc0 );
	B_OUT = mxCreateDoubleMatrix( nks, 1, mxREAL );
	for( i = 0 ; i < nks ; i++ ){
		for( j = 0 ; j < nkc0 ; j++ ){
			double **p = as[0][i][j];

			if( p != NULL ){
				int index = i + j * nks;
				mxArray *q = mxCreateDoubleMatrix( na, na, mxREAL );
				mxSetCell( A_OUT, index, q );
				coefp = (double*)mxGetPr( q );
				memmove( coefp, p[0], na*na*sizeof(double) );
			}
		}
	}
	coefp = (double*)mxGetPr( B_OUT );
	memmove( coefp, bs[0], nks * sizeof(double) );

	return;
}
