Line data Source code
1 : /**
2 : Bojan Nikolic <bojan@bnikolic.co.uk>, <b.nikolic@mrao.cam.ac.uk>
3 :
4 : \file prior_sampler.hxx
5 : Renamed to prior_sampler.h
6 :
7 : */
8 :
9 : #include <functional> // was <boost/bind.hpp>
10 :
11 : #include "prior_sampler.h"
12 :
13 : #include "priors.h"
14 : #include "metro_propose.h"
15 : #include "bnmin_main.h"
16 : #include "mcpoint.h"
17 : #include "mcmonitor.h"
18 : #include "markovchain.h"
19 : #include "nestedsampler.h"
20 :
21 : namespace Minim
22 : {
23 :
24 : using namespace std::placeholders;
25 :
26 81 : CPriorSampler::CPriorSampler(PriorNLikelihood &ml,
27 81 : ModelDesc &md):
28 81 : md(md),
29 81 : ml(ml),
30 81 : mon(NULL)
31 : {
32 81 : }
33 :
34 81 : CPriorSampler::~CPriorSampler()
35 : {
36 81 : }
37 :
38 0 : CSPMetro::CSPMetro(PriorNLikelihood &ml,
39 : ModelDesc &md,
40 : const std::vector<double> &sigmas,
41 0 : unsigned seed):
42 : CPriorSampler(ml,md),
43 0 : prop(new MetroPropose(sigmas, seed))
44 : {
45 :
46 0 : }
47 :
48 0 : CSPMetro::~CSPMetro()
49 : {
50 :
51 0 : }
52 :
53 0 : double CSPMetro::advance(double L,
54 : size_t maxprop)
55 : {
56 0 : const size_t npars=prop->nPars();
57 :
58 0 : std::vector<double> cpoint(npars);
59 0 : md.get(cpoint);
60 :
61 0 : std::vector<double> propose(npars);
62 :
63 : // Current log prior value
64 0 : double clprior = -ml.pprob();
65 :
66 0 : double cllikel = -ml.llprob();
67 :
68 0 : for(size_t sn=0 ; sn<maxprop; ++sn)
69 : {
70 0 : propose=cpoint;
71 0 : prop->displace(propose);
72 0 : md.put(propose);
73 :
74 0 : const double propllikel= -ml.llprob();
75 :
76 0 : if (mon)
77 : {
78 0 : MCPoint p(propose);
79 0 : p.ll=propllikel;
80 0 : mon->propose(p);
81 0 : }
82 :
83 0 : if (propllikel > L)
84 : {
85 0 : const double proplprior = - ml.pprob();
86 :
87 0 : double ldiff = proplprior - clprior;
88 :
89 0 : if (ldiff > 0 or
90 0 : prop->raccept() < exp(ldiff) )
91 : {
92 0 : cpoint=propose;
93 0 : clprior=proplprior;
94 0 : cllikel=propllikel;
95 :
96 0 : if (mon)
97 : {
98 0 : MCPoint p(cpoint);
99 0 : p.ll=cllikel;
100 0 : mon->accept(p);
101 0 : }
102 :
103 : }
104 : }
105 : }
106 : // Store current point in the model
107 0 : md.put(cpoint);
108 :
109 0 : return cllikel;
110 :
111 0 : }
112 :
113 77873505 : double likelihood(ModelDesc &md,
114 : PriorNLikelihood &ml,
115 : const std::vector<double> &x)
116 : {
117 77873505 : md.put(x);
118 77873505 : return ml.llprob();
119 : }
120 :
121 77873505 : double prior(ModelDesc &md,
122 : PriorNLikelihood &ml,
123 : const std::vector<double> &x)
124 : {
125 77873505 : md.put(x);
126 77873505 : return ml.pprob();
127 : }
128 :
129 0 : CSPAdaptive::CSPAdaptive(PriorNLikelihood &ml,
130 : ModelDesc &md,
131 0 : const std::vector<double> &initSigmas):
132 : CPriorSampler(ml,md),
133 0 : sigmas(initSigmas)
134 : {
135 0 : MarkovChain::fx_t flkl=std::bind(likelihood,
136 0 : std::ref(md),
137 0 : std::ref(ml),
138 0 : _1);
139 :
140 0 : MarkovChain::fx_t fprior=std::bind(prior,
141 0 : std::ref(md),
142 0 : std::ref(ml),
143 0 : _1);
144 0 : std::vector<double> ic(sigmas.size());
145 0 : md.get(ic);
146 :
147 0 : c.reset(new InitPntChain(ic,
148 : flkl,
149 : fprior,
150 0 : constrPriorP));
151 0 : }
152 :
153 0 : CSPAdaptive::~CSPAdaptive()
154 : {
155 0 : }
156 :
157 0 : double CSPAdaptive::advance(double /*L*/,
158 : size_t maxprop)
159 : {
160 0 : const size_t qrtr=maxprop/4;
161 0 : const double sf=0.95;
162 0 : const double xf=1.0;
163 :
164 0 : for (size_t j=0; j<sigmas.size(); ++j)
165 0 : sigmas[j]*=xf;
166 :
167 0 : std::vector<double> ic(c->n);
168 0 : md.get(ic);
169 :
170 0 : c->reset(ic);
171 :
172 0 : bool accepted=false;
173 0 : for(size_t i=0; i<maxprop; ++i)
174 : {
175 0 : bool a=normProp(*c,
176 0 : sigmas);
177 0 : accepted=accepted or a;
178 :
179 0 : if (accepted==false and i>qrtr)
180 : {
181 0 : for (size_t j=0; j<sigmas.size(); ++j)
182 0 : sigmas[j]*=sf;
183 : }
184 : }
185 0 : md.put(c->gcx());
186 :
187 : // Note the convention is to return the negative value
188 0 : return -c->gcl();
189 0 : }
190 :
191 :
192 :
193 81 : CSRMSSS::CSRMSSS(PriorNLikelihood &ml,
194 : ModelDesc &md,
195 81 : const std::set<MCPoint> &ss):
196 : CPriorSampler(ml,md),
197 81 : ss(ss)
198 : {
199 81 : }
200 :
201 0 : CSRMSSS::CSRMSSS(PriorNLikelihood &ml,
202 0 : NestedS &s):
203 : CPriorSampler(ml,s),
204 0 : ss(s.g_ss())
205 : {
206 0 : }
207 :
208 162 : CSRMSSS::~CSRMSSS()
209 : {
210 162 : }
211 :
212 81 : void CSRMSSS::initChain(void)
213 : {
214 81 : MarkovChain::fx_t flkl=std::bind(likelihood,
215 81 : std::ref(md),
216 81 : std::ref(ml),
217 81 : _1);
218 :
219 81 : MarkovChain::fx_t fprior=std::bind(prior,
220 81 : std::ref(md),
221 81 : std::ref(ml),
222 81 : _1);
223 81 : std::vector<double> ic(ss.begin()->p.size());
224 81 : md.get(ic);
225 :
226 162 : c.reset(new ILklChain(ic,
227 : flkl,
228 : fprior,
229 81 : constrPriorL));
230 :
231 81 : nprop=0;
232 81 : }
233 :
234 771024 : double CSRMSSS::advance(double L,
235 : size_t maxprop)
236 : {
237 771024 : const double sf=0.1;
238 :
239 771024 : if (not c)
240 81 : initChain();
241 :
242 771024 : const size_t n=c->n;
243 :
244 :
245 771024 : std::vector<double> ic(n);
246 771024 : md.get(ic);
247 771024 : c->reset(ic,
248 : L);
249 :
250 771024 : std::vector<double> cv, eigvals, eigvects;
251 771024 : omoment2(ss, cv);
252 771024 : principalCV(cv,
253 : eigvals, eigvects);
254 :
255 3084096 : for(size_t j=0; j<eigvals.size(); ++j)
256 2313072 : eigvals[j]= pow(eigvals[j],0.5)*sf;
257 :
258 77873424 : for(size_t i=0; i<maxprop; ++i)
259 : {
260 77102400 : std::vector<double> sigmas(n,0);
261 77102400 : sigmas[nprop%n]=eigvals[nprop%n];
262 77102400 : eigenProp(*c,
263 : //sigmas,
264 : eigvals,
265 : eigvects);
266 77102400 : ++nprop;
267 77102400 : }
268 771024 : md.put(c->gcx());
269 :
270 : // Note the convention is to return the negative value
271 1542048 : return -c->gcl();
272 771024 : }
273 :
274 :
275 : }
276 :
277 :
|
