// DLM.C
// This version is for the NSL-book

// Author: (c) Laurenz Wiskott, wiskott@salk.edu, 96-03-11 ;
// read DLMreadme ;

#include "nsl_include.h"

#include "DLMwin.h"
#include "DLMdisplay.h"
#include "DLMSimilarity.h"
#include "DLM.h"

static const int smallLayer  = 0;
static const int smallPatches = 1;
          
DLM::DLM() : NslModel("DLM"),
	similarity("similarity",this),
	recognition("recognition",this),
	correlation("correlation",this),
	w12("w12",this),
	w21("w21",this),
	h1("h1",this),
	h2("h2",this),
	a1("a1",this),
	a2("a2",this),
	display("display",this),

/* For the sake of efficiency some operations such as the modification
   of the dynamic links or the display of several layers takes place
   only after loops single iterations. The correlation for the
   modification is accumulated in each simulation step. */
   
	loops("loops",this), 	// number of cheap simulation steps ;

	gallerySize("gallerySize",this),   // actual size of the simmulated gallery ;

/* Since not all layers 2 can be displayed, only one layer 2 is copied
   into respective 'show' layers. You can either select one special
   layer to be displayed indicated by the variable 'preferredLayer' or
   if this variable does not have a valid value, the layer with the
   highest activity is displayed. */

	preferredLayer("preferredLayer",this),        // number of the layer 2 to show ;
                             // a value of -1 indicates no preferrences ;
	showLayer("showLayer",this),             // actual number of the layer 2 to show ;
                             // depending on preferredLayer and the ;
                             // activity of the layers ;

/* For objects and models there are stored graphs with the precomputed
   Gabor-wavelet transform coefficients.  For the models they are
   taken from a grid of 10 x 10 nodes centered on the faces. For the
   objects the grids cover the whole image plane with 16 x 17 nodes.
   From these stored graphs a subgraph can be selected.  From the 16 x
   17 graph for example a 12 x 12 subgraph will automatically be
   selected if the size of layer 1 is 12 x 12.  In addition one can
   choose the location of the subgraph by the offsets Si1offset, ... ,
   and Sj2offset.  They are given as integer numbers behind the model
   names in the gallery files.  */

	workOnAverage("workOnAverage",this),  // workOnAverage indicates whether to simulate ... ;
                      // .. currently in average instead of the model gallery ;
	avTimeLimit("avTimeLimit",this)    // time for attention dynamics on average model only ;
{
	nsl_string _str;
	for (int i=0; i<gallerySizeMax; i++) {
		strcpy(_str,"shSum");
		sprintf(_str,"%s(%d)",_str,i);
		shSum[i].initialize(_str,this);
		strcpy(_str,"skipModel");
		sprintf(_str,"%s(%d)",_str,i);
		skipModel[i].initialize(_str,this);
	}

	memAlloc();
	makeConn();
}

void DLM::memAlloc()
{
	int frame;      // width of the frame around layer 1
	int i1max; 		// size of the image layer 1, including frame;
	int j1max;
	int i2max;      // size of the model layers 2 ;
	int j2max;
	int i1Rmax;     // size of the projection patches .. ;
	int j1Rmax;     // ... from model layer 2 to image layer 1 ;

	if (smallLayer == 1) {
		frame = 0;          // width of the frame around layer 1
		i1max = 10+2*frame; // size of the image layer 1, including frame;
		j1max = 10+2*frame;
	} else {
		frame = 2;
		i1max = 17+2*frame;
		j1max = 16+2*frame;
	}

	if (smallPatches == 1) {
		i1Rmax= 8;          // size of the projection patches .. ;
		j1Rmax= 8;          // ... from model layer 2 to image layer 1 ;
	} else {
		i1Rmax= i1max-2*frame;
		j1Rmax= j1max-2*frame;
	}

	i2max = 10;         // size of the model layers 2 ;
	j2max = 10;
	
	similarity.memAlloc(frame,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax);
	w12.memAlloc(frame,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax);
	w21.memAlloc(frame,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax);
	correlation.memAlloc(frame,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax);
	a1.memAlloc(i1max,j1max);
	a2.memAlloc(i2max,j2max);
	h1.memAlloc(i1max,j1max);
	h2.memAlloc(i2max,j2max);
	display.memAlloc(frame,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax);
}
void DLM::makeConn()
{
	nslConnect(w12.hInput,h1.hInput);
	nslConnect(a1.sa,h1.sa);
	nslConnect(a2.sa,h2.sa);
	nslConnect(h1.sh,w21.sh);
	nslConnect(h1.sh,correlation.sh1);
	nslConnect(h1.sh,a1.sh);
	nslConnect(h2.shMax,a2.sh);
	for (int i=0; i<gallerySizeMax; i++) {
		nslConnect(similarity.sim[i],w12.sim[i]);
		nslConnect(similarity.sim[i],w21.sim[i]);
		nslConnect(w21.hInput[i],h2.hInput[i]);
		nslConnect(h2.shSum[i],recognition.shSum[i]);
		nslConnect(h2.shSum[i],shSum[i]); // DLM
		nslConnect(h2.sh[i],w12.sh[i]);
		nslConnect(h2.sh[i],correlation.sh2[i]);
		nslConnect(correlation.correlSum[i],w12.correlSum[i]);
		nslConnect(correlation.correlSum[i],w21.correlSum[i]);
		nslConnect(recognition.skipModel[i],w12.skipModel[i]);
		nslConnect(recognition.skipModel[i],w21.skipModel[i]);
		nslConnect(recognition.skipModel[i],correlation.skipModel[i]);
		nslConnect(recognition.skipModel[i],h2.skipModel[i]);
		nslConnect(recognition.skipModel[i],skipModel[i]);
	}
}
void DLM::initSys()
{
  if (gallerySize>=gallerySizeMax) {
    	NSLoutput("WARNING : gallerySize(",gallerySize.elem());
      	NSLoutput(")>=gallerySizeMax(",gallerySizeMax);
		NSLoutput("),\n");
    	NSLoutput("          set gallerySize=gallerySizeMax !!\n");
    	gallerySize=gallerySizeMax;
  }

  loop = (int)loops.elem();

	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}

  showLayer = 0;

  for (int model=0; model<=gallerySize; model++)
    skipModel[model] = 0;
}
void DLM::initRun()
{
	loop = 0;
}

void DLM::simRun()
{
	NSLoutput(":");

	if (loop<=0) {
		loop=(int)(loops.elem());

		// ---- after 'avTimeLimit' switch from average to model simulation -- ;

		if (nslSystem.getSimTime()>avTimeLimit)
			workOnAverage = 0;

		if (workOnAverage == 1) {
			modelLoIndex=0; 
			modelHiIndex=0;
		} else {
			modelLoIndex=1; 
			modelHiIndex=int(gallerySize.elem());
		}

		if (workOnAverage == 1)
			showLayer = 0;
		else if ((preferredLayer>0)&&(preferredLayer<=gallerySize))
			// if a valid preferredLayer is given take it, ... ;
			showLayer = preferredLayer;
		else {
			// ... otherwise take the most active layer ;
			showLayer = 1;
			for (int model=modelLoIndex; model<=modelHiIndex; model++) {
				if (skipModel[model].get_data() == 1) continue;
				if (shSum[model]>shSum[int(showLayer.elem())])
					showLayer = model;
			}
		}

		// ---- print the current recognition status: ----------------------- ;
		// ---- time, showed layer, recognition values ;

		NSLoutput("\nt = ",nslSystem.getSimTime());
		NSLoutput("\nshowLayer = ",showLayer.elem());

	} // if (loop<=0) ;

	loop--;
}

// W

DlmW::DlmW(nsl_string str,NslModule* parent) : 
	DlmCommon(str,parent), 
	loops("loops",this), 	// number of cheap simulation steps ;
	workOnAverage("workOnAverage",this),  
	avTimeLimit("avTimeLimit",this),
	gallerySize("gallerySize",this),  
	lambda_W("lambda_W",this)        
{
	nsl_string _str;
	for (int i=0; i<gallerySizeMax; i++) {
		strcpy(_str,"sim");
		sprintf(_str,"%s(%d)",_str,i);
		sim[i].initialize(_str,this);
		strcpy(_str,"w");
		sprintf(_str,"%s(%d)",_str,i);
		w[i].initialize(_str,this);
		strcpy(_str,"wSums");
		sprintf(_str,"%s(%d)",_str,i);
		wSums[i].initialize(_str,this);
		strcpy(_str,"correlSum");
		sprintf(_str,"%s(%d)",_str,i);
		correlSum[i].initialize(_str,this);
		strcpy(_str,"skipModel");
		sprintf(_str,"%s(%d)",_str,i);
		skipModel[i].initialize(_str,this);
	}
}

// W12

DlmW12::DlmW12(nsl_string str,NslModule* parent) : 
	DlmW(str,parent), 
    normFactor("normFactor",this), 
	hInput("hInput",this)
{
	nsl_string _str;
	for (int i=0; i<gallerySizeMax; i++) {
		strcpy(_str,"sh");
		sprintf(_str,"%s(%d)",_str,i);
		sh[i].initialize(_str,this);
	}
}

void DlmW12::memAlloc(int _frame, int _i1max,int _j1max,int _i2max,int _j2max,int _i1Rmax,int _j1Rmax)
{
	frame = _frame;
	i1max = _i1max;
	j1max = _j1max;
	i2max = _i2max;
	j2max = _j2max;
	i1Rmax = _i1Rmax;
	j1Rmax = _j1Rmax;
	
	normFactor.memAlloc(i1max,j1max);
	hInput.memAlloc(i1max,j1max);

	for (int i=0; i<gallerySizeMax;i++) {
		sim[i].memAlloc(i2max,j2max,i1Rmax,j1Rmax);
		w[i].memAlloc(i1Rmax,j1Rmax,i2max,j2max);
		sh[i].memAlloc(i2max,j2max);
		correlSum[i].memAlloc(i2max,j2max,i1Rmax,j1Rmax);
	}
}

void DlmW12::initRun()
{
	int model,i2,ij2,j2,i1R,i1,j1,j1R,ij1R;
	
	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}

	if (gallerySize>0) {
		// ---- copy 's1' to 'W21' and 'W12' ;
		for (model=1; model<=gallerySize; model++) {
			for (i2=0; i2<i2max; i2++)
				for (j2=0; j2<j2max; j2++)
					for (i1R=0; i1R<i1Rmax; i1R++)
						for (j1R=0; j1R<j1Rmax; j1R++)
							w[model][i1R][j1R][i2][j2] = sim[model][i2][j2][i1R][j1R];
    }

    float stmp;
    if (gallerySize>0)
		for (i2=0; i2<i2max; i2++)
			for (j2=0; j2<j2max; j2++)
				for (i1R=0; i1R<i1Rmax; i1R++)
					for (j1R=0; j1R<j1Rmax; j1R++) {
	  					stmp = 0;
	  					for (model=1; model<=gallerySize; model++)
	    					stmp = NSLmax(stmp,sim[model][i2][j2][i1R][j1R]);
						w[0][i1R][j1R][i2][j2] = stmp;
	}
    // ---- initialize W21Sums ;
    for (model=0; model<=gallerySize; model++)
		wSums[model] = w[model].sum()/(i1Rmax*j1Rmax*i2max*j2max);
  }
}
void DlmW12::simRun()
{
	int model,i2,j2,i1R,i1,j1,j1R;

	hInput=0;
	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model].get_data() == 1) continue;
		for (i2=0; i2<i2max; i2++)
			for (j2=0; j2<j2max; j2++) {
				if (sh[model][i2][j2]<=0.) continue;
				for (i1R=0, i1=i1Index(i2, 0,frame,i1max,i1Rmax,i2max); i1R<i1Rmax; i1R++, i1++)
					for (j1R=0, j1=j1Index(j2, 0,frame,j1max,j1Rmax,j2max); j1R<j1Rmax; j1R++, j1++) 
						// 1 
						hInput[i1][j1] = NSLmax(hInput[i1][j1],
							w[model][i1R][j1R][i2][j2]*sh[model][i2][j2]);
			}
	}

	if (loop<=0) {
		loop=(int)(loops.elem());

		if (workOnAverage == 0) {
			if (lambda_W!=0) {
			// ---- modify weights 'W21' and 'W12' ------------------------- ;
				normFactor = 1.;
      
				for (model=modelLoIndex; model<=modelHiIndex; model++) {
					if (skipModel[model].get_data() == 1) continue;
					for (i2=0; i2<i2max; i2++)
						for (j2=0; j2<j2max; j2++)
							for (i1R=0, i1=i1Index(i2,0,frame,i1max,i1Rmax,i2max); i1R<i1Rmax; i1R++, i1++)
								for (j1R=0, j1=j1Index(j2,0,frame,j1max,j1Rmax,j2max); j1R<j1Rmax; j1R++, j1++) {	
									// --- let weights grow ;
									// 16 
									w[model][i1R][j1R][i2][j2] += nslSystem.getSimDelta()
										*w[model][i1R][j1R][i2][j2]
										*lambda_W.elem()*correlSum[model][i2][j2][i1R][j1R];
		
									// ---- determine normalization factors ;
									if (w[model][i1R][j1R][i2][j2]>sim[model][i2][j2][i1R][j1R]) {
									// 17 
										normFactor[i1][j1] = NSLmin(normFactor[i1][j1],
											sim[model][i2][j2][i1R][j1R]/w[model][i1R][j1R][i2][j2]);
									}
								} // for (j1R=...) ;
				} // for (model=...) ;
      
				// ---- normalize weights ------------------------------------------ ;
      
				for (i2=0; i2<i2max; i2++)
					for (j2=0; j2<j2max; j2++)
						for (i1R=0, i1=i1Index(i2,0,frame,i1max,i1Rmax,i2max); i1R<i1Rmax; i1R++, i1++)
							for (j1R=0, j1=j1Index(j2,0,frame,j1max,j1Rmax,j2max); j1R<j1Rmax; j1R++, j1++) {
								for (model=modelLoIndex; model<=modelHiIndex; model++) {
									// 20 
									w[model][i1R][j1R][i2][j2] *= normFactor[i1][j1];
								}
							}
      
			} // if (lambda_W!=0) ;
      
			// ---- compute normed total weights of the nslConnectivities ---------- ;

			/* 'W21' does not vary for layers not simulated any more
			(skipModel==1) while 'W12' does.  The reason is that there is
			no competition between links going to different models while
			there is competition between links converging to layer 1 from
			different models.  That is the reason why the computation of
			'W21Sums' comes after the continue statement. */

			for (model=modelLoIndex; model<=modelHiIndex; model++) {
				wSums[model] = w[model].sum()/(i1Rmax*j1Rmax*i2max*j2max);
				if (skipModel[model].get_data() == 1) continue;
			}  
		} // if (!workOnAverage) ;
    
		if (nslSystem.getSimTime()>avTimeLimit)
			workOnAverage = 0;

		if (workOnAverage == 1) {
			modelLoIndex=0; 
			modelHiIndex=0;
		} else {
			modelLoIndex=1; 
			modelHiIndex=int(gallerySize.elem());
		}
	} // if (loop<=0)

	loop--;
}

// W21

DlmW21::DlmW21(nsl_string str,NslModule* parent) : 
	DlmW(str,parent), 
	sh("sh",this)        
{
	nsl_string _str;
	for (int i=0; i<gallerySizeMax; i++) {
		strcpy(_str,"hInput");
		sprintf(_str,"%s(%d)",_str,i);
		hInput[i].initialize(_str,this);
		strcpy(_str,"normFactor");
		sprintf(_str,"%s(%d)",_str,i);
		normFactor[i].initialize(_str,this);
	}
}

void DlmW21::memAlloc(int _frame, int _i1max,int _j1max,int _i2max,int _j2max,int _i1Rmax,int _j1Rmax)
{
	frame = _frame;
	i1max = _i1max;
	j1max = _j1max;
	i2max = _i2max;
	j2max = _j2max;
	i1Rmax = _i1Rmax;
	j1Rmax = _j1Rmax;
	
	sh.memAlloc(i1max,j1max);
	wShow.memAlloc(i2max,j2max,i1Rmax,j1Rmax);

	for (int i=0; i<gallerySizeMax; i++) {
		sim[i].memAlloc(i2max,j2max,i1Rmax,j1Rmax);
		w[i].memAlloc(i2max,j2max,i1Rmax,j1Rmax);
		hInput[i].memAlloc(i2max,j2max);
		correlSum[i].memAlloc(i2max,j2max,i1Rmax,j1Rmax);
		normFactor[i].memAlloc(i2max,j2max);
	}
}

void DlmW21::initRun()
{
	int model,i2,ij2,j2,i1R,i1,j1,j1R,ij1R;
	
	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}

	if (gallerySize>0) {
		for (model=1; model<=gallerySize; model++) {
			w[model] = sim[model]; 
		}

    // ---- initialize 'W21', for average layer ;
		float stmp;
		if (gallerySize>0)
			for (i2=0; i2<i2max; i2++)
				for (j2=0; j2<j2max; j2++)
					for (i1R=0; i1R<i1Rmax; i1R++)
						for (j1R=0; j1R<j1Rmax; j1R++) {
	  						stmp = 0;
	  						for (model=1; model<=gallerySize; model++)
	    						stmp = NSLmax(stmp,sim[model][i2][j2][i1R][j1R]);
	  						w[0][i2][j2][i1R][j1R] = stmp;
						}
    // ---- initialize W21Sums ;
		for (model=0; model<=gallerySize; model++) 
			wSums[model] = w[model].sum()/(i1Rmax*j1Rmax*i2max*j2max);
    
	}
}
void DlmW21::simRun()
{
	int model,i2,j2,i1R,i1,j1,j1R;

  // ---- compute layer input 'hInput' propagated via 'W21' ---------- ;

	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model].get_data() == 1) continue;
    		hInput[model]=0;
	}

	for (i1=frame; i1<i1max-frame; i1++)
		for (j1=frame; j1<j1max-frame; j1++) {
			if (sh[i1][j1]<=0.) continue;
			for (model=modelLoIndex; model<=modelHiIndex; model++) {
				if (skipModel[model].get_data() == 1) continue;
				for (i2=0; i2<i2max; i2++) {
					i1R = i1RIndex(i2,i1,frame,i1max,i1Rmax,i2max);
					if ((i1R<0)||(i1R>=i1Rmax)) continue;
					for (j2=0; j2<j2max; j2++) {
						j1R = j1RIndex(j2,j1,frame,j1max,j1Rmax,j2max);
						if ((j1R<0)||(j1R>=j1Rmax)) continue;
						// 2 
						hInput[model][i2][j2] = NSLmax(hInput[model][i2][j2],
							w[model][i2][j2][i1R][j1R]*sh[i1][j1]);
					}
				}
			}
		}

	if (loop<=0) {
		loop=(int)(loops.elem());

		if (workOnAverage == 0) {
			if (lambda_W!=0) {

				// ---- modify weights 'W21' and 'W12' ------------------------- ;
				for (model=modelLoIndex; model<=modelHiIndex; model++) {
					if (skipModel[model].get_data() == 1) continue;
					normFactor[model] = 1.; // *(normFactor[model]) = 1.;	
					for (i2=0; i2<i2max; i2++)
						for (j2=0; j2<j2max; j2++)
							for (i1R=0, i1=i1Index(i2,0,frame,i1max,i1Rmax,i2max); i1R<i1Rmax; i1R++, i1++)
								for (j1R=0, j1=j1Index(j2,0,frame,j1max,j1Rmax,j2max); j1R<j1Rmax; j1R++, j1++) {
									// --- let weights grow ;
									// 15 
									w[model][i2][j2][i1R][j1R] += nslSystem.getSimDelta()
										*w[model][i2][j2][i1R][j1R]
										*lambda_W.elem()*correlSum[model][i2][j2][i1R][j1R];
		
									// ---- determine normalization factors ;
									if (w[model][i2][j2][i1R][j1R]>sim[model][i2][j2][i1R][j1R]) {
										// 18 
										normFactor[model][i2][j2] = NSLmin(normFactor[model][i2][j2],
											sim[model][i2][j2][i1R][j1R]/w[model][i2][j2][i1R][j1R]);
									}
								} // for (j1R=...) ;
				} // for (model=...) ;
      
			// ---- normalize weights ------------------------------------------ ;
      
				for (i2=0; i2<i2max; i2++)
					for (j2=0; j2<j2max; j2++)
						for (i1R=0, i1=i1Index(i2,0,frame,i1max,i1Rmax,i2max); i1R<i1Rmax; i1R++, i1++)
							for (j1R=0, j1=j1Index(j2,0,frame,j1max,j1Rmax,j2max); j1R<j1Rmax; j1R++, j1++) {
								for (model=modelLoIndex; model<=modelHiIndex; model++) {
									// 19 
									w[model][i2][j2][i1R][j1R] *= normFactor[model][i2][j2]; 
								}
							}
			} // if (lambda_W!=0) ;
      
			// ---- compute normed total weights of the nslConnectivities ---------- ;

			/* 'W21' does not vary for layers not simulated any more
			(skipModel==1) while 'W12' does.  The reason is that there is
			no competition between links going to different models while
			there is competition between links converging to layer 1 from
			different models.  That is the reason why the computation of
			'W21Sums' comes after the continue statement. */

			for (model=modelLoIndex; model<=modelHiIndex; model++) {
				if (skipModel[model].get_data() == 1) continue;
				wSums[model] = w[model].sum()/(i1Rmax*j1Rmax*i2max*j2max);
			}
		} // if (!workOnAverage) ;

		if (nslSystem.getSimTime()>avTimeLimit)
			workOnAverage = 0;

		if (workOnAverage == 1) {
			modelLoIndex=0; 
			modelHiIndex=0;
		} else {
			modelLoIndex=1; 
			modelHiIndex=int(gallerySize.elem());
		}
	} // if (loop<=0) 

	loop--;
}

// Correlation

DlmCorrelation::DlmCorrelation(nsl_string str,NslModule* parent) : 
	DlmCommon(str,parent), 
	loops("loops",this), 	// number of cheap simulation steps ;
	correlLIShow("correl21LIShow",this),
	showLayer("showLayer",this),        
	workOnAverage("workOnAverage",this),  
	avTimeLimit("avTimeLimit",this),
	gallerySize("gallerySize",this),  
	sh1("sh",this),	     // maximum activity in the layers 2
	lambda_c("lambda_c",this),               // rate for the leaky correlation integrater ;
	temp()
{
	nsl_string _str;
	for (int i=0; i<gallerySizeMax; i++) {
		strcpy(_str,"sh");
		sprintf(_str,"%s(%d)",_str,i);
		sh2[i].initialize(_str,this);
		strcpy(_str,"correlSum");
		sprintf(_str,"%s(%d)",_str,i);
		correlSum[i].initialize(_str,this);
		strcpy(_str,"correlLI");
		sprintf(_str,"%s(%d)",_str,i);
		correlLI[i].initialize(_str,this);
		strcpy(_str,"skipModel");
		sprintf(_str,"%s(%d)",_str,i);
		skipModel[i].initialize(_str,this);
	}
}

void DlmCorrelation::memAlloc(int _frame, int _i1max,int _j1max,int _i2max,int _j2max,int _i1Rmax,int _j1Rmax)
{
	frame = _frame;
	i1max = _i1max;
	j1max = _j1max;
	i2max = _i2max;
	j2max = _j2max;
	i1Rmax = _i1Rmax;
	j1Rmax = _j1Rmax;

	sh1.memAlloc(i1max,j1max);
	temp.memAlloc(i2max,j2max,i1Rmax,j1Rmax);
	correlLIShow.memAlloc(i2max,j2max,i1Rmax,j1Rmax);

	for (int i=0; i<gallerySizeMax;i++){
		sh2[i].memAlloc(i2max,j2max);
		correlSum[i].memAlloc(i2max,j2max,i1Rmax,j1Rmax);
		correlLI[i].memAlloc(i2max,j2max,i1Rmax,j1Rmax);
	}
}
void DlmCorrelation::initRun()
{
	int model;
	
	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}

  /* The correlation matrices 'correl21LI' are initialized with a very
     small value, so that the net display can work in the
     beginning. */

	for (model=0; model<=gallerySize; model++)
		correlLI[model] = 0.00001;

	for (model=0; model<=gallerySize; model++)
		correlSum[model] = 0;
}
void DlmCorrelation::simRun()
{
	int model,i2,j2,i1R,i1,j1,j1R;
	
	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model] == 1) continue;
		for (i2=0; i2<i2max; i2++)
			for (j2=0; j2<j2max; j2++) {
				if (sh2[model][i2][j2]<=0.) continue;
				for (i1R=0, i1=i1Index(i2, 0,frame,i1max,i1Rmax,i2max); i1R<i1Rmax; i1R++, i1++)
					for (j1R=0, j1=j1Index(j2, 0,frame,j1max,j1Rmax,j2max); j1R<j1Rmax; j1R++, j1++) {
						// 14 
						correlSum[model][i2][j2][i1R][j1R] += sh1[i1][j1]*sh2[model][i2][j2];
					}
			}
	}

	// ---- simulate one time step on 'correl21LI' ------------------------ ;
	// ---- needs valid 'correl21Sum' ;

	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model].get_data() == 1) continue;
//				correlLI[model] = correlLI[model]+nslSystem.getSimDelta()*lambda_c.elem()
//				*(correlSum[model] -loops.elem()* correlLI[model]);

		NSLmult_equal(temp,-loops.elem(),correlLI[model]);
		NSLadd_plus(temp,correlSum[model]);
		NSLmult_equal(temp,lambda_c.elem(),temp);
		NSLmult_plus(correlLI[model],nslSystem.getSimDelta(),temp);
	}

	if (loop<=0) {
		loop=(int)(loops.elem());

		for (model=modelLoIndex; model<=modelHiIndex; model++) {
			if (skipModel[model].get_data() == 1) continue;
			// 21 
			correlSum[model] = 0;
		}   

		if (nslSystem.getSimTime()>avTimeLimit)
			workOnAverage = 0;

		if (workOnAverage == 1) {
			modelLoIndex=0; 
			modelHiIndex=0;
		} else {
			modelLoIndex=1; 
			modelHiIndex=int(gallerySize.elem());
		}

	} // if (loop<=0)
	loop--;
}

// Attention
   
DlmAttention::DlmAttention(nsl_string str,NslModule* parent) : 
	DlmAHCommon(str,parent),
	a("a",this),       // attention blob on layer 1 ;
	sa("sa",this),      // cell activities ;
	aTransE("aTransE",this), // excitatoryly transfered signal ;
	sh("sh",this),	     // maximum activity in the layers 2
	lambda_a("lambda_a",this),       // time constant for the attention dynamics ;
	beta_a("beta_a",this),         // strength of global inhibition for attention blob ;
	kappa_ah("kappa_ah",this),       // effect of running blob on attention blob ;
	alpha_N("alpha_N",this)       // parameter for attention blob initialization ;
{
}

void DlmAttention::memAlloc(int imax,int jmax)
{
	a.memAlloc(imax,jmax);
	sa.memAlloc(imax,jmax);
	aTransE.memAlloc(imax,jmax);
	sh.memAlloc(imax,jmax);
	temp.memAlloc(imax,jmax);
}

void DlmAttention::initRun()
{
	if (attention == 1) {
		a = alpha_N; // a2: alpha_N = 0;
		computeOutputFunc(sa,a,rho); // a2: sa = 0;
		initGauss(g,sigma_g);
	} // if (attention!=0) ;
}
void DlmAttention::simRun()
{
	if (attention!=0) {
		// 10, 12a
		// aTransE = gaussConvolved(g,sa,rho);
		gaussConvolved(aTransE,g,sa);

		// 11, 12b
		// nslDiff(a1,1.,lambda_a * (- a + aTransE  - beta_a*sa.sum() + kappa_ah*sh));
		NSLmult_equal(temp,kappa_ah,sh);
		NSLmult_minus(temp,beta_a,sa.sum());
		NSLadd_plus(temp,aTransE);
		NSLadd_minus(temp,a);
		NSLmult_equal(temp,lambda_a,temp);
		nslDiff(a,1.,temp);

		// 13 
		computeOutputFunc(sa,a,rho);
	} // if (attention!=0)
}

// h

DlmH::DlmH(nsl_string str,NslModule* parent) : 
	DlmAHCommon(str,parent),
	sa("sa",this),      // cell activities ;
	hTransE("hTransE",this),    
	d("d",this),
	beta_h("beta_h",this),     
	beta_ac("beta_ac",this),        // strength of global inhibition compensating ... ;
	kappa_hs("kappa_hs",this),     
	kappa_hh("kappa_hh",this),     
	kappa_ha("kappa_ha",this),       // effect of attention blob on running blob ;
	lambda_p("lambda_p",this),     
	lambda_m("lambda_m",this)
{
}
// ---- Gaussian excitation kernel -------------------------------------- ;

DlmH1::DlmH1(nsl_string str,NslModule* parent) : 
	DlmH(str,parent),
	hInput("hInput",this),
	h("h",this),
	s("s",this),
	sh("sh",this)  
{
}

void DlmH1::memAlloc(int _i1max,int _j1max)
{
	i1max = _i1max;
	j1max = _j1max;
	
	hInput.memAlloc(i1max,j1max);
	hTransE.memAlloc(i1max,j1max);
	h.memAlloc(i1max,j1max);
	s.memAlloc(i1max,j1max);
	d.memAlloc(i1max,j1max);
	sh.memAlloc(i1max,j1max);
	sa.memAlloc(i1max,j1max);
	temp.memAlloc(i1max,j1max);
}

void DlmH1::initRun()
{
	hInput=0;
	h=0;
	s=0;
	sh=0;
	initGauss(g,sigma_g);
}
void DlmH1::simRun()
{
	int i1,j1;
  
	// ---- introduce 'noise' on 'h' if necessary ----------------------- ;

	if (sh.sum()<=0) {
		h[i1max/2][j1max/2]=0.10;
		h[i1max/4][j1max/2]=0.12;
		h[i1max/2][j1max/4]=0.14;
	}

	// 3 
	hTransE = gaussConvolved(g,sh,rho);

	if (attention!=0) {
		// 5a 
		// nslDiff(h,1.,- h + hTransE - beta_h*sh.sum() - kappa_hs*s
		//	+ kappa_hh*hInput + kappa_ha*(sa-beta_ac));

		temp = sa;
		NSLadd_minus(temp,beta_ac);
		NSLmult_equal(temp,kappa_ha,temp);
		NSLmult_plus(temp,kappa_hh,hInput);
		NSLmult_minus(temp,kappa_hs,s);
		NSLmult_minus(temp,beta_h,sh.sum());
		NSLadd_plus(temp,hTransE);
		NSLadd_minus(temp,h);
		nslDiff(h,1.,temp);

	} // if (attention!=0) ;
	else {
		// 5b 
		//nslDiff(h,1., - h + hTransE - beta_h*sh.sum() - kappa_hs*s
		//	+ kappa_hh* hInput);

		NSLmult_equal(temp,kappa_hh,hInput);
		NSLmult_minus(temp,kappa_hs,s);
		NSLmult_minus(temp,beta_h,sh.sum());
		NSLadd_plus(temp,hTransE);
		NSLadd_minus(temp,h);
		nslDiff(h,1.,temp);

	} // if (attention!=0) else ;

	// ---- compute 'sh' and 'sh' ---------------------------------- ;
	// ---- needs valid 'h' and 'h' ;

	// 7a 
	computeOutputFunc(sh,h,rho);

	// ---- simulate one time step on 's' and 's' ------------------- ;
	// ---- needs valid 'h' and 'h' ;

	// 8a 
	d = h - s;
	for (i1=0; i1<i1max; i1++)
		for (j1=0; j1<j1max; j1++)
			// 8b 
			if (d[i1][j1]>0)
				d[i1][j1] = d[i1][j1]*lambda_p.elem();
			else
				d[i1][j1] = d[i1][j1]*lambda_m.elem();
	// 8c 
	nslDiff(s,1.,d);
}

// h

DlmH2::DlmH2(nsl_string str,NslModule* parent) : 
	DlmH(str,parent),
	showLayer("showLayer",this),        
	avTimeLimit("avTimeLimit",this),
	workOnAverage("workOnAverage",this),  
	gallerySize("gallerySize",this),  
	shShow("shShow",this),
	shMax("shMax",this)	     // maximum activity in the layers 2
{
	nsl_string _str;
	for (int i=0; i<gallerySizeMax; i++) {
		strcpy(_str,"shSum");
		sprintf(_str,"%s(%d)",_str,i);
		shSum[i].initialize(_str,this);
		strcpy(_str,"hInput");
		sprintf(_str,"%s(%d)",_str,i);
		hInput[i].initialize(_str,this);
		strcpy(_str,"h");
		sprintf(_str,"%s(%d)",_str,i);
		h[i].initialize(_str,this);
		strcpy(_str,"s");
		sprintf(_str,"%s(%d)",_str,i);
		s[i].initialize(_str,this);
		strcpy(_str,"sh");
		sprintf(_str,"%s(%d)",_str,i);
		sh[i].initialize(_str,this);
		strcpy(_str,"skipModel");
		sprintf(_str,"%s(%d)",_str,i);
		skipModel[i].initialize(_str,this);
	}
}

void DlmH2::memAlloc(int _i2max,int _j2max)
{
	i2max = _i2max;
	j2max = _j2max;
	
	hTransE.memAlloc(i2max,j2max);
	d.memAlloc(i2max,j2max);
	sa.memAlloc(i2max,j2max);
	shMax.memAlloc(i2max,j2max);
	temp.memAlloc(i2max,j2max);
	shShow.memAlloc(i2max,j2max);

	for (int i=0; i<gallerySizeMax; i++){
		shSum[i].memAlloc();
		hInput[i].memAlloc(i2max,j2max);
		h[i].memAlloc(i2max,j2max);
		s[i].memAlloc(i2max,j2max);
		sh[i].memAlloc(i2max,j2max);
	}
}

void DlmH2::initRun()
{
	int model;
	
	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}

	for (model=0; model<=gallerySize; model++) {
		hInput[model] = 0;
		h[model] = 0;
		s[model] = 0;
		sh[model] = 0;
	}
	shMax = 0;
	initGauss(g,sigma_g);
}
void DlmH2::simRun()
{
	int model,i2,j2;

	// 4 
	hTransE = gaussConvolved(g,shMax,rho);

	if (attention!=0) {
		// 5a
		for (model=modelLoIndex; model<=modelHiIndex; model++) {
			if (skipModel[model].get_data() == 1) continue;
			// 6a 
			//nslDiff(h[model],1.0, - h[model] + hTransE
			//	- beta_h*shSum[model] - kappa_hs* s[model]
			//	+ kappa_hh* hInput[model] + kappa_ha*(sa-beta_ac));

			temp = sa;
			NSLadd_minus(temp,beta_ac);
			NSLmult_equal(temp,kappa_ha,temp);
			NSLmult_plus(temp,kappa_hh,hInput[model]);
			NSLmult_minus(temp,kappa_hs,s[model]);
			NSLmult_minus(temp,beta_h,shSum[model]);
			NSLadd_plus(temp,hTransE);
			NSLadd_minus(temp,h[model]);
			nslDiff(h[model],1.,temp);

		}
	} // if (attention!=0) ;
	else {
		// 5b
		for (model=modelLoIndex; model<=modelHiIndex; model++) {
			if (skipModel[model].get_data() == 1) continue;
			// 6b 
			//nslDiff(h[model],1.,- h[model] + hTransE
			//	- beta_h*shSum[model] - kappa_hs* s[model]
			//	 + kappa_hh* hInput[model]);

			NSLmult_equal(temp,kappa_hh,hInput[model]);
			NSLmult_minus(temp,kappa_hs,s[model]);
			NSLmult_minus(temp,beta_h,shSum[model]);
			NSLadd_plus(temp,hTransE);
			NSLadd_minus(temp,h[model]);
			nslDiff(h[model],1.,temp);
		}
	}
  
	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model].get_data() == 1) continue;
		// 7b 
		computeOutputFunc(sh[model],h[model],rho);
	}

	// ---- compute 'shMax' ----------------------------------------- ;
	// ---- needs valid 'sh' ;

	for (i2=0; i2<i2max; i2++)
		for (j2=0; j2<j2max; j2++) {
			shMax[i2][j2] = 0;
			for (model=modelLoIndex; model<=modelHiIndex; model++) {
				if (skipModel[model].get_data() == 1) continue;
				// 7c 
				shMax[i2][j2] = NSLmax(shMax[i2][j2],sh[model][i2][j2]);
			}
		}

	//  shSum = 0;
	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model].get_data() == 1) 
			shSum[model] = 0; // continue;
		else // 7d 
			shSum[model] = sh[model].sum();
	}

	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model].get_data() == 1) continue;
		// 9a 
		d = h[model] - s[model];
		for (i2=0; i2<i2max; i2++)
			for (j2=0; j2<j2max; j2++)
				// 9b 
				if (d[i2][j2]>0)
					d[i2][j2] = d[i2][j2]*lambda_p.elem();
				else
					d[i2][j2] = d[i2][j2]*lambda_m.elem();
		// 9c 
		nslDiff(s[model],1.,d);
	}

	if (nslSystem.getSimTime()>avTimeLimit)
		workOnAverage = 0;

	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}
}

// Recognition

DlmRecognition::DlmRecognition(nsl_string str,NslModule* parent) : 
	DlmCommon(str,parent),
	avTimeLimit("avTimeLimit",this),
	workOnAverage("workOnAverage",this),  
	gallerySize("gallerySize",this),  
	lambda_r("lambda_r",this),                // time constant for the recognition dynamics ;
	r_theta("r_theta",this)                 // threshold for model suppression ;
{
	nsl_string _str;
	for (int i=0; i<gallerySizeMax; i++) {
		strcpy(_str,"shSum");
		sprintf(_str,"%s(%d)",_str,i);
		shSum[i].initialize(_str,this);
		strcpy(_str,"rec");
		sprintf(_str,"%s(%d)",_str,i);
		rec[i].initialize(_str,this);
		strcpy(_str,"skipModel");
		sprintf(_str,"%s(%d)",_str,i);
		skipModel[i].initialize(_str,this);
	}
}

void DlmRecognition::initRun()
{
	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}

	rec[0] = 0.;
	for (int model=1; model<=gallerySize; model++)
		rec[model] = 1.;
}
void DlmRecognition::simRun()
{
	int model;
	float recshSumMax;            // temporary variable ;
	  
	if (workOnAverage == 0) {
		recshSumMax = 0;
		for (model=modelLoIndex; model<=modelHiIndex; model++)
			// 22 
			recshSumMax = NSLmax(recshSumMax,rec[model].elem()*shSum[model].elem());
		// 23 
		for (model=1; model<=gallerySize; model++)
			nslDiff(rec[model],1.,lambda_r*rec[model]^(shSum[model] - recshSumMax));
	} // if (!workOnAverage) ;

  // ---- determine which models to skip ------------------------------- ;
  // ---- needs valid 'rec' ;
  
	for (model=modelLoIndex; model<=modelHiIndex; model++) {
		if (skipModel[model].get_data() == 1) continue;
		if (model>0 && rec[model] <= r_theta) {
			skipModel[model] = 1;
		}
	}
  
    // ---- print the current recognition status: ----------------------- ;
    // ---- time, showed layer, recognition values ;

    for (model=1; model<=gallerySize; model++)
		if (rec[model]>r_theta)
			NSLoutput(" ");
		else
			NSLoutput(" *");
	NSLoutput("_");

	if (nslSystem.getSimTime()>avTimeLimit)
		workOnAverage = 0;

	if (workOnAverage == 1) {
		modelLoIndex=0; 
		modelHiIndex=0;
	} else {
		modelLoIndex=1; 
		modelHiIndex=int(gallerySize.elem());
	}

}

DlmAHCommon::DlmAHCommon(nsl_string str,NslModule* parent) : 
	DlmCommon(str,parent),
	rho("rho",this),            // slope radius of squashing function ;
	attention("attention",this),      // attention!=0 indicates use of  attention dynamics ;
	sigma_g("sigma_g",this),  // Gauss width of excitatory kernel ;
	g("g",this,9),      // Gaussian interaction kernel ; gSize
	temp()
{
}

float DlmAHCommon::initGauss(NslFloat0& sigma_g,float x)
{
	return exp(-x*x/(2*sigma_g.elem()*sigma_g.elem()));
}

void DlmAHCommon::initGauss(NslFloat1& g,NslFloat0& sigma_g)
{
	for (int k=0; k<g.get_imax(); k++) {
		float x = float(k-g.get_imax()/2);
		g[k] = exp(-x*x/(2*sigma_g.elem()*sigma_g.elem()));
	}
}
// ---- convolution with Gauss kernel  ---------------------------------- ;

NslFloat2 DlmAHCommon::gaussConvolved(NslFloat1& g,NslFloat2& s,NslFloat0& rho)
{
	return NSLconv_zero(g.expand_row(1),NSLconv_zero(g.expand_col(1),s));
}

void DlmAHCommon::gaussConvolved(NslFloat2& out,NslFloat1& g,NslFloat2& s)
{
	NSLconv_zero_col_equal(out,g,s);
	NSLconv_zero_row_equal(out,g,out);
}

// ---- output function 'sigma' ----------------------------------- ;

float DlmAHCommon::outputFunc(float h,NslFloat0& rho) {
  float rhoInv = 1./rho.elem();
  if (h <=0.)              return 0.;
  else if (h >=rho) return 1.;
  else                     return sqrt(h*rhoInv); 
}
void DlmAHCommon::computeOutputFunc(NslFloat2& sh, NslFloat2& h, NslFloat0& rho)
{
  int i,j;
  for (i=0; i<sh.get_rows(); i++)
    for (j=0; j<sh.get_cols(); j++)
      sh[i][j] = outputFunc(h[i][j],rho);
}

AslSchemaModel<DLM> _DLM("DLM");