constrained_sample_nI Subroutine

private subroutine constrained_sample_nI(this, contain, renormalization, pos, tgt, prob)

draw a random element from 1:size(this%probs) with the probabilities listed in prob @param[in] constraint pick only elements from constraint @param[out] tgt on return, this is a random number in the sampling range of this @param[out] pos the position of tgt in contain @param[out] prob on return, the probability of picking tgt from constraint

Type Bound

AliasSampler_t

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(in) :: this
integer, intent(in) :: contain(:)
real(kind=dp), intent(in) :: renormalization
integer, intent(out) :: pos
integer, intent(out) :: tgt
real(kind=dp), intent(out) :: prob

Contents

Source Code

constrained_sample_nI

Source Code

    subroutine constrained_sample_nI(this, contain, renormalization, pos, tgt, prob)
        class(AliasSampler_t), intent(in) :: this
        integer, intent(in) :: contain(:)
        real(dp), intent(in) :: renormalization
        integer, intent(out) :: pos, tgt
        real(dp), intent(out) :: prob
        character(*), parameter :: this_routine = 'constrained_sample_nI'

        ASSERT(is_set(contain))
        ASSERT(1 <= contain(1) .and. contain(size(contain)) <= size(this%probs%ptr))
        ASSERT(renormalization .isclose. (sum(this%get_prob(contain))))

        if (near_zero(renormalization)) then
            tgt = 0
            prob = 1.0
            return
        end if

        if (renormalization < redrawing_cutoff) then
        block
            type(CDF_Sampler_t) :: cdf_sampler
            cdf_sampler = CDF_Sampler_t(this%probs%ptr(contain), renormalization)
            call cdf_sampler%sample(pos, prob)
            tgt = contain(pos)
            return
        end block
        end if

        tgt = this%table%sample()
        pos = int(binary_search_int(contain, tgt))
        do while (pos == -1)
            tgt = this%table%sample()
            pos = int(binary_search_int(contain, tgt))
        end do
        prob = this%probs%ptr(tgt) / renormalization
        ASSERT(prob .isclose. (this%probs%ptr(tgt) / sum(this%probs%ptr(contain))))
        ASSERT(contain(pos) == tgt)
    end subroutine constrained_sample_nI