-
Notifications
You must be signed in to change notification settings - Fork 24
/
form_point_mass.jl
148 lines (113 loc) · 6.69 KB
/
form_point_mass.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
export PointMassFormConstraint
import ReactiveMP: default_form_check_strategy, default_prod_constraint, constrain_form
import DomainSets: Domain, infimum, supremum
import Optim
using BayesBase, Distributions, ExponentialFamily
"""
PointMassFormConstraint
One of the form constraint objects. Constraint a message to be in a form of dirac's delta point mass.
By default uses `Optim.jl` package to find argmin of `-logpdf(x)`.
Accepts custom `optimizer` callback which might be used to customise optimisation procedure with different packages
or different arguments for `Optim.jl` package.
# Keyword arguments
- `optimizer`: specifies a callback function for logpdf optimisation. See also: `RxInfer.default_point_mass_form_constraint_optimizer`
- `starting_point`: specifies a callback function for initial optimisation point: See also: `RxInfer.default_point_mass_form_constraint_starting_point`
- `boundaries`: specifies a callback function for determining optimisation boundaries: See also: `RxInfer.default_point_mass_form_constraint_boundaries`
## Custom optimizer callback interface
```julia
# This is an example of the `custom_optimizer` interface
function custom_optimizer(::Type{ Univariate }, ::Type{ Continuous }, constraint::PointMassFormConstraint, distribution)
# should return argmin of the -logpdf(distribution)
end
```
## Custom starting point callback interface
```julia
# This is an example of the `custom_starting_point` interface
function custom_starting_point(::Type{ Univariate }, ::Type{ Continuous }, constraint::PointMassFormConstraint, distribution)
# built-in optimizer expects an array, even for a univariate distribution
return [ 0.0 ]
end
```
## Custom boundaries callback interface
```julia
# This is an example of the `custom_boundaries` interface
function custom_boundaries(::Type{ Univariate }, ::Type{ Continuous }, constraint::PointMassFormConstraint, distribution)
# returns a tuple of `lower` and `upper` boundaries
return (-Inf, Inf)
end
```
"""
struct PointMassFormConstraint{F, P, B} <: AbstractFormConstraint
optimizer :: F
starting_point :: P
boundaries :: B
end
Base.show(io::IO, ::PointMassFormConstraint) = print(io, "PointMassFormConstraint()")
PointMassFormConstraint(;
optimizer = default_point_mass_form_constraint_optimizer,
starting_point = default_point_mass_form_constraint_starting_point,
boundaries = default_point_mass_form_constraint_boundaries
) = PointMassFormConstraint(optimizer, starting_point, boundaries)
ReactiveMP.default_form_check_strategy(::PointMassFormConstraint) = FormConstraintCheckLast()
ReactiveMP.default_prod_constraint(::PointMassFormConstraint) = GenericProd()
call_optimizer(pmconstraint::PointMassFormConstraint, distribution::D) where {D} = pmconstraint.optimizer(variate_form(D), value_support(D), pmconstraint, distribution)
call_boundaries(pmconstraint::PointMassFormConstraint, distribution::D) where {D} = pmconstraint.boundaries(variate_form(D), value_support(D), pmconstraint, distribution)
call_starting_point(pmconstraint::PointMassFormConstraint, distribution::D) where {D} = pmconstraint.starting_point(variate_form(D), value_support(D), pmconstraint, distribution)
ReactiveMP.constrain_form(pmconstraint::PointMassFormConstraint, distribution) = call_optimizer(pmconstraint, distribution)
# There is no need to call the optimizer on a `Distribution` object since they should have a well defined `mode`
ReactiveMP.constrain_form(::PointMassFormConstraint, distribution::Distribution) = PointMass(mode(distribution))
"""
default_point_mass_form_constraint_optimizer(::Type{<:VariateType}, ::Type{<:ValueSupport}, constraint::PointMassFormConstraint, distribution)
Defines a default optimisation procedure for the `PointMassFormConstraint`. By default uses `Optim.jl` package to find argmin of `-logpdf(x)`.
Uses the `starting_point` and `boundaries` callbacks to determine the starting point and boundaries for the optimisation procedure.
"""
function default_point_mass_form_constraint_optimizer end
function default_point_mass_form_constraint_optimizer(::Type{Univariate}, ::Type{Continuous}, constraint::PointMassFormConstraint, distribution)
target = let distribution = distribution
(x) -> -logpdf(distribution, x[1])
end
lower, upper = call_boundaries(constraint, distribution)
result = if isinf(lower) && isinf(upper)
Optim.optimize(target, call_starting_point(constraint, distribution), Optim.LBFGS())
else
Optim.optimize(target, [lower], [upper], call_starting_point(constraint, distribution), Optim.Fminbox(Optim.GradientDescent()))
end
if Optim.converged(result)
return PointMass(Optim.minimizer(result)[1])
else
error("Optimisation procedure for point mass estimation did not converge", result)
end
end
function default_point_mass_form_constraint_optimizer(::Type{Univariate}, ::Type{Discrete}, constraint::PointMassFormConstraint, distribution)
# fetch probvec
p = probvec(distribution)
# create new probvec
p_new = zeros(length(p))
p_new[argmax(p)] = 1
return PointMass(p_new)
end
"""
default_point_mass_form_constraint_boundaries(::Type{<:VariateType}, ::Type{<:ValueSupport}, constraint::PointMassFormConstraint, distribution)
Defines a default boundaries for the `PointMassFormConstraint`. By default simply uses the support of the distribution.
"""
function default_point_mass_form_constraint_boundaries end
function default_point_mass_form_constraint_boundaries(::Type{Univariate}, ::Type{Continuous}, constraint::PointMassFormConstraint, distribution)
return __default_univariate_boundaries(support(distribution))
end
__default_univariate_boundaries(interval::AbstractRange) = (minimum(interval), maximum(interval))
__default_univariate_boundaries(interval::Distributions.RealInterval) = (minimum(interval), maximum(interval))
__default_univariate_boundaries(domain::Domain) = (infimum(domain), supremum(domain))
"""
default_point_mass_form_constraint_starting_point(::Type{<:VariateType}, ::Type{<:ValueSupport}, constraint::PointMassFormConstraint, distribution)
Defines a default starting point for the `PointMassFormConstraint`. By default uses the support of the distribution.
If support is unbounded returns a zero point. Otherwise throws an error.
"""
function default_point_mass_form_constraint_starting_point end
function default_point_mass_form_constraint_starting_point(::Type{Univariate}, ::Type{Continuous}, constraint::PointMassFormConstraint, distribution)
lower, upper = call_boundaries(constraint, distribution)
return if isinf(lower) && isinf(upper)
return zeros(1)
else
error("No default starting point specified for a range [ $(lower), $(upper) ]")
end
end