! This is an example of a one-layer perceptron letter detector. 
! We will train it to discriminate between two letters
! We'll use the feature definitions for letters from Rumelhart & McClelland

PROGRAM Perceptron

USE number_generators
 
implicit none

  integer, parameter    :: N_Feat = 14,  &  ! For single-layer, make sure nodes=features
                           N_Node = 14,  &
			   N_REPS = 100
  real,    parameter    :: Alpha  = 0.01
			   
			   
  real                  :: Target(2, N_Feat), &    ! Input feature representation for letter
                           Weight(N_Node)   , &    ! Weights applied to neural connections
			   Activation       , &    ! Summed activation at THD unit
			   x
			   
  logical               :: THD,               &    ! Boolean threhold device--on or off
                           Gotit(2)=.false.,  & 
			   N_Correct(2)=.false.,&  ! Correct on this trial for type?
			   L_Correct(2)=.false.    ! Correct on last trial for type?
			   		   		
  integer               :: i, j, epoch, targ_type, resp, idum			   


			   
			   
!++++++++BEGIN SIMULATION++++++++++++++++++++++++++++++++++++++
  call RandSeed
  
  
  !call Get_Features(Target)
  call Generate_Vectors(Target(1,:))
  call Generate_Vectors(Target(2,:))
  
  call Initialize_Weights(Weight)

  epoch = 0
  do 
    epoch = epoch+1
    targ_type = 1
    x = ran3(idum)
    if (x < 0.5) targ_type = 2
    N_Correct(targ_type)=.false.
    call THD_Activation(Target(targ_type,:), Weight, Activation, resp)
    
    write(*,*) 'Target/Resp/Activation ', targ_type, resp, Activation
  
    if (resp <> targ_type) then 
      call Update_Weights(Weight, Alpha, Activation, targ_type)
    else
      N_Correct(targ_type) = .true.
    endif
    
    if(N_Correct(1).and.N_Correct(2).and.L_Correct(1).and.L_Correct(2)) exit
    L_Correct(targ_type) = N_Correct(targ_type)
    
  enddo 

  write(*,*) epoch

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

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

!=== get_features ======================================================
!      Reads a table of input patterns (14 features for each letter)

!************************************************************************
subroutine THD_Activation(Target, Weight, Activation, resp)
   integer :: i, resp
   real    :: Target(N_Feat), Weight(N_Node), Activation

 Activation = 0.0
 Activation = dot_product(Target, Weight)
 resp = 1
 if (Activation > 0) resp = 2

end subroutine THD_Activation
!************************************************************************

!************************************************************************
subroutine Update_Weights(Weight, Alpha, Activation, targ_type)
   integer :: i, targ_type
   real    :: Alpha, Weight(N_Node), Activation, true
   
 if (targ_type == 1) then
   true = -1.0
 else
   true = 1.0
 endif
 
 
 do i = 1, N_Node
   Weight(i) = Weight(i) + (Alpha*(true-Activation)*Target(targ_type,i))
 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
   val = ran3(idum)
   X(j) = ran3(idum)
   if (val > 0.5) X(j) = -1*X(j)
   
   !X(j) = gaussian(0.0, 0.2)
  enddo

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

!************************************************************************
subroutine Get_Features(Target)
   integer   :: i, j
   real      :: Target(2,N_Feat)
   character :: dum
   
      open (UNIT = 1, file = 'ltrdefs', status='old' )
      read(1, '(5(/))')
      do i = 1, 2
          read(1, ' (A1, F4.0, 13F5.0) ') dum, (Target(i,j), j=1,14)
          write(*, '(A4, 14F5.0)') dum, Target(i,:)
      enddo
    close(1)
end subroutine Get_Features
!************************************************************************

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

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

END PROGRAM Perceptron