////////////////////////////////////////////////////////////////////////////
//
// Limulus.cc
//
// by Jim Lesko       May 1997
//
// Compiles using G++ in UNIX; I don't know how well it will work for other
// systems.  Try the following settings for good results:
// Max Inh .2   Length 1  ep don't care     -    good for no delta.
// MI .3 length 2.
//
/////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <iostream.h>
#include <stream.h>
#include <math.h>
#include <fstream.h>
#include <stdlib.h>
#include <time.h>

// Default offsets for stimulus
#define XOFFSET 6
#define YOFFSET 6

// Default Eye size
#define MAXELEMS 20

// Some protos -- bad style.  Don't let me catch YOU doing this!
dumpVector(float *v);
saveVector(float *v,char *n,char dot);


// Datatypes:

struct vector {                            // vector: self-explanitory
float element[MAXELEMS];
};

struct matrix {                            // matrix: ditto
float element[MAXELEMS][MAXELEMS];
};

struct stateMatrix {                       // stateMatrix: for ease of passing.
float initialStateMatrix[MAXELEMS][MAXELEMS];
float state[MAXELEMS][MAXELEMS];
float weights[MAXELEMS][MAXELEMS][MAXELEMS][MAXELEMS];
};

struct parameters {                        // parameters: same as for stateMatrix.
float upperLimit;
float lowerLimit;
int lengthConstant;
float epsilon;
float maxStrength;
int iterations;                            // iterations is last of the original
int delta;                                 // settings.
float contrast;
int noise;
int eyetype;
int noiserate;
int xoffset;
int yoffset;
};


// checkXY: takes in a pair of ints and returns 0 if they're > 0 && < MAXELEMS

int checkXY(int i,int j) {
int error=0;
if(i<0 || i>(MAXELEMS-1) || j<0 || j>(MAXELEMS-1)) {
  cout << "\nInput out of bounds.  Pos: "<<i<<","<<j<<".\n"; 
  error = 1;
  }
return error;
}


// zeroMatrix:  Sets state and initialStateMatrix to all zeros.

void zeroMatrix(stateMatrix *s) {
int i,j;

for (j=0;j<MAXELEMS;j++) {

for (i=0;i<MAXELEMS;i++) {
	s->initialStateMatrix[i][j]=0;  }
	s->state[i][j] = s->initialStateMatrix[i][j];  }
}


// initializeStateVector:  Reads in file input and creates initialStateMatrix

void initializeStateVector(stateMatrix *s, parameters *p) {
int i=YOFFSET,j=XOFFSET,error=0;
FILE *source;
char filename[20];
char input;

zeroMatrix(s);    // Start with a matrix of zero's (ala Matlab)

cout<<"\nWhat is the filename of this input:";
cin >> filename;

source = fopen(filename,"r");

if(source == NULL) {
	cout << "File " << filename << "  couldn't be opened.\n";
	exit(-1);
	}

cout << "\nFile opened.\n";

// The following reads in the data line by line; if the line goes past the
// limits of the "retina,"  an error will be returned by checkXY.  Data is
// read in until the EOF is reached.
 
i+=p->yoffset;                            
j+=p->xoffset;
while (error ==0) {
  input=fgetc(source);
  if(input=='\n') { 
      i++;
      j=XOFFSET+p->xoffset;
      error = checkXY(i,j);
  } 
  else {
    if(input == EOF) { error = 1; } else {
      if((input < '0' || input > '1')) { cout << "\n Non-numeric input."; error=1; } 
      else {
	error=checkXY(i,j);
  	s->initialStateMatrix[i][j] = (int)(input-'0');
	j++;
  	}
     }
    }
  }
fclose(source);
}


// reMapStateMatrix:  provides contrast and noise

void reMapStateMatrix (stateMatrix *s,parameters *p) {
  int i,j;

  time_t now;             // use the current time to seed random.
  time(&now);
  srand(now);

  for(i=0;i<MAXELEMS;i++) {
    for(j=0;j<MAXELEMS;j++) {
	s->initialStateMatrix[i][j] = (s->initialStateMatrix[i][j] * (100- (1-p->contrast)* p->upperLimit)) + (100-(p->contrast *100));
	if (p->noise!=0) {
		cout <<"."<<(rand()/(RAND_MAX/p->noise));
		if((rand()/(RAND_MAX/10)) <= p->noiserate) {
		s->initialStateMatrix[i][j] += (rand()/(RAND_MAX/p->noise)*(p->upperLimit/10));
		if(s->initialStateMatrix[i][j] > p->upperLimit)
			s->initialStateMatrix[i][j] = p->upperLimit;
		}
	}
      s->state[i][j] = s->initialStateMatrix[i][j];
    }
  }
}


// initializeParameters:  Takes in user choices.  Depending on what I was doing,
// some lines are commented out and parameters set to constants.  For full user
// input remove all p->x=y; lines and uncomment cin/cout lines.

void initializeParameters(parameters *p)  {
int i,j;
char inChar;

cout << "Lateral Inhibition Demonstration...\n";
cout << "\nEnter Inhibitory stregth:";
cin >> p->maxStrength;
//p->maxStrength = 0;
cout << "\nEnter Length Constant:";
cin >> p->lengthConstant;
//p->lengthConstant = 1;
//cout << "\nEnter epsilon:";
//cin >> p->epsilon;
p->epsilon=0;
p->upperLimit = 100;
p->lowerLimit = 0;
p->iterations=50;
//cout << "\nUse delta rule?";
//cin >> inChar;
//if ((inChar=='y')|(inChar=='Y')) { p->delta=1; }
//else {p->delta=0;}
p->delta = 0;
cout << "\nContrast (0-none - 10-full):";
cin>> i;
p->contrast = (float)i/10;
cout <<"\nNoise rating: (0-10):";
cin >> p->noise;
cout << "\nNoise Rate (0 = never 10 = 100%):";
cin >> p->noiserate;
cout<<"\nEye type:  (1. limulus, 2. Center On, 3. Center Off)";
cout<<"           (4-6: above w/ edge reduction 7: Center On, radius 1):";
cin >> p->eyetype;
//p->eyetype=1;
/*cout << "\nOffset stimulus?";
cin >> inChar;
if (inChar == 'Y' || inChar == 'y') {
	cout << "\nX offset from center:";
	cin >> p->xoffset;
	cout << "\nY offset from center:";
	cin >> p->yoffset;
	}
else { */
	p->xoffset = p->yoffset = 0;
	//}
}


// initializeDisplay:  Simply prints a header to make reading the data easier.
// The name of the procedure is a product of what it did under Anderson's 
// program:  it actually cleared space on the display.

void initializeDisplay()  {
cout << "\n\nROW      | COLUMNS -> \n";
cout << "         0..............5.............10.............15..........19\n";
}


// displayStateMatrix:  Takes in a state matrix and displays it in integral
// values from 0 to 9.  Note that the start and dispChar parameters are artifacts
// of the 2D version of this program.  start was used to determine the displayed
// cell range (which 20 out of MAXELEMS), and dispChar was used as the character
// displayed at the appropriate place (ie, an "x" or a "+" to mark the value).
// the v[][] is the most important parameter, as it is the matrix that will
// be displayed.  Note that rounding occurs here to keep the values integral
// and between 0 and 9; any other values would be difficult to print to STDOUT.

displayStateMatrix (float v[][MAXELEMS], int start, char dispChar) {
	int i,j,k,value;
 	initializeDisplay();	
	
	for(i=start;i<(start+20);i++) {
  	cout << "\n\nrow "<<i<<":  ";	
	if(i<10) { cout << " "; }
	for(j=start;j<(start+20);j++) {
	value = (int)floor ((v[i][j])/10.1);
	if (value != 0) {
	  cout << value << "  "; }
	else {
	  cout << "   "; }
	}}

}


// innerProduct:  Takes two matricies and returns their inner product.
// used in determining inhibitory effects.  

float innerProduct (float m1[][MAXELEMS], float m2[][MAXELEMS]) 
{
        int i,j;
	float sumOfProducts = 0.0;
	for(j=0;j<MAXELEMS;j++) {
       	  for(i=0;i<MAXELEMS;i++) { 
	    sumOfProducts = m1[i][j] * m2[i][j] + sumOfProducts;
	    if(fabs(sumOfProducts) <.00001) {sumOfProducts = 0.0;}
          }
	}
	//cout << "\nSoP: "<<sumOfProducts<<"\t";
  	return sumOfProducts;		
}


// subtractVectors:  Takes two vectors and subracts them.  Passes result in *diff.

subtractVectors (float *v1, float *v2,float *diff) 
{
	int i;
	
       	for (i=0;i<MAXELEMS;i++) { 
             diff [i] = v1[i] - v2[i];
        }
}


// subractMatricies:  Takes two matracies and subtracts them via iterative calls
// to subtractVectors.  Returns the resulting matrix in diff[][].

subtractMatricies (float m1[][MAXELEMS],float m2[][MAXELEMS],float diff[][MAXELEMS])
{
	int i,j;
	for(j=0;j<MAXELEMS;j++) {
	  for(i=0;i<MAXELEMS;i++) {
	    diff [i][j] = m1[i][j] - m2[i][j];
	  }
	}
}


// vectorLength:  Takes a vector and returns the sum of squares of that vector.
// NOTE:  this isn't true vector length; it simply adds up sum of squares for
//        the following procedure, which calculates the matrix "length."

float vectorLength (float v[MAXELEMS])
{
        int i;
	float sumOfSquares=0;
	
	for(i=0;i<MAXELEMS;i++) {
             sumOfSquares = v[i]* v[i] + sumOfSquares;
	}
	return sumOfSquares;
}


// matrixLength:  uses iterative calls to vectorLength to determine the "length"
// of the matrix.

float matrixLength (float m[][MAXELEMS])
{
	int j;
	float sumOfSquares=0;

	for(j=0;j<MAXELEMS;j++)  {
	    sumOfSquares += vectorLength(m[j]);
	}
	sumOfSquares = sqrt(sumOfSquares);
	return sumOfSquares;
}


// Limit state matrix:  takes in the state matrix and parameters, and alters
// the matrix to fit within those parameters.  This is used to keep the simulation
// within the observed biological limits of the Limulus eye.  It should be noted
// that this is rate coding, not actual potentiation.

limitStateMatrix (float m[][MAXELEMS], parameters *p)
{
	int i,j;

	for(j=0;j<MAXELEMS;j++) {
	  for(i=0;i<MAXELEMS;i++) {
	    //cout << " *"<<m[i][j];
	    if(m[i][j] > p->upperLimit) {m[i][j] = p->upperLimit;}
	    if(m[i][j] < p->lowerLimit) {m[i][j] = p->lowerLimit;}
          }
        }
}


// Convergence test:  removed due to the fact that good parameters can keep
// the simulation from converging.  Also, I don't think I got this to work in
// my first set of runs, and never got back to fixing it up once the model
// was working.
/*
convergenceTest(float newMat[][MAXELEMS], stateMatrix *s) 
{
	float difference[MAXELEMS][MAXELEMS];

        subtractVectors (newVec, s->state,difference);
        if (matrixLength (difference) > 1.0) {
            cout << ">>> Convergence Problem.\n";
            cout << "    Length difference between last iterations is";
            cout << vectorLength(difference);
	}
	cout << "c";
}
*/


// computeISM: Compute Inhibited State Matrix:
// Runs through the entire state matrix and figures out what the max. inhibition
// of that cell is.

computeISM(stateMatrix *s, parameters *p)
{
     int i,j;

     // cout<<(1/MAXELEMS);
     for(j=0;j<MAXELEMS;j++)  {
       for(i=0;i<MAXELEMS;i++) {
	 s->state[i][j] = s->initialStateMatrix[i][j] + (float)(1.0/(3.14*p->lengthConstant))*innerProduct(s->weights[i][j],s->initialStateMatrix);
	 limitStateMatrix(s->state,p);
	   }
     }
}


// computeInhibitedStateMatrix:  iterative method of figuring out the inhibitory
// effects of each cell on all others.  (Delta method)

computeInhibitedStateMatrix(stateMatrix *s, parameters *p)
{
     int i,j,count;
     float newStateMatrix[MAXELEMS][MAXELEMS];

     cout << "\nComputing inhibited state matrix...\n";
     for (count=0;count < p->iterations;count++) {  
	if(count%5==0) { cout<<"Iteration "<<count<<"...\n"; }
        for(j=0;j<MAXELEMS;j++) {
	  for(i=0;i<MAXELEMS;i++) {
          newStateMatrix[i][j] = s->state[i][j] + p->epsilon * ( s->initialStateMatrix[i][j] +
	  innerProduct (s->weights[i][j], s->state) - s->state[i][j]);
	}
	}
	 limitStateMatrix (newStateMatrix,p);
         for(i=0;i<MAXELEMS;i++) {
	   for(j=0;j<MAXELEMS;j++) {
	     s->state[i][j] = newStateMatrix[i][j];
	     }
	
	/* if (i == p->iterations) { convergenceTest(newStateMatrix,s); }*/
	 //cout << "\nCompare: "<<s->state[4][4]<<" and "<<newStateMatrix[4][4];
     } 
  } 
}

float square (float num)
{
  return (num*num);
}


// makeWeightMatrx:  takes in the state matrix and parameters,and creates an 
// inhibition matrix for each cell.

void makeWeightMatrix(stateMatrix *s,parameters *p)
{
	int i,j,iw,jw;
	float distFact;
	
	cout << "\nMaking weight matrix...this may take some time.\n";
	cout << "\nType "<<p->eyetype<<".\n";	
	for (j=0;j<MAXELEMS;j++){
	  for (i=0;i<MAXELEMS;i++) {
	    for (iw=0;iw<MAXELEMS;iw++) {
	      for (jw=0;jw<MAXELEMS;jw++) {
		     s->weights[i][j][iw][jw] = -1*p->maxStrength * 
		       exp((float)((-1*sqrt((square(i-iw))+(square(j-jw))))/ p->lengthConstant));
		     if ((p->eyetype == 2) || (p->eyetype == 5)) {
			s->weights[i][j][i][j] = 0;  }

		     if ((p->eyetype == 3) || (p->eyetype == 6)) {
			s->weights[i][j][iw][jw] = s->weights[i][j][i][j]-s->weights[i][j][iw][jw]; }

		     if (p->eyetype == 7) {
		       if (iw>0) 	
			  s->weights[i][j][i-1][j] = 0;
		       if (iw<19)
			  s->weights[i][j][i+1][j] = 0;
		       if (jw>0)
			  s->weights[i][j][i][j-1] = 0;
		       if (jw<19)
			  s->weights[i][j][i][j+1] = 0;
		       s->weights[i][j][i][j] = 0;
		     }
		     
		     if ((p->eyetype > 3) && (p->eyetype < 7)) {
			if (i<5 || i >15 || j<5 || j>15) 
			  s->weights[i][j][iw][jw] *= 2;
			}
		     }
		 }
	   }	
	 }	
}


// dumpWeights:  takes in the state matrix and prints out all the weights
// for each cell.  This is simply a debugging procedure.

void dumpWeights (stateMatrix *s) {
	int i,j;
	for (i=0;i<MAXELEMS;i++) {
	  cout << "\nrow "<<i<<": ";
	  for(j=0;j<MAXELEMS;j++) {
	    cout << " " << s->weights[10][10][i][j];
          }
        }
}


// saveVector: takes in the stateMatrix and prints a (brute-force) ascii listing
// of the values of the matrix.  (As a column vector)

saveVector (stateMatrix *s)
{
 char fileName[12],input;
 int i,j;
 cout << "\n\nSave the result?"; 
 cin >> input;
 if ((input=='y')||(input=='Y')) {
	cout << "\nVery well.  Name the file:";
	cin >> fileName;
	ofstream out(fileName);
	if(!out) {
		cout << "\nUnable to create file.\n";
		}
	else {
		for (i=0;i<MAXELEMS;i++) {
		  for(j=0;j<MAXELEMS;j++) {
		    out << (s->state[i][j]/100) << " ";
		    }
		}
	     }
	}
}


// SaveExcelFile:  Used previously for the 2D version:  saved selected
// weight vectors.

/*
saveExcelFile (stateMatrix *s) 
{
	  int i;
	  char name[30];
	  cout << "\nSaving Excel File... name:";
	  cin >> name;
	  ofstream out(name);
	  if(!out) {
		cout << "\nUnable to create Excel echo.\n";
		}
	  else {
		out <<"Initial State\tWeight0\tWeight15\tWeight42\tFinal State\n";
		for (i=0;i<MAXELEMS;i++) {
			out << s->initialStateVector[i] << "\t";
			out << s->weights[0][i] << "\t";
			out << s->weights[15][i] << "\t";
			out << s->weights[42][i] << "\t";
			out << s->state[i] << "\n";
			}
	}
	cout << "\nDone.\n";
}


// saveMatlabFile:  Saved vectors for Matlab using saveVector.

saveMatlabFile(stateMatrix *s) 
{
	char input;
	cout << "\nSave Matlab Files?";
	cin >> input;
	if ((input=='y')||(input=='Y')) {
		saveVector (s->initialStateVector,"Initial State Vector",'+');
		saveVector (s->stateVector,"Final State Vector",'*');
		saveVector (s->weights[0],"Weight 0 Vector",'o');
		saveVector (s->weights[11],"Weight 11 Vector",'o');
		saveVector (s->weights[62],"Weight 62 Vector",'o');
		}
}

*/


//////////////////////////////////////////////////////////////////////////
//
//  MAIN
//

void main () {
char contin = 'y';
parameters *param;
stateMatrix *s;

while(contin=='y') {
param = new parameters;
s = new stateMatrix;

initializeParameters(param);
initializeStateVector(s,param);
reMapStateMatrix(s,param);
displayStateMatrix (s->state, 0, '+');
makeWeightMatrix (s,param);
if(param->delta == 1) {
  computeInhibitedStateMatrix (s,param);
}
else {
  computeISM (s,param);
  }
displayStateMatrix(s->state,0,'+');
saveVector(s);
/*
saveMatlabFile (s);
saveExcelFile (s);*/

delete s;
delete param;

cout <<"\nRun another simulation?";
cin >> contin;
  }
cout << "\n";
}


// THE END.