! Implementing theories of forgetting: Decay vs. Interference
! Both models are exactly the same except for the force that operates
! on memory over time. With Decay, memory traces simply decay back to 0
! as a function of time. With Interference, old memory traces have features
! shared with new items overwritten with p(lose) (cf. Nairne feature model)
! Decay is sensitive to time, whereas Interference is sensitive to # and similarity
! of items stored (which just happens to correlate w/ time). This simulates the Waugh-
! Norman paradigm presenting 16 items but varying presentation speed. 
! Output files are: decay1, decay4, and interfere

PROGRAM Forget

use number_generators

implicit none

 integer, parameter  :: D = 100,         &  ! Vector dimensionality
                        N = 16,          &  ! # of items to store
		   N_Flip = D/2             ! # of features to bit-flip from parent
			
   real, parameter   :: Lambda = 1.0,    &  ! Exponential rate parameter
                             e = 2.71828,&  ! Euler's number (constant)
		        P_Lose = 0.5        ! Prob of new item overwriting feature share w exisiting item
 
		real :: Prototype(D),    &  
        		Memory(N,D),     &
			Item(N,D),       &
			t_encode(N),     &
			strength, t,     &
			Tick = 0.5          ! Time increment parameter 

	     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)
 
 do i = 1, N
   Item(i,:) = Distort(Prototype, D, N_Flip)
 enddo

! ***************************************
! Running a trial with Decay; 1 item/sec:
! *************************************** 
 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(1,*) strength/D 
 enddo


! ***************************************
! Running a trial with Decay: 4 items/sec:
! *************************************** 
 Tick = Tick/4
 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(2,*) strength/D 
 enddo


! ***************************************
! Same items, but with interference:
! ***************************************
 t = 0
 do i = 1, N
  t = t + tick
  Memory(i,:) = Item(i,:)
  Memory = Interfere(Memory, i)
  
 enddo

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


 close(1)
 close(2)
 close(3)
 

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


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

  open(unit=1, file='decay1', status='replace')
  open(unit=2, file='decay4', status='replace')
  open(unit=3, file='interfere', status='replace')

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

!****************************************************************
function Interfere(Memory, pos)
   integer :: i, j, pos
      real :: Memory(N,D), New_Item(D), Interfere(N,D), x
 
  if (pos==1) then
    Interfere = Memory
    return
  endif
  
  New_Item = Memory(pos,:)

  do i = 1, (pos-1)
    do j = 1, D
      if (Memory(i,j)==New_Item(j)) then 
       x = ran3(idum)
       if (x < P_Lose) Memory(i,j) = 0
      endif
    enddo
 enddo

 Interfere = Memory

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

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


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