! Exemplar classification single-layer Perceptron
! Uses delta rule to learn patterns from two categories
! Input representations are Gaussian densities
! Luce choice rule is used to compute probability from output node activation
! Ouput files: p_corr: p(correct) classification over epochs
!              w_corr: activation at correct out node over epochs
!	     w_incorr: activation at incorrect out node over epochs

PROGRAM Exemplar_Categorization

USE number_generators
 
implicit none

  integer, parameter    :: N_Feat = 100,  &    ! Features in representation of exemplar
                           N_Output = 2,  &    ! Output response categories   
			   N_REPS = 100000
  
  real,    parameter    :: Alpha  = 0.01, &    ! Learning rate parameter for delta rule
                               e  = 2.71828, & ! Euler's number (constant)
			       b  = 2.0        ! Sensitivity parameter for Luce choice rule
			   
  real                  :: Prototype(2,N_Feat), &
                           Target(2,N_Feat),    &     
                           Weight(N_Output, N_Feat),& 
			   Output(N_Output) ,   &
                           Activation(N_Output),& 
			   x
  			   		   		
  integer               :: i, j, epoch, targ_type, resp, idum			   
			   
!++++++++BEGIN SIMULATION++++++++++++++++++++++++++++++++++++++
  call Init_Files
  call RandSeed
  call Generate_Vectors(Prototype(1,:))   ! Parent pattern for category 1
  call Generate_Vectors(Prototype(2,:))   ! Parent pattern for category 2
  call Initialize_Weights(Weight)         ! Start with uniform random weights
 
  do epoch = 1, N_REPS
    ! Flip a coin and pick a category to present exemplar from:
    x = ran3(idum)
    if (x < 0.5) then
      targ_type = 2
      Target(2,:) = Add_Noise(Prototype(2,:), 10, 0.0, 0.2)
    else
      targ_type = 1
      Target(1,:) = Add_Noise(Prototype(1,:), 10, 0.0, 0.2)
    endif
    
    call Classify(Target(targ_type,:), Weight, Activation, targ_type)
    call Update_Weights(Target(targ_Type,:), Weight, Alpha, Activation, targ_type)

  enddo 

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

!++++++++++END SIMULATION++++++++++++++++++++++++++++++++++++++

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

!************************************************************************
subroutine Init_Files

open(unit=1, file='p_corr', status='replace')
open(unit=2, file='w_corr', status='replace') 
open(unit=3, file='w_incorr', status='replace')

end subroutine Init_Files
!************************************************************************


!************************************************************************
subroutine Classify (Target, Weight, Activation, targ_type)
   integer :: i, j, targ_type, correct, incorrect
      real :: Target(N_Feat), Weight(N_Output,N_Feat), Activation(N_Output), prob(N_Output)

  correct = targ_type
  if (correct==1) then 
    incorrect = 2
  else
    incorrect = 1
  endif

  Activation = 0.0
  Activation(correct) = dot_product(Target, Weight(correct,:))
  Activation(incorrect) = dot_product(Target, Weight(incorrect,:))
  
  ! Luce's choice rule: 
  prob(correct) = (e**(b*Activation(correct))) / ((e**(b*Activation(correct)))+(e**(b*Activation(incorrect))))

 write(1,*) prob(correct)
 write(2,*) Activation(correct)
 write(3,*) Activation(incorrect)
 
end subroutine Classify
!************************************************************************

!************************************************************************
subroutine Update_Weights(Target, Weight, Alpha, Activation, targ_type)
   integer :: i, targ_type
   real    :: Alpha, Weight(N_Output,N_Feat), Activation(N_Output), &
              Correct(N_Output), Target(N_Feat)
  
  if (targ_type==1) then
    Correct(1) = 1.0
    Correct(2) = 0.0
  else
    Correct(1) = 0.0
    Correct(2) = 1.0
  endif
 
  ! Delta Update Rule:
  do i = 1, N_Output
    do j = 1, N_Feat
      Weight(i,j) = Weight(i,j) + (Alpha*(Correct(i)-Activation(i))*Target(j))
    enddo    
  enddo

end subroutine Update_Weights
!************************************************************************

!****************************************************************
subroutine Generate_Vectors(X)
  real     :: X(N_Feat)
  real     :: val
 integer   :: i, j

 val = gaussian(0.0, 0.2)
  do j = 1,N_Feat
   X(j) = gaussian(0.0,0.2)
  enddo

end subroutine Generate_Vectors
!****************************************************************


!************************************************************************
subroutine Initialize_Weights(Weight)
   integer :: i, j
   real    :: Weight(N_Output, N_Feat), x, y
   
 do i = 1, N_Output
   do j = 1, N_Feat
	 x = ran3(idum)
	 y = ran3(idum)
	 if (y < 0.5) x = -1.0*x   ! flip bit
   	 Weight(i,j) = x
    enddo
 enddo

end subroutine Initialize_Weights
!************************************************************************

!****************************************************************
function Add_Noise(Vector, N_Flip, Mu, Sigma)
  integer :: N_Flip, i, j, i_dum
  real:: Vector(N_Feat), Temp(N_Feat), Add_Noise(N_Feat), Mu, Sigma

 Temp = Vector
 do i = 1, N_Flip
  j = int(ran3(i_dum)*N_Feat)+1
  Temp(j) = gaussian(Mu, Sigma)  
 enddo
 Add_Noise = Temp

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


END PROGRAM Exemplar_Categorization