You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
#TODO: Implement exponential for any arbitrary range
DistFix{W, F}(mantissa)
end""" triangle(::Type{DistFix{W, F}}, b::Int)Returns a triangle distribution with a range of `b` bits. A triangle distribution over `b` bits refers to the distributions that assigns probability 0 to 0 and increases linearly with probability 2/(2^b)*(2^b - 1)"""functiontriangle(::Type{DistFix{W, F}}, b::Int) where {W,F}
DistFix{W, F}(triangle(DistInt{W}, b))
end""" unit_exponential(::Type{DistFix{W, F}}, beta::Float64; reverse=false)Returns an exponential distribution e^(beta*x) in the range [0, 1) of the specified type.For advanced users: `reverse` is used to control the order in which flips are created. `reverse=false` refers to the order of flips from MSB to LSB. """functionunit_exponential(::Type{DistFix{W, F}}, beta::Float64; reverse=false) where W where F
if!reverse
DistFix{W, F}(vcat([falsefor i in1:W-F], [flip(exp(beta/2^i)/(1+exp(beta/2^i))) for i in1:F]))
elseDistFix{W, F}(vcat([falsefor i in1:W-F], [flip(exp(beta/2^i)/(1+exp(beta/2^i))) for i in F:-1:1][F:-1:1]))
endend""" exponential(::Type{DistFix{W, F}}, beta::Float64, start::Float64, stop::Float64; reverse=false)Returns an exponential distribution e^(beta*x) in the range [start, stop) of the specified type."""functionexponential(::Type{DistFix{W, F}}, beta::Float64, start::Float64, stop::Float64; reverse=false) where W where F
range = stop - start
#TODO: Implement exponential for any arbitrary range@assertispow2(range) "Currently the function 'exponential' only supports range that is a power of 2."
new_beta = beta*range
bits =Int(log2(range)) + F
intermediate =unit_exponential(DistFix{bits+1, bits}, new_beta; reverse = reverse)
bit_vector =vcat([falsefor i in1:W-bits], intermediate.mantissa.number.bits[2:bits+1]...)
DistFix{W, F}(bit_vector) +DistFix{W, F}(start)
end""" laplace(::Type{DistFix{W, F}}, mean::Float64, scale::Float64, start::Float64, stop::Float64)Returns a Laplace distribution with the given mean and scale in the specified range [start, stop) of the specified type.For more information: https://en.wikipedia.org/wiki/Laplace_distribution"""functionlaplace(::Type{DistFix{W, F}}, mean::Float64, scale::Float64, start::Float64, stop::Float64) where {W, F}
@assert scale >0
coinflip =flip(0.5)
beta1 =-1/scale
e1 =exponential(DistFix{W, F}, beta1, mean, stop)
beta2 =1/scale
e2 =exponential(DistFix{W, F}, beta2, start, mean)
ifelse(coinflip, e1, e2)
end""" n_unit_exponentials(::Type{DistFix{W, F}}, betas::Vector{Float64})Returns n exponential distributions over the unit range with interleaved bits for a smaller BDD size."""functionn_unit_exponentials(::Type{DistFix{W, F}}, betas::Vector{Float64}) where {W, F}
DFiP = DistFix{W, F}
l =length(betas)
ans = [vcat([falsefor _ in1:W-F], Vector(undef, F)) for _ in1:l]
for i in F:-1:1for j in1:l
ans[j][i + W - F] =flip(exp(betas[j]/2^i)/(1+exp(betas[j]/2^i)))
endend
[DFiP(i) for i in ans]
end""" geometric(::Type{DistFix{W, F}}, success::Float64, stop::Int)Returns a geometric distribution of the given type with the given success parameter over integers in the range [0, stop).For more details: https://en.wikipedia.org/wiki/Geometric_distribution"""functiongeometric(::Type{DistFix{W, F}}, success::Float64, stop::Int) where {W, F}
@assertispow2(stop)
bits =Int(log2(stop))
@assert W - F > bits
convert(DistFix{W, F}, DistFix{W, 0}(unit_exponential(DistFix{bits+1, bits}, log(1- success)*2^bits).mantissa))
end""" unit_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64; vec_arg=[], constants = [], discrete_bdd=[])Returns a gamma distribution of the form x^α e^βx in the unit interval of the given type. The keyword arguments are there to control the order in which flips are created. """functionunit_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64; vec_arg=[], constants = [], discrete_bdd=[]) where {W, F}
DFiP = DistFix{W, F}
if alpha ==0unit_exponential(DFiP, beta)
elseif alpha ==1
t = (exp(beta*2.0^(-F))*(beta*2.0^(-F) -1) +1)*(1-exp(beta)) / ((1-exp(beta*2.0^(-F)))*(exp(beta) * (beta -1) +1))
coinflip =flip(t)
if (length(vec_arg) !=0)
(Y, Z, U) = vec_arg
else
(Y, Z, U) =n_unit_exponentials(DFiP, [beta, beta, 0.0])
endobserve(U < Y)
final =ifelse(coinflip, Z, Y)
final
else
α = alpha
β = beta
if (length(vec_arg) !=0)
vec_expo = vec_arg
else
discrete_bdd =Vector(undef, α)
constants =gamma_constants(alpha, beta, 1/2^F)
t = (exp(beta*2.0^(-F))*(beta*2.0^(-F) -1) +1)*(1-exp(beta)) / ((1-exp(beta*2.0^(-F)))*(exp(beta) * (beta -1) +1))
f =flip(t)
count =0for i in α:-1:1
l =discrete(DistUInt{max(Int(ceil(log(i))), 1)}, normalize(constants[count +2:count+i+1]))
count = count+i+1
discrete_bdd[α - i +1] = l
end
vec_expo =n_unit_exponentials(DFiP, exponential_for_gamma(α, β))
end
seq =Int(α*(α^2+5)/6)
x1 =unit_gamma(DFiP, alpha-1, beta, vec_arg=vec_expo[1:seq], constants=constants[α +2:length(constants)], discrete_bdd=discrete_bdd[2:α])
x2 = vec_expo[seq +1]
observe(x2 < x1)
discrete_dist_vec =Vector(undef, α)
count = seq+2for i in1:α
x = vec_expo[count]
count+=1for j in1:α - i
observe(vec_expo[count] < x)
count+=1end
discrete_dist_vec[i] = x
end# l = discrete(DistUInt{Int(ceil(log(α)))}, normalize(constants[2:α+1]))
l = discrete_bdd[1]
t =DFiP(0.0)
for i in1:α
t =ifelse(prob_equals(l, DistUInt{Int(ceil(log(α)))}(i-1)), discrete_dist_vec[i], t)
endifelse(flip(constants[1]), x1, t)
endend""" general_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64, ll::Float64, ul::Float64)Returns bitblast distribution for the density function x^α e^(β*x) in the range [ll, ul)"""functiongeneral_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64, ll::Float64, ul::Float64) where {W, F}
@assertispow2(ul - ll)
multiply =Int(log2(ul - ll))
start =DistFix{W, F}(ul)
new_type = DistFix{W, F + multiply}
DistFix{W, F}(unit_gamma(new_type, alpha, beta).mantissa.number.bits) + start
end##################################################### Helper functions for `unit_gamma`####################################################functionnormalize(v)
l =sum(v)
[i/l for i in v]
end"""Function to compute the mixing coefficients of distributions while constructing a `unit_gamma` distribution. α, β: x^α e^βx ϵ: 1/2^b depends on the type DistFix{W, F}"""functiongamma_constants(α::Int, β::Float64, ϵ::Float64)
@syms varint
@syms v2
if α ==0
[]
else
c1 =Float64(sympy.Poly(integrate(varint^α*exp(β*varint), (varint, 0, 1)), varint).coeffs().evalf()[1])
c2 = [Float64(i) for i in sympy.Poly(simplify(v2*integrate(varint^(α-1)*exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2).coeffs()]
p1 =0for i ineachindex(c2)
p1 +=sum_pgp(β, ϵ, length(c2) +1- i) * c2[i]
end
p1 /= c1
c2 = [Float64(i) for i in sympy.Poly(simplify(v2*integrate(varint^(α-1) * (varint - v2) *exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2).coeffs()]
p2 =Vector(undef, α)
for i ineachindex(c2)
p2[i] =sum_pgp(β, ϵ, length(c2) - i) * c2[i]
endvcat([p1], p2, gamma_constants(α-1, β, ϵ))
endend"""Function to compute sum of polynomial geometric series"""functionsum_pgp(β::Float64, ϵ::Float64, p::Int)
if p ==0sum_gp(β, ϵ)
elseif p ==1sum_agp(β, ϵ)
elseif p ==2sum_qgp(β, ϵ)
else
sum =0for i =0:ϵ:1-ϵ
sum += i^p *exp(β*i)
end
sum
endend"""Sum of quadratic geometric serieshttps://www.wolframalpha.com/input?i=sum+%28a*epsilon%29%5E2+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1"""functionsum_qgp(β::Float64, ϵ::Float64)
ans = (1/ϵ -1)^2*exp(β*ϵ*(2+1/ϵ))
ans += (1/ϵ^2)*exp(β)
ans += (2/ϵ -2/ϵ^2+1)*exp(β * (1+ ϵ))
ans -=exp(β*ϵ)
ans -=exp(2*β*ϵ)
ans *= ϵ^2
ans /= (exp(β * ϵ) -1)^3
ans
end"""Sum of arithmetic geometric progressionhttps://www.wolframalpha.com/input?i=sum+%28a*epsilon%29+*+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1"""functionsum_agp(β::Float64, ϵ::Float64)
ans = (1/ϵ -1)*exp(β*(1+ ϵ))
ans -=exp(β)/ϵ
ans +=exp(β*ϵ)
ans *= ϵ
ans /= (exp(β*ϵ) -1)^2
ans
end"""Sum of geometric progressionhttps://www.wolframalpha.com/input?i=sum+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1"""functionsum_gp(β::Float64, ϵ::Float64)
ans = (exp(β) -1) / (exp(β*ϵ) -1)
ans
end"""This function returns parameters of exponential distributions used in unit_gamma"""functionexponential_for_gamma(α::Int, β::Float64)::Vector{Float64}if α ==0
[]
elseif α ==1
[β, β, 0.0]
else
v = []
for i in1:α
v =vcat(vcat([β], zeros(i-1)), v)
endvcat(vcat(exponential_for_gamma(α-1, β), [0.0]), v)
endend############################################## bitblasting any general distribution#############################################""" bitblast(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, numpieces::Int, start::Float64, stop::Float64, blast_strategy=:linear; kwargs...)The function deploys the appropriate bitblasting function based on `blast_strategy`"""# The following is the function for individual piecesfunctionbitblast(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, numpieces::Int,
start::Float64, stop::Float64, blast_strategy=:linear; kwargs...) where {W,F}
if blast_strategy ==:linearbitblast_linear(DistFix{W,F}, dist, numpieces, start, stop; kwargs...)
elseif bitblast_strategy ==:exponentialbitblast_exponential(DistFix{W,F}, dist, numpieces, start, stop; kwargs...)
elseif bitblast_strategy ==:samplebitblast_sample(DistFix{W,F}, dist, numpieces, start, stop; kwargs...)
elseif bitblast_strategy ==:gammaerror("Not implemented yet")
elseerror("Unknown bitblasting strategy: $strategy")
endend""" bitblast_linear(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, numpieces::Int, start::Float64, stop::Float64)Returns a bitblasted representation of type `DistFix{W, F}` of the distribution `dist` using `numpieces` linear pieces in the range [start, stop) """functionbitblast_linear(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution,
numpieces::Int, start::Float64, stop::Float64) where {W,F}
# count bits and pieces@assert-(2^(W-F-1)) <= start < stop <=2^(W-F-1)
f_range_bits =log2((stop - start)*2^float(F))
@assertisinteger(f_range_bits) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two (not $f_range_bits)."@assertispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
intervals_per_piece = (2^Int(f_range_bits))/numpieces
bits_per_piece =Int(log2(intervals_per_piece))
# truncated distribution
dist =truncated(dist, start, stop)
# computing numbers to construct pieces and sew them together
total_prob =0
piece_probs =Vector(undef, numpieces) # prob of each piece
border_probs =Vector(undef, numpieces) # prob of first and last interval in each piece
linear_piece_probs =Vector(undef, numpieces) # prob of each piece if it were linear between end pointsfor i=1:numpieces
firstinter = start + (i-1)*intervals_per_piece/2^float(F)
lastinter = start + (i)*intervals_per_piece/2^float(F)
piece_probs[i] = (cdf(dist, lastinter) -cdf(dist, firstinter))
total_prob += piece_probs[i]
border_probs[i] = [cdf(dist, firstinter +1/2^float(F) ) -cdf(dist, firstinter),
cdf(dist, lastinter) -cdf(dist, lastinter -1/2^float(F) )]
linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2*2^(bits_per_piece)
end# coming up with discrete distribution to create a mixture of pieces
PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
piecechoice =discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces# computing flip probabilities for each flip
slope_flips =Vector(undef, numpieces)
isdecreasing =Vector(undef, numpieces)
for i=numpieces:-1:1iszero(linear_piece_probs[i]) &&iszero(piece_probs[i]) &&continue
a = border_probs[i][1]/linear_piece_probs[i]
isdecreasing[i] = a >1/2^bits_per_piece
if isdecreasing[i]
slope_flips[i] =flip(2-a*2^bits_per_piece)
else
slope_flips[i] =flip(a*2^bits_per_piece)
endend# building each piece with uniform and triangles
unif =uniform(DistFix{W,F}, bits_per_piece)
tria =triangle(DistFix{W,F}, bits_per_piece)
z =nothingfor i=1:numpieces
iszero(linear_piece_probs[i]) &&continue
firstinterval =DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^float(F))
lastinterval =DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^float(F))
linear_dist =if isdecreasing[i]
(ifelse(slope_flips[i],
(firstinterval + unif),
(lastinterval - tria)))
else
firstinterval +ifelse(slope_flips[i], unif, tria)
end
z =ifisnothing(z)
linear_dist
elseifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)
endendreturn z
end""" bitblast_exponential(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, numpieces::Int, start::Float64, stop::Float64)Returns a bitblasted representation of type `DistFix{W, F}` of the distribution `dist` using `numpieces` exponential pieces in the range [start, stop) """functionbitblast_exponential(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution,
numpieces::Int, start::Float64, stop::Float64) where {W,F}
# count bits and pieces@assert-(2^(W-F-1)) <= start < stop <=2^(W-F-1)
f_range_bits =log2((stop - start)*2^float(F))
@assertisinteger(f_range_bits) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two (not $f_range_bits)."@assertispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
intervals_per_piece = (2^Int(f_range_bits))/numpieces
# truncated distribution
d =truncated(dist, start, stop)
# computing numbers to construct pieces and sew them together
total_prob =0
piece_probs =Vector(undef, numpieces)
expo_beta =Vector(undef, numpieces)
for i=1:numpieces
firstinter = start + (i-1)*intervals_per_piece/2.0^(F)
lastinter = start + (i)*intervals_per_piece/2.0^(F)
expo_beta[i] =beta(d, firstinter, lastinter, 2.0^(-F))
piece_probs[i] = (cdf.(d, lastinter) -cdf.(d, firstinter))
total_prob += piece_probs[i]
end# coming up with discrete distribution to create a mixture of pieces
PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
piecechoice =discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces
z =nothingfor i=1:numpieces
expo_dist =exponential(DistFix{W, F}, expo_beta[i], start + (i-1)*intervals_per_piece/2^float(F), start + (i)*intervals_per_piece/2.0^(F))
z =ifisnothing(z)
expo_dist
elseifelse(prob_equals(piecechoice, PieceChoice(i-1)), expo_dist, z)
endendreturn z
end""" bitblast_sample(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, numpieces::Int, start::Float64, stop::Float64, offset::Float64, width::Float64)A modified version of the function bitblast that works with the assumption of lower bits being sampled."""functionbitblast_sample(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution,
numpieces::Int, start::Float64, stop::Float64; offset::Float64=0.0, width::Float64=1/2^float(F)) where {W,F}
# count bits and pieces@assert-(2^(W-F-1)) <= start < stop <=2^(W-F-1)
f_range_bits =log2((stop - start)*2^F)
@assertisinteger(f_range_bits) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two (not $f_range_bits)."@assertispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
intervals_per_piece = (2^Int(f_range_bits))/numpieces
bits_per_piece =Int(log2(intervals_per_piece))
# truncated distribution
dist =truncated(dist, start, stop)
# computing numbers to construct pieces and sew them together
total_prob =0
piece_probs =Vector(undef, numpieces) # prob of each piece
border_probs =Vector(undef, numpieces) # prob of first and last interval in each piece
linear_piece_probs =Vector(undef, numpieces) # prob of each piece if it were linear between end pointsfor i=1:numpieces
firstinter = start + (i-1)*intervals_per_piece/2^F
lastinter = start + (i)*intervals_per_piece/2^F
# Warning: A potential source of terrible runtime
piece_probs[i] =0for j=1:intervals_per_piece
piece_probs[i] +=cdf(dist, firstinter + offset + width + (j-1)/2^F) -cdf(dist, firstinter + offset + (j-1)/2^F)
end
total_prob += piece_probs[i]
border_probs[i] = [cdf(dist, firstinter + offset+width) -cdf(dist, firstinter + offset),
cdf(dist, lastinter -1/2^F + offset+width) -cdf(dist, lastinter -1/2^F + offset)]
linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2*2^(bits_per_piece)
end# coming up with discrete distribution to create a mixture of pieces
PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
piecechoice =discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces# computing flip probabilities for each piece
slope_flips =Vector(undef, numpieces)
isdecreasing =Vector(undef, numpieces)
for i=numpieces:-1:1iszero(linear_piece_probs[i]) &&iszero(piece_probs[i]) &&continue
a = border_probs[i][1]/linear_piece_probs[i]
isdecreasing[i] = a >1/2^bits_per_piece
if isdecreasing[i]
slope_flips[i] =flip(2-a*2^bits_per_piece)
else
slope_flips[i] =flip(a*2^bits_per_piece)
endend# building each piece with uniform and triangles
unif =uniform(DistFix{W,F}, bits_per_piece)
tria =triangle(DistFix{W,F}, bits_per_piece)
z =nothingfor i=1:numpieces
iszero(linear_piece_probs[i]) &&continue
firstinterval =DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^F)
lastinterval =DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^F)
linear_dist =if isdecreasing[i]
(ifelse(slope_flips[i],
(firstinterval + unif),
(lastinterval - tria)))
else
firstinterval +ifelse(slope_flips[i], unif, tria)
end
z =ifisnothing(z)
linear_dist
elseifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)
endendreturn z
end#################################### Helper functions for bitblasting###################################functionbeta(d::ContinuousUnivariateDistribution, start::Float64, stop::Float64, interval_sz::Float64)
prob_start =cdf(d, start + interval_sz) -cdf(d, start)
if prob_start ==0.0
prob_start =eps(0.0)
end
prob_end =cdf(d, stop) -cdf(d, stop - interval_sz)
result = (log(prob_end) -log(prob_start)) / (stop - start - interval_sz)
result
end#Helper function that returns exponentials functionshift_point_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64) where {W, F}
DFiP = DistFix{W, F}
if alpha ==0unit_exponential(DFiP, beta)
else
x1 =shift_point_gamma(DFiP, alpha -1, beta)
x2 =uniform(DFiP, 0.0, 1.0)
observe(ifelse(flip(1/(1+2.0^(-F))), x2 < x1, true))
x1
endend# # TODO: Write tests for the following function# function unit_concave(t::Type{DistFix{W, F}}, beta::Float64) where {W, F}# @assert beta <= 0
The text was updated successfully, but these errors were encountered:
border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
###################################################
############################################
border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
reverse=true refers to LSB to MSB order of flips
Dice.jl/src/dist/number/bitblast.jl
Line 55 in 4ef0d6e
The text was updated successfully, but these errors were encountered: