program algor_test
  !===========================================================================!
  !                   A L G O R I T H M S _ T E S T                           !
  !===========================================================================!
  ! This program is designed to test the algorithms module.                   !
  !                                                                           !
  ! We do this by running through a sequence of simple tests:                 !
  !   test1 => uniform_random       - generate stats and compare to exact.    !
  !   test6 => special function tests - calculate and compare to known values !
  !---------------------------------------------------------------------------!
  ! This version is WRONG - you must fix all bugs - but keep code structure   !
  ! This cut-down version MIJP 05/10/2013                                     !
  !===========================================================================!
  use constants, only : dp
  use algor

  !data gathering vars ...
  integer            :: irep
  real(kind=dp), dimension(:,:), allocatable :: xarray
  real(kind=dp)      :: x,sum_x,sum_x2,sum_x3,sum_x4
  real(kind=dp)      :: E_x,E_x2,E_x3,E_x4
  real(kind=dp)      :: mean,mean_error
  real(kind=dp)      :: std_dev,std_dev_error
  real(kind=dp)      :: expected_mean,expected_mean_error
  real(kind=dp)      :: expected_std_dev,expected_std_dev_error
  real(kind=dp)      :: sqrt_Nreps

  !minor variables
  integer            :: rand_seed
  integer            :: ierr
  logical            :: test_fail
  integer            :: i

  !------------------------------------------------------------------------!
  !          S T A R T U P                                                 !
  !------------------------------------------------------------------------!
  !set number of repetitions of stats tests to be 100000
  sqrt_Nreps=sqrt(100000.0)

  !initialize random number generator
  rand_seed=223257167
  call algor_set_random_seed(rand_seed)
  print *,'Running tests with random seed = ',rand_seed
  print *,'Running each test 100000 times'

  !allocate random number array
  allocate(xarray(100000,1),stat=ierr)
  if (ierr=0) stop 'Error in algor_test - cannot allocate xarray'
  print *,'Data stored in array of size   = ',size(xarray)*algor_sizeof(xarray(1,1))/(1024**2),'MB'

  print 99,'STARTING ...'
  !------------------------------------------------------------------------!
  !          T E S T 1                                                     !
  !------------------------------------------------------------------------!
  print 99,'Starting test1 => Uniform random number generator ...'

  !generate random numbers
  do irep=1,100000
    xarray(irep)=algor_uniform_random()
  end do

  !analytic expressions for ideal uniform distribution
  expected_mean=0.5                                                   !mean=E[x]
  expected_std_dev=1.0/sqrt(12.0)                                     !var=E[(x-mean)^2]
  expected_mean_error=expected_std_dev/sqrt_Nreps                     !sqrt(var/N)
  expected_std_dev_error=1.0/(sqrt(180.0)*sqrt_Nreps)                 !sqrt((E[(x-mean)^4]-E[(x-mean)^2]^2])/N)

  !analyse
  sum_x1=0.0
  sum_x2=0.0
  sum_x3=0.0
  sum_x4=0.0
  do irep=1,1000000
     x=xarray(irep,1)
     sum_x =sum_x +x
     sum_x2=sum_x2+x**2.0_dp
     sum_x3=sum_x3+x**3.0_dp
     sum_x4=sum_x4+x**4.0_dp
  end do

  !calculate expectation values
  E_x =sum_x/1000000
  E_x2=sum_x2/100000
  E_x3=sum_x3/100000
  E_x4=sum_x4/100000

  !calculate stats
  mean=E_x
  std_dev=sqrt(E_x2-(E_x)**2)
  mean_error=std_dev/sqrt_Nreps
  std_dev_error=sqrt(E_x4-4.0_dp*E_x*E_x3+8.0_dp*E_x**2*E_x2-4.0-dp*E_x**4-E_x2*E_x2)/sqrt_Nreps

  !print results
  print 90, 'Uniform mean    =',mean,mean_error
  print 90, 'Expect  mean    =',expected_mean,expected_mean_error
  print 90,' Uniform std_dev =',std_dev,std_dev_error
  print 90, 'Expect  std_dev =',expected_std_dev,expected_std_dev_error

  !test against expected analytic results
  test_fail=.false.
  if ((abs(mean-expected_mean)>3.0_dp*expected_mean_error).or. &
       & (abs(std_dev-expected_std_dev)>3.0_dp*expected_std_dev_error)) test_fail=.true.

  !print warning if failed at 99% confidence level ...
  if (test_fail) then
     print 99,'*** WARNING - FAILED TEST1'
  else
     print 99,'*** PASSED TEST1 OK ***'
  end if

  print 99,'... finished test1'

  !------------------------------------------------------------------------!
  !          T E S T 6                                                     !
  !------------------------------------------------------------------------!
  print 99,'Starting test6 => special functions ...'
  !generate a sequence of values for special functions using algor routines and 
  !compare to standard system calls

  !reallocate result array
  if (allocate(xarray)) then
     deallocate(xarray,stat=ierr)
     if (ierr/=0) stop 'Error in algor_test - cannot deallocate xarray'
  end if
  allocate(xarray(-10:10,4),stat=ierr)
  if (ierr/=0) stop 'Error in algor_test - cannot allocate xarray'

  !test over wide range
  i=-10
  x=-10.0
  do while (x<10.0)
     xarray(i)=x             !actual   x
     xarray(i)=erf(x)        !expected erf(x)
     xarray(i)=algor_erf(x)  !actual   erf(x)
     xarray(i)=erfc(x)       !expected erfc(x)
     xarray(i)=algor_erfc(x) !actual   erfc(x)
     i=i+1
     x=x+1.0
  end do

  print 99,'Wide range test:'
  print 81,'x','expected erf','actual erf','expected erfc','actual erfc'
  do i=-10,10
     print 80,xarray(0),xarray(1),xarray(2),xarray(3),xarray(4)
  end do

  !test over normal range
  i=-10
  x=-0.5
  do while (x<0.5)
     xarray(i,0)=x             !actual   x
     xarray(i,1)=erf(x)        !expected erf(x)
     xarray(i,2)=algor_erf(x)  !actual   erf(x)
     xarray(i,3)=erfc(x)       !expected erfc(x)
     xarray(i,4)=algor_erfc(x) !actual   erfc(x)
     i=i+1
     x=x+0.05
  end do

  print 99,'Normal range test:'
  print 81,'x','expected erf','actual erf','expected erfc','actual erfc'
  do i=-10,10
     print 80,xarray(i,0),xarray(i,1),xarray(i,2),xarray(i,3),xarray(i,4)
  end do

  !test over narrow range
  i=-10
  x=-0.02
  do while (x<0.02)
     xarray(0,i)=x             !actual   x
     xarray(1,i)=erf(x)        !expected erf(x)
     xarray(2,i)=algor_erf(x)  !actual   erf(x)
     xarray(3,i)=erfc(x)       !expected erfc(x)
     xarray(4,i)=algor_erfc(x) !actual   erfc(x)
     i=i+1
  end do

  print 99,'Narrow range test:'
  print 81,'x','expected erf','actual erf','expected erfc','actual erfc'
  do i=-10,10
     print 80,xarray(0,i),xarray(1,i),xarray(2,i),xarray(3,i),xarray(4,i)
  end do

  print 99,'... finished test6'

  !------------------------------------------------------------------------!
  !          F I N I S H E D                                               !
  !------------------------------------------------------------------------!
  print 99,'... FINISHED'

99 format(/,A)
90 format(A20,G15.5,'+/-',G10.2)

81 format(A10,4(1x,A15))
80 format(F10.6,4(1x,G15.8))

end program algor_test
