// DLMdisplay.C

#include "nsl_include.h"
#include "DLMwin.h"
#include "DLMdisplay.h"

inline int i1Index(int i2, int i1R, int frame,int i1max,int i1Rmax,int i2max) 
{return i1R + frame + int((i1max-2*frame-i1Rmax)*i2/float(i2max-1)+0.5);};
inline int j1Index(int j2, int j1R, int frame,int j1max,int j1Rmax,int j2max) 
{return j1R + frame + int((j1max-2*frame-j1Rmax)*j2/float(j2max-1)+0.5);};

DlmDisplay::DlmDisplay(nsl_string str,NslModule* parent) : NslModule(str,parent),
	winT1("traj1",winWidth,winWidth,580,25), // display window ;
	winT2("traj2",winWidth,winWidth,580,320),
	winW("netW",winWidth,winWidth,900,25),
	winC("netC",winWidth,winWidth,900,320),

	netWeightExp("netWeightExp",this), // exponent for the net display computation ;
	gravI2("gravI2",this), // center of gravity for the connections ;
	gravJ2("gravJ2",this),  // center of gravity for the connections ;

	correl21LI("correl21LI",this), // exponent for the net display computation ;
	w21("w21",this),
	Sh1("Sh1",this),     
	Sh2("Sh2",this)
{
}

void DlmDisplay::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;

	gravI2.memAlloc(i2max,j2max); // center of gravity for the connections ;
	gravJ2.memAlloc(i2max,j2max); // center of gravity for the connections ;
			
    correl21LI.memAlloc(i2max,j2max,i1Rmax,j1Rmax);
	w21.memAlloc(i2max,j2max,i1Rmax,j1Rmax);
	Sh1.memAlloc(i1max,j1max);
	Sh2.memAlloc(i2max,j2max);
}

void DlmDisplay::initRun()
{
  // ---- clear trajectory windows ---------------------------------- ;
	winT1.Clear();
	winT2.Clear();

    showNet(correl21LI,frame,winC,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax,
    	netWeightExp,gravI2,gravJ2);

	showNet(w21,frame,winW,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax,netWeightExp,gravI2,gravJ2);
}

void DlmDisplay::simRun()
{
	showNet(correl21LI,frame,winC,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax,
    	netWeightExp,gravI2,gravJ2);

	showNet(w21,frame,winW,i1max,j1max,i2max,j2max,i1Rmax,j1Rmax,
    	netWeightExp,gravI2,gravJ2);

    showTrajectory(Sh1,frame,winT1);
  
//	if (gallerySize.elem()>0)
		showTrajectory(Sh2,0,winT2);
}

// ==== compute center of gravity of layer activity =================== ;

void DlmDisplay::centerOfGravity(NslFloat2& Sh, float centerI, float centerJ) {
  int i,j;
  int imax = Sh.get_imax();
  int jmax = Sh.get_jmax();

  // ---- compute the center of gravity for Sh ;
  float sumI, sumJ, sumW, weight;
  sumI = 0; sumJ = 0; sumW = 0;
  for (i=0; i<imax; i++)
    for (j=0; j<jmax; j++) {
      weight = Sh[i][j];
      sumI += weight*i;
      sumJ += weight*j;
      sumW += weight;
    };

  if (sumW==0.) {
    centerI = -1;
    centerJ = -1;
    return;
  }

  centerI = sumI/sumW;
  centerJ = sumJ/sumW;
}

// ==== trajectory display function =================================== ;

/* In order to follow the course of the blobs, the following function
   displays the center of gravity of the blobs as a sequence of little
   dots. */

void DlmDisplay::showTrajectory(nsl_matrix &Sh, int frame, XMonoWin &win) {

  int imax = Sh.get_imax();
  int jmax = Sh.get_jmax();

#ifdef WIN_FRAME
  frame = 0;
#endif \\ WIN_FRAME ;

  // ---- compute the center of gravity for Sh ;
  float centerI, centerJ;
  centerOfGravity(Sh,centerI,centerJ);  

  // ---- print '_' and return if layer has no activity ;
  if (centerI==-1.) {
    NSLoutput("_"); return;}

  // ---- print '^' and return if layer has no consistent blob, ... ;
  // ... i.e. if the center of gravity does not ly on an active cell ;
  if (Sh.elem(int(centerI+0.5),int(centerJ+0.5))==0.) {
    NSLoutput("^"); return;}

  // ---- plot a dot at the center of gravity position ;
  win.Dot(int((centerJ-frame+0.5)/(jmax-2*frame)*winWidth),
	  int((centerI-frame+0.5)/(imax-2*frame)*winWidth),0);
  win.Flush();
}


// ==== net display function ========================================== ;

/* As an alternative to the matrix display of the connections, the
   following function provides a net display. The net has one node per
   neuron. The location of the nodes represent the receptive fields
   the neurons. Neighbours are connected by edges. The exponent
   netWeightExp determines the relative weight of strong
   links. netWeightExp=1 results in a linear mean over the weights,
   netWeightExp=20 favours very much the strongest links. */

void DlmDisplay::showNet(NslFloat4& M21, int frame, XMonoWin &win, 
	int i1max, int j1max, int i2max, int j2max, int i1Rmax, int j1Rmax, 
	NslFloat0& netWeightExp, NslFloat2& gravI2, NslFloat2& gravJ2) 
{
  int i1R,j1R,i2,j2,i1,j1;
  
#ifdef WIN_FRAME
  frame = 0;
#endif \\ WIN_FRAME ;

  // ---- compute the center of receptive fields for the matrix M21 ;
  double sumI, sumJ, sumW, weight;
  for (i2=0; i2<i2max; i2++)
	for (j2=0; j2<j2max; j2++) {
		sumI = 0; sumJ = 0; sumW = 0;
		for (i1R=0; i1R<i1Rmax; i1R++) 
			for (j1R=0; j1R<j1Rmax; j1R++) {
				i1=i1Index(i2,i1R/j1Rmax,frame,i1max,i1Rmax,i2max);
				j1=j1Index(j2,j1R%j1Rmax,frame,j1max,j1Rmax,j2max);
				if (M21[i2][j2][i1R][j1R]>0) {
					weight = pow(M21[i2][j2][i1R][j1R],netWeightExp.elem());
					sumI += weight*(i1-frame);
					sumJ += weight*(j1-frame);
					sumW += weight;
				}
				if (sumW>0) {
					gravI2[i2][j2] = (sumI/sumW+0.5)/(i1max-2*frame);
					gravJ2[i2][j2] = (sumJ/sumW+0.5)/(j1max-2*frame);
				}
				else {
					gravI2[i2][j2] = -i1max;
					gravJ2[i2][j2] = -j1max;
				}	
			}
  }
  
  // ---- draw the net ;
  win.Clear();
  win.Text(int(gravJ2.elem(0,0)*winWidth),
	   int(gravI2.elem(0,0)*winWidth),"o",0);
  win.Text(int(gravJ2.elem(0,j2max-1)*winWidth),
	   int(gravI2.elem(0,j2max-1)*winWidth),"x",0);
  for (i2=0; i2<i2max; i2++) 
    for (j2=0; j2<j2max; j2++) {
      if (gravJ2.elem(i2,j2)<0) continue;
      if ((i2<i2max-1)&&(gravJ2.elem(i2+1,j2)>0))
	  win.Line(int(gravJ2.elem(i2,j2)*winWidth),
		   int(gravI2.elem(i2,j2)*winWidth),
		   int(gravJ2.elem(i2+1,j2)*winWidth),
		   int(gravI2.elem(i2+1,j2)*winWidth),0);
      if ((j2<j2max-1)&&(gravJ2.elem(i2,j2+1)>0))
	  win.Line(int(gravJ2.elem(i2,j2)*winWidth),
		   int(gravI2.elem(i2,j2)*winWidth),
		   int(gravJ2.elem(i2,j2+1)*winWidth),
		   int(gravI2.elem(i2,j2+1)*winWidth),0);
    };

  // ---- draw a frame around the net ;
  win.Line(         0,         0,         0,winWidth-1,0);
  win.Line(         0,winWidth-1,winWidth-1,winWidth-1,0);
  win.Line(winWidth-1,winWidth-1,winWidth-1,         0,0);
  win.Line(winWidth-1,         0,         0,         0,0);
  win.Flush();

}