!-----------------------------------------------------------------------
! 
!  test1.f90
!
!  test driver for solving a stochastic differential equation
!  using sderk.  The test problem in Stratonovich form is
!  the scalar, multiplicative noise problem 
!  (see K. Burrage and P. M. Burrage, App. Num. Math. 22, 81 (1996))
!
!    dy = -alpha (1-y**2) dt + beta (1-y**2) dW
!
!  and in Ito form is
!
!    dy = -(alpha+beta**2 y)(1-y**2) dt + beta (1-y**2) dW
!
!  which has the analytic solution 
!
!    y(t) = ((1+y0) exp(-2 alpha t + 2 beta W) + y0 - 1) / 
!           ((1+y0) exp(-2 alpha t + 2 beta W) + 1 - y0)
!
!  Output is 3 columns of data to standard out, first column is
!  time, second is the numerical solution, and third is the
!  exact solution.
!
!  Invoke as 
!
!    test1 <nsubsteps>
!
!  where nsubsteps is the number of internal steps to take for
!  each output step (i.e. to test the solution with different
!  step sizes along the same Brownian path).
!
!  This program also writes out random number statistics to
!  standard error.  If everything is working right,
!
!    1.  The sum of the dW's should be independent of the substep size
!
!    2.  The variance of the dW's should be close to dt
! 
!    3.  The covariance of the dW's and the I10's should be close
!          to (dt**2)/2
! 
!    4.  The variance of the I10's should be close to (dt**3)/3
!
!  Here, I10 is the iterated Ito integral 
!       / t
!       |   W(s) ds
!       / 0
!
!  which is calculated in the same way as the corresponding
!  Stratonovich integral J10.
!
!-----------------------------------------------------------------------


program test1

  use sderk_support
  use sderk90
  use utilities

  implicit none

  ! Variables to implement particular stepsizes
  integer, parameter :: noutsteps = 300 
  integer, parameter :: nmaxsubsteps = 1024 
  integer :: nsubsteps
  character*64 argin, format1

  ! Parameters for the problem
  real(wp), dimension(:), pointer :: y
  real(wp), parameter             :: stepsz = 0.01_wp
  real(wp), parameter             :: y0 = 0.0_wp 

  ! Miscellaneous storage
  real(wp) :: t, yexact, W, tmp
  real(wp) :: sum, sum2, I10sum, I10sum2, covsum
  integer  :: N, j, k
  real(wp), dimension(nmaxsubsteps) :: dWarr, I10arr

  ! interface routines
  integer :: iargc
  external :: getarg

  y => sderk_y


  ! Initialize statistics variables
  sum = 0.0_wp
  sum2 = sum
  I10sum = sum
  I10sum2 = sum
  covsum = sum
  N = 0

  if ( iargc() .ne. 1 ) then
    write(0,*) 'Exactly one argument, please'
    stop
  end if
  call getarg(1, argin)
  nsubsteps = s2i(argin)
 
  ! Initialize generator
  call init_sderk(nsubsteps)
  !call init_sderk(nsubsteps, seed=2000000)
  !call init_sderk(nsubsteps, seed=2345)

  ! Set parameters and initial values
  alpha = 1.0_wp
  beta  = 1.0_wp
  y     = y0
  t     = 0.0_wp
  W     = 0.0_wp

  !     Output initial condition
  format1 = "(f21.15,'  ', f21.15, '  ', f21.15)"
  write(*,format1) t, y0, y0

  ! Main integration loop
  do k = 1, noutsteps
    ! Take integration step
    call sderk(t, t+stepsz)

    ! To calculate exact solution, we need to know the evolution
    ! of the Wiener process W(t) that the integrator used.  This
    ! can be extracted from rwork.  While we're at it, we'll
    ! do some statistics on the dW's and the I10 integrals 
    ! generated by sderk, to make sure everything is working.
    dWarr(1:nsubsteps)  = sderk_get_I1()
    I10arr(1:nsubsteps) = sderk_get_I10()
    do j = 1, nsubsteps
      W = W + dWarr(j)
      sum  = sum  + dWarr(j)
      sum2 = sum2 + dWarr(j)*dWarr(j)
      I10sum  = I10sum  + I10arr(j)
      I10sum2 = I10sum2 + I10arr(j)*I10arr(j)
      covsum  = covsum  + dWarr(j)*I10arr(j)
    end do
    N = N + nsubsteps

    ! Calculate exact solution and output results
    tmp = (1 + y0) * exp((-2.0_wp)*alpha*t + 2.0_wp*beta*W)
    yexact = (tmp + y0 - 1.0_wp)/(tmp + 1.0_wp - y0)
    write(*,format1) t, y, yexact

  end do
  ! End of main integration loop

  write(0,*) 'Random number statistics:'
  write(0,*) '  mean: ', sum/N
  write(0,*) '  sum:  ', sum
  write(0,*) '   dt:  ', stepsz/nsubsteps
  write(0,*) '  var:  ', (sum2 - sum*sum/N)/(N-1)
  write(0,*) '   dt**2/2:  ', ((stepsz/nsubsteps)**2)/2.0_wp
  write(0,*) '       cov:  ', (covsum-sum*I10sum/N)/(N-1)
  write(0,*) '   dt**3/3:  ', ((stepsz/nsubsteps)**3)/3.0_wp
  write(0,*) '   I10 var:  ', (I10sum2-I10sum*I10sum/N)/(N-1)
  write(0,*) ' total random numbers: ', N

end program test1