algorithm2_1!(p::Vector{T}, I::Vector{Int}, w::Vector{<:Real}) where {T<:Real}

Fill p with the probabilities that result from normalizing the weights w[I]. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm2_1

algorithm2_1(I::Vector{Int}, w::Vector{<:Real})

Create a vector of probabilities by normalizing the weights selected by I from w. It is assumed that 0 ≤ wᵢ < Inf, i ∈ I.

Mathematically, given:

I ∈ ℕᴺ, 𝐰 ∈ ℝᴰ; I ⊆ {1,…,D}

The iᵗʰ term will be computed as: pᵢ = 𝐰ᵢ / ∑ⱼ 𝐰ⱼ; j ∈ I

See also: algorithm2_1!, algorithm2_2

Examples

julia> I = [1, 5, 2]; w = [5, 4, 3, 2, 1];

julia> algorithm2_1(I, w)
3-element Vector{Float64}:
 0.5
 0.1
 0.4

julia> algorithm2_1([1, 1, 2], Rational.(w))
3-element Vector{Rational{Int64}}:
 5//14
 5//14
 2//7

julia> w[2] = -w[2];

julia> algorithm2_1(I, w)                # Nonsense results if `wᵢ` constraints violated
3-element Vector{Float64}:
  2.5
  0.5
 -2.0

julia> algorithm2_1(I, [5, 4, 3, 2, Inf])
3-element Vector{Float64}:
   0.0
 NaN
   0.0

julia> algorithm2_1(I, [5, NaN, 3, 2, 1])
3-element Vector{Float64}:
 NaN
 NaN
 NaN

algorithm2_1_algorithm3!(p::Vector{T}, I::Vector{Int}, w::Vector{<:Real}, u::Real) where {T<:Real}

Normalize w[I] to probabilities, storing the result in p, then spreading probability mass 0 ≤ u ≤ 1 across the 0 or more elements of w[I] which are equal to zero. If all values of w[I] are zero and u ≠ 0, p will be filled with uniform probability mass. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm2_1_algorithm3

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; u = 0.5;

julia> algorithm2_1_algorithm3!(similar(I, Rational{Int}), I, w, u)
4-element Vector{Rational{Int64}}:
 1//6
 1//4
 1//4
 1//3

algorithm2_1_algorithm3(I::Vector{Int}, w::Vector{<:Real}, u::Real)

Return a vector of probabilities, normalizing the components selected from w by the index set I, then spreading the probability mass 0 ≤ u ≤ 1 across the 0 or more elements which are equal to zero. If all values of w[I] are zero and u ≠ 0, a vector of uniform probability mass is returned. Equivalent to algorithm3!(algorithm2_1(I, w), u) but more efficient.

See also: algorithm2_1_algorithm3!, algorithm2_1, algorithm3

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; u = 0.5;

julia> algorithm2_1_algorithm3(I, w, u)
4-element Vector{Float64}:
 0.16666666666666666
 0.25
 0.25
 0.3333333333333333

julia> algorithm3!(algorithm2_1(I, w), u)
4-element Vector{Float64}:
 0.16666666666666666
 0.25
 0.25
 0.3333333333333333

algorithm2_1_algorithm3!(p::Vector{T}, I::Vector{Int}, w::Vector{<:Real}, u::Real) where {T<:Real}

Normalize w[I] to probabilities, storing the result in p, then spreading probability mass u = r / (1 + r) across the 0 or more elements of w[I] which are equal to zero such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm2_1_algorithm3_ratio, algorithm2_1!, algorithm3_ratio!

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; r = 1.0;

julia> algorithm2_1_algorithm3_ratio!(similar(I, Rational{Int}), I, w, r)
4-element Vector{Rational{Int64}}:
 1//6
 1//4
 1//4
 1//3

algorithm2_1_algorithm3_ratio(I::Vector{Int}, w::Vector{<:Real}, r::Real)

Return a vector of probabilities, normalizing the components selected from w by the index set I, then spreading the probability mass u = r / (1 + r) across the 0 or more elements of w[I] which are equal to zero such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. If all values of w[I] are zero and r ≠ 0, a vector of uniform probability mass is returned. Equivalent to algorithm3_ratio!(algorithm2_1(I, w), u) but more efficient.

See also: algorithm2_1_algorithm3_ratio!, algorithm2_1, algorithm3_ratio

Examples

julia> I = [1, 2, 5, 6]; w = [10, 0, 30, 40, 0, 20]; r = 1.0;

julia> algorithm2_1_algorithm3_ratio(I, w, r)
4-element Vector{Float64}:
 0.16666666666666666
 0.25
 0.25
 0.3333333333333333

julia> algorithm3_ratio!(algorithm2_1(I, w), r)
4-element Vector{Float64}:
 0.16666666666666666
 0.25
 0.25
 0.3333333333333333

algorithm2_2!(p::Vector{T}, Is::NTuple{M, Vector{Int}}, ws::NTuple{M, Vector{<:Real}}) where {T<:Real, M}

Fill p with the probabilities that result from normalizing the element-wise product of weights selected by the index set, Is[m], respective to each weight vector, ws[m]. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm2_2

Examples

julia> Is = ([1,2,3], [4,5,6], [7,8,9]); ws = ([1,2,3], fill(1/6, 6), fill(1//10, 9));

julia> algorithm2_2!(zeros(3), Is, ws)
3-element Vector{Float64}:
 0.16666666666666666
 0.3333333333333333
 0.5

julia> algorithm2_2!(zeros(Rational{Int}, 3), Is, (ws[1], fill(1//6, 6), ws[3]))
3-element Vector{Rational{Int64}}:
 1//6
 1//3
 1//2

algorithm2_2(Is::NTuple{M, Vector{Int}}, ws::NTuple{M, Vector{<:Real}}) where {M}

Compute the product of weights selected by the respective index sets Is, then normalize the resultant weight vector to probabilities.

Mathematically, given:

I₁ ∈ ℕᴺ , 𝐰₁ ∈ ℝᴰ¹, 0 ≤ 𝐰₁ᵢ < Inf, i ∈ I₁

I₂ ∈ ℕᴺ , 𝐰₂ ∈ ℝᴰ², 0 ≤ 𝐰₂ᵢ < Inf, i ∈ I₂

⋮ , ⋮

Iₘ ∈ ℕᴺ , 𝐰ₘ ∈ ℝᴰᵐ, 0 ≤ 𝐰ₘᵢ < Inf, i ∈ Iₘ

The iᵗʰ term will be computed as: pᵢ = ∏ₘ₌₁ᴹ 𝐰ₘ[Iₘ[i]] / ∑ⱼ₌₁ᴺ ∏ₘ₌₁ᴹ 𝐰ₘ[Iₘ[j]]

See also: algorithm2_2!, algorithm2_1

Examples

julia> Is = ([1,2,3], [4,5,6], [7,8,9]); ws = ([1.0, 2.0, 3.0], fill(0.5, 6), fill(0.1, 9));

julia> algorithm2_2(Is, ws)
3-element Vector{Float64}:
 0.16666666666666666
 0.3333333333333333
 0.5

julia> w = ws[1][Is[1]] .* ws[2][Is[2]] .* ws[3][Is[3]]    # unnormalized
3-element Vector{Float64}:
 0.05
 0.1
 0.15000000000000002

julia> w ./= sum(w)                                        # normalized
3-element Vector{Float64}:
 0.16666666666666666
 0.3333333333333333
 0.5

algorithm3!(p::Vector{T}, w::Vector{<:Real}, u::Real) where {T<:Real}

Normalize w to probabilities, storing the result in p, spreading probability mass 0 ≤ u ≤ 1 across the 0 or more elements of w which are equal to zero. If all values of w are zero and u ≠ 0, p will be filled with uniform probability mass. Note that T must be a type which is able to hold the result of inv(one(T)).

Examples

julia> w = [0, 10, 5, 1]; u = 0.5;

julia> algorithm3!(similar(w, Float64), w, u)
4-element Vector{Float64}:
 0.5
 0.3125
 0.15625
 0.03125

algorithm3!(p::Vector{T}, u::Real) where {T<:Real}

Normalize p to probabilities, spreading probability mass u across the 0 or more elements of p which are equal to zero. If all values of w are zero and u ≠ 0, p will be filled with uniform probability mass. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm3, algorithm3_ratio!

Examples

julia> algorithm3!(Rational{Int}[0, 10, 5, 0], 0.5)
4-element Vector{Rational{Int64}}:
 1//4
 1//3
 1//6
 1//4

algorithm3(w::Vector{<:Real}, u::Real)

Return the vector of probabilities created by normalizing w to probabilities, then spreading the probability mass 0 ≤ u ≤ 1 across the 0 or more elements of w which are equal to zero. If all values of w are zero and u ≠ 0, a vector of uniform probability mass is returned.

Mathematically, given:

𝐰 ∈ ℝᴺ, 0 ≤ wᵢ < ∞, u ∈ ℝ, 0 ≤ u ≤ 1, J = {i : 𝐰ᵢ = 0}

pᵢ =
    Case 1: if J ≠ ∅
            u / |J|                     if i ∈ J
            (1 - u) * 𝐰ᵢ / ∑ᵢ₌₁ᴺ 𝐰ᵢ     otherwise
    Case 2: if J = {1,…,N}
            u / (u * N)                 Equivalent to 1/N if 𝐰 ≠ ̲0

See also: algorithm3!, algorithm3_ratio

Examples

julia> algorithm3([0, 10, 5, 0], 0.5)
4-element Vector{Float64}:
 0.25
 0.3333333333333333
 0.16666666666666666
 0.25

julia> algorithm3(Rational{Int}[0, 0, 0], 0.25)    # 𝐰 = ̲0
3-element Vector{Rational{Int64}}:
 1//3
 1//3
 1//3

julia> algorithm3([0, 0], 0.0)                     # 𝐰 = ̲0 and u = 0
2-element Vector{Float64}:
 NaN
 NaN

julia> algorithm3([1, 2, 3], 0.9)                  # in absence of 0's, just normalize
3-element Vector{Float64}:
 0.16666666666666666
 0.3333333333333333
 0.5

algorithm3_ratio!(p::Vector{T}, w::Vector{<:Real}, r::Real) where {T<:Real}

Normalize w to probabilities, storing the result in p, then spread the probability mass u = r / (1 + r) across the 0 or more elements of w such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. Note that T must be a type which is able to hold the result of inv(one(T)).

Examples

julia> w = [0, 10, 5, 1]; r = 1.0;

julia> algorithm3_ratio!(similar(w, Float64), w, r)
4-element Vector{Float64}:
 0.5
 0.3125
 0.15625
 0.03125

algorithm3_ratio!(p::Vector{T}, r::Real) where {T<:Real}

Normalize p to probabilities, then spread the probability mass u = r / (1 + r) across the 0 or more elements of p such that the ratio of the sum of (inititally) zero elements to the sum of the non-zero elements is equal to r. Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm3_ratio, algorithm3!

algorithm3_ratio(w::Vector{<:Real}, r::Real)

Return a vector of probabilities created by normalizing w to probabilities, then spread the probability mass u = r / (1 + r) across the 0 or more elements of w which are equal to zero such that the ratio of the sum of (inititally) zero elements to the sum of non-zero elements is equal to r. If all values of w are zero and r ≠ 0, a vector of uniform probability mass is returned.

Mathematically, given:

𝐰 ∈ ℝᴺ, 0 ≤ wᵢ < ∞, r ∈ ℝ, 0 ≤ r ≤ Inf, J = {i : 𝐰ᵢ = 0}

pᵢ =
    Case 1: if J ≠ ∅
            (r / (1+r)) / |J|                     if i ∈ J
            (1 / (1+r)) * 𝐰ᵢ / ∑ᵢ₌₁ᴺ 𝐰ᵢ           otherwise
    Case 2: if J = {1,…,N}
            r / (r * N)                           Equivalent to 1/N if 𝐰 ≠ ̲0

See also: algorithm3_ratio!, algorithm3

Examples

julia> w = Rational{Int}[1, 0, 3, 0, 5]; r = 3;

julia> p = algorithm3_ratio(w, r)
5-element Vector{Rational{Int64}}:
 1//36
 3//8
 1//12
 3//8
 5//36

julia> r′ = sum(p[findall(iszero, w)]) / sum(p[findall(!iszero, w)]); (r′, r′ == r)
(3//1, true)

julia> algorithm3(w, r / (1 + r))              # Note equivalence
5-element Vector{Rational{Int64}}:
 1//36
 3//8
 1//12
 3//8
 5//36

julia> algorithm3_ratio(w, Inf)                # r = Inf ⟹ u = 1
5-element Vector{Rational{Int64}}:
 0//1
 1//2
 0//1
 1//2
 0//1

algorithm4!(p::Vector{T}, 𝐰₁::Vector{<:Real}, 𝐰₂::Vector{<:Real}) where {T<:Real}

Fill p with the probabilities which result from algorithm4(𝐰₁, 𝐰₂). Note that T must be a type which is able to hold the result of inv(one(T)).

algorithm4!(𝐰₁::Vector{T}, 𝐰₂::Vector{<:Real}) where {T<:Real}

Fill 𝐰₁ with the probabilities which result from algorithm4(𝐰₁, 𝐰₂). Note that T must be a type which is able to hold the result of inv(one(T)).

See also: algorithm4

algorithm4(𝐰₁::Vector{<:Real}, 𝐰₂::Vector{<:Real})

Return a vector of probabilities constructed according to the following algorithm:

Define:

I = {1,…,N}

J₁ = {i: 𝐰₁ᵢ = 0}, I₁′ = {i: 𝐰₁ᵢ ≠ 0} = I ∖ J₁

J₂ = {i: 𝐰₂ᵢ = 0}, I₂′ = {i: 𝐰₂ᵢ ≠ 0} = I ∖ J₂

𝐰₁ ∈ ℝᴺ : initial weights, 0 ≤ 𝐰₁ᵢ < Inf

𝐰₂ ∈ ℝᴺ : augment weights, 0 ≤ 𝐰₂ᵢ < Inf; a value of zero indicates no re-weight

Then:

pᵢ = 𝐰₁ᵢ / ∑ₗ₌₁ᴺ 𝐰₁ₗ, i ∈ I ∖ I₂′

pᵢ = (𝐰₂ᵢ * ∑ₗ 𝐰₁ₗ, l ∈ I₂′) / (∑ₗ₌₁ᴺ 𝐰₂ₗ * ∑ₗ₌₁ᴺ 𝐰₁ₗ), i ∈ I₂′

This algorithm can produce a wide variety of probability vectors as the result of the various combinations of intersections which can be formed from J₁, J₂, I₁′, and I₂′. However, complexity of outputs aside, the motivating concept is quite simple: take a vector of weights, 𝐰₁ and re-weight some subset (I₂′) of those weights using a second set of weights, 𝐰₂, while preserving the proportion of probability mass derived from 𝐰₁. That is, given p = algorithm4(𝐰₁, 𝐰₂), the following relationships are preserved: sum(p[J₂]) ≈ sum(𝐰₁[J₂]) / sum(𝐰₁[I₁′]), sum(w₁[J₂]) / sum(w₁[I₂′]) ≈ sum(p[J₂]) / sum(p[I₂′]).

See also: algorithm4!

Examples

julia> w₁ = [1, 1, 1, 1, 0];

julia> algorithm4(w₁, [2, 1, 3, 4, 0])    # J₁ ∩ I₂′ = ∅
5-element Vector{Float64}:
 0.2
 0.1
 0.30000000000000004
 0.4
 0.0

julia> algorithm4(w₁, [2, 1, 3, 0, 5])    # J₂ = [4] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
 0.13636363636363635
 0.06818181818181818
 0.20454545454545453
 0.25
 0.3409090909090909

julia> w₁ = [1, 1, 1, 0, 0];

julia> algorithm4(w₁, [2, 1, 3, 4, 0])    # J₂ = [5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
 0.2
 0.1
 0.30000000000000004
 0.4
 0.0

julia> w₁ = [1, 1, 0, 1, 0];

julia> algorithm4(w₁, [0, 1, 0, 4, 0])    # J₂ = [1,3,5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
 0.3333333333333333
 0.13333333333333333
 0.0
 0.5333333333333333
 0.0

julia> algorithm4(w₁, [0, 0, 3, 4, 0])    # J₂ = [1,2,5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
 0.3333333333333333
 0.3333333333333333
 0.14285714285714285
 0.19047619047619047
 0.0

julia> algorithm4(w₁, [2, 0, 3, 0, 0])    # J₂ = [2,4,5] not re-weighted; I₂′ re-weighted
5-element Vector{Float64}:
 0.13333333333333333
 0.3333333333333333
 0.2
 0.3333333333333333
 0.0

normalize1!(A::AbstractArray{<:Real})

Normalize the values in A such that sum(A) ≈ 1 and 0 ≤ A[i] ≤ 1 ∀i. This is not quite the L¹-norm, which would require that abs(A[i]) be used. It is assumed that 0 ≤ A[i] < Inf ∀i. Inf values are not handled and will result in NaN's.

See also: normalize1

julia> normalize1!([1.0, 2.0, 3.0])
3-element Vector{Float64}:
 0.16666666666666666
 0.3333333333333333
 0.5

julia> normalize1!([1.0, 2.0, Inf])
3-element Vector{Float64}:
   0.0
   0.0
 NaN

julia> normalize1!([1.0, 2.0, NaN])     # NaN propagates, as expected
3-element Vector{Float64}:
 NaN
 NaN
 NaN

julia> normalize1!([1.0, -2.0, 3.0])    # not the L¹-norm
3-element Vector{Float64}:
  0.5
 -1.0
  1.5

normalize1!(B::AbstractArray{<:Real}, A::AbstractArray{<:Real})

Normalize the values in A such that sum(B) ≈ 1 and 0 ≤ B[i] ≤ 1 ∀i, storing the result in B. It is assumed that A[i] ≥ 0 ∀i.

normalize1(A::AbstractArray{<:Real})

Return an array of equal size which satisfies sum(B) ≈ 1 and 0 ≤ B[i] 1 ∀i. It is assumed that A[i] ≥ 0 ∀i.

See also: normalize1!

sample!(B::AbstractArray, A::AbstractArray)

Draw samples, summing (and potentially reducing) in-place into B. The shape of B determines the extent of reduction performed.

See also: tsample!, vsample!, vtsample!

sample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing])

Draw n_sim samples from the distribution which corresponds to the sum of the independent categorical distributions defined by the probability mass vector(s), A, storing the result in an array of eltype T of sufficient size and dimension as determined by the input, A, the number of categories, n_cat, and the (potential) reduction dimensions, dims.

dims is an optional keyword argument used to specify an in-place sum on the indices of A (if A is an array of arrays). n_cat may be optionally specified, or inferred from the available data; validity indices will be checked regardless.

In the simplest case, A may be a single probability mass vector, for which compatible types are AbstractVector{<:Real}, SparseVector{<:Real, <:Integer}, and AbstractVector{Tuple{<:Integer, <:Real}}, the last of which is simply another representation of a sparse vector.

In the second case, A may be an AbstractArray{V, N} where {V, N} in which V is any of the types specified above for a single probability mass vector. That is, an array of vectors.

In the third case, A may be an AbstractArray{W, M} where {W<:AbstractArray{V, N}, M} where {V, N} in which W is any of the types specified in the second case. In other words, an array of arrays of vectors.

See also: tsample, vsample, vtsample

splitranges(start, stop, chunksize)

Divide the range start:stop into segments, each of size chunksize. The last segment will contain the remainder, (start - stop + 1) % chunksize, if it exists.

splitranges(ur::UnitRange{Int}, chunksize)

Divide the range ur into segments, each of size chunksize.

tsample!(B::AbstractArray, A::AbstractArray; [chunksize=5000])

Identical to sample! except that thread-based parallelism is used to accelerate the computation.

See also: sample!, vsample!, vtsample!

tsample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing], [chunksize=5000])

See sample for full documentation; identical in behavior except that thread-based parallelism is used to accelerate the computation.

The optional chunksize argument provides a fixed upper bound for the number of simulations to be passed to each thread. Depending on number of independent categorical distributions to be sample and number of simulations to be performed, a chunksize of 5000 is likely too conservative. The user is encouraged to try smaller chunk sizes.

See also: sample, vsample, vtsample

vsample!(B::AbstractArray, A::AbstractArray)

Identical to sample! except that except that the underlying Marsaglia sampler uses LoopVectorization.

See also: vtsample!, sample!, tsample!

vsample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing])

See sample for full documentation; identical in behavior except that the underlying Marsaglia sampler uses LoopVectorization. Sometimes, but not always, confers speedup.

See also: vtsample, sample, tsample

tsample!(B::AbstractArray, A::AbstractArray; [chunksize=5000])

Identical to sample! except that thread-based parallelism is used to accelerate the computation. also, the underlying Marsaglia sampler uses LoopVectorization.

See also: vsample!, sample!, tsample!

tsample(::Type{T<:Real}=Int, A::AbstractArray, n_sim::Int; [dims=:], [n_cat=nothing], [chunksize=5000])

See sample for full documentation; identical in behavior except that thread-based parallelism is used to accelerate the computation; also, the underlying Marsaglia sampler uses LoopVectorization.

The optional chunksize argument provides a fixed upper bound for the number of simulations to be passed to each thread. Depending on number of independent categorical distributions to be sample and number of simulations to be performed, a chunksize of 5000 is likely too conservative. The user is encouraged to try smaller chunk sizes.

See also: vsample, sample, tsample