// ***************************************************************************
//			    NSL RETINA MODEL .C
//
//	Retina model description:
//	------------------------ 
//	Quantitative modeling of Responses of Anuran retina: Stimulus Shape and 
//	Size Dependency. In this simulations we test the model's ability to
//	reproduce quantitatively tabulated data on the dependency on stimulus
//	shape and size (Ewert 1976).
//	Input to the model is Light on the receptors (40X40 to simulate the 
//	receptieve field of a single ganglion cell of each type).
//	The model also simulates "simple" Horizontal and Bipolar cells.
// 	The output of the model the temporal firing rate of a  Ganglion cell of 
//    each type (R2, R3, R4). So the graphs given by NSL show the firing rate 
//	of a ganglion cell vs. the time.
//    
//	Notes:
//	-----
// 	No Trailing edge effect since Ewert computed his data with only
//	the response to the leading edge.
// 	No Leaky Integrators for the Ganglion Cells.
//
// 	Goal: Match Ewert's quantitative data on the Toad's retinal ganglion 
// 	cells.
//  	Status: goal achieved, paper submitted to Vision Research.
//
// 	File : "retina.c"
//      Contents: Definition of the neural net and main code. (see "retina.nsl")
// 	Last update: April 1st, 1992 
// 	By: Fernando Corbacho 
//
// ***************************************************************************

#include "nsl_include.h"
#include "Retina.h"

// ___________________________ AVERAGE EXACT METHOD ____________________________

nsl_matrix diff_ae(nsl_matrix& V,nsl_matrix& S,nsl_matrix& Prev, nsl_data& tm)
{

         float dt0,tc,tmp;
         int Vmax = V.get_imax();    // Size of the matrices.
         
	 nsl_matrix term1 (Vmax,Vmax);
	 nsl_matrix term2 (Vmax,Vmax);

	if (NSL_SYSTEM != NULL)
		dt0 = NSL_SYSTEM->getSimDelta();
	else
	{
		cmd_error("*** AVERAGE EXACT METHOD ***");
		cmd_error("Unknown Model");
		return term1;
	}

	if (dt0 == 0 || tm.get_data() == 0)
	{
		cmd_error("*** AVERAGE EXACT METHOD ***");
		cmd_error("delta: ",dt0);
		cmd_error("time constant: ",tm.get_data());
		cmd_error("Bad Integration Parameters");
		return term1;
	}

/* ********************* This is the real Calculation ********************** */

	tc = tm.get_data(); // time constant.
	tmp = exp(-dt0/tc);

	term1 = (nsl_matrix) ( (1 - tmp) * S + tmp * V ) ;
	term2 = (nsl_matrix) ( (S - Prev) * (tmp * (tc + dt0) - tc) / dt0); 
	 
	return (nsl_matrix) (term1 + term2);
}

nsl_matrix newconv(nsl_matrix& a,nsl_matrix& b)
// a is the Mask and b is the input layer.
{
         int saimax = a.get_imax();
         int sajmax = a.get_jmax();
         int sbimax = b.get_imax();
         int sbjmax = b.get_jmax();

        int leftbound = 1;       // 32; for the 72x72
        nsl_matrix c(1,1); // Make this variable size.
			   // c(8,8) for 72x72

        for (int i = 0; i < leftbound; i = i+4)
	{

           for (int j = 0; j < leftbound; j = j+4)
           {    
                num_type val = (num_type) 0;

                for (int m = 0; m < saimax; m++)
                   for (int n = 0; n < sajmax; n++)
                        val = val + a.elem(m,n) * b.elem(i+m,j+n);
                c.elem(i/4,j/4) = val;
			
           }
	}
        return c;
}

nsl_matrix nslGauss(int imax,int jmax,nsl_data& dia,nsl_data& sig,nsl_data& wgt)
{
        int i,j;
        num_type rr,rad;

        nsl_matrix q(imax,jmax);

//      int imax = q.get_imax();
//      int jmax = q.get_jmax();
        int imax2 = imax/2;
        int jmax2 = jmax/2;
             
        num_type rad0 = dia.elem()/2.0;
        num_type sum = 0.0;
        if (dia.elem() == 0.0)
                return q;

        for (i=0; i<imax; i++)
            for (j=0; j<jmax; j++)
            {
                /* rr is in units of 1/2 degrees, same as kernel array */
                rr = (num_type) sqrt((num_type) pow((num_type)(i-imax2),2.0) + 
                        (num_type) pow((num_type)(j-jmax2),2.0));
                /* rad is in units of degrees */
                rad = rr/2.0;
                if (rad < rad0)
                     q[i][j] = (num_type) 
                        exp(-(rad*rad)/(2.0*sig.elem()*sig.elem()));
                sum = sum + q[i][j];
            }

        if (sum == 0.0)
                q = 0;
        else 
                q = q/sum*wgt;

        return q;
}

// To Provide the Blurring edge effect use the following function:
// extern void stimulus(nsl_matrix&,nsl_data&,nsl_data&,nsl_data&,
//                    nsl_data&,nsl_data&,nsl_data&,nsl_data&);

Visin::Visin(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	out("out",this,size1,size2),
	in("in",this,size1,size2)
{
}

void Visin::initRun()
{
	in = 0;
}

// stimulus(r, r0, c0, r1,c1, vx ,vy, flag); // if needed blurred edge effect.
 
void Visin::simRun()	
{
	in = 0;		// needed to reset all values
	in.run();	//Get the stimuli. Defined in the file "retina.nsl"
	out = in;
}

Receptor::Receptor(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	in("in",this,size1,size2), 	
	r("r",this,size1,size2),	// Receptors
	rB("rB",this,size1,size2)	// For sustained response.
{
}

void Receptor::initRun()	// **** Receptors ****
{
	r = 0;
}

void Receptor::simRun()	// **** Receptors ****
{
	float t = NSL_SYSTEM->getSimTime();
	
	rB = r;
// stimulus(r, r0, c0, r1,c1, vx ,vy, flag); // if needed blurred edge effect.
 
	if (t <= 7.0)	      //1.55 in order to stop it in the ERF.
   	{
		r = in;	//Get the stimuli. Defined in the file "retina.nsl"
   	}
  	else
		r = rB;
}

Horizontal::Horizontal(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	hor("hor",this,size1,size2),	// Horizontal cells
	hor_tm("hor_tm",this),		// Horizontal cells time constant.
	H_level("H_level",this)	// Uniform horizontal cell potential
                        //    (0 if the background is bright, 1 if dark.)
{
}

void Horizontal::initRun()	// **** HORIZONTAL CELLS ****
{
	hor = 0;
}

void Horizontal::simRun()	// **** HORIZONTAL CELLS ****
{
	nslDiff(hor,hor_tm,H_level - hor);
}

BipolarOff::BipolarOff(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	r("r",this,size1,size2),	// Receptors
	hor("hor",this,size1,size2),	// Horizontal cells
	hbc("hbc",this,size1,size2),	// Hyperpolarizing Bipolar cells
	pbh("pbh",this,size1,size2)
{
}

void BipolarOff::simRun()	// **** HYPERPOLARIZING BIPOLARS ****
{
	hbc = r - hor;
	pbh = NSLmax(hbc,0);
}

BipolarOn::BipolarOn(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	r("r",this,size1,size2),	// Receptors
	hor("hor",this,size1,size2),	// Horizontal cells
	dbc("dbc",this,size1,size2),	// Depolarizing Bipolar cells
	pbd("pbd",this,size1,size2)
{
}

void BipolarOn::simRun()	// **** DEPOLARIZING BIPOLARS ****
{
	dbc = hor - r;
	pbd = NSLmax(dbc,0);
}

AmacrineOff::AmacrineOff(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	hbc("hbc",this,size1,size2),	// Hyperpolarizing Bipolar cells
	prev("prev",this,size1,size2), // Previous value of the receptor 
	axh("axh",this, size1,size2), 	//  ***** Off-channel Amacrines *****
	ath1("ath1",this,size1,size2), 
	ath2("ath2",this,size1,size2), 
	ath("ath",this,size1,size2),
	old_ath("old_ath",this,size1,size2),	// previous value of the cell.
	axh_tm("axh_tm",this)		// Amacrine cells time constants.
{
}

void AmacrineOff::initRun()	// **** Off-Channel Amacrine Cells ****
{
	prev = 0;
	axh = 0; ath1 = 0; ath2 = 0; ath = 0; 
	dt = NSL_SYSTEM->getSimDelta();
}
void AmacrineOff::simRun()	// **** Off-Channel Amacrine Cells ****
{
	// DIFF.eq(axh,axh_tm) = hbc - axh; 	// Euler method.
        axh = diff_ae(axh,hbc,prev,axh_tm);	// Average Exact Method.	
	prev = hbc;

	old_ath = ath;
	ath1 = 5*(hbc - axh);
	ath2 = exp(-dt/axh_tm.elem())*old_ath;

	ath = NSLmax(ath1,ath2) ;
}

AmacrineOn::AmacrineOn(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	dbc("dbc",this,size1,size2),	// Depolarizing Bipolar cells
	prev("prev",this,size1,size2), // Previous value of the receptor 
	axd("axd",this,size1,size2), 	//  ***** On-channel Amacrines *****
	atd1("atd1",this,size1,size2), 
	atd2("atd2",this,size1,size2),
	atd("atd",this,size1,size2), 
	old_atd("old_atd",this,size1,size2), 	// previous value of the cell.
	axd_tm("axd_tm",this)
{
}

void AmacrineOn::initRun()	// **** On-Channel Amacrine Cells ****
{
	prev = 0;
 	axd = 0; atd1 = 0; atd2 = 0; atd = 0;
	dt = NSL_SYSTEM->getSimDelta();
}

void AmacrineOn::simRun()	// **** On-Channel Amacrine Cells ****
{
	// DIFF.eq(axd,axd_tm) = dbc - axd;
	axd = diff_ae(axd, dbc, prev, axd_tm); 
	prev = dbc;

	old_atd = atd;
	atd1 = 5*(dbc - axd);             // This is because no BIPOLARS.
 	atd2 = exp(-dt/axd_tm.elem())*old_atd;

	atd = NSLmax(atd1,atd2);  
}

R2::R2(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	pbh("pbh",this,size1,size2),
	ath("ath",this,size1,size2),
	atd("atd",this,size1,size2), 
	r2("r2",this,1,1), 
	r2_i("r2_i",this,1,1),	// previous value of r2.
	sust("sust",this,1,1),  // Sustained input to R2 (from the bipolar cells).
	sust_irf("sust_irf",this,1,1),
	sust_erf("sust_erf",this,1,1),
	r2t("r2t",this,1,1),
	r2f("r2f",this,1,1),	// Scaled R2.
	pbh_irf("pbh_irf",this),
	pbh_erf("pbh_erf",this),
	erf("erf",this,size1,size2 ), 	// R2 Excitatory Receptive Field.
	irf("irf",this,size1,size2 ), 	// R2 Inhibitory Receptive Field.
	rf("rf",this,size1,size2 ),		// DOG.
	erf_dia("erf_dia",this),		// diameter of the Gaussian mask.
	irf_dia("irf_dia",this),
	erf_sig("erf_sig",this),		// Standard deviation of the Gaussian.
	irf_sig("irf_sig",this),
	erf_wgt("erf_wgt",this),		// Weight.
	irf_wgt("irf_wgt",this),
	thresh("thresh",this),		// Thresholds for ganglion cells.
	k("k",this),			// scale constant for the ganglion cells
	tm("tm",this),		// Time constants for the ganglion cells.
	trailing("trailing",this),
	temp("temp",this,size1,size2),
	stat("stat",this,3,3)  		 // **** stores the statistics ****
				 // Will contain the Average firing rates
				 // for all ganglion cells.
{
}

void R2::initRun()
{
	r2 = 0; r2t = 0; r2f = 0; 
	
	int size1 = erf.get_imax();
	int size2 = erf.get_jmax();
	
// Construct the Gaussian kernels for the ganglion cells.
	erf = nslGauss(size1,size2,erf_dia,erf_sig,erf_wgt);
	irf = nslGauss(size1,size2,irf_dia,irf_sig,irf_wgt);

	rf = erf - irf;            // DOG for the r2 ganglion cells
}
void R2::simRun()
{
	temp = ath + trailing * atd;   // trailing is the effect of the trailing
				       // edge set to 0 to get Ewert's data.

	// sustained input.
        sust_erf = newconv(erf, pbh_erf * pbh);
	sust_irf = newconv(irf, pbh_irf * pbh);

	sust = sust_erf - sust_irf;	

	r2 = newconv(rf, temp) + sust;    // New convolution and No Leaky Integ.
	
	r2_i = r2;

	r2t = NSLramp(r2);
	r2f = k*r2t;		// Scaling.

	if  (r2f.elem(0,0) >  thresh.elem())  //  ****** STATISTICS ******
	{
		// Compute the average firing rate.
		stat.elem(0,0) = stat.elem(0,0) + r2f.elem(0,0);
		stat.elem(0,1)++;
		stat.elem(0,2) = stat.elem(0,0) / stat.elem(0,1);
	}
}

R3::R3(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	ath("ath",this,size1,size2),
	atd("atd",this,size1,size2), 
	old_ath("old_ath",this,size1,size2),
	old_atd("old_atd",this,size1,size2), 
	r3("r3",this,1,1), 
	r3t("r3t",this,1,1),
	r3f("r3f",this,1,1), 
	erf("erf",this,size1,size2 ),
	irf("irf",this,size1,size2 ),
	erf_dia("erf_dia",this),
	irf_dia("irf_dia",this),
	erf_sig("erf_sig",this),
	irf_sig("irf_sig",this),
	erf_wgt("erf_wgt",this),
	irf_wgt("irf_wgt",this),
	thresh("thresh",this),
	k("k",this),			// outputs.
	tm("tm",this),
	p("p",this),		// Weight of trailing edge for R3.
	all("all",this,size1,size2),	// effect of both amacrine channels.
	old("old",this,size1,size2),	// former value of amacrines.
	stat("stat",this,3,3)  		 // **** stores the statistics ****
				 // Will contain the Average firing rates
				 // for all ganglion cells.
{
}

void R3::initRun()
{
	r3 = 0; r3t = 0; r3f = 0;

	int size1 = erf.get_imax();
	int size2 = erf.get_jmax();
	
	erf = nslGauss(size1,size2,erf_dia,erf_sig,erf_wgt);
	irf = nslGauss(size1,size2,irf_dia,irf_sig,irf_wgt);
}
void R3::simRun()
{
	all = p * atd + ath;
	old = p * old_atd + old_ath;
	
        // DIFF.eq(r3,tm) = -r3 + newconv(erf, all) - newconv(irf, old); 
	r3 = newconv(erf, all) - newconv(irf, old);

	r3t = NSLramp(r3);
	r3f = k*r3t;

	if (r3f.elem(0,0) > thresh.elem())   //  ***** STATISTICS *****
	{
		stat.elem(1,0) = stat.elem(1,0) + r3f.elem(0,0);
		//   stores the summation of the responses.
		stat.elem(1,1)++;
		//   stores the number of time steps .
		stat.elem(1,2) = stat.elem(1,0) / stat.elem(1,1);
		//   stores the Average Firing Rate up to the moment.
	}
}

R4::R4(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	ath("ath",this,size1,size2),
	atd("atd",this,size1,size2), 
	r4("r4",this,1,1), 
	r4t("r4t",this,1,1), 
	r4f("r4f",this,1,1),	// Scaled R4.
	erf_dia("erf_dia",this),
	erf_sig("erf_sig",this),
	erf_wgt("erf_wgt",this),
	erf("erf",this,size1,size2), 
	thresh("thresh",this),
	k("k",this),
	tm("tm",this),
	stat("stat",this,3,3)  		 // **** stores the statistics ****
				 // Will contain the Average firing rates
				 // for all ganglion cells.
{
}

void R4::initRun()
{
	r4 = 0; r4t = 0; r4f = 0;

	int size1 = erf.get_imax();
	int size2 = erf.get_jmax();
	
// Construct the Gaussian kernels for the ganglion cells.
	erf = nslGauss(size1,size2,erf_dia,erf_sig,erf_wgt);
}

void R4::simRun()
{
	// DIFF.eq(r4,r4_tm) = -r4 + newconv(erf,(ath - atd)); 
      r4 = newconv(erf, (ath - atd)); 

	r4t = NSLramp(r4);
	r4f = k*r4t/(r4t+0.2);  // Squashing function.
  
	if (r4f.elem(0,0) > thresh.elem())   //  ***** STATISTICS *****
	{
		stat.elem(2,0) = stat.elem(2,0) + r4f.elem(0,0);
		stat.elem(2,1)++;
		stat.elem(2,2) = stat.elem(2,0) / stat.elem(2,1);
	}
}

Retina::Retina(char* _name, NslModule* parent,int size1,int size2) :
	NslModule(_name, parent),
	visin("visin",this,size1,size2),
	receptor("receptor",this,size1,size2),
	horizontal("horizontal",this,size1,size2),
	bipolarOn("bipolarOn",this,size1,size2),
	bipolarOff("bipolarOff",this,size1,size2),
	amacrineOn("amacrineOn",this,size1,size2),
	amacrineOff("amacrineOff",this,size1,size2),
	r2("r2",this,size1,size2),	
	r3("r3",this,size1,size2),	
	r4("r4",this,size1,size2)
{	
	initModule();
}

void Retina::initModule()
{
	nslConnect(visin.out,receptor.in);
	nslConnect(receptor.r,bipolarOn.r);
	nslConnect(horizontal.hor,bipolarOn.hor);
	nslConnect(receptor.r,bipolarOff.r);
	nslConnect(horizontal.hor,bipolarOff.hor);
	nslConnect(bipolarOn.dbc,amacrineOn.dbc);
	nslConnect(bipolarOff.hbc,amacrineOff.hbc);
	nslConnect(bipolarOff.pbh,r2.pbh);
	nslConnect(amacrineOff.ath,r2.ath);
	nslConnect(amacrineOn.atd,r2.atd);
	nslConnect(amacrineOff.ath,r3.ath);
	nslConnect(amacrineOn.atd,r3.atd);
	nslConnect(amacrineOff.old_ath,r3.old_ath);
	nslConnect(amacrineOn.old_atd,r3.old_atd);
	nslConnect(amacrineOff.ath,r4.ath);
	nslConnect(amacrineOn.atd,r4.atd);
}

RetinaModel::RetinaModel() : NslModel("RetinaModel"),
	retina("retina",this,40,40)
{
}

AslSchemaModel<RetinaModel> _RetinaModel("RetinaModel");