Index: trunk/src/vegas/vegas.nw =================================================================== --- trunk/src/vegas/vegas.nw (revision 8344) +++ trunk/src/vegas/vegas.nw (revision 8345) @@ -1,4886 +1,4909 @@ %% -*- ess-noweb-default-code-mode: f90-mode; noweb-default-code-mode: f90-mode; -*- % WHIZARD code as NOWEB source: VEGAS algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{VEGAS Integration} \label{cha:vegas-integration} @ The backbone integrator of WHIZARD is a object-oriented implemetation of the VEGAS algorithm. <<[[vegas.f90]]>>= <> module vegas <> <> <> <> <> <> <> contains <> end module vegas @ %def vegas <>= use diagnostics use io_units use format_utils, only: write_indent use format_defs, only: FMT_17 use rng_base, only: rng_t use rng_stream, only: rng_stream_t @ @ MPI Module. <>= use mpi_f08 !NODEP! @ \section{Integration modes} \label{sec:integration-modes} @ VEGAS operates in three different modes: [[vegas_mode_importance_only]], [[vegas_mode_importance]] or [[vegas_mode_stratified]]. The default mode is [[vegas_mode_importance]], where the algorithm decides whether if it is possible to use importance sampling or stratified sampling. In low dimensions VEGAS uses strict stratified sampling. <>= integer, parameter, public :: VEGAS_MODE_IMPORTANCE = 0, & & VEGAS_MODE_STRATIFIED = 1, VEGAS_MODE_IMPORTANCE_ONLY = 2 @ %def vegas_mode_importance vegas_mode_stratified vegas_mode_importance_only @ \section{Type: vegas\_func\_t} \label{sec:type:vegas_func_t} We define a abstract [[func]] type which only gives an interface to an [[evaluate]] procedure. The inside of implementation and also the optimization of are not a concern of the [[vegas]] implementation. <>= public :: vegas_func_t <>= type, abstract :: vegas_func_t ! contains procedure(vegas_func_evaluate), deferred, pass, public :: evaluate end type vegas_func_t @ %def vegas_func_t @ The only procedure called in [[vegas]] is [[vegas_func_evaluate]]. It takes a real value [[x]] and returns a real value [[f]]. <>= abstract interface real(default) function vegas_func_evaluate (self, x) result (f) import :: default, vegas_func_t class(vegas_func_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x end function vegas_func_evaluate end interface @ %def vegas_func_evaluate @ \section{Type: vegas\_config\_t} \label{sec:type:vegas_config_t} We store the complete configuration in a transparent container. The [[vegas_config_t]] object inside VEGAS must not be directly accesible. We provide a get method which returns a copy of the [[vegas_config_t]] object. Apart from the options which can be set by the constructor of [[vegas_t]] object, we store the run-time configuration [[n_calls]], [[calls_per_box]], [[n_bins]] and [[n_boxes]]. Those are calculated and set accordingly by VEGAS. <>= public :: vegas_config_t <>= type :: vegas_config_t integer :: n_dim = 0 real(default) :: alpha = 1.5 integer :: n_bins_max = 50 integer :: iterations = 5 integer :: mode = VEGAS_MODE_STRATIFIED integer :: calls_per_box = 0 integer :: n_calls = 0 integer :: n_calls_min = 20 integer :: n_boxes = 1 integer :: n_bins = 1 contains <> end type vegas_config_t @ %def vegas_config_t, n_calls, calls_per_box, n_bins, n_boxes @ Write out the configuration of the grid. <>= procedure, public :: write => vegas_config_write <>= subroutine vegas_config_write (self, unit, indent) class(vegas_config_t), intent(in) :: self integer, intent(in), optional :: unit integer, intent(in), optional :: indent integer :: u, ind u = given_output_unit (unit) ind = 0; if (present (indent)) ind = indent call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Number of dimensions = ", self%n_dim call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Adaption power (alpha) = ", self%alpha call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Max. number of bins (per dim.) = ", self%n_bins_max call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Number of iterations = ", self%iterations call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Mode (stratified or importance) = ", self%mode call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Calls per box = ", self%calls_per_box call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Number of calls = ", self%n_calls call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Min. number of calls = ", self%n_calls_min call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Number of bins = ", self%n_bins call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Number of boxes = ", self%n_boxes end subroutine vegas_config_write @ %def vegas_config_write @ \section{Type: vegas\_grid\_t} \label{sec:type:-vegas_g} We provide a simple and transparent grid container. The container can then later be used, to export the actual grid. <>= public :: vegas_grid_t <>= type :: vegas_grid_t integer :: n_dim = 1 integer :: n_bins = 1 real(default), dimension(:), allocatable :: x_lower real(default), dimension(:), allocatable :: x_upper real(default), dimension(:), allocatable :: delta_x real(default), dimension(:,:), allocatable :: xi contains <> end type vegas_grid_t @ %def vegas_grid_t @ Initialise grid. <>= interface vegas_grid_t module procedure vegas_grid_init end interface vegas_grid_t <>= type(vegas_grid_t) function vegas_grid_init (n_dim, n_bins_max) result (self) integer, intent(in) :: n_dim integer, intent(in) :: n_bins_max self%n_dim = n_dim self%n_bins = 1 allocate (self%x_upper(n_dim), source=1.0_default) allocate (self%x_lower(n_dim), source=0.0_default) allocate (self%delta_x(n_dim), source=1.0_default) allocate (self%xi((n_bins_max + 1), n_dim), source=0.0_default) end function vegas_grid_init @ %def vegas_grid_init @ Output. <>= procedure, public :: write => vegas_grid_write <>= - subroutine vegas_grid_write (self, unit) + subroutine vegas_grid_write (self, unit, pacify) class(vegas_grid_t), intent(in) :: self integer, intent(in), optional :: unit + logical, intent(in), optional :: pacify + logical :: pac integer :: u, i, j + pac = .false.; if (present (pacify)) pac = pacify u = given_output_unit (unit) write (u, descr_fmt) "begin vegas_grid_t" write (u, integer_fmt) "n_dim = ", self%n_dim write (u, integer_fmt) "n_bins = ", self%n_bins write (u, descr_fmt) "begin x_lower" do j = 1, self%n_dim - write (u, double_array_fmt) j, self%x_lower(j) + if (pac) then + write (u, double_array_pac_fmt) j, self%x_lower(j) + else + write (u, double_array_fmt) j, self%x_lower(j) + end if end do write (u, descr_fmt) "end x_lower" write (u, descr_fmt) "begin x_upper" do j = 1, self%n_dim - write (u, double_array_fmt) j, self%x_upper(j) + if (pac) then + write (u, double_array_pac_fmt) j, self%x_upper(j) + else + write (u, double_array_fmt) j, self%x_upper(j) + end if end do write (u, descr_fmt) "end x_upper" write (u, descr_fmt) "begin delta_x" do j = 1, self%n_dim - write (u, double_array_fmt) j, self%delta_x(j) + if (pac) then + write (u, double_array_pac_fmt) j, self%delta_x(j) + else + write (u, double_array_fmt) j, self%delta_x(j) + end if end do write (u, descr_fmt) "end delta_x" write (u, descr_fmt) "begin xi" do j = 1, self%n_dim do i = 1, self%n_bins + 1 - write (u, double_array2_fmt) i, j, self%xi(i, j) + if (pac) then + write (u, double_array2_pac_fmt) i, j, self%xi(i, j) + else + write (u, double_array2_fmt) i, j, self%xi(i, j) + end if end do end do write (u, descr_fmt) "end xi" write (u, descr_fmt) "end vegas_grid_t" end subroutine vegas_grid_write @ %def vegas_grid_write @ Compare two grids, if they match up to an given precision. <>= public :: operator (==) <>= interface operator (==) module procedure vegas_grid_equal end interface operator (==) <>= logical function vegas_grid_equal (grid_a, grid_b) result (yorn) type(vegas_grid_t), intent(in) :: grid_a, grid_b yorn = .true. yorn = yorn .and. (grid_a%n_dim == grid_b%n_dim) yorn = yorn .and. (grid_a%n_bins == grid_b%n_bins) yorn = yorn .and. all (grid_a%x_lower == grid_b%x_lower) yorn = yorn .and. all (grid_a%x_upper == grid_b%x_upper) yorn = yorn .and. all (grid_a%delta_x == grid_b%delta_x) yorn = yorn .and. all (grid_a%xi == grid_b%xi) end function vegas_grid_equal @ %def vegas_grid_equal @ Resize each bin accordingly to its corresponding weight [[w]]. Can be used to resize the grid to a new size of bins or refinement. The procedure expects two arguments, firstly, [[n_bins]] and, secondly, the refinement weights [[w]]. If [[n_bins]] differs from the internally stored one, we resize the grid under consideration of [[w]]. If each element of [[w]] equals one, then the bins are resized preserving their original bin density. Anytime else, we refine the grid accordingly to [[w]]. <>= procedure, private :: resize => vegas_grid_resize <>= subroutine vegas_grid_resize (self, n_bins, w) class(vegas_grid_t), intent(inout) :: self integer, intent(in) :: n_bins real(default), dimension(:, :), intent(in) :: w real(default), dimension(size(self%xi)) :: xi_new integer :: i, j, k real(default) :: pts_per_bin real(default) :: d_width do j = 1, self%n_dim if (self%n_bins /= n_bins) then pts_per_bin = real(self%n_bins, default) / real(n_bins, default) self%n_bins = n_bins else if (all (w(:, j) == 0.)) then call msg_bug ("[VEGAS] grid_resize: resize weights are zero.") end if pts_per_bin = sum(w(:, j)) / self%n_bins end if d_width = 0. k = 0 do i = 2, self%n_bins do while (pts_per_bin > d_width) k = k + 1 d_width = d_width + w(k, j) end do d_width = d_width - pts_per_bin associate (x_upper => self%xi(k + 1, j), x_lower => self%xi(k, j)) xi_new(i) = x_upper - (x_upper - x_lower) * d_width / w(k, j) end associate end do self%xi(:, j) = 0. ! Reset grid explicitly self%xi(2:n_bins, j) = xi_new(2:n_bins) self%xi(n_bins + 1, j) = 1. end do end subroutine vegas_grid_resize @ %def vegas_grid_resize @ Find the probability for a given [[x]] in the unit hypercube. For the case [[n_bins < N_BINARY_SEARCH]], we utilize linear search which is faster for short arrays. Else we make use of a binary search. Furthermore, we calculate the inverse of the probability and invert the result only at the end (saving some time on division). <>= procedure, public :: get_probability => vegas_grid_get_probability <>= function vegas_grid_get_probability (self, x) result (g) class(vegas_grid_t), intent(in) :: self real(default), dimension(:), intent(in) :: x integer, parameter :: N_BINARY_SEARCH = 100 real(default) :: g, y integer :: j, i_lower, i_higher, i_mid g = 1. if (self%n_bins > N_BINARY_SEARCH) then g = binary_search (x) else g = linear_search (x) end if ! Move division to the end, which is more efficient. if (g /= 0) g = 1. / g contains real(default) function linear_search (x) result (g) real(default), dimension(:), intent(in) :: x real(default) :: y integer :: j, i g = 1. ndim: do j = 1, self%n_dim y = (x(j) - self%x_lower(j)) / self%delta_x(j) if (y >= 0. .and. y <= 1.) then do i = 2, self%n_bins + 1 if (self%xi(i, j) > y) then g = g * (self%delta_x(j) * & & self%n_bins * (self%xi(i, j) - self%xi(i - 1, j))) cycle ndim end if end do g = 0 exit ndim else g = 0 exit ndim end if end do ndim end function linear_search real(default) function binary_search (x) result (g) real(default), dimension(:), intent(in) :: x ndim: do j = 1, self%n_dim y = (x(j) - self%x_lower(j)) / self%delta_x(j) if (y >= 0. .and. y <= 1.) then i_lower = 1 i_higher = self%n_bins + 1 search: do if (i_lower >= (i_higher - 1)) then g = g * (self%delta_x(j) * & & self%n_bins * (self%xi(i_higher, j) - self%xi(i_higher - 1, j))) cycle ndim end if i_mid = (i_higher + i_lower) / 2 if (y > self%xi(i_mid, j)) then i_lower = i_mid else i_higher = i_mid end if end do search else g = 0. exit ndim end if end do ndim end function binary_search end function vegas_grid_get_probability @ %def vegas_grid_get_probability @ Broadcast the grid information. As safety measure, we get the actual grid object from VEGAS (correclty allocated, but for non-root unfilled) and broadcast the root object. On success, we set grid into VEGAS. We use the non-blocking broadcast routine, because we have to send quite a bunch of integers and reals. We have to be very careful with [[n_bins]], the number of bins can actually change during different iterations. If we reuse a grid, we have to take that, every grid uses the same [[n_bins]]. We expect, that the number of dimension does not change, which is in principle possible, but will be checked onto in [[vegas_set_grid]]. <>= procedure, public :: broadcast => vegas_grid_broadcast <>= subroutine vegas_grid_broadcast (self) class(vegas_grid_t), intent(inout) :: self integer :: j, ierror type(MPI_Request), dimension(self%n_dim + 4) :: status ! Blocking call MPI_Bcast (self%n_bins, 1, MPI_INTEGER, 0, MPI_COMM_WORLD) ! Non blocking call MPI_Ibcast (self%n_dim, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, status(1)) call MPI_Ibcast (self%x_lower, self%n_dim, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, status(2)) call MPI_Ibcast (self%x_upper, self%n_dim, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, status(3)) call MPI_Ibcast (self%delta_x, self%n_dim, & & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, status(4)) ndim: do j = 1, self%n_dim call MPI_Ibcast (self%xi(1:self%n_bins + 1, j), self%n_bins + 1,& & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, status(4 + j)) end do ndim call MPI_Waitall (self%n_dim + 4, status, MPI_STATUSES_IGNORE) end subroutine vegas_grid_broadcast @ %def vegas_grid_broadcast @ \section{Type: vegas\_result\_t} \label{sec:type:-vegas_r} We store the result of the latest iteration(s) in a transparent container. The [[vegas_result_t]] object inside VEGAS must not be directly accessible. We export the a copy of the result via a get-method of the [[vegas_t]] object. We store latest event weight in [[evt_weight]] and a (possible) evebt weight excess in [[evt_weight_excess]], if the event weight is larger than [[max_abs_f]]. <>= public :: vegas_result_t <>= type :: vegas_result_t integer :: it_start = 0 integer :: it_num = 0 integer :: samples = 0 real(default) :: sum_int_wgtd = 0. real(default) :: sum_wgts real(default) :: sum_chi = 0. real(default) :: chi2 = 0. real(default) :: efficiency = 0. real(default) :: efficiency_pos = 0. real(default) :: efficiency_neg = 0. real(default) :: max_abs_f = 0. real(default) :: max_abs_f_pos = 0. real(default) :: max_abs_f_neg = 0. real(default) :: result = 0. real(default) :: std = 0. real(default) :: evt_weight = 0. real(default) :: evt_weight_excess = 0. contains <> end type vegas_result_t @ %def vegas_results_t @ Write out the current status of the integration result. <>= procedure, public :: write => vegas_result_write <>= subroutine vegas_result_write (self, unit, indent) class(vegas_result_t), intent(in) :: self integer, intent(in), optional :: unit integer, intent(in), optional :: indent integer :: u, ind u = given_output_unit (unit) ind = 0; if (present (indent)) ind = indent call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Start iteration = ", self%it_start call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Iteration number = ", self%it_num call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Sample number = ", self%samples call write_indent (u, ind) write (u, "(2x,A," // FMT_17 //")") & & "Sum of weighted integrals = ", self%sum_int_wgtd call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Sum of weights = ", self%sum_wgts call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Sum of chi = ", self%sum_chi call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "chi2 = ", self%chi2 call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Overall efficiency = ", self%efficiency call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "f-positive efficiency = ", self%efficiency_pos call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "f-negative efficiency = ", self%efficiency_neg call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Maximum absolute overall value = ", self%max_abs_f call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Maximum absolute positive value = ", self%max_abs_f_pos call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Maximum absolute negative value = ", self%max_abs_f_neg call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Integral (of latest iteration) = ", self%result call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Standard deviation = ", self%std write (u, "(2x,A," // FMT_17 // ")") & & "Event weight = ", self%evt_weight write (u, "(2x,A," // FMT_17 // ")") & & "Event weight excess = ", self%evt_weight_excess end subroutine vegas_result_write @ %def vegas_results_write @ Send the result object to specified rank, internally in a non-blocking way. We do not need to handle the event results, because each event result is atomic. <>= procedure, public :: send => vegas_result_send <>= subroutine vegas_result_send (self, receiver, tag) class(vegas_result_t), intent(in) :: self integer, intent(in) :: receiver integer, intent(in) :: tag type(MPI_Request), dimension(13) :: request call MPI_Isend (self%it_start, 1, MPI_INTEGER, receiver, 1 + tag,& & MPI_COMM_WORLD, request(1)) call MPI_Isend (self%it_num, 1, MPI_INTEGER, receiver , 2 + tag,& & MPI_COMM_WORLD, request(2)) call MPI_Isend (self%samples, 1, MPI_INTEGER, receiver, 3 + tag,& & MPI_COMM_WORLD, request(3)) call MPI_Isend (self%sum_int_wgtd, 1, MPI_DOUBLE_PRECISION, receiver, 4 +& & tag, MPI_COMM_WORLD, request(4)) call MPI_Isend (self%sum_wgts, 1, MPI_DOUBLE_PRECISION, receiver, 5 + tag,& & MPI_COMM_WORLD, request(5)) call MPI_Isend (self%sum_chi, 1, MPI_DOUBLE_PRECISION, receiver, 6 + tag,& & MPI_COMM_WORLD, request(6)) call MPI_Isend (self%efficiency, 1, MPI_DOUBLE_PRECISION, receiver, 7 + tag& &, MPI_COMM_WORLD, request(7)) call MPI_Isend (self%efficiency_pos, 1, MPI_DOUBLE_PRECISION, receiver, 8 +& & tag, MPI_COMM_WORLD, request(8)) call MPI_Isend (self%efficiency_neg, 1, MPI_DOUBLE_PRECISION, receiver, 9 +& & tag, MPI_COMM_WORLD, request(9)) call MPI_Isend (self%max_abs_f, 1, MPI_DOUBLE_PRECISION, receiver, 10 + tag& &, MPI_COMM_WORLD, request(10)) call MPI_Isend (self%max_abs_f_pos, 1, MPI_DOUBLE_PRECISION, receiver, 11 +& & tag, MPI_COMM_WORLD, request(10)) call MPI_Isend (self%max_abs_f_neg, 1, MPI_DOUBLE_PRECISION, receiver, 12 +& & tag, MPI_COMM_WORLD, request(11)) call MPI_Isend (self%result, 1, MPI_DOUBLE_PRECISION, receiver, 13 + tag,& & MPI_COMM_WORLD, request(12)) call MPI_Isend (self%std, 1, MPI_DOUBLE_PRECISION, receiver, 14 + tag,& & MPI_COMM_WORLD, request(13)) call MPI_waitall (13, request, MPI_STATUSES_IGNORE) end subroutine vegas_result_send @ %def vegas_result_communicate @ Receive the result object from a specified rank, internally in a non-blocking way. <>= procedure, public :: receive => vegas_result_receive <>= subroutine vegas_result_receive (self, sender, tag) class(vegas_result_t), intent(inout) :: self integer, intent(in) :: sender integer, intent(in) :: tag type(MPI_Request), dimension(13) :: request call MPI_Irecv (self%it_start, 1, MPI_INTEGER, sender, 1 + tag,& & MPI_COMM_WORLD, request(1)) call MPI_Irecv (self%it_num, 1, MPI_INTEGER, sender , 2 + tag,& & MPI_COMM_WORLD, request(2)) call MPI_Irecv (self%samples, 1, MPI_INTEGER, sender, 3 + tag,& & MPI_COMM_WORLD, request(3)) call MPI_Irecv (self%sum_int_wgtd, 1, MPI_DOUBLE_PRECISION, sender, 4 + tag& &, MPI_COMM_WORLD, request(4)) call MPI_Irecv (self%sum_wgts, 1, MPI_DOUBLE_PRECISION, sender, 5 + tag,& & MPI_COMM_WORLD, request(5)) call MPI_Irecv (self%sum_chi, 1, MPI_DOUBLE_PRECISION, sender, 6 + tag,& & MPI_COMM_WORLD, request(6)) call MPI_Irecv (self%efficiency, 1, MPI_DOUBLE_PRECISION, sender, 7 + tag,& & MPI_COMM_WORLD, request(7)) call MPI_Irecv (self%efficiency_pos, 1, MPI_DOUBLE_PRECISION, sender, 8 +& & tag, MPI_COMM_WORLD, request(8)) call MPI_Irecv (self%efficiency_neg, 1, MPI_DOUBLE_PRECISION, sender, 9 +& & tag, MPI_COMM_WORLD, request(9)) call MPI_Irecv (self%max_abs_f, 1, MPI_DOUBLE_PRECISION, sender, 10 + tag,& & MPI_COMM_WORLD, request(10)) call MPI_Irecv (self%max_abs_f_pos, 1, MPI_DOUBLE_PRECISION, sender, 11 +& & tag, MPI_COMM_WORLD, request(10)) call MPI_Irecv (self%max_abs_f_neg, 1, MPI_DOUBLE_PRECISION, sender, 12 +& & tag, MPI_COMM_WORLD, request(11)) call MPI_Irecv (self%result, 1, MPI_DOUBLE_PRECISION, sender, 13 + tag,& & MPI_COMM_WORLD, request(12)) call MPI_Irecv (self%std, 1, MPI_DOUBLE_PRECISION, sender, 14 + tag,& & MPI_COMM_WORLD, request(13)) call MPI_waitall (13, request, MPI_STATUSES_IGNORE) end subroutine vegas_result_receive @ %def vegas_result_receive \section{Type: vegas\_t} \label{sec:type:-vegas_t} The VEGAS object contains the methods for integration and grid resize- and refinement. We store the grid configuration and the (current) result in transparent containers alongside with the actual grid and the distribution. The values of the distribution depend on the chosen mode whether the function value or the variance is stored. The distribution is used after each iteration to refine the grid. <>= public :: vegas_t <>= type :: vegas_t private type(vegas_config_t) :: config real(default) :: hypercube_volume = 0. real(default) :: jacobian = 0. real(default), dimension(:, :), allocatable :: d type(vegas_grid_t) :: grid integer, dimension(:), allocatable :: bin integer, dimension(:), allocatable :: box type(vegas_result_t) :: result contains <> end type vegas_t @ %def vegas_t @ We overload the type constructor of [[vegas_t]] which initialises the mandatory argument [[n_dim]] and allocate memory for the grid. <>= interface vegas_t module procedure vegas_init end interface vegas_t <>= type(vegas_t) function vegas_init (n_dim, alpha, n_bins_max, iterations, mode) result (self) integer, intent(in) :: n_dim integer, intent(in), optional :: n_bins_max real(default), intent(in), optional :: alpha integer, intent(in), optional :: iterations integer, intent(in), optional :: mode self%config%n_dim = n_dim if (present (alpha)) self%config%alpha = alpha if (present (n_bins_max)) self%config%n_bins_max = n_bins_max if (present (iterations)) self%config%iterations = iterations if (present (mode)) self%config%mode = mode self%grid = vegas_grid_t (n_dim, self%config%n_bins_max) allocate (self%d(self%config%n_bins_max, n_dim), source=0.0_default) allocate (self%box(n_dim), source=1) allocate (self%bin(n_dim), source=1) self%config%n_bins = 1 self%config%n_boxes = 1 call self%set_limits (self%grid%x_lower, self%grid%x_upper) call self%reset_grid () call self%reset_result () end function vegas_init @ %def vegas_init @ Finalize the grid. Deallocate grid memory. <>= procedure, public :: final => vegas_final <>= subroutine vegas_final (self) class(vegas_t), intent(inout) :: self deallocate (self%grid%x_upper) deallocate (self%grid%x_lower) deallocate (self%grid%delta_x) deallocate (self%d) deallocate (self%grid%xi) deallocate (self%box) deallocate (self%bin) end subroutine vegas_final @ %def vegas_final \section{Get-/Set-methods} \label{sec:set-get-methods} @ The VEGAS object prohibits direct access from outside. Communication is handle via get- or set-methods. Set the limits of integration. The defaults limits correspong the $n$-dimensionl unit hypercube. \textit{Remark:} After setting the limits, the grid is initialised, again. Previous results are lost due to recalculation of the overall jacobian. <>= procedure, public :: set_limits => vegas_set_limits <>= subroutine vegas_set_limits (self, x_lower, x_upper) class(vegas_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x_lower real(default), dimension(:), intent(in) :: x_upper if (size (x_lower) /= self%config%n_dim & & .or. size (x_upper) /= self%config%n_dim) then write (msg_buffer, "(A, I5, A, I5, A, I5)") & "VEGAS: [set_limits] n_dim of new lower/upper integration limit& & does not match previously set n_dim. ", self%config%n_dim, " =/=& & ", size (x_lower), " =/= ", size (x_upper) call msg_fatal () end if if (any(x_upper < x_lower)) then call msg_fatal ("VEGAS: [set_limits] upper limits are smaller than lower limits.") end if if (any((x_upper - x_lower) > huge(0._default))) then call msg_fatal ("VEGAS: [set_limits] upper and lower limits exceed rendering.") end if self%grid%x_upper = x_upper self%grid%x_lower = x_lower self%grid%delta_x = self%grid%x_upper - self%grid%x_lower self%hypercube_volume = product (self%grid%delta_x) call self%reset_result () end subroutine vegas_set_limits @ %def vegas_set_limits @ Set the number of calls. If the number of calls changed during different passes, we resize the grid preserving the probability density. We should reset the results after changing the number of calls which change the size of the grid and the running mode of VEGAS. But, this is a set method only for the number of calls. <>= procedure, public :: set_calls => vegas_set_n_calls <>= subroutine vegas_set_n_calls (self, n_calls) class(vegas_t), intent(inout) :: self integer, intent(in) :: n_calls if (.not. (n_calls > 0)) then write (msg_buffer, "(A, I5)") & "VEGAS: [set_calls] number of calls must be a positive number. Keep& & number of calls = ", self%config%n_calls call msg_warning () else self%config%n_calls = max (n_calls, self%config%n_calls_min) if (self%config%n_calls /= n_calls) then write (msg_buffer, "(A,I5)") & "VEGAS: [set calls] number of calls is too few, reset to ", self%config%n_calls call msg_warning () end if call self%init_grid () end if end subroutine vegas_set_n_calls @ %def vegas_set_n_calls @ Get the grid object and set [[n_bins]], [[n_dim]] inside grid container. <>= procedure, public :: get_grid => vegas_get_grid <>= type(vegas_grid_t) function vegas_get_grid (self) result (grid) class(vegas_t), intent(in) :: self grid = self%grid grid%n_dim = self%config%n_dim grid%n_bins = self%config%n_bins end function vegas_get_grid @ %def vegas_get_grid @ Set grid. We need a set method for the parallelisation. We do some additional checks before copying the object. Be careful, we do not check on [[n_bins]], because the number of bins can change after setting [[n_calls]]. We remind you, that you will loose all your current progress, if you use set the grid. Hence, it will only be used when compiled with [[MPI]]. <>= procedure, public :: set_grid => vegas_set_grid <>= subroutine vegas_set_grid (self, grid) class(vegas_t), intent(inout) :: self type(vegas_grid_t), intent(in) :: grid integer :: j, rank logical :: success call MPI_Comm_rank (MPI_COMM_WORLD, rank) success = .true. success = (success .and. (grid%n_dim .eq. self%config%n_dim)) success = (success .and. all (grid%x_lower .eq. self%grid%x_lower)) success = (success .and. all (grid%x_upper .eq. self%grid%x_upper)) success = (success .and. all (grid%delta_x .eq. self%grid%delta_x)) if (success) then self%config%n_bins = grid%n_bins do j = 1, self%config%n_dim self%grid%xi(1, j) = 0._default self%grid%xi(2:self%config%n_bins, j) = grid%xi(2:grid%n_bins, j) self%grid%xi(self%config%n_bins + 1, j) = 1._default end do else call msg_bug ("VEGAS: set grid: boundary conditions do not match.") end if end subroutine vegas_set_grid @ %def vegas_set_grid @ We check if it is senseful to parallelize the actual grid. In simple, this means that [[n_boxes]] has to be larger than 2. With the result that we could have an actual superimposed stratified grid. In advance, we can give the size of communicator [[n_size]] and check whether we have enough boxes to distribute. <>= procedure, public :: is_parallelizable => vegas_is_parallelizable <>= elemental logical function vegas_is_parallelizable (self, opt_n_size) result (flag) class(vegas_t), intent(in) :: self integer, intent(in), optional :: opt_n_size integer :: n_size n_size = 2 if (present (opt_n_size)) n_size = opt_n_size flag = (self%config%n_boxes**floor (self%config%n_dim / 2.) >= n_size) end function vegas_is_parallelizable @ %def vegas_is_parallelizable @ Get the config object. <>= procedure, public :: get_config => vegas_get_config <>= subroutine vegas_get_config (self, config) class(vegas_t), intent(in) :: self type(vegas_config_t), intent(out) :: config config = self%config end subroutine vegas_get_config @ %def vegas_get_config @ Set non-runtime dependent configuration. It will no be possible to change [[n_bins_max]]. <>= procedure, public :: set_config => vegas_set_config <>= subroutine vegas_set_config (self, config) class(vegas_t), intent(inout) :: self class(vegas_config_t), intent(in) :: config self%config%alpha = config%alpha self%config%iterations = config%iterations self%config%mode = config%mode self%config%n_calls_min = config%n_calls_min end subroutine vegas_set_config @ %def vegas_set_config @ Get the result object. <>= procedure, public :: get_result => vegas_get_result <>= type(vegas_result_t) function vegas_get_result (self) result (result) class(vegas_t), intent(in) :: self result = self%result end function vegas_get_result @ %def vegas_get_result @ Set the result object. Be reminded, that you will loose your current results, if you are not careful! Hence, it will only be avaible during usage with [[MPI]]. <>= procedure, public :: set_result => vegas_set_result <>= subroutine vegas_set_result (self, result) class(vegas_t), intent(inout) :: self type(vegas_result_t), intent(in) :: result self%result = result end subroutine vegas_set_result @ %def vegas_set_result @ Get (actual) number of calls. <>= procedure, public :: get_calls => vegas_get_n_calls <>= elemental real(default) function vegas_get_n_calls (self) result (n_calls) class(vegas_t), intent(in) :: self n_calls = self%config%n_calls end function vegas_get_n_calls @ %def vegas_get_n_calls @ Get the cumulative result of the integration. Recalculate weighted average of the integration. <>= procedure, public :: get_integral => vegas_get_integral <>= elemental real(default) function vegas_get_integral (self) result (integral) class(vegas_t), intent(in) :: self integral = 0. if (self%result%sum_wgts > 0.) then integral = self%result%sum_int_wgtd / self%result%sum_wgts end if end function vegas_get_integral @ %def vegas_get_integral @ Get the cumulative variance of the integration. Recalculate the variance. <>= procedure, public :: get_variance => vegas_get_variance <>= elemental real(default) function vegas_get_variance (self) result (variance) class(vegas_t), intent(in) :: self variance = 0. if (self%result%sum_wgts > 0.) then variance = 1.0 / self%result%sum_wgts end if end function vegas_get_variance @ %def vegas_get_variance @ Get efficiency. <>= procedure, public :: get_efficiency => vegas_get_efficiency <>= elemental real(default) function vegas_get_efficiency (self) result (efficiency) class(vegas_t), intent(in) :: self efficiency = 0. if (self%result%efficiency > 0. ) then efficiency = self%result%efficiency end if end function vegas_get_efficiency @ %def vegas_get_efficiency @ Get [[f_max]]. <>= procedure, public :: get_max_abs_f => vegas_get_max_abs_f <>= elemental real(default) function vegas_get_max_abs_f (self) result (max_abs_f) class(vegas_t), intent(in) :: self max_abs_f = 0. if (self%result%max_abs_f > 0.) then max_abs_f = self%result%max_abs_f end if end function vegas_get_max_abs_f @ %def vegas_get_max_abs_f @ Get [[f_max_pos]]. <>= procedure, public :: get_max_abs_f_pos => vegas_get_max_abs_f_pos <>= elemental real(default) function vegas_get_max_abs_f_pos (self) result (max_abs_f) class(vegas_t), intent(in) :: self max_abs_f = 0. if (self%result%max_abs_f_pos > 0.) then max_abs_f = self%result%max_abs_f_pos end if end function vegas_get_max_abs_f_pos @ %def vegas_get_max_abs_f_pos @ Get [[f_max_neg]]. <>= procedure, public :: get_max_abs_f_neg => vegas_get_max_abs_f_neg <>= elemental real(default) function vegas_get_max_abs_f_neg (self) result (max_abs_f) class(vegas_t), intent(in) :: self max_abs_f = 0. if (self%result%max_abs_f_neg > 0.) then max_abs_f = self%result%max_abs_f_neg end if end function vegas_get_max_abs_f_neg @ %def vegas_get_max_abs_f_neg @ Get event weight and excess. <>= procedure, public :: get_evt_weight => vegas_get_evt_weight procedure, public :: get_evt_weight_excess => vegas_get_evt_weight_excess <>= real(default) function vegas_get_evt_weight (self) result (evt_weight) class(vegas_t), intent(in) :: self evt_weight = self%result%evt_weight end function vegas_get_evt_weight real(default) function vegas_get_evt_weight_excess (self) result (evt_weight_excess) class(vegas_t), intent(in) :: self evt_weight_excess = self%result%evt_weight_excess end function vegas_get_evt_weight_excess @ %def vegas_get_evt_weight, vegas_get_evt_weight_excess @ Get and set distribution. Beware! This method is hideous as it allows to manipulate the algorithm at its very core. <>= procedure, public :: get_distribution => vegas_get_distribution procedure, public :: set_distribution => vegas_set_distribution <>= function vegas_get_distribution (self) result (d) class(vegas_t), intent(in) :: self real(default), dimension(:, :), allocatable :: d d = self%d end function vegas_get_distribution subroutine vegas_set_distribution (self, d) class(vegas_t), intent(inout) :: self real(default), dimension(:, :), intent(in) :: d if (size (d, dim = 2) /= self%config%n_dim) then call msg_bug ("[VEGAS] set_distribution: new distribution has wrong size of dimension") end if if (size (d, dim = 1) /= self%config%n_bins_max) then call msg_bug ("[VEGAS] set_distribution: new distribution has wrong number of bins") end if self%d = d end subroutine vegas_set_distribution @ %def vegas_set_distribution, vegas_get_distribution @ Send distribution to specified rank, internally in a non-blocking way. We send the complete array of [[d]], not just the actually used part. <>= procedure, public :: send_distribution => vegas_send_distribution <>= subroutine vegas_send_distribution (self, receiver, tag) class(vegas_t), intent(in) :: self integer, intent(in) :: receiver integer, intent(in) :: tag integer :: j type(MPI_Request), dimension(self%config%n_dim + 2) :: request call MPI_Isend (self%bin, self%config%n_dim, MPI_INTEGER, receiver, tag + 1& &, MPI_COMM_WORLD, request(1)) call MPI_Isend (self%box, self%config%n_dim, MPI_INTEGER, receiver, tag + 2& &, MPI_COMM_WORLD, request(2)) do j = 1, self%config%n_dim call MPI_Isend (self%d(:, j), self%config%n_bins_max,& & MPI_DOUBLE_PRECISION, receiver, tag + j + 2, MPI_COMM_WORLD,& & request(j + 2)) end do call MPI_Waitall (self%config%n_dim, request, MPI_STATUSES_IGNORE) end subroutine vegas_send_distribution @ %def vegas_send_distribution @ Receive distribution from specified rank, internally in a non-blocking way. <>= procedure, public :: receive_distribution => vegas_receive_distribution <>= subroutine vegas_receive_distribution (self, sender, tag) class(vegas_t), intent(inout) :: self integer, intent(in) :: sender integer, intent(in) :: tag integer :: j type(MPI_Request), dimension(self%config%n_dim + 2) :: request call MPI_Irecv (self%bin, self%config%n_dim, MPI_INTEGER, sender, tag + 1& &, MPI_COMM_WORLD, request(1)) call MPI_Irecv (self%box, self%config%n_dim, MPI_INTEGER, sender, tag + 2& &, MPI_COMM_WORLD, request(2)) do j = 1, self%config%n_dim call MPI_Irecv (self%d(:, j), self%config%n_bins_max,& & MPI_DOUBLE_PRECISION, sender, tag + j + 2, MPI_COMM_WORLD,& & request(j + 2)) end do call MPI_Waitall (self%config%n_dim, request, MPI_STATUSES_IGNORE) end subroutine vegas_receive_distribution @ %def vegas_receive_distribution \section{Grid resize- and refinement} \label{sec:grid-resize-refin} Before integration the grid itself must be initialised. Given the number of [[n_calls]] and [[n_dim]] we prepare the grid for the integration. The grid is binned according to the VEGAS mode and [[n_calls]]. If the mode is not set to [[vegas_importance_only]], the grid is divided in to equal boxes. We try for 2 calls per box \begin{equation} boxes = \sqrt[n_{dim}]{\frac{calls}{2}}. \end{equation} If the numbers of boxes exceeds the number of bins, which is the case for low dimensions, the algorithm switches to stratified sampling. Otherwise, we are still using importance sampling, but keep the boxes for book keeping. If the number of bins changes from the previous invocation, bins are expanded or contracted accordingly, while preserving bin density. <>= procedure, private :: init_grid => vegas_init_grid <>= subroutine vegas_init_grid (self) class(vegas_t), intent(inout) :: self integer :: n_bins, n_boxes, box_per_bin, n_total_boxes real(default), dimension(:, :), allocatable :: w n_bins = self%config%n_bins_max n_boxes = 1 if (self%config%mode /= VEGAS_MODE_IMPORTANCE_ONLY) then ! We try for 2 calls per box n_boxes = max (floor ((self%config%n_calls / 2.)**(1. / self%config%n_dim)), 1) self%config%mode = VEGAS_MODE_IMPORTANCE if (2 * n_boxes >= self%config%n_bins_max) then ! if n_bins/box < 2 box_per_bin = max (n_boxes / self%config%n_bins_max, 1) n_bins = min (n_boxes / box_per_bin, self%config%n_bins_max) n_boxes = box_per_bin * n_bins self%config%mode = VEGAS_MODE_STRATIFIED end if end if n_total_boxes = n_boxes**self%config%n_dim self%config%calls_per_box = max (floor (real (self%config%n_calls) / n_total_boxes), 2) self%config%n_calls = self%config%calls_per_box * n_total_boxes ! Total volume of x-space/(average n_calls per bin) self%jacobian = self%hypercube_volume * real(n_bins, default)& &**self%config%n_dim / real(self%config%n_calls, default) self%config%n_boxes = n_boxes if (n_bins /= self%config%n_bins) then allocate (w(self%config%n_bins, self%config%n_dim), source=1.0_default) call self%grid%resize (n_bins, w) self%config%n_bins = n_bins end if end subroutine vegas_init_grid @ %def vegas_init_grid @ Reset the cumulative result, and efficiency and max. grid values. <>= procedure, public :: reset_result => vegas_reset_result <>= subroutine vegas_reset_result (self) class(vegas_t), intent(inout) :: self self%result%sum_int_wgtd = 0. self%result%sum_wgts = 0. self%result%sum_chi = 0. self%result%it_num = 0 self%result%samples = 0 self%result%chi2 = 0 self%result%efficiency = 0. self%result%efficiency_pos = 0. self%result%efficiency_neg = 0. self%result%max_abs_f = 0. self%result%max_abs_f_pos = 0. self%result%max_abs_f_neg = 0. end subroutine vegas_reset_result @ %def vegas_reset_results @ Reset the grid. Purge the adapted grid and the distribution. Furthermore, reset the results. The maximal size of the grid remains. Note: Handle [[vegas_reset_grid]] with great care! Instead of reusing an old object, create a new one. <>= procedure, public :: reset_grid => vegas_reset_grid <>= subroutine vegas_reset_grid (self) class(vegas_t), intent(inout) :: self self%config%n_bins = 1 self%d = 0._default self%grid%xi = 0._default self%grid%xi(1, :) = 0._default self%grid%xi(2, :) = 1._default call self%reset_result () end subroutine vegas_reset_grid @ %def vegas_reset_grid @ Refine the grid to match the distribution [[d]]. Average the distribution over neighbouring bins, then contract or expand the bins. The averaging dampens high fluctuations amog the integrand or the variance. We make the type-bound procedure public accessible because the multi-channel integration refines each grid after integration over all grids. <>= procedure, public :: refine => vegas_refine_grid <>= subroutine vegas_refine_grid (self) class(vegas_t), intent(inout) :: self integer :: j real(default), dimension(self%config%n_bins, self%config%n_dim) :: w ndim: do j = 1, self%config%n_dim call average_distribution (self%config%n_bins, self%d(:self%config& &%n_bins, j), self%config%alpha, w(:, j)) end do ndim call self%grid%resize (self%config%n_bins, w) contains <> end subroutine vegas_refine_grid @ %def vegas_refine_grid @ We average the collected values [[d]] of the (sq.) weighted [[f]] over neighbouring bins. The averaged [[d]] are then agian damped by a logarithm to enhance numerical stability. The results are then the refinement weights [[w]]. We have to take care of the special case where we have a very low sampling rate. In those cases we can not be sure that the distribution is satisfying filled, although we have already averaged over neighbouring bins. This will lead to a squashing of the unfilled bins and such the boundaries of those will be pushed together. We circumvent this problem by setting those unfilled bins to the smallest representable value of a default real. The problem becomes very annoying in the multi-channel formualae where have to look up via binary search the corresponding probability of [[x]] and if the width is zero, the point will be neglected. <>= subroutine average_distribution (n_bins, d, alpha, w) integer, intent(in) :: n_bins real(default), dimension(:), intent(inout) :: d real(default), intent(in) :: alpha real(default), dimension(n_bins), intent(out) :: w if (n_bins > 2) then d(1) = (d(1) + d(2)) / 2.0_default d(2:n_bins - 1) = (d(1:n_bins - 2) + d(2:n_bins - 1) + d(3:n_bins)) /& & 3.0_default d(n_bins) = d(n_bins - 1) + d(n_bins) / 2.0_default end if w = 1.0_default if (.not. all (d < tiny (1.0_default))) then d = d / sum (d) where (d < tiny (1.0_default)) d = tiny (1.0_default) end where where (d /= 1.0_default) w = ((d - 1.) / log(d))**alpha elsewhere ! Analytic limes for d -> 1 w = 1.0_default end where end if end subroutine average_distribution @ %def average_distribution @ \section{Integration} \label{sec:integration} Integrate [[func]], in the previous set bounds [[x_lower]] to [[x_upper]], with [[n_calls]]. Use results from previous invocations of [[integrate]] with [[opt_reset_result = .false.]] and better with subsequent calls. Before we walk through the hybercube, we initialise the grid (at a central position). We step through the (equidistant) boxes which ensure we do not miss any place in the n-dim. hypercube. In each box we sample [[calls_per_box]] random points and transform them to bin coordinates. The total integral and the total (sample) variance over each box $i$ is then calculated by \begin{align*} E(I)_{i} = \sum_{j}^{\text{calls per box}} I_{i, j}, \\ V(I)_{i} = \text{calls per box} \frac{\sum_{j}^{\text{calls per box}}} I_{i, j}^{2} - (\sum_{j}^{\text{calls per box}} I_{i, j})**2 \frac{\text{calls per box}}{\text{calls per box} - 1}. \end{align*} The stratification of the $n$-dimensional hybercube allows a simple parallelisation of the algorithm (R. Kreckel, "Parallelization of adaptive MC integrators", Computer Physics Communications, vol. 106, no. 3, pp. 258–266, Nov. 1997.). We have to ensure that all boxes are sampled, but the number of boxes to distribute is too large. We allow each thread to sample a fraction $r$ of all boxes $k$ such that $r \ll k$. Furthermore, we constrain that the number of process $p$ is much smaller than $r$. The overall constraint is \begin{equation} p \ll r \ll k. \end{equation} We divide the intgeration into a parallel and a perpendicular subspace. The number of parallel dimensions is $D_{\parallel} = \lfloor \frac{D}{2} \rfloor$. <>= procedure, public :: integrate => vegas_integrate <>= subroutine vegas_integrate (self, func, rng, iterations, opt_reset_result,& & opt_refine_grid, opt_verbose, result, abserr) class(vegas_t), intent(inout) :: self class(vegas_func_t), intent(inout) :: func class(rng_t), intent(inout) :: rng integer, intent(in), optional :: iterations logical, intent(in), optional :: opt_reset_result logical, intent(in), optional :: opt_refine_grid logical, intent(in), optional :: opt_verbose real(default), optional, intent(out) :: result, abserr integer :: it, j, k real(default), dimension(self%config%n_dim) :: x real(default) :: fval, fval_sq, bin_volume real(default) :: fval_box, fval_sq_box real(default) :: total_integral, total_sq_integral, total_variance, chi, wgt real(default) :: cumulative_int, cumulative_std real(default) :: sum_abs_f_pos, max_abs_f_pos real(default) :: sum_abs_f_neg, max_abs_f_neg logical :: reset_result = .true. logical :: refine_grid = .true. logical :: verbose = .false. <> if (present (iterations)) self%config%iterations = iterations if (present (opt_reset_result)) reset_result = opt_reset_result if (present (opt_refine_grid)) refine_grid = opt_refine_grid if (present (opt_verbose)) verbose = opt_verbose <> if (verbose) then call msg_message ("Results: [it, calls, integral, error, chi^2, eff.]") end if iteration: do it = 1, self%config%iterations <> loop_over_par_boxes: do while (box_success) loop_over_perp_boxes: do while (box_success) fval_box = 0._default fval_sq_box = 0._default do k = 1, self%config%calls_per_box call self%random_point (rng, x, bin_volume) ! Call the function, yeah, call it... fval = self%jacobian * bin_volume * func%evaluate (x) fval_sq = fval**2 fval_box = fval_box + fval fval_sq_box = fval_sq_box + fval_sq if (fval > 0.) then max_abs_f_pos = max(abs (fval), max_abs_f_pos) sum_abs_f_pos = sum_abs_f_pos + abs(fval) else max_abs_f_neg = max(abs (fval), max_abs_f_neg) sum_abs_f_neg = sum_abs_f_neg + abs(fval) end if if (self%config%mode /= VEGAS_MODE_STRATIFIED) then call self%accumulate_distribution (fval_sq) end if end do fval_sq_box = sqrt (fval_sq_box * self%config%calls_per_box) ! (a - b) * (a + b) = a**2 - b**2 fval_sq_box = (fval_sq_box - fval_box) * (fval_sq_box + fval_box) if (fval_sq_box <= 0.0) fval_sq_box = tiny (1.0_default) total_integral = total_integral + fval_box total_sq_integral = total_sq_integral + fval_sq_box if (self%config%mode == VEGAS_MODE_STRATIFIED) then call self%accumulate_distribution (fval_sq_box) end if call increment_box_coord (self%box(n_dim_par + 1:self%config& &%n_dim), box_success) end do loop_over_perp_boxes shift: do k = 1, n_size call increment_box_coord (self%box(1:n_dim_par), box_success) if (.not. box_success) exit shift end do shift <> end do loop_over_par_boxes <> associate (result => self%result) ! Compute final results for this iterations total_variance = total_sq_integral / (self%config%calls_per_box - 1.) ! Ensure variance is always positive and larger than zero. if (total_variance < tiny (1._default) / epsilon (1._default) & & * max (total_integral**2, 1._default)) then total_variance = tiny (1._default) / epsilon (1._default) & & * max (total_integral**2, 1._default) end if wgt = 1. / total_variance total_sq_integral = total_integral**2 result%result = total_integral result%std = sqrt (total_variance) result%samples = result%samples + 1 if (result%samples == 1) then result%chi2 = 0._default else chi = total_integral if (result%sum_wgts > 0) then chi = chi - result%sum_int_wgtd / result%sum_wgts end if result%chi2 = result%chi2 * (result%samples - 2.0_default) result%chi2 = (wgt / (1._default + (wgt / result%sum_wgts))) & & * chi**2 result%chi2 = result%chi2 / (result%samples - 1._default) end if result%sum_wgts = result%sum_wgts + wgt result%sum_int_wgtd = result%sum_int_wgtd + (total_integral * wgt) result%sum_chi = result%sum_chi + (total_sq_integral * wgt) cumulative_int = result%sum_int_wgtd / result%sum_wgts cumulative_std = sqrt (1. / result%sum_wgts) end associate call calculate_efficiency () if (verbose) then - write (msg_buffer, "(I0,1x,I0,1x, 4(E16.8E4,1x))") & + write (msg_buffer, "(I0,1x,I0,1x, 4(E24.16E4,1x))") & & it, self%config%n_calls, cumulative_int, cumulative_std, & & self%result%chi2, self%result%efficiency call msg_message () end if if (refine_grid) call self%refine () end do iteration if (present(result)) result = cumulative_int if (present(abserr)) abserr = abs(cumulative_std) contains <> end subroutine vegas_integrate @ %def vegas_integrate @ Calculate the extras here. We define \begin{equation} \operatorname*{max}_{x} w(x) = \frac{f(x)}{p(x)} \Delta_{\text{jac}}. \end{equation} In the implementation we have to factor out [[n_calls]] in the jacobian. Also, during event generation. <>= subroutine calculate_efficiency () self%result%max_abs_f_pos = self%config%n_calls * max_abs_f_pos self%result%max_abs_f_neg = self%config%n_calls * max_abs_f_neg self%result%max_abs_f = & & max (self%result%max_abs_f_pos, self%result%max_abs_f_neg) self%result%efficiency_pos = 0. if (max_abs_f_pos > 0.) then self%result%efficiency_pos = & & sum_abs_f_pos / max_abs_f_pos end if self%result%efficiency_neg = 0. if (max_abs_f_neg > 0.) then self%result%efficiency_neg = & & sum_abs_f_neg / max_abs_f_neg end if self%result%efficiency = 0. if (self%result%max_abs_f > 0.) then self%result%efficiency = (sum_abs_f_pos + sum_abs_f_neg) & & / self%result%max_abs_f end if end subroutine calculate_efficiency @ %def calculate_efficiency @ We define additional chunk, which will be used later on for inserting MPI/general MPI code. The code can is then removed by additional noweb filter if not compiled with the correspondig compiler flag. Overall variables, some additionally introduced due to the MPI parallelization and needed in sequentiell verison. <>= integer :: n_size integer :: n_dim_par logical :: box_success ! MPI-specific variables below @ Overall initialization. <>= call self%init_grid () if (reset_result) call self%reset_result () self%result%it_start = self%result%it_num cumulative_int = 0. cumulative_std = 0. n_size = 1 n_dim_par = floor (self%config%n_dim / 2.) @ Reset all last-iteration results before sampling. <>= self%result%it_num = self%result%it_start + it self%d = 0. self%box = 1 self%bin = 1 total_integral = 0. total_sq_integral = 0. total_variance = 0. sum_abs_f_pos = 0. max_abs_f_pos = 0. sum_abs_f_neg = 0. max_abs_f_neg = 0. box_success = .true. select type (rng) type is (rng_stream_t) call rng%next_substream () end select @ Pacify output by defining empty chunk (nothing to do here). <>= @ <>= @ Increment the box coordinates by 1. If we reach the largest value for the current axis (starting with the largest dimension number), we reset the counter to 1 and increment the next axis counter by 1. And so on, until we reach the maximum counter value of the axis with the lowest dimension, then we set [[success]] to false and the box coord is set to 1. <>= subroutine increment_box_coord (box, success) integer, dimension(:), intent(inout) :: box logical, intent(out) :: success integer :: j success = .true. do j = size (box), 1, -1 box(j) = box(j) + 1 if (box(j) <= self%config%n_boxes) return box(j) = 1 end do success = .false. end subroutine increment_box_coord @ %def increment_box_coord @ We parallelize [[VEGAS]] in simple forward manner. The hyper-cube is dissambled in to equidistant boxes in which we sample the integrand [[calls_per_box]] times. The workload of calculating those boxes is distributed along the worker. The number of dimensions which will be parallelised are $\lfloor \frac{D}{2} \rfloor$, such MPI Variables for [[vegas_integrate]]. We have to duplicate all buffers for [[MPI_Ireduce]], because we cannot use the same send or receive buffer. We temporarily store a (empty) grid, before communicating. <>= integer :: rank type(vegas_grid_t) :: grid @ MPI procedure-specific initialization. <>= call MPI_Comm_size (MPI_COMM_WORLD, n_size) call MPI_Comm_rank (MPI_COMM_WORLD, rank) @ Pre-sampling communication. We make a copy of the (actual) grid, which is unfilled when non-root. The actual grid is then broadcasted among the workers and inserted into the [[VEGAS]] object. <>= if (self%is_parallelizable ()) then grid = self%get_grid () call grid%broadcast () call self%set_grid (grid) end if @ Start index of the boxes for different ranks. If the random number generator is RngStream, we can advance the current stream in such a way, that we will getting matching numbers. Iff [[n_boxes]] is larger than 2, otherwise parallelization is useless. <>= if (self%is_parallelizable ()) then do k = 1, rank call increment_box_coord (self%box(1:n_dim_par), box_success) if (.not. box_success) exit end do select type (rng) type is (rng_stream_t) call rng%advance_state (self%config%n_dim * self%config%calls_per_box& & * self%config%n_boxes**(self%config%n_dim - n_dim_par) * rank) end select end if @ Increment [[n_size]] times the box coordinates. <>= if (self%is_parallelizable ()) then select type (rng) type is (rng_stream_t) call rng%advance_state (self%config%n_dim * self%config%calls_per_box& & * self%config%n_boxes**(self%config%n_dim - n_dim_par) * (n_size - 1)) end select end if @ Call to [[vegas_integrate_collect]]. <>= if (self%is_parallelizable ()) then call vegas_integrate_collect () if (rank /= 0) cycle iteration end if @ Reduce (in an non-blocking fashion) all sampled information via [[MPI_SUM]] or [[MPI_MAX]]. <>= subroutine vegas_integrate_collect () real(default) :: root_total_integral, root_total_sq_integral real(default) :: root_sum_abs_f_pos, root_max_abs_f_pos real(default) :: root_sum_abs_f_neg, root_max_abs_f_neg real(default), dimension(self%config%n_bins_max, self%config%n_dim) :: root_d type(MPI_Request), dimension(self%config%n_dim + 6) :: status root_d = 0._default root_sum_abs_f_pos = 0._default root_sum_abs_f_neg = 0._default root_max_abs_f_pos = 0._default root_sum_abs_f_neg = 0._default root_total_integral = 0._default root_total_sq_integral = 0._default call MPI_Ireduce (sum_abs_f_pos, root_sum_abs_f_pos, 1, MPI_DOUBLE_PRECISION,& & MPI_SUM, 0, MPI_COMM_WORLD, status(1)) call MPI_Ireduce (sum_abs_f_neg, root_sum_abs_f_neg, 1, MPI_DOUBLE_PRECISION,& & MPI_SUM, 0, MPI_COMM_WORLD, status(2)) call MPI_Ireduce (max_abs_f_pos, root_max_abs_f_pos, 1, MPI_DOUBLE_PRECISION,& & MPI_MAX, 0, MPI_COMM_WORLD, status(3)) call MPI_Ireduce (max_abs_f_neg, root_max_abs_f_neg, 1, MPI_DOUBLE_PRECISION,& & MPI_MAX, 0, MPI_COMM_WORLD, status(4)) call MPI_Ireduce (total_integral, root_total_integral, 1, MPI_DOUBLE_PRECISION,& & MPI_SUM, 0, MPI_COMM_WORLD, status(5)) call MPI_Ireduce (total_sq_integral, root_total_sq_integral, 1,& & MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, status(6)) do j = 1, self%config%n_dim call MPI_Ireduce (self%d(1:self%config%n_bins, j), root_d(1:self%config& &%n_bins, j), self%config%n_bins, MPI_DOUBLE_PRECISION, MPI_SUM, 0,& & MPI_COMM_WORLD, status(6 + j)) end do call MPI_Waitall (self%config%n_dim + 6, status, MPI_STATUSES_IGNORE) if (rank == 0) sum_abs_f_pos = root_sum_abs_f_pos if (rank == 0) sum_abs_f_neg = root_sum_abs_f_neg if (rank == 0) max_abs_f_pos = root_max_abs_f_pos if (rank == 0) max_abs_f_neg = root_max_abs_f_neg if (rank == 0) total_integral = root_total_integral if (rank == 0) total_sq_integral = root_total_sq_integral if (rank == 0) self%d = root_d end subroutine vegas_integrate_collect @ %def vegas_integrate_collect @ Obtain a random point inside the $n$-dimensional hypercube, transform onto the correct interval and calculate the bin volume. The additional factor [[n_bins]] is already applied to the [[jacobian]] (per dimension). <>= procedure, private :: random_point => vegas_random_point <>= subroutine vegas_random_point (self, rng, x, bin_volume) class(vegas_t), intent(inout) :: self class(rng_t), intent(inout) :: rng real(default), dimension(self%config%n_dim), intent(out) :: x real(default), intent(out) :: bin_volume integer :: j real(default) :: r, y, z, bin_width bin_volume = 1. ndim: do j = 1, self%config%n_dim call rng%generate (r) z = ((self%box(j) - 1 + r) / self%config%n_boxes) * self%config%n_bins + 1 self%bin(j) = max (min (int (z), self%config%n_bins), 1) if (self%bin(j) == 1) then bin_width = self%grid%xi(2, j) y = (z - self%bin(j)) * bin_width else bin_width = self%grid%xi(self%bin(j) + 1, j) - self%grid%xi(self%bin(j), j) y = self%grid%xi(self%bin(j), j) + (z - self%bin(j)) * bin_width end if x(j) = self%grid%x_lower(j) + y * self%grid%delta_x(j) bin_volume = bin_volume * bin_width end do ndim end subroutine vegas_random_point @ %def vegas_random_point @ Obtain a random point inside the $n$-dimensional hyper-cube. We neglect stratification and generate the random point in the most simple way. Hence, we do not need to know in which box we are actually sampling. This is useful for only for event generation. <>= procedure, private :: simple_random_point => vegas_simple_random_point <>= subroutine vegas_simple_random_point (self, rng, x, bin_volume) class(vegas_t), intent(inout) :: self class(rng_t), intent(inout) :: rng real(default), dimension(self%config%n_dim), intent(out) :: x real(default), intent(out) :: bin_volume integer :: j, k real(default) :: r, y, z, bin_width bin_volume = 1. ndim: do j = 1, self%config%n_dim call rng%generate (r) z = r * self%config%n_bins + 1 k = max (min (int (z), self%config%n_bins), 1) if (k == 1) then bin_width = self%grid%xi(2, j) y = (z - 1) * bin_width else bin_width = self%grid%xi(k + 1, j) - self%grid%xi(k, j) y = self%grid%xi(k, j) + (z - k) * bin_width end if x(j) = self%grid%x_lower(j) + y * self%grid%delta_x(j) bin_volume = bin_volume * bin_width end do ndim end subroutine vegas_simple_random_point @ %def vegas_simple_random_point @ <>= procedure, private :: accumulate_distribution => vegas_accumulate_distribution <>= subroutine vegas_accumulate_distribution (self, y) class(vegas_t), intent(inout) :: self real(default), intent(in) :: y integer :: j do j = 1, self%config%n_dim self%d(self%bin(j), j) = self%d(self%bin(j), j) + y end do end subroutine vegas_accumulate_distribution @ %def vegas_accumulate_distribution @ Generate weighted random event. The weight given by the overall jacobian \begin{equation} \Delta_{\text{jac}} = \prod_{j=1}^{d} \left( x_j^+ - x_j^- \right) \frac{N_{\text{bins}}^d}{N_{\text{calls}}} \end{equation} includes the overall non-changing factors $\frac{1}{N_{\text{calls}}}$-factor (divisions are expensive) and $N_{\text{bins}}^{d}$, the latter combined with [[bin_volume]] gives rise to the probability, see [[vegas_init_grid]] for details. We have to factor out $N_{\text{calls}}$ to retrieve the correct weight. <>= procedure :: generate_weighted => vegas_generate_weighted_event <>= subroutine vegas_generate_weighted_event (self, func, rng, x) class(vegas_t), intent(inout) :: self class(vegas_func_t), intent(inout) :: func class(rng_t), intent(inout) :: rng real(default), dimension(self%config%n_dim), intent(inout) :: x real(default) :: bin_volume call self%simple_random_point (rng, x, bin_volume) ! Cancel n_calls from jacobian with n_calls self%result%evt_weight = self%config%n_calls * self%jacobian * bin_volume & & * func%evaluate (x) end subroutine vegas_generate_weighted_event @ %def vegas_generate_weighted_event @ Generate random event. We accept on the rate \begin{equation} \frac{|f(x)|}{\underset{x}{\max} |f(x)|}. \end{equation} We keep separate maximum weights for positive and negative weights, and use them, accordingly. In the case of unweighted event generation, if the current weight exceeds the the maximum weight, we update the maximum, accordingly. <>= procedure, public :: generate_unweighted=> vegas_generate_unweighted_event <>= subroutine vegas_generate_unweighted_event (self, func, rng, x) class(vegas_t), intent(inout) :: self class(vegas_func_t), intent(inout) :: func class(rng_t), intent(inout) :: rng real(default), dimension(self%config%n_dim), intent(out) :: x real(default) :: bin_volume real(default) :: max_abs_f real(default) :: r associate (result => self%result) generate: do call self%generate_weighted (func, rng, x) max_abs_f = merge (result%max_abs_f_pos, result%max_abs_f_neg, & & result%evt_weight > 0.) if (result%evt_weight > max_abs_f) then result%evt_weight_excess = & & result%evt_weight / max_abs_f - 1._default exit generate end if call rng%generate (r) ! Do not use division, because max_abs_f could be zero. if (max_abs_f * r <= abs(result%evt_weight)) then exit generate end if end do generate end associate end subroutine vegas_generate_unweighted_event @ %def vegas_random_event \section{I/0 operation} \label{sec:i0-operation} @ Write grid to file. We use the original VAMP formater. <>= character(len=*), parameter, private :: & descr_fmt = "(1X,A)", & integer_fmt = "(1X,A18,1X,I15)", & integer_array_fmt = "(1X,I18,1X,I15)", & logical_fmt = "(1X,A18,1X,L1)", & - double_fmt = "(1X,A18,1X,E16.8E4)", & - double_array_fmt = "(1X,I18,1X,E16.8E4)", & - double_array2_fmt = "(1X,2(1X,I8),1X,E16.8E4)" + double_fmt = "(1X,A18,1X,E24.16E4)", & + double_array_fmt = "(1X,I18,1X,E24.16E4)", & + double_array_pac_fmt = "(1X,I18,1X,E16.8E4)", & + double_array2_fmt = "(1X,2(1X,I8),1X,E24.16E4)", & + double_array2_pac_fmt = "(1X,2(1X,I8),1X,E16.8E4)" @ %def descr_fmt integer_fmt integer_array_fmt logical_fmt @ %def double_fmt double_array_fmt double_array2_fmt <>= procedure, public :: write_grid => vegas_write_grid <>= subroutine vegas_write_grid (self, unit) class(vegas_t), intent(in) :: self integer, intent(in), optional :: unit integer :: u integer :: i, j u = given_output_unit (unit) write (u, descr_fmt) "begin type(vegas_t)" write (u, integer_fmt) "n_dim =", self%config%n_dim write (u, integer_fmt) "n_bins_max =", self%config%n_bins_max write (u, double_fmt) "alpha =", self%config%alpha write (u, integer_fmt) "iterations =", self%config%iterations write (u, integer_fmt) "mode =", self%config%mode write (u, integer_fmt) "calls_per_box =", self%config%calls_per_box write (u, integer_fmt) "n_calls =", self%config%n_calls write (u, integer_fmt) "n_calls_min =", self%config%n_calls_min write (u, integer_fmt) "n_boxes =", self%config%n_boxes write (u, integer_fmt) "n_bins =", self%config%n_bins write (u, integer_fmt) "it_start =", self%result%it_start write (u, integer_fmt) "it_num =", self%result%it_num write (u, integer_fmt) "samples =", self%result%samples write (u, double_fmt) "sum_int_wgtd =", self%result%sum_int_wgtd write (u, double_fmt) "sum_wgts =", self%result%sum_wgts write (u, double_fmt) "sum_chi =", self%result%sum_chi write (u, double_fmt) "chi2 =", self%result%chi2 write (u, double_fmt) "efficiency =", self%result%efficiency write (u, double_fmt) "efficiency =", self%result%efficiency_pos write (u, double_fmt) "efficiency =", self%result%efficiency_neg write (u, double_fmt) "max_abs_f =", self%result%max_abs_f write (u, double_fmt) "max_abs_f_pos =", self%result%max_abs_f_pos write (u, double_fmt) "max_abs_f_neg =", self%result%max_abs_f_neg write (u, double_fmt) "result =", self%result%result write (u, double_fmt) "std =", self%result%std write (u, double_fmt) "hypercube_volume =", self%hypercube_volume write (u, double_fmt) "jacobian =", self%jacobian write (u, descr_fmt) "begin x_lower" do j = 1, self%config%n_dim write (u, double_array_fmt) j, self%grid%x_lower(j) end do write (u, descr_fmt) "end x_lower" write (u, descr_fmt) "begin x_upper" do j = 1, self%config%n_dim write (u, double_array_fmt) j, self%grid%x_upper(j) end do write (u, descr_fmt) "end x_upper" write (u, descr_fmt) "begin delta_x" do j = 1, self%config%n_dim write (u, double_array_fmt) j, self%grid%delta_x(j) end do write (u, descr_fmt) "end delta_x" write (u, integer_fmt) "n_bins =", self%config%n_bins write (u, descr_fmt) "begin bin" do j = 1, self%config%n_dim write (u, integer_array_fmt) j, self%bin(j) end do write (u, descr_fmt) "end n_bin" write (u, integer_fmt) "n_boxes =", self%config%n_boxes write (u, descr_fmt) "begin box" do j = 1, self%config%n_dim write (u, integer_array_fmt) j, self%box(j) end do write (u, descr_fmt) "end box" write (u, descr_fmt) "begin d" do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max write (u, double_array2_fmt) i, j, self%d(i, j) end do end do write (u, descr_fmt) "end d" write (u, descr_fmt) "begin xi" do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max + 1 write (u, double_array2_fmt) i, j, self%grid%xi(i, j) end do end do write (u, descr_fmt) "end xi" write (u, descr_fmt) "end type(vegas_t)" end subroutine vegas_write_grid @ %def vegas_write_grid @ Read grid configuration from file. <>= procedure, public :: read_grid => vegas_read_grid <>= subroutine vegas_read_grid (self, unit) class(vegas_t), intent(out) :: self integer, intent(in) :: unit integer :: i, j character(len=80) :: buffer integer :: ibuffer, jbuffer read (unit, descr_fmt) buffer read (unit, integer_fmt) buffer, ibuffer read (unit, integer_fmt) buffer, jbuffer select type(self) type is (vegas_t) self = vegas_t (n_dim = ibuffer, n_bins_max = jbuffer) end select read (unit, double_fmt) buffer, self%config%alpha read (unit, integer_fmt) buffer, self%config%iterations read (unit, integer_fmt) buffer, self%config%mode read (unit, integer_fmt) buffer, self%config%calls_per_box read (unit, integer_fmt) buffer, self%config%n_calls read (unit, integer_fmt) buffer, self%config%n_calls_min read (unit, integer_fmt) buffer, self%config%n_boxes read (unit, integer_fmt) buffer, self%config%n_bins self%grid%n_bins = self%config%n_bins read (unit, integer_fmt) buffer, self%result%it_start read (unit, integer_fmt) buffer, self%result%it_num read (unit, integer_fmt) buffer, self%result%samples read (unit, double_fmt) buffer, self%result%sum_int_wgtd read (unit, double_fmt) buffer, self%result%sum_wgts read (unit, double_fmt) buffer, self%result%sum_chi read (unit, double_fmt) buffer, self%result%chi2 read (unit, double_fmt) buffer, self%result%efficiency read (unit, double_fmt) buffer, self%result%efficiency_pos read (unit, double_fmt) buffer, self%result%efficiency_neg read (unit, double_fmt) buffer, self%result%max_abs_f read (unit, double_fmt) buffer, self%result%max_abs_f_pos read (unit, double_fmt) buffer, self%result%max_abs_f_neg read (unit, double_fmt) buffer, self%result%result read (unit, double_fmt) buffer, self%result%std read (unit, double_fmt) buffer, self%hypercube_volume read (unit, double_fmt) buffer, self%jacobian read (unit, descr_fmt) buffer do j = 1, self%config%n_dim read (unit, double_array_fmt) jbuffer, self%grid%x_lower(j) end do read (unit, descr_fmt) buffer read (unit, descr_fmt) buffer do j = 1, self%config%n_dim read (unit, double_array_fmt) jbuffer, self%grid%x_upper(j) end do read (unit, descr_fmt) buffer read (unit, descr_fmt) buffer do j = 1, self%config%n_dim read (unit, double_array_fmt) jbuffer, self%grid%delta_x(j) end do read (unit, descr_fmt) buffer read (unit, integer_fmt) buffer, self%config%n_bins read (unit, descr_fmt) buffer do j = 1, self%config%n_dim read (unit, integer_array_fmt) jbuffer, self%bin(j) end do read (unit, descr_fmt) buffer read (unit, integer_fmt) buffer, self%config%n_boxes read (unit, descr_fmt) buffer do j = 1, self%config%n_dim read (unit, integer_array_fmt) jbuffer, self%box(j) end do read (unit, descr_fmt) buffer read (unit, descr_fmt) buffer do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max read (unit, double_array2_fmt) ibuffer, jbuffer, self%d(i, j) end do end do read (unit, descr_fmt) buffer read (unit, descr_fmt) buffer do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max + 1 read (unit, double_array2_fmt) ibuffer, jbuffer, self%grid%xi(i, j) end do end do read (unit, descr_fmt) buffer read (unit, descr_fmt) buffer end subroutine vegas_read_grid @ %def vegas_read_grid Read and write a grid from an unformatted file. <>= procedure :: write_binary_grid => vegas_write_binary_grid procedure :: read_binary_grid => vegas_read_binary_grid <>= subroutine vegas_write_binary_grid (self, unit) class(vegas_t), intent(in) :: self integer, intent(in) :: unit integer :: i, j write (unit) self%config%n_dim write (unit) self%config%n_bins_max write (unit) self%config%alpha write (unit) self%config%iterations write (unit) self%config%mode write (unit) self%config%calls_per_box write (unit) self%config%n_calls write (unit) self%config%n_calls_min write (unit) self%config%n_boxes write (unit) self%config%n_bins write (unit) self%result%it_start write (unit) self%result%it_num write (unit) self%result%samples write (unit) self%result%sum_int_wgtd write (unit) self%result%sum_wgts write (unit) self%result%sum_chi write (unit) self%result%chi2 write (unit) self%result%efficiency write (unit) self%result%efficiency_pos write (unit) self%result%efficiency_neg write (unit) self%result%max_abs_f write (unit) self%result%max_abs_f_pos write (unit) self%result%max_abs_f_neg write (unit) self%result%result write (unit) self%result%std write (unit) self%hypercube_volume write (unit) self%jacobian do j = 1, self%config%n_dim write (unit) j, self%grid%x_lower(j) end do do j = 1, self%config%n_dim write (unit) j, self%grid%x_upper(j) end do do j = 1, self%config%n_dim write (unit) j, self%grid%delta_x(j) end do write (unit) self%config%n_bins do j = 1, self%config%n_dim write (unit) j, self%bin(j) end do write (unit) self%config%n_boxes do j = 1, self%config%n_dim write (unit) j, self%box(j) end do do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max write (unit) i, j, self%d(i, j) end do end do do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max + 1 write (unit) i, j, self%grid%xi(i, j) end do end do end subroutine vegas_write_binary_grid subroutine vegas_read_binary_grid (self, unit) class(vegas_t), intent(out) :: self integer, intent(in) :: unit integer :: i, j integer :: ibuffer, jbuffer read (unit) ibuffer read (unit) jbuffer select type(self) type is (vegas_t) self = vegas_t (n_dim = ibuffer, n_bins_max = jbuffer) end select read (unit) self%config%alpha read (unit) self%config%iterations read (unit) self%config%mode read (unit) self%config%calls_per_box read (unit) self%config%n_calls read (unit) self%config%n_calls_min read (unit) self%config%n_boxes read (unit) self%config%n_bins self%grid%n_bins = self%config%n_bins read (unit) self%result%it_start read (unit) self%result%it_num read (unit) self%result%samples read (unit) self%result%sum_int_wgtd read (unit) self%result%sum_wgts read (unit) self%result%sum_chi read (unit) self%result%chi2 read (unit) self%result%efficiency read (unit) self%result%efficiency_pos read (unit) self%result%efficiency_neg read (unit) self%result%max_abs_f read (unit) self%result%max_abs_f_pos read (unit) self%result%max_abs_f_neg read (unit) self%result%result read (unit) self%result%std read (unit) self%hypercube_volume read (unit) self%jacobian do j = 1, self%config%n_dim read (unit) jbuffer, self%grid%x_lower(j) end do do j = 1, self%config%n_dim read (unit) jbuffer, self%grid%x_upper(j) end do do j = 1, self%config%n_dim read (unit) jbuffer, self%grid%delta_x(j) end do read (unit) self%config%n_bins do j = 1, self%config%n_dim read (unit) jbuffer, self%bin(j) end do read (unit) self%config%n_boxes do j = 1, self%config%n_dim read (unit) jbuffer, self%box(j) end do do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max read (unit) ibuffer, jbuffer, self%d(i, j) end do end do do j = 1, self%config%n_dim do i = 1, self%config%n_bins_max + 1 read (unit) ibuffer, jbuffer, self%grid%xi(i, j) end do end do end subroutine vegas_read_binary_grid @ %def vegas_write_binary_grid, vegas_read_binary_grid \section{Unit tests} \label{sec:unit-tests} Test module, followed by the corresponding implementation module. <<[[vegas_ut.f90]]>>= <> module vegas_ut use unit_tests use vegas_uti <> <> contains <> end module vegas_ut @ %def vegas_ut @ <<[[vegas_uti.f90]]>>= <> module vegas_uti <> use io_units use constants, only: pi use format_defs, only: FMT_10, FMT_12 use rng_base use rng_stream use vegas <> <> <> contains <> end module vegas_uti @ %def vegas_uti @ API: driver for the unit tests below. <>= public :: vegas_test <>= subroutine vegas_test (u, results) integer, intent(in) :: u type(test_results_t), intent(inout) :: results <> end subroutine vegas_test @ %def vegas_test @ \subsubsection{Test function} \label{sec:test-function} We use the example from the Monte Carlo Examples of the GSL library \begin{equation} I = \int_{-pi}^{+pi} {dk_x/(2 pi)} \int_{-pi}^{+pi} {dk_y/(2 pi)} \int_{-pi}^{+pi} {dk_z/(2 pi)} 1 / (1 - cos(k_x)cos(k_y)cos(k_z)). \end{equation} The integral is reduced to region (0,0,0) $\rightarrow$ ($\pi$, $\pi$, $\pi$) and multiplied by 8. <>= type, extends (vegas_func_t) :: vegas_test_func_t ! contains <> end type vegas_test_func_t @ %def vegas_test_func_t @ Evaluate the integrand. <>= procedure, public :: evaluate => vegas_test_func_evaluate <>= real(default) function vegas_test_func_evaluate (self, x) result (f) class(vegas_test_func_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x f = 1.0 / (pi**3) f = f / ( 1.0 - cos (x(1)) * cos (x(2)) * cos (x(3))) end function vegas_test_func_evaluate @ %def vegas_test_func_evaluate @ The second test function is the normalised n-dim.\@ gaussian distribution. <>= type, extends (vegas_func_t) :: vegas_gaussian_test_func_t ! contains <> end type vegas_gaussian_test_func_t @ %def vegas_gaussian_test_func_t @ Evaluate the integrand. <>= procedure, public :: evaluate => vegas_gaussian_evaluate <>= real(default) function vegas_gaussian_evaluate (self, x) result (f) class(vegas_gaussian_test_func_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x real(default), parameter :: inv_sqrt_pi = 1._default / sqrt(pi) f = inv_sqrt_pi**size (x) f = f * exp (- dot_product(x, x)) end function vegas_gaussian_evaluate @ %def vegas_gaussian_evaluate @ The third test function is a three-dimensional polynomial function which factories. The function is defined in such a way that the integral in the unit range is normalised to zero. \begin{equation} f(x) = - \frac{8}{3} (x + 1)*(y-1)*z \end{equation} <>= type, extends (vegas_func_t) :: vegas_polynomial_func_t ! contains <> end type vegas_polynomial_func_t @ %def vegas_polynomial_func_t <>= procedure, public :: evaluate => vegas_polynomial_evaluate <>= real(default) function vegas_polynomial_evaluate (self, x) result (f) class(vegas_polynomial_func_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x f = - 8. / 3. * (x(1) + 1.) * (x(2) - 1.) * x(3) end function vegas_polynomial_evaluate @ %def vegas_polynomial_evaluate @ \subsubsection{MC Integrator check} \label{sec:mc-integrator-check} Initialise the VEGAS MC integrator and call to [[vegas_init_grid]] for the initialisation of the grid. <>= call test (vegas_1, "vegas_1", "VEGAS initialisation and& & grid preparation", u, results) <>= public :: vegas_1 <>= subroutine vegas_1 (u) integer, intent(in) :: u type(vegas_t) :: mc_integrator class(rng_t), allocatable :: rng class(vegas_func_t), allocatable :: func real(default), dimension(3), parameter :: x_lower = 0., & x_upper = pi real(default) :: result, abserr write (u, "(A)") "* Test output: vegas_1" write (u, "(A)") "* Purpose: initialise the VEGAS MC integrator and the grid" write (u, "(A)") write (u, "(A)") "* Initialise random number generator (default seed)" write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_dim = 3" write (u, "(A)") allocate (vegas_test_func_t :: func) mc_integrator = vegas_t (3) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 10000" write (u, "(A)") call mc_integrator%set_limits (x_lower, x_upper) call mc_integrator%set_calls (10000) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 10000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 2000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (2000) call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () call rng%final () deallocate (rng) end subroutine vegas_1 @ %def vegas_1 @ \subsubsection{Configuration and result check} \label{sec:conf-result-check} Initialise the MC integrator. Get and write the config object, also the (empty) result object. <>= call test (vegas_2, "vegas_2", "VEGAS configuration and result object", u, results) <>= public :: vegas_2 <>= subroutine vegas_2 (u) integer, intent(in) :: u type(vegas_t) :: mc_integrator type(vegas_config_t) :: mc_integrator_config type(vegas_result_t) :: mc_integrator_result write (u, "(A)") "* Test output: vegas_2" write (u, "(A)") "* Purpose: use transparent containers for& & configuration and result." write (u, "(A)") write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_dim = 10" write (u, "(A)") mc_integrator = vegas_t (10) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 10000 (Importance Sampling)" write (u, "(A)") call mc_integrator%set_calls (10000) write (u, "(A)") write (u, "(A)") "* Get VEGAS config object and write out" write (u, "(A)") call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(A)") write (u, "(A)") "* Get VEGAS empty result object and write out" write (u, "(A)") mc_integrator_result = mc_integrator%get_result () call mc_integrator_result%write (u) write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () end subroutine vegas_2 @ %def vegas_2 @ \subsubsection{Grid check} \label{sec:conf-result-check} Initialise the MC integrator. Get and write the config object. Integrate the gaussian distribution. Get and write the result object. Before and after integration get the grid object and output both. Repeat with different number of dimensions. <>= call test (vegas_3, "vegas_3", "VEGAS integration of multi-dimensional gaussian", u, results) <>= public :: vegas_3 <>= subroutine vegas_3 (u) integer, intent(in) :: u type(vegas_t) :: mc_integrator class(rng_t), allocatable :: rng class(vegas_func_t), allocatable :: func real(default), dimension(3), parameter :: x_lower_3 = -10._default, & x_upper_3 = 10._default type(vegas_config_t) :: mc_integrator_config type(vegas_grid_t) :: mc_integrator_grid type(vegas_result_t) :: mc_integrator_result real(default) :: result, abserr write (u, "(A)") "* Test output: vegas_3" write (u, "(A)") "* Purpose: Integrate gaussian distribution." write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_dim = 3" write (u, "(A)") allocate (vegas_gaussian_test_func_t :: func) mc_integrator = vegas_t (3) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 10000" write (u, "(A)") call mc_integrator%set_limits (x_lower_3, x_upper_3) call mc_integrator%set_calls (10000) write (u, "(A)") write (u, "(A)") "* Get VEGAS config object and write out" write (u, "(A)") call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(A)") write (u, "(A)") "* Get VEGAS grid object and write out" write (u, "(A)") mc_integrator_grid = mc_integrator%get_grid () - call mc_integrator_grid%write (u) + call mc_integrator_grid%write (u, pacify = .true.) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 20000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 2000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (2000) call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Get VEGAS result object and write out" write (u, "(A)") mc_integrator_result = mc_integrator%get_result () call mc_integrator_result%write (u) write (u, "(A)") write (u, "(A)") "* Get VEGAS grid object and write out" write (u, "(A)") mc_integrator_grid = mc_integrator%get_grid () - call mc_integrator_grid%write (u) + call mc_integrator_grid%write (u, pacify = .true.) write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () end subroutine vegas_3 @ %def vegas_3 \subsubsection{Three-dimensional integration with polynomial function} \label{sec:conf-result-check} Initialise the MC integrator. Get and write the config object. Integrate the factorisable polynomial function. Get and write the result object. Repeat with different number of dimensions. <>= call test (vegas_4, "vegas_4", "VEGAS integration of three& &-dimensional factorisable polynomial function", u, results) <>= public :: vegas_4 <>= subroutine vegas_4 (u) integer, intent(in) :: u type(vegas_t) :: mc_integrator class(rng_t), allocatable :: rng class(vegas_func_t), allocatable :: func real(default), dimension(3), parameter :: x_lower_3 = 0._default, & x_upper_3 = 1._default type(vegas_config_t) :: mc_integrator_config type(vegas_result_t) :: mc_integrator_result real(default) :: result, abserr write (u, "(A)") "* Test output: vegas_4" write (u, "(A)") "* Purpose: Integrate gaussian distribution." write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_dim = 3" write (u, "(A)") allocate (vegas_polynomial_func_t :: func) mc_integrator = vegas_t (3) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 2000" write (u, "(A)") call mc_integrator%set_limits (x_lower_3, x_upper_3) call mc_integrator%set_calls (2000) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 2000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(A)") write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 20000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (20000) call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(A)") write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () end subroutine vegas_4 @ %def vegas_4 @ \subsubsection{Event generation} Initialise the MC integrator. Integrate the gaussian distribution. Get and write the result object. Finally, generate events in accordance to the adapted grid and print them out. <>= call test (vegas_5, "vegas_5", "VEGAS integration and event& & generation of multi-dimensional gaussian", u, results) <>= public :: vegas_5 <>= subroutine vegas_5 (u) integer, intent(in) :: u type(vegas_t) :: mc_integrator class(rng_t), allocatable :: rng class(vegas_func_t), allocatable :: func real(default), dimension(1), parameter :: x_lower_1 = -10._default, & x_upper_1 = 10._default type(vegas_config_t) :: mc_integrator_config type(vegas_result_t) :: mc_integrator_result integer :: i, u_event real(default), dimension(1) :: event, mean, delta, M2 real(default) :: result, abserr write (u, "(A)") "* Test output: vegas_5" write (u, "(A)") "* Purpose: Integrate gaussian distribution." write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_dim = 1" write (u, "(A)") allocate (vegas_gaussian_test_func_t :: func) mc_integrator = vegas_t (1) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 20000" write (u, "(A)") call mc_integrator%set_limits (x_lower_1, x_upper_1) call mc_integrator%set_calls (20000) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, opt_verbose=.true., result=result, abserr=abserr) call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") & & "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 2000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (2000) call mc_integrator%integrate (func, rng, 3, opt_verbose=.true., result=result, abserr=abserr) call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") & & "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Generate 10000 events based on the adaptation and& & calculate mean and variance" write (u, "(A)") mean = 0._default M2 = 0._default do i = 1, 10000 call mc_integrator%generate_unweighted (func, rng, event) delta = event - mean mean = mean + delta / i M2 = M2 + delta * (event - mean) end do write (u, "(2X,A)") "Result:" write (u, "(4X,A," // FMT_12 //")") & & "mean = ", mean write (u, "(4X,A," // FMT_12 //")") & & "(sample) std. dev. = ", sqrt (M2 / (9999)) write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () end subroutine vegas_5 @ %def vegas_5 @ \subsubsection{Grid I/O} \label{sec:grid-io} Initialise the MC integrator. Get and write the config object. Integrate the factorisable polynomial function. Get and write the result object. Write grid to file and start with fresh grid. <>= call test (vegas_6, "vegas_6", "VEGAS integrate and write grid, & & read grid and continue", u, results) <>= public :: vegas_6 <>= subroutine vegas_6 (u) integer, intent(in) :: u type(vegas_t) :: mc_integrator class(rng_t), allocatable :: rng class(vegas_func_t), allocatable :: func real(default), dimension(3), parameter :: x_lower_3 = 0._default, & x_upper_3 = 1._default type(vegas_config_t) :: mc_integrator_config type(vegas_result_t) :: mc_integrator_result real(default) :: result, abserr integer :: unit write (u, "(A)") "* Test output: vegas_6" write (u, "(A)") "* Purpose: Write and read grid, and continue." write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_dim = 3" write (u, "(A)") allocate (vegas_polynomial_func_t :: func) mc_integrator = vegas_t (3) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 2000" write (u, "(A)") call mc_integrator%set_limits (x_lower_3, x_upper_3) call mc_integrator%set_calls (2000) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 2000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(A)") write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Write grid to file vegas_io.grid" write (u, "(A)") unit = free_unit () open (unit, file = "vegas_io.grid", & action = "write", status = "replace") call mc_integrator%write_grid (unit) close (unit) write (u, "(A)") write (u, "(A)") "* Read grid from file vegas_io.grid" write (u, "(A)") call mc_integrator%final () open (unit, file = "vegas_io.grid", & action = "read", status = "old") call mc_integrator%read_grid (unit) close (unit) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 20000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (20000) call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) call mc_integrator%get_config (mc_integrator_config) call mc_integrator_config%write (u) write (u, "(A)") write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () end subroutine vegas_6 @ %def vegas_6 @ \clearpage %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{VAMP2} \label{sec:vamp2} We concentrate all configuration and run-time data in a derived-type, such that, [[mci_t]] can spwan each time a distinctive MCI VEGAS integrator object. <<[[vamp2.f90]]>>= <> module vamp2 <> <> use io_units use format_utils, only: pac_fmt use format_utils, only: write_separator, write_indent use format_defs, only: FMT_17 use diagnostics use rng_base use rng_stream, only: rng_stream_t use vegas <> <> <> <> <> <> contains <> end module vamp2 @ %def vamp2 <>= @ <>= use mpi_f08 !NODEP! @ \subsection{Type: vamp2\_func\_t} \label{sec:vamp2-func} We extend [[vegas_func_t]] with the multi-channel weights and the [[vegas_grid_t]], such that, the overall multi-channel weight can be calculated by the function itself. We add an additional logicial [[valid_x]], if it is set to [[.false.]], we do not compute weighted function and just set the weighted integrand to zero. This behavior is in particular very useful, if a mapping is prohibited or fails. Or in the case of WHIZARD, a phase cut is applied. <>= public :: vamp2_func_t <>= type, abstract, extends(vegas_func_t) :: vamp2_func_t integer :: current_channel = 0 integer :: n_dim = 0 integer :: n_channel = 0 integer :: n_calls = 0 logical :: valid_x = .false. real(default), dimension(:, :), allocatable :: xi real(default), dimension(:), allocatable :: det real(default), dimension(:), allocatable :: wi real(default), dimension(:), allocatable :: gi type(vegas_grid_t), dimension(:), allocatable :: grids real(default) :: g = 0. contains <> end type vamp2_func_t @ %def vamp2_func_t @ Init. <>= procedure, public :: init => vamp2_func_init <>= subroutine vamp2_func_init (self, n_dim, n_channel) class(vamp2_func_t), intent(out) :: self integer, intent(in) :: n_dim integer, intent(in) :: n_channel self%n_dim = n_dim self%n_channel = n_channel allocate (self%xi(n_dim, n_channel), source=0._default) allocate (self%det(n_channel), source=1._default) allocate (self%wi(n_channel), source=0._default) allocate (self%gi(n_channel), source=0._default) allocate (self%grids(n_channel)) end subroutine vamp2_func_init @ %def vamp2_func_init @ Set current channel. <>= procedure, public :: set_channel => vamp2_func_set_channel <>= subroutine vamp2_func_set_channel (self, channel) class(vamp2_func_t), intent(inout) :: self integer, intent(in) :: channel self%current_channel = channel end subroutine vamp2_func_set_channel @ %def vamp2_func_set_channel @ Get number of function calls for which $f \neq 0$. <>= procedure, public :: get_n_calls => vamp2_func_get_n_calls <>= integer function vamp2_func_get_n_calls (self) result (n_calls) class(vamp2_func_t), intent(in) :: self n_calls = self%n_calls end function vamp2_func_get_n_calls @ %def vamp2_func_get_func_calls @ Reset number of calls. <>= procedure, public :: reset_n_calls => vamp2_func_reset_n_calls <>= subroutine vamp2_func_reset_n_calls (self) class(vamp2_func_t), intent(inout) :: self self%n_calls = 0 end subroutine vamp2_func_reset_n_calls @ %def vamp2_func_reset_n_calls @ Evaluate mappings. We defer this method to be implemented by the user. The result must be written to [[xi]] and [[det]]. The mapping is defined by $\phi : U \rightarrow M$. We map $x \in M$ to the different mappings of the hypercube $U_{i}$, such that $x_{i} \in U_{i}$. The mapping should determine, whether [[x]] is a valid point, e.g. can be mapped, or is restricted otherwise. <>= procedure(vamp2_func_evaluate_maps), deferred :: evaluate_maps <>= abstract interface subroutine vamp2_func_evaluate_maps (self, x) import :: vamp2_func_t, default class(vamp2_func_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x end subroutine vamp2_func_evaluate_maps end interface @ %def vamp2_evaluate_func @ Evaluate channel weights. The calling procedure must handle the case of a vanishing overall probability density where either a channel weight or a channel probability vanishes. <>= procedure, private :: evaluate_weight => vamp2_func_evaluate_weight <>= subroutine vamp2_func_evaluate_weight (self) class(vamp2_func_t), intent(inout) :: self integer :: ch self%g = 0. self%gi = 0. !$OMP PARALLEL DO PRIVATE(ch) SHARED(self) do ch = 1, self%n_channel if (self%wi(ch) /= 0) then self%gi(ch) = self%grids(ch)%get_probability (self%xi(:, ch)) end if end do !$OMP END PARALLEL DO if (self%gi(self%current_channel) /= 0) then do ch = 1, self%n_channel if (self%wi(ch) /= 0 .and. self%det(ch) /= 0) then self%g = self%g + self%wi(ch) * self%gi(ch) / self%det(ch) end if end do self%g = self%g / self%gi(self%current_channel) end if end subroutine vamp2_func_evaluate_weight @ %def vamp2_func_evaluate_weight @ Evaluate function at [[x]]. We call this procedure in [[vamp2_func_evaluate]]. <>= procedure(vamp2_func_evaluate_func), deferred :: evaluate_func <>= abstract interface real(default) function vamp2_func_evaluate_func (self, x) result (f) import :: vamp2_func_t, default class(vamp2_func_t), intent(in) :: self real(default), dimension(:), intent(in) :: x end function vamp2_func_evaluate_func end interface @ %def vamp2_func_evaluate_func <>= procedure, public :: evaluate => vamp2_func_evaluate <>= real(default) function vamp2_func_evaluate (self, x) result (f) class(vamp2_func_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x call self%evaluate_maps (x) f = 0. self%gi = 0. self%g = 1 if (self%valid_x) then call self%evaluate_weight () if (self%g /= 0) then f = self%evaluate_func (x) / self%g self%n_calls = self%n_calls + 1 end if end if end function vamp2_func_evaluate @ %def vamp2_func_evaluate \subsection{Type: vamp2\_config\_t} \label{sec:vamp2-config} This is a transparent container which incorporates and extends the definitions in [[vegas_config]]. The parent object can then be used to parametrise the VEGAS grids directly, where the new parameters are exclusively used in the multi-channel implementation of VAMP2. [[n_calls_min]] is calculated by [[n_calls_min_per_channel]] and [[n_channel]]. The channels weights (and the result [[n_calls]] for each channel) are calculated regarding [[n_calls_threshold]]. <>= public :: vamp2_config_t <>= type, extends(vegas_config_t) :: vamp2_config_t integer :: n_channel = 0 integer :: n_calls_min_per_channel = 20 integer :: n_calls_threshold = 10 integer :: n_chains = 0 logical :: stratified = .true. logical :: equivalences = .false. real(default) :: beta = 0.5_default real(default) :: accuracy_goal = 0._default real(default) :: error_goal = 0._default real(default) :: rel_error_goal = 0._default contains <> end type vamp2_config_t @ %def vamp2_config_t @ Write. <>= procedure, public :: write => vamp2_config_write <>= subroutine vamp2_config_write (self, unit, indent) class(vamp2_config_t), intent(in) :: self integer, intent(in), optional :: unit integer, intent(in), optional :: indent integer :: u, ind u = given_output_unit (unit) ind = 0; if (present (indent)) ind = indent call self%vegas_config_t%write (unit, indent) call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Number of channels = ", self%n_channel call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Min. number of calls per channel (setting calls) = ", & & self%n_calls_min_per_channel call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Threshold number of calls (adapting weights) = ", & & self%n_calls_threshold call write_indent (u, ind) write (u, "(2x,A,I0)") & & "Number of chains = ", self%n_chains call write_indent (u, ind) write (u, "(2x,A,L1)") & & "Stratified = ", self%stratified call write_indent (u, ind) write (u, "(2x,A,L1)") & & "Equivalences = ", self%equivalences call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "Adaption power (beta) = ", self%beta if (self%accuracy_goal > 0) then call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "accuracy_goal = ", self%accuracy_goal end if if (self%error_goal > 0) then call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "error_goal = ", self%error_goal end if if (self%rel_error_goal > 0) then call write_indent (u, ind) write (u, "(2x,A," // FMT_17 // ")") & & "rel_error_goal = ", self%rel_error_goal end if end subroutine vamp2_config_write @ %def vamp2_config_write @ \subsection{Type: vamp2\_result\_t} \label{sec:vamp2-result} This is a transparent container which incorporates and extends the definitions of [[vegas_result_t]]. <>= public :: vamp2_result_t <>= type, extends(vegas_result_t) :: vamp2_result_t contains <> end type vamp2_result_t @ %def vamp2_result_t @ Output. <>= procedure, public :: write => vamp2_result_write <>= subroutine vamp2_result_write (self, unit, indent) class(vamp2_result_t), intent(in) :: self integer, intent(in), optional :: unit integer, intent(in), optional :: indent integer :: u, ind u = given_output_unit (unit) ind = 0; if (present (indent)) ind = indent call self%vegas_result_t%write (unit, indent) end subroutine vamp2_result_write @ %def vamp2_result_write @ \subsection{Type: vamp2\_equivalences\_t} \label{sec:vamp2-eqv} <>= integer, parameter, public :: & VEQ_IDENTITY = 0, VEQ_INVERT = 1, VEQ_SYMMETRIC = 2, VEQ_INVARIANT = 3 @ @ Channel equivalences. Store retrieving and sourcing channel. <>= type :: vamp2_equi_t integer :: ch integer :: ch_src integer, dimension(:), allocatable :: perm integer, dimension(:), allocatable :: mode contains <> end type vamp2_equi_t @ %def vamp2_equi_t @ Write equivalence. <>= procedure :: write => vamp2_equi_write <>= subroutine vamp2_equi_write (self, unit, indent) class(vamp2_equi_t), intent(in) :: self integer, intent(in), optional :: unit integer, intent(in), optional :: indent integer :: u, ind u = given_output_unit (unit) ind = 0; if (present (indent)) ind = indent call write_indent (u, ind) write (u, "(2(A,1X,I0))") "src:", self%ch_src, "-> dest:", self%ch call write_indent (u, ind) write (u, "(A,99(1X,I0))") "Perm: ", self%perm call write_indent (u, ind) write (u, "(A,99(1X,I0))") "Mode: ", self%mode end subroutine vamp2_equi_write @ %def vamp2_equi_write @ <>= public :: vamp2_equivalences_t <>= type :: vamp2_equivalences_t private integer :: n_eqv = 0 integer :: n_channel = 0 integer :: n_dim = 0 type(vamp2_equi_t), dimension(:), allocatable :: eqv integer, dimension(:), allocatable :: map integer, dimension(:), allocatable :: multiplicity integer, dimension(:), allocatable :: symmetry logical, dimension(:), allocatable :: independent integer, dimension(:), allocatable :: equivalent_to_ch logical, dimension(:, :), allocatable :: dim_is_invariant contains <> end type vamp2_equivalences_t @ %def vamp2_equivalences_t @ Constructor. <>= interface vamp2_equivalences_t module procedure vamp2_equivalences_init end interface vamp2_equivalences_t <>= type(vamp2_equivalences_t) function vamp2_equivalences_init (& n_eqv, n_channel, n_dim) result (eqv) integer, intent(in) :: n_eqv, n_channel, n_dim eqv%n_eqv = n_eqv eqv%n_channel = n_channel eqv%n_dim = n_dim allocate (eqv%eqv(n_eqv)) allocate (eqv%map(n_channel), source = 0) allocate (eqv%multiplicity(n_channel), source = 0) allocate (eqv%symmetry(n_channel), source = 0) allocate (eqv%independent(n_channel), source = .true.) allocate (eqv%equivalent_to_ch(n_channel), source = 0) allocate (eqv%dim_is_invariant(n_dim, n_channel), source = .false.) end function vamp2_equivalences_init @ %def vamp2_equivlences_init @ Write equivalences. <>= procedure :: write => vamp2_equivalences_write <>= subroutine vamp2_equivalences_write (self, unit, indent) class(vamp2_equivalences_t), intent(in) :: self integer, intent(in), optional :: unit integer, intent(in), optional :: indent integer :: u, ind, i_eqv, ch u = given_output_unit (unit) ind = 0; if (present (indent)) ind = indent write (u, "(A)") "Inequivalent channels:" if (allocated (self%independent)) then do ch = 1, self%n_channel if (self%independent(ch)) then write (u, "(2X,A,1x,I0,A,4x,A,I0,4x,A,I0,4x,A,999(L1))") & "Channel", ch, ":", & "Mult. = ", self%multiplicity(ch), & "Symm. = ", self%symmetry(ch), & "Invar.: ", self%dim_is_invariant(:, ch) end if end do else write (u, "(A)") "[not allocated]" end if write (u, "(A)") "Equivalence list:" if (allocated (self%eqv)) then do i_eqv = 1, self%n_eqv write (u, "(2X,A,1X,I0)") "i_eqv:", i_eqv call self%eqv(i_eqv)%write (unit, indent = ind + 4) end do else write (u, "(A)") "[not allocated]" end if end subroutine vamp2_equivalences_write @ %def vamp2_equivalences_write @ Is allocated. <>= procedure, public :: is_allocated => vamp2_equivalences_is_allocated <>= logical function vamp2_equivalences_is_allocated (self) result (yorn) class(vamp2_equivalences_t), intent(in) :: self yorn = allocated (self%eqv) end function vamp2_equivalences_is_allocated @ %def vamp2_equivalences_is_allocated @ Get source channel and destination channel for given equivalence. <>= procedure, public :: get_channels => vamp2_equivalences_get_channels <>= subroutine vamp2_equivalences_get_channels (eqv, i_eqv, dest, src) class(vamp2_equivalences_t), intent(in) :: eqv integer, intent(in) :: i_eqv integer, intent(out) :: dest, src dest = eqv%eqv(i_eqv)%ch src = eqv%eqv(i_eqv)%ch_src end subroutine vamp2_equivalences_get_channels @ %def vamp2_equivalences_get_channels @ <>= procedure, public :: get_mode => vamp2_equivalences_get_mode procedure, public :: get_perm => vamp2_equivalences_get_perm <>= function vamp2_equivalences_get_mode (eqv, i_eqv) result (mode) class(vamp2_equivalences_t), intent(in) :: eqv integer, intent(in) :: i_eqv integer, dimension(:), allocatable :: mode mode = eqv%eqv(i_eqv)%mode end function vamp2_equivalences_get_mode function vamp2_equivalences_get_perm (eqv, i_eqv) result (perm) class(vamp2_equivalences_t), intent(in) :: eqv integer, intent(in) :: i_eqv integer, dimension(:), allocatable :: perm perm = eqv%eqv(i_eqv)%perm end function vamp2_equivalences_get_perm @ %def vamp2_equivalences_get_perm, vamp2_equivalences_get_mode @ <>= procedure, public :: set_equivalence => vamp2_equivalences_set_equivalence <>= subroutine vamp2_equivalences_set_equivalence & (eqv, i_eqv, dest, src, perm, mode) class(vamp2_equivalences_t), intent(inout) :: eqv integer, intent(in) :: i_eqv integer, intent(in) :: dest, src integer, dimension(:), intent(in) :: perm, mode integer :: i if (dest < 1 .or. dest > eqv%n_channel) call msg_bug & ("VAMP2: set_equivalences: destination channel out of range.") if (src < 1 .or. src > eqv%n_channel) call msg_bug & ("VAMP2: set_equivalences: source channel out of range.") if (size(perm) /= eqv%n_dim) call msg_bug & ("VAMP2: set_equivalences: size(perm) does not match n_dim.") if (size(mode) /= eqv%n_dim) call msg_bug & ("VAMP2: set_equivalences: size(mode) does not match n_dim.") eqv%eqv(i_eqv)%ch = dest eqv%eqv(i_eqv)%ch_src = src allocate (eqv%eqv(i_eqv)%perm (size (perm))) do i = 1, size (perm) eqv%eqv(i_eqv)%perm(i) = perm(i) end do allocate (eqv%eqv(i_eqv)%mode (size (mode))) do i = 1, size (mode) eqv%eqv(i_eqv)%mode(i) = mode(i) end do end subroutine vamp2_equivalences_set_equivalence @ %def vamp2_equivalences_set_equivalence @ Freeze equivalences. <>= procedure, public :: freeze => vamp2_equivalences_freeze <>= subroutine vamp2_equivalences_freeze (self) class(vamp2_equivalences_t), intent(inout) :: self integer :: i_eqv, ch, upper, lower ch = 0 do i_eqv = 1, self%n_eqv if (ch /= self%eqv(i_eqv)%ch) then ch = self%eqv(i_eqv)%ch self%map(ch) = i_eqv end if end do do ch = 1, self%n_channel lower = self%map(ch) if (ch == self%n_channel) then upper = self%n_eqv else upper = self%map(ch + 1) - 1 end if associate (eqv => self%eqv, n_eqv => size (self%eqv(lower:upper))) if (.not. all(eqv(lower:upper)%ch == ch) .or. & eqv(lower)%ch_src > ch) then do i_eqv = lower, upper call self%eqv(i_eqv)%write () end do call msg_bug ("VAMP2: vamp2_equivalences_freeze: & &equivalence order is not correct.") end if self%symmetry(ch) = count (eqv(lower:upper)%ch_src == ch) if (mod (n_eqv, self%symmetry(ch)) /= 0) then do i_eqv = lower, upper call self%eqv(i_eqv)%write () end do call msg_bug ("VAMP2: vamp2_equivalences_freeze: & &permutation count is not correct.") end if self%multiplicity(ch) = n_eqv / self%symmetry(ch) self%independent(ch) = all (eqv(lower:upper)%ch_src >= ch) self%equivalent_to_ch(ch) = eqv(lower)%ch_src self%dim_is_invariant(:, ch) = eqv(lower)%mode == VEQ_INVARIANT end associate end do end subroutine vamp2_equivalences_freeze @ %def vamp2_equivalences_freeze @ \subsection{Type: vamp2\_t} \label{sec:vamp2-t} <>= public :: vamp2_t <>= type :: vamp2_t private type(vamp2_config_t) :: config type(vegas_t), dimension(:), allocatable :: integrator integer, dimension(:), allocatable :: chain real(default), dimension(:), allocatable :: weight real(default), dimension(:), allocatable :: integral real(default), dimension(:), allocatable :: variance real(default), dimension(:), allocatable :: efficiency type(vamp2_result_t) :: result type(vamp2_equivalences_t) :: equivalences logical :: event_prepared real(default), dimension(:), allocatable :: event_weight contains <> end type vamp2_t <>= interface vamp2_t module procedure vamp2_init end interface vamp2_t @ %def vamp2_t @ Constructor. <>= type(vamp2_t) function vamp2_init (n_channel, n_dim, alpha, beta, n_bins_max,& & n_calls_min_per_channel, iterations, mode) result (self) integer, intent(in) :: n_channel integer, intent(in) :: n_dim integer, intent(in), optional :: n_bins_max integer, intent(in), optional :: n_calls_min_per_channel real(default), intent(in), optional :: alpha real(default), intent(in), optional :: beta integer, intent(in), optional :: iterations integer, intent(in), optional :: mode integer :: ch self%config%n_dim = n_dim self%config%n_channel = n_channel if (present (n_bins_max)) self%config%n_bins_max = n_bins_max if (present (n_calls_min_per_channel)) self%config%n_calls_min_per_channel = n_calls_min_per_channel if (present (alpha)) self%config%alpha = alpha if (present (beta)) self%config%beta = beta if (present (iterations)) self%config%iterations = iterations if (present (mode)) self%config%mode = mode allocate (self%chain(n_channel), source=0) allocate (self%integrator(n_channel)) allocate (self%weight(n_channel), source=0._default) do ch = 1, n_channel self%integrator(ch) = vegas_t (n_dim, alpha, n_bins_max, 1, mode) end do self%weight = 1. / self%config%n_channel call self%reset_result () allocate (self%event_weight(self%config%n_channel), source = 0._default) self%event_prepared = .false. end function vamp2_init @ %def vamp2_init <>= procedure, public :: final => vamp2_final <>= subroutine vamp2_final (self) class(vamp2_t), intent(inout) :: self integer :: ch do ch = 1, self%config%n_channel call self%integrator(ch)%final () end do end subroutine vamp2_final @ %def vamp2_final @ Output. <>= procedure, public :: write => vamp2_write <>= subroutine vamp2_write (self, unit, indent) class(vamp2_t), intent(in) :: self integer, intent(in), optional :: unit integer, intent(in), optional :: indent integer :: u, ind, ch u = given_output_unit (unit) ind = 0; if (present (indent)) ind = indent call write_indent (u, ind) write (u, "(A)") "VAMP2: VEGAS AMPlified 2" call write_indent (u, ind) call self%config%write (unit, indent) call self%result%write (unit, indent) end subroutine vamp2_write @ %def vamp2_write @ Get the config object. <>= procedure, public :: get_config => vamp2_get_config <>= subroutine vamp2_get_config (self, config) class(vamp2_t), intent(in) :: self type(vamp2_config_t), intent(out) :: config config = self%config end subroutine vamp2_get_config @ %def vamp2_get_config @ Set non-runtime dependent configuration. It will no be possible to change [[n_bins_max]]. <>= procedure, public :: set_config => vamp2_set_config <>= subroutine vamp2_set_config (self, config) class(vamp2_t), intent(inout) :: self class(vamp2_config_t), intent(in) :: config integer :: ch self%config%equivalences = config%equivalences self%config%n_calls_min_per_channel = config%n_calls_min_per_channel self%config%n_calls_threshold = config%n_calls_threshold self%config%n_calls_min = config%n_calls_min self%config%beta = config%beta self%config%accuracy_goal = config%accuracy_goal self%config%error_goal = config%error_goal self%config%rel_error_goal = config%rel_error_goal do ch = 1, self%config%n_channel call self%integrator(ch)%set_config (config) end do end subroutine vamp2_set_config @ %def vamp2_set_config @ Set the overall number of calls. The number of calls each channel is scaled by the channel weights \begin{equation} N_i = \alpha_i N. \end{equation} <>= procedure, public :: set_calls => vamp2_set_n_calls <>= subroutine vamp2_set_n_calls (self, n_calls) class(vamp2_t), intent(inout) :: self integer, intent(in) :: n_calls integer :: ch self%config%n_calls_min = self%config%n_calls_min_per_channel & & * self%config%n_channel self%config%n_calls = max(n_calls, self%config%n_calls_min) if (self%config%n_calls > n_calls) then write (msg_buffer, "(A,I0)") "VAMP2: [set_calls] number of calls too few,& & reset to = ", self%config%n_calls call msg_message () end if do ch = 1, self%config%n_channel call self%integrator(ch)%set_calls (max (nint (self%config%n_calls *& & self%weight(ch)), self%config%n_calls_min_per_channel)) end do end subroutine vamp2_set_n_calls @ %def vamp2_set_n_calls @ Set limits. We only support same limits for all channels. <>= procedure, public :: set_limits => vamp2_set_limits <>= subroutine vamp2_set_limits (self, x_upper, x_lower) class(vamp2_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x_upper real(default), dimension(:), intent(in) :: x_lower integer :: ch do ch = 1, self%config%n_channel call self%integrator(ch)%set_limits (x_upper, x_lower) end do end subroutine vamp2_set_limits @ %def vamp2_set_limits @ Set [[n_chains]] and the (actual) chains. [[chain]] must have size [[n_channels]] and each elements must store an index to a corresponding chain. This means, that channels with equal index correspond to the same chain, and we refer to those as chained weights, where we average the contributions of the chained weights in [[vamp2_adapt_weights]]. <>= procedure, public :: set_chain => vamp2_set_chain <>= subroutine vamp2_set_chain (self, n_chains, chain) class(vamp2_t), intent(inout) :: self integer, intent(in) :: n_chains integer, dimension(:), intent(in) :: chain if (size (chain) /= self%config%n_channel) then call msg_bug ("VAMP2: set chain: size of chain array does not match n_channel.") else call msg_message ("VAMP2: set chain: use chained weights.") end if self%config%n_chains = n_chains self%chain = chain end subroutine vamp2_set_chain @ %def vamp2_set_chain @ Set channel equivalences. <>= procedure, public :: set_equivalences => vamp2_set_equivalences <>= subroutine vamp2_set_equivalences (self, equivalences) class(vamp2_t), intent(inout) :: self type(vamp2_equivalences_t), intent(in) :: equivalences self%equivalences = equivalences end subroutine vamp2_set_equivalences @ %def vamp2_set_equivalences @ Get [[n_calls]] calculated by [[VEGAS]]. <>= procedure, public :: get_n_calls => vamp2_get_n_calls <>= elemental real(default) function vamp2_get_n_calls (self) result (n_calls) class(vamp2_t), intent(in) :: self n_calls = sum (self%integrator%get_calls ()) end function vamp2_get_n_calls @ %def vamp2_get_n_calls @ Get the cumulative result of the integration. Recalculate weighted average of the integration. <>= procedure, public :: get_integral => vamp2_get_integral <>= elemental real(default) function vamp2_get_integral (self) result (integral) class(vamp2_t), intent(in) :: self integral = 0. if (self%result%sum_wgts > 0.) then integral = self%result%sum_int_wgtd / self%result%sum_wgts end if end function vamp2_get_integral @ %def vamp2_get_integral @ Get the cumulative variance of the integration. Recalculate the variance. <>= procedure, public :: get_variance => vamp2_get_variance <>= elemental real(default) function vamp2_get_variance (self) result (variance) class(vamp2_t), intent(in) :: self variance = 0. if (self%result%sum_wgts > 0.) then variance = 1.0 / self%result%sum_wgts end if end function vamp2_get_variance @ %def vamp2_get_variance @ Get efficiency. <>= procedure, public :: get_efficiency => vamp2_get_efficiency <>= elemental real(default) function vamp2_get_efficiency (self) result (efficiency) class(vamp2_t), intent(in) :: self efficiency = 0. if (self%result%efficiency > 0.) then efficiency = self%result%efficiency end if end function vamp2_get_efficiency @ %def vamp2_get_efficiency @ Get event weight and event weight excess. <>= procedure :: get_evt_weight => vamp2_get_evt_weight procedure :: get_evt_weight_excess => vamp2_get_evt_weight_excess <>= real(default) function vamp2_get_evt_weight (self) result (evt_weight) class(vamp2_t), intent(in) :: self evt_weight = self%result%evt_weight end function vamp2_get_evt_weight real(default) function vamp2_get_evt_weight_excess (self) result (evt_weight_excess) class(vamp2_t), intent(in) :: self evt_weight_excess = self%result%evt_weight_excess end function vamp2_get_evt_weight_excess @ %def vamp2_get_evt_weight, vamp2_get_evt_weight_excess @ Get procedure to retrieve channel-th grid. <>= procedure :: get_grid => vamp2_get_grid <>= type(vegas_grid_t) function vamp2_get_grid (self, channel) result (grid) class(vamp2_t), intent(in) :: self integer, intent(in) :: channel if (channel < 1 .or. channel > self%config%n_channel) & call msg_bug ("VAMP2: vamp2_get_grid: channel index < 1 or > n_channel.") grid = self%integrator(channel)%get_grid () end function vamp2_get_grid @ %def vamp2_get_grid @ Adapt. We adapt the weights due the contribution of variances with $\beta > 0$. \begin{equation} \alpha_i = \frac{\alpha_i V_i^\beta}{\sum_i \alpha_i V_i^\beta} \end{equation} If [[n_calls_threshold]] is set, we rescale the weights in such a way, that the [[n_calls]] for each channel are greater than [[n_calls_threshold]]. We calculate the distance of the weights to the [[weight_min]] and reset those weights which are less than [[weight_mins]] to this value. The other values are accordingly resized to fit the boundary condition of the partition of unity. <>= procedure, private :: adapt_weights => vamp2_adapt_weights <>= subroutine vamp2_adapt_weights (self) class(vamp2_t), intent(inout) :: self integer :: n_weights_underflow real(default) :: weight_min, sum_weights_underflow self%weight = self%weight * self%integrator%get_variance ()**self%config%beta if (sum (self%weight) == 0) self%weight = real(self%config%n_calls, default) if (self%config%n_chains > 0) then call chain_weights () end if self%weight = self%weight / sum(self%weight) if (self%config%n_calls_threshold /= 0) then weight_min = real(self%config%n_calls_threshold, default) & & / self%config%n_calls sum_weights_underflow = sum (self%weight, self%weight < weight_min) n_weights_underflow = count (self%weight < weight_min) where (self%weight < weight_min) self%weight = weight_min elsewhere self%weight = self%weight * (1. - n_weights_underflow * weight_min) & & / (1. - sum_weights_underflow) end where end if call self%set_calls (self%config%n_calls) contains <> end subroutine vamp2_adapt_weights @ %def vamp2_adapt_weights @ We average the weights over their respective chain members. <>= subroutine chain_weights () integer :: ch real(default) :: average do ch = 1, self%config%n_chains average = max (sum (self%weight, self%chain == ch), 0._default) if (average /= 0) then average = average / count (self%chain == ch) where (self%chain == ch) self%weight = average end where end if end do end subroutine chain_weights @ %def chain_weights <>= procedure, private :: apply_equivalences => vamp2_apply_equivalences <>= subroutine vamp2_apply_equivalences (self) class(vamp2_t), intent(inout) :: self integer :: ch, ch_src, j, j_src, i_eqv real(default), dimension(:, :, :), allocatable :: d real(default), dimension(:, :), allocatable :: d_src integer, dimension(:), allocatable :: mode, perm if (.not. self%equivalences%is_allocated ()) then call msg_bug ("VAMP2: vamp2_apply_equivalences: & &cannot apply not-allocated equivalences.") end if allocate (d(self%config%n_bins_max, self%config%n_dim, & self%config%n_channel), source=0._default) associate (eqv => self%equivalences, nb => self%config%n_bins_max) do i_eqv = 1, self%equivalences%n_eqv call eqv%get_channels (i_eqv, ch, ch_src) d_src = self%integrator(ch_src)%get_distribution () mode = eqv%get_mode (i_eqv) perm = eqv%get_perm (i_eqv) do j = 1, self%config%n_dim select case (mode (j)) case (VEQ_IDENTITY) d(:, j, ch) = d(:, j, ch) + & d_src(:, perm(j)) case (VEQ_INVERT) d(:, j, ch) = d(:, j, ch) + & d_src(nb:1:-1, perm(j)) case (VEQ_SYMMETRIC) d(:, j, ch) = d(:, j, ch) + & d_src(:, perm(j)) / 2. + & d_src(nb:1:-1, perm(j)) / 2. case (VEQ_INVARIANT) d(:, j, ch) = 1._default end select end do end do end associate do ch = 1, self%config%n_channel call self%integrator(ch)%set_distribution (d(:, :, ch)) end do end subroutine vamp2_apply_equivalences @ %def vamp2_apply_equivalences @ Reset the cumulative result. <>= procedure, public :: reset_result => vamp2_reset_result <>= subroutine vamp2_reset_result (self) class(vamp2_t), intent(inout) :: self self%result%sum_int_wgtd = 0. self%result%sum_wgts = 0. self%result%sum_chi = 0. self%result%it_num = 0 self%result%samples = 0 self%result%chi2 = 0 self%result%efficiency = 0. end subroutine vamp2_reset_result @ %def vamp2_reset_result @ Integrate. We integrate each channel separately and combine the results \begin{align} I & = \sum_i \alpha_i I_i, \\ \sigma^2 & = \sum_i \alpha_i^2 \sigma^2_i. \end{align} Although, the (population) variance is given by \begin{equation} \begin{split} \sigma^2 & = \frac{1}{N} \left( \sum_i \alpha_i I^2_i - I^2 \right) \\ & = \frac{1}{N - 1} \left( \sum_i \left( N_i \sigma^2_i + I^2_i \right) -I^2 \right) \\ & = \frac{1}{N - 1} \left( \sum_i \alpha_i \sigma^2_i + \alpha_i I^2_i - I^2 \right), \end{split} \end{equation} where we used $\sigma^2_i = \frac{1}{N} \left( \langle I^2_i \rangle - \langle I_i \rangle^2 \right)$, we use the approximation for numeric stability. The population variance relates to sample variance \begin{equation} s^2 = \frac{n}{n - 1} \sigma^2, \end{equation} which gives an unbiased error estimate. Beside those adaption to multichannel, the overall processing of [[total_integral]], [[total_sq_integral]] and [[total_variance]] is the same as in [[vegas_integrate]]. <>= procedure, public :: integrate => vamp2_integrate <>= subroutine vamp2_integrate (self, func, rng, iterations, opt_reset_result,& & opt_refine_grid, opt_adapt_weight, opt_verbose, result, abserr) class(vamp2_t), intent(inout) :: self class(vamp2_func_t), intent(inout) :: func class(rng_t), intent(inout) :: rng integer, intent(in), optional :: iterations logical, intent(in), optional :: opt_reset_result logical, intent(in), optional :: opt_refine_grid logical, intent(in), optional :: opt_adapt_weight logical, intent(in), optional :: opt_verbose real(default), optional, intent(out) :: result, abserr integer :: it, ch real(default) :: total_integral, total_sq_integral, total_variance, chi, wgt real(default) :: cumulative_int, cumulative_std logical :: reset_result = .true. logical :: adapt_weight = .true. logical :: refine_grid = .true. logical :: verbose = .false. <> if (present (iterations)) self%config%iterations = iterations if (present (opt_reset_result)) reset_result = opt_reset_result if (present (opt_adapt_weight)) adapt_weight = opt_adapt_weight if (present (opt_refine_grid)) refine_grid = opt_refine_grid if (present (opt_verbose)) verbose = opt_verbose <> if (verbose) then call msg_message ("Results: [it, calls, integral, error, chi^2, eff.]") end if iteration: do it = 1, self%config%iterations <> do ch = 1, self%config%n_channel func%wi(ch) = self%weight(ch) func%grids(ch) = self%integrator(ch)%get_grid () end do channel: do ch = 1, self%config%n_channel <> call func%set_channel (ch) call self%integrator(ch)%integrate ( & & func, rng, iterations, opt_refine_grid = .false., opt_verbose = verbose) end do channel <> total_integral = dot_product (self%weight, self%integrator%get_integral ()) total_sq_integral = dot_product (self%weight, self%integrator%get_integral ()**2) total_variance = self%config%n_calls * dot_product (self%weight**2, self%integrator%get_variance ()) associate (result => self%result) ! a**2 - b**2 = (a - b) * (a + b) total_variance = sqrt (total_variance + total_sq_integral) total_variance = 1. / self%config%n_calls * & & (total_variance + total_integral) * (total_variance - total_integral) ! Ensure variance is always positive and larger than zero if (total_variance < tiny (1._default) / epsilon (1._default) * max (total_integral**2, 1._default)) then total_variance = tiny (1._default) / epsilon (1._default) * max (total_integral**2, 1._default) end if wgt = 1. / total_variance result%result = total_integral result%std = sqrt (total_variance) result%samples = result%samples + 1 if (result%samples == 1) then result%chi2 = 0._default else chi = total_integral if (result%sum_wgts > 0) chi = chi - result%sum_int_wgtd / result%sum_wgts result%chi2 = result%chi2 * (result%samples - 2.0_default) result%chi2 = (wgt / (1._default + (wgt / result%sum_wgts))) & & * chi**2 result%chi2 = result%chi2 / (result%samples - 1._default) end if result%sum_wgts = result%sum_wgts + wgt result%sum_int_wgtd = result%sum_int_wgtd + (total_integral * wgt) result%sum_chi = result%sum_chi + (total_sq_integral * wgt) cumulative_int = result%sum_int_wgtd / result%sum_wgts cumulative_std = sqrt (1. / result%sum_wgts) call calculate_efficiency () if (verbose) then - write (msg_buffer, "(I0,1x,I0,1x, 4(E16.8E4,1x))") & + write (msg_buffer, "(I0,1x,I0,1x, 4(E24.16E4,1x))") & & it, self%config%n_calls, cumulative_int, cumulative_std, & & self%result%chi2, self%result%efficiency call msg_message () end if end associate if (adapt_weight) then call self%adapt_weights () end if if (refine_grid) then if (self%config%equivalences .and. self%equivalences%is_allocated ()) then call self%apply_equivalences () end if do ch = 1, self%config%n_channel call self%integrator(ch)%refine () end do end if end do iteration if (present (result)) result = cumulative_int if (present (abserr)) abserr = abs (cumulative_std) <> end subroutine vamp2_integrate @ %def vamp2_integrate @ <>= contains subroutine calculate_efficiency () self%result%max_abs_f = dot_product (self%weight, & & self%integrator%get_max_abs_f ()) self%result%max_abs_f_pos = dot_product (self%weight, & & self%integrator%get_max_abs_f_pos ()) self%result%max_abs_f_neg = dot_product (self%weight, & & self%integrator%get_max_abs_f_neg ()) self%result%efficiency = 0. if (self%result%max_abs_f > 0.) then self%result%efficiency = & & dot_product (self%weight * self%integrator%get_max_abs_f (), & & self%integrator%get_efficiency ()) / self%result%max_abs_f ! TODO pos. or. negative efficiency would be very nice. end if end subroutine calculate_efficiency @ %def calculate_efficiency @ We define additional chunks, which we use to insert parallel/MPI code. <>= @ <>= cumulative_int = 0. cumulative_std = 0. if (reset_result) call self%reset_result () @ <>= total_integral = 0._default total_sq_integral = 0._default total_variance = 0._default @ <>= @ <>= @ @ Distribute workers up in chunks of [[n_size]]. <>= integer function map_channel_to_worker (channel, n_size) result (worker) integer, intent(in) :: channel integer, intent(in) :: n_size worker = mod (channel, n_size) end function map_channel_to_worker @ %def map_channel_to_rank <>= type(vegas_grid_t) :: grid type(MPI_Request) :: status integer :: rank, n_size, worker @ <>= call MPI_Comm_rank (MPI_COMM_WORLD, rank) call MPI_Comm_size (MPI_COMM_WORLD, n_size) @ Broadcast all a-priori weights. After setting the weights, we have to update the number of calls in each channel. Afterwards, we can collect the number of channels, which are not parallelized by [[VEGAS]] itself, [[n_channel_non_parallel]]. <>= call MPI_Ibcast (self%weight, self%config%n_channel, MPI_DOUBLE_PRECISION, 0,& & MPI_COMM_WORLD, status) do ch = 1, self%config%n_channel grid = self%integrator(ch)%get_grid () call grid%broadcast () call self%integrator(ch)%set_grid (grid) end do call MPI_Wait (status, MPI_STATUS_IGNORE) call self%set_calls (self%config%n_calls) @ We check on the parallelization state of the current [[VEGAS]] integrator. If [[VEGAS]] can not be parallelized on lowest level, we map the current channel to a rank and calculate the channel on that rank. On all other worker we just enhance the random-generator (when supported), see [[vegas_integrate]] for the details on the random-generator handling. <>= if (.not. self%integrator(ch)%is_parallelizable ()) then worker = map_channel_to_worker (ch, n_size) if (rank /= worker) then select type (rng) type is (rng_stream_t) call rng%next_substream () end select cycle channel end if else call MPI_Barrier (MPI_COMM_WORLD) end if @ Collect results, the actual communication is done inside the different objects. <>= call vamp2_integrate_collect () <>= subroutine vamp2_integrate_collect () type(vegas_result_t) :: result integer :: root_n_calls integer :: worker do ch = 1, self%config%n_channel if (self%integrator(ch)%is_parallelizable ()) cycle worker = map_channel_to_worker (ch, n_size) result = self%integrator(ch)%get_result () if (rank == 0) then if (worker /= 0) then call result%receive (worker, ch) call self%integrator(ch)%receive_distribution (worker, ch) call self%integrator(ch)%set_result (result) end if else if (rank == worker) then call result%send (0, ch) call self%integrator(ch)%send_distribution (0, ch) end if end if end do select type (func) class is (vamp2_func_t) call MPI_reduce (func%n_calls, root_n_calls, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD) if (rank == 0) then func%n_calls = root_n_calls else call func%reset_n_calls () end if end select end subroutine vamp2_integrate_collect @ %def vegas_integrate_collect @ Skip results analyze if non-root, after waiting for all processes to reach the barrier. <>= call MPI_barrier (MPI_COMM_WORLD) if (rank /= 0) cycle iteration @ Generate event from multi-channel weight $w(x) = f(x) / g(x)$. We select a channel using the a-priori weights and $f_{i}^{\text{max}}$, to flatten possible unbalanced channel weight(s). An additional rescale factor [[opt_event_rescale]] is applied to [[f_max]], iff set. <>= procedure, public :: generate_weighted => vamp2_generate_weighted_event <>= subroutine vamp2_generate_weighted_event (& self, func, rng, x) class(vamp2_t), intent(inout) :: self class(vamp2_func_t), intent(inout) :: func class(rng_t), intent(inout) :: rng real(default), dimension(self%config%n_dim), intent(out) :: x integer :: ch, i real(default) :: r if (.not. self%event_prepared) then call prepare_event () end if call rng%generate (r) nchannel: do ch = 1, self%config%n_channel r = r - self%event_weight(ch) if (r <= 0._default) exit nchannel end do nchannel ch = min (ch, self%config%n_channel) call func%set_channel (ch) call self%integrator(ch)%generate_weighted (func, rng, x) ! Norm weight by f_max, hidden in event_weight(ch), else by 1 self%result%evt_weight = self%integrator(ch)%get_evt_weight () & * self%weight(ch) / self%event_weight(ch) contains <> end subroutine vamp2_generate_weighted_event @ %def vamp2_generate_weighted_event @ Generate unweighted events. After selecting a channel $ch$ by the acceptance $r$ \begin{equation*} r > \operatorname*{argmax}_{ch} \sum_{i = 1}^{ch} \alpha_i, \end{equation*} we try for an event from the previously selected channel. If the event is rejected, we also reject the selected channel. <>= procedure, public :: generate_unweighted => vamp2_generate_unweighted_event <>= subroutine vamp2_generate_unweighted_event ( & & self, func, rng, x, opt_event_rescale) class(vamp2_t), intent(inout) :: self class(vamp2_func_t), intent(inout) :: func class(rng_t), intent(inout) :: rng real(default), dimension(self%config%n_dim), intent(out) :: x real(default), intent(in), optional :: opt_event_rescale integer :: ch, i real(default) :: r, max_abs_f, event_rescale event_rescale = 1._default if (present (opt_event_rescale)) then event_rescale = opt_event_rescale end if if (.not. self%event_prepared) then call prepare_event () end if generate: do call rng%generate (r) nchannel: do ch = 1, self%config%n_channel r = r - self%event_weight(ch) if (r <= 0._default) exit nchannel end do nchannel ch = min (ch, self%config%n_channel) call func%set_channel (ch) call self%integrator(ch)%generate_weighted (func, rng, x) self%result%evt_weight = self%integrator(ch)%get_evt_weight () max_abs_f = merge ( & self%integrator(ch)%get_max_abs_f_pos (), & self%integrator(ch)%get_max_abs_f_neg (), & self%result%evt_weight > 0.) self%result%evt_weight_excess = 0._default if (self%result%evt_weight > max_abs_f) then self%result%evt_weight_excess = self%result%evt_weight / max_abs_f - 1._default exit generate end if call rng%generate (r) ! Do not use division, because max_abs_f could be zero. if (event_rescale * max_abs_f * r <= abs(self%result%evt_weight)) then exit generate end if end do generate contains <> end subroutine vamp2_generate_unweighted_event @ %def vamp2_generate_event Prepare event generation. We have to set the channel weights and the grids for the integrand's object. We use an ansatz proposed by T. Ohl in the original VAMP code where we do not have to accept on \begin{equation*} \frac{w_i(x)}{\operatorname*{max}_{i, x} w_i(x)}, \end{equation*} after we have selected a channel by the weights $\alpha_i$. But rather, we use a more efficient way where we rescale the channel weights $\alpha_i$ \begin{equation*} \alpha_i \rightarrow \alpha_i \frac{\operatorname*{max}_x w_i(x)}{\operatorname*{max}_{i, x} w_i(x)}. \end{equation*} The overall magic is to insert a "1" and to move the uneasy part into the channel selection, such that we can generate events likewise in the single channel mode. We generate an unweighted event by \begin{equation*} \frac{w_i(x)}{\operatorname*{max}_{x} w_i{x}}, \end{equation*} after we have selected a channel by the rescaled event channel weights. The overall normalization $\operatorname*{max}_{i, x}$ is not needed because we normalize the event channel weights to one and therefore the overall normalization cancels. <>= subroutine prepare_event () integer :: i self%event_prepared = .false. do i = 1, self%config%n_channel func%wi(i) = self%weight(i) func%grids(i) = self%integrator(i)%get_grid () end do if (any (self%integrator%get_max_abs_f () > 0)) then self%event_weight = self%weight * self%integrator%get_max_abs_f () else self%event_weight = self%weight end if self%event_weight = self%event_weight / sum (self%event_weight) self%event_prepared = .true. end subroutine prepare_event @ %def prepare_event @ Write grids to unit. <>= character(len=*), parameter, private :: & descr_fmt = "(1X,A)", & integer_fmt = "(1X,A18,1X,I15)", & integer_array_fmt = "(1X,I18,1X,I15)", & logical_fmt = "(1X,A18,1X,L1)", & - double_fmt = "(1X,A18,1X,E16.8E4)", & - double_array_fmt = "(1X,I18,1X,E16.8E4)", & - double_array2_fmt = "(1X,2(1X,I8),1X,E16.8E4)" + double_fmt = "(1X,A18,1X,E24.16E4)", & + double_array_fmt = "(1X,I18,1X,E24.16E4)", & + double_array_pac_fmt = "(1X,I18,1X,E16.8E4)", & + double_array2_fmt = "(1X,2(1X,I8),1X,E24.16E4)", & + double_array2_pac_fmt = "(1X,2(1X,I8),1X,E16.8E4)" @ %def descr_fmt integer_fmt integer_array_fmt logical_fmt @ %def double_fmt double_array_fmt double_array2_fmt <>= procedure, public :: write_grids => vamp2_write_grids <>= subroutine vamp2_write_grids (self, unit) class(vamp2_t), intent(in) :: self integer, intent(in), optional :: unit integer :: u integer :: ch u = given_output_unit (unit) write (u, descr_fmt) "begin type(vamp2_t)" write (u, integer_fmt) "n_channel =", self%config%n_channel write (u, integer_fmt) "n_dim =", self%config%n_dim write (u, integer_fmt) "n_calls_min_ch =", self%config%n_calls_min_per_channel write (u, integer_fmt) "n_calls_thres =", self%config%n_calls_threshold write (u, integer_fmt) "n_chains =", self%config%n_chains write (u, logical_fmt) "stratified =", self%config%stratified write (u, double_fmt) "alpha =", self%config%alpha write (u, double_fmt) "beta =", self%config%beta write (u, integer_fmt) "n_bins_max =", self%config%n_bins_max write (u, integer_fmt) "iterations =", self%config%iterations write (u, integer_fmt) "n_calls =", self%config%n_calls write (u, integer_fmt) "it_start =", self%result%it_start write (u, integer_fmt) "it_num =", self%result%it_num write (u, integer_fmt) "samples =", self%result%samples write (u, double_fmt) "sum_int_wgtd =", self%result%sum_int_wgtd write (u, double_fmt) "sum_wgts =", self%result%sum_wgts write (u, double_fmt) "sum_chi =", self%result%sum_chi write (u, double_fmt) "chi2 =", self%result%chi2 write (u, double_fmt) "efficiency =", self%result%efficiency write (u, double_fmt) "efficiency_pos =", self%result%efficiency_pos write (u, double_fmt) "efficiency_neg =", self%result%efficiency_neg write (u, double_fmt) "max_abs_f =", self%result%max_abs_f write (u, double_fmt) "max_abs_f_pos =", self%result%max_abs_f_pos write (u, double_fmt) "max_abs_f_neg =", self%result%max_abs_f_neg write (u, double_fmt) "result =", self%result%result write (u, double_fmt) "std =", self%result%std write (u, descr_fmt) "begin weight" do ch = 1, self%config%n_channel write (u, double_array_fmt) ch, self%weight(ch) end do write (u, descr_fmt) "end weight" if (self%config%n_chains > 0) then write (u, descr_fmt) "begin chain" do ch = 1, self%config%n_channel write (u, integer_array_fmt) ch, self%chain(ch) end do write (u, descr_fmt) "end chain" end if write (u, descr_fmt) "begin integrator" do ch = 1, self%config%n_channel call self%integrator(ch)%write_grid (unit) end do write (u, descr_fmt) "end integrator" write (u, descr_fmt) "end type(vamp2_t)" end subroutine vamp2_write_grids @ %def vamp2_write_grids @ Read grids from unit. <>= procedure, public :: read_grids => vamp2_read_grids <>= subroutine vamp2_read_grids (self, unit) class(vamp2_t), intent(out) :: self integer, intent(in), optional :: unit integer :: u integer :: ibuffer, jbuffer, ch character(len=80) :: buffer read (unit, descr_fmt) buffer read (unit, integer_fmt) buffer, ibuffer read (unit, integer_fmt) buffer, jbuffer select type (self) type is (vamp2_t) self = vamp2_t (n_channel = ibuffer, n_dim = jbuffer) end select read (unit, integer_fmt) buffer, self%config%n_calls_min_per_channel read (unit, integer_fmt) buffer, self%config%n_calls_threshold read (unit, integer_fmt) buffer, self%config%n_chains read (unit, logical_fmt) buffer, self%config%stratified read (unit, double_fmt) buffer, self%config%alpha read (unit, double_fmt) buffer, self%config%beta read (unit, integer_fmt) buffer, self%config%n_bins_max read (unit, integer_fmt) buffer, self%config%iterations read (unit, integer_fmt) buffer, self%config%n_calls read (unit, integer_fmt) buffer, self%result%it_start read (unit, integer_fmt) buffer, self%result%it_num read (unit, integer_fmt) buffer, self%result%samples read (unit, double_fmt) buffer, self%result%sum_int_wgtd read (unit, double_fmt) buffer, self%result%sum_wgts read (unit, double_fmt) buffer, self%result%sum_chi read (unit, double_fmt) buffer, self%result%chi2 read (unit, double_fmt) buffer, self%result%efficiency read (unit, double_fmt) buffer, self%result%efficiency_pos read (unit, double_fmt) buffer, self%result%efficiency_neg read (unit, double_fmt) buffer, self%result%max_abs_f read (unit, double_fmt) buffer, self%result%max_abs_f_pos read (unit, double_fmt) buffer, self%result%max_abs_f_neg read (unit, double_fmt) buffer, self%result%result read (unit, double_fmt) buffer, self%result%std read (unit, descr_fmt) buffer do ch = 1, self%config%n_channel read (unit, double_array_fmt) ibuffer, self%weight(ch) end do read (unit, descr_fmt) buffer if (self%config%n_chains > 0) then read (unit, descr_fmt) buffer do ch = 1, self%config%n_channel read (unit, integer_array_fmt) ibuffer, self%chain(ch) end do read (unit, descr_fmt) buffer end if read (unit, descr_fmt) buffer do ch = 1, self%config%n_channel call self%integrator(ch)%read_grid (unit) end do read (unit, descr_fmt) buffer read (unit, descr_fmt) buffer end subroutine vamp2_read_grids @ %def vamp2_read_grids @ Read and write grids from an unformatted file. <>= procedure :: write_binary_grids => vamp2_write_binary_grids procedure :: read_binary_grids => vamp2_read_binary_grids <>= subroutine vamp2_write_binary_grids (self, unit) class(vamp2_t), intent(in) :: self integer, intent(in) :: unit integer :: ch write (unit) write (unit) self%config%n_channel write (unit) self%config%n_dim write (unit) self%config%n_calls_min_per_channel write (unit) self%config%n_calls_threshold write (unit) self%config%n_chains write (unit) self%config%stratified write (unit) self%config%alpha write (unit) self%config%beta write (unit) self%config%n_bins_max write (unit) self%config%iterations write (unit) self%config%n_calls write (unit) self%result%it_start write (unit) self%result%it_num write (unit) self%result%samples write (unit) self%result%sum_int_wgtd write (unit) self%result%sum_wgts write (unit) self%result%sum_chi write (unit) self%result%chi2 write (unit) self%result%efficiency write (unit) self%result%efficiency_pos write (unit) self%result%efficiency_neg write (unit) self%result%max_abs_f write (unit) self%result%max_abs_f_pos write (unit) self%result%max_abs_f_neg write (unit) self%result%result write (unit) self%result%std do ch = 1, self%config%n_channel write (unit) ch, self%weight(ch) end do if (self%config%n_chains > 0) then do ch = 1, self%config%n_channel write (unit) ch, self%chain(ch) end do end if do ch = 1, self%config%n_channel call self%integrator(ch)%write_binary_grid (unit) end do end subroutine vamp2_write_binary_grids subroutine vamp2_read_binary_grids (self, unit) class(vamp2_t), intent(out) :: self integer, intent(in) :: unit integer :: ch, ibuffer, jbuffer read (unit) read (unit) ibuffer read (unit) jbuffer select type (self) type is (vamp2_t) self = vamp2_t (n_channel = ibuffer, n_dim = jbuffer) end select read (unit) self%config%n_calls_min_per_channel read (unit) self%config%n_calls_threshold read (unit) self%config%n_chains read (unit) self%config%stratified read (unit) self%config%alpha read (unit) self%config%beta read (unit) self%config%n_bins_max read (unit) self%config%iterations read (unit) self%config%n_calls read (unit) self%result%it_start read (unit) self%result%it_num read (unit) self%result%samples read (unit) self%result%sum_int_wgtd read (unit) self%result%sum_wgts read (unit) self%result%sum_chi read (unit) self%result%chi2 read (unit) self%result%efficiency read (unit) self%result%efficiency_pos read (unit) self%result%efficiency_neg read (unit) self%result%max_abs_f read (unit) self%result%max_abs_f_pos read (unit) self%result%max_abs_f_neg read (unit) self%result%result read (unit) self%result%std do ch = 1, self%config%n_channel read (unit) ibuffer, self%weight(ch) end do if (self%config%n_chains > 0) then do ch = 1, self%config%n_channel read (unit) ibuffer, self%chain(ch) end do end if do ch = 1, self%config%n_channel call self%integrator(ch)%read_binary_grid (unit) end do end subroutine vamp2_read_binary_grids @ %def vamp2_write_binary_grids, vamp2_read_binary_grids @ \section{Unit tests} \label{sec:unit-tests} Test module, followed by the corresponding implementation module. <<[[vamp2_ut.f90]]>>= <> module vamp2_ut use unit_tests use vamp2_uti <> <> contains <> end module vamp2_ut @ %def vamp2_ut @ <<[[vamp2_uti.f90]]>>= <> module vamp2_uti <> use io_units use constants, only: pi use numeric_utils, only: nearly_equal use format_defs, only: FMT_12 use rng_base use rng_stream use vegas, only: vegas_func_t, vegas_grid_t, operator(==) use vamp2 <> <> <> contains <> end module vamp2_uti @ %def vamp2_uti @ API: driver for the unit tests below. <>= public :: vamp2_test <>= subroutine vamp2_test (u, results) integer, intent(in) :: u type(test_results_t), intent(inout) :: results <> end subroutine vamp2_test @ %def vamp2_test @ \subsubsection{Test function} \label{sec:test-function} We use the example from the Monte Carlo Examples of the GSL library \begin{equation} I = \int_{-pi}^{+pi} {dk_x/(2 pi)} \int_{-pi}^{+pi} {dk_y/(2 pi)} \int_{-pi}^{+pi} {dk_z/(2 pi)} 1 / (1 - cos(k_x)cos(k_y)cos(k_z)). \end{equation} The integral is reduced to region (0,0,0) $\rightarrow$ ($\pi$, $\pi$, $\pi$) and multiplied by 8. <>= type, extends (vamp2_func_t) :: vamp2_test_func_t ! contains <> end type vamp2_test_func_t @ %def vegas_test_func_t @ <>= procedure, public :: evaluate_maps => vamp2_test_func_evaluate_maps <>= subroutine vamp2_test_func_evaluate_maps (self, x) class(vamp2_test_func_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x self%xi(:, 1) = x self%det(1) = 1 self%valid_x = .true. end subroutine vamp2_test_func_evaluate_maps @ %def vamp2_test_func_evaluate_maps @ Evaluate the integrand. <>= procedure, public :: evaluate_func => vamp2_test_func_evaluate <>= real(default) function vamp2_test_func_evaluate (self, x) result (f) class(vamp2_test_func_t), intent(in) :: self real(default), dimension(:), intent(in) :: x f = 1.0 / (pi**3) f = f / ( 1.0 - cos (x(1)) * cos (x(2)) * cos (x(3))) end function vamp2_test_func_evaluate @ %def vamp2_test_func_evaluate @ The second test function implements \begin{equation} f(\vec{x}) = 4 \sin^{2}(\pi x_{1})\sin^{2}(\pi x_{2}) + 2\sin^2(\pi v), \end{equation} where \begin{align} x = u^{v} & y = u^{1 - v} \\ u = xy & v = \frac{1}{2} \left( 1 + \frac{\log(x/y}{\log(xy)} \right). \end{align} The jacobian is $\frac{\partial (x, y)}{\partial (u, v)}$. <>= type, extends(vamp2_func_t) :: vamp2_test_func_2_t ! contains <> end type vamp2_test_func_2_t @ %def vamp2_test_func_2_t @ Evaluate maps. <>= procedure :: evaluate_maps => vamp2_test_func_2_evaluate_maps <>= subroutine vamp2_test_func_2_evaluate_maps (self, x) class(vamp2_test_func_2_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x select case (self%current_channel) case (1) self%xi(:, 1) = x self%xi(1, 2) = x(1) * x(2) self%xi(2, 2) = 0.5 * ( 1. + log(x(1) / x(2)) / log(x(1) * x(2))) case (2) self%xi(1, 1) = x(1)**x(2) self%xi(2, 1) = x(1)**(1. - x(2)) self%xi(:, 2) = x end select self%det(1) = 1. self%det(2) = abs (log(self%xi(1, 2))) self%valid_x = .true. end subroutine vamp2_test_func_2_evaluate_maps @ %def vamp2_test_func_2_evaluate_maps @ Evaluate func. <>= procedure :: evaluate_func => vamp2_test_func_2_evaluate_func <>= real(default) function vamp2_test_func_2_evaluate_func (self, x) result (f) class(vamp2_test_func_2_t), intent(in) :: self real(default), dimension(:), intent(in) :: x f = 4. * sin(pi * self%xi(1, 1))**2 * sin(pi * self%xi(2, 1))**2 & + 2. * sin(pi * self%xi(2, 2))**2 end function vamp2_test_func_2_evaluate_func @ %def vamp2_test_func_2_evaluate_func @ The third test function implements \begin{equation} f(\vec{x}) = 5 x_{1}^4 + 5 (1 - x_{1})^4, \end{equation} where \begin{equation} x_1 = u^{1 / 5} \quad \vee \quad x_1 = 1 - v^{1 / 5} \end{equation} The jacobians are $\frac{\partial x_1}{\partial u} = \frac{1}{5} u^{-\frac{4}{5}}$ and $\frac{\partial x_1}{\partial v} = \frac{1}{5} v^{-\frac{4}{5}}$. <>= type, extends(vamp2_func_t) :: vamp2_test_func_3_t ! contains <> end type vamp2_test_func_3_t @ %def vamp2_test_func_3_t @ Evaluate maps. <>= procedure :: evaluate_maps => vamp2_test_func_3_evaluate_maps <>= subroutine vamp2_test_func_3_evaluate_maps (self, x) class(vamp2_test_func_3_t), intent(inout) :: self real(default), dimension(:), intent(in) :: x real(default) :: u, v, xx select case (self%current_channel) case (1) u = x(1) xx = u**0.2_default v = (1 - xx)**5._default case (2) v = x(1) xx = 1 - v**0.2_default u = xx**5._default end select self%det(1) = 0.2_default * u**(-0.8_default) self%det(2) = 0.2_default * v**(-0.8_default) self%xi(:, 1) = [u] self%xi(:, 2) = [v] self%valid_x = .true. end subroutine vamp2_test_func_3_evaluate_maps @ %def vamp2_test_func_3_evaluate_maps @ Evaluate func. <>= procedure :: evaluate_func => vamp2_test_func_3_evaluate_func <>= real(default) function vamp2_test_func_3_evaluate_func (self, x) result (f) class(vamp2_test_func_3_t), intent(in) :: self real(default), dimension(:), intent(in) :: x real(default) :: xx select case (self%current_channel) case (1) xx = x(1)**0.2_default case (2) xx = 1 - x(1)**0.2_default end select f = 5 * xx**4 + 5 * (1 - xx)**4 end function vamp2_test_func_3_evaluate_func @ %def vamp2_test_func_3_evaluate_func @ \subsubsection{MC Integrator check} \label{sec:mc-integrator-check} We reproduce the first test case of VEGAS. Initialise the VAMP2 MC integrator and call to [[vamp2_init_grid]] for the initialisation of the grid. <>= call test (vamp2_1, "vamp2_1", "VAMP2 initialisation and& & grid preparation", u, results) <>= public :: vamp2_1 <>= subroutine vamp2_1 (u) integer, intent(in) :: u type(vamp2_t) :: mc_integrator class(rng_t), allocatable :: rng class(vamp2_func_t), allocatable :: func real(default), dimension(3), parameter :: x_lower = 0., & x_upper = pi real(default) :: result, abserr write (u, "(A)") "* Test output: vamp2_1" write (u, "(A)") "* Purpose: initialise the VAMP2 MC integrator and the grid" write (u, "(A)") write (u, "(A)") "* Initialise random number generator (default seed)" write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_channel = 1 and n_dim = 3" write (u, "(A)") allocate (vamp2_test_func_t :: func) call func%init (n_dim = 3, n_channel = 1) mc_integrator = vamp2_t (1, 3) call mc_integrator%write (u) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 10000" write (u, "(A)") call mc_integrator%set_limits (x_lower, x_upper) call mc_integrator%set_calls (10000) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 10000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 2000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (2000) call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () call rng%final () deallocate (rng) end subroutine vamp2_1 @ %def vamp2_1 @ Integrate a function with two-dimensional argument and two channels. <>= call test (vamp2_2, "vamp2_2", "VAMP2 intgeration of two-dimensional & & function with two channels", u, results) <>= public :: vamp2_2 <>= subroutine vamp2_2 (u) integer, intent(in) :: u type(vamp2_t) :: mc_integrator class(rng_t), allocatable :: rng class(vamp2_func_t), allocatable :: func real(default), dimension(2), parameter :: x_lower = 0., & x_upper = 1. real(default) :: result, abserr write (u, "(A)") "* Test output: vamp2_2" write (u, "(A)") "* Purpose: intgeration of two-dimensional & & function with two channels" write (u, "(A)") write (u, "(A)") "* Initialise random number generator (default seed)" write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_channel = 1 and n_dim = 3" write (u, "(A)") allocate (vamp2_test_func_2_t :: func) call func%init (n_dim = 2, n_channel = 2) mc_integrator = vamp2_t (2, 2) call mc_integrator%write (u) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 10000" write (u, "(A)") call mc_integrator%set_limits (x_lower, x_upper) call mc_integrator%set_calls (1000) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 10000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, opt_verbose = .true., result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 2000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (200) call mc_integrator%integrate (func, rng, 3, opt_verbose = .true., result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () call rng%final () deallocate (rng) end subroutine vamp2_2 @ %def vamp2_2 @ Integrate a function with two-dimensional argument and two channels. <>= call test (vamp2_3, "vamp2_3", "VAMP2 intgeration of two-dimensional & & function with two channels", u, results) <>= public :: vamp2_3 <>= subroutine vamp2_3 (u) integer, intent(in) :: u type(vamp2_t) :: mc_integrator class(rng_t), allocatable :: rng class(vamp2_func_t), allocatable :: func real(default), dimension(2), parameter :: x_lower = 0., & x_upper = 1. real(default) :: result, abserr integer :: unit write (u, "(A)") "* Test output: vamp2_3" write (u, "(A)") "* Purpose: intgeration of two-dimensional & & function with two channels" write (u, "(A)") write (u, "(A)") "* Initialise random number generator (default seed)" write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_channel = 1 and n_dim = 3" write (u, "(A)") allocate (vamp2_test_func_2_t :: func) call func%init (n_dim = 2, n_channel = 2) mc_integrator = vamp2_t (2, 2) call mc_integrator%write (u) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 20000" write (u, "(A)") call mc_integrator%set_limits (x_lower, x_upper) call mc_integrator%set_calls (20000) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 20000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Write grid to file vamp2_3.grids" write (u, "(A)") unit = free_unit () open (unit, file = "vamp2_3.grids", & action = "write", status = "replace") call mc_integrator%write_grids (unit) close (unit) write (u, "(A)") write (u, "(A)") "* Read grid from file vamp2_3.grids" write (u, "(A)") call mc_integrator%final () unit = free_unit () open (unit, file = "vamp2_3.grids", & action = "read", status = "old") call mc_integrator%read_grids (unit) close (unit) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 5000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (5000) call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () call rng%final () deallocate (rng) end subroutine vamp2_3 @ %def vamp2_3 @ Integrate a function with two-dimensional argument and two channels. Use chained weights, although we average over each weight itself. <>= call test (vamp2_4, "vamp2_4", "VAMP2 intgeration of two-dimensional & & function with two channels with chains", u, results) <>= public :: vamp2_4 <>= subroutine vamp2_4 (u) integer, intent(in) :: u type(vamp2_t) :: mc_integrator class(rng_t), allocatable :: rng class(vamp2_func_t), allocatable :: func real(default), dimension(2), parameter :: x_lower = 0., & x_upper = 1. real(default) :: result, abserr integer :: unit write (u, "(A)") "* Test output: vamp2_4" write (u, "(A)") "* Purpose: intgeration of two-dimensional & & function with two channels with chains" write (u, "(A)") write (u, "(A)") "* Initialise random number generator (default seed)" write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_channel = 2 and n_dim = 2" write (u, "(A)") allocate (vamp2_test_func_2_t :: func) call func%init (n_dim = 2, n_channel = 2) mc_integrator = vamp2_t (2, 2) call mc_integrator%write (u) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 20000 and set chains" write (u, "(A)") call mc_integrator%set_limits (x_lower, x_upper) call mc_integrator%set_calls (20000) call mc_integrator%set_chain (2, [1, 2]) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 10000 (Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Write grid to file vamp2_4.grids" write (u, "(A)") unit = free_unit () open (unit, file = "vamp2_4.grids", & action = "write", status = "replace") call mc_integrator%write_grids (unit) close (unit) write (u, "(A)") write (u, "(A)") "* Read grid from file vamp2_4.grids" write (u, "(A)") call mc_integrator%final () unit = free_unit () open (unit, file = "vamp2_4.grids", & action = "read", status = "old") call mc_integrator%read_grids (unit) close (unit) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 5000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (5000) call mc_integrator%integrate (func, rng, 3, result=result, abserr=abserr) write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ")") "Result: ", result, " +/- ", abserr write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () call rng%final () deallocate (rng) end subroutine vamp2_4 @ %def vamp2_4 @ <>= call test (vamp2_5, "vamp2_5", "VAMP2 intgeration of two-dimensional & & function with two channels with equivalences", u, results) <>= public :: vamp2_5 <>= subroutine vamp2_5 (u) integer, intent(in) :: u type(vamp2_t) :: mc_integrator class(rng_t), allocatable :: rng class(vamp2_func_t), allocatable :: func real(default), dimension(1), parameter :: x_lower = 0., & x_upper = 1. real(default) :: result, abserr integer :: unit type(vamp2_config_t) :: config type(vamp2_equivalences_t) :: eqv type(vegas_grid_t), dimension(2) :: grid write (u, "(A)") "* Test output: vamp2_5" write (u, "(A)") "* Purpose: intgeration of two-dimensional & & function with two channels with equivalences" write (u, "(A)") write (u, "(A)") "* Initialise random number generator (default seed)" write (u, "(A)") allocate (rng_stream_t :: rng) call rng%init () call rng%write (u) write (u, "(A)") write (u, "(A)") "* Initialise MC integrator with n_channel = 2 and n_dim = 1" write (u, "(A)") allocate (vamp2_test_func_3_t :: func) call func%init (n_dim = 1, n_channel = 2) config%equivalences = .true. mc_integrator = vamp2_t (n_channel = 2, n_dim = 1) call mc_integrator%set_config (config) call mc_integrator%write (u) write (u, "(A)") write (u, "(A)") "* Initialise grid with n_calls = 20000 and set chains" write (u, "(A)") call mc_integrator%set_limits (x_lower, x_upper) call mc_integrator%set_calls (20000) write (u, "(A)") write (u, "(A)") "* Initialise equivalences" write (u, "(A)") eqv = vamp2_equivalences_t (n_eqv = 4, n_channel = 2, n_dim = 1) call eqv%set_equivalence & (i_eqv = 1, dest = 2, src = 1, perm = [1], mode = [VEQ_IDENTITY]) call eqv%set_equivalence & (i_eqv = 2, dest = 1, src = 2, perm = [1], mode = [VEQ_IDENTITY]) call eqv%set_equivalence & (i_eqv = 3, dest = 1, src = 1, perm = [1], mode = [VEQ_IDENTITY]) call eqv%set_equivalence & (i_eqv = 4, dest = 2, src = 2, perm = [1], mode = [VEQ_IDENTITY]) call eqv%write (u) call mc_integrator%set_equivalences (eqv) write (u, "(A)") write (u, "(A)") & "* Integrate with n_it = 3 and n_calls = 10000 (Grid-only Adaptation)" write (u, "(A)") call mc_integrator%integrate (func, rng, 3, & opt_adapt_weight = .false., result=result, abserr=abserr) if (nearly_equal & (result, 2.000_default, rel_smallness = 0.003_default)) then write (u, "(2x,A)") "Result: 2.000 [ok]" else write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ",A)") & "Result: ", result, " +/- ", abserr, " [not ok]" end if write (u, "(A)") write (u, "(A)") "* Compare the grids of both channels" write (u, "(A)") grid(1) = mc_integrator%get_grid(channel = 1) grid(2) = mc_integrator%get_grid(channel = 2) write (u, "(2X,A,1X,L1)") "Equal grids =", (grid(1) == grid(2)) write (u, "(A)") write (u, "(A)") "* Write grid to file vamp2_5.grids" write (u, "(A)") unit = free_unit () open (unit, file = "vamp2_5.grids", & action = "write", status = "replace") call mc_integrator%write_grids (unit) close (unit) write (u, "(A)") write (u, "(A)") "* Integrate with n_it = 3 and n_calls = 5000 (Precision)" write (u, "(A)") call mc_integrator%set_calls (5000) call mc_integrator%integrate (func, rng, 3, opt_adapt_weight = .false., & opt_refine_grid = .false., result=result, abserr=abserr) if (nearly_equal & (result, 2.000_default, rel_smallness = 0.002_default)) then write (u, "(2x,A)") "Result: 2.000 [ok]" else write (u, "(2x,A," // FMT_12 // ",A," // FMT_12 // ",A)") & "Result: ", result, " +/- ", abserr, " [not ok]" end if write (u, "(A)") write (u, "(A)") "* Cleanup" call mc_integrator%final () call rng%final () deallocate (rng) end subroutine vamp2_5