test2 Program

    program test2

    use iso_fortran_env
    use numerical_differentiation_module
    use numdiff_kinds_module, only: wp
    use pyplot_module

    implicit none

    integer,parameter :: n = 1 !! number of variables
    integer,parameter :: m = 1 !! number of functions
    real(wp),dimension(n),parameter :: x     = 1.0_wp !! point at which to compute the derivative
    real(wp),dimension(n),parameter :: xlow  = -100000.0_wp !! bounds not really needed
    real(wp),dimension(n),parameter :: xhigh =  100000.0_wp !!
    integer,parameter :: perturb_mode  = 1 !! absolute step
    integer,parameter :: cache_size    = 0 !! `0` indicates not to use cache
    integer,parameter :: sparsity_mode = 1 !! assume dense
    integer,dimension(*),parameter  :: methods = [1,3,10,21,36,500,600,700,800]
                                                !! array of method IDs:
                                                !! [1,3,10,21,36,500,600,700,800]
                                                !! forward + all the central diff ones
    integer,parameter :: exp_star = -ceiling(log10(epsilon(1.0_wp))) !! exponent to start with
    integer,parameter :: exp_stop = -2 !! exponent to end with
    integer,parameter :: exp_step = -1 !! ecponent step
    integer,parameter :: exp_scale = 5 !! number of substeps from one to the next

    type(numdiff_type)                  :: my_prob          !! main class to compute the derivatives
    integer                             :: i                !! method counter
    integer                             :: j                !! counter
    real(wp),dimension(:),allocatable   :: jac              !! jacobian
    integer                             :: func_evals       !! function evaluation counter
    character(len=:),allocatable        :: error_msg        !! error message string
    real(wp),dimension(n)               :: dpert            !! perturbation step size
    integer                             :: ipert            !! perturbation step size counter
    real(wp)                            :: error            !! diff from true derivative
    integer                             :: num_dperts       !! number of dperts to test
    integer                             :: num_methods      !! number of methods to test
    real(wp),dimension(:),allocatable   :: results_dpert    !! results array - dpert
    real(wp),dimension(:),allocatable   :: results_errors   !! results array - errors
    type(pyplot)                        :: plt              !! for plotting the results
    character(len=:),allocatable        :: formula          !! finite diff formula for the plot legend
    character(len=:),allocatable        :: name             !! finite diff name for the plot legend
    integer                             :: idx              !! index in results arrays
    character(len=:),allocatable        :: real_kind_str    !! real kind for the plot title
    real(wp),dimension(3)               :: color            !! line color array
    real(wp),dimension(2)               :: ylim             !! plot y limit array

    ! size the arrays:
    num_methods = size(methods)
    num_dperts = size([(i, i = exp_star*exp_scale, exp_stop*exp_scale, exp_step)])
    allocate(results_errors(num_dperts)); results_errors = -huge(1.0_wp)

    ! for plot title:
    select case(wp)
        case(REAL32);  real_kind_str = '[Single Precision]'
        case(REAL64);  real_kind_str = '[Double Precision]'
        case(REAL128); real_kind_str = '[Quad Precision]'
        case default; error stop 'Invalid real kind'
    end select

    ! initialize the plot:
    call plt%initialize(grid            = .true.,&
                        figsize         = [20,10],&
                        axes_labelsize  = 30, &
                        xtick_labelsize = 30, &
                        ytick_labelsize = 30, &
                        font_size       = 30, &
                        xlabel          = 'Finite Difference Perturbation Step Size $h$',&
                        ylabel          = 'Finite Difference Derivative Error',&
                        title           = 'Derivative of $x + \sin(x)$ at $x=1$ '//real_kind_str,&
                        legend          = .true., &
                        legend_fontsize = 10,&
                        usetex          = .true.)

    ! try different finite diff methods
    do j=1,size(methods)

        idx = 0
        i = methods(j)  ! method id
        call get_finite_diff_formula(i,formula,name)

        ! cycle through perturbation step sizes:
        do ipert = exp_star*exp_scale, exp_stop*exp_scale, exp_step

            idx = idx + 1 ! index for arrays
            dpert = 10.0_wp ** (-ipert/real(exp_scale,wp)) ! compute perturbation step size
            if (j==1) then
                if (.not. allocated(results_dpert)) allocate(results_dpert(0))
                results_dpert = [results_dpert, dpert(1)] ! save dpert
            end if

            func_evals = 0
            call my_prob%destroy()
            call my_prob%initialize(n,m,xlow,xhigh,perturb_mode,dpert,&
                                    problem_func=my_func,&
                                    sparsity_mode=sparsity_mode,&
                                    jacobian_method=i,&
                                    partition_sparsity_pattern=.false.,&
                                    cache_size=cache_size)
            if (my_prob%failed()) then
                call my_prob%get_error_status(error_msg=error_msg)
                error stop error_msg
            end if

            call my_prob%compute_jacobian(x,jac)
            if (my_prob%failed()) then
                call my_prob%get_error_status(error_msg=error_msg)
                error stop error_msg
            end if
            error = abs(deriv(x(1)) - jac(1))
            !write(output_unit,'(I5,1X,*(F27.16))') i , dpert, error
            if (error<=epsilon(1.0_wp)) error = 0.0_wp
            results_errors(idx) = error ! save result

        end do

        ! line color for the plot:
        select case (j)
        case(1); color = [1.0_wp, 0.0_wp, 0.0_wp] ! red
        case(2); color = [0.0_wp, 1.0_wp, 0.0_wp] ! green
        case default
            ! blue gradient for the others
            color = [real(j-2,wp)/num_methods, &
                     real(j-2,wp)/num_methods, &
                     0.9_wp]
        end select
        ylim = [10.0_wp**(ceiling(log10(epsilon(1.0_wp)))), 1.0_wp]

        ! plot for this method:
        call plt%add_plot(results_dpert,results_errors,&
                          xscale='log', yscale='log',&
                          label = formula,linestyle='.-',markersize=5,linewidth=2, &
                          color = color,&
                          xlim = ylim,&
                          ylim = ylim)

    end do

    ! save plot:
    call plt%savefig('results '//real_kind_str//'.pdf', pyfile='results.py')

contains

    function func(x)
        !! Problem function
        real(wp), intent(in) :: x
        real(wp) :: func
        func = x + sin(x)
    end function func
    function deriv(x)
        !! Problem function true derivative
        real(wp), intent(in) :: x
        real(wp) :: deriv
        deriv = 1.0_wp + cos(x)
    end function deriv

    subroutine my_func(me,x,f,funcs_to_compute)
        !! Problem function interface for numdiff
        class(numdiff_type),intent(inout) :: me
        real(wp),dimension(:),intent(in)  :: x
        real(wp),dimension(:),intent(out) :: f
        integer,dimension(:),intent(in)   :: funcs_to_compute
        if (any(funcs_to_compute==1)) f(1) = func(x(1))
        func_evals = func_evals + 1
    end subroutine my_func

!*******************************************************************************
    end program test2