/* ------------------------------------------------------------ */
/*                                                              */
/*     REM  Memory Version january 12   1996                    */
/*                                                              */
/*     Mark Steyvers                                            */ 
/*     msteyver@falstaff.ucs.indiana.edu                        */
/*     http://EZinfo.ucs.indiana.edu/~msteyver/home.html        */
/*                                                              */
/*     PAIRS   -  LIST LENGTH + STRENGTH +                      */
/*                LIST STRENGTH + FREQUENCY + ROC               */
/*                                                              */
/*     Note: this is NOT the program used for the simulations   */
/*           reported in the technical report; this version is  */
/*           a stripped down version of the original program    */
/* ------------------------------------------------------------ */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>         
/* the random generator (from numerical recipes) */
#include "ran1.c"

/* file handling */
#define  FILENM     "results.txt"
/* append to existing file ? */
#define  APPEND     0

/* random seed */
long     IDUM    = -3;

/* number of data points */
#define  NPOINTS    20000

/* the parameters for the model */
#define  W       20

#define  SystemG   0.4
#define  LowG      0.325
#define  HighG     0.45
float    G[ 3 ] = { LowG , SystemG , HighG  };

#define  CO      0.7    /* probability storing correct */
#define  US      0.04   /* probability storage attempt per time unit */

#define  TINCR   0.0    /* criterion for recognition in incremental storage */

#define  T1     10      /* use T1 when strong, T2 when weak */
#define  T2     7       /* in normal op. t1=10 t2=7;
						   for incremental storage, set t1=2 and t2=1 !!!!!!!!  */

#define  REPnew 7       /* storage attempts in incremental storage */
#define  REPold 3

#define  NLEN   4
#define  NSTR   4
#define  NCRI   11
#define  NFRE   4

int   LEN[ NLEN ] = { 4 , 10 , 20 , 40 };
         
/* criteria used in log scale for ZROC plots */
float CRI[ NCRI ] = { -1.0 , -0.8 , -0.7 , -0.2 , 0.0 ,
         		      0.2  , 0.7 , 1.0 , 1.25 , 1.75 , 2.5 };
                                               
/* what kind of simulations are we running? */
#define  DOLENGTH    1
#define  DOSTRENGTH  1
#define  DOFREQ      1

#define  DOROC         1   /* calculate ZROC curves ? */
#define  DOPAIRS       0   /* do we do pairs ? */
#define  DOINCREMENT   0   /* incremental storage ? */
#define  SHOWZSCORES   0
                                                      
/* these are the default settings for length, strength and frequency when
   one of these are held constant */
                                                      
#define  FIXED_LEN    2  /* look up length in length array          */
#define  FIXED_STR    0  /* 0 : strength= t1; 1 : strength= t1&t2 (mixed) ; 2 : strength= t2    */
#define  FIXED_CRI    4  /* look up criterion in criterion array    */
#define  FIXED_FRE    2  /* 0 = lowg ; 1= mixed 2=high */

/* some constants */
#define  MAXK   50
#define  NMAX   40
#define  IMAX   100
#define  LMAX   100
#define  MAXPOS 40

float    MINV = -5.0;
float    MAXV = +5.0;
#define  STEPSV 5000
float    PI   = 3.1415927;

/* strings used by the program */
#define  NPARAM 12
struct PARAM
{
  char *word;
} SPARAM[ NPARAM ] =
{
  "Length          ",
  "Strength        ",
  "Frequency       ",
  "Criterion Used  ",
  "Criterion ROC   ",
  "W               ",
  "CO              ",
  "US              ",
  "G               ",
  "REPnew          ",
  "REPold          ",
  "TINCR           ",
};

#define  ENDWHAT 13
#define  NWHAT   ENDWHAT
struct   WHAT
{
  char *word;
  int  SWI;
} SWHAT[ NWHAT ] =
{
  "Hits                   " , 1 ,
  "False Alarms           " , 1 ,
  "Dprime                 " , 1 ,
  "Slope                  " , 1 , 
  "List average Ln(O)     " , 1 ,
  "Max Odds               " , 1 ,
  "P( F crit for list   ) " , 1 ,
  "#traces                " , 1 ,
  "P( both in 1 trace )   " , 1 ,
  "P( second in new )     " , 1 ,
  "P( second with other ) " , 1 ,
  "max images in 1 trace  " , 1 
};

struct TAG
{
  char *word;
} TAGS[ 4 ] =

{
  "               targets ",
  "           distractors ",
  "              matching ",
  "           mismatching ",
};

/* declarations of the global arrays */
float  GF[ MAXK ][ 4 ];
float    INF[ NMAX*2 ][ NWHAT  ];
float INFCRI[ NMAX*2 ][ NCRI   ];
float  SLOPE[ NLEN + NSTR + NFRE ]; 

/* declarations for the structured arrays   */
struct  VEC
{
   int  W1[ W ];
   int  W2[ W ];
   int  PO[ MAXPOS ];
   int  DR;
   int  TI;
};

struct  VEC LEX[ LMAX ];         /* lexical array */
struct  VEC LEXCOP;
struct  VEC IMA[ IMAX + NMAX ];  /* images array or memory traces */
struct  VEC TEST;                /* probe array */
           
/* structure for statistics bookkeeping */
struct  STA
{
   float SUM;
   float SUS;
   float AVE;
   float STD;
   long  N;
};

struct  STA    STATS[ NWHAT ][ NLEN + NSTR + NFRE ][ 2 ];
struct  STA STATSCRI[ NCRI  ][ NLEN + NSTR + NFRE ][ 4 ];

/* declarations of the global variables */

FILE *OUTD,*OUTFREQ;
int  o,ULEN,UFRE,USTR,whv,iv;
long n,NRUN,ntraces,testpos;
long cycle,tcycle;
long SAVEIDUM,storewhere;
char *comment;
float dprime,INTERCEPT,PMCC;

/* declarations of the function used */
void  checkparams();
void  initstats();
void  initf();
int   getf();
void  openfile();
void  myclosefile();
void  makelexicon();
void  dorecognition();
void  dorecognitionstorage();
void  doincrementstats();
void  writecycle();
void  calcstats();
float powm( float w1, int w2 );
float conv( float Pagg );
void  makeincrementlexicon();
void  regression();
void  calcslopes();
void  procstats( int whv, 
		 int iv,
		 int nstart,
		 int nend,
		 int td    );
int   wh( int whv, int iv );
int   convert( int   ustr );
void  statslength();
void  statsstrength();
void  statsfreq();
void  showparams();
void  showstats( int whv , int nv );
void  dorun();

main()
{
   SAVEIDUM = IDUM;

   checkparams();
   initf();
  
   openfile();
   initstats();   

   cycle  = 0;
   tcycle = 0;
   if (DOLENGTH) for (o=0; o<NLEN; o++) 
		     tcycle += NPOINTS / LEN[ o ];
   if (DOSTRENGTH)   tcycle += (NSTR-1) * NPOINTS / LEN[ FIXED_LEN ];
   if (DOFREQ)       tcycle += (NFRE-1) * NPOINTS / LEN[ FIXED_LEN ];
     
   /* do a length simulation */
   if (DOLENGTH)
   for (o=0; o<NLEN; o++)
   {
     ULEN = LEN[ o ];
     USTR = FIXED_STR;
     UFRE = FIXED_FRE;

     whv  = 0;
     iv   = o;
     
     NRUN = NPOINTS / ULEN;

     for (n=0; n<NRUN; n++)
     {
       ++cycle;
       writecycle();
       dorun();               
       /* separate the statistics for olds and news: 
          the test list 0..Ulen-1 are the old items and
          the test list ulen..(ulen*2)-1 are the new
          items */
       procstats( whv , iv ,  0    ,   ULEN-1   , 0 );
       procstats( whv , iv ,  ULEN , (ULEN*2)-1 , 1 );
     }
   }
                
   /* do a strength simulation */
   if (DOSTRENGTH)
   for (o=0; o<NSTR-1; o++)
   {
     ULEN = LEN[ FIXED_LEN ];
     USTR = o;
     UFRE = FIXED_FRE;

     whv  = 1;
     iv   = convert( o );    
     NRUN = NPOINTS / ULEN;

     for (n=0; n<NRUN; n++)
     {
       ++cycle;
       writecycle();
       dorun();
       
       if (USTR==0) /* pure strong */
       {
         procstats( whv , iv ,  0   , ULEN-1 , 0 );
         procstats( whv , iv , ULEN , (ULEN*2)-1 , 1 );
       }
       if (USTR==2) /* pure weak */
       {
         procstats( whv , iv ,  0 , ULEN-1 , 0 );
         procstats( whv , iv , ULEN , (ULEN*2)-1 , 1 );
       }
       if (USTR==1)  /* mixed */
       {
	   procstats( whv , iv   , 0 , (ULEN / 2)-1   , 0 );
	   procstats( whv , iv+1 , (ULEN / 2) , ULEN-1 , 0 );
           procstats( whv , iv   , ULEN   , ULEN + (ULEN/2)-1 , 1 );
           procstats( whv , iv+1 , ULEN+(ULEN/2) , (ULEN*2)-1 , 1 );
       }
       
     }
   }
        
   /* do a frequency simulation */
   if (DOFREQ)
   for (o=0; o< (NFRE-1); o++)
   {
     ULEN = LEN[ FIXED_LEN ];
     USTR = FIXED_STR;
     UFRE = o;
    
     whv  = 2;
     iv   = convert( o );
     
     NRUN = NPOINTS / ULEN;

     for (n=0; n<NRUN; n++)
     {
       ++cycle;
       writecycle();
       dorun();
       
       if (UFRE==0)  /* pure low */
       {
         procstats( whv , iv ,  0   , ULEN-1 , 0 );
         procstats( whv , iv , ULEN , (ULEN*2)-1 , 1 );
       }
       if (UFRE==2)  /* pure high */
       {
         procstats( whv , iv ,  0 , ULEN-1 , 0 );
         procstats( whv , iv , ULEN , (ULEN*2)-1 , 1 );
       }
       if (UFRE==1)  /* mixed */
       {
	   procstats( whv , iv   , 0 , (ULEN / 2)-1   , 0 );
	   procstats( whv , iv+1 , (ULEN / 2) , ULEN-1 , 0 );
           procstats( whv , iv   , ULEN   , ULEN + (ULEN/2)-1 , 1 );
           procstats( whv , iv+1 , ULEN+(ULEN/2) , (ULEN*2)-1 , 1 );
       }
       
     }
   }
    
   /* process statistics */
   calcstats();
   if (DOROC) regression();   
                       
    
   if (DOLENGTH)   statslength();
   if (DOSTRENGTH) statsstrength();
   if (DOFREQ)     statsfreq();

   myclosefile();
}
              
/* get a Z-score */
float conv( float pagg )
{
  float normf,onestep,paggc,zscore,zf;
  int   zindex;


  normf  = 1.0 / ( sqrt( 2.0 * PI ));
  paggc  = 0.0;
  zindex = 0;

  onestep = (MAXV - MINV) / (1.0*STEPSV);

  while ((paggc < pagg) && (zindex < STEPSV))
  {
    ++zindex;
    zscore = (1.0*zindex) / (1.0*STEPSV);
    zscore = zscore * (MAXV-MINV) + MINV;
    zf     = normf * exp( -0.5 * zscore * zscore );

    paggc  += zf * onestep;
  }

  zscore = (1.0*zindex) / (1.0*STEPSV);
  zscore = zscore * (MAXV-MINV) + MINV;

  return zscore;
}
     
/* regression analysis for SRN or ZROC slope estimation */
void regression()
{
  float sum1,sum2,sumsqr1,sumsqr2,value,std1,std2,av1,av2;
  int   ni1,ni2,ir,i;

  for (i=0; i< NLEN+NSTR+NFRE; i++)
  {
    sum1    = 0.0;
    sumsqr1 = 0.0;
    ni1     = 0;

    sum2    = 0.0;
    sumsqr2 = 0.0;
    ni2     = 0;
 
    for (ir=0; ir < NCRI; ir++)
    {
      value    = STATSCRI[ ir ][ i ][ 3 ].AVE;
      sum1    += value;
      sumsqr1 += value * value;
      ++ni1;

      value    = STATSCRI[ ir ][ i ][ 2 ].AVE;
      sum2    += value;
      sumsqr2 += value * value;
      ++ni2;
    }

    if ((ni1>0) && (ni2>0))
    {
      std1 = sqrt( (1.0*ni1*sumsqr1 - sum1*sum1) / (1.0*ni1*ni1) );
      av1  = sum1 / (1.0*ni1);

      std2 = sqrt( (1.0*ni2*sumsqr2 - sum2*sum2) / (1.0*ni2*ni2) );
      av2  = sum2 / (1.0*ni2);

      if (std1*std2 != 0)
      {
        PMCC = 0.0;
        for (ir=0; ir< NCRI; ir++) 
          PMCC += ( STATSCRI[ ir ][ i ][ 3 ].AVE - av1 ) * 
                  ( STATSCRI[ ir ][ i ][ 2 ].AVE - av2 ) /
            
                  ( std1 * std2 );

        PMCC        = PMCC / (1.0*NCRI);
        SLOPE[ i ]  = PMCC * ( std2 / std1 );
        INTERCEPT   = av2 - SLOPE[ i ] * av1;
      } else SLOPE[ i ] = 0.0;
    } else SLOPE[ i ]  = 0.0;
  }
}
        
/* show the parameters of the simulation */
void showparams()
{
   char *stars = "***************************";
   int  i,j;

   fprintf( OUTD , "\n" );
   fprintf( OUTD , stars );
   fprintf( OUTD , " "   );
   fprintf( OUTD , comment );
   fprintf( OUTD , " "     );
   fprintf( OUTD , stars , "\n" );
   fprintf( OUTD , "\n\n" );
   fprintf( OUTD , "# points " );
   fprintf( OUTD , "%d\n" , NPOINTS );
   fprintf( OUTD , "seed     " );
   fprintf( OUTD , "%ld\n" , SAVEIDUM );
   fprintf( OUTD , "LowG SystemG HighG: %1.3f %1.3f %1.3f\n" , LowG , SystemG , HighG );
   if (DOPAIRS) fprintf( OUTD , "PAIRS\n" ); else
		fprintf( OUTD , "SINGLES\n" );
   if (!DOINCREMENT) fprintf( OUTD , "NO INCREMENTAL STORAGE\n" ); else
                     fprintf( OUTD , "INCREMENTAL STORAGE\n" );

   for (i=0; i<NPARAM; i++)
   {
     fprintf( OUTD , "%s" , SPARAM[ i ] );
     
     if (i==0)     
     {
       if (whv==0) 
       {
	 for (j=0; j<NLEN; j++) fprintf( OUTD , "%d \t",LEN[ j ] );
       }  else 
			        fprintf( OUTD , "%d"   ,LEN[ FIXED_LEN ] );
     };

     if (i==1)
     {
       if (whv==1) 
       {
	 fprintf( OUTD , "%d \t %d \t %d \t %d" , T1 , T1 , T2 , T2 );
       } else
       {
	 if (FIXED_STR==0) fprintf( OUTD , "%d" , T1 ); else
			   fprintf( OUTD , "%d" , T2 );
       };
     }; 
     if (i==2)
     {
	if (whv==2) 
        { fprintf( OUTD , "%1.3f\t%1.3f\t%1.3f",G[ 0 ] , G[ 1 ] , G[ 2 ] );
	} else
        {
	  if (FIXED_FRE==0) fprintf( OUTD , "low" ); else
			    fprintf( OUTD , "high" );
	}
     }
     if (i==3) fprintf( OUTD , "%1.2f" , CRI[ FIXED_CRI ] );
     if (i==4) 
     {
	for (j=0; j< NCRI; j++) 
	{
	   fprintf( OUTD , "%1.2f" , CRI[ j ] );
	   fprintf( OUTD , " " );
	}
     }  
     if (i==5) fprintf( OUTD , "%d"    , W  );
     if (i==6)  fprintf( OUTD , "%1.2f" , CO );
     if (i==7)  fprintf( OUTD , "%1.2f" , US );
     if (i==8) fprintf( OUTD , "%1.3f" , G[ 1 ]  );
     if (i==9) fprintf( OUTD , "%d"    , REPnew );
     if (i==10) fprintf( OUTD , "%d"    , REPold );
     if (i==11) fprintf( OUTD , "%1.3f" , TINCR  );
  
     fprintf( OUTD , "\n" );

   }
   fprintf( OUTD , "\n" );
}
         
/* show the results of the simulation */
void showstats( int whv, int nv )
{
   int i,j,iv,ir,offset;

   if (1 > 0)
   {
     fprintf( OUTD , "AVERAGES\n" );
     for (i=0; i< NWHAT; i++)
     if ((SWHAT[i].SWI) || (i>=ENDWHAT))
     {
       if (i<ENDWHAT) fprintf( OUTD , SWHAT[ i ].word );
       if ((i==ENDWHAT) || (i==ENDWHAT)) fprintf( OUTD , "\n" );
       if  ((i==0) || (i==1) || (i==2) || (i==3) || (i>=ENDWHAT) ||
            ((i>=9) && (i<=12)) || ((i>=18) && (i<=ENDWHAT))  )
       {
	 for (iv=0; iv<nv; iv++)
	 {
	   if ((i==0) || (i==1))  fprintf( OUTD , "%1.3f\t" , 
	       STATS[ i ][ wh( whv , iv ) ][ i ].AVE );
	   if (i==2)              fprintf( OUTD , "%1.3f\t" ,
	       STATS[ i ][ wh( whv , iv ) ][ 1 ].AVE );
	   if ((i==3) && (DOROC)) fprintf( OUTD , "%1.3f\t" ,
	      SLOPE[ wh( whv , iv ) ] );
           if (i>=ENDWHAT)        fprintf( OUTD , "%1.3f\t" ,
               STATS[ i ][ wh( whv , iv ) ][ 1 ].AVE );
           if ((i>=9) && (i<=12)) fprintf( OUTD , "%1.3f\t" ,
               STATS[ i ][ wh( whv , iv ) ][ 1 ].AVE );
           if ((i>=18) && (i<ENDWHAT)) fprintf( OUTD , "%4.3f\t" ,
               STATS[ i ][ wh( whv , iv ) ][ 1 ].AVE );
	 }
	 fprintf( OUTD , "\n" );
       }
	 else
       {
	 fprintf( OUTD , "\n" );
	 if ((i>12) && (i<18)) offset = 2; else offset = 0;
         for (j=0; j<2; j++)
	 {
	    fprintf( OUTD , TAGS[ j+offset ].word );
	    for (iv=0; iv<nv; iv++)
	    {
	      fprintf( OUTD , "%1.3f\t" , 
		 STATS[ i ][ wh( whv , iv ) ][ j ].AVE );
	    }
	    fprintf( OUTD , "\n" ); 
	 }
       }  
     } 
   }

   if ((DOROC) && (SHOWZSCORES))
   {
     fprintf( OUTD , "\n" );
     for (ir=0; ir < NCRI; ir++)
     {
       fprintf( OUTD , "crit= " );
       fprintf( OUTD , "%2.2f\t" , CRI[ ir ] );
       
       for (iv=0; iv<nv; iv++)
       {
         fprintf( OUTD , "%1.3f  %1.3f\t" , STATSCRI[ ir ][ iv ][ 3 ].AVE ,
                                            STATSCRI[ ir ][ iv ][ 2 ].AVE  );
       }
       fprintf( OUTD , "\n" );
     }

     fprintf( OUTD , "\n" );
     for (ir=0; ir < NCRI; ir++)
     {
       fprintf( OUTD , "crit= " );
       fprintf( OUTD , "%2.2f\t" , CRI[ ir ] );
       
       for (iv=0; iv<nv; iv++)
       {
         fprintf( OUTD , "%1.3f  %1.3f\t" , STATSCRI[ ir ][ iv ][ 1 ].AVE ,
                                            STATSCRI[ ir ][ iv ][ 0 ].AVE  );
       }
       fprintf( OUTD , "\n" );
     }
   }

}
                
/* the following procedures are used for statistics calculations */
int wh( int whv , int iv )
{
    int offs;

    if (whv==0) offs = 0;
    if (whv==1) offs = NLEN;
    if (whv==2) offs = NLEN + NSTR;
    offs += iv;
    if ((offs <0) || (offs > NLEN+NSTR+NFRE-1 ))
    {
      printf( "something is wrong with iv\n" );
      exit(0);
    }
    return (offs);
}

void statslength()
{
   comment = "  LENGTH  "; 
   whv     = 0;

   for (iv=0; iv<NLEN; iv++)
   {
     dprime = conv( STATS[ 0 ][ wh( whv , iv ) ][ 0 ].AVE ) -
              conv( STATS[ 1 ][ wh( whv , iv ) ][ 1 ].AVE );
     STATS[ 2 ][ wh( whv , iv ) ][ 1 ].AVE = dprime;
   }

   showparams();
   showstats( whv , NLEN );
}

void statsstrength()
{
   comment = " STRENGTH ";
   whv     = 1;

   for (iv=0; iv<NSTR; iv++)
   {
     dprime = conv( STATS[ 0 ][ wh( whv , iv ) ][ 0 ].AVE ) -
              conv( STATS[ 1 ][ wh( whv , iv ) ][ 1 ].AVE );
     STATS[ 2 ][ wh( whv , iv ) ][ 1 ].AVE = dprime;
   }

   showparams();
   showstats( whv , NSTR );
}

void statsfreq()
{
   comment = "   FREQ   ";
   whv     = 2;

   for (iv=0; iv<NFRE; iv++)
   {
     dprime = conv( STATS[ 0 ][ wh( whv , iv ) ][ 0 ].AVE ) -
              conv( STATS[ 1 ][ wh( whv , iv ) ][ 1 ].AVE );
     STATS[ 2 ][ wh( whv , iv ) ][ 1 ].AVE = dprime;
   }

   showparams();
   showstats( whv , NFRE );
}

void writecycle()
{
  float percen;

  if (1 > 0)
  {
     if (tcycle >= 10)
     {
       if ((cycle % (tcycle / 10)) == 0)
       {
	  percen = 100.0 * cycle / (tcycle * 1.0);
	  printf( "Percentage done : " );
	  printf( "%3.0f\n" , percen );
	  fflush( stdout );
       }
     }
   }
}

int convert( int ust )
{
  int vout;

  vout = ust;
  if (ust==2) vout = 3; 
  return (vout);
}
               
               
void procstats( int whv, 
		int iv,
		int nstart,
		int nend,
		int td    )
{
  float sumin,susin;
  long  nin;
  int   f,i,cri_now;
  float value;

  for (f=0; f<NCRI; f++)
  {
	sumin = STATSCRI[ f ][ wh( whv , iv ) ][ td ].SUM;
	susin = STATSCRI[ f ][ wh( whv , iv ) ][ td ].SUS;
	  nin = STATSCRI[ f ][ wh( whv , iv ) ][ td ].N;

	for (i=nstart; i<=nend; i++)
	{
	   value = INFCRI[ i ][ f ];
	
	   if (value != -1000.0)
	   {
	      sumin += value;
	      susin += value*value;
	      ++nin;
	   }
	}
     
	STATSCRI[ f ][ wh( whv , iv ) ][ td ].SUM = sumin;
	STATSCRI[ f ][ wh( whv , iv ) ][ td ].SUS = susin;
	STATSCRI[ f ][ wh( whv , iv ) ][ td ].N   = nin;

  }

  for (f=0; f<NWHAT; f++)
  {
	sumin = STATS[ f ][ wh( whv , iv ) ][ td ].SUM;
	susin = STATS[ f ][ wh( whv , iv ) ][ td ].SUS;
	  nin = STATS[ f ][ wh( whv , iv ) ][ td ].N;

	for (i=nstart; i<=nend; i++)
	{
	   value = INF[ i ][ f ];
	
	   if (value != -1000.0)
	   {
	      sumin += value;
	      susin += value*value;
	      ++nin;
	   }
	}
     
	STATS[ f ][ wh( whv , iv ) ][ td ].SUM = sumin;
	STATS[ f ][ wh( whv , iv ) ][ td ].SUS = susin;
	STATS[ f ][ wh( whv , iv ) ][ td ].N   = nin;   
  }
}

void initstats()
{
   int f,i,j,k,l;

   for (f=0;f<NWHAT;f++)
     for (i=0;i<NLEN+NSTR+NFRE;i++)
       for (l=0;l<2;l++)
       {
	  STATS[ f ][ i ][ l ].SUM = 0;
	  STATS[ f ][ i ][ l ].SUS = 0;
	  STATS[ f ][ i ][ l ].AVE = 0;
	  STATS[ f ][ i ][ l ].STD = 0;
	  STATS[ f ][ i ][ l ].N   = 0;
       }

    for (f=0;f<NCRI;f++)
      for (i=0;i<NLEN+NSTR+NFRE;i++)
	for (l=0; l<4;l++)
	{
	  STATSCRI[ f ][ i ][ l ].SUM = 0;
	  STATSCRI[ f ][ i ][ l ].SUS = 0;
	  STATSCRI[ f ][ i ][ l ].AVE = 0;
	  STATSCRI[ f ][ i ][ l ].STD = 0;
	  STATSCRI[ f ][ i ][ l ].N   = 0;
	}
    
    for (i=0;i<NLEN+NSTR+NFRE;i++)
    {
       SLOPE[ i ] = 0;
    }
}

void calcstats()
{
   int f,i,j,k,l;
   float sumin,susin,ave,std;
   long  nin;

   for (f=0;f<NWHAT;f++)
     for (i=0;i<NLEN+NSTR+NFRE;i++)
       for (l=0;l<2;l++)
       {
	 sumin = STATS[ f ][ i ][ l ].SUM;
	 susin = STATS[ f ][ i ][ l ].SUS;
	 nin   = STATS[ f ][ i ][ l ].N;
	 ave   = 0.0;
	 std   = 0.0;
	 if (nin != 0)
	 {
	   ave = sumin / (1.0*nin);
	   std = sqrt( (susin * nin - sumin*sumin) /
		       (1.0*nin*nin)
		     );
	 }  
			 
	 STATS[ f ][ i ][ l ].AVE = ave;
	 STATS[ f ][ i ][ l ].STD = std;
       }
   
   for (f=0;f<NCRI;f++)
     for (i=0;i<NLEN+NSTR+NFRE;i++)
       for (l=0;l<2;l++)
       {
	 sumin = STATSCRI[ f ][ i ][ l ].SUM;
	 susin = STATSCRI[ f ][ i ][ l ].SUS;
	 nin   = STATSCRI[ f ][ i ][ l ].N;
	 ave   = 0.0;
	 std   = 1.0;
	 if (nin != 0)
	 {
	   ave = sumin / (1.0*nin);
	   std = sqrt( (susin * nin - sumin*sumin) /
					(1.0*nin*nin)
		      );
	 }  
			 
	 STATSCRI[ f ][ i ][ l ].AVE = ave;
	 STATSCRI[ f ][ i ][ l ].STD = std;

         STATSCRI[ f ][ i ][ l+2 ].AVE = conv( ave );
       }
}
      
/* to speed up sampling features, the P(particular feature) is precalculated */      
void initf()
{
   float probg;
   int   i,useg;

   for (useg=0; useg<3; useg++)
   {
     probg = G[ useg ];
     for (i=0; i<MAXK; ++i)
     {
	GF[ i ][ useg ] = probg;
	probg = probg * ( 1.0 - G[ useg ] );
     }
   }
}
   
/* sample a feature value according to g */
int  getf( int useg )
{
   int    kindex;
   float  drawr;
   float  sumg;

   do
   {
      kindex = 0;
      sumg   = 0.0;
      drawr  = ran1( &IDUM );

      while ((sumg < drawr) && (kindex < MAXK))
	 sumg = sumg + GF[ kindex++ ][ useg ];

   }
   while ( kindex >= MAXK );

   return kindex;
}
          
/* file handling */
void openfile()
{
   int error = 0;

   if (APPEND)
   {
	 if (( OUTD = fopen( FILENM , "at" )) == NULL)
	 {
		error = 1;
		OUTD = fopen( FILENM , "wt" );
	 }
   }
	 else
   {
	 OUTD = fopen( FILENM , "wt" );
   }

   if (error==1)
   {
	  printf( "WARNING: attempt to append to non-existing file\n" );
   }
}

void myclosefile()
{
   fclose( OUTD );
}

void checkparams()
{
   int i;

   for (i=0;i<NLEN;++i)
	 if (LEN[ i ] > NMAX) 
     {
        printf( "error in initializing length array\n" );
        exit(0);
     }      
   
}

/* make a lexicon of LIST-LENGTH*2 items. Later, we will store the first half, so the
first half of the lexicon are old items and the seconf half new items */
void makelexicon()
{
  long  i,ulmax,uimax,lexpos;
  int   j,dr,fr,time,t,value,useg;
  
  uimax = ULEN;
  ulmax = ULEN * 2;

  for (i=0;i<ulmax;i++)
  {
    if (UFRE==0)                 useg = 0; /* low frequency */
    if (UFRE==1)                           /* mixed frequency */
    {
       if (i < ULEN / 2)         useg = 0; else
       if (i < ULEN    )         useg = 2; else
       if (i < ULEN + ULEN / 2 ) useg = 0; else
                                 useg = 2;
    }
    if (UFRE==2)                 useg = 2; /* high frequnecy */

    for (j=0;j<W;j++)
    {
      LEX[ i ].W1[ j ] = getf( useg );
      if (DOPAIRS) LEX[ i ].W2[ j ] = getf( useg );      
    }
  }

  for (i=0; i < ULEN; i++)
  {
     lexpos = i;

     IMA[ i ].PO[ 0 ] = 1;
     IMA[ i ].PO[ 1 ] = lexpos;

     /* clean storage array */
     for (j=0; j<W; j++) 
     {
       IMA[ i ].W1[ j ] = 0;
       if (DOPAIRS) IMA[ i ].W2[ j ] = 0;
     }

     /* copy the lexical item */
     for (j=0; j<W; j++)
     {
       LEXCOP.W1[ j ] = LEX[ lexpos ].W1[ j ];
       if (DOPAIRS) LEXCOP.W2[ j ] = LEX[ lexpos ].W2[ j ];
     }


     /* do the storage */
     if (USTR==0) time = T1; 
     if (USTR==2) time = T2;
     if (USTR==1)
     {
       if (i < (ULEN / 2)) time = T1; else
                           time = T2;
     }

     /* first word */
     for (j=0; j<W; j++)
       for (t=0; t<time; t++)
       {
         if ((IMA[ i ].W1[ j ] == 0) && (ran1( &IDUM ) < US ))
	 {
            value = LEXCOP.W1[ j ];
            if (ran1( &IDUM ) > CO) value = getf( 1 );
            IMA[ i ].W1[ j ] = value;
	 }
       }
     
     /* second word */
     if (DOPAIRS)
     for (j=0; j<W; j++)
       for (t=0; t<time; t++)
       {
         if ((IMA[ i ].W2[ j ] == 0) && (ran1( &IDUM ) < US ))
	 {
            value = LEXCOP.W2[ j ];
            if (ran1( &IDUM ) > CO) value = getf( 1 );
            IMA[ i ].W2[ j ] = value;
	 }
       }
  } 
}

/* this is only used in the incremental storage version */
void makeincrementlexicon()
{
  long  i,ulmax,uimax,lexpos,totaldraws,q;
  int   j,dr,fr,time,t,value,useg;
 
  /* note: ulmax-1 is upper boundary of lexicon array */
  ulmax = ULEN * 2;
  uimax = ULEN * 2;
 
  totaldraws = 0;

  for (i=0;i<ulmax;i++)
  {
    if (UFRE==0)                 useg = 0;
    if (UFRE==1)
    {
       if (i < ULEN / 2)         useg = 0; else
       if (i < ULEN    )         useg = 2; else
       if (i < ULEN + ULEN / 2 ) useg = 0; else
                                 useg = 2;
    }
    if (UFRE==2)                 useg = 2;

    for (j=0;j<W;j++)
    {
      LEX[ i ].W1[ j ] = getf( useg );
      if (DOPAIRS) LEX[ i ].W2[ j ] = getf( useg );
    }

    if (i>=ULEN) dr = 0;
          else
    {
       if  (USTR==0)                    dr=T1;
       if ((USTR==1) && (i<  ULEN / 2)) dr=T1;
       if ((USTR==1) && (i>= ULEN / 2)) dr=T2;
       if  (USTR==2)                    dr=T2; 
    }

    totaldraws += dr;

    LEX[ i ].DR = dr;
  }

  /* clean storage array */
  for (i=0; i<uimax; i++)
  {
     IMA[ i ].PO[ 0 ] = 0;
     for (j=0; j<W; j++) 
     {
       IMA[ i ].W1[ j ] = 0;
       if (DOPAIRS) IMA[ i ].W2[ j ] = 0;
     }
  }

  ntraces = 0;
  for (q=0; q < totaldraws; q++)
  {
     do
     {
       lexpos = ran1( &IDUM ) * ULEN; 
     } while (LEX[ lexpos ].DR == 0);

     --LEX[ lexpos ].DR; 

     /* where to put the image from lexpos ? */     

     if (ntraces > 0) 
     {
       testpos = lexpos;
       dorecognitionstorage();
       time = REPold;
       if (storewhere==ntraces) 
       { 
          ++ntraces;
          time = REPnew;
	}
     }
       else   
     {
       time = REPnew;
       storewhere = 0;
       ++ntraces;
     }

     i = storewhere;

     /* do the storage */
     ++IMA[ i ].PO[ 0 ];
     if (IMA[ i ].PO[ 0 ] == MAXPOS)
     {
       printf( "\n\n\n\ warning !!!! storing too many images in one trace\n" );
       exit(0);
     }

     IMA[ i ].PO[ IMA[ i ].PO[ 0 ] ] = lexpos;
     IMA[ i ].TI = time;

     /* first word */
     for (j=0; j<W; j++)
       for (t=0; t<time; t++)
       {
         if 
            ((IMA[ i ].W1[ j ] == 0) && 
             (ran1( &IDUM ) < US ))
	 {
            value = LEX[ lexpos ].W1[ j ];
            if (ran1( &IDUM ) > CO) value = getf( 1 );
            IMA[ i ].W1[ j ] = value;
	 }
       }
     
     /* second word */
     if (DOPAIRS)
     for (j=0; j<W; j++)
       for (t=0; t<time; t++)
       {
         if 
            ((IMA[ i ].W2[ j ] == 0) && 
             (ran1( &IDUM ) < US ))
	 {
            value = LEX[ lexpos ].W2[ j ];
            if (ran1( &IDUM ) > CO) value = getf( 1 );
            IMA[ i ].W2[ j ] = value;
	 }
       }
  }

  
}

/* this is only used in the incremental storage version */
void dorecognitionstorage()
{
  float gamma,F1,F2;
  long  uimax,k,q;
  long  kl,bestpos;
  int   j,mis1,mis2,misc,mat1,mat2,matc,f,mis,mat;
  int   check,bestmis,bestnf;
  float odds,odds1,odds2,firstt1,firstt2,secondt1,secondt2;
  float lnodds,best,b,gterm,sumodds,ni;

  
  uimax = ntraces;

  b = 0.0;
  
  for (j=0; j < W; j++) 
  {
    TEST.W1[ j ] = LEX[ testpos ].W1[ j ];
    if (DOPAIRS) TEST.W2[ j ] = LEX[ testpos ].W2[ j ];
           else  TEST.W2[ j ] = 0;
  }

  /* go through all the test items */

  best    = -100.0;
  bestpos = -1;
  bestmis = 100;

  sumodds = 0.0;
  ni      = 0.0;

  for (k=0; k<uimax; k++)
  {
      misc = 0;
      matc = 0;
      mis1 = 0;
      mat1 = 0;
      mis2 = 0;
      mat2 = 0;
    
      odds1 = 1.0;
      odds2 = 1.0;

      for (j=0; j < W; j++)
      {
         if ( IMA[ k ].W1[ j ] != 0 )
         {
            if  ( TEST.W1[ j ] == IMA[ k ].W1[ j ] ) 
	    {
             
               gterm = powm( 1.0 - SystemG , TEST.W1[ j ] - 1 );
               odds1 = odds1 * ( CO + ( 1.0 - CO ) * SystemG * gterm ) / 
                       ( b * CO + ( 1.0 - b * CO ) * SystemG * gterm);
               ++mat1;
            } else
	    {
              
               odds1 = odds1 * ( 1.0 - CO ) / ( 1.0 - b * CO );
       	       ++mis1;
	    }
         }

         if ((DOPAIRS) && ( IMA[ k ].W2[ j ] != 0 ))
	 {
           if  ( TEST.W2[ j ] == IMA[ k ].W2[ j ] ) 
           {
             
              gterm = powm( 1.0 - SystemG , TEST.W2[ j ] - 1 );
              odds1 = odds1 * ( CO + ( 1.0 - CO ) * SystemG * gterm ) / 
                              ( b * CO + ( 1.0 - b * CO ) * SystemG * gterm);
              ++mat1;
            } else
	    {
             
              odds1 = odds1 * ( 1.0 - CO ) / ( 1.0 - b * CO );
      	      ++mis1;
	    }           
         }
      }   

      odds = odds1;
      mat  = mat1 + matc;
      mis  = mis1 + misc;
      
      if (odds>0.0) 
      {
         lnodds = log( odds );
   
         if ( 1.0 > 0.0 ) 
	 {
            sumodds += odds;
            ni      += 1.0;
         }

         if (lnodds > best) 
         {
           best    = lnodds;
           bestpos = k;
           bestmis = mis;
           bestnf  = mat+mis;
	 }
      }

    }

    storewhere = ntraces;
    if ((ni>0.0) && (log(sumodds/ni) > TINCR)) storewhere = bestpos;
}
      
/* main routine for calculating the odds */
void dorecognition()
{
  float gamma,F1,F2,listodds;
  long  i,uimax,ni,nilist,k,q,ulist;
  int   j,mis1,mis2,misc,mat1,mat2,matc,f,mis,mat,nfeatp,targetp,order;
  int   check,myc;
  float sumodds,odds,odds1,odds2,firstt1,firstt2,secondt1,secondt2;
  float lnodds,best,floatodds,floatlog,matchln,mismatchln;
  float matchodds,mismatchodds,Pd,Ps,gterm,odds3,odds4,term1;
  int   intlog,intodds,fp,fv;
  float b,codds,odds1b,odds2b;
  int    ch;
  
  uimax = ULEN;

  if (DOINCREMENT) uimax = ntraces;
  if (DOINCREMENT) ulist = ntraces;
              else ulist = ULEN;

  b = 0.0;  

  for (i=0; i < ULEN*2; i++)
  {
    /* construct the probe in array TEST. 
       In all simulations (in this program), we will 
       probe with a single item, so word-1 is put in TEST.W1 and
       since there is no word-2, TEST.W1 is empty */
    for (j=0; j < W; j++) 
    {
       TEST.W1[ j ] = LEX[ i ].W1[ j ];
       TEST.W2[ j ] = 0;
    }

    /* go through all the test items */

    
    sumodds = 0.0;
    listodds = 0.0;
    ni      = 0;
    nilist  = 0;
    nfeatp  = 0;
    best    = -100.0;
    targetp = 0;

    for (k=0; k<uimax; k++)
    {
      odds1  = 1.0;
      odds2  = 1.0;
      
      misc = 0;
      matc = 0;
      mis1 = 0;
      mat1 = 0;
      mis2 = 0;
      mat2 = 0;

      for (j=0; j < W; j++)
      {   
         /* only compare filled entries */
         if ( IMA[ k ].W1[ j ] != 0 )
         {
            if  ( TEST.W1[ j ] == IMA[ k ].W1[ j ] ) 
	        {
              fv    = TEST.W1[ j ];             
              /* odds calculation as described in technical report */
              gterm = powm( 1.0 - SystemG , fv - 1 );
              odds1  = odds1  * ( CO + ( 1.0 - CO ) * SystemG * gterm ) / 
                       ( SystemG*gterm);
              ++mat1;
            } else
	        {
              odds1  = odds1 *( 1.0 - CO );
	          ++mis1;
	        }
         }

         /* in the pairs simulation, the probe items will be singles
         	and the stored items will be doubles. So if we probe with
         	X and A-B is stored, then the odds of comparison X-A and
         	X-B will be calculated. The odds of all these comparisons,
         	2N in total, will be added and divided by 2N. The alignment
         	problem is only partly addressed by this approach. We do
         	not magically align X to X if X-Y is stored but we do align
         	the features of X to the features of X and the features of Y.
         */
          
         if ((DOPAIRS) && ( IMA[ k ].W2[ j ] != 0 ))
	     {
           if  ( TEST.W1[ j ] == IMA[ k ].W2[ j ] ) 
           {
              fv    = TEST.W1[ j ];
              gterm = powm( 1.0 - SystemG , fv - 1 );
              odds2 = odds2 * ( CO + ( 1.0 - CO ) * SystemG * gterm ) / 
                 ( SystemG * gterm);
              ++mat2;
           } else
	       {
             odds2  = odds2  * ( 1.0 - CO );
             ++mis2;
	       }           
         }
      }   

      if (!DOPAIRS)
      {
        odds  = odds1;
        mat   = mat1;
        mis   = mis1;
      } else
      {
        odds  = odds1 + odds2;
        mat   = mat1  + mat2;
        mis   = mat1  + mis2;
      }

      
      lnodds = log( odds );
         
     
      sumodds += odds;
      ++ni;
      if ((lnodds > best) && (k < ulist)) best = lnodds;
      if ((i<ULEN) && (k==i) && (!DOINCREMENT)) ++targetp;
            
      if (k<ulist) ++nilist;
      if (k<ulist) listodds += odds;
    }
                      
    /* -1000.0 is a flag for 'no information available'.
       In the statistics procedures, these will not count in the
       averages */
    for (f=0; f<NWHAT; f++)    INF[ i ][ f ] = -1000.0;
    for (f=0; f<NCRI ; f++) INFCRI[ i ][ f ] = -1000.0;
                                                                              
    /* ni is the number of single or double comparisons made */
    if (ni==0) /* could happen if we have a threshold (not in this program) */
    {
      lnodds   = -100;
      odds     = 0.0;
      ulist    = 1;
      nilist   = 0;
      listodds = 1.0;
    }
      else
    {
      if (DOPAIRS) ni = ni * 2.0;
      odds     = sumodds / (ni*1.0);
      lnodds   = log( odds );
    }

      INF[ i ][ 0 ] = 0.0;
      INF[ i ][ 1 ] = 0.0;
      
      if (( lnodds > CRI[ FIXED_CRI ]) && ( i <  ULEN )) INF[ i ][ 0 ] = 1.0;
      if (( lnodds > CRI[ FIXED_CRI ]) && ( i >= ULEN )) INF[ i ][ 1 ] = 1.0;

      if (nilist>0) INF[ i ][ 4 ] = log( listodds / (1.0*nilist));
      INF[ i ][ 5 ] = best;

      if (i>=ULEN) INF[ i ][ 6 ] = nilist / ( 1.0 * ulist );
      if (i<ULEN)  INF[ i ][ 7 ] = nfeatp / ( 1.0 * ulist );
      if (i<ULEN)  INF[ i ][ 8 ] = targetp;
 
      for (f=0; f<NCRI; f++)
      {
        INFCRI[ i ][ f ] = 0.0;
        if (lnodds > CRI[ f ]) INFCRI[ i ][ f ] = 1.0;
      }
    }
  
}


float powm( float w1, int w2 )
{
  float poww;
  int   wc;

  poww = 1.0;
  wc   = w2;
  while (wc != 0)
  {
    --wc;
    poww = poww * w1;
  }

  return poww;
}
                                    
/* calculate some statistics for the incremental storage version */
void doincrementstats()
{
  long i,j,k;
  int  indiffnew,indiffother,inone,prevj,firstother;
  int  tracesin[ 5 ];
  int  maxtraces,numbert;

  for (i=0; i<ULEN*2; i++ ) INF[ i ][ 18 ] = ntraces;    

  tracesin[ 1 ] = 0;
  tracesin[ 2 ] = 0;
  tracesin[ 3 ] = 0;
  tracesin[ 4 ] = 0;

  maxtraces     = 0;

  for (i=0; i<ULEN; i++)
  {
    inone       = 0;
    indiffnew   = 0;
    indiffother = 0;
    firstother  = 0;

    prevj       = -1;

    for (j=0; j<ntraces; j++)
    {
      numbert = IMA[ j ].PO[ 0 ];
      if (numbert > maxtraces) maxtraces = numbert;
 
      for (k=1; k<=numbert; k++)
      {
        if (IMA[ j ].PO[ k ]==i) 
	{
           if (prevj == -1)
	   {
             prevj = j;
             if (k>1) firstother=1; 
           } 
             else 
               if (prevj == j )
	       {
                 inone = 1;
               }
                 else
                   if (prevj != j)
	           {
                      if (k==1) indiffnew   = 1; else
                                indiffother = 1;
                   }
	}
      }       
    }

    if (inone==1)       ++tracesin[ 1 ];
    if (indiffnew==1)   ++tracesin[ 2 ];
    if (indiffother==1) ++tracesin[ 3 ];
    if (firstother==1)  ++tracesin[ 4 ];
  }

  INF[ i ][ 19 ] = tracesin[ 1 ] / (1.0*ULEN); /* both in one     */
  INF[ i ][ 20 ] = tracesin[ 2 ] / (1.0*ULEN); /* second in new   */
  INF[ i ][ 21 ] = tracesin[ 3 ] / (1.0*ULEN); /* second in other */
  INF[ i ][ 22 ] = firstother;
}
     
/* main running routine */
void dorun()
{
  int f,i;

  if (DOINCREMENT) makeincrementlexicon();
             else  makelexicon();

  dorecognition();
   
  if (DOINCREMENT) doincrementstats();
}