
		////////////////////////////////////////////////////////
		/////  Neocognitron               by K. Fukushima  /////
		/////           4 staged network                   /////
		/////   Edge extractors only (plane 0--15) in US1  /////
		/////   Inhibitory surround in connection ac       /////
		/////          New thinning-out process            /////
		////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#ifdef NEOCOG_MATLAB	// by H.Shouno 2009/02/20
#include "mex.h"	// to use mx-Functions. Setup MEX environment appropriately by MATLAB manual.
// initialyze routine may become a risk of memory leak in MATLAB application.
#define MALLOC(X)	mxMalloc( (X) )
#define PRINTF		mexPrintf
#define STATIC
#define EXIT		mxAssert(0, "Logical Error in Neocog.c")
#else
#define MALLOC(X)	malloc( (X) )
#define PRINTF		printf
#define	STATIC		static
#define EXIT		exit(1)
#endif // NEOCOG_MATLAB

////  Prototype
void main(void);
double ***getp(			// memory for 1-D pointer array **u[L]
	int		L
	);
double *get1D(			// get memory for 1-D array m[N]
	int		N
	);
int *get1Di(			// get memory for 1-D int array m[N]
	int		N
	);
double **get2D(			// get memory for 2-D array u[N1][N2]
	int		N1,
	int		N2
	);
int **get2Di(			// get memory for 2-D int array m[N1][N2]
	int		N1,
	int		N2
	);
unsigned char **get2Dc(	// get memory for 2-D unsigned char array ch[N1][N2]
	int		N1,
	int		N2
	);
double ***geta3D(		// get addable memory for 3-D array u[K][N1][N2]
	int		K,			// initial array size
	int		N1,
	int		N2,
	int		Kmax		// maximum possible array size: u[Kmax][N1][N2]
	);
void add3D(				// add memory to 3-D array u[K][N1][N2]
	double	***u,		// (in) u = (double ***)malloc(sizeof(double **) * Kmax)
	int		Knew,		// new array size
	int		N1,
	int		N2,
	int		Kmax		// maximum possiboe array size
	);
double ****geta4D1(		// get addable memory for 1-D pointer array ***a[]
	int		Kmax		// maximum possible pointer array size: ***a[Kmax]
	);
void add4D2(			// add memory to 2-D pointer array **a[K][Kc]
	double	****a,		// (in) a = (double ****)malloc(sizeof(double *) * Kmax)
	int		Knew,		// new pointer array size: **a[Knew][Kc]
	int		Kc,
	int		Kmax		// maximum possible array size: **a[Kmax][Kc]
	);
void Sres(				// response of S-layer
	double	***us,		// (out) response of S-cells: us[Ks][Ns][Ns]
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	****as,		// (in) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	**cs,		// (in) exc. connections to V-cells: cs[Na][Na]
	double	*bs,		// (in) inh. connections to S-cells: bs[Ks]
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Ks,			// numberof cell-planes in S-layer
	int		Na,			// size of the connections (as & cs): (na*2) or (na*2+1)
	int		NcMAX,		// max of Nc (and Ng) among all stages
	int		NsMAX,		// max of Ns among all stages
	double	theta		// threshold of an S-cell
	);
void Cres(				// response of C-layer
	double	***uc,		// (out) response of C-cells: uc[K][Nc][Nc]
	double	***us,		// (in) response of S-cells: us[K][Ns][Ns]
	double	**ac,		// (in) exc. connections to C-cells: ac[Nd][Nd]
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		mc,			// ratio of thinning-out
	int		K,			// number of cell-planes in S- and C-layers (K=Ks=Kc)
	int		Nd,			// size of the connections (ac): (nd*2) or (nd*2+1)
	int		NsMAX		// max of Ns among all stages
	);
int CresL(				// response of C-layer without lateral inhibition
	double	***uc,		// (out) response of C-cells: uc[Kc][1][1]
	double	***us,		// (in) response of S-cells: us[K][Ns][Ns]
	int		*cat,		// (in) category of the S-plane: cat[Ks]
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Ks,			// number of cell-planes in S-layers (Ks)
	int		Kc			// number of cell-planes in C-layers (Kc)
	);
void onoff(				// response of on- & off-center cell layers
	double	***uc,		// (out) response of on- & off-center cells: uc[2][Nc][Nc]
	double	**us,		// (in) response of receptor cells: us[Ns][Ns]
	double	**cg,		// (in) connections for on-center: cg[Ng][Ng]
	int		Ns,			// size of the receptor layer: (ns*2) or (ns*2+1)
	int		Nc,			// size of the on- & off-center cell layer: (nc*2) or (nc*2+1)
	int		Ng			// size of the connections (ag): (ng*2) or (ng*2+1)
	);
void make_cg(			// make connections for on-center
	double	**cg,		// (out) Connections: cg[Ng][Ng]
	int		Ng,			// size of connection cg: (ng*2) or (ng*2+1)
	double	Ag,			// radius of the off-periphery (connect. area)
	double	Ago,		// radius of the on-center component
	double	gamcg		// spatial decay at x=Ag
	);
void make_ac(			// make center surround connections ac
	double	**ac,		// (out) Connections: ac[Na][Na]
	int		Na,			// size of the connection ac: (na*2) or (na*2+1)
	double	A,			// radius of the off-periphery (connect. area)
	double	Ao,			// radius of the on-center component
	double	gamma,		// spatial decay of exc. connect. at x=Ao
	double	gammao,		// spatial decay of inh. connect. at x=Ao
	double	gammap,		// spatial decay of inh. connect. at x=A
	double	co,			// value of ac[][] at the center
	double	pn			// ratio of negative to positibe components
	);
void make_c(			// make connections (spatial weight)
	double	**c,		// (out) Connections: c[Na][Na]
	int		Na,			// size of the connection A: (na*2) or (na*2+1)
	double	A,			// radius of the connectable area
	double	gamma,		// spatial decay at x=A
	double  co			// value of c[][] at the center
	);
void angle(				// draw an angle of ang degree
	double	**u,		// (out) generated pattern: u[N][N]
	int		N,			// size of the layer
	int		K,			// number of different orientations
	int		k,			// sequence number of the orientation of the edge
	double	ang,		// angle (in degree) of the bend
	double	x0,			// x offset of the center of the pattern
	double	y0			// y offset of the center of the pattern
	);
void edge(
	double	**u,		// (i/o) generated pattern: u[N][N]
	int		N,			// size of the layer
	double	xs,
	double	ys,			// edge from (xs,ys) to (xe,ye)
	double	xe,
	double	ye
	);
int Sorg(				// self-organization of S-layer
	double	****as,		// (out) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	*bs,		// (out) inh. connections to S-cells: bs[Ks], new Ks
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**cs,		// (in) spatial weight: cs[Na][Na]
	double	**asini,	// (in) exc. initial connections to S-cells: as[Na][Na]
	double	**cl,		// (in) connect. for lat.inh. from us to usini: cl[Nl][Nl]
	double	***us,		// (work) response of S-cells: us[Ks][Ns][Ns], old Ks
	int		Nc,			// size of the C-layer:  (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer:  (ns*2) or (ns*2+1)
	int		Ks,			// number of cell-planes in S-layer, old value
	int		Na,			// size of connections (as & cs):  (na*2) or (na*2+1)
	int		Nasini,		// size of connect. asini:  (nasini*2) or (nasini*2+1)
	int		Nl,			// size of connection cl (lat. inh. from us to usini):  Nl=nl*2+1
	double	D,			// radius of competition area
	double	qs,			// parameter determining the speed of learning
	double	theta,		// threshold of an S-cell
	int		Ksmax,		// maximum possible Ks
	int		KsMAX,		// max of KsMAX among all stages
	int		NcMAX,		// max of Nc (and Ng) among all stages
	int		NsMAX,		// max of Ns among all stages
	double	thrini		// threshold for cutting usini
	);
void delta(				// strengthen connections
	double	****as,		// (out) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	*bs,		// (out) inh. connections to S-cells: bs[Ks]
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**cs,		// (in) spatial weight: cs[Na][Na]
	int		Na,			// size of the connections (as & cs): (na*2) or (na*2+1)
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		Kc,			// numberof cell-planes in C-layer
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		kseed,		// cell-plane of the seed cell
	int		iseed,		// x-coordinate of the seed cell
	int		jseed,		// y-coordinate of the seed cell
	int		KsMAX,		// max of KsMAX among all stages, because KcMAX <= KsMAX
	double	qs			// parameter determining the speed of learning
	);
void seed(				// select seed cells
	int		*iseed,		// (out) x-coordinate of the k-th cell-plane: iseed[Ks]
	int		*jseed,		// (out) y-coordinate of the k-th cell-plane: jseed[Ks]
	double	***us,		// (in) response of S-cells: us[Ks][Ns][Ns]
	double	**wmax,		// (work) working area: wmax[Ns][Ns]
	int		**kwmax,	// (work) working area: kwmax[Ns][Ns]
	int		Ns,			// size of the S-layer
	int		Ks,			// number of cell-planes in S-layer
	int		NsMAX,		// max of Ns among all stages
	double	D			// radius of competition area
	);
void Sresini(			// response of seed-selecting plane
	double	**usini,	// (out) response of seed-selecting plane: usini[Ns][Ns]
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**asini,	// (in) exc. initial connections to S-cells: as[Na][Na]
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Na,			// size of the connections (asini, as & cs): (na*2) or (na*2+1)
	int		NcMAX,		// max of Nc (and Ng) among all stages
	double	thrini		// threshold for cutting usini
	);
void seedini(			// select seed cells
	int		*iseed0,	// (out) x-coordinate of the k-th cell-plane: iseed0[1]
	int		*jseed0,	// (out) y-coordinate of the k-th cell-plane: jseed0[1]
	double	**usini,	// (in) response of S-cells: usini[Ns][Ns]
	int		Ns			// size of the S-layer
	);
void lat(				// lateral inhibition from uin to uout
	double	**ut,		// (i/o) cell-plane to be inhibited: ut[N][N]
	double	**uf,		// (in) inhibiting cell-plane: uf[N][N]
	double	**cl,		// (in) connections for lateral inhibition: cl[Nl][Nl]
	int		N,			// size of the layer
	int		Nl			// size of connection cl:  Nl=nl*2+1
	);
int SorgL(				// self-organization of S-layer (new cell-plane if winner is not from correct category)
	double	****as,		// (out) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	*bs,		// (out) inh. connections to S-cells: bs[Ks], new Ks
	int		*cat,		// (out) category of the S-plane: cat[Ks]
	int		Cat,		// (in) category of the training pattern
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**cs,		// (in) spatial weight: cs[Na][Na]
	double	**asini,	// (in) exc. initial connections to S-cells: as[Na][Na]
	double	***us,		// (work) response of S-cells: us[Ks][Ns][Ns], old Ks
	int		Nc,			// size of the C-layer:  (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer:  (ns*2) or (ns*2+1)
	int		Ks,			// number of cell-planes in S-layer, old value
	int		Na,			// size of the connections (as & cs): (na*2) or (na*2+1)
	int		Nasini,		// size of the connect. asini:  (nasini*2) or (nasini*2+1)
	double	qs,			// parameter determining the speed of learning
	double	theta,		// threshold of an S-cell
	int		Ksmax,		// maximum possible Ks
	int		KsMAX,		// max of KsMAX among all stages, because KcMAX <= KsMAX
	int		NcMAX,		// max of Nc (and Ng) among all stages
	int		NsMAX,		// max of Ns among all stages
	FILE	*fpT		// (in) file pointer for debugging.  If ==NULL, no debugging
	);
void rwcon(				// read/write connections.  Also get memores as and bs, in the read mode
	FILE	*fpC,		// (in) file pointer
	char	*fnameC,	// (in) file name
	double	*****as,	// (i/o) exc. connect. to S-cells: as[L][Ks][Kc][Na][Na]
	double	**bs,		// (i/o) inh. connect. to S-cells: bs[L][Ks]
	int		**cat0,		// (i/o) category of the cell-plane UsL: cat0[1][KsL]=cat[Ks]
	int		*Na,		// (in) size of connect. as
	int		*Kc,		// (i/o) no. of C-cell-planes: Kc[L+1]
	int		*Ks,		// (i/o) no. of S-cell-planes: Ks[L]
	int		*Mas,		// (i/o) no. of active as[L][Ks][Kc]: Mas[L]
	int		L,			// number of stages
	int		learn		// switch:  if (learn>0) write, else read.  If (learn<0) up to (-learn-1)
	);
void printu(
	FILE	*fp,		// (in) file pointer
	double	***u,		// (in) layer to be printed: u[K][N][N]
	int		N,			// size of a cell-plane: (n*2+1) or (n*2)
	int		K,			// number of cell-planes
	char	title[],	// title to be printed
	int		l,			// stage number to be printed in the title. if (l<0) no print
	int		KMAX		// >0: maximum possible K, Initialize only, No printing.  ==0: print response
	);
void outH(				// hexadecinmal output of the response of a layer
	FILE	*fp,		// (in) file pointer
	double	***u,		// (in) layer to be printed: u[K][N][N]
	int		N,			// size of a cell-plane: (n*2) or (n*2+1)
	int		K,			// number of cell-planes
	char	title[],	// title to be printed
	int		l,			// stage number in the title. if (l==-1) no number; if (l==-2) no title
	int		KMAX		// >0: maximum possible K, Initialize only, No output.  ==0: write output
	);
void outHhead(			// write header of hexadecinmal output file fpN
	FILE	*fpN,		// (in) file pointer
	char	fnameNh[],	// header of fnameN (fpNh)
	int		Npages		// no. of pages. >0: all data in one ps file.  =0: separate eps files.
	);
void printuadd(			// print sum of k1-th to k2-th cell-planes
	FILE	*fp,		// (in) file pointer
	double	***u,		// (in) layer to be printed: u[K][N][N]
	int		N,			// size of a cell-plane: (n*2) or (n*2+1)
	char	title[],	// title to be printed
	int		l,			// stage number to be printed in the title. if (l<0) no print
	int		k1,			// added cell-planes from k1
	int		k2			// added cell-planes to k2
	);
void pat(				// get a pattern (from ETL1 file or memory)
	double	**u,		// (out) generated pattern: u[N][N]
	int		N,			// size of the layer:  N*N = 64*64 or 65*65
	double	sat,		// saturation factor
	int		t,			// sequence number of the pattern to be generated
	int		T0,			// starting no. of patterns to be read from file
	int		TT,			// total no. of patterns to be read from file
	int		TTmax,		// maximum possible TT
	char	fnameETL[]	// name of ETL1 file (fpETL)
	);
void score(				// make confusion matrix
	int		**mat,		// (i/o) confusion matrix: mat[KcL][KcL+1]
	int		recCat,		// result of recognition
	int		Cat,		// category of the input pattern
	int		KcL			// number of categories
	);
void printsc(			// print confusion matrix
	FILE	*fp,		// (in) file pointer
	int		**mat,		// (in) confusion matrix: mat[KcL][KcL+1]
	int		KcL			// number of categories
	);
void ErrOut(			// hexadecinal output of erroneously recognized patterns
	double	***u,		// (in) layer to be printed: u[1][N][N]
	int		N,			// size of a cell-plane: (n*2) or (n*2+1)
	char	fnameE[],	// file name (fpE), erroneously recognized patterns
	char	fnameEh[],	// header of fnameE (fpEh)
	int		iend		// if (iend==1), no more output (close file)
	);

// *************************************************************************

#ifndef NEOCOG_MATLAB

void main(void) {

  ////  Parameters

	int		L = 4;										// number of stages

		// 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		mc[]	= {          2,    2,    2,    1};	// Cres: step of C-cells,  mc[L]

	double	sat0	=  0.5;								// pat: saturation, input pat U0
	double	Ag		=     3.3;							// radius of connect. ag (off-surround)
	double	Ago		=     1.5;							// radius of on-center of ag
	double	As[]	= {     3.4,  3.4,  3.4,  4.9   };	// radius of connect. as, As[L]
	double	Aco[]	= {        3.4,  3.4,  4.4,  0.0};	// radius of on center in connect. ac, Ac[L]
	double	Ac[]	= {        9.4,  7.4, Aco[2],0.0};	// radius of connect. ac, Ac[L]
	double	D[]		= {     0.0,  3.1,  3.1,  0.0   };	// radius of competition area,  D[L]
	double	Asini[]	= {     0.0,As[1],As[2],As[3]   };	// make_c(asini): radius of connect.: Asini[L]
	double	Al[]	= {     0.0, D[1], D[2],  0.0   };	// radius of lat.inh. cl: Nl=((int)Al)*2+1, Al[L]

	double	gamcg	=  0.05;								// make_cg(cg): spatial decay at x=Ag
	double	gamS[]	= {     0.7,  0.7,  0.7,  0.6   };	// make_c(cs): spatial decay, gammaS,  gamS[L]
	double	coS[]	= {     1.0,  1.0,  1.0,  1.0   };	// make_c(cs): value at the center,  coS[L]
	double	gamSini[]={  0.0,gamS[1],gamS[2],gamS[3]};	// make_c(asini): spatial decay of asini,  gamSini[L]
	double	coSini[]= {   0.0001,0.0001,0.0001,0.0001};	// make_c(asini): value at the center of asini,  coSini[L]
	double	gamC[]	= {        0.7,  0.7,  0.7,  1.0};	// make_ac(ac): spatial decay of exc. center at x=Ao,  gamC[L]
	double	gamCo[]	= {        0.12, 0.12, 0.12, 0.0};	// make_ac(ac): spatial decay of inh. connect. at x=Ao
	double	gamCp[]	= {        0.8,  0.8,  0.8,  0.0};	// make_ac(ac): spatial decay of inh. connect. at x=A
	double	coC[]	= {        4.0,  4.0,  4.0,  1.0};	// make_ac(ac): value at the center,  coC[L]
	double	pnC[]	= {        1.0,  1.0,  1.0,  0.0};	// make_ac(ac): ratio of negative to positive components
	double	gaml[]	= {     0.0,  0.7,  0.7,  1.0   };	// make_c: spatial decay, lat.inh. cl,  gaml[L]
	double	col[]	= {     0.0,  1.0,  1.0,  1.0   };	// make_c: strength of lat.inh. cl,  col[L]

	double	thetal[]= {    0.00, 0.66, 0.67, 0.30   };	// Sres: threshold of S-cells (learning),  theta[L]
	double	thetar[]= {    0.55, 0.51, 0.58, 0.30   };	// Sres: threshold of S-cells (recognition),  theta[L]

	double	qs[]	= {   100.0,100.0,100.0,100.0   };	// speed of learning,  qs[L]
	double	thrini[]= {     0.0,  0.3,  0.3,  0.0   };	// Sresini: threshold for cutting usini,  thrini[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]

	int		learn = 1;						// stage from which learning start (if learn=0, test)
//	int		learn = 0;						// stage from which learning start (if learn=0, test)
	int		TL0 = 0;						// starting no. of learning patterns
	int		TL = 400;						// no. of learning pat for each stage
//	int		TL = 3000;						// no. of learning pat for each stage
//	int		TTL = 0;						// test for TTL learned patterns (if TTL=0, skip; if TTL=TL, all learned patterns)
	int		TTL = TL;						// test for TTL learned patterns (if TTL=0, skip; if TTL=TL, all learned patterns)
	int		TT0 = 0;						// starting no. of test patterns
	int		TT = 400;						// no. of tests for unlearned patterns
//	int		TT = 3000;						// no. of tests for unlearned patterns
	int		repeat[]= {       1,    1,    1,    9   };	// no. of max. repetition in the learning,  repeat[L]

	int		iprl[]	= { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};	// print control in the learning
//	int		iprl[]	= { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};	// print control in the learning
	int		iprt[]	= { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};	// print control in the test (==2: printuadd)
//	int		iprt[]	= { 1, 1, 3, 3, 3, 3, 3, 3, 1, 1};	// print control in the test (==2: printuadd)
	int		ifpN[]	= {0, 19};	// output to fpN, tt == ifpN[0] thru ifpN[1],   (ifpN[1]-ifpN[0]<100)
	int		ifpNall = 1;		// =1: all data in one ps file.  =0: separate eps files.

	char	fnamea[] = "res.dat";				// name of output file (fpa: response)
	char	fnameC[] = "connect.dat";			// name of i/o file (fpC: connections)
	char	fnameL[] = "PatLearn.dat";			// name of learning pattern file (fpL)
	char	fnameT[] = "PatTest.dat";			// name of test pattern file (fpT)
	char	fnameNhead[] = "Header.dat";		// header of fnameN and fnameNall (fpNh)
	char	fnameN[] = "out  .eps";				// name of eps (hexadecimal) output file (fpN: response)
	char	fnameNall[] = "outall.ps";			// file name (fpN), all responses together
	char	fnameEh[] = "ErrHead.dat";			// header of fnameE
	char	fnameE[] = "outErr.ps";				// file name, erroneously recognized patterns


   // memories for parameters
	int		Ng;			// size of on-center connect. cg[Ng][Ng]
	int		*Na;		// size of connect. as & cs:  Na=((int)(As+0.5))*2,  Na[L]
	int		*Nd;		// size of connect. ac:  Nd=((int)(Ac+0.5))*2,  Nd[L]
	int		*Nasini;	// size of connect. asini: Nasini=((int)(Asini+0.5))*2, Nasini[L]
	int		*Nl;		// size of lat.inh. cl:  Nl=((int)Al)*2+1,  Nl[L]
	int		*Mas;		// no. of active as[L][Ks][Kc]: Mas[L]

   // memories for fixed connections
	double	**cg;		// connect. of on-center: cg[Ng][Ng]
	double	***cs;		// exc. connect. to V-cells: cs[L][Na][Na]
	double	***ac;		// exc. connect. to C-cells: ac[L][Nd][Nd]
	double	***asini;	// exc. initial connect. to S-cells: asini[L][Nasini][Nasini]
	double	***cl;		// lat.inh. from us to usini: cl[L][Nl][Nl]

   // memories for responses -- fixed
	double	***u0;		// response of input layer: u0[N0][N0]
	int		**mat;		// confusion matrix: mat[KcL][KcL+1]

   // memories for responses -- expandable
	double	****us;		// response of S-cells: us[L][Ks][Ns][Ns]
	double	****uc;		// response of C-cells: uc[L][Kc][Nc][Nc]
	int		*cat;		// category of the cell-plane in UsL: cat[KsL]

   // memories for connections -- expandable
	double	*****as;	// exc. connect. to S-cells: as[L][Ks][Kc][Na][Na]
	double	**bs;		// inh. connect. to S-cells: bs[L][Ks]


   // working
	int		NsMAX;				// NsMAX=max(Ns[l]) for l=0,1,..,L-1
	int		KsMAX;				// KsMAX=max(Ksmax[l]) for l=0,1,..,L-1
	int		NcMAX;				// NcMAX=max(Nc[l]) for l=0,1,..,L
	FILE	*fpa = NULL;		// file pointer (response)
	FILE	*fpN = NULL;		// file pointer (response, hexadecimal)
	FILE	*fpC = NULL;		// file pointer (connections)
	int		l, ll, t, tt, k, Cat, recCat, i, j, Npages, learn1,
			Iprt, KsL, repeatL, TTmax;


   // open file "fnamea"
    if ((fpa = fopen(fnamea, "w")) == NULL) {
		PRINTF("file open error fpa\n");
		EXIT;
    }

	/// Check file-output condition
	Iprt = 0;
	for (l = 0; l < 2 + L * 2; l++) {
		if (iprt[l] > 0)  Iprt = 1;
	}									// if Iprt=0, no file-output at each time


  /// calculate connection sizes N depending on Ns and Nc
	// get array and set parameters

	Ng = ((int) Ag) * 2 + 1;	// size of on-center connect. cg[Ng][Ng]

	if (((N0 + Nc[0]) % 2) == 0) {		// alligned
		Ng = ((int) Ag) * 2 + 1;		// size of on-center connect. cg[Ng][Ng]
	} else {							// staggered
		Ng = ((int) (Ag + 0.5)) * 2;
	}

	Na = get1Di(L);				// size of connect. as & cs: Na[L]
	Nasini = get1Di(L);			// size of init. connect. asini: Nasini[L]
	for (l = 0; l < L; l++) {
		if ((Nc[l] + Ns[l]) % 2 == 0) {					// alligned (Na is odd)
			Na[l] = ((int) As[l]) * 2 + 1;
			Nasini[l] = ((int) Asini[l]) * 2 + 1;
		} else {										// staggered (Na is even)
			Na[l] = ((int) (As[l] + 0.5)) * 2;
			Nasini[l] = ((int) (Asini[l] + 0.5)) * 2;	// approxinate
			// 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;
		}
	}

	Nd = get1Di(L);				// size of connect. ac: Nd[L]
	for (l = 0; l < L; l++) {
		if ( ((mc[l] == 2) && (Ns[l] % 2 == 1))
					|| ((mc[l] == 1) && ((Ns[l] + Nc[l+1]) % 2 == 0)) ) {	// alligned (Nd is odd)
			Nd[l] = ((int) Ac[l]) * 2 + 1;
		} else {														// staggered (Nd is even)
			Nd[l] = ((int) (Ac[l] + 0.5)) * 2;			// approximate
			// Exactly, Nd should be: (nd-0.5)^2 + 0.5^2 = A^2;  nd = 0.5 + sqrt(A^2 - 0.5^2);  Nd = nd * 2;
		}
	}

	Nl = get1Di(L);				// size of lat.inh. cl: Nl[L]
	for (l = 0; l < L; l++)  Nl[l] = ((int) Al[l]) * 2 + 1;

	Mas = get1Di(L);			// no. of active as[L][Ks][Kc]: Mas[L]


  ////  get memory for arrays

  ///  fixed connections -- get memory & calculate
	cg = get2D(Ng, Ng);
	make_cg(cg, Ng, Ag, Ago, gamcg);

	cs = getp(L);
	for (l = 0; l < L; l++) {
		cs[l] = get2D(Na[l], Na[l]);
		make_c(cs[l], Na[l], As[l], gamS[l], coS[l]);
	}

	ac = getp(L);
	for (l = 0; l < L - 1; l++) {
		ac[l] = get2D(Nd[l], Nd[l]);
		make_ac(ac[l], Nd[l], Ac[l], Aco[l], gamC[l], gamCo[l], gamCp[l], coC[l], pnC[l]);
	}

	as = (double *****)getp(L);		// actual memories will be allocated later
	bs = (double **)getp(L);		// actual memories will be allocated later

	if (learn > 0) {
		learn1 = learn - 1;
		if (learn1 == 0)  learn1 = 1;

		asini = getp(L);
		cl = getp(L);
		for (l = learn1; l < L; l++) {
			asini[l] = get2D(Nasini[l], Nasini[l]);
			make_c(asini[l], Nasini[l], Asini[l], gamSini[l], coSini[l]);

			cl[l] = get2D(Nl[l], Nl[l]);
			make_c(cl[l], Nl[l], Al[l], gaml[l], col[l]);
		}
	}

   // memories for responses -- fixed
	u0 = geta3D(1, N0, N0, 1);
	mat = get2Di(Kc[L], Kc[L] + 1);

   // memories for responses -- memory for arrays will be allocated later
	us = (double ****)getp(L);
	uc = (double ****)getp(L+1);	// memory(l>0) will be allocated after learning of us
	uc[0] = geta3D(Kc[0], Nc[0], Nc[0], Kc[0]);

   // calculate parameters
	NsMAX = 0;
	NcMAX = Nc[L];
	for (l = 0; l < L; l++) {
		if (Ns[l] > NsMAX)  NsMAX = Ns[l];
		if (Nc[l] > NcMAX)  NcMAX = Nc[l];
	}

	TTmax = TL;
	if (TT > TTmax)  TTmax = TT;

  ////  start learning
	if (learn > 0) {
   /// get memories for learning (including test after learning)
		KsMAX = 0;			// to be recalculated after finishing the learning
		for (l = 0; l < L; l++) {
			if (Ksmax[l] > KsMAX)  KsMAX = Ksmax[l];
		}
		printu(fpa, u0, N0, 1, " ", 0, KsMAX);		// Initialize printu().  No output.
		if (learn > 1) {
		  ////  read connections up to the (learn-1)th stage
		   // Also get memories for as, bs, Mas[Na][Na] array.
			rwcon(fpC, fnameC, as, bs, &cat, Na, Kc, Ks, Mas, L, -learn);

		  /// get memories for responses, for stages that finised learning
			for (l = 0; l < learn - 1; l++) {
				us[l] = geta3D(Ks[l], Ns[l], Ns[l], Ks[l]);
				uc[l+1] = geta3D(Kc[l+1], Nc[l+1], Nc[l+1], Kc[l+1]);
			}
		}

	// memories for variable connections -- expandable
		for (l = learn - 1; l < L; l++) {
			as[l] = geta4D1(Ksmax[l]);	// memory will be allocated in Sorg
		}

	// memories for variable connections -- fixed to the maximum possible values
		for (l = learn - 1; l < L; l++) {
			bs[l] = get1D(Ksmax[l]);
		}

	// memories for responses -- fixed to the maximum possible values
		cat = get1Di(Ksmax[L-1]);

	// memories for responses -- expandable
		for (l = learn - 1; l < L; l++) {
			us[l] = geta3D(Ks[l], Ns[l], Ns[l], Ksmax[l]);	// may be expanded after Sorg
		}

   ///  start learning
		if (learn == 1) {
		   // Train the first layer (Us1, i.e., ll=0) to detect edges
			add4D2(as[0], Ks[0], Kc[0], Ks[0]);	// create pointer array as[0][Ks[0]][Kc[0]]
			for (k = 0; k < Ks[0]; k++) {
				angle(u0[0], N0, Ks[0], k, 180.0, 0.5, -0.5);	// generate a straight edge
				onoff(uc[0], u0[0], cg, N0, Nc[0], Ng);
				delta(as[0], bs[0], uc[0], cs[0], Na[0], Nc[0], Kc[0], Ns[0],
						k, Ns[0]/2, Ns[0]/2, KsMAX, qs[0]);
			}
			uc[1] = geta3D(Kc[1], Nc[1], Nc[1], Kc[1]);	// get memory for uc
		}

		KsL = 0;
		repeatL = 0;
		for (ll = learn1; ll < L; ll++) {	   // Learning of Us(ll+1)
			Ks[ll] = 0;
			for (t = 0; t < TL * repeat[ll]; t++) {			// repeat while (Ks[L-1]!=KsL)
				pat(u0[0], N0, sat0, t % TL, TL0, TL, TTmax, fnameL);
				if (ll == 1  &&  iprl[0] == 1) {
					fprintf(fpa, "\nLEARNING FOR LAYER %d,   time=%d\n", ll + 1, t);
					printu(fpa, u0, N0, 1, "U", 0, 0);
				}
				onoff(uc[0], u0[0], cg, N0, Nc[0], Ng);
				if (ll == 1  &&  iprl[1] == 1) {
					printu(fpa, uc[0], Nc[0], 2, "UG", -1, 0);
				}
				for (l = 0; l < ll; l++) {
					Sres(us[l], uc[l], as[l], cs[l], bs[l], Nc[l], Kc[l], Ns[l], Ks[l],
							 Na[l], NcMAX, NsMAX, thetar[l]);
					if (l == ll-1  &&  iprl[l*2+2] == 1) {
						fprintf(fpa, "\nLEARNING FOR LAYER %d,   time=%d\n", ll + 1, t);
						printu(fpa, us[l], Ns[l], Ks[l], "Us", l+1, 0);
					}
					Cres(uc[l+1], us[l], ac[l], Ns[l], Nc[l+1], mc[l], Kc[l+1], Nd[l], NsMAX);
					if (l == ll-1  &&  iprl[l*2+3] == 1) {
						fprintf(fpa, "\nLEARNING FOR LAYER %d,   time=%d\n", ll + 1, t);
						printu(fpa, uc[l+1], Nc[l+1], Kc[l+1], "Uc", l+1, 0);
					}
				}
				if (ll != L - 1) {
					Ks[ll] = Sorg(as[ll], bs[ll], uc[ll], cs[ll], asini[ll], cl[ll], us[ll],
							Nc[ll], Kc[ll], Ns[ll], Ks[ll], Na[ll], Nasini[ll], Nl[ll],
							D[ll], qs[ll], thetal[ll], Ksmax[ll], KsMAX, NcMAX, NsMAX, thrini[ll]);
					Kc[ll+1] = Ks[ll];
				} else {
					Cat = t % Kc[L];
					Ks[ll] = SorgL(as[ll], bs[ll], cat, Cat, uc[ll], cs[ll], asini[ll], us[ll],
							Nc[ll], Kc[ll], Ns[ll], Ks[ll], Na[ll], Nasini[ll],
							qs[ll], thetal[ll], Ksmax[ll], KsMAX, NcMAX, NsMAX, fpa);
				}
				add3D(us[ll], Ks[ll], Ns[ll], Ns[ll], Ksmax[ll]);		// add memory to us
				PRINTF("   Learned:    layer=%d,    t=%d,    Ks=%d\n", ll+1, t, Ks[ll]);
				fprintf(fpa, "   Learned:    layer=%d,    t=%d,    Ks=%d\n", ll+1, t, Ks[ll]);
				if ((ll == L - 1) && (t % TL == TL - 1)) {
					repeatL++;
					if (Ks[L - 1] == KsL)  break;		// no error in the last cycle
					KsL = Ks[L - 1];
				}
			}
			uc[ll+1] = geta3D(Kc[ll+1], Nc[ll+1], Nc[ll+1], Kc[ll+1]);	// get memory for uc
			PRINTF(" Finished learning of layer %d,    Ks=%d\n\n", ll+1, Ks[ll]);
			fprintf(fpa, " Finished learning of layer %d,    Ks=%d\n\n", ll+1, Ks[ll]);
			PRINTF("#################################################################\n\n");
			fprintf(fpa, "################################################################################\n\n");
		}
		PRINTF("  Number of repetition in the learning for layer L:  repeatL=%d\n", repeatL);
		fprintf(fpa, "  Number of repetition in the learning for layer L:  repeatL=%d\n", repeatL);
		if (Ks[L - 1] != KsL) {
			PRINTF("    more repetition might be necessary\n");
			fprintf(fpa, "    more repetition might be necessary\n");
		}
		PRINTF("\n");
		fprintf(fpa, "\n");
	}


  ////  read/write connections  ====>  get Mas arrays [Na][Na] when read
   // Also get memories for as, bs and cat, when read (recognition)
	rwcon(fpC, fnameC, as, bs, &cat, Na, Kc, Ks, Mas, L, learn); 	// if (learn>0) write, else read

  /// get memories for responses, in case of recognition without learning
	KsMAX = 0;							// recalculate KsMAX
	for (l = 0; l < L; l++) {
		if (Ks[l] > KsMAX)  KsMAX = Ks[l];
	}
	outH(fpN, u0, N0, 1, " ", 0, KsMAX);		// Initialize outH().  No output.
	if (learn == 0) {
		for (l = 0; l < L; l++) {
			us[l] = geta3D(Ks[l], Ns[l], Ns[l], Ks[l]);
			uc[l+1] = geta3D(Kc[l+1], Nc[l+1], Nc[l+1], Kc[l+1]);
		}
		printu(fpa, u0, N0, 1, " ", 0, KsMAX);		// Initialize printu().  No output.
	}

  ////  test recognition
	fprintf(fpa, "Scale of the network:\n");
	for (l = 0; l < L; l++) {
		fprintf(fpa, "    Stage %d:        Ks = %d,      Mas = %d\n", l+1, Ks[l], Mas[l]);
	}
	PRINTF("\n");

	for (tt = -TTL; tt < TT; tt++) {		// test TTL+TT times

		if (tt < 0) {
			t = tt + TTL;					// test for learned patterns

			PRINTF("Test response for LEARNED pattern      t=%d\n", t);
			if (Iprt != 0) {		// file output of responses of some layers
				fprintf(fpa, "\nTest response for LEARNED pattern       t=%d\n", t);
			}

		} else {
			t = tt;							// test for unlearned patterns from TT0

			PRINTF("Test response for TEST pattern       t=%d\n", t);
			if (Iprt != 0) {		// file output of responses of some layers
				fprintf(fpa, "\nTest response for TEST pattern      t=%d\n", t);
			}

		}


	   // open file "fnameN" (hexadecimal output)
		if (ifpNall != 1) {
			if (tt >= ifpN[0]  &&  tt <= ifpN[1]) {
				fnameN[3] = '0' + ((tt - ifpN[0]) / 10);		// fnameN[] = "out  .eps"
				fnameN[4] = '0' + ((tt - ifpN[0]) % 10);		// fnameN[] = "out  .eps"
				if ((fpN = fopen(fnameN, "w")) == NULL) {
					PRINTF("file open error fpN\n");
					EXIT;
				}
				Npages = 0;							// ==0:  write header (eps) to fpN
			}
		} else {
			if (tt == ifpN[0]) {
				if ((fpN = fopen(fnameNall, "w")) == NULL) {
					PRINTF("file open error fpN\n");
					EXIT;
				}
				Npages = ifpN[1] - ifpN[0] + 1;			// >0:  write header (ps) to fpN
			}
		}
	   // write header (and/or %%Page ) to fpN
		if (fpN != NULL) {
			outHhead(fpN, fnameNhead, Npages);
		}


		if (tt < 0) {						// learned patterns
			pat(u0[0], N0, sat0, t, TL0, TL, TTmax, fnameL);
		} else {							// unlearned patterns
			pat(u0[0], N0, sat0, t, TT0, TT, TTmax, fnameT);
		}
		if (iprt[0] == 1)  printu(fpa, u0, N0, 1, "U", 0, 0);
		if (fpN != NULL)  outH(fpN, u0, N0, 1, "U", 0, 0);

		Cat = t % Kc[L];

		onoff(uc[0], u0[0], cg, N0, Nc[0], Ng);
		if (iprt[1] == 1)  printu(fpa, uc[0], Nc[0], 2, "UG", -1, 0);
		if (fpN != NULL)  outH(fpN, uc[0], Nc[0], 2, "UG", -1, 0);

		for (l = 0; l < L; l++) {
			Sres(us[l], uc[l], as[l], cs[l], bs[l], Nc[l], Kc[l], Ns[l], Ks[l],
					Na[l], NcMAX, NsMAX, thetar[l]);
			if (iprt[l*2+2] / 2 == 1)  printuadd(fpa, us[l], Ns[l], "Us", l+1, 0, Ks[l]-1);
			if (iprt[l*2+2] % 2 == 1)  printu(fpa, us[l], Ns[l], Ks[l], "Us", l+1, 0);

			if (l != L - 1) {
				Cres(uc[l+1], us[l], ac[l], Ns[l], Nc[l+1], mc[l], Kc[l+1], Nd[l], NsMAX);
			} else {
				recCat = CresL(uc[l+1], us[l], cat, Ns[l], Ks[l], Kc[l+1]);
			}
			if (iprt[l*2+3] / 2 == 1)  printuadd(fpa, uc[l+1], Nc[l+1], "Uc", l+1, 0, Kc[l+1]-1);
			if (iprt[l*2+3] % 2 == 1)  printu(fpa, uc[l+1], Nc[l+1], Kc[l+1], "Uc", l+1, 0);
			if (fpN != NULL)  outH(fpN, uc[l+1], Nc[l+1], Kc[l+1], "UC", l+1, 0);
		}
		if (fpN != NULL) {
			fprintf(fpN, "\ndisplay\n");			// "display"
			if (ifpNall != 1) {
				fclose(fpN);
				fpN = NULL;
			} else {
				fprintf(fpN, "\nshowpage\n");		// "showpage"
				if (tt == ifpN[1]) {
					fclose(fpN);
					fpN = NULL;
				}
			}
		}

		if (recCat != Cat) {
			PRINTF("             Erroneous response:   recCat=%d   (Input Cat=%d)\n", recCat, Cat);
			if (Iprt == 0) {		// no file-output of response of any layer
				fprintf(fpa,
					"Test  t=%d     Erroneous response:   recCat=%d   (Input Cat=%d)\n",
					tt, recCat, Cat);
			} else {
				fprintf(fpa,
					"             Erroneous response:   recCat=%d   (Input Cat=%d)\n",
					recCat, Cat);
			}
			ErrOut(u0, N0, fnameE, fnameEh, 0);
		}

		score(mat, recCat, Cat, Kc[L]);
		if (tt < 0) {
			if (tt == -1) {			// TTL > 0
				fprintf(fpa, "\n\nRecognition rate for Learned patterns:");
				PRINTF("\n\nRecognition rate for Learned patterns:\n");
				printsc(fpa, mat, Kc[L]);
				for (i = 0; i < Kc[L]; i++) {
					for (j = 0; j <= Kc[L]; j++) {
						mat[i][j] = 0;
					}
				}
			}
		} else if ( (tt + 1 == TT) || ((tt + 1) % 1000 == 0) ) {
			fprintf(fpa, "\n\nRecognition rate for Unlearned patterns  (from t=%d to t=%d):", TT0, t);
			PRINTF("\n\nRecognition rate for Unlearned patterns  (from t=%d to t=%d):\n", TT0, t);
			printsc(fpa, mat, Kc[L]);
		}
	}

	ErrOut(u0, N0, fnameE, fnameEh, 1);		// close file fpE
	if (fpa != NULL)  fclose(fpa);
}													//  end of main

#endif // NEOCOG_MATLAB

/*************************************************************************/

////  Function Definition

double ***getp(			// memory for 1-D pointer array **u[L]
	int		L
	)
{
	double	***u;

	if( L <= 0 )	return NULL;
	u = (double ***)MALLOC(sizeof(double **) * L);
	if (u == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}
	return  u;
}



double *get1D(			// get memory for 1-D array m[N]
	int		N
	)
{
	double		*u;

	if( N <= 0 )	return NULL;
	u = (double *)MALLOC(sizeof(double) * N);
	if (u == NULL)  PRINTF("Memory allocation error!\n");
	return  u;
}



int *get1Di(			// get memory for 1-D int array m[N]
	int		N
	)
{
	int		*m;

	if( N <= 0 )	return NULL;
	m = (int *)MALLOC(sizeof(int) * N);
	if (m == NULL)  PRINTF("Memory allocation error!\n");
	return  m;
}



double **get2D(			// get memory for 2-D array u[N1][N2]
	int		N1,
	int		N2
	)
{
	int		i;
	double	**u;

	if( N1 * N2 <= 0 )	return NULL;
	u = (double **)MALLOC(sizeof(double *) * N1);
	if (u == NULL)  PRINTF("Memory allocation error!\n");

	u[0] = (double *)MALLOC(sizeof(double) * N1 * N2);
	if (u[0] == NULL)  PRINTF("Memory allocation error!\n");

	for (i = 1; i < N1; i++) {
		u[i] = u[i-1] + N2;
	}

	return  u;
}



int **get2Di(			// get memory for 2-D int array m[N1][N2]
	int		N1,
	int		N2
	)
{
	int		i;
	int		**m;

	if( N1 * N2 <= 0 )	return NULL;
	m = (int **)MALLOC(sizeof(int *) * N1);
	if (m == NULL)  PRINTF("Memory allocation error!\n");

	m[0] = (int *)MALLOC(sizeof(int) * N1 * N2);
	if (m[0] == NULL)  PRINTF("Memory allocation error!\n");

	for (i = 1; i < N1; i++) {
		m[i] = m[i-1] + N2;
	}

	return  m;
}



unsigned char **get2Dc(	// get memory for 2-D unsigned char array ch[N1][N2]
	int		N1,
	int		N2
	)
{
	int		i;
	unsigned char	**ch;

	if( N1 * N2 <= 0 )	return NULL;
	ch = (unsigned char **)MALLOC(sizeof(unsigned char *) * N1);
	if (ch == NULL)  PRINTF("Memory allocation error!\n");

	ch[0] = (unsigned char *)MALLOC(sizeof(unsigned char) * N1 * N2);
	if (ch[0] == NULL)  PRINTF("Memory allocation error!\n");

	for (i = 1; i < N1; i++) {
		ch[i] = ch[i-1] + N2;
	}

	return  ch;
}



double ***geta3D(		// get addable memory for 3-D array u[K][N1][N2]
	int		K,			// initial array size
	int		N1,
	int		N2,
	int		Kmax		// maximum possible array size: u[Kmax][N1][N2]
	)
{
	int		i, j;
	double	***u;

	if( Kmax * N1 * N2 <= 0 )	return NULL;
	if (K > Kmax) {
		PRINTF("Cannot get array larger than Kmax:  K = %d", K);
		EXIT;
	}

	u = (double ***)MALLOC(sizeof(double **) * Kmax);
	if (u == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}

	u[0] = (double **)MALLOC(sizeof(double *) * K * N1);
	if (u[0] == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}
	for (i = 1; i < K; i++) {
		u[i] = u[i-1] + N1;
	}
	for (i = K; i < Kmax; i++) {
		u[i] = NULL;
	}

	u[0][0] = (double *)MALLOC(sizeof(double) * K * N1 * N2);
	if (u[0][0] == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}
	for (i = 1; i < K; i++) {
		u[i][0] = u[i-1][0] + N1 * N2;
	}
	for (i = 0; i < K; i++) {
		for (j = 1; j < N1; j++) {
			u[i][j] = u[i][j-1] + N2;
		}
	}

	return  u;
}



void add3D(				// add memory to 3-D array u[K][N1][N2]
	double	***u,		// (in) u = (double ***)malloc(sizeof(double **) * Kmax)
	int		Knew,		// new array size
	int		N1,
	int		N2,
	int		Kmax		// maximum possiboe array size
	)
{
	int		i, j;
	int		Kold;		// previous array size
	int		Kadd;		// size of the added memory

	for (Kold = 0; Kold < Kmax; Kold++) {
		if (u[Kold] == NULL)  break;
	}
	Kadd = Knew - Kold;
	if (Kadd <= 0)  return;

	if (Knew > Kmax) {
		PRINTF("Cannot get array larger than Kmax:  Knew = %d", Knew);
		EXIT;
	}

	u[Kold] = (double **)MALLOC(sizeof(double *) * Kadd * N1);
	if (u[Kold] == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}
	for (i = Kold + 1; i < Knew; i++) {
		u[i] = u[i-1] + N1;
	}

	u[Kold][0] = (double *)MALLOC(sizeof(double) * Kadd * N1 * N2);
	if (u[Kold][0] == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}
	for (i = Kold + 1; i < Knew; i++) {
		u[i][0] = u[i-1][0] + N1 * N2;
	}
	for (i = Kold; i < Knew; i++) {
		for (j = 1; j < N1; j++) {
			u[i][j] = u[i][j-1] + N2;
		}
	}
}



double ****geta4D1(		// get addable memory for 1-D pointer array ***a[]
	int		Kmax		// maximum possible pointer array size: ***a[Kmax]
	)
{
	int		i;
	double	****a;

	a = (double ****)MALLOC(sizeof(double *) * Kmax);
	if (a == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}
	for (i = 0; i < Kmax; i++) {
		a[i] = NULL;
	}
	return  a;
}



void add4D2(			// add memory to 2-D pointer array **a[K][Kc]
	double	****a,		// (in) a = (double ****)malloc(sizeof(double *) * Kmax)
	int		Knew,		// new pointer array size: **a[Knew][Kc]
	int		Kc,
	int		Kmax		// maximum possible array size: **a[Kmax][Kc]
	)
{
	int		i, j;
	int		Kold;			// previous array size
	int		Kadd;			// added memory size

	for (Kold = 0; Kold < Kmax; Kold++) {
		if (a[Kold] == NULL)  break;
	}
	Kadd = Knew - Kold;
	if (Kadd <= 0)  return;

	if (Knew > Kmax) {
		PRINTF("Cannot get array larger than Kmax:  Knew = %d", Knew);
		EXIT;
	}

	a[Kold] = (double ***)MALLOC(sizeof(double *) * Kc * Kadd);
	if (a[Kold] == NULL) {
		PRINTF("Memory allocation error!\n");
		EXIT;
	}
	for (i = Kold + 1; i < Knew; i++) {
		a[i] = a[i-1] + Kc;
	}

	for (i = Kold; i < Knew; i++) {
		for (j = 0; j < Kc; j++) {
			a[i][j] = NULL;
		}
	}
}



void Sres(				// response of S-layer
	double	***us,		// (out) response of S-cells: us[Ks][Ns][Ns]
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	****as,		// (in) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	**cs,		// (in) exc. connections to V-cells: cs[Na][Na]
	double	*bs,		// (in) inh. connections to S-cells: bs[Ks]
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Ks,			// numberof cell-planes in S-layer
	int		Na,			// size of the connections (as & cs): (na*2) or (na*2+1)
	int		NcMAX,		// max of Nc (and Ng) among all stages
	int		NsMAX,		// max of Ns among all stages
	double	theta		// threshold of an S-cell
	)
{
	STATIC int		init = 1;	// flag indicationg initial entry
	STATIC int		*nsmax;		// (work) upper bound of the loop: nsmax[NcMAX]
	STATIC int		*nsmin;		// (work) lower bound of the loop: nsmin[NcMAX]
	STATIC int		*namax;		// (work) upper bound of the conection: namax[NcMAX]
	STATIC double	**vs2;		// (work) response**2 of V-cells: vs2[NsMAX][NsMAX]
	int		ns2 = Ns - 1;
	int		nc2 = Nc - 1;
	int		na2 = Na - 1;
	int		na0 = (na2 - nc2 + ns2) / 2;	// na0=an-cn+sn; sn=ns2*0.5, cn=nc2*0.5, an=na2*0.5
	double	theta1 = theta / (1 - theta);
	int		is, js, ks, ic, jc, kc, ia, ja, ic1;
	double	uu, uu2, aas, vv, thv;

	if (init != 0) {
		nsmax = get1Di(NcMAX);
		nsmin = get1Di(NcMAX);
		namax = get1Di(NcMAX);
		vs2 = get2D(NsMAX, NsMAX);
		init = 0;
	}
	for (ic = 0; ic < Nc; ic++) {
		nsmax[ic] = ic + na0;					// = ic - cn + an + sn;
		if (nsmax[ic] > ns2)  nsmax[ic] = ns2;
		ic1 = nc2 - ic;
		nsmin[ic1] = ns2 - nsmax[ic];
		namax[ic1] = ic1 - nsmin[ic1] + na0;	// =(ic-cn) - (nsmin[ic]-sn) + an;
	}
	for (ks = 0; ks < Ks; ks++) {				// claer us & usma
		for (is = 0; is < Ns; is++) {
			for (js = 0; js < Ns; js++) {
				us[ks][is][js] = 0.0;
			}
		}
	}
	for (is = 0; is < Ns; is++) {				// claer vs2
		for (js = 0; js < Ns; js++) {
			vs2[is][js] = 0.0;
		}
	}
	for (kc = 0; kc < Kc; kc++) {	// send signals from uc to us & vs2
		for (ic = 0; ic < Nc; ic++) {
			for (jc = 0; jc < Nc; jc++) {
				if ((uu = uc[kc][ic][jc]) > 0.0) {
					uu2 = uu * uu;
					for (is = nsmin[ic], ia = namax[ic]; is <= nsmax[ic]; is++, ia--) {
						for (js = nsmin[jc], ja = namax[jc]; js <= nsmax[jc]; js++, ja--) {
							vs2[is][js] += cs[ia][ja] * uu2;
						}
					}
					for (ks = 0; ks < Ks; ks++) {
						if (as[ks][kc] != NULL) {		// sarting adress in the memory as
							for (is = nsmin[ic], ia = namax[ic]; is <= nsmax[ic]; is++, ia--) {
								for (js = nsmin[jc], ja = namax[jc]; js <= nsmax[jc]; js++, ja--) {
									if ((aas = as[ks][kc][ia][ja]) > 0.0) {
										us[ks][is][js] += aas * uu;
									}
								}
							}
						}
					}
				}
			}
		}
	}
	for (is = 0; is < Ns; is++) {
		for (js = 0; js < Ns; js++) {
			if ((vv = vs2[is][js])> 0.0) {
				thv = theta * sqrt(vv);
				for (ks = 0; ks < Ks; ks++) {
					if ((uu = us[ks][is][js]) > 0.0) {
						uu = (1.0 + uu)/(1.0 + thv * bs[ks]);
						if (uu > 1.0) {
							us[ks][is][js] = theta1 * (uu - 1.0);
						} else {
							us[ks][is][js] = 0.0;
						}
					}
				}
			}
		}
	}
}



void Cres(				// response of C-layer
	double	***uc,		// (out) response of C-cells: uc[K][Nc][Nc]
	double	***us,		// (in) response of S-cells: us[K][Ns][Ns]
	double	**ac,		// (in) exc. connections to C-cells: ac[Nd][Nd]
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		mc,			// ratio of thinning-out
	int		K,			// number of cell-planes in S- and C-layers (K=Ks=Kc)
	int		Nd,			// size of the connections (ac): (nd*2) or (nd*2+1)
	int		NsMAX		// max of Ns among all stages
	)
{
	STATIC int	init = 1;	// flag indicationg initial entry
	STATIC int	*ncmax;		// (work) upper bound of the loop: ncmax[NsMAX]
	STATIC int	*ncmin;		// (work) lower bound of the loop: ncmin[NsMAX]
	STATIC int	*namax;		// (work) upper bound of the conection: namax[NsMAX]
	int		nc2 = Nc - 1;			// cn=nc2*0.5
	int		ns2 = Ns - 1;			// sn=ns2*0.5
	int		nd2 = Nd - 1;			// dn=nd2*0.5
	int		nd0 = (nc2 * mc - ns2 + nd2) / 2;	// =dn0=cn*mc-sn+dn=(nc2*mc-ns2+nd2)*0.5
	int		k, ic, jc, is, js, ia, ja, is1;
	double	uu;

	if (init != 0) {
		ncmax = get1Di(NsMAX);
		ncmin = get1Di(NsMAX);
		namax = get1Di(NsMAX);
		init = 0;
	}
	for (is = 0; is < Ns; is++) {
		ncmax[is] = (is + nd0) / mc;	// =floor((is-sn+dn+cn*mc)/mc)
		if ((is + nd0) < 0)  ncmax[is]--;	// not necessary, except for Ns>>Nc
		if (ncmax[is] > nc2)  ncmax[is] = nc2;
		is1 = ns2 - is;
		ncmin[is1] = nc2 - ncmax[is];
		namax[is1] = is1 - ncmin[is1] * mc + nd0;
									// = (is1-sn) - (ncmin[is1]-cn)*mc + dn;
									// = is1 - ncmin[is1]*mc + cn*mc-sn+dn;
									// = is1 - ncmin[is1]*mc + dn0;
	}
	for (k = 0; k < K; k++) {
		for (ic = 0; ic < Nc; ic++) {		// claer uc
			for (jc = 0; jc < Nc; jc++) {
				uc[k][ic][jc] = 0.0;
			}
		}
		for (is = 0; is < Ns; is++) {
			for (js = 0; js < Ns; js++) {
				if ((uu = us[k][is][js]) > 0.0) {
					for (ic = ncmin[is], ia = namax[is]; ic <= ncmax[is]; ic++, ia -= mc) {
						for (jc = ncmin[js], ja = namax[js]; jc <= ncmax[js]; jc++, ja -= mc) {
							uc[k][ic][jc] += ac[ia][ja] * uu;
						}
					}
				}
			}
		}
		for (ic = 0; ic < Nc; ic++) {			// calculate the response of uc
			for (jc = 0; jc < Nc; jc++) {
				if ((uu = uc[k][ic][jc]) > 0.0) {
					uc[k][ic][jc] = uu / (1.0 + uu);
				} else {
					uc[k][ic][jc] = 0.0;
				}
			}
		}
	}
}



int CresL(				// response of C-layer without lateral inhibition
	double	***uc,		// (out) response of C-cells: uc[Kc][1][1]
	double	***us,		// (in) response of S-cells: us[K][Ns][Ns]
	int		*cat,		// (in) category of the S-plane: cat[Ks]
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Ks,			// number of cell-planes in S-layers (Ks)
	int		Kc			// number of cell-planes in C-layers (Kc)
	)
{
	int		is, js, ks, kc, Cat;
	double	Usmax;

	for (kc = 0; kc < Kc; kc++) {
		uc[kc][0][0] = 0.0;
	}

	Cat = -1;
	Usmax = 0.0;
	for (ks = 0; ks < Ks; ks++) {
		for (is = 0; is < Ns; is++) {
			for (js = 0; js < Ns; js++) {
				if (us[ks][is][js] > Usmax) {
					Usmax = us[ks][is][js];
					Cat = cat[ks];
				}
			}
		}
	}

//	PRINTF("Cat=%d,  Usmax=%f\n", Cat, Usmax);

	if (Cat >= 0) {
		uc[Cat][0][0] = Usmax;
	}
	return Cat;
}



void onoff(				// response of on- & off-center cell layers
	double	***uc,		// (out) response of on- & off-center cells: uc[2][Nc][Nc]
	double	**us,		// (in) response of receptor cells: us[Ns][Ns]
	double	**cg,		// (in) connections for on-center: cg[Ng][Ng]
	int		Ns,			// size of the receptor layer: (ns*2) or (ns*2+1)
	int		Nc,			// size of the on- & off-center cell layer: (nc*2) or (nc*2+1)
	int		Ng			// size of the connections (ag): (ng*2) or (ng*2+1)
	)
{
	STATIC int	init = 1;	// flag indicationg initial entry
	STATIC int	*ncmax;		// (work) upper bound of the loop: ncmax[Ns]
	STATIC int	*ncmin;		// (work) lower bound of the loop: ncmin[Ns]
	STATIC int	*namax;		// (work) upper bound of the conection: namax[Ns]
	int		ns2 = Ns - 1;
	int		nc2 = Nc - 1;
	int		ng2 = Ng - 1;
	int		ng0 = (ng2 - ns2 + nc2) / 2;	// ng0=gn-sn+cn; sn=ns2*0.5, cn=nc2*0.5, gn=ng2*0.5
	int		ic, jc, is, js, ia, ja, is1;
	double	uu;

	if (init != 0) {
		ncmax = get1Di(Ns);
		ncmin = get1Di(Ns);
		namax = get1Di(Ns);
		init = 0;
	}
	for (is = 0; is < Ns; is++) {
		ncmax[is] = is + ng0;
		if (ncmax[is] > nc2)  ncmax[is] = nc2;
		is1 = ns2 - is;
		ncmin[is1] = nc2 - ncmax[is];
		namax[is1] = is1 - ncmin[is1] + ng0;	// =(is-sn) - (ncmin[is]-cn) + gn;
	}
	for (ic = 0; ic < Nc; ic++) {
		for (jc = 0; jc < Nc; jc++) {
			uc[0][ic][jc] = 0.0;
		}
	}
	for (is = 0; is < Ns; is++) {
		for (js = 0; js < Ns; js++) {
			if ((uu = us[is][js]) > 0.0) {
				for (ic = ncmin[is], ia = namax[is]; ic <= ncmax[is]; ic++, ia--) {
					for (jc = ncmin[js], ja = namax[js]; jc <= ncmax[js]; jc++, ja--) {
						uc[0][ic][jc] += cg[ia][ja] * uu;
					}
				}
			}
		}
	}
	for (ic = 0; ic < Nc; ic++) {
		for (jc = 0; jc < Nc; jc++) {
			if ((uu = uc[0][ic][jc]) > 0.0001) {
				uc[1][ic][jc] = 0.0;
			} else if (uu < -0.0001) {
				uc[0][ic][jc] = 0.0;
				uc[1][ic][jc] = -uu;
			} else {
				uc[0][ic][jc] = 0.0;
				uc[1][ic][jc] = 0.0;
			}
		}
	}
}



void make_cg(			// make connections for on-center
	double	**cg,		// (out) Connections: cg[Ng][Ng]
	int		Ng,			// size of connection cg: (ng*2) or (ng*2+1)
	double	Ag,			// radius of the off-periphery (connect. area)
	double	Ago,		// radius of the on-center component
	double	gamcg		// spatial decay at x=Ag
	)
{
	int		i, j;
	double	x, y, y2, r2, Spos, Sneg, SSpos, SSneg, cc;
	double	gn = (Ng - 1) * 0.5;
	double	Ag2 = Ag * Ag;
	double	Ago2 = Ago * Ago;
	double	C = 1 / Ag2;
	double	Co = 1 / Ago2;

	Spos = 0.0;
	Sneg = 0.0;
	for (i = 0; i < Ng; i++) {
		y = i - gn;
		y2 = y * y;
		for (j = 0; j < Ng; j++) {
			x = j - gn;
			r2 = y2 + x * x;
			if (r2 <= Ag2) {
				cc = Co * pow(gamcg, r2 / Ago2) - C * pow(gamcg, r2 / Ag2);
				cg[i][j] = cc;
				if (cc > 0.0)  Spos += cc;
				else           Sneg -= cc;
			} else {
				cg[i][j] = 0.0;
			}
		}
	}
	SSpos = 1.0 / Spos;
	SSneg = 1.0 / Sneg;
	for (i = 0; i < Ng; i++) {
		for (j = 0; j < Ng; j++) {
			if (cg[i][j] > 0.0) {
				cg[i][j] *= SSpos;
			} else {
				cg[i][j] *= SSneg;
			}
		}
	}
}



void make_ac(			// make center surround connections ac
	double	**ac,		// (out) Connections: ac[Na][Na]
	int		Na,			// size of the connection ac: (na*2) or (na*2+1)
	double	A,			// radius of the off-periphery (connect. area)
	double	Ao,			// radius of the on-center component
	double	gamma,		// spatial decay of exc. connect. at x=Ao
	double	gammao,		// spatial decay of inh. connect. at x=Ao
	double	gammap,		// spatial decay of inh. connect. at x=A
	double	co,			// value of ac[][] at the center
	double	pn			// ratio of negative to positibe components
	)
{
	int		i, j;
	double	x, y, y2, r2, Spos, Sneg, SSneg, cc, C, Co, CCo, CC;
	double	an = (Na - 1) * 0.5;
	double	A2 = A * A;
	double	Ao2 = Ao * Ao;

	if (Ao2 < A2) {
		C = 1 / A2;
		Co = 1 / Ao2;
		CCo = co * Co / (Co - C);
		CC = co * C / (Co - C);
	}

	Spos = 0.0;
	Sneg = 0.0;
	for (i = 0; i < Na; i++) {
		y = i - an;
		y2 = y * y;
		for (j = 0; j < Na; j++) {
			x = j - an;
			r2 = y2 + x * x;
			if (r2 <= Ao2) {
				cc = co * pow(gamma, r2 / Ao2);
				ac[i][j] = cc;
				Spos += cc;
			} else if (r2 <= A2) {
				cc = CCo * pow(gammao, r2 / Ao2) - CC * pow(gammap, r2 / A2);
				ac[i][j] = cc;
				if (cc > 0.0)  Spos += cc;
				else           Sneg -= cc;
			} else {
				ac[i][j] = 0.0;
			}
		}
	}
	if (Ao2 < A2) {
		SSneg = pn * Spos / Sneg;
		for (i = 0; i < Na; i++) {
			for (j = 0; j < Na; j++) {
				if (ac[i][j] < 0.0) {
					ac[i][j] *= SSneg;
				}
			}
		}
	}
}



void make_c(			// make connections (spatial weight)
	double	**c,		// (out) Connections: c[Na][Na]
	int		Na,			// size of the connection c: (na*2) or (na*2+1)
	double	A,			// radius of the connectable area
	double	gamma,		// spatial decay at x=A
	double  co			// value of c[][] at the center
	)
{
	int		i, j;
	double	x, y, y2, r2;
	double	an = (Na - 1) * 0.5;
	double	A2 = A * A;

	for (i = 0; i < Na; i++) {
		y = i - an;
		y2 = y * y;
		for (j = 0; j < Na; j++) {
			x = j - an;
			r2 = y2 + x * x;
			if (r2 <= A2) {
				c[i][j] = co * pow(gamma, r2 / A2);
			} else {
				c[i][j] = 0.0;
			}
		}
	}
}



void angle(				// draw an angle of ang degree
	double	**u,		// (out) generated pattern: u[N][N]
	int		N,			// size of the layer
	int		K,			// number of different orientations
	int		k,			// sequence number of the orientation of the edge
	double	ang,		// angle (in degree) of the bend
	double	x0,			// x offset of the center of the pattern
	double	y0			// y offset of the center of the pattern
	)
{
	int		i, j;
	double	rad, Rcosa, Rsina, Rcosb, Rsinb, x1, y1, x2, y2;
	double	pi = 3.141592653589793;
	double	radc = k * (2.0 * pi / K);
	double	radang2 = ang * pi / 360.0;		// ang/2 in radian
	double	rada = radc + radang2;
	double	radb = radc - radang2;
	double	radstep = pi * 0.0625;
	double	R = ((double) N) * 0.43;

	Rcosa = R * cos(rada);
	Rsina = R * sin(rada);
	Rcosb = R * cos(radb);
	Rsinb = R * sin(radb);

	for (i = 0; i < N; i++) {			// clear u[][]
		for (j = 0; j < N; j++) {
			u[i][j] = 0.0;
		}
	}

	edge(u, N, Rcosa + x0, Rsina + y0, x0, y0);
	edge(u, N, x0, y0, Rcosb + x0, Rsinb + y0);
	x1 = Rcosb + x0;
	y1 = Rsinb + y0;
	for (rad = radb; rad < rada; rad += radstep) {
		x2 = R * cos(rad) + x0;
		y2 = R * sin(rad) + y0;
		edge(u, N, x1, y1, x2, y2);
		x1 = x2;
		y1 = y2;
	}
	edge(u, N, x1, y1, Rcosa + x0, Rsina + y0);

	for (i = 0; i < N; i++) {			// eliminate rounding error
		for (j = 0; j < N; j++) {
			if (u[i][j] < 0.0001)  u[i][j] = 0.0;
			if (u[i][j] > 0.9999)  u[i][j] = 1.0;
		}
	}
}



void edge(
	double	**u,		// (i/o) generated pattern: u[N][N]
	int		N,			// size of the layer
	double	xs,
	double	ys,			// edge from (xs,ys) to (xe,ye)
	double	xe,
	double	ye
	)
{
	int		i, j, ia, ii;
	double	center, xpol, x1, y1, x2, y2, p, x, y, xa, ya, w;

	center = ((double) N) / 2.0;
	if (fabs(xe - xs) < 0.0001)  return;
	if (xs < xe) {
		xpol = 1.0;
		x1 = center + xs;
		y1 = center - ys;		// i increases from up to down
		x2 = center + xe;
		y2 = center - ye;
	} else {
		xpol = -1.0;
		x1 = center + xe;
		y1 = center - ye;
		x2 = center + xs;
		y2 = center - ys;
	}
	p = (y2 - y1) / (x2 - x1);
	if (x1 < 0.0) {
		y1 = y1 - p * x1;
		x1 = 0.0;
	}
	if (x2 > N) {
		y2 = y2 - p * (x2 - N);
		x2 = N;
	}
	if (fabs(p) < 0.0001)  p = 0.0;
	x = x1;
	y = y1;
	while (x < x2) {
		j = (int) x;
		i = (int) floor(y);
		if ((i == y) && (p < 0.0))  i = i - 1;
		xa = j + 1;
		if (xa > x2)  xa = x2;
		ya = y + p * (xa - x);
		ia = (int) floor(ya);
		if (ia > i) {
			ya = i + 1;
			xa = x + (ya - y) / p;
		} else if (ia < i) {							// p<0.0
			ya = i;
			xa = x + (ya - y) / p;
		}
		w = xpol * (xa - x);
		if (i >= 0) {
			if (i < N) {
				u[i][j] += w * ((y + ya) / 2.0 - i);
				for (ii = 0; ii < i; ii++) {
					u[ii][j] += w;
				}
			} else {
				for (ii = 0; ii < N; ii++) {
					u[ii][j] += w;
				}
			}
		}
		x = xa;
		y = ya;
	}
}



int Sorg(				// self-organization of S-layer
	double	****as,		// (out) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	*bs,		// (out) inh. connections to S-cells: bs[Ks], new Ks
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**cs,		// (in) spatial weight: cs[Na][Na]
	double	**asini,	// (in) exc. initial connections to S-cells: as[Na][Na]
	double	**cl,		// (in) connect. for lat.inh. from us to usini: cl[Nl][Nl]
	double	***us,		// (work) response of S-cells: us[Ks][Ns][Ns], old Ks
	int		Nc,			// size of the C-layer:  (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer:  (ns*2) or (ns*2+1)
	int		Ks,			// number of cell-planes in S-layer, old value
	int		Na,			// size of connections (as & cs):  (na*2) or (na*2+1)
	int		Nasini,		// size of connect. asini:  (nasini*2) or (nasini*2+1)
	int		Nl,			// size of connection cl (lat. inh. from us to usini):  Nl=nl*2+1
	double	D,			// radius of competition area
	double	qs,			// parameter determining the speed of learning
	double	theta,		// threshold of an S-cell
	int		Ksmax,		// maximum possible Ks
	int		KsMAX,		// max of KsMAX among all stages
	int		NcMAX,		// max of Nc (and Ng) among all stages
	int		NsMAX,		// max of Ns among all stages
	double	thrini		// threshold for cutting usini
	)
{
	STATIC int		init = 1;	// flag indicationg initial entry
	STATIC double	**usini;	// response of seed-selecting plane: usini[NsMAX][NsMAX]
	STATIC double	**work;		// for wmax[][] & us[Ks][][]: work[NsMAX][NsMAX]
	STATIC int		**kwmax;	// kwmax[Ns][Ns].  Also used as a print-buffer when fpT!=NULL
	STATIC int		*iseed;		// iseed[Ks], old Ks
	STATIC int		*jseed;		// jseed[Ks], old Ks
	int		ks, iseed0, jseed0;

	if (init != 0) {
		usini = get2D(NsMAX, NsMAX);
		work = get2D(NsMAX, NsMAX);
		kwmax = get2Di(NsMAX, NsMAX);
		iseed = get1Di(KsMAX);
		jseed = get1Di(KsMAX);
		init = 0;
	}

	if (Ks > 0) {
		Sres(us, uc, as, cs, bs, Nc, Kc, Ns, Ks, Na, NcMAX, NsMAX, theta);
										// Sres, using old connections
		seed(iseed, jseed, us, work, kwmax, Ns, Ks, NsMAX, D);		// work=wmax

		for (ks = 0; ks < Ks; ks++) {
			if (iseed[ks] >= 0) {
				delta(as, bs, uc, cs, Na, Nc, Kc, Ns,
						ks, iseed[ks], jseed[ks], KsMAX, qs);
								// reinforce connect. for old cell-planes
								// memory for as[] may be increased by delta()
			}
		}

		Sres(us, uc, as, cs, bs, Nc, Kc, Ns, Ks, Na, NcMAX, NsMAX, theta);
										// Sres, using renewed connections
	}
	Sresini(usini, uc, asini, Nc, Kc, Ns, Nasini, NcMAX, thrini);
	for (ks = 0; ks < Ks; ks++) {
		lat(usini, us[ks], cl, Ns, Nl);			// lat.inh. from us to usini
	}
	for (ks = Ks; ks < Ksmax; ks++) {
		seedini(&iseed0, &jseed0, usini, Ns);	// location of seed = max.-output cell
		if (iseed0 == -1)  break;
		add4D2(as, ks + 1, Kc, Ksmax);		// add Kc elements to pointer array **as
		delta(as, bs, uc, cs, Na, Nc, Kc, Ns, ks, iseed0, jseed0, KsMAX, qs);
								// reinforce connect. for new cell-planes
								// memory for as[] may be increased by delta()
								// memory for bs[] must have been prepared beforehand
		Sres(&work, uc, &as[ks], cs, &bs[ks], Nc, Kc, Ns, 1, Na, NcMAX, NsMAX, theta);
											// Sres, for the new cell-plane
		lat(usini, work, cl, Ns, Nl);		// lat.inh. from new us to usini
	}
	return  ks;
}



void delta(				// strengthen connections
	double	****as,		// (out) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	*bs,		// (out) inh. connections to S-cells: bs[Ks]
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**cs,		// (in) spatial weight: cs[Na][Na]
	int		Na,			// size of the connections (as & cs): (na*2) or (na*2+1)
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		Kc,			// numberof cell-planes in C-layer
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		kseed,		// cell-plane of the seed cell
	int		iseed,		// x-coordinate of the seed cell
	int		jseed,		// y-coordinate of the seed cell
	int		KsMAX,		// max of KsMAX among all stages, because KcMAX <= KsMAX
	double	qs			// parameter determining the speed of learning
	)
{
	STATIC int	init = 1;
	STATIC int	*ucmax;		// ucmax[KsMAX], because KcMAX <= KsMAX
	int		i, j, ic, jc, kc, icmin, icmax, jcmin, jcmax, iamin, jamin, icmin0, jcmin0;
	double	**ma;
	int		m_inc;
	double	uu, b2, cc, aa;
	int 	ns2 = Ns - 1;						// sn=ns2*0.5
	int		nc2 = Nc - 1;						// cn=nc2*0.5
	int		na2 = Na - 1;						// an=na2*0.5
	int		na0 = (nc2 - ns2 + na2) / 2;		// na0=cn-sn+an
	int		na1 = (nc2 - ns2 - na2) / 2;		// na1=cn-sn-an

	if (init != 0) {
		ucmax = get1Di(KsMAX);			// ucmax[KsMAX], because KcMAX <= KsMAX
	}

	icmax = iseed + na0;				// = (iseed - sn) + an + cn;
	if (icmax > nc2)  icmax = nc2;
	icmin = iseed + na1;				// = (iseed - sn) - an + cn;
	icmin0 = icmin;
	if (icmin < 0)  icmin = 0;
	iamin = icmin - icmin0;
	jcmax = jseed + na0;				// = (jseed - sn) + an + cn;
	if (jcmax > nc2)  jcmax = nc2;
	jcmin = jseed + na1;				// = (jseed - sn) - an + cn;
	jcmin0 = jcmin;
	if (jcmin < 0)  jcmin = 0;
	jamin = jcmin - jcmin0;

	m_inc = 0;
	for (kc = 0; kc < Kc; kc++) {
		for (ic = icmin, i = iamin; ic <= icmax; ic++, i++) {
			for (jc = jcmin, j = jamin; jc <= jcmax; jc++, j++) {
				if (uc[kc][ic][jc] > 0.0) {		// if at least one positive cell exists
					ucmax[kc] = 1;
					if (as[kseed][kc] == NULL) {
						m_inc++;
					}
					goto skip;
				}
			}
		}
		ucmax[kc] = 0;
skip:	;
	}
	ma = get2D(Na * m_inc, Na);			// get m_inc arrays [Na][Na]

	for (kc = 0; kc < Kc; kc++) {
		if (ucmax[kc] > 0) {
			if (as[kseed][kc] == NULL) {
				as[kseed][kc] = ma;
				ma += Na;
				for (i = 0; i < Na; i++) {
					for (j = 0; j < Na; j++) {
						as[kseed][kc][i][j] = 0.0;
					}
				}
			}
			for (ic = icmin, i = iamin; ic <= icmax; ic++, i++) {
				for (jc = jcmin, j = jamin; jc <= jcmax; jc++, j++) {
					if ((uu = uc[kc][ic][jc]) > 0.0) {
						if ((cc = cs[i][j]) > 0.0) {
							as[kseed][kc][i][j] += qs * cc * uu;
						}
					}
				}
			}
		}
	}
	b2 = 0.0;
	for (kc = 0; kc < Kc; kc++) {
		if (as[kseed][kc] != NULL) {
			for (i = 0; i < Na; i++) {
				for (j = 0; j < Na; j++) {
					if ((aa = as[kseed][kc][i][j]) > 0.0) {
						if ((cc = cs[i][j]) > 0.0) {
							b2 += aa * aa / cc;
						}
					}
				}
			}
		}
	}
	bs[kseed] = sqrt(b2);
}



void seed(				// select seed cells
	int		*iseed,		// (out) x-coordinate of the k-th cell-plane: iseed[Ks]
	int		*jseed,		// (out) y-coordinate of the k-th cell-plane: jseed[Ks]
	double	***us,		// (in) response of S-cells: us[Ks][Ns][Ns]
	double	**wmax,		// (work) working area: wmax[Ns][Ns]
	int		**kwmax,	// (work) working area: kwmax[Ns][Ns]
	int		Ns,			// size of the S-layer
	int		Ks,			// number of cell-planes in S-layer
	int		NsMAX,		// max of Ns among all stages
	double	D			// radius of competition area
	)
{
	STATIC int	init = 1;	// flag indicationg initial entry
	STATIC int	*nmax;		// (work) upper bound of the competition: nmax[NsMAX]
	STATIC int	*nmin;		// (work) lower bound of the competition: nmin[Ns]
	double	d2 = D * D;		// D**2 in the scale of this layer
	int		id = (int) D;
	int		k, i, j, ii, jj, iii, jjj, i2, r2, ks;
	double	wm;

	if (init != 0) {
		nmax = get1Di(Ns);
		nmin = get1Di(Ns);
		init = 0;
	}
	for (i = 0; i < Ns; i++) {
		nmin[i] = i - id;
		if (nmin[i] < 0)  nmin[i] = 0;
		nmax[i] = i + id;
		if (nmax[i] >= Ns)  nmax[i] = Ns - 1;
	}

	for (k = 0; k < Ks; k++) {		// reset iseed & jseed
		iseed[k] = -1;
		jseed[k] = -1;
	}
	for (i = 0; i < Ns; i++) {		// competition among k at each location
		for (j = 0; j < Ns; j++) {
			wmax[i][j] = 0.0;
			kwmax[i][j] = -1;
			for (k = 0; k < Ks; k++) {
				if (us[k][i][j] > wmax[i][j]) {
					wmax[i][j] = us[k][i][j];
					kwmax[i][j] = k;
				}
			}
		}
	}
	for (i = 0; i < Ns; i++) {		// find a winner in each hypercolumn
		for (j = 0; j < Ns; j++) {
			if ((wm = wmax[i][j]) > 0.0) {
				for (ii = nmin[i]; ii <= nmax[i]; ii++) {
					iii = ii - i;
					i2 = iii * iii;
					for (jj = nmin[j]; jj <= nmax[j]; jj++) {
						jjj = jj - j;
						r2 = i2 + jjj * jjj;
						if (r2 <= d2 && r2 != 0.0) {
							if (wmax[ii][jj] > wm) {		// defeated
								kwmax[i][j] = -1;
								goto skip;
							} else if (wmax[ii][jj] == wm) {	// tie
								wmax[ii][jj] = wm * 0.9999;
							}
						}
					}
				}
			}
			skip:	;
		}
	}
	for (i = 0; i < Ns; i++) {		// find a winner for each cell-plane
		for (j = 0; j < Ns; j++) {
			if ((ks = kwmax[i][j]) >= 0) {
				if ((iseed[ks] < 0) ||	// first appearance of a winner in ks
						(wmax[i][j] > wmax[iseed[ks]][jseed[ks]])) {
					iseed[ks] = i;		// replace the seed
					jseed[ks] = j;
				}
			}
		}
	}
}



void Sresini(			// response of seed-selecting plane
	double	**usini,	// (out) response of seed-selecting plane: usini[Ns][Ns]
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**asini,	// (in) exc. initial connections to S-cells: as[Na][Na]
	int		Nc,			// size of the C-layer: (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer: (ns*2) or (ns*2+1)
	int		Na,			// size of the connections (asini, as & cs): (na*2) or (na*2+1)
	int		NcMAX,		// max of Nc (and Ng) among all stages
	double	thrini		// threshold for cutting usini
	)
{
	STATIC int	init = 1;	// flag indicationg initial entry
	STATIC int	*nsmax;		// (work) upper bound of the loop: nsmax[NcMAX]
	STATIC int	*nsmin;		// (work) lower bound of the loop: nsmin[NcMAX]
	STATIC int	*namax;		// (work) upper bound of the conection: namax[NcMAX]
	int		ns2 = Ns - 1;
	int		nc2 = Nc - 1;
	int		na2 = Na - 1;
	int		na0 = (na2 - nc2 + ns2) / 2;	// na0=an-cn+sn; sn=ns2*0.5, cn=nc2*0.5, an=na2*0.5
	int		is, js, ic, jc, kc, ia, ja, ic1;
	double	uu, usimax, usicut;

	if (init != 0) {
		nsmax = get1Di(NcMAX);
		nsmin = get1Di(NcMAX);
		namax = get1Di(NcMAX);
		init = 0;
	}
	for (ic = 0; ic < Nc; ic++) {
		nsmax[ic] = ic + na0;					// = ic - cn + an + sn;
		if (nsmax[ic] > ns2)  nsmax[ic] = ns2;
		ic1 = nc2 - ic;
		nsmin[ic1] = ns2 - nsmax[ic];
		namax[ic1] = ic1 - nsmin[ic1] + na0;	// =(ic-cn) - (nsmin[ic]-sn) + an;
	}
	for (is = 0; is < Ns; is++) {
		for (js = 0; js < Ns; js++) {
			usini[is][js] = 0.0;
		}
	}
	for (kc = 0; kc < Kc; kc++) {
		for (ic = 0; ic < Nc; ic++) {
			for (jc = 0; jc < Nc; jc++) {
				if ((uu = uc[kc][ic][jc]) > 0.0) {
					for (is = nsmin[ic], ia = namax[ic]; is <= nsmax[ic]; is++, ia--) {
						for (js = nsmin[jc], ja = namax[jc]; js <= nsmax[jc]; js++, ja--) {
							usini[is][js] += asini[ia][ja] * uu;
						}
					}
				}
			}
		}
	}
	usimax = 0.0;
	for (is = 0; is < Ns; is++) {
		for (js = 0; js < Ns; js++) {
			if (usini[is][js] > usimax)  usimax = usini[is][js];
		}
	}
	usicut = usimax * thrini;
	for (is = 0; is < Ns; is++) {
		for (js = 0; js < Ns; js++) {
			usini[is][js] -= usicut;
			if (usini[is][js] < 0.0)  usini[is][js] = 0.0;
		}
	}
}



void seedini(			// select seed cells
	int		*iseed0,	// (out) x-coordinate of the k-th cell-plane: iseed0[1]
	int		*jseed0,	// (out) y-coordinate of the k-th cell-plane: jseed0[1]
	double	**usini,	// (in) response of S-cells: usini[Ns][Ns]
	int		Ns			// size of the S-layer
	)
{
	int		i, j;
	double	umax;

	umax = 0.0;
	iseed0[0] = -1;
	jseed0[0] = -1;
	for (i = 0; i < Ns; i++) {		// find maximum output
		for (j = 0; j < Ns; j++) {
			if (usini[i][j] > umax) {
				umax = usini[i][j];
				iseed0[0] = i;
				jseed0[0] = j;
			}
		}
	}
}



void lat(				// lateral inhibition from uin to uout
	double	**ut,		// (i/o) cell-plane to be inhibited: ut[N][N]
	double	**uf,		// (in) inhibiting cell-plane: uf[N][N]
	double	**cl,		// (in) connections for lateral inhibition: cl[Nl][Nl]
	int		N,			// size of the layer
	int		Nl			// size of connection cl:  Nl=nl*2+1
	)
						//  Nl=nl*2+1,  nl=Al
{
	int		i, j, ii, jj, iii, jjj;
	double	uu;
	int		nl = Nl / 2;		// Nl=nl*2+1

	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			if ((uu = ut[i][j]) > 0.0) {
				for (ii = -nl; ii <= nl; ii++) {
					iii = i + ii;
					if (iii >= 0 && iii < N) {
						for (jj = -nl; jj <= nl; jj++) {
							jjj = j + jj;
							if (jjj >= 0 && jjj < N) {
								uu -= cl[ii + nl][jj + nl] * uf[iii][jjj];
							}
						}
					}
				}
				if (uu < 0.0)  uu = 0.0;
				ut[i][j] = uu;
			}
		}
	}
}



int SorgL(				// self-organization of S-layer (new cell-plane if winner is not from correct category)
	double	****as,		// (out) exc. connections to S-cells: as[Ks][Kc][Na][Na]
	double	*bs,		// (out) inh. connections to S-cells: bs[Ks], new Ks
	int		*cat,		// (out) category of the S-plane: cat[Ks]
	int		Cat,		// (in) category of the training pattern
	double	***uc,		// (in) response of C-cells: uc[Kc][Nc][Nc]
	double	**cs,		// (in) spatial weight: cs[Na][Na]
	double	**asini,	// (in) exc. initial connections to S-cells: as[Na][Na]
	double	***us,		// (work) response of S-cells: us[Ks][Ns][Ns], old Ks
	int		Nc,			// size of the C-layer:  (nc*2) or (nc*2+1)
	int		Kc,			// number of cell-planes in C-layer
	int		Ns,			// size of the S-layer:  (ns*2) or (ns*2+1)
	int		Ks,			// number of cell-planes in S-layer, old value
	int		Na,			// size of the connections (as & cs): (na*2) or (na*2+1)
	int		Nasini,		// size of the connect. asini:  (nasini*2) or (nasini*2+1)
	double	qs,			// parameter determining the speed of learning
	double	theta,		// threshold of an S-cell
	int		Ksmax,		// maximum possible Ks
	int		KsMAX,		// max of KsMAX among all stages, because KcMAX <= KsMAX
	int		NcMAX,		// max of Nc (and Ng) among all stages
	int		NsMAX,		// max of Ns among all stages
	FILE	*fpT		// (in) file pointer for debugging.  If ==NULL, no debugging
	)
{
	STATIC int		init = 1;	// flag indicationg initial entry
	STATIC double	**usini;	// (work) response of seed-selecting plane: usini[Ns][Ns]
	int		ks, iseed0, jseed0, kseed0, k, i, j;
	double	usMax;		// max value of us[k][i][j]
	int		kwin;		// k that gives usMax

	if (init != 0) {
		usini = get2D(Ns, Ns);
		init = 0;
	}
	kseed0 = -1;
	if (Ks > 0) {
		ks = Ks;
		Sres(us, uc, as, cs, bs, Nc, Kc, Ns, Ks, Na, NcMAX, NsMAX, theta);
										// Sres, using old connections
		usMax = 0.0;
		kwin = -1;
		iseed0 = -1;
		jseed0 = -1;
		for (k = 0; k < Ks; k++) {
			for (i = 0; i < Ns; i++) {
				for (j = 0; j < Ns; j++) {
					if (us[k][i][j] > usMax) {
						usMax = us[k][i][j];
						kwin = k;
						iseed0 = i;
						jseed0 = j;
					}
				}
			}
		}
		if ((kwin >= 0) && (cat[kwin] == Cat)) {
			kseed0 = kwin;							// else kseed0 = -1 (unchanged)
		}

//TEST
		if (fpT != NULL) {
			if ((kwin >= 0) && (cat[kwin] != Cat)) {
				PRINTF("        Erroneous response:   recCat=%d   (Input Cat=%d)\n",
					cat[kwin], Cat);
				fprintf(fpT, "        Erroneous response:   recCat=%d   (Input Cat=%d)\n",
					cat[kwin], Cat);
			}
		}
//TEST

	}

	if (kseed0 < 0) {
		ks = Ks + 1;
		Sresini(usini, uc, asini, Nc, Kc, Ns, Nasini, NcMAX, 0.0);
		seedini(&iseed0, &jseed0, usini, Ns);	// location of seed = max.-output cell
		if (iseed0 == -1) {
			return Ks;
		}

		cat[Ks] = Cat;
		add4D2(as, Ks + 1, Kc, Ksmax);		// add Kc elements to pointer array **as
		kseed0 = Ks;						// cell-plane to be reinforced

//TEST
		if (fpT != NULL) {
			PRINTF("        New cell-plane created (next)\n");
			fprintf(fpT, "        New cell-plane created (next)\n");
		}
//TEST

	}

	delta(as, bs, uc, cs, Na, Nc, Kc, Ns, kseed0, iseed0, jseed0, KsMAX, qs);
				// if (0 <= kseed0 < Ks): reinforce connect. for old cell-planes
				// if (kseed0 = Ks): reinforce connect. for new cell-planes
				//                   memory for bs[] must be prepared beforehand
				// memory for as[] may be increased by delta()

	return  ks;
}



void rwcon(				// read/write connections.  Also get memores as and bs, in the read mode
	FILE	*fpC,		// (in) file pointer
	char	*fnameC,	// (in) file name
	double	*****as,	// (i/o) exc. connect. to S-cells: as[L][Ks][Kc][Na][Na]
	double	**bs,		// (i/o) inh. connect. to S-cells: bs[L][Ks]
	int		**cat0,		// (i/o) category of the cell-plane UsL: cat0[1][KsL]=cat[Ks]
	int		*Na,		// (in) size of connect. as
	int		*Kc,		// (i/o) no. of C-cell-planes: Kc[L+1]
	int		*Ks,		// (i/o) no. of S-cell-planes: Ks[L]
	int		*Mas,		// (i/o) no. of active as[L][Ks][Kc]: Mas[L]
	int		L,			// number of stages
	int		learn		// switch:  if (learn>0) write, else read.  If (learn<0) up to (-learn-1)
	)
{
	double	**ma;
	int		l, isize, Na2, ks, kc, m, Lread;

	if (learn > 0) {							// write
		fpC = fopen(fnameC, "wb");
		if (fpC == NULL) {
			PRINTF("file open error fpC\n");
			EXIT;
		}

		isize = fwrite(Kc, sizeof(int), L+1, fpC);
		if (isize != L+1)  goto error;

		isize = fwrite(Ks, sizeof(int), L, fpC);
		if (isize != L)  goto error;

		for (l = 0; l < L; l++) {
			Mas[l] = 0;
			for (ks = 0; ks < Ks[l]; ks++) {
				for (kc = 0; kc < Kc[l]; kc++) {
					if (as[l][ks][kc] != NULL) {
						Mas[l]++;
					}
				}
			}
		}
		isize = fwrite(Mas, sizeof(int), L, fpC);
		if (isize != L)  goto error;

		for (l = 0; l < L; l++) {
			Na2 = Na[l] * Na[l];
			for (ks = 0; ks < Ks[l]; ks++) {
				for (kc = 0; kc < Kc[l]; kc++) {
					if ((ma = as[l][ks][kc]) != NULL) {
						isize = fwrite(&ks, sizeof(int), 1, fpC);
						if (isize != 1)  goto error;
						isize = fwrite(&kc, sizeof(int), 1, fpC);
						if (isize != 1)  goto error;
						isize = fwrite(ma[0], sizeof(double), Na2, fpC);
						if (isize != Na2)  goto error;
					}
				}
			}

			isize = fwrite(bs[l], sizeof(double), Ks[l], fpC);
			if (isize != Ks[l])  goto error;
		}

		isize = fwrite(cat0[0], sizeof(int), Ks[L-1], fpC);
		if (isize != Ks[L-1])  goto error;
	} else if (learn != -1) {	// if (learn==0, or <=-2), read
		if (learn == 0) {		// if (learn==0), read all stages
			Lread = L;
		} else {				// if (learn<=-2), read up to (-lean-1)th stage
			Lread = -learn - 1;
		}

		fpC = fopen(fnameC, "rb");
		if (fpC == NULL) {
			PRINTF("file open error fpC\n");
			EXIT;
		}

		isize = fread(Kc, sizeof(int), L+1, fpC);
		if (isize != L+1)  goto error;

		isize = fread(Ks, sizeof(int), L, fpC);
		if (isize != L)  goto error;

		isize = fread(Mas, sizeof(int), L, fpC);
		if (isize != L)  goto error;

		for (l = 0; l < Lread; l++) {
			as[l] = geta4D1(Ks[l]);
			add4D2(as[l], Ks[l], Kc[l], Ks[l]);	// get Kc[l]*Ks[l] elements to pointer array
			Na2 = Na[l] * Na[l];
			for (ks = 0; ks < Ks[l]; ks++) {
				for (kc = 0; kc < Kc[l]; kc++) {
					as[l][ks][kc] = NULL;
				}
			}
			ma = get2D(Na[l] * Mas[l], Na[l]);			// get Mas arrays [Na][Na]
			for (m = 0; m < Mas[l]; m++) {
				isize = fread(&ks, sizeof(int), 1, fpC);
				if (isize != 1)  goto error;
				isize = fread(&kc, sizeof(int), 1, fpC);
				if (isize != 1)  goto error;
				as[l][ks][kc] = ma;
				isize = fread(ma[0], sizeof(double), Na2, fpC);
				if (isize != Na2)  goto error;
				ma += Na[l];
			}

			bs[l] = get1D(Ks[l]);
			isize = fread(bs[l], sizeof(double), Ks[l], fpC);
			if (isize != Ks[l])  goto error;
		}
		if (learn == 0) {			// Namely, Lread == L
			cat0[0] = get1Di(Ks[L-1]);
			isize = fread(cat0[0], sizeof(int), Ks[L-1], fpC);
			if (isize != Ks[L-1])  goto error;
		}
	}

	PRINTF("Scale of the network:\n");
	for (l = 0; l < L; l++) {
		PRINTF("    Stage %d:        Ks = %d,      Mas = %d\n", l+1, Ks[l], Mas[l]);
	}
	PRINTF("\n");
	fclose(fpC);
	return;

  error:	PRINTF("i/o error in fpC:       isize=%d\n", isize);
	fclose(fpC);
	return;
}



void printu(
	FILE	*fp,		// (in) file pointer
	double	***u,		// (in) layer to be printed: u[K][N][N]
	int		N,			// size of a cell-plane: (n*2+1) or (n*2)
	int		K,			// number of cell-planes
	char	title[],	// title to be printed
	int		l,			// stage number to be printed in the title. if (l<0) no print
	int		KMAX		// >0: maximum possible K, Initialize only, No printing.  ==0: print response
	)
{
	STATIC int		Kmax;
	STATIC double	*umax;		// max response of cell-planes: umax[KMAX]
	int		k, i, j, iuu, line, n, odd, j0;
	double	amax, anorm, uu;
	char	ch[2] = " ";

	if (KMAX > 0) {				// Initialize only, No printing.
		umax = get1D(KMAX);
		Kmax = KMAX;

	} else {					// Print response

		if (K > Kmax) {				// Unusual case, which should not occur!
			PRINTF("printu() should have been initialized with a larger KMAX!");
			umax = get1D(K);
			Kmax = K;
		}

		n = N / 2;		// N=(n*2+1) or (n*2)
		odd = N % 2;
		amax = 0.0;
		for (k = 0; k < K; k++) {
			umax[k] = 0.0;
			for (i = 0; i < N; i++) {
				for (j = 0; j < N; j++) {
					if (u[k][i][j] > umax[k]) {
						umax[k] = u[k][i][j];
					}
				}
			}
			if (umax[k] > amax)  amax = umax[k];
		}
		if (l >= 0) {
			fprintf(fp, "\n%s%d      max = %.9f        %d cell-planes\n", title, l, amax, K);
		} else {
			fprintf(fp, "\n%s      max = %.9f        %d cell-planes\n", title, amax, K);
		}
		anorm = (amax > 0.0) ?  10.0 / amax : 0.0;

		if (N == 1) {
			if (amax > 0.0) {
				fprintf(fp, "\n   ");
				for (k = 0; k < K; k++)  fprintf(fp, "%3d", k);
				fprintf(fp, "\n   ");
				for (k = 0; k < K; k++) {
					uu = u[k][0][0] * anorm;
					if (uu == 0.0)  ch[0] = ' ';
					else if (u[k][0][0] == amax)  ch[0] = 'X';
					else if (uu < 0.0)  ch[0] = '-';
					else if (uu < 0.01)  ch[0] = '_';
					else  ch[0] = '0' + ((int) uu);
					fprintf(fp, "  %1s", ch);
				}
				fprintf(fp, "\n");
			}
			return;
		}

		for (k = 0; k < K; k++) {
			if (umax[k] > 0.0) {
				if (K != 1) {
					if (l >= 0) {
						fprintf(fp, "\n %s%d(k=%d)  (max of this cell-plane) = %.9f\n", title, l, k, umax[k]);
					} else {
						fprintf(fp, "\n %s(k=%d)  (max of this cell-plane) = %.9f\n", title, k, umax[k]);
					}
				}

				fprintf(fp, "   ");
				for (j = 0; j < N; j++) {
					j0 = abs(j - n);
					if ((odd == 0) && (j >= n)) {
						j0++;
					}
					fprintf(fp, "%1d", j0 % 10);
				}
				fprintf(fp, "\n");

				fprintf(fp, "   ");
				for (j = 0; j < N; j++)  fprintf(fp, "*");
				fprintf(fp, "\n");

				for (i = 0; i < N; i++) {
					line = abs(i - n);
					if ((odd == 0) && (i >= n)) {
						line++;
					}
					fprintf(fp, "%2d*", line);
					for (j = 0; j < N; j++) {
						uu = u[k][i][j] * anorm;
						iuu = (int) uu;
						if (u[k][i][j] == amax)  iuu = 10;		// in case of rounding error

						if (uu == 0.0) {
							if ( (abs(j0) % 10 == 0) || (line % 10 == 0) )  ch[0] = ',';
							else  ch[0] = '.';
						} else if (uu < 0.0) {
							ch[0] = '-';
						} else if (uu < 0.01) {
							ch[0] = '_';
						} else if (iuu >= 10) {
							ch[0] = 'X';
						} else {
							ch[0] = '0' + iuu;
						}
						fprintf(fp, "%1s", ch);
					}
					fprintf(fp, "*%1d\n", line % 10);
				}

				fprintf(fp, "   ");
				for (j = 0; j < N; j++)  fprintf(fp, "*");
				fprintf(fp, "\n");

				fprintf(fp, "   ");
				for (j = 0; j < N; j++) {
					j0 = abs(j - n);
					if ((odd == 0) && (j >= n)) {
						j0++;
					}
					fprintf(fp, "%1d", j0 % 10);
				}
				fprintf(fp, "\n");
			}
		}
	}
}



void outH(				// hexadecinmal output of the response of a layer
	FILE	*fp,		// (in) file pointer
	double	***u,		// (in) layer to be printed: u[K][N][N]
	int		N,			// size of a cell-plane: (n*2) or (n*2+1)
	int		K,			// number of cell-planes
	char	title[],	// title to be printed
	int		l,			// stage number in the title. if (l==-1) no number; if (l==-2) no title
	int		KMAX		// >0: maximum possible K, Initialize only, No output.  ==0: write output
	)
{
	STATIC int		Kmax;
	STATIC double	*umax;		// max response of cell-planes: umax[KMAX]
	int		k, i, j, iuu, odd;
	double	amax, anorm, uu;

	if (KMAX > 0) {				// Initialize only, No output.
		umax = get1D(KMAX);
		Kmax = KMAX;
	} else {					// Write output

		if (K > Kmax) {				// Unusual case, which should not occur!
			PRINTF("outH() should have been initialized with a larger KMAX!");
			umax = get1D(K);
			Kmax = K;
		}

		if (l >= 0) {
			fprintf(fp, "/%s%1dplnumber %d def\n", title, l, K);
															// "/??plnumber K def"
			fprintf(fp, "/%s%1dpattern [\n", title, l);		// "/??pattern ["
		} else if (l == -1) {
			fprintf(fp, "/%splnumber %d def\n", title, K);
															// "/?plnumber K def"
			fprintf(fp, "/%spattern [\n", title);			// "/?pattern ["
		}

		odd = N % 2;
		amax = 0.0;
		for (k = 0; k < K; k++) {
			umax[k] = 0.0;
			for (i = 0; i < N; i++) {
				for (j = 0; j < N; j++) {
					if (u[k][i][j] > umax[k]) {
						umax[k] = u[k][i][j];
					}
				}
			}
			if (umax[k] > amax)  amax = umax[k];
		}
		anorm = (amax > 0.0) ?  16.0 / amax : 0.0;

		for (k = 0; k < K; k++) {
			if (umax[k] > 0.0) {
				fprintf(fp, "<\n");
				for (i = 0; i < N; i++) {
					for (j = 0; j < N; j++) {
						uu = u[k][i][j] * anorm;
						iuu = 15 - ((int) uu);
						if (iuu > 15) {
							iuu = 15;
						} else if (iuu < 0) {
							iuu = 0;
						}
						fprintf(fp, "%1X", iuu);
					}
					if (odd == 1) {
						fprintf(fp, "F\n");
					} else {
						fprintf(fp, "\n");
					}
				}
				fprintf(fp, ">\n");
			} else {
				fprintf(fp, "<FF>\n");
			}
		}
		if (l != -2) {
			fprintf(fp, "] def\n");
		}
	}
}



void outHhead(			// write header of hexadecinmal output file fpN
	FILE	*fpN,		// (in) file pointer
	char	fnameNh[],	// header of fnameN (fpNh)
	int		Npages		// no. of pages. >0: all data in one ps file.  =0: separate eps files.
	)
{
	static int	ipage = 0;
	FILE	*fpNh = NULL;	// file pointer (header of fpN file)
	char	linebuf[101];	// line buffer for copying data to fpN

	if (ipage == 0) {		// if ((ifpNall==0) || (ifpNall==1 && ipage==0))
	   // open file "fnameNh"
		if ((fpNh = fopen(fnameNh, "r")) == NULL) {
			PRINTF("file open error fpNh\n");
			EXIT;
		}
		fgets(linebuf, 101, fpNh);						// read 1 line and discard
		if (Npages > 0) {								// ps file
			fprintf(fpN, "%%!PS-Adobe-3.0\n");			// "%!PS-Adobe-3.0\n"
		} else {										// eps file
			fprintf(fpN, "%%!PS-Adobe-3.0 EPSF-3.0\n");	// "%!PS-Adobe-3.0 EPSF-3.0\n"
		}
		fgets(linebuf, 101, fpNh);		// copy 1 line (%%BoundingBox) from fpNh to fpN
			if (fputs(linebuf, fpN) != 0)  PRINTF("file write error fpN\n");
		if (Npages > 0) {							// ps file
			fprintf(fpN, "%%%%Pages: %d\n", Npages);	// "%%Pages: ?"
		}
		while (fgets(linebuf, 101, fpNh) != NULL) {		// copy fpNh to fpN
			if (fputs(linebuf, fpN) != 0)  PRINTF("file write error fpN\n");
		}
		fclose(fpNh);
	}
	if (Npages > 0) {								// ps file
		ipage++;
		fprintf(fpN, "\n%%%%Page: %d %d\n", ipage, ipage);	// "%%Page: ? ?"
	}
}



void printuadd(			// print sum of k1-th to k2-th cell-planes
	FILE	*fp,		// (in) file pointer
	double	***u,		// (in) layer to be printed: u[K][N][N]
	int		N,			// size of a cell-plane: (n*2) or (n*2+1)
	char	title[],	// title to be printed
	int		l,			// stage number to be printed in the title. if (l<0) no print
	int		k1,			// added cell-planes from k1
	int		k2			// added cell-planes to k2
	)
{
	int		k, i, j, iuu, line, n, odd, j0;
	double	amax, anorm, uu, uuu;
	char	ch[2] = " ";

	n = N / 2;		// N=(n*2+1) or (n*2)
	odd = N % 2;
	amax = 0.0;
	for (i = 0; i < N; i++) {
		for (j = 0; j < N; j++) {
			uu = 0.0;
			for (k = k1; k <= k2; k++) {
				uu += u[k][i][j];
			}
			if (amax < uu)  amax = uu;
		}
	}

	if (l >= 0) {
		fprintf(fp, "\n%s%d        Sum of the %d-th to the %d-th cell-planes      max = %.9f\n",
													title, l, k1, k2, amax);
	} else {
		fprintf(fp, "\n%s        Sum of the %d-th to the %d-th cell-planes      max = %.9f\n",
													title, k1, k2, amax);
	}

	if (amax > 0.0) {
		anorm = 10.0 / amax;
		fprintf(fp, "   ");
		for (j = 0; j < N; j++) {
			j0 = abs(j - n);
			if ((odd == 0) && (j >= n)) {
				j0++;
			}
			fprintf(fp, "%1d", j0 % 10);
		}
		fprintf(fp, "\n");

		fprintf(fp, "   ");
		for (j = 0; j < N; j++)  fprintf(fp, "*");
		fprintf(fp, "\n");

		for (i = 0; i < N; i++) {
			line = abs(i - n);
			if ((odd == 0) && (i >= n)) {
				line++;
			}
			fprintf(fp, "%2d*", line);
			for (j = 0; j < N; j++) {
				uuu = 0.0;
				for (k = k1; k <= k2; k++) {
					uuu += u[k][i][j];
				}
				uu = uuu * anorm;
				iuu = (int) uu;

				if (uuu == amax) {
					iuu = 10;		// in case of rounding error
				}
				if (uu == 0.0) {
					if ( (abs(j0) % 10 == 0) || (line % 10 == 0) )  ch[0] = ',';
					else  ch[0] = '.';
				} else if (uu < 0.0) {
					ch[0] = '-';
				} else if (uu < 0.01) {
					ch[0] = '_';
				} else if (iuu >= 10) {
					ch[0] = 'X';
				} else {
					ch[0] = '0' + iuu;
				}
				fprintf(fp, "%1s", ch);
			}
			fprintf(fp, "*%1d\n", line % 10);
		}

		fprintf(fp, "   ");
		for (j = 0; j < N; j++)  fprintf(fp, "*");
		fprintf(fp, "\n");

		fprintf(fp, "   ");
		for (j = 0; j < N; j++) {
			j0 = abs(j - n);
			if ((odd == 0) && (j >= n)) {
				j0++;
			}
			fprintf(fp, "%1d", j0 % 10);
		}
		fprintf(fp, "\n");
	}
}



void pat(				// get a pattern (from ETL1 file or memory)
	double	**u,		// (out) generated pattern: u[N][N]
	int		N,			// size of the layer:  N*N = 64*64 or 65*65
	double	sat,		// saturation factor
	int		t,			// sequence number of the pattern to be generated
	int		T0,			// starting no. of patterns to be read from file
	int		TT,			// total no. of patterns to be read from file
	int		TTmax,		// maximum possible TT
	char	fnameETL[]	// name of ETL1 file (fpETL)
	)
{
	int		Lpat = 2052;		// Logical record length / pat for ETL1
	STATIC int		init = 1;
	static int		T0read = -1;
	static int		TTread = 0;
	STATIC unsigned char	**pat;	// read buffer/binary data for patterns,  pat[TTmax][Lpat]
	static char		fnameOld[] = "                         ";

	FILE	*fpETL;		// file pointer (ETL1)
	int		i, j, j2, m, pp, isize, ipat;
	double	pmax, norm, uu, sat1;

   // get memory for pattern buffers
	if (init == 1) {
		pat = get2Dc(TTmax, Lpat);		// read buffer,  pat[TTmax][Lpat]
		init = 0;
	}

   // read patterns from ETL1
	if ( (strcmp(fnameETL, fnameOld) != 0) || (T0 != T0read) || (TT > TTread) ) {
		if ((fpETL = fopen(fnameETL, "rb")) == NULL) {
			PRINTF("file open error fpETL\n");
			EXIT;
		}
		if (fseek(fpETL, sizeof(char) * Lpat * T0, SEEK_SET) != 0) {
			PRINTF("seek error in reading ETL1\n");
			fclose(fpETL);
			EXIT;
		}
		for (ipat = 0; ipat < TT; ipat++) {
			isize = fread(pat[ipat], sizeof(unsigned char), Lpat, fpETL);
			if (isize != Lpat) {
				PRINTF("read error in ETL1 (TL):   isize=%d,   ipat=%d\n", isize, ipat + T0);
				fclose(fpETL);
				EXIT;
			}
		}
		fclose(fpETL);
		fpETL = NULL;
		T0read = T0;
		TTread = TT;
		strcpy(fnameOld, fnameETL);
	}

   // decode a pattern from binary data
	pmax = 0.0;
	for (j = 0; j < N; j++) {
		u[0][j] = 0.0;
	}
	m = 32;
	for (i = 1; i < 64; i++) {
		for (j = 0; j < 32; j++) {
			j2 = j * 2;
			pp = pat[t][m];
			u[i][j2] = (double) (pp / 16);
			if (u[i][j2] > pmax)  pmax = u[i][j2];
			u[i][j2 + 1] = (double) (pp % 16);
			if (u[i][j2 + 1] > pmax)  pmax = u[i][j2 + 1];
			m++;
		}
		for (j = 64; j < N; j++) {
			u[i][j] = 0.0;
		}
	}
	for (i = 64; i < N; i++) {
		for (j = 0; j < N; j++) {
			u[i][j] = 0.0;
		}
	}

	if (sat == 0.0) {					// no saturation
		norm = 1.0 / pmax;
		for (i = 0; i < N; i++) {
			for (j = 0; j < N; j++) {
				u[i][j] *= norm;
			}
		}
	} else {
		sat1 = sat + 1.0;
		norm = sat * pmax;
		for (i = 0; i < N; i++) {
			for (j = 0; j < N; j++) {
				uu = u[i][j];
				u[i][j] = sat1 * uu / (uu + norm);
			}
		}
	}
}



void score(				// make confusion matrix
	int		**mat,		// (i/o) confusion matrix: mat[KcL][KcL+1]
	int		recCat,		// result of recognition
	int		Cat,		// category of the input pattern
	int		KcL			// number of categories
	)
{
	static int	init = 1;
	int		i, j, res;
	double	amax;

	if (init == 1) {
		init = 0;
		for (i = 0; i < KcL; i++) {
			for (j = 0; j <= KcL; j++) {
				mat[i][j] = 0;
			}
		}
	}
	amax = 0.0;
	if (recCat < 0) {		// rejected
		res = KcL;
	} else {
		res = recCat;
	}
	mat[Cat][res]++;
}



void printsc(			// print confusion matrix
	FILE	*fp,		// (in) file pointer
	int		**mat,		// (in) confusion matrix: mat[KcL][KcL+1]
	int		KcL			// number of categories
	)
{
	int		i, j, sum, cor, Sum;
	double	rate;

	Sum = 0;
	cor = 0;
	fprintf(fp, "\nConfusion Matrix\n");
	fprintf(fp, "   out=");
	for (j = 0; j < KcL; j++) {
		fprintf(fp, "%5d", j);
	}
	fprintf(fp, "  reject  Rec.rate\n");
	for (i = 0; i < KcL; i++) {
		sum = 0;
		fprintf(fp, " in=%d  ", i);
		for (j = 0; j <= KcL; j++) {
			fprintf(fp, "%5d", mat[i][j]);
			sum += mat[i][j];
		}
		rate = (sum > 0)  ?  (mat[i][i] * 100.0) / sum  :  0.0;
		fprintf(fp, "     %7.3f\n", rate);
		Sum += sum;
		cor += mat[i][i];
	}
	rate = (Sum > 0)  ?  (cor * 100.0) / Sum  :  0.0;
	fprintf(fp, "   Recognition rate = %7.3f\n\n\n", rate);
	PRINTF("   Recognition rate = %7.3f\n\n\n", rate);
}



void ErrOut(			// hexadecinal output of erroneously recognized patterns
	double	***u,		// (in) layer to be printed: u[1][N][N]
	int		N,			// size of a cell-plane: (n*2) or (n*2+1)
	char	fnameE[],	// file name (fpE), erroneously recognized patterns
	char	fnameEh[],	// header of fnameE (fpEh)
	int		iend		// if (iend==1), no more output (close file)
	)
{
	static int	ipat = 0;
	static int	ipage = 0;
	static FILE	*fpE = NULL;		// file pointer (erroneously recognized patterns, hexadecimal)
	static FILE	*fpEh = NULL;		// file pointer (header of fpE file)
	int		Mpat = 10 * 10;			// maximum number of patterns in a page
	int		i, mpat, mpat1, mpat10;
	char	linebuf[101];			// line buffer for copying data to fpE

	if (iend != 1) {				// Not the end
		if (ipat == 0) {
			if (ipage == 0) {		// write header
			   // open file "fnameE" (hexadecimal output, erroneously recognized patterns)
				if ((fpE = fopen(fnameE, "w")) == NULL) {
					PRINTF("file open error fpN\n");
					EXIT;
				}
			   // copy fpEh (header of fpE) to fpE
				if ((fpEh = fopen(fnameEh, "r")) == NULL) {		// open file "fnameEh"
					PRINTF("file open error fpEh\n");
					EXIT;
				}
				while (fgets(linebuf, 101, fpEh) != NULL) {		// copy fpEh to fpE
					if (fputs(linebuf, fpE) != 0)  PRINTF("file write error fpE\n");
				}
				fclose(fpEh);
			}
			ipage++;
			fprintf(fpE, "\n%%%%Page: %d %d\n", ipage, ipage);	// "%%Page: ? ?"
			fprintf(fpE, "/pattern [\n");				// "/pattern ["
		}
		outH(fpE, u, N, 1, " ", -2, 0);
		ipat++;
		if (ipat == Mpat) {
			fprintf(fpE, "] def\n");					// "] def"
			fprintf(fpE, "/Npat %d def\n", ipat);		// "/Npat ? def"
			fprintf(fpE, "display\n");					// "display"
			fprintf(fpE, "showpage\n");					// "showpage"
			ipat = 0;
		}
	} else {						// (iend == 1)
		if ((ipage != 0) || (ipat != 0)) {		// at least one pattern has been written
			if (ipat < Mpat) {
				mpat = Mpat - ipat;
				mpat1 = mpat % 10;
				mpat10 = mpat / 10;
				if (mpat1 != 0) {
					for (i = 0; i < mpat1; i++) {
						fprintf(fpE, "<FF>");
					}
					fprintf(fpE, "\n");
				}
				for (i = 0; i < mpat10; i++) {
					fprintf(fpE, "<FF><FF><FF><FF><FF><FF><FF><FF><FF><FF>\n");
				}
				fprintf(fpE, "] def\n");						// "] def"
				fprintf(fpE, "/Npat %d def\n", ipat);		// "/Npat ? def"
				fprintf(fpE, "display\n");						// "display"
				fprintf(fpE, "showpage\n");						// "showpage"
			}
			if (fpE != NULL)  fclose(fpE);
		}
	}
}



// ############################################################################

