gmisclib.mcmc

1 2 """An extension of mcmc that includes new stepping algorithms. 3 """ 4 5 from __future__ import with_statement 6 7 import math 8 import random 9 import numpy 10 11 import die 12 import mcmc 13 from mcmc import Debug 14 import gpk_lsq 15 16 SIGFAC = 3.0 17 18 NotGoodPosition = mcmc.NotGoodPosition 19

20 -def N_maximum(a):

21 return a[numpy.argmax(a)].item()

22 23

24 -def _parab_interp_guts(xylist):

25 ND = 3 26 assert len(xylist) >= 3 27 28 yy = numpy.zeros((len(xylist), 1)) 29 a = numpy.zeros((len(xylist), ND)) 30 for (i, (x, y)) in enumerate(xylist): 31 # Order of parameter is a Constant, linear slope, then quadratic. 32 yy.itemset(i, 0, y) 33 a.itemset(i, 0, 1.0) 34 a.itemset(i, 1, x) 35 a.itemset(i, 2, x*x) 36 37 lls = gpk_lsq.linear_least_squares(a, yy, copy=False) 38 39 try: 40 sv = lls.sv() 41 x = lls.x() 42 except (ValueError, numpy.linalg.linalg.LinAlgError), ex: 43 raise mcmc.NoBoot, "stepParab: %s in _compute_move" % str(ex) 44 if not (sv[-1] > 1e-6*sv[0]): 45 raise mcmc.NoBoot, "Sv ratio too extreme: %g/%g" % (sv[-1], sv[0]) 46 c = x.item(0,0) 47 b = x.item(1,0) 48 a = x.item(2,0) 49 # print 'c=', c 50 # print 'cc=', cc 51 return (a, b, c)

52 53

54 -def _parab_interp(points):

55 xlist = [x for (x, pos) in points] 56 tmp = [pos.logp_nocompute() for (x, pos) in points] 57 tmax = max(tmp) 58 ylist = [q-tmax for q in tmp] 59 # print 'XY=', (x0, y0), (x1, y1), (x2, y2) 60 61 # print 'xy=', (x0, y0), (x1, y1), (x2, y2) 62 a, b, c = _parab_interp_guts(zip(xlist, ylist)) 63 # print 'abc=', a, b, c 64 mn = min(xlist) 65 mx = max(xlist) 66 w = mx - mn 67 if a >= 0.0: 68 raise mcmc.NoBoot, "parab: positive curvature" 69 else: 70 # This looks like a nice parabola with a maximum. We step 71 # into the vicinity of the maximum. 72 xmin = -b/(2*a) 73 sigma = math.sqrt(-SIGFAC/(2*a)) 74 75 TOOCLOSE = 0.02 76 TOOFAR = 2.0 77 if xmin < -TOOFAR*w: 78 xmin = -TOOFAR*w 79 elif xmin > TOOFAR*w: 80 xmin = TOOFAR*w 81 if sigma > TOOFAR*w: 82 sigma = TOOFAR*w 83 elif sigma < 2*TOOCLOSE*w: 84 sigma = 2*TOOCLOSE*w 85 86 xtmp = random.normalvariate(xmin, sigma) 87 while min([abs(xtmp-x) for x in xlist]) < TOOCLOSE*w: 88 xtmp = random.normalvariate(xtmp, sigma) 89 # The minimum is too nearly a replication of an existing point, 90 # so we try again to make sure we get some new information. 91 return xtmp

92 93

94 -def test():

95 a, b, c = _parab_interp_guts([(0.0, 1.0), (1.0, 0.0), (2.0, 1.0)]) 96 assert abs(a-1.0) < 0.001 97 assert abs(-b/(2*a) - 1.0) < 0.001 98 assert abs(c-1.0) < 0.001

99 100 101

102 -def find_closest_p(v, vp):

103 """Searches a list of (v,p) pairs and finds the one whose v 104 is closest to the first argument. Returns 105 (v,p) of the closest pair. 106 """ 107 assert len(vp) > 0 108 vc, pc = vp[0] 109 for (vi, pi) in vp[1:]: 110 if abs(vi-v) < abs(vc-v): 111 vc = vi 112 pc = pi 113 return (vc, pc)

114 115 116

117 -def _fast_adjust(moveB, moveV, vfac, aV, aB, accepted):

118 """@note: this breaks the Markov assumption.""" 119 SLOP = 5.0 120 # If the moveB completely dominates moveV, then we can apply our knowledge of 121 # whether the step was accepted or not to adjusting the size of the Bootstrap step. 122 if numpy.greater(numpy.absolute(moveB), SLOP*numpy.absolute(moveV)).all(): 123 aB.inctry(accepted) 124 # Likewise if moveV dominates: 125 if numpy.greater(numpy.absolute(moveV), SLOP*numpy.absolute(moveB)).all(): 126 aV.inctry(accepted) 127 # Otherwise, we assume that the bootstrap steps are fairly successful, 128 # so the bootstrap size is approximately correct, and we 129 # adjust the stepV step size to approximately match the stepBoot step size. 130 if numpy.greater(numpy.absolute(moveB), SLOP*numpy.absolute(moveV)/vfac).all(): 131 aV.inctry(1) # moveV is too small: pretend stepV succeeded so it gets bigger. 132 if numpy.greater(numpy.absolute(moveV)/vfac, SLOP*numpy.absolute(moveB)).all(): 133 aV.inctry(0) # moveV is too big: pretend stepV failed so it gets smaller.

134 135 136

137 -def _has_enough(points, dP):

138 if len(points) < 3: 139 return False 140 if len(points) > 10: # Something is wrong. Stop wasting time on a parabolic approximation. 141 return True 142 lp = [pos.logp_nocompute() for (x, pos) in points] 143 spread = max(lp) - min(lp) 144 return spread > 3*dP/math.sqrt(len(points)-2)

145 146

147 -class BootStepper(mcmc.BootStepper):

148 - def __init__(self, lop, v, strategy=mcmc.BootStepper.SSAUTO, 149 maxArchSize=None, parallelSizeDiv=1):

150 mcmc.BootStepper.__init__(self, lop, v, strategy=strategy, 151 maxArchSize=maxArchSize, 152 parallelSizeDiv=parallelSizeDiv)

153 154

155 - def step(self):

156 mcmc.stepper.step(self) 157 Wboot = max(self.archive.distinct_count()-2, 0) 158 WV = self.np 159 Wmixed = math.sqrt(Wboot * WV) 160 # Parabolic steps are a bootstrap method, so they needs a large archive. 161 # But, Step_parab is not really Markov, so we only want to do it 162 if self.archive.sorted: 163 Wparab = max(self.archive.distinct_count()-2, 0) 164 else: 165 Wparab = 0 166 Wprobe = (Wboot+WV+Wmixed+1.5*Wparab) * self.F 167 W = float(Wboot + Wmixed + Wparab + WV + Wprobe) 168 Pboot = Wboot/W 169 Pmixed = Pboot + Wmixed/W 170 Pparab = Pmixed + Wparab/W 171 PV = Pparab + WV/W 172 173 again = True 174 while again: 175 P = random.random() 176 try: 177 if P < Pboot: 178 accepted = self.step_boot() 179 elif P < Pmixed: 180 accepted = self.step_mixed() 181 elif P < Pparab: 182 # Note! This violates MCMC requirements! 183 # The resulting probability distribution will not be correct 184 # when step_parab() is used! 185 accepted = self.step_parab() 186 elif P < PV: 187 accepted = self.stepV() 188 else: 189 accepted = self.step_probe() 190 except mcmc.NoBoot, x: 191 if Debug > 0: 192 die.info('NoBoot: %s' % str(x)) 193 again = True 194 else: 195 again = False 196 return accepted

197 198

199 - def step_mixed(self):

200 vfac = math.exp(random.random()*math.log(1e-5)) 201 self.steptype = 'step_mixed' 202 if len(self.archive) <= 2: 203 # print 'NoBoot', len(self.archive) 204 raise mcmc.NoBoot, "mixed short archive" 205 p1 = self.archive.choose() 206 p2 = self.archive.choose() 207 if p1.uid() == p2.uid(): 208 # print 'NoBoot dup' 209 raise mcmc.NoBoot, "mixed duplicate" 210 211 vsV = self.aV.vs() * vfac 212 if self.archive.distinct_count() >= min(5, self.np_eff): 213 vsV *= (self.archive.variance()/self.v.diagonal())**(1./4.) 214 215 moveB = (p1.vec() - p2.vec()) * self.aB.vs() 216 moveV = vsV * self.V.sample() 217 move = moveB + moveV 218 try: 219 tmp = self.current().new(move) 220 delta = tmp.logp() - self.current().logp() 221 except mcmc.NotGoodPosition: 222 die.warn('NotGoodPosition') 223 # Should I call self.aM.inctry(accepted) ? 224 return 0 225 if self.acceptable(delta): 226 accepted = 1 227 self._set_current(tmp) 228 if Debug>2: 229 die.info('StepMixed: Accepted logp=%g' % tmp.logp_nocompute()) 230 else: 231 self._set_failed( tmp ) 232 accepted = 0 233 self._set_current(self.current()) 234 if Debug>2: 235 die.info('StepMixed: Rejected logp=%g vs. %g, T=%g' 236 % (tmp.logp_nocompute(), self.current().logp_nocompute(), self.acceptable.T())) 237 238 # Only do this after a reset or when we don't care about breaking Markov assumptions: 239 if self.archive.sorted: 240 _fast_adjust(moveB, moveV, vfac, self.aV, self.aB, accepted) 241 return accepted

242 243 244

245 - def step_parab(self):

246 self.steptype = 'step_parab' 247 # The parabolic fit will be degenerate if we pick vs0 too close 248 # to either of the existing points. 249 if len(self.archive) <= 2: 250 raise mcmc.NoBoot, "parab short archive" 251 252 p1 = self.current() 253 p2 = self.archive.choose() 254 if p1.uid() == p2.uid(): 255 raise mcmc.NoBoot, "parab duplicate" 256 points = [(0.0, p1), (1.0, p2)] 257 pp = [] 258 vsB = self.aB.vs() 259 while True: 260 vs0 = random.normalvariate(0.0, vsB) 261 # Don't pick a point too close to the existing points. 262 tooclose = False 263 for (x, pos) in points: 264 if abs(vs0-x) < 0.5/len(points): 265 tooclose = True 266 break 267 if tooclose: 268 continue 269 270 vbase, pbase = find_closest_p(vs0, points) 271 move = (p2.vec() - p1.vec()) * (vs0-vbase) 272 # This is an approximation to a MCMC step, at least when the 273 # archive is large and well-annealed. So, we might as well add 274 # it to the archive. 275 try: 276 tmp = pbase.new(move) 277 pp.append( (tmp, tmp.logp()) ) 278 except mcmc.NotGoodPosition: 279 die.warn('NotGoodPosition2a at %.3f' % vs0) 280 return 0 281 282 points.append( (vs0, tmp) ) 283 if _has_enough(points, self.acceptable.jitter()): 284 break 285 vsB *= 1.2 286 if Debug > 0: 287 die.info("step_parab: %d points at %s" % (len(points), [q for (q,pos) in points])) 288 289 # Now, we take one of the (potentially many) new points and treat that as 290 # a MCMC step and add it to the archive. We don't want to add them all, 291 # because we might build up strong linear structures in the archive and 292 # make it more likely that the optimization gets stuck in a subspace. 293 tmp, tmp_logp = random.choice(pp) 294 try: 295 delta = tmp_logp - self.current().logp() 296 except mcmc.NotGoodPosition: 297 die.warn('NotGoodPosition2a at %.3f' % vs0) 298 return 0 299 if self.acceptable(delta): 300 if Debug > 2: 301 die.info('StepParab Accepted1 at %.3f, logp=%g' % (vs0, tmp_logp)) 302 self._set_current(tmp) 303 else: 304 self._set_failed( tmp ) 305 self._set_current(self.current()) 306 # Return from here? No. It still can be used to 307 # define the interpolation parabola, even if logP 308 # isn't particularly good. 309 310 try: 311 vsnew = _parab_interp(points) 312 except mcmc.NotGoodPosition: 313 die.warn('NotGoodPosition preparing for _parab_interp().') 314 return 0 315 vbase, pbase = find_closest_p(vsnew, points) 316 move = (p2.vec() - p1.vec()) * (vsnew-vbase) 317 318 try: 319 tmp = pbase.new(move) 320 delta = tmp.logp() - self.current().logp() 321 except mcmc.NotGoodPosition: 322 die.warn('NotGoodPosition') 323 accepted = 0 324 if self.acceptable(delta): 325 if Debug > 2: 326 die.info('StepParab Accepted2 at %.3f, logp=%g' % (vsnew, tmp.logp_nocompute())) 327 accepted = 1 328 self._set_current(tmp) 329 else: 330 self._set_failed( tmp ) 331 accepted = 0 332 self._set_current(self.current()) 333 if Debug>2: 334 die.info('StepParab: Rejected logp=%g vs. %g, T=%g' 335 % (tmp.logp_nocompute(), self.current().logp_nocompute(), self.acceptable.T())) 336 return accepted

337 338 339 340

341 -def bootstepper(logp, x, v, c=None, strategy=BootStepper.SSAUTO, fixer=None, repeatable=True):

342 """This is (essentially) another interface to the class constructor. 343 It's really there for backwards compatibility. 344 """ 345 pd = mcmc.problem_definition_F(logp_fcn=logp, c=c, fixer=fixer) 346 position_constructor = [mcmc.position_nonrepeatable, mcmc.position_repeatable][repeatable] 347 return BootStepper(mcmc.make_list_of_positions(x, position_constructor, pd), v, strategy=strategy)

348 349 diag_variance = mcmc.diag_variance 350 stepper = mcmc.stepper 351 problem_definition = mcmc.problem_definition 352 problem_definition_F = mcmc.problem_definition_F 353 problem_definition = mcmc.problem_definition 354 position_repeatable = mcmc.position_repeatable 355 position_nonrepeatable = mcmc.position_nonrepeatable 356

357 -def test2d(stepper):

358 import dictops 359 start = numpy.array([9, 1]) 360 c = numpy.array([1.0, 0.1]) 361 V = numpy.identity(len(c), numpy.float) 362 # V = numpy.array([[ 20.69626808, 20.6904984 ], [ 20.6904984, 20.69477235]]) 363 364 x = stepper(mcmc._logp1, start, V, c=c) 365 v = numpy.zeros((x.np, x.np)) 366 n_per_steptype = dictops.dict_of_averages() 367 for i in range(1000): 368 accepted = x.step() 369 n_per_steptype.add(x.steptype, accepted) 370 lpsum = 0.0 371 lps2 = 0.0 372 na = 0 373 psum = numpy.zeros((x.np,)) 374 N = 30000 375 nreset = 0 376 nsorted = 0 377 rid = x.reset_id() 378 for i in range(N): 379 accepted = x.step() 380 na += accepted 381 n_per_steptype.add(x.steptype, accepted) 382 nsorted += x.archive.sorted 383 if x.reset_id() != rid: 384 rid = x.reset_id() 385 nreset += 1 386 lp = x.current().logp() 387 lpsum += lp 388 lps2 += lp**2 389 p = x.current().vec() 390 numpy.add(psum, p, psum) 391 v += numpy.outer(p, p) 392 for steptype in n_per_steptype.keys(): 393 print 'Step %s has %.0f successes out of %.0f trials: %.2f' % ( 394 steptype, n_per_steptype.get_both(steptype)[0], 395 n_per_steptype.get_both(steptype)[1], 396 n_per_steptype.get_avg(steptype) 397 ) 398 # assert N*x.PBootLim**1.5 < nboot < N*x.PBootLim**0.7 399 assert nsorted < N//30 400 assert N//8 < na < N//2 401 assert nreset < 30 402 lpsum /= N 403 assert abs(lpsum+0.5*x.np) < 0.15 404 lps2 /= N-1 405 lpvar = lps2 - lpsum**2 406 assert abs(lpvar-0.5*x.np) < 1.0 407 numpy.divide(psum, N, psum) 408 assert numpy.alltrue(numpy.less(numpy.absolute(psum), 6.0/numpy.sqrt(c))) 409 numpy.divide(v, N-1, v) 410 for i in range(x.np): 411 for j in range(x.np): 412 if i == j: 413 # print '2*v[ii]*c=', 2*v[i,j]*c[i]**2 414 assert abs(math.log(2*v[i,j]*c[i]*c[j])) < 0.1 415 else: 416 assert abs(v[i,j]) < 20/(c[i]*c[j])

417 418 419 if __name__ == '__main__': 420 # mcmc.test(stepper=bootstepper) 421 Debug = 3 422 test() 423 mcmc.test_(stepper=bootstepper) 424

Source Code for Module gmisclib.mcmc_big