halo_optimizer Program

    program halo_optimizer

    use halo_kinds_module,          only: wp
    use argv_module,                only: argv
    use simulated_annealing_module, only: simulated_annealing_type
    use json_module,                only: json_file
    use config_file_module,         only: config_file
    use halo_module

    implicit none

    type(simulated_annealing_type) :: sa
    real(wp), dimension(2) :: x !! [period, et_ref]
    real(wp) :: f
    integer :: iunit, istat

    integer,parameter :: n = 2 !! the 2 optimization variables are `et_ref` and `period` in the config file
    character(len=*),parameter :: default_base_config_file = './examples/example_sparse.json' !! default if not specified
    !character(len=*),parameter :: base_script_file = './run-sa.sh' ! the script to run for the function evaluation

    integer :: ns   !! number of cycles.
    integer :: nt  !! number of iterations before temperature reduction.
    integer :: neps  !! number of final function values used to decide upon termination.
    integer :: n_resets  !! number of times to run the main loop
    real(wp),dimension(:),allocatable :: lb !! lower bounds (size `n`)
    real(wp),dimension(:),allocatable :: ub !! upper bounds (size `n`)
    real(wp),dimension(:),allocatable :: vm !! size `n`
    real(wp) :: fopt, t, rt, period
    real(wp),dimension(n) :: xopt
    integer  :: nacc, nfcnev, ier
    logical :: found
    type(config_file) :: config
    integer :: epoch_mode, year, month, day, hour, minute
    real(wp) :: sec, et_ref
    character(len=:),allocatable :: base_config_file

    write(*,*) ''
    write(*,*) '----------------------'
    write(*,*) 'Initializing...'
    write(*,*) '----------------------'
    write(*,*) ''

    if (command_argument_count() == 0 ) then
        base_config_file = default_base_config_file
    else
        ! can pass this in as the first command line argument
        base_config_file = argv(1)
        write(*,*) 'Using base config file: '//base_config_file
    end if

    call config%open(base_config_file)
    call config%json%print()

    open(newunit = iunit, file = 'function_evals.csv', status='REPLACE', iostat=istat)
    write(iunit,'(A20,1X,2(A30,1X),A8)') 'CONFIG', 'PERIOD', 'ET', 'OBJFCN'

    ! initial guess:
    !x = [1.5872714606_wp, 757339200.0_wp] ! period, et [2024-01-01]
    call config%get('period', period) ! required
    call read_epoch(config, epoch_mode, year, month, day, hour, minute, sec, et_ref)
    x = [period, et_ref]

    ! defaults for optional arguments:
    t        = 5.0_wp
    lb       = [1.40_wp, x(2)-7*day2sec]  ! * vary epoch by +/- 7 days
    ub       = [1.60_wp, x(2)+7*day2sec]
    vm       = [0.2_wp, 7*day2sec]
    ! set some small values for testing...
    ns       = 3      ! number of cycles.
    nt       = 3      ! number of iterations before temperature reduction.
    neps     = 2      ! number of final function values used to decide upon termination.
    n_resets = 1      ! number of times to run the main loop
    rt       = 0.5_wp ! temperature reduction factor

    ! if these are specified in the config file then use those instead:
    call config%get('sa_t',       t,        found)
    call config%get('sa_lb',      lb,       found)
    call config%get('sa_ub',      ub,       found)
    call config%get('sa_vm',      vm,       found)
    call config%get('sa_ns',      ns,       found)
    call config%get('sa_nt',      nt,       found)
    call config%get('sa_neps',    neps,     found)
    call config%get('sa_n_resets',n_resets, found)
    call config%get('sa_rt',      rt,       found)

    write(*,*) ''
    write(*,*) '----------------------'
    write(*,*) 'Running SA optimizer..'
    write(*,*) '----------------------'
    write(*,*) ''

    call sa%initialize(fcn,n,lb,ub,&
                       n_resets = n_resets, &
                       ns       = ns, &
                       nt       = nt, &
                       neps     = neps, &
                       optimal_f_specified = .true.,&
                       optimal_f     = 0.0_wp,& ! 0 is the desired solution (no eclipses) if it's found, we can stop.
                       optimal_f_tol = epsilon(1.0_wp)) ! since this is really an int, this doesn't matter much
    call sa%optimize(x, rt, t, vm, xopt, fopt, nacc, nfcnev, ier)

    write(*,*) ''
    write(*,*) '----------------------'
    write(*,*) 'Done with SA'
    write(*,*) '----------------------'
    write(*,*) ''

    write(*,*) ''
    write(*,*) '=============='
    write(*,*) '   SOLUTION   '
    write(*,*) '=============='
    write(*,*) ''
    write(*,*) 'xopt = ', xopt
    write(*,*) 'fopt = ', fopt
    write(*,*) 'ier  = ', ier
    write(*,*) ''

    !TODO: maybe make a new config for the solution ?

    close(iunit, iostat=istat) ! function eval file
    call config%close()

    contains

        subroutine fcn(me, x, f, istat)

            !! function for the SA algorithm.

            class(simulated_annealing_type),intent(inout) :: me
            real(wp), dimension(:), intent(in) :: x !! [period, et_ref]
            real(wp), intent(out) :: f
            integer,intent(out) :: istat

            integer,save :: istep = 0 !! number of function evaluations

            character(len=5) :: istep_str
            character(len=:),allocatable :: script !! script to run
            character(len=:),allocatable :: step_config !! config file for this step
            character(len=:),allocatable :: eclipse_file !! eclipse file
            type(json_file) :: json !! for loading the eclipse json file
            real(wp),dimension(:),allocatable :: et_vec !! ephemeris time of eclipse violations
            type(my_solver_type) :: solver !! just so we can read the base config and generate the case name

            istep = istep + 1
            write(istep_str,'(I5)') istep

            ! create the script
            ! read the original script and modify with the inputs, then save the new one
            step_config = './STEP='//trim(adjustl(istep_str))//'.json'
            call config%json%update('period', x(1), found)
            call config%json%update('et_ref', x(2), found)
            ! make sure the output settings are right:
            call config%json%update('initial_guess_from_file', '', found)
            call config%json%update('generate_plots',                    .false., found)
            call config%json%update('generate_trajectory_files',         .false., found)
            call config%json%update('generate_guess_and_solution_files', .false., found)
            call config%json%update('generate_kernel',                   .false., found)
            call config%json%update('generate_defect_file',              .false., found)
            call config%json%update('generate_json_trajectory_file',     .false., found)
            call config%json%update('run_pyvista_script',                .false., found)
            call config%json%update('generate_eclipse_files', .true., found)  ! only need these
            call config%json%update('eclipse_filetype', 2, found)             !
            call config%json%remove('year')
            call config%json%remove('month')
            call config%json%remove('day')
            call config%json%remove('hour')
            call config%json%remove('minute')
            call config%json%remove('sec')
            call config%json%print(step_config)

            ! have to read the config so that get_case_name will return the right string:
            write(*,*) 'reading base config file: '//base_config_file
            call solver%read_config_file(step_config)

            ! the env var is used by the script to indicate the config file to use:
            !call execute_command_line('export HALO_CONFIG_FILE='//step_config//'; '//base_script_file)
            ! new way: call the solver as a function:
            call halo_solver_main(step_config,debug=.false.)

            ! read the eclipse function output file
            eclipse_file = './eclipse_'//solver%mission%get_case_name()//'.json'
            call json%load(eclipse_file)
            call json%get('et', et_vec, found)
            if (.not. found) error stop 'error reading eclipse file. et not found.'
            call json%destroy()

            ! f is just the count the 1-hr epochs where the phi is < 0
            !TODO:: could use unique() to remove duplicate et's. just count the unique ets (can be zero)
            f = real(size(et_vec), wp)

            istat = 0 ! no problems

            write(*,*) '======================'
            write(*,*) ' Latest f: ', f
            write(*,*) '======================'
            write(iunit,'(A20,1X,2(E30.16,1X),I8)') step_config, x, int(f)

        end subroutine fcn

    end program halo_optimizer