module test_dll

  ! StatLab module to test dll build

  ! Crown Copyright (c) $shortdate: 2016-03-23 $ Defence Ordnance Safety Group

  ! $id: test_dll.f90 Sam Ellis 2016-03-23 1847cb3b7cee $

  ! .. preprocessing ..
  ! _short, _dll, _short_output and _no_terminal are assumed undefined unless defined in the compile command
  ! they are included here commented out to allow compilation in development
!#define _short
!#define _dll
!#define _cref
!#define _short_output
#define _no_terminal
#ifdef _no_lab_debug_filters
#define _entry_variables               = 0
#define _local_variables               = 0
#define _exit_variables                = 0
#else
#define _entry_variables               = 9
#define _local_variables               = 9
#define _exit_variables                = 9
#endif
#define _hide                          = 0
#define _kullback_mdis                 = 0

  ! .. used modules ..
!  use kinds
!  use ehlib
!  use module_ids, only : id_step,contingency_tables_id
!
!  use constants
!  use logarithms
!  use acm_toms_643

  ! .. implicit statement ..
  implicit none

  ! .. default attribute ..
  private

  ! .. public types ..
  public exception_t

  ! .. public module procedures ..
  public kullback_mdis
  public fisher_exact

  ! .. exception ids ..
  integer         (kind = ki4), parameter :: kullback_mdis_id =  1 * id_step
  integer         (kind = ki4), parameter :: fisher_exact_id  =  2 * id_step

  ! .. bound error module variables ..
#ifdef _short_output
  integer         (kind = ki4) :: log_status_m
#endif
  logical         (kind = kl4) :: bound_error_m
  type         (bound_error_t) :: bound_error_tm
  type           (exception_t) :: exception_m

  contains

  subroutine kullback_mdis(n_column,n_row,contingency_table,result,exception)

    ! routine to calculate Kullback's Minimum Discriminant Information Statistic

    ! .. scalar arguments ..
    integer         (kind = ki4), intent(in)    :: n_column,n_row
    real            (kind = kr8), intent(out)   :: result

    ! .. array arguments ..
    integer         (kind = ki4), intent(in),    dimension(n_row,n_column) :: contingency_table

    ! .. type arguments ..
    type           (exception_t), intent(inout) :: exception

    ! input
    ! n_column          - number of columns (1 < n_column)
    ! n_row             - number of rows (1 < n_row)
    ! contingency_table - array representing contingency table

    ! output
    ! result            - Kullback's Minimum Discriminant Information Statistic

    ! .. automatic arrays ..
    real            (kind = kr8), dimension(n_row,n_column) :: c    ! copy of contingency table
    real            (kind = kr8), dimension(1,n_column)     :: nr_  ! column sums
    real            (kind = kr8), dimension(n_row,1)        :: n_c  ! row sums

    ! ..local scalars ..
    real            (kind = kr8) :: n__                             ! overall sum

    ! ..
#if _hide == 0
#if _kullback_mdis > 0
#if _entry_variables > 0
    write (666,fmt = '(1x,a)')         '**** kullback_mdis: entry ****'
#if _entry_variables > 1
    write (666,fmt = '(1x,a)')         '**** kullback_mdis: scalars ****'
    write (666,fmt = '(1x,a,i12)')     'n_column',n_column
    write (666,fmt = '(1x,a,i12)')     'n_row   ',n_row
#if _entry_variables > 5
    write (666,fmt = '(1x,a)')         '**** kullback_mdis: scalar arrays ****'
    write (666,fmt = '(1x,a)')         'contingency_table(:,:)'
    write (666,fmt = '(1x,10i12)')      contingency_table(:,:)
#endif
#endif
#endif
#endif

    ! initialize
    result = nan

    ! check arguments
    call initialize_bound_error(bound_error_tm,contingency_tables_id,'contingency_tables',kullback_mdis_id,'kullback_mdis')
    bound_error_m = check_bound(bound_error_tm,1,1,'n_column',n_column,lt = 2)
    bound_error_m = check_bound(bound_error_tm,2,2,'n_row',n_row,lt = 2)
    bound_error_m = check_bound(bound_error_tm,3,3,'minval(contingency_table(:,:))',minval(contingency_table(:,:)),lt = 0)
    bound_error_m = check_bound(bound_error_tm,3,4,'maxval(contingency_table(:,:))',maxval(contingency_table(:,:)),lt = 1)

    A_if: if (bound_error_set(bound_error_tm)) then
      call raise_exception(exception,argument_error_p + error_p,5,bound_error_tm, &
        details_string = 'invalid arguments to kullback_mdis')

    else A_if

      ! copy contingency table into local array replacing zeros by tiny
      where (contingency_table > 0)
        c = real(contingency_table,kind = kr8)
      else where
        c = one_tiny
      end where

#if _kullback_mdis > 0
#if _local_variables > 0
      write (666,fmt = '(1x,a)')         '**** kullback_mdis: local ****'
#if _local_variables > 5
      write (666,fmt = '(1x,a)')         '**** kullback_mdis: scalar arrays ****'
      write (666,fmt = '(1x,a)')         'c(:,:)'
      write (666,fmt = '(1x,5es25.16)')   c(:,:)
#endif
#endif
#endif

      ! calculate marginal totals and total
      nr_(1,:) = sum(c,dim = 1)
      n_c(:,1) = sum(c,dim = 2)
      n__      = sum(c)

#if _kullback_mdis > 0
#if _local_variables > 0
      write (666,fmt = '(1x,a)')         '**** kullback_mdis: local ****'
#if _local_variables > 1
      write (666,fmt = '(1x,a)')         '**** kullback_mdis: scalars ****'
      write (666,fmt = '(1x,a,es25.16)') 'n__     ',n__
#if _local_variables > 5
      write (666,fmt = '(1x,a)')         '**** kullback_mdis: scalar arrays ****'
      write (666,fmt = '(1x,a)')         'nr_(:,1)'
      write (666,fmt = '(1x,5es25.16)')   nr_(:,1)
      write (666,fmt = '(1x,a)')         'n_c(1,:)'
      write (666,fmt = '(1x,5es25.16)')   n_c(1,:)
#endif
#endif
#endif
#endif

      ! calculate result
      result = 2.0_kr8 * sum(c * ln(c / (matmul(n_c,nr_) / n__)))

    end if A_if

#if _kullback_mdis > 0
#if _exit_variables > 0
    write (666,fmt = '(1x,a)')         '**** kullback_mdis: exit ****'
#if _exit_variables > 1
    write (666,fmt = '(1x,a)')         '**** kullback_mdis: scalars ****'
    write (666,fmt = '(1x,a,es25.16)') 'result  ',result
#endif
#endif
#endif
#endif
  end subroutine kullback_mdis

#ifdef _short
  subroutine z_kullback_mdis(n_column,n_row,contingency_table,result,error_code)

    ! short kullback_mdis

#ifdef _dll
    !DEC$ ATTRIBUTES DLLEXPORT :: z_kullback_mdis
    !DEC$ ATTRIBUTES ALIAS     : 'z_kullback_mdis' :: z_kullback_mdis
#ifdef _cref
    !DEC$ ATTRIBUTES C         :: z_kullback_mdis
#else
    !DEC$ ATTRIBUTES STDCALL   :: z_kullback_mdis
#endif
    !DEC$ ATTRIBUTES REFERENCE :: z_kullback_mdis
#endif

    ! .. scalar arguments ..
    integer         (kind = ki4), intent(in)    :: n_column,n_row
    integer         (kind = ki4), intent(out)   :: error_code
    real            (kind = kr8), intent(out)   :: result

    ! .. array arguments ..
    integer         (kind = ki4), intent(in),    dimension(n_column,n_row) :: contingency_table

    ! ..
#if _hide == 0

#ifdef _short_output
    ! turn on log output (gives both errors and warnings)
    log_status_m = restart_exception_log_unit(exception_m)
#endif
    write (666,'(a)') 'into z_kullback_mdis'

    ! call full routine
    call kullback_mdis(n_column,n_row,transpose(contingency_table),result,exception_m)

    ! clear exception
    call clear_exception(exception_m,error_code)

#endif
  end subroutine z_kullback_mdis

  subroutine z_scale(i4_in,r4_in,r8_in,i4_out,r4_out,r8_out,error_code)

    ! short double

#ifdef _dll
    !DEC$ ATTRIBUTES DLLEXPORT :: z_scale
    !DEC$ ATTRIBUTES ALIAS     : 'z_scale' :: z_scale
#ifdef _cref
    !DEC$ ATTRIBUTES C         :: z_scale
#else
    !DEC$ ATTRIBUTES STDCALL   :: z_scale
#endif
    !DEC$ ATTRIBUTES REFERENCE :: z_scale
#endif

    ! .. scalar arguments ..
    integer         (kind = ki4), intent(in)    :: i4_in
    real            (kind = kr4), intent(in)    :: r4_in
    real            (kind = kr8), intent(in)    :: r8_in
    integer         (kind = ki4), intent(out)   :: i4_out
    real            (kind = kr4), intent(out)   :: r4_out
    real            (kind = kr8), intent(out)   :: r8_out
    integer         (kind = ki4), intent(out)   :: error_code

    ! ..

    i4_out     = i4_in * 2
    r4_out     = r4_in * 3.0_kr4
    r8_out     = r8_in * 4.0_kr8
    error_code = 0

  end subroutine z_scale
#endif

#undef _private_code
#undef _hide
#undef _entry_variables
#undef _local_variables
#undef _exit_variables

end module test_dll
