/* wr.c */


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

# include "random.h"
# include "gen_sim.h"
# include "nn.h"
# include "vector.h"



# define tau1  3.
# define tau2  5.0

long stp = 0;	

float sim_time, noise_fac;

Layer       input,
            pot1, pot2,
            f1, f2;

SpikeLayer  out1, out2;

# define L_SIZE11   10.0       /* FWHM in columns (float) */
# define M_SIZE11   11         /* Kernel dimension (int)  */
# define FM_SIZE11  (2*M_SIZE11+1)

# define L_SIZE12   10.0       /* FWHM in columns (float) */
# define M_SIZE12    8         /* Kernel dimension (int)  */
# define FM_SIZE12  (2*M_SIZE12+1)

# define L_SIZE21   10.0       /* FWHM in columns (float) */
# define M_SIZE21    8        /* Kernel dimension (int)  */
# define FM_SIZE21  (2*M_SIZE21+1)


UniKernel   kernel11,
            kernel12,
            kernel21;  

Layer       link11, 
            link12,
            link21; 



# define barlength  25
# define barskip     0    /* 5  */
# define barsigma    7
# define BARINITOFFS 14.

double  yy1, yy2;
int bardirection = 1;




SwitchValue santi = OFF;
SwitchValue scent = OFF;

SliderValue sI1     = 85;
SliderValue sI2     = 85;
SliderValue sI     = 85;
SliderValue snoise = 20;
SliderValue sJ11  = 100;
SliderValue sJ12  = 40;
SliderValue sJ21  = 600;
SliderValue sspeed = 0;

BEGIN_DISPLAY

 SWITCH( "anti", santi)
 SWITCH( "center", scent )

 SLIDER( "Signal Input", sI, 0, 1000 )

 SLIDER( "E ", sI1, -200, 200 )
 SLIDER( "I ", sI2, -200, 200 )

 SLIDER( "noise", snoise, 0, 1000 )

 SLIDER( "J11", sJ11, 0, 500)
 SLIDER( "J12", sJ12, 0, 300)
 SLIDER( "J21", sJ21, 0, 1000)
 SLIDER( "speed", sspeed, 0, 1000);



      
 WINDOW("Signale")

  IMAGE( " input ", AR, AC, input, LAYER, xsize, ysize, 0.0, 2.1, 1)
  IMAGE( "  pot1  ", AR, NC, pot1, LAYER, xsize, ysize, -.5, 1.0, 1)
  IMAGE( "  out1  ", NR, C0, out1, SPIKE_LAYER, xsize, ysize, 0.0, 1.0, 1)


 WINDOW("Graphs")

  GRAPH( "pot1", AR, AC, pot1, MATRIX, xsize, ysize, 32, 32, -.5, 1.5)
  GRAPH( "pot2", NR, AC, pot2, MATRIX, xsize, ysize, 32, 32, -.5, 1.5)
  
WINDOW("Inhibition")

  IMAGE( " input ", AR, AC, input, LAYER, xsize, ysize, 0.0, 2.1, 1)
  IMAGE( "  pot2  ", AR, NC, pot2, LAYER, xsize, ysize, -.5, 1.0, 1)
  IMAGE( "  out2  ", NR, C0, out2, SPIKE_LAYER, xsize, ysize, 0.0, 1.0, 1)   


WINDOW("Links")

  IMAGE( "link11", AR, AC, link11, LAYER, xsize, ysize, 0.0, 1.0, 1)
  IMAGE( "link12", NR, AC, link12, LAYER, xsize, ysize, 0.0, 1.0, 1)
  IMAGE( "link21", NR, AC, link21, LAYER, xsize, ysize, 0.0, 1.0, 1)

WINDOW("kernels")

  IMAGE( " k11", AR, AC, kernel11, CONSTANT LAYER, 
                 FM_SIZE11, FM_SIZE11, 0.0, 1., 5)
  IMAGE( " k12", NR, AC, kernel12, CONSTANT LAYER, 
                 FM_SIZE12, FM_SIZE12, 0.0, 1., 5)
  IMAGE( " k21", NR, AC, kernel21, CONSTANT LAYER, 
                 FM_SIZE21, FM_SIZE21, 0.0, 1., 5)

END_DISPLAY


/* NO_OUTPUT */

BEGIN_OUTPUT

 OUTFILE("phi1")

  SET_SAVE_FILE_FLAG( THISFILE, ASCII, ON)
  SAVE_VARIABLE( "phi1 (pot2)", pot1, MATRIX, xsize, ysize,  SKIP | GRID , 
                   TSkip(2), Grid(26, xsize, 100, 32, ysize, 100) ) 
                   
 OUTFILE("phi2")

  SET_SAVE_FILE_FLAG( THISFILE, ASCII, ON)
  SAVE_VARIABLE( "phi2 (pot2)", pot1, MATRIX, xsize, ysize,  SKIP | GRID , 
                   TSkip(2), Grid(38, xsize, 100, 32, ysize, 100) ) 


/* 
 OUTFILE("phi")

  SET_SAVE_FILE_FLAG( THISFILE, ASCII, ON)
  SAVE_VARIABLE( "pot1", pot1, MATRIX, xsize, ysize,  SKIP | GRID , 
                   TSkip(5), Grid(0, xsize, 16, 0, ysize, 16) )
*/                   
                   
END_OUTPUT







static void init_bars(centerflag)
int centerflag;
{
  if (centerflag)  /* center */
  {
    yy1 = yy2 = ysize/2;
  } 
  else
  {
    yy1 = BARINITOFFS;
    if (santi) 
      yy2 = ysize-BARINITOFFS;
    else 
      yy2 = BARINITOFFS;    
    bardirection = 1;
  }
}

static void move_bars()
{
  if (scent) 
  {
    init_bars(1);
    return;
  }

  if (yy1 > ysize-BARINITOFFS ||
      yy1 < BARINITOFFS)
    bardirection *= -1;
    
  yy1 += .001*bardirection*sspeed;
  if (santi) 
    yy2 -= .001*bardirection*sspeed;
  else 
    yy2 += .001*bardirection*sspeed;
    
/*   printf("yy1: %f yy2: %f\n", yy1,yy2 ); */
}



static void smooth_bars( out )
Matrix out;
{
  int i, j, s1, s2, s3, s4;
  static double fac=0;
  double h;
  
  if (fac==0) fac = -.5/(float)(barsigma*barsigma);
  
  s2 = (xsize - barskip)/2;
  s1 = s2-barlength;
  s3 = (xsize + barskip)/2;
  s4 = s3 + barlength;
 
/*   printf("%d %d %d %d \n", s1 , s2, s3, s4 ); */
 
  for (j = 0; j<ysize; j++)
  {
    h = elem( out, j, s1, xsize) = gaussian( fac * (yy1-j)*(yy1-j));
    for (i=s1+1; i<s2; i++)
      elem( out, j, i, xsize) = h;

    h = elem( out, j, s3, xsize) = gaussian( fac * (yy2-j)*(yy2-j));
    for (i=s3+1; i<s4; i++)
      elem( out, j, i, xsize) = h;      
      
  }

}



int main_init()
{
  int i;

  randomize(  time(NULL) % 32000  ); 

  input = Get_Layer();
  pot1  = Get_Layer();
  f1    = Get_Layer();
  out1  = Get_SpikeLayer();
  pot2  = Get_Layer();
  f2    = Get_Layer();
  out2  = Get_SpikeLayer();
    
  kernel11 = Get_UniKernel( FM_SIZE11, FM_SIZE11 );
  kernel12 = Get_UniKernel( FM_SIZE12, FM_SIZE12 );
  kernel21 = Get_UniKernel( FM_SIZE21, FM_SIZE21 );

  Set_Circ_Func_Uni_Kernel( kernel11, FM_SIZE11, FM_SIZE11, gaussian, 
                            1., L_SIZE11, 0. );
  Set_Circ_Func_Uni_Kernel( kernel12, FM_SIZE12, FM_SIZE12, gaussian, 
                            1., L_SIZE12, 0. );
  Set_Circ_Func_Uni_Kernel( kernel21, FM_SIZE21, FM_SIZE21, gaussian, 
                            1., L_SIZE21, 0. );

  link11 = Get_Layer();
  link12 = Get_Layer();
  link21 = Get_Layer();

  SET_STEPSIZE(0.5);
  noise_fac = sqrt(24.0/step_size);
  
}

int init( )
{
  register int i,j;

  stp = 0;
  
  Clear_Layer(input);
  init_bars( scent );
  smooth_bars( input );
  
  Clear_Layer(pot1);
  Clear_SpikeLayer(out1);
  
  Clear_Layer(pot2);
  Clear_SpikeLayer(out2);  
  
}

int step()
{
  register int i,j,k;
    
  if (stp >= 36050) 
    exit (0);
  
  /********************/
  /* compute stimulus */
  /********************/
  
  move_bars();
  smooth_bars( input );

  /********************/
  /* compute ddynamics */
  /********************/

  /* excit. units */

  for (i=0; i<ysize; i++)
  {
    for (j=0; j<xsize; j++)
    {
      leaky_integrate( tau1, elem( pot1, i, j, xsize) , 
           0.001*(sI1 + sI*gauss_noise()*elem( input , i, j, xsize) 
              + sJ11*elem( link11,i,j, xsize)
              - sJ12*elem( link12,i,j, xsize)
              + (snoise*noise_fac)*(equal_noise() - 0.5) ) ) ; 
      elem( f1, i, j, xsize) = RAMP( elem( pot1, i, j, xsize) );
      elem( out1, i, j, xsize) = PROB_FIRE( elem( f1, i, j, xsize) );
      
    } /* END j */
  } /* END i */


  /* inhib. units */

  for (i=0; i<ysize; i++)
  {
    for (j=0; j<xsize; j++)
    {
      leaky_integrate( tau2, elem( pot2, i, j, xsize) , 
           0.001*( sI2 +  sJ21*elem( link21,i,j, xsize)
                 + (snoise*noise_fac)*(equal_noise() - 0.5) ) ) ; 
      elem( f2, i, j, xsize) = RAMP( elem( pot2, i, j, xsize) );
      elem( out2, i, j, xsize) = PROB_FIRE( elem( f2, i, j, xsize) );
    } /* END j */
  } /* END i */

  bConvolute_2d_Uni( out1, kernel11, xsize, ysize, FM_SIZE11, FM_SIZE11, link11);
  bConvolute_2d_Uni( out1, kernel21, xsize, ysize, FM_SIZE21, FM_SIZE21, link21);
  bConvolute_2d_Uni( out2, kernel12, xsize, ysize, FM_SIZE12, FM_SIZE12, link12);

  stp++;

} /* END of step() */