! Implementation of our Decay model from class. 
! Memory traces decay as an exponential function of time
! This simulates the Waugh-Normal paradigm

PROGRAM Decay_Model

use number_generators

implicit none

 integer, parameter  :: D = 100,         &  ! Vector dimensionality
                        N = 16              ! # of items to store
			
   real, parameter   :: Lambda = 1.0,    &  ! Exponential rate parameter
                             e = 2.71828,&  ! Euler's number (constant)
                        Tick   = 0.025      ! Time increment parameter
 
		real :: Prototype(D),    &  
        		Memory(N,D),     &
			Item(N,D),       &
			t_encode(N),     &
			strength, t

	     integer :: i, j, idum

!++++++++BEGIN SIMULATION++++++++++++++++++++++++++++++++++++++
 call RandSeed                 ! initializing the r-num generator
 call Init_Output              ! initializing files for time-series output
 
 Prototype = Random_Vector(D)  ! Create parent item
 
 do i = 1, N
   Item(i,:) = Distort(Prototype, D, 50)  
 enddo

 t = 0
 do i = 1, N
  t = t + tick
  Memory(i,:) = Item(i,:)
  t_encode(i) = t
  Memory = Decay(Memory, t_encode, t, i)   
 enddo

 do i = 1, N
 strength = 0
  do j = 1, D
    strength = strength + abs(Memory(i,j))
  enddo
  write(*,*) strength/D 
 enddo

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


!****************************************************************
function Decay(Memory, t_encode, time, pos)
   integer :: i, j, pos
      real :: time, t_delta, t_encode(N), Memory(N,D), t(N), Decay(N,D)
 
  do i = 1, pos
    t_delta = time - t_encode(i)
    Memory(i,:) = Memory(i,:) * (e**(-1*(t_delta*Lambda)))
  enddo

  Decay = Memory

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

!****************************************************************
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
!****************************************************************

END PROGRAM Decay_Model