Coverage for /usr/local/lib/python3.11/dist-packages/grond/optimisers/highscore/optimiser.py: 87%
453 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-27 13:30 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-27 13:30 +0000
1from __future__ import print_function
2import math
3import os.path as op
4import os
5import logging
6import time
7import numpy as num
8from collections import OrderedDict
10from pyrocko.guts import StringChoice, Int, Float, Object, List
11from pyrocko.guts_array import Array
13from grond.meta import GrondError, Forbidden, has_get_plot_classes
14from grond.problems.base import ModelHistory
15from grond.optimisers.base import Optimiser, OptimiserConfig, BadProblem, \
16 OptimiserStatus
18guts_prefix = 'grond'
20logger = logging.getLogger('grond.optimisers.highscore.optimiser')
23def nextpow2(i):
24 return 2**int(math.ceil(math.log(i)/math.log(2.)))
27def excentricity_compensated_probabilities(xs, sbx, factor):
28 inonflat = num.where(sbx != 0.0)[0]
29 scale = num.zeros_like(sbx)
30 scale[inonflat] = 1.0 / (sbx[inonflat] * (factor if factor != 0. else 1.0))
31 distances_sqr_all = num.sum(
32 ((xs[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
33 scale[num.newaxis, num.newaxis, :])**2, axis=2)
34 probabilities = 1.0 / num.sum(distances_sqr_all < 1.0, axis=1)
35 # print(num.sort(num.sum(distances_sqr_all < 1.0, axis=1)))
36 probabilities /= num.sum(probabilities)
37 return probabilities
40def excentricity_compensated_choice(xs, sbx, factor, rstate):
41 probabilities = excentricity_compensated_probabilities(
42 xs, sbx, factor)
43 r = rstate.random_sample()
44 ichoice = num.searchsorted(num.cumsum(probabilities), r)
45 ichoice = min(ichoice, xs.shape[0]-1)
46 return ichoice
49def local_std(xs):
50 ssbx = num.sort(xs, axis=0)
51 dssbx = num.diff(ssbx, axis=0)
52 mdssbx = num.median(dssbx, axis=0)
53 return mdssbx * dssbx.shape[0] / 2.6
56class SamplerDistributionChoice(StringChoice):
57 choices = ['multivariate_normal', 'normal']
60class StandardDeviationEstimatorChoice(StringChoice):
61 choices = [
62 'median_density_single_chain',
63 'standard_deviation_all_chains',
64 'standard_deviation_single_chain']
67class SamplerStartingPointChoice(StringChoice):
68 choices = ['excentricity_compensated', 'random', 'mean']
71class BootstrapTypeChoice(StringChoice):
72 choices = ['bayesian', 'classic']
75def fnone(i):
76 return i if i is not None else -1
79class Sample(Object):
81 '''Sample model with context about how it was generated.'''
83 model = Array.T(shape=(None,), dtype=num.float64, serialize_as='list')
84 iphase = Int.T(optional=True)
85 ichain_base = Int.T(optional=True)
86 ilink_base = Int.T(optional=True)
87 imodel_base = Int.T(optional=True)
89 def preconstrain(self, problem):
90 self.model = problem.preconstrain(self.model)
92 def pack_context(self):
93 i = num.zeros(4, dtype=int)
94 i[:] = (
95 fnone(self.iphase),
96 fnone(self.ichain_base),
97 fnone(self.ilink_base),
98 fnone(self.imodel_base))
100 return i
103class SamplerPhase(Object):
104 niterations = Int.T(
105 help='Number of iteration for this phase.')
106 ntries_preconstrain_limit = Int.T(
107 default=1000,
108 help='Tries to find a valid preconstrained sample.')
109 seed = Int.T(
110 optional=True,
111 help='Random state seed.')
113 def __init__(self, *args, **kwargs):
114 Object.__init__(self, *args, **kwargs)
115 self._rstate = None
117 def get_rstate(self):
118 if self._rstate is None:
119 self._rstate = num.random.RandomState(self.seed)
121 return self._rstate
123 def get_raw_sample(self, problem, iiter, chains):
124 raise NotImplementedError
126 def get_sample(self, problem, iiter, chains):
127 assert 0 <= iiter < self.niterations
129 ntries_preconstrain = 0
130 for ntries_preconstrain in range(self.ntries_preconstrain_limit):
131 try:
132 sample = self.get_raw_sample(problem, iiter, chains)
133 sample.preconstrain(problem)
134 return sample
136 except Forbidden:
137 pass
139 raise GrondError(
140 'could not find any suitable candidate sample within %i tries' % (
141 self.ntries_preconstrain_limit))
144class InjectionSamplerPhase(SamplerPhase):
145 xs_inject = Array.T(
146 dtype=num.float64, shape=(None, None),
147 help='Array with the reference model.')
149 def get_raw_sample(self, problem, iiter, chains):
150 return Sample(model=self.xs_inject[iiter, :])
153class UniformSamplerPhase(SamplerPhase):
155 def get_raw_sample(self, problem, iiter, chains):
156 xbounds = problem.get_parameter_bounds()
157 return Sample(model=problem.random_uniform(xbounds, self.get_rstate()))
160class DirectedSamplerPhase(SamplerPhase):
161 scatter_scale = Float.T(
162 optional=True,
163 help='Scales search radius around the current `highscore` models')
164 scatter_scale_begin = Float.T(
165 optional=True,
166 help='Scaling factor at beginning of the phase.')
167 scatter_scale_end = Float.T(
168 optional=True,
169 help='Scaling factor at the end of the directed phase.')
170 starting_point = SamplerStartingPointChoice.T(
171 default='excentricity_compensated',
172 help='Tunes to the center value of the sampler distribution.'
173 'May increase the likelihood to draw a highscore member model'
174 ' off-center to the mean value')
176 sampler_distribution = SamplerDistributionChoice.T(
177 default='normal',
178 help='Distribution new models are drawn from.')
180 standard_deviation_estimator = StandardDeviationEstimatorChoice.T(
181 default='median_density_single_chain')
183 ntries_sample_limit = Int.T(default=1000)
185 def get_scatter_scale_factor(self, iiter):
186 s = self.scatter_scale
187 sa = self.scatter_scale_begin
188 sb = self.scatter_scale_end
190 assert s is None or (sa is None and sb is None)
192 if sa != sb:
193 tb = float(self.niterations-1)
194 tau = tb/(math.log(sa) - math.log(sb))
195 t0 = math.log(sa) * tau
196 t = float(iiter)
197 return num.exp(-(t-t0) / tau)
199 else:
200 return s or 1.0
202 def get_raw_sample(self, problem, iiter, chains):
203 rstate = self.get_rstate()
204 factor = self.get_scatter_scale_factor(iiter)
205 npar = problem.nparameters
206 pnames = problem.parameter_names
207 xbounds = problem.get_parameter_bounds()
209 ilink_choice = None
210 ichain_choice = num.argmin(chains.accept_sum)
212 if self.starting_point == 'excentricity_compensated':
213 models = chains.models(ichain_choice)
214 ilink_choice = excentricity_compensated_choice(
215 models,
216 chains.standard_deviation_models(
217 ichain_choice, self.standard_deviation_estimator),
218 2., rstate)
220 xchoice = chains.model(ichain_choice, ilink_choice)
222 elif self.starting_point == 'random':
223 ilink_choice = rstate.randint(0, chains.nlinks)
224 xchoice = chains.model(ichain_choice, ilink_choice)
226 elif self.starting_point == 'mean':
227 xchoice = chains.mean_model(ichain_choice)
229 else:
230 assert False, 'invalid starting_point choice: %s' % (
231 self.starting_point)
233 ntries_sample = 0
234 if self.sampler_distribution == 'normal':
235 x = num.zeros(npar, dtype=float)
236 sx = chains.standard_deviation_models(
237 ichain_choice, self.standard_deviation_estimator)
239 for ipar in range(npar):
240 ntries = 0
241 while True:
242 if xbounds[ipar, 0] == xbounds[ipar, 1]:
243 v = xbounds[ipar, 0]
244 break
246 if sx[ipar] > 0.:
247 v = rstate.normal(
248 xchoice[ipar],
249 factor*sx[ipar])
250 else:
251 v = xchoice[ipar]
253 if xbounds[ipar, 0] <= v and \
254 v <= xbounds[ipar, 1]:
256 break
258 if ntries > self.ntries_sample_limit:
259 logger.warning(
260 'failed to produce a suitable '
261 'candidate sample from normal '
262 'distribution for parameter \'%s\''
263 '- drawing from uniform instead.' %
264 pnames[ipar])
265 v = rstate.uniform(xbounds[ipar, 0],
266 xbounds[ipar, 1])
267 break
269 ntries += 1
271 x[ipar] = v
273 elif self.sampler_distribution == 'multivariate_normal':
274 ok_mask_sum = num.zeros(npar, dtype=int)
275 while True:
276 ntries_sample += 1
277 xcandi = rstate.multivariate_normal(
278 xchoice, factor**2 * chains.covariance_models(
279 ichain_choice))
281 ok_mask = num.logical_and(
282 xbounds[:, 0] <= xcandi, xcandi <= xbounds[:, 1])
284 if num.all(ok_mask):
285 break
287 ok_mask_sum += ok_mask
289 if ntries_sample > self.ntries_sample_limit:
290 logger.warning(
291 'failed to produce a suitable candidate '
292 'sample from multivariate normal '
293 'distribution, (%s) - drawing from uniform instead' %
294 ', '.join('%s:%i' % xx for xx in
295 zip(pnames, ok_mask_sum)))
296 xbounds = problem.get_parameter_bounds()
297 xcandi = problem.random_uniform(xbounds, rstate)
298 break
300 x = xcandi
302 imodel_base = None
303 if ilink_choice is not None:
304 imodel_base = chains.imodel(ichain_choice, ilink_choice)
306 return Sample(
307 model=x,
308 ichain_base=ichain_choice,
309 ilink_base=ilink_choice,
310 imodel_base=imodel_base)
313def make_bayesian_weights(nbootstrap, nmisfits,
314 type='bayesian', rstate=None):
315 ws = num.zeros((nbootstrap, nmisfits))
316 if rstate is None:
317 rstate = num.random.RandomState()
319 for ibootstrap in range(nbootstrap):
320 if type == 'classic':
321 ii = rstate.randint(0, nmisfits, size=nmisfits)
322 ws[ibootstrap, :] = num.histogram(
323 ii, nmisfits, (-0.5, nmisfits - 0.5))[0]
324 elif type == 'bayesian':
325 f = rstate.uniform(0., 1., size=nmisfits+1)
326 f[0] = 0.
327 f[-1] = 1.
328 f = num.sort(f)
329 g = f[1:] - f[:-1]
330 ws[ibootstrap, :] = g * nmisfits
331 else:
332 assert False
333 return ws
336class Chains(object):
337 def __init__(
338 self, problem, history, nchains, nlinks_cap):
340 self.problem = problem
341 self.history = history
342 self.nchains = nchains
343 self.nlinks_cap = nlinks_cap
344 self.chains_m = num.zeros(
345 (self.nchains, nlinks_cap), dtype=float)
346 self.chains_i = num.zeros(
347 (self.nchains, nlinks_cap), dtype=int)
348 self.nlinks = 0
349 self.nread = 0
351 self.accept_sum = num.zeros(self.nchains, dtype=int)
352 self._acceptance_history = num.zeros(
353 (self.nchains, 1024), dtype=bool)
355 history.add_listener(self)
357 def goto(self, n=None):
358 if n is None:
359 n = self.history.nmodels
361 n = min(self.history.nmodels, n)
363 assert self.nread <= n
365 while self.nread < n:
366 nread = self.nread
367 gbms = self.history.bootstrap_misfits[nread, :]
369 self.chains_m[:, self.nlinks] = gbms
370 self.chains_i[:, self.nlinks] = nread
371 nbootstrap = self.chains_m.shape[0]
373 self.nlinks += 1
374 chains_m = self.chains_m
375 chains_i = self.chains_i
377 for ichain in range(nbootstrap):
378 isort = num.argsort(chains_m[ichain, :self.nlinks])
379 chains_m[ichain, :self.nlinks] = chains_m[ichain, isort]
380 chains_i[ichain, :self.nlinks] = chains_i[ichain, isort]
382 if self.nlinks == self.nlinks_cap:
383 accept = (chains_i[:, self.nlinks_cap-1] != nread) \
384 .astype(bool)
385 self.nlinks -= 1
386 else:
387 accept = num.ones(self.nchains, dtype=bool)
389 self._append_acceptance(accept)
390 self.accept_sum += accept
391 self.nread += 1
393 def load(self):
394 return self.goto()
396 def extend(self, ioffset, n, models, misfits, sampler_contexts):
397 self.goto(ioffset + n)
399 def indices(self, ichain):
400 if ichain is not None:
401 return self.chains_i[ichain, :self.nlinks]
402 else:
403 return self.chains_i[:, :self.nlinks].ravel()
405 def models(self, ichain=None):
406 return self.history.models[self.indices(ichain), :]
408 def model(self, ichain, ilink):
409 return self.history.models[self.chains_i[ichain, ilink], :]
411 def imodel(self, ichain, ilink):
412 return self.chains_i[ichain, ilink]
414 def misfits(self, ichain=0):
415 return self.chains_m[ichain, :self.nlinks]
417 def misfit(self, ichain, ilink):
418 assert ilink < self.nlinks
419 return self.chains_m[ichain, ilink]
421 def mean_model(self, ichain=None):
422 xs = self.models(ichain)
423 return num.mean(xs, axis=0)
425 def best_model(self, ichain=0):
426 xs = self.models(ichain)
427 return xs[0]
429 def best_model_misfit(self, ichain=0):
430 return self.chains_m[ichain, 0]
432 def standard_deviation_models(self, ichain, estimator):
433 if estimator == 'median_density_single_chain':
434 xs = self.models(ichain)
435 return local_std(xs)
436 elif estimator == 'standard_deviation_all_chains':
437 bxs = self.models()
438 return num.std(bxs, axis=0)
439 elif estimator == 'standard_deviation_single_chain':
440 xs = self.models(ichain)
441 return num.std(xs, axis=0)
442 else:
443 assert False, 'invalid standard_deviation_estimator choice'
445 def covariance_models(self, ichain):
446 xs = self.models(ichain)
447 return num.cov(xs.T)
449 @property
450 def acceptance_history(self):
451 return self._acceptance_history[:, :self.nread]
453 def _append_acceptance(self, acceptance):
454 if self.nread >= self._acceptance_history.shape[1]:
455 new_buf = num.zeros(
456 (self.nchains, nextpow2(self.nread+1)), dtype=bool)
457 new_buf[:, :self._acceptance_history.shape[1]] = \
458 self._acceptance_history
459 self._acceptance_history = new_buf
460 self._acceptance_history[:, self.nread] = acceptance
463@has_get_plot_classes
464class HighScoreOptimiser(Optimiser):
465 '''Monte-Carlo-based directed search optimisation with bootstrap.'''
467 sampler_phases = List.T(SamplerPhase.T())
468 chain_length_factor = Float.T(default=8.)
469 nbootstrap = Int.T(default=100)
470 bootstrap_type = BootstrapTypeChoice.T(default='bayesian')
471 bootstrap_seed = Int.T(default=23)
473 SPARKS = u'\u2581\u2582\u2583\u2584\u2585\u2586\u2587\u2588'
474 ACCEPTANCE_AVG_LEN = 100
476 def __init__(self, **kwargs):
477 Optimiser.__init__(self, **kwargs)
478 self._bootstrap_weights = None
479 self._bootstrap_residuals = None
480 self._correlated_weights = None
481 self._status_chains = None
482 self._rstate_bootstrap = None
484 def get_rstate_bootstrap(self):
485 if self._rstate_bootstrap is None:
486 self._rstate_bootstrap = num.random.RandomState(
487 self.bootstrap_seed)
489 return self._rstate_bootstrap
491 def init_bootstraps(self, problem):
492 self.init_bootstrap_weights(problem)
493 self.init_bootstrap_residuals(problem)
495 def init_bootstrap_weights(self, problem):
496 logger.info('Initializing Bayesian bootstrap weights.')
498 nmisfits_w = sum(
499 t.nmisfits for t in problem.targets if t.can_bootstrap_weights)
501 ws = make_bayesian_weights(
502 self.nbootstrap,
503 nmisfits=nmisfits_w,
504 rstate=self.get_rstate_bootstrap())
506 imf = 0
507 for t in problem.targets:
508 if t.can_bootstrap_weights:
509 t.set_bootstrap_weights(ws[:, imf:imf+t.nmisfits])
510 imf += t.nmisfits
511 else:
512 t.set_bootstrap_weights(
513 num.ones((self.nbootstrap, t.nmisfits)))
515 def init_bootstrap_residuals(self, problem):
516 logger.info('Initializing Bayesian bootstrap residuals.')
518 for t in problem.targets:
519 if t.can_bootstrap_residuals:
520 t.init_bootstrap_residuals(
521 self.nbootstrap, rstate=self.get_rstate_bootstrap())
522 else:
523 t.set_bootstrap_residuals(
524 num.zeros((self.nbootstrap, t.nmisfits)))
526 def get_bootstrap_weights(self, problem):
527 if self._bootstrap_weights is None:
528 try:
529 problem.targets[0].get_bootstrap_weights()
530 except Exception:
531 self.init_bootstraps(problem)
533 bootstrap_weights = num.hstack(
534 [t.get_bootstrap_weights()
535 for t in problem.targets])
537 self._bootstrap_weights = num.vstack((
538 num.ones((1, problem.nmisfits)),
539 bootstrap_weights))
541 return self._bootstrap_weights
543 def get_bootstrap_residuals(self, problem):
544 if self._bootstrap_residuals is None:
545 try:
546 problem.targets[0].get_bootstrap_residuals()
547 except Exception:
548 self.init_bootstraps(problem)
550 bootstrap_residuals = num.hstack(
551 [t.get_bootstrap_residuals()
552 for t in problem.targets])
554 self._bootstrap_residuals = num.vstack((
555 num.zeros((1, problem.nmisfits)),
556 bootstrap_residuals))
558 return self._bootstrap_residuals
560 def get_correlated_weights(self, problem):
561 if self._correlated_weights is None:
562 corr = dict()
563 misfit_idx = num.cumsum(
564 [0.] + [t.nmisfits for t in problem.targets], dtype=int)
566 for it, target in enumerate(problem.targets):
567 weights = target.get_correlated_weights()
568 if weights is None:
569 continue
570 corr[misfit_idx[it]] = weights
572 self._correlated_weights = corr
574 return self._correlated_weights
576 @property
577 def nchains(self):
578 return self.nbootstrap + 1
580 def chains(self, problem, history):
581 nlinks_cap = int(round(
582 self.chain_length_factor * problem.nparameters + 1))
584 return Chains(
585 problem, history,
586 nchains=self.nchains, nlinks_cap=nlinks_cap)
588 def get_sampler_phase(self, iiter):
589 niter = 0
590 for iphase, phase in enumerate(self.sampler_phases):
591 if iiter < niter + phase.niterations:
592 return iphase, phase, iiter - niter
594 niter += phase.niterations
596 assert False, 'sample out of bounds'
598 def log_progress(self, problem, iiter, niter, phase, iiter_phase):
599 t = time.time()
600 if self._tlog_last < t - 10. \
601 or iiter_phase == 0 \
602 or iiter_phase == phase.niterations - 1:
604 logger.info(
605 '%s at %i/%i (%s, %i/%i)' % (
606 problem.name,
607 iiter+1, niter,
608 phase.__class__.__name__, iiter_phase, phase.niterations))
610 self._tlog_last = t
612 def optimise(self, problem, rundir=None):
613 if rundir is not None:
614 self.dump(filename=op.join(rundir, 'optimiser.yaml'))
616 history = ModelHistory(problem,
617 nchains=self.nchains,
618 path=rundir, mode='w')
619 chains = self.chains(problem, history)
621 niter = self.niterations
622 isbad_mask = None
623 self._tlog_last = 0
624 for iiter in range(niter):
625 iphase, phase, iiter_phase = self.get_sampler_phase(iiter)
626 self.log_progress(problem, iiter, niter, phase, iiter_phase)
628 sample = phase.get_sample(problem, iiter_phase, chains)
629 sample.iphase = iphase
631 if isbad_mask is not None and num.any(isbad_mask):
632 isok_mask = num.logical_not(isbad_mask)
633 else:
634 isok_mask = None
636 misfits = problem.misfits(sample.model, mask=isok_mask)
638 bootstrap_misfits = problem.combine_misfits(
639 misfits,
640 extra_weights=self.get_bootstrap_weights(problem),
641 extra_residuals=self.get_bootstrap_residuals(problem),
642 extra_correlated_weights=self.get_correlated_weights(problem))
643 isbad_mask_new = num.isnan(misfits[:, 0])
644 if isbad_mask is not None and num.any(
645 isbad_mask != isbad_mask_new):
647 errmess = [
648 'Problem "%s": inconsistency in data availability'
649 ' at iteration %i' %
650 (problem.name, iiter)]
652 for target, isbad_new, isbad in zip(
653 problem.targets, isbad_mask_new, isbad_mask):
655 if isbad_new != isbad:
656 errmess.append(' %s, %s -> %s' % (
657 target.string_id(), isbad, isbad_new))
659 raise BadProblem('\n'.join(errmess))
661 isbad_mask = isbad_mask_new
663 if num.all(isbad_mask):
664 raise BadProblem(
665 'Problem "%s": all target misfit values are NaN.'
666 % problem.name)
668 if num.all(num.isnan(bootstrap_misfits)):
669 raise BadProblem(
670 'Problem "%s": all misfits are NaN.' % problem.name)
672 history.append(
673 sample.model, misfits,
674 bootstrap_misfits,
675 sample.pack_context())
677 @property
678 def niterations(self):
679 return sum([ph.niterations for ph in self.sampler_phases])
681 def get_status(self, history):
682 if self._status_chains is None:
683 self._status_chains = self.chains(history.problem, history)
685 self._status_chains.goto(history.nmodels)
687 chains = self._status_chains
688 problem = history.problem
690 row_names = [p.name_nogroups for p in problem.parameters]
691 row_names.append('Misfit')
693 def colum_array(data):
694 arr = num.full(len(row_names), fill_value=num.nan)
695 arr[:data.size] = data
696 return arr
698 phase = self.get_sampler_phase(history.nmodels-1)[1]
700 bs_mean = colum_array(chains.mean_model(ichain=None))
701 bs_std = colum_array(chains.standard_deviation_models(
702 ichain=None, estimator='standard_deviation_all_chains'))
704 glob_mean = colum_array(chains.mean_model(ichain=0))
705 glob_mean[-1] = num.mean(chains.misfits(ichain=0))
707 glob_std = colum_array(chains.standard_deviation_models(
708 ichain=0, estimator='standard_deviation_single_chain'))
709 glob_std[-1] = num.std(chains.misfits(ichain=0))
711 glob_best = colum_array(chains.best_model(ichain=0))
712 glob_best[-1] = chains.best_model_misfit()
714 glob_misfits = chains.misfits(ichain=0)
716 acceptance_latest = chains.acceptance_history[
717 :, -min(chains.acceptance_history.shape[1], self.ACCEPTANCE_AVG_LEN):] # noqa
718 acceptance_avg = acceptance_latest.mean(axis=1)
720 def spark_plot(data, bins):
721 hist, _ = num.histogram(data, bins)
722 hist_max = num.max(hist)
723 if hist_max == 0.0:
724 hist_max = 1.0
725 hist = hist / hist_max
726 vec = num.digitize(hist, num.linspace(0., 1., len(self.SPARKS)))
727 return ''.join([self.SPARKS[b-1] for b in vec])
729 return OptimiserStatus(
730 row_names=row_names,
731 column_data=OrderedDict(
732 zip(['BS mean', 'BS std',
733 'Glob mean', 'Glob std', 'Glob best'],
734 [bs_mean, bs_std, glob_mean, glob_std, glob_best])),
735 extra_header= # noqa
736 u'Optimiser phase: {phase}, exploring {nchains} BS chains\n' # noqa
737 u'Global chain misfit distribution: \u2080{mf_dist}\xb9\n'
738 u'Acceptance rate distribution: \u2080{acceptance}'
739 u'\u2081\u2080\u2080\ufe6a (Median {acceptance_med:.1f}%)'
740 .format(
741 phase=phase.__class__.__name__,
742 nchains=chains.nchains,
743 mf_dist=spark_plot(
744 glob_misfits, num.linspace(0., 1., 25)),
745 acceptance=spark_plot(
746 acceptance_avg,
747 num.linspace(0., 1., 25)),
748 acceptance_med=num.median(acceptance_avg) * 100.
749 ))
751 def get_movie_maker(
752 self, problem, history, xpar_name, ypar_name, movie_filename):
754 from . import plot
755 return plot.HighScoreOptimiserPlot(
756 self, problem, history, xpar_name, ypar_name, movie_filename)
758 @classmethod
759 def get_plot_classes(cls):
760 from .plot import HighScoreAcceptancePlot
761 plots = Optimiser.get_plot_classes()
762 plots.append(HighScoreAcceptancePlot)
763 return plots
766class HighScoreOptimiserConfig(OptimiserConfig):
768 sampler_phases = List.T(
769 SamplerPhase.T(),
770 default=[UniformSamplerPhase(niterations=1000),
771 DirectedSamplerPhase(niterations=5000)],
772 help='Stages of the sampler: Start with uniform sampling of the model'
773 ' model space and narrow down through directed sampling.')
774 chain_length_factor = Float.T(
775 default=8.,
776 help='Controls the length of each chain: '
777 'chain_length_factor * nparameters + 1')
778 nbootstrap = Int.T(
779 default=100,
780 help='Number of bootstrap realisations to be tracked simultaneously in'
781 ' the optimisation.')
783 def get_optimiser(self):
784 return HighScoreOptimiser(
785 sampler_phases=list(self.sampler_phases),
786 chain_length_factor=self.chain_length_factor,
787 nbootstrap=self.nbootstrap)
790def load_optimiser_history(dirname, problem):
791 fn = op.join(dirname, 'accepted')
792 with open(fn, 'r') as f:
793 nmodels = os.fstat(f.fileno()).st_size // (problem.nbootstrap+1)
794 data1 = num.fromfile(
795 f,
796 dtype='<i1',
797 count=nmodels*(problem.nbootstrap+1)).astype(bool)
799 accepted = data1.reshape((nmodels, problem.nbootstrap+1))
801 fn = op.join(dirname, 'choices')
802 with open(fn, 'r') as f:
803 data2 = num.fromfile(
804 f,
805 dtype='<i8',
806 count=nmodels*2).astype(num.int64)
808 ibootstrap_choices, imodel_choices = data2.reshape((nmodels, 2)).T
809 return ibootstrap_choices, imodel_choices, accepted
812__all__ = '''
813 SamplerDistributionChoice
814 StandardDeviationEstimatorChoice
815 SamplerPhase
816 InjectionSamplerPhase
817 UniformSamplerPhase
818 DirectedSamplerPhase
819 Chains
820 HighScoreOptimiserConfig
821 HighScoreOptimiser
822'''.split()