!/ ====================================================================== BEGIN FILE =====
!/ **                                 T E S T _ I R I S                                 **
!/ =======================================================================================
!/ **                                                                                   **
!/ **  Copyright (c) 2016, Stephen W. Soliday                                           **
!/ **                      stephen.soliday@trncmp.org                                   **
!/ **                      http://research.trncmp.org                                   **
!/ **                                                                                   **
!/ **  -------------------------------------------------------------------------------  **
!/ **                                                                                   **
!/ **  This program is free software: you can redistribute it and/or modify it under    **
!/ **  the terms of the GNU General Public License as published by the Free Software    **
!/ **  Foundation, either version 3 of the License, or (at your option)                 **
!/ **  any later version.                                                               **
!/ **                                                                                   **
!/ **  This program is distributed in the hope that it will be useful, but WITHOUT      **
!/ **  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS    **
!/ **  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.   **
!/ **                                                                                   **
!/ **  You should have received a copy of the GNU General Public License along with     **
!/ **  this program. If not, see <http://www.gnu.org/licenses/>.                        **
!/ **                                                                                   **
!/ =======================================================================================
!
!> @brief   Test Cluster
!!
!! @author  Stephen W. Soliday
!! @date    2016-04-20
!!
!! @details This program will use the means and inverse covariances in a configuration
!!          file to classify a table of data and then compare the labels to known values.
!!          This code has been extracted and simplified from a FORTRAN module I wrote
!!          that supports both supervised and unsupervised learning of Gaussian clusters.
!!          That module is part of a research project, that is on going. At the end I
!!          may release the entire Machine Learning library that it is part of.
!!
!!  The configuration file begins with two integers on the first line: number of classes
!!  and number of variables. For the Iris data this is 3 4
!!  the next line contains the means for the first cluster.
!!  the next four lines contain the lower triangle of the inverse covariance for the first
!!  cluster. This is repeated twice more for the second and third clusters.
!!
!!  The data file begins with two integers on the first line: number of samples and
!!  number of variables. For the Iris data this is 150 4
!!  each line starts with an integer label followed by floating-point variables.
!!
!!  This program will report then number of correct matches, and list the incorrect line:
!!
!!  Example:
!!           3 != sample(73): 2 |  6.30  2.50  4.90  1.50 
!!
!!     This indicates that sample number 73 which has the values (6.30  2.50  4.90  1.50)
!!     should have been labeled 2 but was incorrectly labeled 3.
!
!/ =======================================================================================
module test_iris
  !/ -------------------------------------------------------------------------------------
  use iso_fortran_env
  implicit none

  integer,  parameter :: dp=kind(1.0d0)
  real(dp), parameter :: N_ZERO=1.0d0

  !/ =====================================================================================
  type :: m_center
     !/ ----------------------------------------------------------------------------------
     integer               :: nVar = 0      !< number of variables
     real(dp), allocatable :: mean(:)       !< cluster centers
     real(dp), allocatable :: inv_cov(:,:)  !< inverse covariance work array
   contains
     procedure :: metric => compute_mcenter_metric
  end type m_center

  !/ =====================================================================================
  type :: m_class
     !/ ----------------------------------------------------------------------------------
     integer                     :: nClass = 0   !< number of classes
     integer                     :: nVar   = 0   !< number of variables
     type(m_center), allocatable :: cluster(:)   !< array of m_centers
     real(dp),       allocatable :: temp_row(:)  !< work space
   contains
     procedure, private :: label_single
     procedure, private :: label_list
     generic            :: label => label_single, label_list
  end type m_class


  interface create
     module procedure :: create_mcenter
     module procedure :: create_mclass
  end interface create

  interface delete
     module procedure :: delete_mcenter
     module procedure :: delete_mclass
  end interface delete




  !/ =====================================================================================
contains !/**                     P R O C E D U R E   S E C T I O N                     **
  !/ =====================================================================================




  !/ =====================================================================================
  !> @brief Mahalanobis distance metric.
  !! @param[in] A metric.
  !! @param[in] x first vector.
  !! @param[in] y second vector.
  !! @return square of the Mahalanobis distance metric.
  !!
  !! Compute the Mahalanobis distance metric. A is a square positive semi-definite matrix.
  !! Typically this will be the inverse of a covariance matrix. If A is the Identity
  !! matrix then this metric reduces to Euclidean distance.
  !!
  !! @note this function returns the square for the distance. In most comparisons for 
  !! relative distance calling SQRT increases the computation time but does not buy you
  !! anything. If a < b then sqrt(a) < sqrt(b) for all a,b in REAL
  !/ -------------------------------------------------------------------------------------
  function Mahalanobis( A, x, y ) result( d )
    !/ -----------------------------------------------------------------------------------
    real(dp)             :: d
    real(dp), intent(in) :: A(:,:)
    real(dp), intent(in) :: x(:)
    real(dp), intent(in) :: y(:)
    !/ -----------------------------------------------------------------------------------
    integer  :: i,j,n
    real(dp) :: di, dj
    n = size(x)
    d = N_ZERO
    do i=1,n
       di = x(i)-y(i)
       d  = d + A(i,i)*di*di
       do j=i+1,n
          dj = x(j)-y(j)
          d  = d + 2.0d0*A(j,i)*di*dj
       end do
    end do
  end function Mahalanobis




  !/ =====================================================================================
  !> @brief Create.
  !! @param[in,out] mc     reference to this object.
  !! @param[in]     unit   file unit.
  !! @param[in]     status optional return status when using unit.
  !!
  !! Create a cluster center. The configuration for this center will be read from a
  !! file unit. If status is present and unit was used, then a return status may be
  !! provided.
  !!
  !! @note vars and unit may not be used together. Using status with vars has no effect.
  !/ -------------------------------------------------------------------------------------
  subroutine create_mcenter( mc, nv, unit, status )
    !/ -----------------------------------------------------------------------------------
    class(m_center), intent(inout) :: mc
    integer,         intent(in)    :: nv
    integer,         intent(in)    :: unit
    integer,         intent(out)   :: status
    !/ -----------------------------------------------------------------------------------
    integer :: i, j

    status = 0

    if ( allocated( mc%mean ) ) deallocate( mc%mean )
    allocate( mc%mean(nv) )

    if ( allocated( mc%inv_cov ) ) deallocate( mc%inv_cov )
    allocate( mc%inv_cov(nv,nv) )

    read(unit,*,end=800,err=900) ( mc%mean(i), i=1,nv )

    do i=1,nv
       read(unit,*,end=800,err=900) ( mc%inv_cov(i,j), j=1,i )
       do j=1,i
          if ( i.ne.j ) then
             mc%inv_cov(j,i) = mc%inv_cov(i,j)
          end if
       end do
    end do

    goto 1000

    !/ -----------------------------------------------------------------------------------
800 continue
    write(*,'(A)') '**E** create_mcenter: Premature EOF encountered'
    status = 1
    goto 1000

900 continue
    write(*,'(A)') '**E** create_mcenter: I/O error occurred'
    status = 2

1000 continue
  end subroutine create_mcenter

  !/ =====================================================================================
  !> @brief Delete.
  !! @param[in,out] mc reference to this object.
  !!
  !! Deallocate allocations.
  !/ -------------------------------------------------------------------------------------
  subroutine delete_mcenter( mc )
    !/ -----------------------------------------------------------------------------------
    class(m_center), intent(inout) :: mc
    !/ -----------------------------------------------------------------------------------
    mc%nVar = 0
    if ( allocated( mc%mean ) )    deallocate( mc%mean )
    if ( allocated( mc%inv_cov ) ) deallocate( mc%inv_cov )
  end subroutine delete_mcenter

  !/ =====================================================================================
  !> @brief Metric.
  !! @param[in,out] this reference to this object.
  !! @param[in]     tmean reference to a test vector.
  !! @return the test metric.
  !!
  !! Compute a similarity match between the test vector and this center's means.
  !/ -------------------------------------------------------------------------------------
  function compute_mcenter_metric( this, tmean ) result( M )
    !/ -----------------------------------------------------------------------------------
    real(dp)                       :: M
    class(m_center), intent(inout) :: this
    real(dp),        intent(in)    :: tmean(:)
    !/ -----------------------------------------------------------------------------------
    ! M = Euclid( this%inv_cov, this%mean, tmean )
    M = Mahalanobis( this%inv_cov, this%mean, tmean )
  end function compute_mcenter_metric








  !/ =====================================================================================
  !> @brief Create.
  !! @param[in,out] mc      reference to this object.
  !! @param[in]     fspc    path to configuration file.
  !! @param[in]     status  optional return status when using unit.
  !!
  !! Create a classifier. If status is present and file was used, then a return status
  !! may be provided.
  !/ -------------------------------------------------------------------------------------
  subroutine create_mclass( mc, fspc, status )
    !/ -----------------------------------------------------------------------------------
    class(m_class),             intent(inout) :: mc
    character(len=*),           intent(in)    :: fspc
    integer,          optional, intent(out)   :: status
    !/ -----------------------------------------------------------------------------------
    integer :: i, nc, nv, un, ios, rv

    rv = 0

    open( NEWUNIT=un, FILE=fspc, ACTION='READ', STATUS='OLD', IOSTAT=ios )
    if ( 0.ne.ios ) then
       write(*,110) fspc, ios
       rv=ios
       goto 1100
    end if

    !/ -----------------------------------------------------------------------------------

    if ( allocated( mc%cluster ) ) then
       do i=1,mc%nClass
          call delete( mc%cluster(i) )
       end do
       deallocate( mc%cluster )
    end if

    read(un,*,end=800,err=900) nc, nv

    allocate( mc%cluster(nc) )

    mc%nClass = nc
    mc%nVar   = nv

    do i=1,nc
       call create( mc%cluster(i), mc%nVar, un, ios )
       if ( 0.ne.ios ) then
          write(*,120) i
          rv=ios
          goto 1000
       end if
    end do

    write(*,130) nc, nv

    if ( allocated( mc%temp_row ) ) deallocate( mc%temp_row )
    allocate( mc%temp_row(mc%nVar) )

    goto 1000
    !/ -----------------------------------------------------------------------------------

110 format( 'Cannot open file:[',A,'] for reading, errno:',I0)
120 format( 'create_mclass: Cluster(',I0,') failed to read')
130 format( 'Configured ',I0,' classes of ',I0,' variables' )

800 continue
    write(*,'(A)') '**E** create_mclass: Premature EOF encountered'
    rv = 1
    goto 1000

900 continue
    write(*,'(A)') '**E** create_mclass: I/O error occurred'
    rv = 2

1000 continue
    close(un)

1100 continue
    if ( present( status ) ) then
       status = rv
    end if
  end subroutine create_mclass

  !/ =====================================================================================
  !> @brief Delete.
  !! @param[in,out] mc reference to this object.
  !!
  !! Remove any allocations and zero out buffers.
  !/ -------------------------------------------------------------------------------------
  subroutine delete_mclass( mc )
    !/ -----------------------------------------------------------------------------------
    class(m_class), intent(inout) :: mc
    !/ -----------------------------------------------------------------------------------
    integer :: i

    if ( allocated( mc%cluster ) ) then
       do i=1,mc%nClass
          call delete( mc%cluster(i) )
       end do
       deallocate( mc%cluster )
    end if

    if ( allocated( mc%temp_row ) ) deallocate( mc%temp_row )

  end subroutine delete_mclass

  !/ =====================================================================================
  !> @brief Label.
  !! @param[in,out] this reference to this object.
  !! @param[out]    label class that data fits.     
  !! @param[in]     data  sample of data to be classified.
  !!
  !! Assign label the class that data best fits.
  !/ -------------------------------------------------------------------------------------
  subroutine label_single( this, label, data )
    !/ -----------------------------------------------------------------------------------
    class(m_class), intent(inout) :: this
    integer,        intent(out)   :: label
    real(dp),       intent(in)    :: data(:)
    !/ -----------------------------------------------------------------------------------
    real(dp) :: metric, min_metric
    integer  :: c
    label      = 1

    min_metric = this%cluster(1)%metric( data )

    do c=2,this%nClass
       metric = this%cluster(c)%metric( data )
       if ( metric.lt.min_metric ) then
          min_metric = metric
          label      = c
       end if
    end do
  end subroutine label_single

  !/ =====================================================================================
  !> @brief Label.
  !! @param[in,out] this reference to this object.
  !! @param[out]    label array of that data fits.     
  !! @param[in]     data  sample of data to be classified.
  !!
  !! Assign a class to each label that corresponds to a row of data. Data is sample rows,
  !! and nvar columns.
  !/ -------------------------------------------------------------------------------------
  subroutine label_list( this, label, data )
    !/ -----------------------------------------------------------------------------------
    class(m_class), intent(inout) :: this
    integer,        intent(inout) :: label(:)
    real(dp),       intent(in)    :: data(:,:)
    !/ -----------------------------------------------------------------------------------
    integer :: s, ns, v, nv
    ns=size( data, dim=1 )
    nv=size( data, dim=2 )
    do s=1,ns
       do v=1,nv
          this%temp_row(v) = data(s,v)
       end do
       call label_single( this, label(s), this%temp_row )
    end do
  end subroutine label_list








  !/ =====================================================================================
  !> @brief Read Iris Data.
  !! @param[in,out] true_label allocatable array of labels.
  !! @param[in,out] data allocatable array for the data.
  !! @param[in]     file path to the labeled iris data.
  !! @param[out]    status 0 if no errors.
  !!
  !! Read the iris data and allocate the arrays.
  !! The first line is the number of samples and the number of variables.
  !! Each subsequent line has one integer label, followed by real variables.
  !/ -------------------------------------------------------------------------------------
  subroutine read_data( true_label, data, fspc, status )
    !/ -----------------------------------------------------------------------------------
    integer,  allocatable, intent(inout) :: true_label(:)
    real(dp), allocatable, intent(inout) :: data(:,:)
    character(len=*),      intent(in)    :: fspc
    integer,               intent(out)   :: status
    !/ -----------------------------------------------------------------------------------
    integer :: is, iv, ns, nv, un, ios 
    status = 0

    open( NEWUNIT=un, FILE=fspc, ACTION='READ', STATUS='OLD', IOSTAT=ios )
    if ( 0.ne.ios ) then
       write(*,110) fspc, ios
       status=ios
       goto 1100
    end if

    if ( allocated( true_label ) ) deallocate( true_label )
    if ( allocated( data ) )       deallocate( data )

    !/ -----------------------------------------------------------------------------------

    read(un,*,end=800,err=900) ns, nv

    allocate( true_label(ns) )
    allocate( data(ns,nv) )

    do is=1,ns
       read(un,*,end=800,err=900) true_label(is), ( data(is,iv), iv=1,nv )
    end do

    write(*,120) ns, nv

    goto 1000
    !/ -----------------------------------------------------------------------------------
110 format( 'Cannot open file:[',A,'] for reading, errno:',I0)
120 format( 'Read ',I0,' samples of ',I0,' variables' )

800 continue
    write(*,'(A)') '**E** read_data: Premature EOF encountered'
    status = 1
    goto 1000

900 continue
    write(*,'(A)') '**E** read_data: I/O error occurred'
    status = 2

1000 continue
    close(un)
1100 continue
  end subroutine read_data




end module test_iris








!/ =======================================================================================
program main
  !/ -------------------------------------------------------------------------------------
  use test_iris
  implicit none

  character(len=*), parameter :: CFG_FILE = 'iris.cfg'
  character(len=*), parameter :: DAT_FILE = 'iris.labeled'

  integer,  allocatable :: true_label(:)
  integer,  allocatable :: test_label(:)
  real(dp), allocatable :: data(:,:)
  integer               :: i, j, n, m, ios, match

  type(m_class)         :: classifier

  !/ -------------------------------------------------------------------------------------

  write (*,*)

  call create( classifier, CFG_FILE, status=ios )
  if ( 0.ne.ios ) then
     write(*,110) CFG_FILE
     goto 1000
  end if

  !/ -------------------------------------------------------------------------------------

  call read_data( true_label, data, DAT_FILE, status=ios )
  if ( 0.ne.ios ) then
     write(*,120) ios
     goto 1000
  end if

  n = size(true_label)
  allocate( test_label(n) )
  m = size( data, dim=2 )

  !/ -------------------------------------------------------------------------------------

  call classifier%label( test_label, data )

  write (*,*)

  match = 0
  do i=1,n
     if ( true_label(i).eq.test_label(i) ) then
        match = match + 1
     else
        write(*,130) test_label(i),i,true_label(i), ( data(i,j), j=1,m )
     end if
  end do

  write (*,*)
  write (*,140) match,n, ( real(match*100) / real(n) )
  write (*,*)

  !/ -------------------------------------------------------------------------------------

110 format( 'Failed to read classifier configuration: ',A)
120 format( 'Cannot read iris data, errno:',I0)
130 format( '  ',I0,' != sample(',I0,'): ',I0,' |',100(' ',F5.2) )
140 format( 'Matched ',I0,' out of ',I0,' (',F5.2,' percent)' )

  call delete( classifier )

1000 continue

  if ( allocated( data ) )       deallocate( data )
  if ( allocated( true_label ) ) deallocate( true_label )
  if ( allocated( test_label ) ) deallocate( test_label )

end program main

!/ =======================================================================================
!/ **                                 T E S T _ I R I S                                 **
!/ =========================================================================== END FILE ==
