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

1# https://pyrocko.org/grond - GPLv3 

2# 

3# The Grond Developers, 21st Century 

4import math 

5import os.path as op 

6import os 

7import logging 

8import time 

9import numpy as num 

10from collections import OrderedDict 

11 

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

13from pyrocko.guts_array import Array 

14 

15from grond.meta import GrondError, Forbidden, has_get_plot_classes 

16from grond.problems.base import ModelHistory 

17from grond.optimisers.base import Optimiser, OptimiserConfig, BadProblem, \ 

18 OptimiserStatus 

19 

20guts_prefix = 'grond' 

21 

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

23 

24 

25def nextpow2(i): 

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

27 

28 

29def excentricity_compensated_probabilities(xs, sbx, factor): 

30 inonflat = num.where(sbx != 0.0)[0] 

31 scale = num.zeros_like(sbx) 

32 scale[inonflat] = 1.0 / (sbx[inonflat] * (factor if factor != 0. else 1.0)) 

33 distances_sqr_all = num.sum( 

34 ((xs[num.newaxis, :, :] - xs[:, num.newaxis, :]) * 

35 scale[num.newaxis, num.newaxis, :])**2, axis=2) 

36 probabilities = 1.0 / num.sum(distances_sqr_all < 1.0, axis=1) 

37 # print(num.sort(num.sum(distances_sqr_all < 1.0, axis=1))) 

38 probabilities /= num.sum(probabilities) 

39 return probabilities 

40 

41 

42def excentricity_compensated_choice(xs, sbx, factor, rstate): 

43 probabilities = excentricity_compensated_probabilities( 

44 xs, sbx, factor) 

45 r = rstate.random_sample() 

46 ichoice = num.searchsorted(num.cumsum(probabilities), r) 

47 ichoice = min(ichoice, xs.shape[0]-1) 

48 return ichoice 

49 

50 

51def local_std(xs): 

52 ssbx = num.sort(xs, axis=0) 

53 dssbx = num.diff(ssbx, axis=0) 

54 mdssbx = num.median(dssbx, axis=0) 

55 return mdssbx * dssbx.shape[0] / 2.6 

56 

57 

58class SamplerDistributionChoice(StringChoice): 

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

60 

61 

62class StandardDeviationEstimatorChoice(StringChoice): 

63 choices = [ 

64 'median_density_single_chain', 

65 'standard_deviation_all_chains', 

66 'standard_deviation_single_chain'] 

67 

68 

69class SamplerStartingPointChoice(StringChoice): 

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

71 

72 

73class BootstrapTypeChoice(StringChoice): 

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

75 

76 

77def fnone(i): 

78 return i if i is not None else -1 

79 

80 

81class Sample(Object): 

82 

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

84 

85 model = Array.T(shape=(None,), dtype=num.float64, serialize_as='list') 

86 iphase = Int.T(optional=True) 

87 ichain_base = Int.T(optional=True) 

88 ilink_base = Int.T(optional=True) 

89 imodel_base = Int.T(optional=True) 

90 

91 def preconstrain(self, problem): 

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

93 

94 def pack_context(self): 

95 i = num.zeros(4, dtype=int) 

96 i[:] = ( 

97 fnone(self.iphase), 

98 fnone(self.ichain_base), 

99 fnone(self.ilink_base), 

100 fnone(self.imodel_base)) 

101 

102 return i 

103 

104 

105class SamplerPhase(Object): 

106 niterations = Int.T( 

107 help='Number of iteration for this phase.') 

108 ntries_preconstrain_limit = Int.T( 

109 default=1000, 

110 help='Tries to find a valid preconstrained sample.') 

111 seed = Int.T( 

112 optional=True, 

113 help='Random state seed.') 

114 

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

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

117 self._rstate = None 

118 

119 def get_rstate(self): 

120 if self._rstate is None: 

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

122 

123 return self._rstate 

124 

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

126 raise NotImplementedError 

127 

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

129 assert 0 <= iiter < self.niterations 

130 

131 ntries_preconstrain = 0 

132 for ntries_preconstrain in range(self.ntries_preconstrain_limit): 

133 try: 

134 sample = self.get_raw_sample(problem, iiter, chains) 

135 sample.preconstrain(problem) 

136 return sample 

137 

138 except Forbidden: 

139 pass 

140 

141 raise GrondError( 

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

143 self.ntries_preconstrain_limit)) 

144 

145 

146class InjectionSamplerPhase(SamplerPhase): 

147 xs_inject = Array.T( 

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

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

150 

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

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

153 

154 

155class UniformSamplerPhase(SamplerPhase): 

156 

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

158 xbounds = problem.get_parameter_bounds() 

159 return Sample(model=problem.random_uniform(xbounds, self.get_rstate())) 

160 

161 

162class DirectedSamplerPhase(SamplerPhase): 

163 scatter_scale = Float.T( 

164 optional=True, 

165 help='Scales search radius around the current `highscore` models') 

166 scatter_scale_begin = Float.T( 

167 optional=True, 

168 help='Scaling factor at beginning of the phase.') 

169 scatter_scale_end = Float.T( 

170 optional=True, 

171 help='Scaling factor at the end of the directed phase.') 

172 starting_point = SamplerStartingPointChoice.T( 

173 default='excentricity_compensated', 

174 help='Tunes to the center value of the sampler distribution.' 

175 'May increase the likelihood to draw a highscore member model' 

176 ' off-center to the mean value') 

177 

178 sampler_distribution = SamplerDistributionChoice.T( 

179 default='normal', 

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

181 

182 standard_deviation_estimator = StandardDeviationEstimatorChoice.T( 

183 default='median_density_single_chain') 

184 

185 ntries_sample_limit = Int.T(default=1000) 

186 

187 def get_scatter_scale_factor(self, iiter): 

188 s = self.scatter_scale 

189 sa = self.scatter_scale_begin 

190 sb = self.scatter_scale_end 

191 

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

193 

194 if sa != sb: 

195 tb = float(self.niterations-1) 

196 tau = tb/(math.log(sa) - math.log(sb)) 

197 t0 = math.log(sa) * tau 

198 t = float(iiter) 

199 return num.exp(-(t-t0) / tau) 

200 

201 else: 

202 return s or 1.0 

203 

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

205 rstate = self.get_rstate() 

206 factor = self.get_scatter_scale_factor(iiter) 

207 npar = problem.nparameters 

208 pnames = problem.parameter_names 

209 xbounds = problem.get_parameter_bounds() 

210 

211 ilink_choice = None 

212 ichain_choice = num.argmin(chains.accept_sum) 

213 

214 if self.starting_point == 'excentricity_compensated': 

215 models = chains.models(ichain_choice) 

216 ilink_choice = excentricity_compensated_choice( 

217 models, 

218 chains.standard_deviation_models( 

219 ichain_choice, self.standard_deviation_estimator), 

220 2., rstate) 

221 

222 xchoice = chains.model(ichain_choice, ilink_choice) 

223 

224 elif self.starting_point == 'random': 

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

226 xchoice = chains.model(ichain_choice, ilink_choice) 

227 

228 elif self.starting_point == 'mean': 

229 xchoice = chains.mean_model(ichain_choice) 

230 

231 else: 

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

233 self.starting_point) 

234 

235 ntries_sample = 0 

236 if self.sampler_distribution == 'normal': 

237 x = num.zeros(npar, dtype=float) 

238 sx = chains.standard_deviation_models( 

239 ichain_choice, self.standard_deviation_estimator) 

240 

241 for ipar in range(npar): 

242 ntries = 0 

243 while True: 

244 if xbounds[ipar, 0] == xbounds[ipar, 1]: 

245 v = xbounds[ipar, 0] 

246 break 

247 

248 if sx[ipar] > 0.: 

249 v = rstate.normal( 

250 xchoice[ipar], 

251 factor*sx[ipar]) 

252 else: 

253 v = xchoice[ipar] 

254 

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

256 v <= xbounds[ipar, 1]: 

257 

258 break 

259 

260 if ntries > self.ntries_sample_limit: 

261 logger.warning( 

262 'failed to produce a suitable ' 

263 'candidate sample from normal ' 

264 'distribution for parameter \'%s\'' 

265 '- drawing from uniform instead.' % 

266 pnames[ipar]) 

267 v = rstate.uniform(xbounds[ipar, 0], 

268 xbounds[ipar, 1]) 

269 break 

270 

271 ntries += 1 

272 

273 x[ipar] = v 

274 

275 elif self.sampler_distribution == 'multivariate_normal': 

276 ok_mask_sum = num.zeros(npar, dtype=int) 

277 while True: 

278 ntries_sample += 1 

279 xcandi = rstate.multivariate_normal( 

280 xchoice, factor**2 * chains.covariance_models( 

281 ichain_choice)) 

282 

283 ok_mask = num.logical_and( 

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

285 

286 if num.all(ok_mask): 

287 break 

288 

289 ok_mask_sum += ok_mask 

290 

291 if ntries_sample > self.ntries_sample_limit: 

292 logger.warning( 

293 'failed to produce a suitable candidate ' 

294 'sample from multivariate normal ' 

295 'distribution, (%s) - drawing from uniform instead' % 

296 ', '.join('%s:%i' % xx for xx in 

297 zip(pnames, ok_mask_sum))) 

298 xbounds = problem.get_parameter_bounds() 

299 xcandi = problem.random_uniform(xbounds, rstate) 

300 break 

301 

302 x = xcandi 

303 

304 imodel_base = None 

305 if ilink_choice is not None: 

306 imodel_base = chains.imodel(ichain_choice, ilink_choice) 

307 

308 return Sample( 

309 model=x, 

310 ichain_base=ichain_choice, 

311 ilink_base=ilink_choice, 

312 imodel_base=imodel_base) 

313 

314 

315def make_bayesian_weights(nbootstrap, nmisfits, 

316 type='bayesian', rstate=None): 

317 ws = num.zeros((nbootstrap, nmisfits)) 

318 if rstate is None: 

319 rstate = num.random.RandomState() 

320 

321 for ibootstrap in range(nbootstrap): 

322 if type == 'classic': 

323 ii = rstate.randint(0, nmisfits, size=nmisfits) 

324 ws[ibootstrap, :] = num.histogram( 

325 ii, nmisfits, (-0.5, nmisfits - 0.5))[0] 

326 elif type == 'bayesian': 

327 f = rstate.uniform(0., 1., size=nmisfits+1) 

328 f[0] = 0. 

329 f[-1] = 1. 

330 f = num.sort(f) 

331 g = f[1:] - f[:-1] 

332 ws[ibootstrap, :] = g * nmisfits 

333 else: 

334 assert False 

335 return ws 

336 

337 

338class Chains(object): 

339 def __init__( 

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

341 

342 self.problem = problem 

343 self.history = history 

344 self.nchains = nchains 

345 self.nlinks_cap = nlinks_cap 

346 self.chains_m = num.zeros( 

347 (self.nchains, nlinks_cap), dtype=float) 

348 self.chains_i = num.zeros( 

349 (self.nchains, nlinks_cap), dtype=int) 

350 self.nlinks = 0 

351 self.nread = 0 

352 

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

354 self._acceptance_history = num.zeros( 

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

356 

357 history.add_listener(self) 

358 

359 def goto(self, n=None): 

360 if n is None: 

361 n = self.history.nmodels 

362 

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

364 

365 assert self.nread <= n 

366 

367 while self.nread < n: 

368 nread = self.nread 

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

370 

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

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

373 nbootstrap = self.chains_m.shape[0] 

374 

375 self.nlinks += 1 

376 chains_m = self.chains_m 

377 chains_i = self.chains_i 

378 

379 for ichain in range(nbootstrap): 

380 isort = num.argsort(chains_m[ichain, :self.nlinks]) 

381 chains_m[ichain, :self.nlinks] = chains_m[ichain, isort] 

382 chains_i[ichain, :self.nlinks] = chains_i[ichain, isort] 

383 

384 if self.nlinks == self.nlinks_cap: 

385 accept = (chains_i[:, self.nlinks_cap-1] != nread) \ 

386 .astype(bool) 

387 self.nlinks -= 1 

388 else: 

389 accept = num.ones(self.nchains, dtype=bool) 

390 

391 self._append_acceptance(accept) 

392 self.accept_sum += accept 

393 self.nread += 1 

394 

395 def load(self): 

396 return self.goto() 

397 

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

399 self.goto(ioffset + n) 

400 

401 def indices(self, ichain): 

402 if ichain is not None: 

403 return self.chains_i[ichain, :self.nlinks] 

404 else: 

405 return self.chains_i[:, :self.nlinks].ravel() 

406 

407 def models(self, ichain=None): 

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

409 

410 def model(self, ichain, ilink): 

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

412 

413 def imodel(self, ichain, ilink): 

414 return self.chains_i[ichain, ilink] 

415 

416 def misfits(self, ichain=0): 

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

418 

419 def misfit(self, ichain, ilink): 

420 assert ilink < self.nlinks 

421 return self.chains_m[ichain, ilink] 

422 

423 def mean_model(self, ichain=None): 

424 xs = self.models(ichain) 

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

426 

427 def best_model(self, ichain=0): 

428 xs = self.models(ichain) 

429 return xs[0] 

430 

431 def best_model_misfit(self, ichain=0): 

432 return self.chains_m[ichain, 0] 

433 

434 def standard_deviation_models(self, ichain, estimator): 

435 if estimator == 'median_density_single_chain': 

436 xs = self.models(ichain) 

437 return local_std(xs) 

438 elif estimator == 'standard_deviation_all_chains': 

439 bxs = self.models() 

440 return num.std(bxs, axis=0) 

441 elif estimator == 'standard_deviation_single_chain': 

442 xs = self.models(ichain) 

443 return num.std(xs, axis=0) 

444 else: 

445 assert False, 'invalid standard_deviation_estimator choice' 

446 

447 def covariance_models(self, ichain): 

448 xs = self.models(ichain) 

449 return num.cov(xs.T) 

450 

451 @property 

452 def acceptance_history(self): 

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

454 

455 def _append_acceptance(self, acceptance): 

456 if self.nread >= self._acceptance_history.shape[1]: 

457 new_buf = num.zeros( 

458 (self.nchains, nextpow2(self.nread+1)), dtype=bool) 

459 new_buf[:, :self._acceptance_history.shape[1]] = \ 

460 self._acceptance_history 

461 self._acceptance_history = new_buf 

462 self._acceptance_history[:, self.nread] = acceptance 

463 

464 

465@has_get_plot_classes 

466class HighScoreOptimiser(Optimiser): 

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

468 

469 sampler_phases = List.T(SamplerPhase.T()) 

470 chain_length_factor = Float.T(default=8.) 

471 nbootstrap = Int.T(default=100) 

472 bootstrap_type = BootstrapTypeChoice.T(default='bayesian') 

473 bootstrap_seed = Int.T(default=23) 

474 

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

476 ACCEPTANCE_AVG_LEN = 100 

477 

478 def __init__(self, **kwargs): 

479 Optimiser.__init__(self, **kwargs) 

480 self._bootstrap_weights = None 

481 self._bootstrap_residuals = None 

482 self._correlated_weights = None 

483 self._status_chains = None 

484 self._rstate_bootstrap = None 

485 

486 def get_rstate_bootstrap(self): 

487 if self._rstate_bootstrap is None: 

488 self._rstate_bootstrap = num.random.RandomState( 

489 self.bootstrap_seed) 

490 

491 return self._rstate_bootstrap 

492 

493 def init_bootstraps(self, problem): 

494 self.init_bootstrap_weights(problem) 

495 self.init_bootstrap_residuals(problem) 

496 

497 def init_bootstrap_weights(self, problem): 

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

499 

500 nmisfits_w = sum( 

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

502 

503 ws = make_bayesian_weights( 

504 self.nbootstrap, 

505 nmisfits=nmisfits_w, 

506 rstate=self.get_rstate_bootstrap()) 

507 

508 imf = 0 

509 for t in problem.targets: 

510 if t.can_bootstrap_weights: 

511 t.set_bootstrap_weights(ws[:, imf:imf+t.nmisfits]) 

512 imf += t.nmisfits 

513 else: 

514 t.set_bootstrap_weights( 

515 num.ones((self.nbootstrap, t.nmisfits))) 

516 

517 def init_bootstrap_residuals(self, problem): 

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

519 

520 for t in problem.targets: 

521 if t.can_bootstrap_residuals: 

522 t.init_bootstrap_residuals( 

523 self.nbootstrap, rstate=self.get_rstate_bootstrap()) 

524 else: 

525 t.set_bootstrap_residuals( 

526 num.zeros((self.nbootstrap, t.nmisfits))) 

527 

528 def get_bootstrap_weights(self, problem): 

529 if self._bootstrap_weights is None: 

530 try: 

531 problem.targets[0].get_bootstrap_weights() 

532 except Exception: 

533 self.init_bootstraps(problem) 

534 

535 bootstrap_weights = num.hstack( 

536 [t.get_bootstrap_weights() 

537 for t in problem.targets]) 

538 

539 self._bootstrap_weights = num.vstack(( 

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

541 bootstrap_weights)) 

542 

543 return self._bootstrap_weights 

544 

545 def get_bootstrap_residuals(self, problem): 

546 if self._bootstrap_residuals is None: 

547 try: 

548 problem.targets[0].get_bootstrap_residuals() 

549 except Exception: 

550 self.init_bootstraps(problem) 

551 

552 bootstrap_residuals = num.hstack( 

553 [t.get_bootstrap_residuals() 

554 for t in problem.targets]) 

555 

556 self._bootstrap_residuals = num.vstack(( 

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

558 bootstrap_residuals)) 

559 

560 return self._bootstrap_residuals 

561 

562 def get_correlated_weights(self, problem): 

563 if self._correlated_weights is None: 

564 corr = dict() 

565 misfit_idx = num.cumsum( 

566 [0.] + [t.nmisfits for t in problem.targets], dtype=int) 

567 

568 for it, target in enumerate(problem.targets): 

569 weights = target.get_correlated_weights() 

570 if weights is None: 

571 continue 

572 corr[misfit_idx[it]] = weights 

573 

574 self._correlated_weights = corr 

575 

576 return self._correlated_weights 

577 

578 @property 

579 def nchains(self): 

580 return self.nbootstrap + 1 

581 

582 def chains(self, problem, history): 

583 nlinks_cap = int(round( 

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

585 

586 return Chains( 

587 problem, history, 

588 nchains=self.nchains, nlinks_cap=nlinks_cap) 

589 

590 def get_sampler_phase(self, iiter): 

591 niter = 0 

592 for iphase, phase in enumerate(self.sampler_phases): 

593 if iiter < niter + phase.niterations: 

594 return iphase, phase, iiter - niter 

595 

596 niter += phase.niterations 

597 

598 assert False, 'sample out of bounds' 

599 

600 def log_progress(self, problem, iiter, niter, phase, iiter_phase): 

601 t = time.time() 

602 if self._tlog_last < t - 10. \ 

603 or iiter_phase == 0 \ 

604 or iiter_phase == phase.niterations - 1: 

605 

606 logger.info( 

607 '%s at %i/%i (%s, %i/%i)' % ( 

608 problem.name, 

609 iiter+1, niter, 

610 phase.__class__.__name__, iiter_phase, phase.niterations)) 

611 

612 self._tlog_last = t 

613 

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

615 if rundir is not None: 

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

617 

618 history = ModelHistory(problem, 

619 nchains=self.nchains, 

620 path=rundir, mode='w') 

621 chains = self.chains(problem, history) 

622 

623 niter = self.niterations 

624 isbad_mask = None 

625 self._tlog_last = 0 

626 for iiter in range(niter): 

627 iphase, phase, iiter_phase = self.get_sampler_phase(iiter) 

628 self.log_progress(problem, iiter, niter, phase, iiter_phase) 

629 

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

631 sample.iphase = iphase 

632 

633 if isbad_mask is not None and num.any(isbad_mask): 

634 isok_mask = num.logical_not(isbad_mask) 

635 else: 

636 isok_mask = None 

637 

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

639 

640 bootstrap_misfits = problem.combine_misfits( 

641 misfits, 

642 extra_weights=self.get_bootstrap_weights(problem), 

643 extra_residuals=self.get_bootstrap_residuals(problem), 

644 extra_correlated_weights=self.get_correlated_weights(problem)) 

645 isbad_mask_new = num.isnan(misfits[:, 0]) 

646 if isbad_mask is not None and num.any( 

647 isbad_mask != isbad_mask_new): 

648 

649 errmess = [ 

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

651 ' at iteration %i' % 

652 (problem.name, iiter)] 

653 

654 imisfit = 0 

655 for target in problem.targets: 

656 for itarget_misfit in range(target.nmisfits): 

657 i = imisfit + itarget_misfit 

658 isbad = isbad_mask[i] 

659 isbad_new = isbad_mask_new[i] 

660 if isbad_new != isbad: 

661 errmess.append(' %s, %i, %s -> %s' % ( 

662 target.string_id(), itarget_misfit, 

663 isbad, isbad_new)) 

664 

665 imisfit += target.nmisfits 

666 

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

668 

669 isbad_mask = isbad_mask_new 

670 

671 if num.all(isbad_mask): 

672 raise BadProblem( 

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

674 % problem.name) 

675 

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

677 raise BadProblem( 

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

679 

680 history.append( 

681 sample.model, misfits, 

682 bootstrap_misfits, 

683 sample.pack_context()) 

684 

685 @property 

686 def niterations(self): 

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

688 

689 def get_status(self, history): 

690 if self._status_chains is None: 

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

692 

693 self._status_chains.goto(history.nmodels) 

694 

695 chains = self._status_chains 

696 problem = history.problem 

697 

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

699 row_names.append('Misfit') 

700 

701 def colum_array(data): 

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

703 arr[:data.size] = data 

704 return arr 

705 

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

707 

708 bs_mean = colum_array(chains.mean_model(ichain=None)) 

709 bs_std = colum_array(chains.standard_deviation_models( 

710 ichain=None, estimator='standard_deviation_all_chains')) 

711 

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

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

714 

715 glob_std = colum_array(chains.standard_deviation_models( 

716 ichain=0, estimator='standard_deviation_single_chain')) 

717 glob_std[-1] = num.std(chains.misfits(ichain=0)) 

718 

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

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

721 

722 glob_misfits = chains.misfits(ichain=0) 

723 

724 acceptance_latest = chains.acceptance_history[ 

725 :, -min(chains.acceptance_history.shape[1], self.ACCEPTANCE_AVG_LEN):] # noqa 

726 acceptance_avg = acceptance_latest.mean(axis=1) 

727 

728 def spark_plot(data, bins): 

729 hist, _ = num.histogram(data, bins) 

730 hist_max = num.max(hist) 

731 if hist_max == 0.0: 

732 hist_max = 1.0 

733 hist = hist / hist_max 

734 vec = num.digitize(hist, num.linspace(0., 1., len(self.SPARKS))) 

735 return ''.join([self.SPARKS[b-1] for b in vec]) 

736 

737 return OptimiserStatus( 

738 row_names=row_names, 

739 column_data=OrderedDict( 

740 zip(['BS mean', 'BS std', 

741 'Glob mean', 'Glob std', 'Glob best'], 

742 [bs_mean, bs_std, glob_mean, glob_std, glob_best])), 

743 extra_header= # noqa 

744 u'Optimiser phase: {phase}, exploring {nchains} BS chains\n' # noqa 

745 u'Global chain misfit distribution: \u2080{mf_dist}\xb9\n' 

746 u'Acceptance rate distribution: \u2080{acceptance}' 

747 u'\u2081\u2080\u2080\ufe6a (Median {acceptance_med:.1f}%)' 

748 .format( 

749 phase=phase.__class__.__name__, 

750 nchains=chains.nchains, 

751 mf_dist=spark_plot( 

752 glob_misfits, num.linspace(0., 1., 25)), 

753 acceptance=spark_plot(