! Using Minerva to illustrate a process dissociation between recognition and 
! classification with delay from a single memory stores
! Output is p(correct) for each task as a function of test 
! delay (files: 'class' and 'recog')

PROGRAM Minerva2

use number_generators

implicit none


 integer, parameter  :: N_Feat = 20, &       ! # of features in each field
                        N_Cats = 4,  &       ! # of categories
		     N_Per_Cat = 8,  &       ! # exemplars per category
			   Exp = 3,  &       ! # power of similarity metric
		       N_Items = N_Per_Cat*N_Cats, &  ! total # of exemplars to be learned
    	             N_Distort = 5,  &       ! # of features to flip for exemplars
		    Crit_Recog = 3.0,&       ! Recognition criterion
		    Crit_Cat   = .1, &       ! Min Categorization criterion 
		         SUBS  = 2000        ! # of subjects (replications) to simulate
			   
   real              ::      F = 0.4   ! probability of forgettig a feature

   real  :: Memory(N_Feat*2, N_Items), &
            Trace(N_Feat*2, N_Items),  &
            Prototype(N_Feat, N_Cats), &
	    Label(N_Feat, N_Cats),     &
	    Probe(N_Feat),             &
            Echo(N_Feat),              &
	    Activation,                &
	    Correct(2),                &
	    intensity             
	    
 integer :: i, j, idum, pos, sample, time, Cat=1

!++++++++BEGIN SIMULATION++++++++++++++++++++++++++++++++++++++
 call RandSeed                 ! initializing the r-num generator
 call Init_Output              ! initializing files for time-series output
 
 DO time = 1, 6                ! How many time points? 
   
   F = F +0.05                 ! At each time point, p(forget) increases by this amount
   Correct = 0
  
   DO sample = 1, SUBS         ! For each time point, simulate this many subjects (new stimuli and memories)
      pos = 0
       ! Generating Prototypes and Labels for N_Cats categories: 
      do i = 1, N_Cats
       Prototype(:,i) = Random_Vector(N_Feat)
       Label(:,i) = Random_Vector(N_Feat)
      enddo
      Memory = 0
    
      ! BUILDING MEMORY:   
      do i = 1, N_Cats
        do j = 1, N_Per_Cat
	  pos = pos+1
	  Memory(1:N_Feat,pos) = Distort(Prototype(:,i), N_Feat, N_Distort)
	  Memory(N_Feat+1:N_Feat*2,pos) = Label(1:N_Feat,i)	
	enddo
      enddo
      
    call Forget(F,Memory)
    
    ! CLASSIFYING A PATTERN:
    probe = Memory(1:N_Feat, 2)      ! Using  instance of Category 1
    Correct(1) = Correct(1) + Classify(Probe, Label, Memory, Cat)           
    
    ! RECOGNIZING A PATTERN: 
    intensity = Echo_Intensity(Probe, Memory)
    if (intensity > Crit_Recog) Correct(2) = Correct(2) + 1
    
  ENDDO ! sample
   
   do i = 1, 2
      write(i,*) time, ' ', Correct(i)/SUBS   
   enddo
    
 ENDDO ! time
 
 close(1)
 close(2)
 
!++++++++END SIMULATION++++++++++++++++++++++++++++++++++++++++

 CONTAINS
! Subroutines follow this statement
!=================================================================

!****************************************************************
subroutine Init_Output

  open(unit=1, file='class', status='replace')
  open(unit=2, file='recog', status='replace')

end subroutine Init_Output
!****************************************************************

!****************************************************************
function Random_Vector(D)
 integer  :: D, i
 real     :: Random_Vector(D), x

do i = 1, D
  x = ran3(idum)  ! random number 0 <= x <= 1 from uniform 
  if  (x >= 0.5) then
    Random_Vector(i) = 1    
  else
    Random_Vector(i) = -1
  endif
enddo
end function Random_Vector
!****************************************************************

!****************************************************************
function Distort(Vector, D, n_flip)
  integer :: i, j, D, n_flip, List(n_flip)
  real    :: Vector(D), Distort(D), x

 List = Rand_List(n_flip, D)        ! get a list of n_flip random numbers from 1..D
 Distort = Vector

do i = 1, n_flip
 Distort(List(i)) = -1.0 * Vector(List(i))   ! flipping bit
enddo

end function
!****************************************************************

!*****************************************************************
subroutine Forget(F, Matrix)
  integer  :: i, j
  real     :: F, Matrix(N_Feat*2,N_Items), x
  
 do i = 1, N_Items
   do j = 1, N_Feat*2
      x = ran3(idum) 
      if (x < F) Matrix(j,i) = 0   ! Forget the feature with a probability of F
   enddo
 enddo

end subroutine Forget
!*****************************************************************

!*****************************************************************
function Echo_Intensity(Probe, Memory)
   integer :: i, j 
   real    :: Echo_Intensity, Probe(N_Feat), Echo(N_Feat), Memory(N_Feat*2, N_Items)

 Echo = 0
 Echo_Intensity = 0
 
 do i = 1, N_Items
   Activation = Vector_Cosine(Probe, Memory(1:N_Feat,i), N_Feat) ** Exp
   Echo = Echo + (Memory(1:N_Feat,i)*Activation)
 enddo

 do i = 1, N_Feat
   Echo_Intensity = Echo_Intensity + abs(Echo(i))
 enddo

end function Echo_Intensity
!*****************************************************************

!*****************************************************************
function Classify(Probe, Label, Memory, Correct_Label)
  integer :: i, j, Correct_Label, Assigned_Label, Classify, loc
  real    :: Probe(N_Feat), Label(N_Feat,N_Cats), Memory(N_Feat*2,N_Items), &
             Similarity, Max 

 Classify = 0
 echo = 0
 Max = 0.0
 
 do i = 1, N_Items
  Activation = Vector_Cosine(Probe, Memory(1:N_Feat,i), N_Feat) ** Exp
  Echo = Echo + (Memory(N_Feat+1:N_Feat*2,i)*Activation) 
 enddo
 
 do i = 1, N_Cats
   Similarity = Vector_Cosine(Echo, Label(:,i), N_Feat)
   if (Similarity > Max) then 
      Max = Similarity
      loc = i
   endif
 enddo

 Assigned_Label = loc
 if ((Correct_Label == Assigned_Label).and. (Max>Crit_Cat)) Classify = 1         ! if we pick the right label, give us a point

end function
!*****************************************************************



!************************************************
!     SOME MATRIX ALGEBRA TOOLS:
!
!************************************************
!****************************************************************
function Vector_Length (Vector, n)
   integer :: i, n 
   real    :: Vector_Length, Vector(n), SS

 SS = 0.0
  SS = dot_product(Vector, Vector)
 Vector_Length = sqrt(SS)
end function Vector_Length 
!****************************************************************


!****************************************************************
subroutine Normalize (Vector_In, n)
   integer             :: i, n
   real                :: Vector_In(n), Vect_Length, Test_Val
  
  Test_Val = dot_product(Vector_In, Vector_In)
  if (Test_Val == 0.0) then 
    return
  endif
   
  Vect_Length = 0.0 
  Vect_Length = Vector_Length(Vector_In, n)
  Vector_In = Vector_In / Vect_Length
end subroutine Normalize
!****************************************************************


!****************************************************************
function Vector_Cosine (Vector1, Vector2, n)
   
   real      :: Vector1(n), Vector2(n)
   real      :: Vector_Cosine
   integer :: n, i
 
 call Normalize(Vector1, n) 
 call Normalize(Vector2, n) 
 Vector_Cosine = dot_product(Vector1, Vector2)
 
end function Vector_Cosine
!****************************************************************


END PROGRAM Minerva2