Coverage for /usr/local/lib/python3.11/dist-packages/grond/optimisers/highscore/optimiser.py: 87%

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 

9 

10from pyrocko.guts import StringChoice, Int, Float, Object, List 

11from pyrocko.guts_array import Array 

12 

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 

17 

18guts_prefix = 'grond' 

19 

20logger = logging.getLogger('grond.optimisers.highscore.optimiser') 

21 

22 

23def nextpow2(i): 

24 return 2**int(math.ceil(math.log(i)/math.log(2.))) 

25 

26 

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 

38 

39 

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 

47 

48 

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 

54 

55 

56class SamplerDistributionChoice(StringChoice): 

57 choices = ['multivariate_normal', 'normal'] 

58 

59 

60class StandardDeviationEstimatorChoice(StringChoice): 

61 choices = [ 

62 'median_density_single_chain', 

63 'standard_deviation_all_chains', 

64 'standard_deviation_single_chain'] 

65 

66 

67class SamplerStartingPointChoice(StringChoice): 

68 choices = ['excentricity_compensated', 'random', 'mean'] 

69 

70 

71class BootstrapTypeChoice(StringChoice): 

72 choices = ['bayesian', 'classic'] 

73 

74 

75def fnone(i): 

76 return i if i is not None else -1 

77 

78 

79class Sample(Object): 

80 

81 '''Sample model with context about how it was generated.''' 

82 

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) 

88 

89 def preconstrain(self, problem): 

90 self.model = problem.preconstrain(self.model) 

91 

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)) 

99 

100 return i 

101 

102 

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.') 

112 

113 def __init__(self, *args, **kwargs): 

114 Object.__init__(self, *args, **kwargs) 

115 self._rstate = None 

116 

117 def get_rstate(self): 

118 if self._rstate is None: 

119 self._rstate = num.random.RandomState(self.seed) 

120 

121 return self._rstate 

122 

123 def get_raw_sample(self, problem, iiter, chains): 

124 raise NotImplementedError 

125 

126 def get_sample(self, problem, iiter, chains): 

127 assert 0 <= iiter < self.niterations 

128 

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 

135 

136 except Forbidden: 

137 pass 

138 

139 raise GrondError( 

140 'could not find any suitable candidate sample within %i tries' % ( 

141 self.ntries_preconstrain_limit)) 

142 

143 

144class InjectionSamplerPhase(SamplerPhase): 

145 xs_inject = Array.T( 

146 dtype=num.float64, shape=(None, None), 

147 help='Array with the reference model.') 

148 

149 def get_raw_sample(self, problem, iiter, chains): 

150 return Sample(model=self.xs_inject[iiter, :]) 

151 

152 

153class UniformSamplerPhase(SamplerPhase): 

154 

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())) 

158 

159 

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') 

175 

176 sampler_distribution = SamplerDistributionChoice.T( 

177 default='normal', 

178 help='Distribution new models are drawn from.') 

179 

180 standard_deviation_estimator = StandardDeviationEstimatorChoice.T( 

181 default='median_density_single_chain') 

182 

183 ntries_sample_limit = Int.T(default=1000) 

184 

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 

189 

190 assert s is None or (sa is None and sb is None) 

191 

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) 

198 

199 else: 

200 return s or 1.0 

201 

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() 

208 

209 ilink_choice = None 

210 ichain_choice = num.argmin(chains.accept_sum) 

211 

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) 

219 

220 xchoice = chains.model(ichain_choice, ilink_choice) 

221 

222 elif self.starting_point == 'random': 

223 ilink_choice = rstate.randint(0, chains.nlinks) 

224 xchoice = chains.model(ichain_choice, ilink_choice) 

225 

226 elif self.starting_point == 'mean': 

227 xchoice = chains.mean_model(ichain_choice) 

228 

229 else: 

230 assert False, 'invalid starting_point choice: %s' % ( 

231 self.starting_point) 

232 

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) 

238 

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 

245 

246 if sx[ipar] > 0.: 

247 v = rstate.normal( 

248 xchoice[ipar], 

249 factor*sx[ipar]) 

250 else: 

251 v = xchoice[ipar] 

252 

253 if xbounds[ipar, 0] <= v and \ 

254 v <= xbounds[ipar, 1]: 

255 

256 break 

257 

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 

268 

269 ntries += 1 

270 

271 x[ipar] = v 

272 

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)) 

280 

281 ok_mask = num.logical_and( 

282 xbounds[:, 0] <= xcandi, xcandi <= xbounds[:, 1]) 

283 

284 if num.all(ok_mask): 

285 break 

286 

287 ok_mask_sum += ok_mask 

288 

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 

299 

300 x = xcandi 

301 

302 imodel_base = None 

303 if ilink_choice is not None: 

304 imodel_base = chains.imodel(ichain_choice, ilink_choice) 

305 

306 return Sample( 

307 model=x, 

308 ichain_base=ichain_choice, 

309 ilink_base=ilink_choice, 

310 imodel_base=imodel_base) 

311 

312 

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() 

318 

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 

334 

335 

336class Chains(object): 

337 def __init__( 

338 self, problem, history, nchains, nlinks_cap): 

339 

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 

350 

351 self.accept_sum = num.zeros(self.nchains, dtype=int) 

352 self._acceptance_history = num.zeros( 

353 (self.nchains, 1024), dtype=bool) 

354 

355 history.add_listener(self) 

356 

357 def goto(self, n=None): 

358 if n is None: 

359 n = self.history.nmodels 

360 

361 n = min(self.history.nmodels, n) 

362 

363 assert self.nread <= n 

364 

365 while self.nread < n: 

366 nread = self.nread 

367 gbms = self.history.bootstrap_misfits[nread, :] 

368 

369 self.chains_m[:, self.nlinks] = gbms 

370 self.chains_i[:, self.nlinks] = nread 

371 nbootstrap = self.chains_m.shape[0] 

372 

373 self.nlinks += 1 

374 chains_m = self.chains_m 

375 chains_i = self.chains_i 

376 

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] 

381 

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) 

388 

389 self._append_acceptance(accept) 

390 self.accept_sum += accept 

391 self.nread += 1 

392 

393 def load(self): 

394 return self.goto() 

395 

396 def extend(self, ioffset, n, models, misfits, sampler_contexts): 

397 self.goto(ioffset + n) 

398 

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() 

404 

405 def models(self, ichain=None): 

406 return self.history.models[self.indices(ichain), :] 

407 

408 def model(self, ichain, ilink): 

409 return self.history.models[self.chains_i[ichain, ilink], :] 

410 

411 def imodel(self, ichain, ilink): 

412 return self.chains_i[ichain, ilink] 

413 

414 def misfits(self, ichain=0): 

415 return self.chains_m[ichain, :self.nlinks] 

416 

417 def misfit(self, ichain, ilink): 

418 assert ilink < self.nlinks 

419 return self.chains_m[ichain, ilink] 

420 

421 def mean_model(self, ichain=None): 

422 xs = self.models(ichain) 

423 return num.mean(xs, axis=0) 

424 

425 def best_model(self, ichain=0): 

426 xs = self.models(ichain) 

427 return xs[0] 

428 

429 def best_model_misfit(self, ichain=0): 

430 return self.chains_m[ichain, 0] 

431 

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' 

444 

445 def covariance_models(self, ichain): 

446 xs = self.models(ichain) 

447 return num.cov(xs.T) 

448 

449 @property 

450 def acceptance_history(self): 

451 return self._acceptance_history[:, :self.nread] 

452 

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 

461 

462 

463@has_get_plot_classes 

464class HighScoreOptimiser(Optimiser): 

465 '''Monte-Carlo-based directed search optimisation with bootstrap.''' 

466 

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) 

472 

473 SPARKS = u'\u2581\u2582\u2583\u2584\u2585\u2586\u2587\u2588' 

474 ACCEPTANCE_AVG_LEN = 100 

475 

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 

483 

484 def get_rstate_bootstrap(self): 

485 if self._rstate_bootstrap is None: 

486 self._rstate_bootstrap = num.random.RandomState( 

487 self.bootstrap_seed) 

488 

489 return self._rstate_bootstrap 

490 

491 def init_bootstraps(self, problem): 

492 self.init_bootstrap_weights(problem) 

493 self.init_bootstrap_residuals(problem) 

494 

495 def init_bootstrap_weights(self, problem): 

496 logger.info('Initializing Bayesian bootstrap weights.') 

497 

498 nmisfits_w = sum( 

499 t.nmisfits for t in problem.targets if t.can_bootstrap_weights) 

500 

501 ws = make_bayesian_weights( 

502 self.nbootstrap, 

503 nmisfits=nmisfits_w, 

504 rstate=self.get_rstate_bootstrap()) 

505 

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))) 

514 

515 def init_bootstrap_residuals(self, problem): 

516 logger.info('Initializing Bayesian bootstrap residuals.') 

517 

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))) 

525 

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) 

532 

533 bootstrap_weights = num.hstack( 

534 [t.get_bootstrap_weights() 

535 for t in problem.targets]) 

536 

537 self._bootstrap_weights = num.vstack(( 

538 num.ones((1, problem.nmisfits)), 

539 bootstrap_weights)) 

540 

541 return self._bootstrap_weights 

542 

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) 

549 

550 bootstrap_residuals = num.hstack( 

551 [t.get_bootstrap_residuals() 

552 for t in problem.targets]) 

553 

554 self._bootstrap_residuals = num.vstack(( 

555 num.zeros((1, problem.nmisfits)), 

556 bootstrap_residuals)) 

557 

558 return self._bootstrap_residuals 

559 

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) 

565 

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 

571 

572 self._correlated_weights = corr 

573 

574 return self._correlated_weights 

575 

576 @property 

577 def nchains(self): 

578 return self.nbootstrap + 1 

579 

580 def chains(self, problem, history): 

581 nlinks_cap = int(round( 

582 self.chain_length_factor * problem.nparameters + 1)) 

583 

584 return Chains( 

585 problem, history, 

586 nchains=self.nchains, nlinks_cap=nlinks_cap) 

587 

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 

593 

594 niter += phase.niterations 

595 

596 assert False, 'sample out of bounds' 

597 

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: 

603 

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)) 

609 

610 self._tlog_last = t 

611 

612 def optimise(self, problem, rundir=None): 

613 if rundir is not None: 

614 self.dump(filename=op.join(rundir, 'optimiser.yaml')) 

615 

616 history = ModelHistory(problem, 

617 nchains=self.nchains, 

618 path=rundir, mode='w') 

619 chains = self.chains(problem, history) 

620 

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) 

627 

628 sample = phase.get_sample(problem, iiter_phase, chains) 

629 sample.iphase = iphase 

630 

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 

635 

636 misfits = problem.misfits(sample.model, mask=isok_mask) 

637 

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): 

646 

647 errmess = [ 

648 'Problem "%s": inconsistency in data availability' 

649 ' at iteration %i' % 

650 (problem.name, iiter)] 

651 

652 for target, isbad_new, isbad in zip( 

653 problem.targets, isbad_mask_new, isbad_mask): 

654 

655 if isbad_new != isbad: 

656 errmess.append(' %s, %s -> %s' % ( 

657 target.string_id(), isbad, isbad_new)) 

658 

659 raise BadProblem('\n'.join(errmess)) 

660 

661 isbad_mask = isbad_mask_new 

662 

663 if num.all(isbad_mask): 

664 raise BadProblem( 

665 'Problem "%s": all target misfit values are NaN.' 

666 % problem.name) 

667 

668 if num.all(num.isnan(bootstrap_misfits)): 

669 raise BadProblem( 

670 'Problem "%s": all misfits are NaN.' % problem.name) 

671 

672 history.append( 

673 sample.model, misfits, 

674 bootstrap_misfits, 

675 sample.pack_context()) 

676 

677 @property 

678 def niterations(self): 

679 return sum([ph.niterations for ph in self.sampler_phases]) 

680 

681 def get_status(self, history): 

682 if self._status_chains is None: 

683 self._status_chains = self.chains(history.problem, history) 

684 

685 self._status_chains.goto(history.nmodels) 

686 

687 chains = self._status_chains 

688 problem = history.problem 

689 

690 row_names = [p.name_nogroups for p in problem.parameters] 

691 row_names.append('Misfit') 

692 

693 def colum_array(data): 

694 arr = num.full(len(row_names), fill_value=num.nan) 

695 arr[:data.size] = data 

696 return arr 

697 

698 phase = self.get_sampler_phase(history.nmodels-1)[1] 

699 

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')) 

703 

704 glob_mean = colum_array(chains.mean_model(ichain=0)) 

705 glob_mean[-1] = num.mean(chains.misfits(ichain=0)) 

706 

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)) 

710 

711 glob_best = colum_array(chains.best_model(ichain=0)) 

712 glob_best[-1] = chains.best_model_misfit() 

713 

714 glob_misfits = chains.misfits(ichain=0) 

715 

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) 

719 

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]) 

728 

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 )) 

750 

751 def get_movie_maker(