/*
 * PatLoad.c
 *
 *  Created on: 2009/02/01
 *      Author: H.Shouno
 *  This is a wrapper program for pat() 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 U_OUT		lhs[0] // output matrix
#define FNAME_IN	rhs[0] // file name
#define USIZE_IN	rhs[1] // Pattern size 64 or 65
#define TSTART_IN	rhs[2] // start of read out
#define TEND_IN		rhs[3] // end of read out
#define SAT_IN		rhs[4] // satruation factor


#define MAXSTRING 65536
#define UDEFAULT 65
#define NDIMS 3
#define TTMAX 400

static const int Lpat = 2052; // unit record length for ETL1

double* PatRead( int psize, double sat, int start, int end, char *fname );

// Gateway function
void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] )
{
	int i, size, start, end;
	int T_max = TTMAX;
	int dims[NDIMS];
	double sat;
	char fname[MAXSTRING];
	double udat[UDEFAULT*UDEFAULT], *umat[UDEFAULT];
	struct stat st;
	double *datap;


	umat[0] = &udat[0];
	for( i = 1 ; i < UDEFAULT ; i++ ){
		umat[i] = umat[i-1] + UDEFAULT;
	}

	if( nrhs < 1 ){
		mexErrMsgTxt( "PatTest requires 4 input arguments: filename, [pattern size = 64/65], [start number], and [end number]" );
	}
	if( nlhs != 1 ){
		mexErrMsgTxt( "PatTest requires 1 output argument: pattern matrix" );
	}

	if( mxIsChar(FNAME_IN) != 1 ){
		mexErrMsgTxt( "File name must be a string." );
	}
	mxGetString( FNAME_IN, fname, MAXSTRING-1 );
	if( stat( fname, &st ) != 0 ){
		mexPrintf( "File '%s' does not exist. ", fname );
		mexErrMsgTxt( "stat() was failed." );
	}
	T_max = st.st_size / Lpat; // max pattern numbers for ETL1 data
	mexPrintf( "%s contains %d patterns\n", fname, T_max );

	if( nrhs >= 2 )		size = (int)mxGetScalar( USIZE_IN );
	else			size = UDEFAULT;
	if( size != 64 && size != 65 ){
		mexErrMsgTxt( "Pattern size should be 64 or 65. This is specification of Fukushima's Neocognitron." );
	}
	if( size == 64 ){
		mexWarnMsgTxt( "Loading Pattern may be something wrong. I recommend to use this parameter as 65" );
	}
	if( nrhs >= 3 )		start = (int)mxGetScalar( TSTART_IN ); // How can I validate the argument?
	else			start = 0;
	if( start < 0 || start > T_max ){
		mexErrMsgTxt( "Start number is wrong" );
	}
	if( nrhs >= 4 )		end = (int)mxGetScalar( TEND_IN );  // How can I validate the argument?
	else			end = T_max;
	if( end < 0 || end > T_max || start > end ){
		mexErrMsgTxt( "End/Start number is wrong" );
	}
	if( nrhs >= 5 )		sat = mxGetScalar( SAT_IN );
	else			sat = 0.5;

	// U_OUT = mxCreateDoubleMatrix( size, size, mxREAL );
	dims[0] = size;
	dims[1] = size;
	dims[2] = end-start;
	U_OUT = mxCreateNumericArray( NDIMS, dims, mxDOUBLE_CLASS, mxREAL );
	if( U_OUT == NULL ){
		mexErrMsgTxt( "mxCreateNumericArray() is failed." );
	}
    
	datap = (double*)mxGetPr( U_OUT );
	/*
	for( t = start ; t < end ; t++ ){
		pat( umat, size, sat, t, start, end, T_max, fname );
		// mxSetPr( U_OUT, udat );
		memmove( datap, udat, size*size*sizeof(double) );
		datap += size*size;
	}
	*/
	memmove( datap, PatRead( size, sat, start, end, fname ),
		 size*size*(end-start)*sizeof(double) );

	return;
}



double* PatRead( int psize, double sat, int start, int end, char *fname )
{
	int pnums;
	int isize, t, i, T_max;
	unsigned char *upat, *up;
	double *dpat, *dp, pmax;
	FILE* ifp;
	struct stat st;
	const int dataoff = 32;
	const int W = 64;
	const int H = 63;

	if( stat( fname, &st ) != 0 ){
		mexPrintf( "File '%s' does not exist. ", fname );
		mexErrMsgTxt( "stat() was failed." );
	}
	if( (ifp = fopen( fname, "rb" ) ) == NULL ){
		mexPrintf( "Cannot open %s", fname );
		mexErrMsgTxt( "File open failure" );
	}
	T_max = st.st_size / Lpat; // max pattern numbers for ETL1 data
	upat = mxCalloc( st.st_size, sizeof(unsigned char) );
	isize = fread( upat, Lpat, T_max, ifp );
	if( isize != T_max ){
		mexPrintf("read error in ETL1 (TL):   isize=%d\n", isize );
		fclose(ifp);
		mexErrMsgTxt( "Read Error Occured" );
	}
	fclose(ifp);

	pnums = end-start;
	dpat = mxCalloc( pnums*psize*psize, sizeof(double) );
	for( t = 0 ; t < pnums ; t++ ){
		up = &upat[(start+t)*Lpat] + dataoff;	
		dp = &dpat[t*psize*psize];
		pmax = 0.0;

		for( i = 0 ; i < (W/2)*H ; i++ ){
			int y = i / (W/2) + 1;
			int xl = 2 * ( i % (W/2) );
			int xr = xl+1;

			dp[y*psize+xl] = (double)( ((*up) >> 4) & 0x0f );
			dp[y*psize+xr] = (double)( (*up) & 0x0f );
			if( dp[y*psize+xl] > pmax )	pmax = dp[y*psize+xl];					
			if( dp[y*psize+xr] > pmax )	pmax = dp[y*psize+xr];
			up++;
		}
		if( pmax > 0.0 ){
			if( sat == 0.0 ){
				double norm = 1.0 / pmax;
				for( i = 0 ; i < psize*psize ; i++ ){
					dp[i] *= norm;
				}
			}
			else{
				double sat1 = sat + 1.0;
				double norm = sat * pmax;
				for( i = 0 ; i < psize*psize ; i++ ){
					dp[i] = sat1 * dp[i] / (dp[i] + norm);
				}
			}
		}
	}
	return dpat;
}
