!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! fit_model.f90
!
! You are welcome to use this code for your own research but please 
! properly cite and attibute:
!
! Author: 
! Cameron Bracken, cameron.bracken@gmail.com
!
! Citation: 
! C. Bracken, K. Holman, B. Rajagopalan, H. Moradkhani, A Bayesian 
! hierarchical approach to multivariate nonstationary hydrologic frequency 
! analysis, Water Resources Research, 2017, Under Review. 
!
! This is the main executable program for running the sampler. 
! It has a number of options.
!
! usage: cmdline [OPTIONS]
! 
! options:
! 
!   -h, --help                print this message
!   -s, --sample-file file    name of the file to save samples [default=samples]
!   -l, --lp-file file        name of the file to save samples [default=none]
!   -i, -iterations n         number of iterations to run [default=20000]
!   -r, -refresh n            print after each refresh iterations [default=1]
!   --no-sample-output        do not write out any samples (useful for model selection).
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program fit
    use model
    use mfu_sampler

    implicit none

    double precision, allocatable:: x(:)
    integer:: i, iterations, warmup, refresh, iflen, ier, narg
    double precision:: t1, t2, total1, total2, lp
    character(len=1000):: sample_file, lp_file, iterations_c, arg
    logical:: write_lp = .false., skip_next_arg = .false., no_sample_output = .false.

    ! defaults
    iterations = 20000
    warmup = 10000
    sample_file = 'samples'
    lp_file = ''
    write_lp = .false.
    refresh = 1

    narg = command_argument_count() 

    if(narg>0)then
      !loop across options
      do i=1,narg
        if(skip_next_arg) then
          skip_next_arg = .false.
          cycle
        end if
        call get_command_argument(i, arg)
          select case(adjustl(trim(arg)))

            case("--help","-h")
              call get_command_argument(i+1, arg)
              skip_next_arg = .false.
              call print_help()

            case("--sample-file","-s")
              call get_command_argument(i+1, arg)
              skip_next_arg = .true.
              sample_file = arg 

            case("--lp-file","-l")
              call get_command_argument(i+1, arg)
              skip_next_arg = .true.
              lp_file = arg 
              write_lp = .true.

            case("--iterations","-i")
              call get_command_argument(i+1, arg)
              skip_next_arg = .true.
              read(arg,'(I10)') iterations

            case("--warmup","-w")
              call get_command_argument(i+1, arg)
              skip_next_arg = .true.
              read(arg,'(I10)') warmup

            case("--refresh","-r")
              call get_command_argument(i+1, arg)
              skip_next_arg = .true.
              read(arg,'(I10)') refresh

            case("--no-sample-output")
              call get_command_argument(i+1, arg)
              skip_next_arg = .false.
              no_sample_output = .true.

            case default
              write(*,*)"Option ",adjustl(trim(arg))," unknown"
              call print_help()

         end select
      end do
    end if

    write(*,*)
    write(*,*)"Sample file is: ", trim(sample_file)
    if(write_lp) write(*,*)"log posterior file is: ", trim(lp_file)
    write(*,*)"Iterations: ", iterations

    write(*,*)"Refresh: ", refresh
    write(*,*)

    open(88,file=trim(sample_file))
    if(write_lp) open(99,file=trim(lp_file))

    !
    call init_random_seed()

    ! read data into global variables (yeah I know)
    call read_data()
    call mfu_init()

    !write(*,*)'Target Ratio:', target_ratio

    allocate(x(n_pars))
    
    ! defined in read_data subroutine 
    x = x0
    lp = lp_fun(x)
    write(*,*)'Initial lp: ',lp

    call cpu_time(total1)

    write(*,*)'Sampling started.'
    do i = 1, iterations

        if(i .eq. 1) call cpu_time(t1)

        call mfu_sample(x, i, warmup)
        lp = lp_fun(x)

        if(.not. no_sample_output)write(88,*) x
        if(write_lp) write(99,*) lp

        if(modulo(i,refresh) .eq. 0)then
            call cpu_time(t2)
            write(*,'(a10,i6,a4,i6,a1,f8.3,a6,f15.2)') 'Iteration ', i,' of ', iterations,',', t2-t1, 's, lp:', lp
            call cpu_time(t1)
        end if

    end do

    call cpu_time(total2)

    write(*,*) 'Computation time: '
    write(*,'(f15.1,a8)') total2 - total1, ' seconds'

    close(88)
    close(99) 

    call exit(0)

end program 


subroutine print_help()
    write(*,*) 'usage: cmdline [OPTIONS]'
    write(*,*) ''
    write(*,*) 'options:'
    write(*,*) ''
    write(*,*) '  -h, --help                print this message'
    write(*,*) '  -s, --sample-file file    name of the file to save samples [default=samples]'
    write(*,*) '  -l, --lp-file file        name of the file to save samples [default=none]'
    write(*,*) '  -i, -iterations n         number of iterations to run [default=20000]'
    write(*,*) '  -r, -refresh n            print after each refresh iterations [default=1]'
    write(*,*) '  --no-sample-output        do not write out any samples (useful for model selection).'
    call exit(0)
end subroutine print_help