/*
 * OnOff.c
 *
 *  Created on: 2009/02/10
 *      Author: H.Shouno
 *  This is a wrapper program for onoff() 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 on/off matricies
#define CG_OUT		lhs[1] // output coefficient if you need. (optional)
#define US_IN		rhs[0] // Input from ETL1 pattern
#define UCSIZE_IN	rhs[1] // ouput uc-plane size
#define CG_IN		rhs[2] // filter coefficient if you set. (optional)

//
// available expressions...
// Uc = OnOff( Us );
// Uc = OnOff( Us, Ucsize );
// Uc = OnOff( Us, Ucsize, CG );
// [Uc CGout] = OnOff( Us );
// etc... 
//
//


#define UC_DEFAULT_SIZE 71
#define OUT_PLANES 2 // on and off planes
#define UCDIMS 3

#define AG 3.3 // from Neocog.c
#define AGO 1.5
#define GamCG 0.05
#define CG_DEFAULT_SIZE ((int)(2*(AG)+1))

//
// Gateway function
//
void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] )
{
	int n0, nc0, ng, i, k, n, m;
	int ucdims[UCDIMS];
	double *inp, *outp, *coefp;
	double **cg, **uc[2], **us;
	double default_cg[CG_DEFAULT_SIZE*CG_DEFAULT_SIZE];
	double *default_cg_mat[CG_DEFAULT_SIZE];

	if( nrhs < 1 )		mexErrMsgTxt( "OnOff() requires input parameter(s): OnOff( Us, [UcSize], [CoefMatrix])" );

	if( nlhs < 1 )		mexErrMsgTxt( "OnOff() requires ouput parameter(s): [Uc [CoefMatrix] ] = OnOff( Us, ... )" );

	// right side
	m = mxGetM( US_IN );	n = mxGetN( US_IN );	
	if( !mxIsDouble( US_IN ) || n != m ){
		mexErrMsgTxt( "Input must be a double valued square matrix" );
	}
	n0 = n;
	if( n0 < 0 )		mexErrMsgTxt( "Us size is illegal." );
	inp = (double*)mxGetPr( US_IN );
	us = (double**)mxCalloc( sizeof(double*), n0 );
	us[0] = inp;
	for( n = 1 ; n < n0 ; n++ )	us[n] = us[n-1] + n0;

	if( nrhs >= 2 )		nc0 = mxGetScalar( UCSIZE_IN );
	else			nc0 = UC_DEFAULT_SIZE;
	if( nc0 < 0 )		mexErrMsgTxt( "Uc size is illegal." );

	if( nrhs >= 3 ){ // Coefficients are assigned by CG_IN
		int cgn, cgm;
		cgn = mxGetN( CG_IN );	cgm = mxGetM( CG_IN );
		if( !mxIsDouble( CG_IN ) || cgn != cgm ){
			mexErrMsgTxt( "Coefficient Matrix must be a double valued square matrix" );
		}
		ng = cgn;
		if( ng < 0 )		mexErrMsgTxt( "CG size is illegal." );
		cg = mxCalloc( sizeof(double*), ng );
		cg[0] = (double*)mxGetPr( CG_IN );
		for( i = 1 ; i < cgn ; i++ ){
			cg[i] = cg[i-1] + cgm;
		}
	}
	else{  	// No 3rd argument. use default coefficient set
		ng = CG_DEFAULT_SIZE;	// size of on-center connect. cg[ng][ng]; staggered version is not implemented.
		default_cg_mat[0] = &default_cg[0];
		for( i = 1 ; i < ng ; i++ ){
			default_cg_mat[i] = default_cg_mat[i-1] + ng;
		}
		cg = default_cg_mat;
		make_cg( cg, ng, AG, AGO, GamCG );
	}

	ucdims[0] = ucdims[1] = nc0;
	ucdims[2] = OUT_PLANES;
	UC_OUT = mxCreateNumericArray( UCDIMS, ucdims, mxDOUBLE_CLASS, mxREAL );
	outp = (double*)mxGetPr( UC_OUT );
	for( k = 0 ; k < OUT_PLANES ; k++ ){
		uc[k] = mxCalloc( sizeof(double*), nc0 );
		uc[k][0] = &outp[k*nc0*nc0];
		for( n = 1 ; n < nc0 ; n++ ){
			uc[k][n] = uc[k][n-1] + nc0;
		}
	}
	if( nlhs >= 2 ){
		CG_OUT = mxCreateDoubleMatrix( ng, ng, mxREAL );
		coefp = (double*)mxGetPr( CG_OUT );
		memmove( coefp, cg[0], sizeof(double)*ng*ng );
	}
    
	onoff( uc, us, cg, n0, nc0, ng );

	return;
}
