aliasSampling Module

This module contains the utility to use alias table lookup on lists, requiring to precompute biases but making the lookup O(1)


Uses


Variables

Type Visibility Attributes Name Initial
real(kind=dp), private, parameter :: redrawing_cutoff = 0.1_dp

If we draw from constrained subsets of a precomputed probability distributions we can use two different algorithms: 1. We just redraw until we sample an element from our subset 2. We reconstruct probability distributions for our narrower subset. If the sum of probabilities of our subset is larger than redrawing_cutoff, we take the first method, otherwise the second.

Derived Types

type, public ::  AliasTable_t

This class implements the Walker-Vose Alias Sampling method.

Read more…

Components

Type Visibility Attributes Name Initial
type(shared_array_real_t), private :: bias

this is the table of bias
type(shared_array_int32_t), private :: alias

this is the lookup table for the resulting random number

Type-Bound Procedures

procedure, public :: setup => setup_AliasTable_t

constructor
procedure, private :: init => init_AliasTable_t

only compute the data, without allocation
procedure, public :: finalize => finalize_AliasTable_t

destructor - final would be suited better, but is not supported by all compilers
procedure, public :: sample => sample_AliasTable_t

get a random value from the alias table
procedure, public :: get_memory_demand => get_memory_demand_AliasTable_t

type, public ::  AliasSampler_t

Alias Sampler class

Read more…

Components

Type Visibility Attributes Name Initial
type(AliasTable_t), private :: table

alias table used for sampling
type(shared_array_real_t), private :: probs

the probabilities

Type-Bound Procedures

procedure, public :: setup => setup_AliasSampler_t

constructor
procedure, private :: init => init_AliasSampler_t

only compute the data, i.e. bias tables etc., without allocation
procedure, private :: init_probs => init_probs_AliasSampler_t

initialize the probabilities
procedure, public :: finalize => finalize_AliasSampler_t

destructor
procedure, public :: sample

get a random element and the generation probability
generic, public :: constrained_sample => constrained_sample_fast, constrained_sample_nI

get a random element from a constrained set and its normalized! generation probability
procedure, private :: constrained_sample_fast
procedure, private :: constrained_sample_nI
procedure, public :: get_prob => get_prob_AliasSampler_t

get the probability to produce a given value
procedure, public :: constrained_getProb

get the probability to draw a given value from a constrained set
procedure, public :: get_memory_demand => get_memory_demand_AliasSampler

type, public ::  AliasSampler_3D_t

Components

Type Visibility Attributes Name Initial
type(AliasSampler_t), private, allocatable :: samplerArray(:,:,:)
type(shared_array_real_t), private :: allProbs
type(shared_array_real_t), private :: allBiasTable
type(shared_array_int32_t), private :: allAliasTable

Type-Bound Procedures

procedure, public :: shared_alloc => setupSamplerArray_3D
procedure, public :: setup_entry => setupEntry_3D
procedure, public :: finalize => samplerArrayDestructor_3D
procedure, public :: sample => aSample_3D
generic, public :: constrained_sample => constrained_sample_3D_nI, constrained_sample_3D_fast
procedure, private :: constrained_sample_3D_fast
procedure, private :: constrained_sample_3D_nI
procedure, public :: get_prob => aGetProb_3D
procedure, public :: constrained_getProb => constrained_get_prob_3D
procedure, public :: get_memory_demand => get_memory_demand_AliasSampler_3D_t

type, public ::  AliasSampler_2D_t

Components

Type Visibility Attributes Name Initial
type(AliasSampler_3D_t), private :: alias_sampler

Type-Bound Procedures

procedure, public :: shared_alloc => setupSamplerArray_2D
procedure, public :: setup_entry => setupEntry_2D
procedure, public :: finalize => samplerArrayDestructor_2D
procedure, public :: sample => aSample_2D
generic, public :: constrained_sample => constrained_sample_2D_nI, constrained_sample_2D_fast
procedure, private :: constrained_sample_2D_fast
procedure, private :: constrained_sample_2D_nI
procedure, public :: get_prob => aGetProb_2D
procedure, public :: constrained_getProb => constrained_get_prob_2D
procedure, public :: get_memory_demand => get_memory_demand_AliasSampler_2D_t

type, public ::  AliasSampler_1D_t

Components

Type Visibility Attributes Name Initial
type(AliasSampler_3D_t), private :: alias_sampler

Type-Bound Procedures

procedure, public :: shared_alloc => setupSamplerArray_1D
procedure, public :: setup_entry => setupEntry_1D
procedure, public :: finalize => samplerArrayDestructor_1D
procedure, public :: sample => aSample_1D
generic, public :: constrained_sample => constrained_sample_1D_nI, constrained_sample_1D_fast
procedure, private :: constrained_sample_1D_fast
procedure, private :: constrained_sample_1D_nI
procedure, public :: get_prob => aGetProb_1D
procedure, public :: constrained_getProb => constrained_get_prob_1D
procedure, public :: get_memory_demand => get_memory_demand_AliasSampler_1D_t

Functions

private elemental function get_memory_demand_AliasTable_t(this) result(res)

Arguments

Type IntentOptional Attributes Name
class(AliasTable_t), intent(in) :: this

Return Value type(ByteSize_t)

private function sample_AliasTable_t(this) result(ind)

Draw a random number from an alias table created with the corresponding probabilities aliasTable object

Arguments

Type IntentOptional Attributes Name
class(AliasTable_t), intent(in) :: this

Return Value integer

private elemental function get_memory_demand_AliasSampler(this) result(res)

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(in) :: this

Return Value type(ByteSize_t)

private elemental function get_prob_AliasSampler_t(this, tgt) result(prob)

Returns the probability to draw tgt from this sampler

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(in) :: this
integer, intent(in) :: tgt

the number for which we request the probability of sampling

Return Value real(kind=dp)

the probability of drawing tgt with the sample routine

private pure function constrained_getProb(this, contain, renorm, tgt) result(prob)

Returns the probability to draw tgt from this sampler (has to be a set, i.e. unique and ordered) the probability of drawing anything from an empty sampler is 0

Arguments

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

the number for which we request the probability of sampling

Return Value real(kind=dp)

This loosened threshhold might be a good idea if the renormalization was calculated via the complement, i.e. \Sum{ p_i } {i \in D_i} = 1 - \Sum{ p_i } {i \notin D_i}

the probability of picking tgt from constraint

private elemental function get_memory_demand_AliasSampler_1D_t(this) result(res)

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(in) :: this

Return Value type(ByteSize_t)

private elemental function aGetProb_1D(this, iEntry, tgt) result(prob)

Returns the probability to draw tgt from the sampler with index iEntry

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(in) :: this
integer, intent(in) :: iEntry

index of the sampler to use
integer, intent(in) :: tgt

the number for which we request the probability of sampling

Return Value real(kind=dp)

private pure function constrained_get_prob_1D(this, i, contain, renorm, tgt) result(prob)

Returns the probability to draw tgt from the sampler with index iEntry

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(in) :: this
integer, intent(in) :: i

Index of the sampler to use
integer, intent(in) :: contain(:)
real(kind=dp), intent(in) :: renorm
integer, intent(in) :: tgt

the number for which we request the probability of sampling

Return Value real(kind=dp)

private elemental function get_memory_demand_AliasSampler_2D_t(this) result(res)

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(in) :: this

Return Value type(ByteSize_t)

private elemental function aGetProb_2D(this, i, j, tgt) result(prob)

Returns the probability to draw tgt from the sampler with index iEntry

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(in) :: this
integer, intent(in) :: i
integer, intent(in) :: j
integer, intent(in) :: tgt

the number for which we request the probability of sampling

Return Value real(kind=dp)

private pure function constrained_get_prob_2D(this, i, j, contain, renorm, tgt) result(prob)

Returns the probability to draw tgt from the sampler with index iEntry

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(in) :: this
integer, intent(in) :: i

Index of the sampler to use
integer, intent(in) :: j

Index of the sampler to use Index of the sampler to use
integer, intent(in) :: contain(:)
real(kind=dp), intent(in) :: renorm
integer, intent(in) :: tgt

the number for which we request the probability of sampling

Return Value real(kind=dp)

private elemental function get_memory_demand_AliasSampler_3D_t(this) result(res)

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(in) :: this

Return Value type(ByteSize_t)

private elemental function aGetProb_3D(this, i, j, k, tgt) result(prob)

Returns the probability to draw tgt from the sampler with index iEntry

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(in) :: this
integer, intent(in) :: i

Index of the sampler to use
integer, intent(in) :: j

Index of the sampler to use Index of the sampler to use
integer, intent(in) :: k

Index of the sampler to use Index of the sampler to use Index of the sampler to use
integer, intent(in) :: tgt

the number for which we request the probability of sampling

Return Value real(kind=dp)

private pure function constrained_get_prob_3D(this, i, j, k, contain, renorm, tgt) result(prob)

Returns the probability to draw tgt from the sampler with index iEntry

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(in) :: this
integer, intent(in) :: i

Index of the sampler to use
integer, intent(in) :: j

Index of the sampler to use Index of the sampler to use
integer, intent(in) :: k

Index of the sampler to use Index of the sampler to use Index of the sampler to use
integer, intent(in) :: contain(:)
real(kind=dp), intent(in) :: renorm
integer, intent(in) :: tgt

the number for which we request the probability of sampling

Return Value real(kind=dp)

public elemental function do_direct_calculation(normalization)

Evaluate if a normalization has to be calculated directly.

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: normalization

Return Value logical


Subroutines

private subroutine setup_AliasTable_t(this, rank_with_info, arr)

pseudo-constructor for alias tables

Arguments

Type IntentOptional Attributes Name
class(AliasTable_t), intent(inout) :: this
integer, intent(in) :: rank_with_info

The intra-node rank that contains the weights to be used all other arr of all other ranks are ignored (and can be allocated with size 0).
real(kind=dp), intent(in) :: arr(:)

arr array containing the (not necessarily normalized) probabilities we want to use for sampling

private subroutine init_AliasTable_t(this, rank_with_info, arr)

Set the bias and alias values for each value in range want to use for sampling

Arguments

Type IntentOptional Attributes Name
class(AliasTable_t), intent(inout) :: this
integer, intent(in) :: rank_with_info

The intra-node rank that contains the weights to be used all other arr of all other ranks are ignored (and can be allocated with size 0).
real(kind=dp), intent(in) :: arr(:)

private subroutine finalize_AliasTable_t(this)

clear the memory used by the alias table

Arguments

Type IntentOptional Attributes Name
class(AliasTable_t), intent(inout) :: this

private subroutine setup_AliasSampler_t(this, rank_with_info, arr)

allocate the resources of this and load the probability distribution from arr into this want to use for sampling

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(inout) :: this
integer, intent(in) :: rank_with_info

The intra-node rank that contains the weights to be used all other arr of all other ranks are ignored (and can be allocated with size 0).
real(kind=dp), intent(in) :: arr(:)

private subroutine init_AliasSampler_t(this, rank_with_info, arr)

load the probability distribution from arr into this we only use this in the sampler array, but fortran has no friend classes, so its public want to use for sampling

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(inout) :: this
integer, intent(in) :: rank_with_info

The intra-node rank that contains the weights to be used all other arr of all other ranks are ignored (and can be allocated with size 0).
real(kind=dp), intent(in) :: arr(:)

private subroutine init_probs_AliasSampler_t(this, rank_with_info, arr)

load the probability distribution from arr into this%probs want to use for sampling

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(inout) :: this
integer, intent(in) :: rank_with_info

The intra-node rank that contains the weights to be used all other arr of all other ranks are ignored (and can be allocated with size 0).
real(kind=dp), intent(in) :: arr(:)

private subroutine finalize_AliasSampler_t(this)

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(inout) :: this

private subroutine sample(this, tgt, prob)

draw a random element from 1:size(this%probs) with the probabilities listed in prob

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_t), intent(in) :: this
integer, intent(out) :: tgt

on return, this is a random number in the sampling range of this
real(kind=dp), intent(out) :: prob

on return, the probability of picking tgt

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

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

the position of tgt in contain
integer, intent(out) :: tgt

the position of tgt in contain on return, this is a random number in the sampling range of this
real(kind=dp), intent(out) :: prob

on return, the probability of picking tgt from constraint

private subroutine constrained_sample_fast(this, contain, contain_ilut, renormalization, pos, val, prob)

draw a random element from 1:size(this%probs) with the probabilities listed in prob

Arguments

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

the position of tgt in contain
integer, intent(out) :: val

the position of tgt in contain
real(kind=dp), intent(out) :: prob

on return, the probability of picking tgt from constraint

public subroutine clear_sampler_array(arr)

call the destructor on all elements of an array, then deallocate it. This is for intrinsic arrays, the sampler array class has its own deallocate routine.

Arguments

Type IntentOptional Attributes Name
type(AliasSampler_t), intent(inout), allocatable :: arr(:)

private subroutine setupSamplerArray_1D(this, nEntries, entrySize, name)

Setup an array of samplers using a single shared resource (split into parts associated with one of them each). This only does the allocation.

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t) :: this
integer, intent(in) :: nEntries

number of samplers to initialise
integer, intent(in) :: entrySize

number of samplers to initialise number of values per sampler
character(len=*), intent(in) :: name

private subroutine setupEntry_1D(this, iEntry, rank_with_info, arr)

Initialise one sampler of an array

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(inout) :: this
integer, intent(in) :: iEntry

index of the entry to initialize
integer, intent(in) :: rank_with_info

index of the entry to initialize
real(kind=dp), intent(in) :: arr(:)

private subroutine samplerArrayDestructor_1D(this)

Deallocate an array of samplers

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(inout) :: this

private subroutine aSample_1D(this, iEntry, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(in) :: this
integer, intent(in) :: iEntry

The index of the sampler.
integer, intent(out) :: tgt

The sampled value tgt.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt.

private subroutine constrained_sample_1D_nI(this, i, contain, renorm, pos, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob while adherring to constraints

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: contain(:)

The constraint in nI format.
real(kind=dp), intent(in) :: renorm

The renormalization. (i.e. sum(this%get_prob(… contain…))
integer, intent(out) :: pos

The sampled value tgt and its position pos in `contain.
integer, intent(out) :: tgt

The sampled value tgt and its position pos in `contain.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt from contain

private subroutine constrained_sample_1D_fast(this, i, contain, contain_ilut, renormalization, pos, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob while adherring to constraints

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_1D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: contain(:)

The constraint in nI format.
integer(kind=n_int), intent(in) :: contain_ilut(0:)

The constraint in ilut (bitmask) format
real(kind=dp), intent(in) :: renormalization

The renormalization. (i.e. sum(this%get_prob(… contain…))
integer, intent(out) :: pos

The sampled value tgt and its position pos in `contain.
integer, intent(out) :: tgt

The sampled value tgt and its position pos in `contain.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt from contain

private subroutine setupSamplerArray_2D(this, dims, entrySize, name)

Setup an array of samplers using a single shared resource (split into parts associated with one of them each). This only does the allocation.

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t) :: this
integer, intent(in) :: dims(2)
integer, intent(in) :: entrySize

number of values per sampler
character(len=*), intent(in) :: name

private subroutine setupEntry_2D(this, i, j, rank_with_info, arr)

Initialise one sampler of an array

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(inout) :: this
integer, intent(in) :: i
integer, intent(in) :: j
integer, intent(in) :: rank_with_info
real(kind=dp), intent(in) :: arr(:)

private subroutine samplerArrayDestructor_2D(this)

Deallocate an array of samplers

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(inout) :: this

private subroutine aSample_2D(this, i, j, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: j

The index of the sampler.
integer, intent(out) :: tgt

The sampled value tgt.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt.

private subroutine constrained_sample_2D_nI(this, i, j, contain, renorm, pos, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob while adherring to constraints

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: j

The index of the sampler.
integer, intent(in) :: contain(:)

The constraint in nI format.
real(kind=dp), intent(in) :: renorm

The renormalization. (i.e. sum(this%get_prob(… contain…))
integer, intent(out) :: pos

The sampled value tgt and its position pos in `contain.
integer, intent(out) :: tgt

The sampled value tgt and its position pos in `contain.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt from contain

private subroutine constrained_sample_2D_fast(this, i, j, contain, contain_ilut, renormalization, pos, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob while adherring to constraints

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_2D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: j

The index of the sampler.
integer, intent(in) :: contain(:)

The constraint in nI format.
integer(kind=n_int), intent(in) :: contain_ilut(0:)

The constraint in ilut (bitmask) format
real(kind=dp), intent(in) :: renormalization

The renormalization. (i.e. sum(this%get_prob(… contain…))
integer, intent(out) :: pos

The sampled value tgt and its position pos in `contain.
integer, intent(out) :: tgt

The sampled value tgt and its position pos in `contain.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt from contain

private subroutine setupSamplerArray_3D(this, dims, entry_size, name)

Setup an array of samplers using a single shared resource (split into parts associated with one of them each). This only does the allocation.

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(inout) :: this
integer, intent(in) :: dims(3)
integer, intent(in) :: entry_size

number of values per sampler
character(len=*), intent(in) :: name

private subroutine setupEntry_3D(this, i, j, k, rank_with_info, arr)

@brief Initialise one sampler of an array

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(inout) :: this
integer, intent(in) :: i

Index of the entry to initialize
integer, intent(in) :: j

Index of the entry to initialize Index of the entry to initialize
integer, intent(in) :: k

Index of the entry to initialize Index of the entry to initialize Index of the entry to initialize
integer, intent(in) :: rank_with_info

The intra-node rank that contains the weights to be used all other arr of all other ranks are ignored (and can be allocated with size 0).
real(kind=dp), intent(in) :: arr(:)

private subroutine samplerArrayDestructor_3D(this)

@brief Deallocate an array of samplers

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(inout) :: this

private subroutine aSample_3D(this, i, j, k, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: j

The index of the sampler.
integer, intent(in) :: k

The index of the sampler.
integer, intent(out) :: tgt

The sampled value tgt.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt.

private subroutine constrained_sample_3D_nI(this, i, j, k, contain, renorm, pos, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob while adherring to constraints

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: j

The index of the sampler.
integer, intent(in) :: k

The index of the sampler.
integer, intent(in) :: contain(:)

The constraint in nI format.
real(kind=dp), intent(in) :: renorm

The renormalization. (i.e. sum(this%get_prob(… contain…))
integer, intent(out) :: pos

The sampled value tgt and its position pos in `contain.
integer, intent(out) :: tgt

The sampled value tgt and its position pos in `contain.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt from contain

private subroutine constrained_sample_3D_fast(this, i, j, k, contain, contain_ilut, renormalization, pos, tgt, prob)

Draw a random element from 1:size(this%probs) with the probabilities listed in prob while adherring to constraints

Arguments

Type IntentOptional Attributes Name
class(AliasSampler_3D_t), intent(in) :: this
integer, intent(in) :: i

The index of the sampler.
integer, intent(in) :: j

The index of the sampler.
integer, intent(in) :: k

The index of the sampler.
integer, intent(in) :: contain(:)

The constraint in nI format.
integer(kind=n_int), intent(in) :: contain_ilut(0:)

The constraint in ilut (bitmask) format
real(kind=dp), intent(in) :: renormalization

The renormalization. (i.e. sum(this%get_prob(… contain…))
integer, intent(out) :: pos

The sampled value tgt and its position pos in `contain.
integer, intent(out) :: tgt

The sampled value tgt and its position pos in `contain.
real(kind=dp), intent(out) :: prob

The probability of sampling tgt from contain