[frames] | no frames]

# Source Code for Module gmisclib.mcmc_big

```  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)
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()
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
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()
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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Thu Sep 22 04:25:13 2011 http://epydoc.sourceforge.net