/* CR.c  Carandini/Ringach model for orientation tuning
    i.e. a 1-dimensional neural field with DOG kernel */


# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# include <time.h>

# include "random.h"
# include "nn.h"


# define N1  131

# define STEPSIZE 1.

# define SIGMA0  30.
# define SIGMA11  7.
# define SIGMA12 30.


float tau1  = 10.0;
float teta1 = 0.00;

Vector      inp, phi, f1,
            link11, link12;

Matrix      J11, J12;




static BaseType kernel(BaseType x);
static BaseType kernel(x)
BaseType x;
{
  return (BaseType)( exp( - .5 * x * x ) / sqrt(2.*M_PI) ) ;
}

static BaseType Input(BaseType x);
static BaseType Input(x)
BaseType x;
{
  return (BaseType)( exp( - .5 * x * x ) ) ;
}



SwitchValue sspikes = OFF;

SliderValue sI    = 85;
SliderValue sJ11  = 10;
SliderValue sJ12  = 60;

SliderValue steta1 = 0;
SliderValue stau1  = 100;
SliderValue sp     = 100;


BEGIN_DISPLAY

 SWITCH( "spikes", sspikes )

 SLIDER( "Input", sI, -500, 10000 )
 SLIDER( "J11", sJ11, 0, 500)
 SLIDER( "J12", sJ12, 0, 500)

 SLIDER( "theta", steta1, 0, 100)
 SLIDER( "tau", stau1, 0, 500)
 SLIDER( "power", sp, 0, 500)


 WINDOW("Signale")

  RASTER( "phi", AR, AC, phi, VECTOR, N1, 0, 0.0, 1.0, MAX(128/N1, 1) )
  RASTER( "f(phi)", NR, AC, f1, VECTOR, N1, 0, 0.0, 1.0, MAX(128/N1, 1) )

 WINDOW( "Matrizen")

  IMAGE( "J11", AR, AC, J11, CONSTANT MATRIX, N1, N1, 0.0, 1.0, 1)
  IMAGE( "J12", AR, NC, J12, CONSTANT MATRIX, N1, N1, 0.0, 1.0, 1)

 WINDOW("Links")

  RASTER( "link11", AR, AC, link11, VECTOR, N1, 0, 0.0, 1.0, MAX(128/N1, 1) )
  RASTER( "link12", NR, AC, link12, VECTOR, N1, 0, 0.0, 1.0, MAX(128/N1, 1) )

 WINDOW("Functions")

  FUNCTION( "phi", AR, AC, phi, VECTOR, N1, 0, 0, 0, 0, 0., 1.0 )
  FUNCTION( "f(phi)", NR, AC, f1, VECTOR, N1, 0, 0, 0, 0, 0., 1.0 )

 WINDOW("Plot phi-phi") 
 
  PLOT( "phi", AR, AC, phi, VECTOR, N1, 0, 0, 0, 0., 1.0, 
                       phi, VECTOR, N1, 0, 0, 0, 0., 1.0, 2 )
 
 WINDOW("Graphs")
 
  GRAPH( "phi",    R0, C0, phi, VECTOR, N1, 0, 0, 0, 0., 1.0 ) 
  GRAPH( "f(phi)", NR, C0, f1,   VECTOR, N1, 0, 0, 0, 0., 1.0 ) 

END_DISPLAY


/* NO_OUTPUT */

BEGIN_OUTPUT

 OUTFILE("phi")

  SET_SAVE_FILE_FLAG( THISFILE, ASCII, ON)
  SAVE_VARIABLE( "phi", phi, VECTOR, N1, 0,  SKIP | GRID , 
                   TSkip(1), Grid(0, N1, 4, 0,0,0) )
                   
 OUTFILE("f1")

  SET_SAVE_FILE_FLAG( THISFILE, ASCII, ON)
  SAVE_VARIABLE( "f1", f1, VECTOR, N1, 0,  SKIP | GRID , 
                   TSkip(1), Grid(0, N1, 4, 0,0,0) )

END_OUTPUT





int gaussinput(n, v, w)
int n;
Vector v;
float w;
{
  int i;
  float k1,k2;
  k1 = n/2.;
  k2 = 1.0/w;
  for (i=0;i<n;i++)
    v[i] = Input( k2*(i-k1));
}


int main_init()
{

  randomize(  time(NULL) % 32000  );
  SET_STEPSIZE(STEPSIZE);

  inp   = Get_Vector(N1);
  phi   = Get_Vector(N1);
  f1    = Get_Vector(N1);

  link11 = Get_Vector(N1);
  link12 = Get_Vector(N1);

  J11 = Get_Matrix(N1, N1);
  J12 = Get_Matrix(N1, N1);

  Make_Func_Band_Matrix( N1, J11, kernel, 1., SIGMA11);
  Make_Func_Band_Matrix( N1, J12, kernel, 1., SIGMA12); 
}



init()
{
  gaussinput(N1, inp, SIGMA0 );
  Clear_Vector( N1, phi );
  Clear_Vector( N1, f1);
  Clear_Vector(N1, link11);
  Clear_Vector(N1, link12);
}



step()
{
  long i,j,k;
  float h1;

  tau1 = .1*stau1;
  teta1 = .01*steta1;

  for (j=0; j<N1; j++)
  {
    leaky_integrate( tau1, phi[j], 
                     0.001*( sI   * inp[j]   
                          + sJ11 * link11[j] 
                          - sJ12 * link12[j] ) );
   h1 = f1[j] = TLIN( phi[j] - teta1 );
   
   if (sp != 100)
     if( h1 > 0)                    
       f1[j] = pow( h1, .01 * sp );
      
    if (sspikes)
      f1[j] = bool_noise( f1[j] ); 

  }

  Mult(N1, N1, J11, f1, link11);
  Mult(N1, N1, J12, f1, link12);
}