/*
 * ReadConnection.c
 *
 *  Created on: 2009/03/04
 *      Author: H.Shouno
 *  This is a wrapper program for read part of the rwcon() 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 CELL matrix a{nKs, nKc}
#define B_OUT	lhs[1] // output B coeeficients b[Ks]
#define FNAME_IN	rhs[0] // file name
#define LAYER_IN	rhs[1] // Layer number.

#define MAXSTRING 65536
#define KSMAX 512 // upper limit of Us Planes
#define LMAX 32   // upper limit of stages.
#define LH   4    // recongition layer (assumed in Neocog.c) : Todo 

//
// Size information is lacked in  'connect.dat'. Thus we should assume size for reading connection.
// Todo: size information should be embedded in data format
//
static	double	As[]	= {     3.4,  3.4,  3.4,  4.9   };	// radius of connect. as, As[L]
//
//                            UG S1 C1 S2 C2 S3 C3 S4 C4
//
static	int	Nc[]	= {   71,   37,   21,   13,    1};	// size of Uc(or U0) cell-plane, Nc[L+1]
static	int	Ns[]	= {      68,   38,   22,    5   };	// size of Us cell-planes,  Ns[L]


// Gateway function
void mexFunction( int nlhs, mxArray *lhs[], int nrhs, const mxArray *rhs[] )
{
	char fname[MAXSTRING];
	double ****a[LH];
	double *b[LH], *coefp;
	int Mas[LH], *cat;
	int Na[LMAX];           // size of coefficients a[]
	int Kc[LMAX], Ks[LMAX]; // number of Uc[], and Us[] planes.
	int i, j, l, ll;
	FILE *ifp; // may be dummy, since fopen() and fclose() pair procedure is carried out.
	struct stat st;

	if( nlhs != 2 )		mexErrMsgTxt( "ReadConnect() requires 2 outputs: [A B] = ReadConnect( fname, layer_num )" );
	if( nrhs != 2 )		mexErrMsgTxt( "ReadConnect() requires 2 inputs: [A B] = ReadConnect( fname, layer_num )" );

	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." );
	}
	l = (int)mxGetScalar( LAYER_IN );
	if( l < 0 || l >= LH )	mexErrMsgTxt( "ReadConnect(): layer number 0 <= l < LH (recognition layer)" );
	ll = -l-2;
	
	for( i = 0 ; i <= l ; i++ ){
		// from Neocog.c
		if ((Nc[i] + Ns[i]) % 2 == 0) {					// alligned (Na is odd)
			Na[i] = ((int) As[i]) * 2 + 1;
		} else {							// staggered (Na is even)
			Na[i] = ((int) (As[i] + 0.5)) * 2;
			// Exactly, Na should be: (na-0.5)^2 + 0.5^2 = A^2;  na = 0.5 + sqrt(A^2 - 0.5^2);  Na = na * 2;
		}
	}

	rwcon( ifp, fname, a, b, &cat, Na, Kc, Ks, Mas, LH, ll );

	A_OUT = mxCreateCellMatrix( Ks[l], Kc[l] );
	for( i = 0 ; i < Ks[l] ; i++ ){
		for( j = 0 ; j < Kc[l] ; j++ ){
			double **p = a[l][i][j];
			if( p != NULL ){
				int index = i + j * Ks[l];
				mxArray *q = mxCreateDoubleMatrix( Na[l], Na[l], mxREAL );

				mxSetCell( A_OUT, index, q );
				coefp = (double*)mxGetPr( q );
				memmove( coefp, p[0], Na[l]*Na[l]*sizeof(double) );
			}
		}
	}

	B_OUT = mxCreateDoubleMatrix( Ks[l], 1, mxREAL );
	coefp = (double*)mxGetPr( B_OUT );
	memmove( coefp, b[l], Ks[l] * sizeof(double) );

	return;
}
