#include <retina.h>

//********************************************************************************************
//******************  The Retina functions   *************************************************


Retina::Retina(double stepp)
{
  step=stepp; input_luminosity_range=255.0; pixels_per_degree=1.0;
  sizeX=0; sizeY=0;
  //Just for safety, really...
  radScheme=0; outerPlexiformLayer=0; contrastGainControl=0;
  done=false; adjusted=false;
}


void Retina::xml_ize()
{

  // General parameters

  xmlClass::add_param(step,"temporal-step__sec");
  xmlClass::add_param( input_luminosity_range , "input-luminosity-range");
  xmlClass::add_param( pixels_per_degree , "pixels-per-degree");

  // Log-polar scheme

  xmlClass::add_child(radScheme,"");
  xmlRegister<1, Retinal_LogPolarScheme>(radScheme , "log-polar-scheme");

  // OPL schemes

  xmlClass::add_child(outerPlexiformLayer,"outer-plexiform-layer");
  xmlRegister<1, LinearOPL > ( outerPlexiformLayer, "linear-version" );

  xmlRegister<1, UndershootOPL > ( outerPlexiformLayer, "undershoot-version" );
  // Obsolete... could be left here only for compatibility of old definition files...but... who cares ;-)
  // xmlRegister<1, UndershootOPL > ( outerPlexiformLayer, "adaptation-version" );

  xmlRegister<1, GainControlOPL > (outerPlexiformLayer, "luminance-gain-control-version" );

  // Gain control scheme

  xmlClass::add_child ( contrastGainControl, "" );
  xmlRegister<1, FastGainControl >( contrastGainControl , "contrast-gain-control" );
  xmlRegister<1, FastGainControl_custom >( contrastGainControl , "contrast-gain-control-custom" );

  // Ganglion layers

  xmlClass::add_child ( ganglionLayers, "" );
  xmlRegister<2, GanglionLayer> (ganglionLayers, xml_names ("" , "ganglion-layer") );

  // Obsolete... could be left here only for compatibility of old definition files...but... who cares ;-).
  // xmlRegister<2, GanglionLayer> (ganglionLayers, "parvocellular-ganglion-layer");
  // xmlRegister<2, GanglionLayer> (ganglionLayers, "magnocellular-ganglion-layer");
}


Retina::~Retina()
{
  if(radScheme) delete radScheme;
  if(outerPlexiformLayer) delete outerPlexiformLayer;
  if(contrastGainControl) delete contrastGainControl;
}

//Return minimum size that an input image should have, given the blur parameters required for the retina in its successive layers, for the filtering to be exact at least in central pixel od the retina.
int Retina::get_minimumSize()
{
  double sumSigmaSquare=0.0;
  for(uint chann=0;chann<ganglionLayers.size();++chann)
    sumSigmaSquare=max(sumSigmaSquare,pow( ganglionLayers[chann]->sigmaPool * pixels_per_degree , 2.0 ) );
  if(outerPlexiformLayer)
    sumSigmaSquare+=(pow(outerPlexiformLayer->sigmaCenter * pixels_per_degree ,2.0)+pow(outerPlexiformLayer->sigmaSurround * pixels_per_degree ,2.0));
  if(contrastGainControl)
    sumSigmaSquare+=(pow(contrastGainControl->sigmaSurround * pixels_per_degree ,2.0));

  //taking 3.5*sigma as the 'zero' boundary for Gaussian kernels
  return (int)ceil(7*sqrt(sumSigmaSquare));
}



//This function converts to PIXELS all parameters that were set in DEGREES.
void Retina::degToPix()
{
  if(radScheme)
    radScheme->degToPix( this->pixels_per_degree );

  for(uint chann=0;chann<ganglionLayers.size();++chann)
  {
    if(ganglionLayers[chann]->spikingChannel)
    {
      Retinal_deg_obj* ret_chann = dynamic_cast<Retinal_deg_obj*>( ganglionLayers[chann]->spikingChannel);
      if(!ret_chann)
        cerr<<"ERROR in Retina::degToPix(): the created spiking channel does not inherit from a Retinal_deg_obj, as it should."<<endl;
      ret_chann->degToPix( this->pixels_per_degree );
    }
  }
}




//Helper function: adjust ONE size of ONE spiking channel as a function of input image dimension.
//output code: 0=no prob, 1=too big, 2=auto (in both cases, automatic fit to input image)
int Retina::adjustChannelSize( int* chanSize, int imSize )
{
  int ret=0;
  if( *chanSize > imSize )
    ret=1;
  else if( *chanSize == 0 )
    ret=2;
  if( ret )
    *chanSize = imSize ;
  return ret ;
}


//adjusting sizeX and sizeY, that is, the size of all maps that are going to be allocated by sub-objects of the retina, as a function of the sizes of the input images that are going to be fed to the retina...
void Retina::adjustSizes(int imX, int imY)
{

  //(0) Setting parameters expressed in retinal degrees, back to good old PIXELS !
  degToPix();



  //(1) FIX sizeX and sizeY from the spiking channels' size specifications

  int nbSC=0;

  for(uint chann=0;chann<ganglionLayers.size();++chann)
  {
    if(ganglionLayers[chann]->spikingChannel)
    {
      nbSC++;

      int status = max( adjustChannelSize( ganglionLayers[chann]->spikingChannel->sizeX , imX ),
                        adjustChannelSize( ganglionLayers[chann]->spikingChannel->sizeY , imY ) ) ;

      if ( status==2 )
        cout<<"From Retina::adjustSizes(): A spiking array <"
        << xmlFC::find_xmlObject_base( *(ganglionLayers[chann]->spikingChannel) )->xmlName
        <<"> is being automatically set to fit the size of input image."<<endl;
      else if ( status==1 )
        cout<<"Warning from Retina::adjustSizes(): A spiking array <"
        << xmlFC::find_xmlObject_base( *(ganglionLayers[chann]->spikingChannel) )->xmlName
        <<"> is being automatically REDUCED to fit the size of input image."<<endl;


      //Take the changes in cosideration at the level of the XML saving file, in degrees.
      Retinal_deg_obj* ret_chann = dynamic_cast<Retinal_deg_obj*>( ganglionLayers[chann]->spikingChannel);
      if(!ret_chann)
        cerr<<"ERROR: the created spiking channel does not inherit from a Retinal_deg_obj, as it should."<<endl;
      ret_chann->pixToDeg( this->pixels_per_degree );


      //Adapt retinal sizes sizeX and sizeY in consequence:
      sizeX=max ( sizeX,*(ganglionLayers[chann]->spikingChannel->sizeX) );
      sizeY=max ( sizeY,*(ganglionLayers[chann]->spikingChannel->sizeY) );

    }
  }



  //(1bis) If the retina has no spiking channel, then sizeX and sizeY are just defined by the input image!
  if(!nbSC)
    sizeX=imX, sizeY=imY;


  //(2) 'Delimitating' the maximum zone in the input image where retinal filtering can be exact
  int minSizeForBlurs=get_minimumSize();

  if((imX<minSizeForBlurs)||(imY<minSizeForBlurs))
    cerr<<"WARNING in Retina::adjustSizes(): Given the blur parameters required for the retina in its successive layers, the input image should be at least of size "<<minSizeForBlurs<<" for the filtering to be exact at least at central pixel. It is not the case presently."<<endl;

  sizeX=min(sizeX+minSizeForBlurs, imX); sizeY=min(sizeY+minSizeForBlurs, imY);

  adjusted=true;
}



void Retina::allocateValues(bool verbose, double initValue)
{
  if(!adjusted)
    cerr<<"ERROR in Retina::allocateValues(): use adjustSizes(imX, imY) first!"<<endl;

  double initVal = (initValue==HUGE_VAL)?
                   input_luminosity_range / 2.0 :
                   initValue ;

  if(verbose)
  {
    cout<<endl<<"The retina is currently allocating all its elements!"<<endl;
    cout<<"Size of its maps was fixed at: ("<<sizeX<<","<<sizeY<<")."<<endl;
    cout<<"Initialization of the retina (OPL only) by a uniform screen of luminance "<<initVal<<"."<<endl;

    cout<<endl<<"(1) OUTER PLEXIFORM LAYER:"<<endl;
  }



  if(outerPlexiformLayer)
  {
    outerPlexiformLayer->fatherRetina=this;
    outerPlexiformLayer->allocateValues(sizeX, sizeY, step, verbose, initVal );
  }
  else
  {
    cerr<<"Error! Your retina does not have an outer plexiform layer, which is its minimal element!"<<endl;
    exit(0);
  }



  if(verbose){cout<<endl<<"(2) CONTRAST GAIN CONTROL MECHANISM:"<<endl;}
  if(contrastGainControl)
  {
    contrastGainControl->fatherRetina=this;
    contrastGainControl->allocateValues(sizeX, sizeY, step, verbose);
  }
  else if(verbose){cout<<"No Contrast Gain Control mechanism was defined."<<endl;}



  if(verbose){cout<<endl<<"(3) SPIKING LAYERS OF GANGLION CELLS:"<<endl;}
  if(ganglionLayers.size())
  {
    int nbUnits=0, nbSpkCh=0, nbCh=ganglionLayers.size();
    for(int chann=0;chann<nbCh;++chann)
    {
      ganglionLayers[chann]->fatherRetina=this;
      ganglionLayers[chann]->allocateValues(sizeX, sizeY, step, verbose);
      if(ganglionLayers[chann]->spikingChannel)
      {
        nbSpkCh++;
        nbUnits+=ganglionLayers[chann]->spikingChannel->units.size();
      }
    }
    if(verbose){cout<<nbCh<<" layers of ganglion cells were defined, amongst which "<<nbSpkCh<<" are spiking channels, for a total of "<<nbUnits<<" spiking ganglion cells."<<endl;}
  }
  else if(verbose){cout<<"No layer of ganglion cells were defined. Retinal output will be an analogic signal corresponding to the potential of a layer of bipolar cells."<<endl;}



  done=true;
}

void Retina::feedInput(const CImg<double>& image)
{
  if(!done){adjustSizes(image.dimx(),image.dimy()); allocateValues();}

  //if the input image should be bigger than necessary, reduce it:
  if((image.dimx()==sizeX)&(image.dimy()==sizeY))
  {
    input.const_targets=&image;
    input.targets=0;
  }
  else
  {
    stockInput = image.get_crop((image.dimx()-sizeX)/2,(image.dimy()-sizeY)/2,(image.dimx()-sizeX)/2+sizeX-1,(image.dimy()-sizeY)/2+sizeY-1);
    input.targets=&stockInput;
  }

  outerPlexiformLayer->feedInput(input.read());

  if(contrastGainControl)
    contrastGainControl->feedInput(outerPlexiformLayer->outputSignal->read());

  for(uint chann=0;chann<ganglionLayers.size();++chann)
    if(contrastGainControl)
      ganglionLayers[chann]->feedInput(contrastGainControl->bipolarSignal->read());
    else
      ganglionLayers[chann]->feedInput(outerPlexiformLayer->outputSignal->read());
}

void Retina::tempStep()
{
  if(!done)
    cerr<<"error in Retina::tempStep(): feedInput first!"<<endl;
  outerPlexiformLayer->tempStep();
  if(contrastGainControl)
    contrastGainControl->tempStep();
  for(uint chann=0;chann<ganglionLayers.size();++chann)
    ganglionLayers[chann]->tempStep();
}






//*****************************************************************************************************
//******************Functions for OPL classes LinearOPL , UndershootOPL and  GainControlOPL  **********





//******************LinearOPL ******************:

LinearOPL::LinearOPL()
{
  sigmaCenter=0.0;
  tauCenter=0.01;
  nCenter=0;
  sigmaSurround=1.0;
  tauSurround=0.01;
  wSurround=0.0;
  leakyHeatEq=0;
  ampOPL=1.0; relative_ampOPL=1.0;

  excitCells = new RetinaMapFilter();
  inhibCells = new RetinaMapFilter();
  undershoot = 0 ;
  outputSignal=0, receptorSignal=0, horizontalSignal=0;
}

void LinearOPL::xml_ize()
{
  xmlClass::add_param(sigmaCenter,"center-sigma__deg");
  xmlClass::add_param(tauCenter,"center-tau__sec");
  xmlClass::add_param(nCenter,"center-n__uint");
  xmlClass::add_param(sigmaSurround,"surround-sigma__deg");
  xmlClass::add_param(tauSurround,"surround-tau__sec");
  xmlClass::add_param(wSurround,"opl-relative-weight");
  xmlClass::add_param(leakyHeatEq,"leaky-heat-equation");
  xmlClass::add_param(relative_ampOPL,"opl-amplification");

  //The undershoot is now added as a subobject:
  xmlClass::add_child( undershoot , "" ) ;
  xmlRegister<1,LinearOPL::Undershoot> ( undershoot , "undershoot" ) ;
}

void LinearOPL::Undershoot::xml_ize()
{
  xmlClass::add_param(wAdap,"relative-weight");
  xmlClass::add_param(tauAdap,"tau__sec");
}


LinearOPL::~LinearOPL()
{
  delete excitCells;
  delete inhibCells;
  if((wSurround)&&(outputSignal)) delete outputSignal;
  if(undershoot) delete undershoot ;
}



void LinearOPL::allocateValues(int sx, int sy, double st, bool verbose, double initValue)
{
  if(verbose)
  {
    cout<<"Allocating values for the LINEAR version LinearOPL."<<endl;
    if(leakyHeatEq)
      cout<<"Low-pass filter for each layer of cells is a non-separable leaky heat equation."<<endl;
    else
      cout<<"Low-pass filter for each layer of cells is a separable spatial_gaussian / temporal_negative_exponential."<<endl;
    if(undershoot)
      cout<<"A linear undershoot (transient) has been added."<<endl;
  }

  //relative_ampOPL must have values of order [1-10] for the OPL output to have typical values of order 1.
  if ( (relative_ampOPL<1) || (relative_ampOPL>20) )
    cerr<<"WARNING from the LinearOPL object: given your range of luminosities fed as input to the retina, the coeficient LinearOPL::opl-amplification is not well suited. It is a RELATIVE factor with respect to input luminosity on the retina. It should always be set between 1 and 20 for the output range of the OPL to have values of order 1. CHANGE IT in your retina definition file!"<<endl;

  ampOPL = relative_ampOPL / fatherRetina->input_luminosity_range ;

  excitCells->set_sizeX(sx).set_sizeY(sy).set_step(st);
  if(leakyHeatEq)
  {
    excitCells -> leakyHeatEquation( tauCenter , sigmaCenter * fatherRetina->pixels_per_degree );
    if ( nCenter != 0 )
      cerr<<"WARNING from the LinearOPL object: You have chosen a 'leaky heat equation filter', so parameter center-n__uint is inefficient ! It will be considered as having value zero (simple exponential filter)."<<endl;
  }
  else
  {
    excitCells -> ExpCascade(tauCenter , nCenter);
    excitCells -> set_sigmaPool(sigmaCenter * fatherRetina->pixels_per_degree );
  }
  (*excitCells) *= ampOPL;

  //We set the initial input of the filters at half the retinal input_luminosity_range
  excitCells -> set_initial_input_value ( initValue );
  excitCells -> allocateValues();
  receptorSignal = excitCells;


  if(undershoot)
  {
    BaseRecFilter adapTr( st );
    adapTr.Exp ( undershoot->tauAdap ).warningVerbose=false;
    adapTr*= undershoot->wAdap ;

    (*undershoot) .set_sizeX(sx).set_sizeY(sy).set_step(st) ;
    (*undershoot) -= adapTr ;

    undershoot -> set_initial_input_value ( excitCells->read()(sx/2,sy/2) );
    undershoot -> allocateValues();
    receptorSignal = undershoot ;
  }


  if(wSurround)
  {
    inhibCells -> set_sizeX(sx).set_sizeY(sy).set_step(st);
    if(leakyHeatEq)
      inhibCells -> leakyHeatEquation( tauSurround , sigmaSurround * fatherRetina->pixels_per_degree );
    else
    {
      inhibCells -> Exp(tauSurround);
      inhibCells -> set_sigmaPool(sigmaSurround * fatherRetina->pixels_per_degree );
    }
    (*inhibCells) *= wSurround;

    //We set the initial input of the filters at half the retinal input_luminosity_range
    //NOTA : receptorSignal = undershoot if it is present, and = excitCells otherwise.
    inhibCells -> set_initial_input_value ( receptorSignal->read()(sx/2,sy/2) );
    inhibCells -> allocateValues();
    horizontalSignal = inhibCells;

    stock_output = CImg<double>(sx,sy,1,1,0.0);
    outputSignal = new CImgReader;
    outputSignal -> targets= &stock_output;
  }
  else
    outputSignal = receptorSignal ;
}

void LinearOPL::feedInput(const CImg<double>& image)
{
  excitCells -> feedInput( image );
  if( undershoot )
    undershoot -> feedInput ( excitCells -> read() );
  //NOTA : receptorSignal = undershoot if it is present, and = excitCells otherwise.
  if( wSurround )
    inhibCells -> feedInput( receptorSignal -> read() );
}

void LinearOPL::tempStep()
{
  excitCells -> tempStep();
  if ( undershoot )
    undershoot -> tempStep();
  if(wSurround)
  {
    inhibCells -> tempStep();
    stock_output = receptorSignal -> read() - inhibCells -> read();
  }
}




//******************UndershootOPL ******************:
//************************ OBSOLETE ! **************


UndershootOPL::UndershootOPL() : wAdap(0.0), tauAdap(1.0), adapType(0)
{}

void UndershootOPL::xml_ize()
{
  LinearOPL::xml_ize();
  xmlClass::add_param(wAdap,"undershoot-relative-weight");
  xmlClass::add_param(tauAdap,"undershoot-tau__sec");

  // Obsolete... left here only for compatibility of the previous XML versions
  // xmlClass::add_param(wAdap,"adap-relative-weight");
  // xmlClass::add_param(tauAdap,"adap-tau__sec");
  // xmlClass::add_param(adapType,"adap-type");
}



void UndershootOPL::allocateValues(int sx, int sy, double st, bool verbose, double initValue)

{
  if(verbose)
  {
    cout<<"Allocating values for the linear version UndershootOPL."<<endl;
    if(leakyHeatEq)
      cout<<"Low-pass filter for each layer of cells is a non-separable leaky heat equation."<<endl;
    else
      cout<<"Low-pass filter for each layer of cells is a separable spatial_gaussian / temporal_negative_exponential."<<endl;
  }
  LinearOPL::allocateValues(sx,sy,st,verbose, initValue);

  BaseRecFilter adapTr( st );
  adapTr.Exp ( tauAdap ).warningVerbose=false;
  adapTr*= wAdap;

  afterAdap.set_sizeX(sx).set_sizeY(sy).set_step(st) ;
  afterAdap -= adapTr ;

  receptorSignal = & afterAdap ;

}

void UndershootOPL::feedInput(const CImg<double>& image)
{

  excitCells -> feedInput( image );
  afterAdap . feedInput ( excitCells -> read() );
  if(wSurround)
    inhibCells -> feedInput ( afterAdap . read() );
}

void UndershootOPL::tempStep()
{
  excitCells -> tempStep();
  afterAdap . tempStep();
  if(wSurround)
  {
    inhibCells -> tempStep();
    stock_output = afterAdap . read() - inhibCells -> read();
  }
}





//****************** GainControlOPL ******************:

GainControlOPL::GainControlOPL()
{
  ampPhoto=1.0;
  sigmaCenter=0.0;
  tauPhoto=0.01;
  sigmaSurround=1.0;
  tauSurround=0.01;
  wSurround=0.0;
  gLeak=10.0;
  ampFeedback=200.0;
  nernstFeedback=0.0;

  phototransFilter=RetinaMapFilter();
  excitCells=CondCellMap();
  controlCond=RetinaMapFilter();
  outputSignal=0, receptorSignal=0, horizontalSignal=0;
}

void GainControlOPL::xml_ize()
{
  xmlClass::add_param(ampPhoto,"input-amplification-phototransduction");
  xmlClass::add_param(sigmaCenter,"sigma-receptors__deg");
  xmlClass::add_param(tauPhoto,"tau-phototransduction__sec");
  xmlClass::add_param(sigmaSurround,"sigma-horizontal-cells__deg");
  xmlClass::add_param(tauSurround,"additional-tau-horizontal-cells__sec");
  xmlClass::add_param(wSurround,"relative-weight-surround-center");
  xmlClass::add_param(gLeak,"inert-leak-in-receptors");
  xmlClass::add_param(ampFeedback,"horizontal-feedback-amplification");
  xmlClass::add_param(nernstFeedback,"horizontal-feedback-nernst-potential");
}


GainControlOPL::~GainControlOPL()
{
  if((wSurround)&&(outputSignal)) delete outputSignal;
}



void GainControlOPL::allocateValues(int sx, int sy, double st, bool verbose, double initValue)
{
  if(verbose)
  {
    cout<<"Allocating values for a CONDUCTANCE-BASED OPL (luminance gain-control)."<<endl;
    cout<<"Very experimental element. the initValue method does not work (amongst other things!)"<<endl;
  }

  phototransFilter.set_sizeX(sx).set_sizeY(sy).set_step(st);
  phototransFilter.Exp(tauPhoto); phototransFilter.set_sigmaPool(sigmaCenter * fatherRetina->pixels_per_degree ); phototransFilter*=ampPhoto;

  excitCells.set_sizeX(sx).set_sizeY(sy).set_step(st);
  excitCells.set_g_leak(gLeak);
  int index=excitCells.new_ionicChannel(nernstFeedback);
  excitCells.inputSynapse(index).linearAmplification(ampFeedback);
  excitCells.allocateValues();
  receptorSignal=&excitCells;

  controlCond.set_sizeX(sx).set_sizeY(sy).set_step(st);
  controlCond.Exp(tauSurround); controlCond.set_sigmaPool(sigmaSurround * fatherRetina->pixels_per_degree );
  controlCond.inputSynapse().rectification(0.0);
  controlCond.allocateValues();
  horizontalSignal=&controlCond;

  if(wSurround)
  {
    stock_output=CImg<double>(sx,sy);
    outputSignal=new CImgReader;
    outputSignal->targets=&stock_output;
  }
  else
    outputSignal=&excitCells;
}

void GainControlOPL::feedInput(const CImg<double>& image)
{
  excitCells.feedCurrent(image);//the input current corresponds to port 0...
  excitCells.feedCond(controlCond.read(),1,true);//the feedback conductance corresponds to port 1...
  controlCond.feedInput(excitCells.read());
}

void GainControlOPL::tempStep()
{
  excitCells.tempStep();
  controlCond.tempStep();
  if(wSurround)
    stock_output=excitCells.read()-wSurround*controlCond.read();
}





//****************** FastGainControl:  ******************
//*******************************************************



FastGainControl::FastGainControl() :  	vThresh(0.0), Einh(0.0), exponent(2.0), symmetrize(1),
    inhibCells(0), wTransient(0), tauTransient(0.03)
{
  ampInputCurrent=100.0;
  sigmaSurround=1.0;
  tauSurround=0.01;
  gLeak=10.0;
  ampFeedback=100.0;

  excitCells=CondCellMap();
  controlCond=RetinaMapFilter();
  bipolarSignal=0, adaptationSignal=0;
}


void FastGainControl::xml_ize()
{
  xmlClass::add_param(ampInputCurrent,"opl-amplification__Hz");
  xmlClass::add_param(gLeak,"bipolar-inert-leaks__Hz");

  xmlClass::add_param(sigmaSurround,"adaptation-sigma__deg");
  xmlClass::add_param(tauSurround,"adaptation-tau__sec");
  xmlClass::add_param(ampFeedback,"adaptation-feedback-amplification__Hz");

  // Obsolete... left here only for compatibility of the previous XML versions
  // xmlClass::add_param(sigmaSurround,"amacrine-sigma__deg");
  // xmlClass::add_param(tauSurround,"amacrine-tau__sec");
  // xmlClass::add_param(ampFeedback,"amacrine-feedback-amplification__Hz");
}


void FastGainControl::allocateValues(int sx, int sy, double st, bool verbose)
{

  if(verbose)
  {cout<<"Allocating values for a conductance-based contrast-gain-control."<<endl;}

  // *********** (1) LAYER OF BIPOLAR CELLS:

  excitCells.set_sizeX(sx).set_sizeY(sy).set_step(st);
  excitCells.set_g_leak(gLeak);

  // (a) Input current in port number 0.
  excitCells.inputSynapse().linearAmplification(ampInputCurrent);

  // (b) Adapting inhibitory (or leak, if Einh=0) conductance in port number 1 -- providing CGC.
  // Einh=0 should be the law if an ON/OFF symmetry is present.
  excitCells.new_ionicChannel( Einh );

  // (c) Optional linear transient, if wTransient!=0 :
  if(wTransient)
    excitCells.new_currentInput();

  // (d) Allocating the retinal map. Pointing the reader on the map.
  excitCells.allocateValues();
  bipolarSignal=&excitCells;



  // *************** (2) ADAPTING CONDUCTANCES IN THE MEMBRANE

  controlCond.set_sizeX(sx).set_sizeY(sy).set_step(st);

  // (a) Synaptic transmission ===> Parameters for this stage can only be fixed through the 'CUSTOM' version.

  controlCond.inputSynapse().rectification( vThresh , ampFeedback );
  controlCond.inputSynapse().set_power_exp ( exponent ) ;
  if (symmetrize)
  {
    controlCond.new_currentInput();
    controlCond.inputSynapse(1).rectification( -vThresh , -ampFeedback ); //The simultaneous 'OFF'-bipolar-control
    controlCond.inputSynapse(1).set_power_exp ( exponent ) ;
  }

  // (b) Spatio-temporal averaging properties

  controlCond.Exp(tauSurround);

  if(fatherRetina) 	//If the CGC scheme is attached to a retina, use the retina's conversion from degrees to pixels.
    controlCond.set_sigmaPool( sigmaSurround * fatherRetina->pixels_per_degree );
  else 			//If the CGC scheme is tested alone, take the default conversion 1pix = 1deg.
    controlCond.set_sigmaPool( sigmaSurround );


  // (c) Allocating the retinal map. Pointing the reader on the map.
  controlCond.allocateValues();
  adaptationSignal=&controlCond;


  // ***************** (3) OPTIONAL ADDTIONAL LINEAR HIGH-PASS
  if(wTransient)
  {
    inhibCells = new RetinaMapFilter() ;
    inhibCells -> set_sizeX(sx).set_sizeY(sy).set_step(st);
    inhibCells -> Exp(tauTransient) ;
    (*inhibCells) *= wTransient ;
    inhibCells -> allocateValues() ;
  }

}

void FastGainControl::feedInput(const CImg<double>& image)
{
  // *******  (1) Inputs to bipolar cells

  // input current in port 0:
  excitCells.feedCurrent(image);
  //feedback conductance corresponds to port 1:
  excitCells.feedCond(controlCond.read(),1,true);
  //Optionally, linear feedback corresponds to port 2:
  if(wTransient)
    excitCells.feedCurrent(inhibCells->read(),2);

  // *******  (2) Inputs to adapting conductance:
  controlCond.feedInput(excitCells.read(),0);
  if(symmetrize)
    controlCond.feedInput(excitCells.read(),1); //Input from 'ON' and 'OFF' bipolar cells.

  // *******  (3) OPTIONALLY, Input to the supplementary transient (short-range amacrine cells)
  if (inhibCells)
    inhibCells->feedInput( excitCells.read() );
}

void FastGainControl::tempStep()
{
  excitCells.tempStep() ;
  controlCond.tempStep() ;
  if(inhibCells)
    inhibCells->tempStep() ;
}


//****************** FastGainControl_custom ******************:


void FastGainControl_custom::xml_ize()
{
  FastGainControl::xml_ize();
  xmlClass::add_param(exponent,"adaptation-powerlaw-exponent");
  xmlClass::add_param(Einh,"adaptation-E-inh");
  xmlClass::add_param(symmetrize,"adaptation-ONOFF-symmetry");
  xmlClass::add_param(vThresh,"bipolar-to-adaptation-threshold");

  xmlClass::add_param( wTransient , "transient-relative-weight");
  xmlClass::add_param( tauTransient , "transient-tau__sec");

  // Obsolete... left here only for compatibility of the previous XML versions
  // xmlClass::add_param(exponent,"amacrine-powerlaw-exponent");
  // xmlClass::add_param(Einh,"amacrine-E-inh");
  // xmlClass::add_param(symmetrize,"amacrine-ONOFF-symmetry");
  // xmlClass::add_param(vThresh,"bipolar-to-amacrine-threshold");
}







//*********************************************************************************************************
//******************Functions for class GanglionLayer, and derived class MagnoGanglionLayer ***************


//parvocellular cells...

GanglionLayer::GanglionLayer()
{
  sign=1.0;
  vThresh=0.0;
  valThresh=70.0;
  ampInputCurrent=100.0;
  tauTransient=0.04;
  wTransient=0.75;
  sigmaPool=0.0;

  inputCurrent=RetinaMapFilter(); supplementaryTransient=RetinaMapFilter(); spikingChannel=0;
}

void GanglionLayer::xml_ize()

{
  xmlClass::add_param(sign,"sign");
  xmlClass::add_param(tauTransient,"transient-tau__sec");
  xmlClass::add_param(wTransient,"transient-relative-weight");
  xmlClass::add_param(vThresh,"bipolar-linear-threshold");
  xmlClass::add_param(valThresh,"value-at-linear-threshold__Hz");
  xmlClass::add_param(ampInputCurrent,"bipolar-amplification__Hz");
  xmlClass::add_param(sigmaPool,"sigma-pool__deg");

  xmlClass::add_child(spikingChannel,"spiking-channel");
  xmlRegister<1, Retinal_sqSpikingChannel> (spikingChannel , "square-spiking-channel");
  xmlRegister<1, Retinal_rdSpikingChannel> (spikingChannel , "circular-spiking-channel");
}

GanglionLayer::~GanglionLayer()
{
  if(spikingChannel) delete spikingChannel;
};

void GanglionLayer::allocateValues(int sx, int sy, double st, bool verbose)
{
  sign=(sign<0)? -1.0:1.0;
  inputCurrent.set_sizeX(sx).set_sizeY(sy).set_step(st);
  supplementaryTransient.set_sizeX(sx).set_sizeY(sy).set_step(st);

  BaseRecFilter amacrineCuttingInGanglionar(st);
  amacrineCuttingInGanglionar.Exp(tauTransient).warningVerbose=false;
  amacrineCuttingInGanglionar*=wTransient;

  //This is what we had before!!: transient AFTER the synaptic rectif/pooling
  //if(wTransient) inputCurrent-=amacrineCuttingInGanglionar;

  //This is what we propose now: transient BEFORE the synaptic rectif/pooling
  if(wTransient)
    supplementaryTransient-=amacrineCuttingInGanglionar;
  supplementaryTransient.allocateValues();

  //Check the corresponding change in functions feedInput() and tempStep() !

  inputCurrent.inputSynapse().rectification(sign*vThresh,sign*ampInputCurrent,valThresh);
  inputCurrent.set_sigmaPool(sigmaPool * fatherRetina->pixels_per_degree );
  inputCurrent.allocateValues();

  if(spikingChannel)
  {
    spikingChannel->set_step(st);
    spikingChannel->allocateValues();
  }

}

void GanglionLayer::feedInput(const CImg<double>& image)
{

  //This is what we had before!!: transient AFTER the synaptic rectif/pooling
  //inputCurrent.feedInput(image);

  //This is what we propose now: transient BEFORE the synaptic rectif/pooling
  if(wTransient)
  {
    supplementaryTransient.feedInput( image );
    inputCurrent.feedInput( supplementaryTransient.read() );
  }
  else
    inputCurrent.feedInput( image );
  if(spikingChannel) spikingChannel->feedInput(inputCurrent.read());
}

void GanglionLayer::tempStep()
{
  if(wTransient)
    supplementaryTransient.tempStep();
  inputCurrent.tempStep();
}




/*******************************************************************/
/****************  UGLY OVERRIDES        ***************************/


//overriding of the 'LogPolarScheme' functions
void Retinal_LogPolarScheme::xml_ize()
{
  xmlClass::add_param( R0_deg , "fovea-radius__deg" );
  xmlClass::add_param( K_deg , "scaling-factor-outside-fovea__inv-deg" );
}

void Retinal_LogPolarScheme::degToPix( double pixPerDeg )
{
  LogPolarScheme::R0 = R0_deg * pixPerDeg ;
  LogPolarScheme::K = K_deg / pixPerDeg ;
}

void Retinal_LogPolarScheme::pixToDeg( double pixPerDeg )
{
  R0_deg = LogPolarScheme::R0 / pixPerDeg ;
  K_deg = LogPolarScheme::K * pixPerDeg ;
}


//overriding of the 'rdSpikingChannel' functions
Retinal_rdSpikingChannel::Retinal_rdSpikingChannel() : diameter_deg(0), foveaDensity_deg(1.0)
{}

void Retinal_rdSpikingChannel::xml_ize()
{
  BaseSpikingChannel::xml_ize();
  xmlClass::add_param(diameter_deg,"diameter__deg");
  xmlClass::add_param(foveaDensity_deg,"fovea-density__inv-deg");
}

void Retinal_rdSpikingChannel::degToPix( double pixPerDeg )
{
  rdSpikingChannel::diameter = (int)floor (diameter_deg * pixPerDeg ) ;
  rdSpikingChannel::foveaDensity = foveaDensity_deg / pixPerDeg ;
}

void Retinal_rdSpikingChannel::pixToDeg( double pixPerDeg )
{
  diameter_deg = rdSpikingChannel::diameter / pixPerDeg ;
  foveaDensity_deg = rdSpikingChannel::foveaDensity * pixPerDeg ;
}


//overriding of the 'sqSpikingChannel' functions
Retinal_sqSpikingChannel::Retinal_sqSpikingChannel() : sizeX_deg(0), sizeY_deg(0), uniformDensity_deg(1.0)
{}

void Retinal_sqSpikingChannel::xml_ize()
{
  BaseSpikingChannel::xml_ize();
  xmlClass::add_param(sizeX_deg,"size-x__deg");
  xmlClass::add_param(sizeY_deg,"size-y__deg");
  xmlClass::add_param(uniformDensity_deg,"uniform-density__inv-deg");
}

void Retinal_sqSpikingChannel::degToPix( double pixPerDeg )
{
  *sqSpikingChannel::sizeX = (int)floor( sizeX_deg * pixPerDeg );
  *sqSpikingChannel::sizeY = (int)floor( sizeY_deg * pixPerDeg );
  sqSpikingChannel::uniformDensity = uniformDensity_deg / pixPerDeg ;
}


void Retinal_sqSpikingChannel::pixToDeg( double pixPerDeg )
{
  sizeX_deg = *sqSpikingChannel::sizeX / pixPerDeg ;
  sizeY_deg = *sqSpikingChannel::sizeY / pixPerDeg ;
  uniformDensity_deg = sqSpikingChannel::uniformDensity * pixPerDeg ;
}
