Source code for grond.problems.base

'''
Base classes for Grond's problem definition and the model history container.

Common behaviour of all source models offered by Grond is implemented here.
Source model specific details are implemented in the respective submodules.
'''

import numpy as num
import math
import copy
import logging
import os.path as op
import os
import time
import struct
import threading

from pyrocko import gf, util, guts
from pyrocko.guts import Object, String, List, Dict, Int

from grond.meta import ADict, Parameter, GrondError, xjoin, Forbidden, \
    StringID, has_get_plot_classes
from ..targets import MisfitResult, MisfitTarget, TargetGroup, \
    WaveformMisfitTarget, SatelliteMisfitTarget, GNSSCampaignMisfitTarget

from grond import stats

from grond.version import __version__

guts_prefix = 'grond'
logger = logging.getLogger('grond.problems.base')
km = 1e3
as_km = dict(scale_factor=km, scale_unit='km')

g_rstate = num.random.RandomState()


def nextpow2(i):
    return 2**int(math.ceil(math.log(i) / math.log(2.)))


def nextcapacity(i):
    return int(math.ceil(i / 1024) * 1024)


def correlated_weights(values, weight_matrix):
    '''
    Applies correlated weights to values

    The resulting weighed values have to be squared! Check out
    :meth:`Problem.combine_misfits` for more information.

    :param values: Misfits or norms as :class:`numpy.Array`
    :param weight: Weight matrix, commonly the inverse of covariance matrix

    :returns: :class:`numpy.Array` weighted values
    '''
    return num.matmul(values, weight_matrix)


[docs]class ProblemConfig(Object): ''' Base class for config section defining the objective function setup. Factory for :py:class:`Problem` objects. ''' name_template = String.T() norm_exponent = Int.T(default=2) nthreads = Int.T(default=1)
[docs] def get_problem(self, event, target_groups, targets): ''' Instantiate the problem with a given event and targets. :returns: :py:class:`Problem` object ''' raise NotImplementedError
[docs]@has_get_plot_classes class Problem(Object): ''' Base class for objective function setup. Defines the *problem* to be solved by the optimiser. ''' name = String.T() ranges = Dict.T(String.T(), gf.Range.T()) dependants = List.T(Parameter.T()) norm_exponent = Int.T(default=2) base_source = gf.Source.T(optional=True) targets = List.T(MisfitTarget.T()) target_groups = List.T(TargetGroup.T()) grond_version = String.T(optional=True) nthreads = Int.T(default=1) def __init__(self, **kwargs): Object.__init__(self, **kwargs) if self.grond_version is None: self.grond_version = __version__ self._target_weights = None self._engine = None self._family_mask = None self._rstate_manager = None if hasattr(self, 'problem_waveform_parameters') and self.has_waveforms: self.problem_parameters =\ self.problem_parameters + self.problem_waveform_parameters unused_parameters = [] for p in self.problem_parameters: if p.optional and p._name not in self.ranges.keys(): unused_parameters.append(p) for p in unused_parameters: self.problem_parameters.remove(p) self.check() @classmethod def get_plot_classes(cls): from . import plot return plot.get_plot_classes() def check(self): paths = set() for grp in self.target_groups: if grp.path == 'all': continue if grp.path in paths: raise ValueError('Path %s defined more than once! In %s' % (grp.path, grp.__class__.__name__)) paths.add(grp.path) logger.debug('TargetGroup check OK.') def copy(self): o = copy.copy(self) o._target_weights = None return o def set_target_parameter_values(self, x): nprob = len(self.problem_parameters) for target in self.targets: target.set_parameter_values(x[nprob:nprob + target.nparameters]) nprob += target.nparameters def get_parameter_dict(self, model, group=None): params = [(p.name, model[ip]) for ip, p in enumerate(self.parameters) if group in p.groups or group is None] return ADict(params) def get_parameter_array(self, d): arr = num.zeros(self.nparameters, dtype=num.float) for ip, p in enumerate(self.parameters): if p.name in d.keys(): arr[ip] = d[p.name] return arr def get_parameter_index(self, param_name): return {k.name: ik for ik, k in enumerate(self.parameters)}[param_name] def get_rstate_manager(self): if self._rstate_manager is None: self._rstate_manager = RandomStateManager() return self._rstate_manager def dump_problem_info(self, dirname): fn = op.join(dirname, 'problem.yaml') util.ensuredirs(fn) guts.dump(self, filename=fn) def dump_problem_data( self, dirname, x, misfits, chains=None, sampler_context=None): fn = op.join(dirname, 'models') if not isinstance(x, num.ndarray): x = num.array(x) with open(fn, 'ab') as f: x.astype('<f8').tofile(f) fn = op.join(dirname, 'misfits') with open(fn, 'ab') as f: misfits.astype('<f8').tofile(f) if chains is not None: fn = op.join(dirname, 'chains') with open(fn, 'ab') as f: chains.astype('<f8').tofile(f) if sampler_context is not None: fn = op.join(dirname, 'choices') with open(fn, 'ab') as f: num.array(sampler_context, dtype='<i8').tofile(f) fn = op.join(dirname, 'rstate') self.get_rstate_manager().save_state(fn) def name_to_index(self, name): pnames = [p.name for p in self.combined] return pnames.index(name) @property def parameters(self): target_parameters = [] for target in self.targets: target_parameters.extend(target.target_parameters) return self.problem_parameters + target_parameters @property def parameter_names(self): return [p.name for p in self.combined] @property def dependant_names(self): return [p.name for p in self.dependants] @property def nparameters(self): return len(self.parameters) @property def ntargets(self): return len(self.targets) @property def nwaveform_targets(self): return len(self.waveform_targets) @property def nsatellite_targets(self): return len(self.satellite_targets) @property def ngnss_targets(self): return len(self.gnss_targets) @property def nmisfits(self): nmisfits = 0 for target in self.targets: nmisfits += target.nmisfits return nmisfits @property def ndependants(self): return len(self.dependants) @property def ncombined(self): return len(self.parameters) + len(self.dependants) @property def combined(self): return self.parameters + self.dependants @property def satellite_targets(self): return [t for t in self.targets if isinstance(t, SatelliteMisfitTarget)] @property def gnss_targets(self): return [t for t in self.targets if isinstance(t, GNSSCampaignMisfitTarget)] @property def waveform_targets(self): return [t for t in self.targets if isinstance(t, WaveformMisfitTarget)] @property def has_satellite(self): if self.satellite_targets: return True return False @property def has_waveforms(self): if self.waveform_targets: return True return False def set_engine(self, engine): self._engine = engine def get_engine(self): return self._engine def get_source(self, x): raise NotImplementedError def pack(self, source): raise NotImplementedError def get_gf_store_ids(self): return tuple(set([t.store_id for t in self.targets])) def get_gf_store(self, target): if not isinstance(target, str): target = target.store_id if self.get_engine() is None: raise GrondError('Cannot get GF Store, modelling is not set up.') return self.get_engine().get_store(target) def random_uniform(self, xbounds, rstate, fixed_magnitude=None): if fixed_magnitude is not None: raise GrondError( 'Setting fixed magnitude in random model generation not ' 'supported for this type of problem.') x = rstate.uniform(0., 1., self.nparameters) x *= (xbounds[:, 1] - xbounds[:, 0]) x += xbounds[:, 0] return x def preconstrain(self, x, optimizer=None): return x def extract(self, xs, i): if xs.ndim == 1: return self.extract(xs[num.newaxis, :], i)[0] if i < self.nparameters: return xs[:, i] else: return self.make_dependant( xs, self.dependants[i - self.nparameters].name) def get_target_weights(self): if self._target_weights is None: self._target_weights = num.concatenate( [target.get_combined_weight() for target in self.targets]) return self._target_weights def get_target_residuals(self): pass def inter_family_weights(self, ns): exp, root = self.get_norm_functions() family, nfamilies = self.get_family_mask() ws = num.zeros(self.nmisfits) for ifamily in range(nfamilies): mask = family == ifamily ws[mask] = 1.0 / root(num.nansum(exp(ns[mask]))) return ws
[docs] def inter_family_weights2(self, ns): ''' :param ns: 2D array with normalization factors ``ns[imodel, itarget]`` :returns: 2D array ``weights[imodel, itarget]`` ''' exp, root = self.get_norm_functions() family, nfamilies = self.get_family_mask() ws = num.zeros(ns.shape) for ifamily in range(nfamilies): mask = family == ifamily ws[:, mask] = (1.0 / root( num.nansum(exp(ns[:, mask]), axis=1)))[:, num.newaxis] return ws
def get_reference_model(self): model = num.zeros(self.nparameters) model_source_params = self.pack(self.base_source) model[:model_source_params.size] = model_source_params return model def get_parameter_bounds(self): out = [] for p in self.problem_parameters: r = self.ranges[p.name] out.append((r.start, r.stop)) for target in self.targets: for p in target.target_parameters: r = target.target_ranges[p.name_nogroups] out.append((r.start, r.stop)) return num.array(out, dtype=num.float) def get_dependant_bounds(self): return num.zeros((0, 2)) def get_combined_bounds(self): return num.vstack(( self.get_parameter_bounds(), self.get_dependant_bounds())) def raise_invalid_norm_exponent(self): raise GrondError('Invalid norm exponent: %f' % self.norm_exponent) def get_norm_functions(self): if self.norm_exponent == 2: def sqr(x): return x**2 return sqr, num.sqrt elif self.norm_exponent == 1: def noop(x): return x return noop, num.abs else: self.raise_invalid_norm_exponent()
[docs] def combine_misfits( self, misfits, extra_weights=None, extra_residuals=None, extra_correlated_weights=dict(), get_contributions=False): ''' Combine misfit contributions (residuals) to global or bootstrap misfits :param misfits: 3D array ``misfits[imodel, iresidual, 0]`` are the misfit contributions (residuals) ``misfits[imodel, iresidual, 1]`` are the normalisation contributions. It is also possible to give the misfit and normalisation contributions for a single model as ``misfits[iresidual, 0]`` and misfits[iresidual, 1]`` in which case, the first dimension (imodel) of the result will be stipped off. :param extra_weights: if given, 2D array of extra weights to be applied to the contributions, indexed as ``extra_weights[ibootstrap, iresidual]``. :param extra_residuals: if given, 2D array of perturbations to be added to the residuals, indexed as ``extra_residuals[ibootstrap, iresidual]``. :param extra_correlated_weights: if a dictionary of ``imisfit: correlated weight matrix`` is passed a correlated weight matrix is applied to the misfit and normalisation values. `imisfit` is the starting index in the misfits vector the correlated weight matrix applies to. :param get_contributions: get the weighted and perturbed contributions (don't do the sum). :returns: if no *extra_weights* or *extra_residuals* are given, a 1D array indexed as ``misfits[imodel]`` containing the global misfit for each model is returned, otherwise a 2D array ``misfits[imodel, ibootstrap]`` with the misfit for every model and weighting/residual set is returned. ''' if misfits.ndim == 2: misfits = misfits[num.newaxis, :, :] return self.combine_misfits( misfits, extra_weights, extra_residuals, extra_correlated_weights, get_contributions)[0, ...] if extra_weights is None and extra_residuals is None: return self.combine_misfits( misfits, False, False, extra_correlated_weights, get_contributions)[:, 0] assert misfits.ndim == 3 assert not num.any(extra_weights) or extra_weights.ndim == 2 assert not num.any(extra_residuals) or extra_residuals.ndim == 2 if self.norm_exponent != 2 and extra_correlated_weights: raise GrondError('Correlated weights can only be used ' ' with norm_exponent=2') exp, root = self.get_norm_functions() nmodels = misfits.shape[0] nmisfits = misfits.shape[1] # noqa mf = misfits[:, num.newaxis, :, :].copy() if num.any(extra_residuals): mf = mf + extra_residuals[num.newaxis, :, :, num.newaxis] res = mf[..., 0] norms = mf[..., 1] for imisfit, corr_weight_mat in extra_correlated_weights.items(): jmisfit = imisfit + corr_weight_mat.shape[0] for imodel in range(nmodels): corr_res = res[imodel, :, imisfit:jmisfit] corr_norms = norms[imodel, :, imisfit:jmisfit] res[imodel, :, imisfit:jmisfit] = \ correlated_weights(corr_res, corr_weight_mat) norms[imodel, :, imisfit:jmisfit] = \ correlated_weights(corr_norms, corr_weight_mat) # Apply normalization family weights (these weights depend on # on just calculated correlated norms!) weights_fam = \ self.inter_family_weights2(norms[:, 0, :])[:, num.newaxis, :] weights_fam = exp(weights_fam) res = exp(res) norms = exp(norms) res *= weights_fam norms *= weights_fam weights_tar = self.get_target_weights()[num.newaxis, num.newaxis, :] if num.any(extra_weights): weights_tar = weights_tar * extra_weights[num.newaxis, :, :] weights_tar = exp(weights_tar) res = res * weights_tar norms = norms * weights_tar if get_contributions: return res / num.nansum(norms, axis=2)[:, :, num.newaxis] result = root( num.nansum(res, axis=2) / num.nansum(norms, axis=2)) assert result[result < 0].size == 0 return result
def make_family_mask(self): family_names = set() families = num.zeros(self.nmisfits, dtype=num.int) idx = 0 for itarget, target in enumerate(self.targets): family_names.add(target.normalisation_family) families[idx:idx + target.nmisfits] = len(family_names) - 1 idx += target.nmisfits return families, len(family_names) def get_family_mask(self): if self._family_mask is None: self._family_mask = self.make_family_mask() return self._family_mask def evaluate(self, x, mask=None, result_mode='full', targets=None, nthreads=1): source = self.get_source(x) engine = self.get_engine() self.set_target_parameter_values(x) if mask is not None and targets is not None: raise ValueError('Mask cannot be defined with targets set.') targets = targets if targets is not None else self.targets for target in targets: target.set_result_mode(result_mode) modelling_targets = [] t2m_map = {} for itarget, target in enumerate(targets): t2m_map[target] = target.prepare_modelling(engine, source, targets) if mask is None or mask[itarget]: modelling_targets.extend(t2m_map[target]) u2m_map = {} for imtarget, mtarget in enumerate(modelling_targets): if mtarget not in u2m_map: u2m_map[mtarget] = [] u2m_map[mtarget].append(imtarget) modelling_targets_unique = list(u2m_map.keys()) resp = engine.process(source, modelling_targets_unique, nthreads=nthreads) modelling_results_unique = list(resp.results_list[0]) modelling_results = [None] * len(modelling_targets) for mtarget, mresult in zip( modelling_targets_unique, modelling_results_unique): for itarget in u2m_map[mtarget]: modelling_results[itarget] = mresult imt = 0 results = [] for itarget, target in enumerate(targets): nmt_this = len(t2m_map[target]) if mask is None or mask[itarget]: result = target.finalize_modelling( engine, source, t2m_map[target], modelling_results[imt:imt + nmt_this]) imt += nmt_this else: result = gf.SeismosizerError( 'target was excluded from modelling') results.append(result) return results def misfits(self, x, mask=None, nthreads=1): results = self.evaluate( x, mask=mask, result_mode='sparse', nthreads=nthreads) misfits = num.full((self.nmisfits, 2), num.nan) imisfit = 0 for target, result in zip(self.targets, results): if isinstance(result, MisfitResult): misfits[imisfit:imisfit + target.nmisfits, :] = result.misfits imisfit += target.nmisfits return misfits def forward(self, x): source = self.get_source(x) engine = self.get_engine() plain_targets = [] for target in self.targets: plain_targets.extend(target.get_plain_targets(engine, source)) resp = engine.process(source, plain_targets) results = [] for target, result in zip(plain_targets, resp.results_list[0]): if isinstance(result, gf.SeismosizerError): logger.debug( '%s.%s.%s.%s: %s' % (target.codes + (str(result),))) else: results.append(result) return results def get_random_model(self, ntries_limit=100): xbounds = self.get_parameter_bounds() for _ in range(ntries_limit): x = self.random_uniform(xbounds, rstate=g_rstate) try: return self.preconstrain(x) except Forbidden: pass raise GrondError( 'Could not find any suitable candidate sample within %i tries' % ( ntries_limit))
[docs]class ProblemInfoNotAvailable(GrondError): pass
[docs]class ProblemDataNotAvailable(GrondError): pass
[docs]class NoSuchAttribute(GrondError): pass
[docs]class InvalidAttributeName(GrondError): pass
[docs]class ModelHistory(object): ''' Write, read and follow sequences of models produced in an optimisation run. :param problem: :class:`grond.Problem` instance :param path: path to rundir, defaults to None :type path: str, optional :param mode: open mode, 'r': read, 'w': write :type mode: str, optional ''' nmodels_capacity_min = 1024 def __init__(self, problem, nchains=None, path=None, mode='r'): self.mode = mode self.problem = problem self.path = path self.nchains = nchains self._models_buffer = None self._misfits_buffer = None self._bootstraps_buffer = None self._sample_contexts_buffer = None self._sorted_misfit_idx = {} self.models = None self.misfits = None self.bootstrap_misfits = None self.sampler_contexts = None self.nmodels_capacity = self.nmodels_capacity_min self.listeners = [] self._attributes = {} if mode == 'r': self.load() @staticmethod def verify_rundir(rundir): _rundir_files = ('misfits', 'models') if not op.exists(rundir): raise ProblemDataNotAvailable( 'Directory does not exist: %s' % rundir) for f in _rundir_files: if not op.exists(op.join(rundir, f)): raise ProblemDataNotAvailable('File not found: %s' % f)
[docs] @classmethod def follow(cls, path, nchains=None, wait=20.): ''' Start following a rundir (constructor). :param path: the path to follow, a grond rundir :type path: str, optional :param wait: wait time until the folder become alive :type wait: number in seconds, optional :returns: A :py:class:`ModelHistory` instance ''' start_watch = time.time() while (time.time() - start_watch) < wait: try: cls.verify_rundir(path) problem = load_problem_info(path) return cls(problem, nchains=nchains, path=path, mode='r') except (ProblemDataNotAvailable, OSError): time.sleep(.25)
@property def nmodels(self): if self.models is None: return 0 else: return self.models.shape[0] @nmodels.setter def nmodels(self, nmodels_new): assert 0 <= nmodels_new <= self.nmodels self.models = self._models_buffer[:nmodels_new, :] self.misfits = self._misfits_buffer[:nmodels_new, :, :] if self.nchains is not None: self.bootstrap_misfits = self._bootstraps_buffer[:nmodels_new, :, :] # noqa if self._sample_contexts_buffer is not None: self.sampler_contexts = self._sample_contexts_buffer[:nmodels_new, :] # noqa @property def nmodels_capacity(self): if self._models_buffer is None: return 0 else: return self._models_buffer.shape[0] @nmodels_capacity.setter def nmodels_capacity(self, nmodels_capacity_new): if self.nmodels_capacity != nmodels_capacity_new: models_buffer = num.zeros( (nmodels_capacity_new, self.problem.nparameters), dtype=num.float) misfits_buffer = num.zeros( (nmodels_capacity_new, self.problem.nmisfits, 2), dtype=num.float) sample_contexts_buffer = num.zeros( (nmodels_capacity_new, 4), dtype=num.int) sample_contexts_buffer.fill(-1) if self.nchains is not None: bootstraps_buffer = num.zeros( (nmodels_capacity_new, self.nchains), dtype=num.float) ncopy = min(self.nmodels, nmodels_capacity_new) if self._models_buffer is not None: models_buffer[:ncopy, :] = \ self._models_buffer[:ncopy, :] misfits_buffer[:ncopy, :, :] = \ self._misfits_buffer[:ncopy, :, :] sample_contexts_buffer[:ncopy, :] = \ self._sample_contexts_buffer[:ncopy, :] self._models_buffer = models_buffer self._misfits_buffer = misfits_buffer self._sample_contexts_buffer = sample_contexts_buffer if self.nchains is not None: if self._bootstraps_buffer is not None: bootstraps_buffer[:ncopy, :] = \ self._bootstraps_buffer[:ncopy, :] self._bootstraps_buffer = bootstraps_buffer def clear(self): assert self.mode != 'r', 'History is read-only, cannot clear.' self.nmodels = 0 self.nmodels_capacity = self.nmodels_capacity_min def extend( self, models, misfits, bootstrap_misfits=None, sampler_contexts=None): nmodels = self.nmodels n = models.shape[0] nmodels_capacity_want = max( self.nmodels_capacity_min, nextpow2(nmodels + n)) if nmodels_capacity_want != self.nmodels_capacity: self.nmodels_capacity = nmodels_capacity_want self._models_buffer[nmodels:nmodels + n, :] = models self._misfits_buffer[nmodels:nmodels + n, :, :] = misfits self.models = self._models_buffer[:nmodels + n, :] self.misfits = self._misfits_buffer[:nmodels + n, :, :] if bootstrap_misfits is not None: self._bootstraps_buffer[nmodels:nmodels + n, :] = bootstrap_misfits self.bootstrap_misfits = self._bootstraps_buffer[:nmodels + n, :] if sampler_contexts is not None: self._sample_contexts_buffer[nmodels:nmodels + n, :] \ = sampler_contexts self.sampler_contexts = self._sample_contexts_buffer[ :nmodels + n, :] if self.path and self.mode == 'w': for i in range(n): self.problem.dump_problem_data( self.path, models[i, :], misfits[i, :, :], bootstrap_misfits[i, :] if bootstrap_misfits is not None else None, sampler_contexts[i, :] if sampler_contexts is not None else None) self._sorted_misfit_idx.clear() self.emit('extend', nmodels, n, models, misfits, sampler_contexts) def append( self, model, misfits, bootstrap_misfits=None, sampler_context=None): if bootstrap_misfits is not None: bootstrap_misfits = bootstrap_misfits[num.newaxis, :] if sampler_context is not None: sampler_context = sampler_context[num.newaxis, :] return self.extend( model[num.newaxis, :], misfits[num.newaxis, :, :], bootstrap_misfits, sampler_context) def load(self): self.mode = 'r' self.verify_rundir(self.path) models, misfits, bootstraps, sampler_contexts = load_problem_data( self.path, self.problem, nchains=self.nchains) self.extend(models, misfits, bootstraps, sampler_contexts)
[docs] def update(self): ''' Update history from path ''' nmodels_available = get_nmodels(self.path, self.problem) if self.nmodels == nmodels_available: return try: new_models, new_misfits, new_bootstraps, new_sampler_contexts = \ load_problem_data( self.path, self.problem, nmodels_skip=self.nmodels, nchains=self.nchains) except ValueError: return self.extend( new_models, new_misfits, new_bootstraps, new_sampler_contexts)
[docs] def add_listener(self, listener): ''' Add a listener to the history The listening class can implement the following methods: * ``extend`` ''' self.listeners.append(listener)
def emit(self, event_name, *args, **kwargs): for listener in self.listeners: slot = getattr(listener, event_name, None) if callable(slot): slot(*args, **kwargs) @property def attribute_names(self): apath = op.join(self.path, 'attributes') if not os.path.exists(apath): return [] return [fn for fn in os.listdir(apath) if StringID.regex.match(fn)] def get_attribute(self, name): if name not in self._attributes: if name not in self.attribute_names: raise NoSuchAttribute(name) path = op.join(self.path, 'attributes', name) with open(path, 'rb') as f: self._attributes[name] = num.fromfile( f, dtype='<i4', count=self.nmodels).astype(num.int) assert self._attributes[name].shape == (self.nmodels,) return self._attributes[name] def set_attribute(self, name, attribute): if not StringID.regex.match(name): raise InvalidAttributeName(name) attribute = attribute.astype(num.int) assert attribute.shape == (self.nmodels,) apath = op.join(self.path, 'attributes') if not os.path.exists(apath): os.mkdir(apath) path = op.join(apath, name) with open(path, 'wb') as f: attribute.astype('<i4').tofile(f) self._attributes[name] = attribute def ensure_bootstrap_misfits(self, optimiser): if self.bootstrap_misfits is None: problem = self.problem self.bootstrap_misfits = problem.combine_misfits( self.misfits, extra_weights=optimiser.get_bootstrap_weights(problem), extra_residuals=optimiser.get_bootstrap_residuals(problem)) def imodels_by_cluster(self, cluster_attribute): if cluster_attribute is None: return [(-1, 100.0, num.arange(self.nmodels))] by_cluster = [] try: iclusters = self.get_attribute(cluster_attribute) iclusters_avail = num.unique(iclusters) for icluster in iclusters_avail: imodels = num.where(iclusters == icluster)[0] by_cluster.append( (icluster, (100.0 * imodels.size) / self.nmodels, imodels)) if by_cluster and by_cluster[0][0] == -1: by_cluster.append(by_cluster.pop(0)) except NoSuchAttribute: logger.warn( 'Attribute %s not set in run %s.\n' ' Skipping model retrieval by clusters.' % ( cluster_attribute, self.problem.name)) return by_cluster def models_by_cluster(self, cluster_attribute): if cluster_attribute is None: return [(-1, 100.0, self.models)] return [ (icluster, percentage, self.models[imodels]) for (icluster, percentage, imodels) in self.imodels_by_cluster(cluster_attribute)] def mean_sources_by_cluster(self, cluster_attribute): return [ (icluster, percentage, stats.get_mean_source(self.problem, models)) for (icluster, percentage, models) in self.models_by_cluster(cluster_attribute)] def get_sorted_misfits_idx(self, chain=0): if chain not in self._sorted_misfit_idx.keys(): self._sorted_misfit_idx[chain] = num.argsort( self.bootstrap_misfits[:, chain]) return self._sorted_misfit_idx[chain] def get_sorted_misfits(self, chain=0): isort = self.get_sorted_misfits_idx(chain) return self.bootstrap_misfits[:, chain][isort] def get_sorted_models(self, chain=0): isort = self.get_sorted_misfits_idx(chain=0) return self.models[isort, :] def get_sorted_primary_misfits(self): return self.get_sorted_misfits(chain=0) def get_sorted_primary_models(self): return self.get_sorted_models(chain=0) def get_best_model(self, chain=0): return self.get_sorted_models(chain)[0, ...] def get_best_misfit(self, chain=0): return self.get_sorted_misfits(chain)[0] def get_mean_model(self): return num.mean(self.models, axis=0) def get_mean_misfit(self, chain=0): return num.mean(self.bootstrap_misfits[:, chain]) def get_best_source(self, chain=0): return self.problem.get_source(self.get_best_model(chain)) def get_mean_source(self, chain=0): return self.problem.get_source(self.get_mean_model()) def get_chain_misfits(self, chain=0): return self.bootstrap_misfits[:, chain] def get_primary_chain_misfits(self): return self.get_chain_misfits(chain=0)
class RandomStateManager(object): MAX_LEN = 64 save_struct = struct.Struct('%ds7s2496sqqd' % MAX_LEN) def __init__(self): self.rstates = {} self.lock = threading.Lock() def get_rstate(self, name, seed=None): assert len(name) <= self.MAX_LEN if name not in self.rstates: self.rstates[name] = num.random.RandomState(seed) return self.rstates[name] @property def nstates(self): return len(self.rstates) def save_state(self, fname): with self.lock: with open(fname, 'wb') as f: for name, rstate in self.rstates.items(): s, arr, pos, has_gauss, chached_gauss = rstate.get_state() f.write( self.save_struct.pack( name.encode(), s.encode(), arr.tobytes(), pos, has_gauss, chached_gauss)) def load_state(self, fname): with self.lock: with open(fname, 'rb') as f: while True: buff = f.read(self.save_struct.size) if not buff: break name, s, arr, pos, has_gauss, chached_gauss = \ self.save_struct.unpack(buff) name = name.replace(b'\x00', b'').decode() s = s.replace(b'\x00', b'').decode() arr = num.frombuffer(arr, dtype=num.uint32) rstate = num.random.RandomState() rstate.set_state((s, arr, pos, has_gauss, chached_gauss)) self.rstates[name] = rstate def get_nmodels(dirname, problem): fn = op.join(dirname, 'models') with open(fn, 'r') as f: nmodels1 = os.fstat(f.fileno()).st_size // (problem.nparameters * 8) fn = op.join(dirname, 'misfits') with open(fn, 'r') as f: nmodels2 = os.fstat(f.fileno()).st_size // (problem.nmisfits * 2 * 8) return min(nmodels1, nmodels2) def load_problem_info_and_data(dirname, subset=None, nchains=None): problem = load_problem_info(dirname) models, misfits, bootstraps, sampler_contexts = load_problem_data( xjoin(dirname, subset), problem, nchains=nchains) return problem, models, misfits, bootstraps, sampler_contexts def load_optimiser_info(dirname): fn = op.join(dirname, 'optimiser.yaml') return guts.load(filename=fn) def load_problem_info(dirname): try: fn = op.join(dirname, 'problem.yaml') return guts.load(filename=fn) except OSError as e: logger.debug(e) raise ProblemInfoNotAvailable( 'No problem info available (%s).' % dirname) def load_problem_data(dirname, problem, nmodels_skip=0, nchains=None): def get_chains_fn(): for fn in (op.join(dirname, 'bootstraps'), op.join(dirname, 'chains')): if op.exists(fn): return fn return False try: nmodels = get_nmodels(dirname, problem) - nmodels_skip fn = op.join(dirname, 'models') offset = nmodels_skip * problem.nparameters * 8 models = num.memmap( fn, dtype='<f8', offset=offset, shape=(nmodels, problem.nparameters))\ .astype(num.float, copy=False) fn = op.join(dirname, 'misfits') offset = nmodels_skip * problem.nmisfits * 2 * 8 misfits = num.memmap( fn, dtype='<f8', offset=offset, shape=(nmodels, problem.nmisfits, 2))\ .astype(num.float, copy=False) chains = None fn = get_chains_fn() if fn and nchains is not None: offset = nmodels_skip * nchains * 8 chains = num.memmap( fn, dtype='<f8', offset=offset, shape=(nmodels, nchains))\ .astype(num.float, copy=False) sampler_contexts = None fn = op.join(dirname, 'choices') if op.exists(fn): offset = nmodels_skip * 4 * 8 sampler_contexts = num.memmap( fn, dtype='<i8', offset=offset, shape=(nmodels, 4))\ .astype(num.int, copy=False) fn = op.join(dirname, 'rstate') if op.exists(fn): problem.get_rstate_manager().load_state(fn) except OSError as e: logger.debug(str(e)) raise ProblemDataNotAvailable( 'No problem data available (%s).' % dirname) return models, misfits, chains, sampler_contexts __all__ = ''' ProblemConfig Problem ModelHistory ProblemInfoNotAvailable ProblemDataNotAvailable load_problem_info load_problem_info_and_data InvalidAttributeName NoSuchAttribute '''.split()