|
1 | | - |
2 | | -export For |
3 | | -using Random |
4 | | -import Base |
5 | | - |
6 | | -struct For{T, F, I} <: AbstractProductMeasure |
7 | | - f::F |
8 | | - inds::I |
9 | | - |
10 | | - @inline function For(f::F, inds::I) where {F,I<:Tuple} |
11 | | - T = Core.Compiler.return_type(f, Tuple{_eltype.(inds)...}) |
12 | | - new{T,instance_type(f),I}(f, inds) |
13 | | - end |
14 | | - |
15 | | - For{T,F,I}(f::F, inds::I) where {T,F,I} = new{T,F,I}(f,inds) |
16 | | -end |
17 | | - |
18 | | -# For(f, gen::Base.Generator) = ProductMeasure(Base.Generator(f ∘ gen.f, gen.iter)) |
19 | | - |
20 | | -@inline function logdensity_def(d::For{T,F,I}, x::AbstractVector{X}) where {X,T,F,I<:Tuple{<:AbstractVector}} |
21 | | - ℓ = zero(float(Core.Compiler.return_type(logdensity_def, Tuple{T,X}))) |
22 | | - @inbounds for j in eachindex(x) |
23 | | - ℓ += logdensity_def(d.f(j), x[j]) |
24 | | - end |
25 | | - ℓ |
26 | | -end |
27 | | - |
28 | | -function logdensity_def(d::For, x::AbstractVector) |
29 | | - sum(eachindex(x)) do i |
30 | | - @inbounds logdensity_def(d.f(getindex.(d.inds,i)...), x[i]) |
31 | | - end |
32 | | -end |
33 | | - |
34 | | -function logdensity_def(d::For{T,F,I}, x::AbstractArray{X}) where {T,F,I,X} |
35 | | - ℓ = zero(float(Core.Compiler.return_type(logdensity_def, Tuple{T,X}))) |
36 | | - |
37 | | - @inbounds for j in CartesianIndices(x) |
38 | | - i = (getindex(ind, j) for ind in d.inds) |
39 | | - ℓ += logdensity_def(d.f(i...), x[j]) |
40 | | - end |
41 | | - ℓ |
42 | | -end |
43 | | - |
44 | | -function logdensity_def(d::For{T,F,I}, x) where {N,T,F,I<:NTuple{N,<:Base.Generator}} |
45 | | - sum(zip(x, d.inds...)) do (xⱼ, dⱼ...) |
46 | | - logdensity_def(d.f(dⱼ...), xⱼ) |
47 | | - end |
48 | | -end |
49 | | - |
50 | | -function logdensity_def(d::For{T,F,I}, x::AbstractVector) where {N,T,F,I<:NTuple{N,<:Base.Generator}} |
51 | | - sum(zip(x, d.inds...)) do (xⱼ, dⱼ...) |
52 | | - logdensity_def(d.f(dⱼ...), xⱼ) |
53 | | - end |
54 | | -end |
55 | | - |
56 | | -function marginals(d::For{T,F,I}) where {T,F,I} |
57 | | - f(x...) = d.f(x...)::T |
58 | | - mappedarray(f, d.inds...) |
59 | | -end |
60 | | - |
61 | | -function marginals(d::For{T,F,I}) where {N,T,F,I<:NTuple{N,<:Base.Generator}} |
62 | | - f(x...) = d.f(x...)::T |
63 | | - Iterators.map(f, d.inds...) |
64 | | -end |
65 | | - |
66 | | -@inline function basemeasure(d::For{T,F,I}) where {T,F,I} |
67 | | - B = Core.Compiler.return_type(basemeasure, Tuple{T}) |
68 | | - _basemeasure(d, B, static(Base.issingletontype(B))) |
69 | | -end |
70 | | - |
71 | | -@inline function _basemeasure(d::For{T,F,I}, ::Type{B}, ::True) where {T,F,I,B} |
72 | | - instance(B) ^ axes(d.inds) |
73 | | -end |
74 | | - |
75 | | -@inline function _basemeasure(d::For{T,F,I}, ::Type{B}, ::False) where {T,F,I,B<:AbstractMeasure} |
76 | | - f = basekleisli(d.f) |
77 | | - _For(B, f, d.inds) |
78 | | -end |
79 | | - |
80 | | -@inline function _basemeasure(d::For{T,F,I}, ::Type{B}, ::False) where {T,F,I,B} |
81 | | - productmeasure(basemeasure.(marginals(d))) |
82 | | -end |
83 | | - |
84 | | -function _basemeasure(d::For{T,F,I}, ::Type{B}, ::True) where {N,T<:AbstractMeasure,F,I<:NTuple{N,<:Base.Generator},B} |
85 | | - return instance(B) ^ minimum(length, d.inds) |
86 | | -end |
87 | | - |
88 | | -function _basemeasure(d::For{T,F,I}, ::Type{B}, ::False) where {N,T<:AbstractMeasure,F,I<:NTuple{N,<:Base.Generator},B} |
89 | | - f = basekleisli(d.f) |
90 | | - _For(B, f, d.inds) |
91 | | -end |
92 | | - |
93 | | -function Pretty.tile(d::For{T}) where {T} |
94 | | - result = Pretty.literal("For{") |
95 | | - result *= Pretty.tile(T) |
96 | | - result *= Pretty.literal("}") |
97 | | - result *= Pretty.list_layout( |
98 | | - [ |
99 | | - Pretty.literal(func_string(d.f, Tuple{_eltype.(d.inds)...})), |
100 | | - Pretty.tile.(d.inds)... |
101 | | - ] |
102 | | - ) |
103 | | -end |
104 | | - |
105 | | -function _For(::Type{T}, f::F, inds::I) where {T,F,I} |
106 | | - For{T,F,I}(f,inds) |
107 | | -end |
108 | | - |
109 | | - |
110 | | -""" |
111 | | - For(f, base...) |
112 | | -
|
113 | | -`For` provides a convenient way to construct a `ProductMeasure`. There are |
114 | | -several options for the `base`. With Julia's `do` notation, this can look very |
115 | | -similar to a standard `for` loop, while maintaining semantics structure that's |
116 | | -easier to work with. |
117 | | -
|
118 | | ------------- |
119 | | -
|
120 | | -# `For(f, base::Int...)` |
121 | | -
|
122 | | -When one or several `Int` values are passed for `base`, the result is treated as |
123 | | -depending on `CartesianIndices(base)`. |
124 | | -
|
125 | | -``` |
126 | | -julia> For(3) do λ Exponential(λ) end |> marginals |
127 | | -3-element mappedarray(MeasureBase.var"#17#18"{var"#15#16"}(var"#15#16"()), ::CartesianIndices{1, Tuple{Base.OneTo{Int64}}}) with eltype Exponential{(:λ,), Tuple{Int64}}: |
128 | | - Exponential(λ = 1,) |
129 | | - Exponential(λ = 2,) |
130 | | - Exponential(λ = 3,) |
131 | | -``` |
132 | | -
|
133 | | -``` |
134 | | -julia> For(4,3) do μ,σ Normal(μ,σ) end |> marginals |
135 | | -4×3 mappedarray(MeasureBase.var"#17#18"{var"#11#12"}(var"#11#12"()), ::CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}) with eltype Normal{(:μ, :σ), Tuple{Int64, Int64}}: |
136 | | - Normal(μ = 1, σ = 1) Normal(μ = 1, σ = 2) Normal(μ = 1, σ = 3) |
137 | | - Normal(μ = 2, σ = 1) Normal(μ = 2, σ = 2) Normal(μ = 2, σ = 3) |
138 | | - Normal(μ = 3, σ = 1) Normal(μ = 3, σ = 2) Normal(μ = 3, σ = 3) |
139 | | - Normal(μ = 4, σ = 1) Normal(μ = 4, σ = 2) Normal(μ = 4, σ = 3) |
140 | | -``` |
141 | | -
|
142 | | -------- |
143 | | -
|
144 | | -# `For(f, base::AbstractArray...)`` |
145 | | -
|
146 | | -In this case, `base` behaves as if the arrays are `zip`ped together before |
147 | | -applying the map. |
148 | | -
|
149 | | -``` |
150 | | -julia> For(randn(3)) do x Exponential(x) end |> marginals |
151 | | -3-element mappedarray(x->Main.Exponential(x), ::Vector{Float64}) with eltype Exponential{(:λ,), Tuple{Float64}}: |
152 | | - Exponential(λ = -0.268256,) |
153 | | - Exponential(λ = 1.53044,) |
154 | | - Exponential(λ = -1.08839,) |
155 | | -``` |
156 | | -
|
157 | | -``` |
158 | | -julia> For(1:3, 1:3) do μ,σ Normal(μ,σ) end |> marginals |
159 | | -3-element mappedarray((:μ, :σ)->Main.Normal(μ, σ), ::UnitRange{Int64}, ::UnitRange{Int64}) with eltype Normal{(:μ, :σ), Tuple{Int64, Int64}}: |
160 | | - Normal(μ = 1, σ = 1) |
161 | | - Normal(μ = 2, σ = 2) |
162 | | - Normal(μ = 3, σ = 3) |
163 | | -``` |
164 | | -
|
165 | | ----- |
166 | | -
|
167 | | -# `For(f, base::Base.Generator)` |
168 | | -
|
169 | | -For `Generator`s, the function maps over the values of the generator: |
170 | | -
|
171 | | -``` |
172 | | -julia> For(eachrow(rand(4,2))) do x Normal(x[1], x[2]) end |> marginals |> collect |
173 | | -4-element Vector{Normal{(:μ, :σ), Tuple{Float64, Float64}}}: |
174 | | - Normal(μ = 0.255024, σ = 0.570142) |
175 | | - Normal(μ = 0.970706, σ = 0.0776745) |
176 | | - Normal(μ = 0.731491, σ = 0.505837) |
177 | | - Normal(μ = 0.563112, σ = 0.98307) |
178 | | -``` |
179 | | -
|
180 | | -""" |
181 | | -For(f, inds...) = For(f, inds) |
182 | | -For(f, n::Integer) = For(f, Base.OneTo(n)) |
183 | | -For(f, inds::Integer...) = For(i -> f(Tuple(i)...), Base.CartesianIndices(inds)) |
184 | | -# For(f, inds::Base.Generator) = productmeasure(mymap(f, inds)) |
185 | | - |
186 | | -function Random.rand!(rng::AbstractRNG, d::For, x) |
187 | | - mar = marginals(d) |
188 | | - @inbounds for (dⱼ, j) in zip(marginals(d), eachindex(x)) |
189 | | - x[j] = rand(rng,dⱼ) |
190 | | - end |
191 | | - return x |
192 | | -end |
193 | | - |
0 commit comments