Compare commits

...

14 Commits

Author SHA1 Message Date
braunfuss 371f7a85c3 call problem rstate manager in target balancing analyser 1 year ago
braunfuss e7adb168ea Merge branch 'dev-rectangularv' of https://git.pyrocko.org/pyrocko/grond into dev-continue-velocity 1 year ago
miili 1aa6210ff2 docs: added continue CLI and mentions in text 1 year ago
miili 54bb4532c7 parameter: `optional=True` now is the default 1 year ago
miili 43e4ac6715 docs: added velocity and optional parameters 1 year ago
miili db29303a3a rectangular: adding velocity parameter 1 year ago
miili 9b9600744e problem: adding Parameter optionality 1 year ago
miili f7d7eb0c6a RectangularSource: adding velocit param 1 year ago
miili 36c2e78bf9 dataset: fix kite import 1 year ago
miili 4857406ff5 docs: upd 1 year ago
miili 100fbe07c2 removing debug statements 1 year ago
miili 016ef5e78b problem: Added `RandomStateManager` 1 year ago
miili 5e196b4ca8 grond continue continued 2 years ago
miili b5c0bf1a6d adding `grond continue` 2 years ago
  1. 13
      CHANGELOG.md
  2. 2
      docs/Makefile
  3. 1
      docs/source/cli/index.rst
  4. 6
      docs/source/config/problems/index.rst
  5. 23
      docs/source/overview/index.rst
  6. 20
      docs/source/quickstart/index.rst
  7. 4
      src/analysers/base.py
  8. 2
      src/analysers/target_balancing/analyser.py
  9. 46
      src/apps/grond.py
  10. 87
      src/core.py
  11. 1
      src/data/snippets/problem_rectangular.gronf
  12. 3
      src/dataset.py
  13. 8
      src/environment.py
  14. 3
      src/meta.py
  15. 2
      src/monitor.py
  16. 41
      src/optimisers/highscore/optimiser.py
  17. 117
      src/problems/base.py
  18. 23
      src/problems/rectangular/problem.py
  19. 57
      test/test_rstate_manager.py

13
CHANGELOG.md

@ -9,6 +9,13 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
### Added
- `grond go` utilises `--threads` arguments for various tasks.
- InSAR plots improved, now plotting 'best' and 'mean' // Improvements
- `grond continue` subcommand to continue aborted and re-configured runs
- Added optional keyword for `Parameters`, defaults to `optional=True`
### Changed
- `RectangularSource` now has a `velocity` parameter
- `RandomStateManager` is now centrally attached to `Problem`.
### Fixed
- Improvements in documentation and examples.
@ -81,7 +88,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
- Transparent event loading and checking.
- Noise analyser: target groups are now handled independently. Each group now
uses its own threshold in weeding mode.
- Improved error handling (`grond check`, instrument responses,
- Improved error handling (`grond check`, instrument responses,
- Only exclude waveform targets when `distance_min` constraint is given in
`problem_config`.
- Improved method chapter in documentation.
@ -110,7 +117,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
- New problem config sections in `grond init`.
### Changed
- Documentation of problems configurations are now centralised at
- Documentation of problems configurations are now centralised at
`src/data/snippets`.
- Output of `grond init list`.
@ -146,7 +153,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
multiple misfits.
- Optimiser can now be configured to yield exactly reproducible results by
providing seed values for all random number generators involved.
- Plots `sequence` and `fits_waveform`: layout as single plot figures by
- Plots `sequence` and `fits_waveform`: layout as single plot figures by
default.
### Fixed

2
docs/Makefile

@ -8,7 +8,7 @@ SPHINXPROJ = grond
SOURCEDIR = source
BUILDDIR = build
GRONDCOMMANDS = scenario init events check go forward harvest plot movie \
GRONDCOMMANDS = scenario init events check go continue forward harvest plot movie \
export report diff qc-polarization upgrade-config version
# Put it first so that "make" without argument is like "make help".

1
docs/source/cli/index.rst

@ -24,6 +24,7 @@ subcommands.
events
check
go
continue
forward
harvest
plot

6
docs/source/config/problems/index.rst

@ -53,7 +53,7 @@ problem configurations:
defines the norm of combining several `normalization_family` in the global misfit. This integer value is 1 or larger. Please find here more information on the global `misfit calculation in Grond`_.
``ranges``
defines the bounds of individual and specific source model parameters. See the details for the source ranges of different problems in the sections below.
defines the bounds of individual and specific source model parameters. See the details for the source ranges of different problems in the sections below. Ranges are optional, if a parameter range is not given, the parameter will be fixed to a default value.
An example for the configuration of a rectangular fault problem is given here:
@ -254,6 +254,10 @@ For the source parameter configuration, please note that the last three paramete
``nucleation_y``
relative along-dip position of the rupture nucleation point on the fault to the centre location. This parameter may range from -1 to 1. With 0 being in the centre, -1 being at the top fault edge, 1 at the bottom fault edge, and 0.5 is half-way between centroid and bottom fault edge.
``velocity``
speed of rupture front m/s. this parameter is optional, default is 3500 m/s.
.. _volume_point:

23
docs/source/overview/index.rst

@ -232,13 +232,13 @@ Before running the optimisation, you may want to check your dataset and configur
::
grond check <configfile> <eventname>
grond check <configfile> [eventname]
Now, you may start the optimization for a given event using
::
grond go <configfile> <eventname>
grond go <configfile> [eventname]
During the optimisation, results are aggregated in an output directory, referred to as `<rundir>` in the configuration and documentation.
@ -254,6 +254,7 @@ During the optimisation, results are aggregated in an output directory, referred
│ ├── laquila2009_joint.grun
│ │ ├── ... # some bookkeeping yaml-files
│ │ ├── optimiser.yaml
│ │ ├── rstate
│ │ ├── models
│ │ ├── misfits
│ │ └── harvest
@ -267,6 +268,14 @@ During the optimisation, results are aggregated in an output directory, referred
You find detailed information on the misfit configuration and model space sampling in the section :doc:`/config/optimisers/index`.
.. note ::
An interrupted optimisation can be prolonged and continued with:
::
grond continue <configfile> [eventname]
Results and visualisation
-------------------------
@ -320,6 +329,16 @@ The results can be exported in various ways by running the subcommand
See the command reference of :option:`grond export` for more details.
.. tip::
An existing grond optimisation can be extended, e.g. by adding iterations to reach a better convergence. Add iterations
or a ``SamplerPhase`` in the corresponding configuration file (for details, please check the sections :doc:`/config/optimisers/index` and :doc:`/config/index`) and use:
::
grond continue <configfile> [eventname]
Terminology
-----------

20
docs/source/quickstart/index.rst

@ -37,20 +37,34 @@ With the following few commands, we will let Grond
*Note: precomputed Green's functions (GF) needed to run this example will be
downloaded from the internet.*
**Check the data setup and configuration**
**Checking data setup and configuration**
.. code-block :: sh
cd my_first_project
grond check config/scenario.gronf
**Start the optimisation**
**Starting optimisation**
.. code-block :: sh
grond go config/scenario.gronf
**Plot the results in a report**
.. note ::
Interrupted or prolonged optimisation can be continued with:
.. code-block :: sh
grond go config/scenario.gronf
**Continuing an aborted/modified optimisation**
.. code-block :: sh
grond continue config/scenario.gronf
**Plotting the results in a report**
.. code-block :: sh

4
src/analysers/base.py

@ -10,6 +10,10 @@ class Analyser(object):
def analyse(self, problem, ds):
pass
def get_rstate(self, problem):
return problem.get_rstate_manager().get_rstate(
self.__class__.__name__)
class AnalyserConfig(Object):

2
src/analysers/target_balancing/analyser.py

@ -69,7 +69,7 @@ class TargetBalancingAnalyser(Analyser):
xbounds = wproblem.get_parameter_bounds()
misfits = num.zeros((self.niter, wproblem.ntargets, 2))
rstate = num.random.RandomState(123)
rstate = problem.get_rstate(problem, seed=123)
isbad_mask = None

46
src/apps/grond.py

@ -46,6 +46,8 @@ subcommand_descriptions = {
'events': 'print available event names for given configuration',
'check': 'check data and configuration',
'go': 'run Grond optimisation',
'continue': 'continue a Grond optimisation in case of an'
' interruption or more iterations are needed',
'forward': 'run forward modelling',
'harvest': 'manually run harvesting',
'cluster': 'run cluster analysis on result ensemble',
@ -69,6 +71,7 @@ subcommand_usages = {
'events': 'events <configfile>',
'check': 'check <configfile> <eventnames> ... [options]',
'go': 'go <configfile> <eventnames> ... [options]',
'continue': 'continue <configfile> <eventnames> ... [options]',
'forward': (
'forward <rundir> [options]',
'forward <configfile> <eventnames> ... [options]'),
@ -120,6 +123,7 @@ Subcommands:
events %(events)s
check %(check)s
go %(go)s
continue %(continue)s
forward %(forward)s
harvest %(harvest)s
cluster %(cluster)s
@ -763,6 +767,48 @@ def command_go(args):
die(str(e))
def command_continue(args):
from grond.environment import Environment
def setup(parser):
parser.add_option(
'--no-preserve', dest='no_preserve', action='store_true',
help='do not preserve old rundir')
parser.add_option(
'--status', dest='status', default='state',
type='choice', choices=['state', 'quiet'],
help='status output selection (choices: state, quiet, default: '
'state)')
parser.add_option(
'--parallel', dest='nparallel', type=int, default=1,
help='set number of events to process in parallel, '
'if set to more than one, --status=quiet is implied.')
parser.add_option(
'--threads', dest='nthreads', type=int, default=1,
help='set number of threads per process (default: 1). '
'Set to 0 to use all available cores.')
parser, options, args = cl_parse('continue', args, setup)
try:
env = Environment(args)
status = options.status
if options.nparallel != 1:
status = 'quiet'
grond.continue_run(
env,
preserve=~bool(options.no_preserve),
status=status,
nparallel=options.nparallel,
nthreads=options.nthreads)
except grond.GrondError as e:
die(str(e))
def command_forward(args):
from grond.environment import Environment

87
src/core.py

@ -10,7 +10,7 @@ import math
import os.path as op
import numpy as num
from pyrocko.guts import Object, String, Float, List
from pyrocko.guts import Object, String, Float, List, clone
from pyrocko import gf, trace, guts, util, weeding
from pyrocko import parimap, model, marker as pmarker
@ -430,7 +430,23 @@ def go(environment,
nparallel=1, status='state', nthreads=0):
g_data = (environment, force, preserve,
status, nparallel, nthreads)
status, nparallel, nthreads, False)
g_state[id(g_data)] = g_data
nevents = environment.nevents_selected
for x in parimap.parimap(
process_event,
range(environment.nevents_selected),
[id(g_data)] * nevents,
nprocs=nparallel):
pass
def continue_run(environment, force=False, preserve=False,
nparallel=1, status='state', nthreads=0):
g_data = (environment, force, preserve,
status, nparallel, nthreads, True)
g_state[id(g_data)] = g_data
nevents = environment.nevents_selected
@ -445,12 +461,12 @@ def go(environment,
def process_event(ievent, g_data_id):
environment, force, preserve, status, nparallel, nthreads = \
env, force, preserve, status, nparallel, nthreads, continue_run = \
g_state[g_data_id]
config = environment.get_config()
event_name = environment.get_selected_event_names()[ievent]
nevents = environment.nevents_selected
config = env.get_config()
event_name = env.get_selected_event_names()[ievent]
nevents = env.nevents_selected
tstart = time.time()
ds = config.get_dataset(event_name)
@ -466,31 +482,46 @@ def process_event(ievent, g_data_id):
rundir = expand_template(
config.rundir_template,
dict(problem_name=problem.name))
environment.set_rundir_path(rundir)
env.set_rundir_path(rundir)
nold_rundirs = len(glob.glob(rundir + '*'))
if op.exists(rundir):
if op.exists(rundir) and not continue_run:
if preserve:
nold_rundirs = len(glob.glob(rundir + '*'))
shutil.move(rundir, rundir+'-old-%d' % (nold_rundirs))
shutil.move(rundir, rundir+'-old-%d' % nold_rundirs)
elif force:
shutil.rmtree(rundir)
else:
logger.warn('Skipping problem "%s": rundir already exists: %s' %
(problem.name, rundir))
logger.warn('Skipping problem "%s": rundir already exists: %s',
problem.name, rundir)
return
util.ensuredir(rundir)
if op.exists(rundir) and continue_run:
logger.info(
'Continuing event %i / %i', ievent + 1, nevents)
logger.info(
'Starting event %i / %i' % (ievent+1, nevents))
env_old = Environment(rundir)
history = env_old.get_history()
targets = env_old.get_problem().targets
for target in targets:
target.set_dataset(ds)
logger.info('Rundir: %s' % rundir)
problem.targets = targets
logger.info('Analysing problem "%s".' % problem.name)
if preserve:
shutil.copytree(rundir, rundir+'-old-%d' % nold_rundirs)
for analyser_conf in config.analyser_configs:
analyser = analyser_conf.get_analyser()
analyser.analyse(problem, ds)
elif not op.exists(rundir) and continue_run:
logger.warn('Cannot find rundir %s to continue...', rundir)
return
else:
logger.info(
'Starting event %i / %i', ievent+1, nevents)
history = None
util.ensuredir(rundir)
logger.info('Rundir: %s', rundir)
basepath = config.get_basepath()
config.change_basepath(rundir)
@ -500,7 +531,13 @@ def process_event(ievent, g_data_id):
optimiser = config.optimiser_config.get_optimiser()
optimiser.set_nthreads(nthreads)
optimiser.init_bootstraps(problem)
if not continue_run:
logger.info('Analysing problem "%s".', problem.name)
for analyser_conf in config.analyser_configs:
analyser = analyser_conf.get_analyser()
analyser.analyse(problem, ds)
optimiser.init_bootstraps(problem)
problem.dump_problem_info(rundir)
monitor = None
@ -524,7 +561,8 @@ def process_event(ievent, g_data_id):
optimiser.optimise(
problem,
rundir=rundir)
rundir=rundir,
history=history)
harvest(rundir, problem, force=True)
@ -540,10 +578,10 @@ def process_event(ievent, g_data_id):
tstop = time.time()
logger.info(
'Stop %i / %i (%g min)' % (ievent+1, nevents, (tstop - tstart)/60.))
'Stop %i / %i (%g min)', ievent+1, nevents, (tstop - tstart)/60.)
logger.info(
'Done with problem "%s", rundir is "%s".' % (problem.name, rundir))
'Done with problem "%s", rundir is "%s".', problem.name, rundir)
class ParameterStats(Object):
@ -779,6 +817,7 @@ __all__ = '''
harvest
cluster
go
continue_run
get_event_names
check
export

1
src/data/snippets/problem_rectangular.gronf

@ -49,3 +49,4 @@ problem_config: !grond.RectangularProblemConfig
time: '-15. .. 10. | add'
nucleation_x: '-1. .. 1.'
nucleation_y: '-1. .. 1.'
velocity: '1700. .. 3800.'

3
src/dataset.py

@ -316,9 +316,8 @@ class Dataset(object):
' please install from https://pyrocko.org.')
logger.debug('Loading kite scene from "%s"...' % filename)
scene = Scene()
scene = Scene.load(filename)
scene._log.setLevel(logger.level)
scene.load(filename)
try:
self.get_kite_scene(scene.meta.scene_id)

8
src/environment.py

@ -6,7 +6,7 @@ import os
from grond.config import read_config, write_config
from grond import meta, run_info
from grond.problems.base import load_optimiser_info, load_problem_info, \
ModelHistory
ModelHistory, RandomStateManager
op = os.path
@ -264,14 +264,14 @@ class Environment(object):
def get_history(self, subset=None):
if subset not in self._histories:
self._histories[subset] = \
history = \
ModelHistory(
self.get_problem(),
nchains=self.get_optimiser().nchains,
path=meta.xjoin(self.get_rundir_path(), subset))
self._histories[subset].ensure_bootstrap_misfits(
self.get_optimiser())
history.ensure_bootstrap_misfits(self.get_optimiser())
self._histories[subset] = history
return self._histories[subset]

3
src/meta.py

@ -4,7 +4,7 @@ import re
import numpy as num
import os.path as op
from string import Template
from pyrocko.guts import Object, String, Float, Unicode, StringPattern
from pyrocko.guts import Object, String, Float, Unicode, StringPattern, Bool
from pyrocko import util
guts_prefix = 'grond'
@ -161,6 +161,7 @@ class Parameter(Object):
scale_factor = Float.T(default=1., optional=True)
scale_unit = Unicode.T(optional=True)
label = Unicode.T(optional=True)
optional = Bool.T(default=True, optional=True)
def __init__(self, *args, **kwargs):
if len(args) >= 1:

2
src/monitor.py

@ -86,7 +86,7 @@ class GrondMonitor(threading.Thread):
self._tm = None
def run(self):
logger.info('Waiting to follow environment %s...' % self.rundir)
logger.info('Waiting to follow environment %s...', self.rundir)
env = Environment.discover(self.rundir)
if env is None:
logger.error('Could not attach to Grond environment.')

41
src/optimisers/highscore/optimiser.py

@ -114,9 +114,10 @@ class SamplerPhase(Object):
Object.__init__(self, *args, **kwargs)
self._rstate = None
def get_rstate(self):
def get_rstate(self, problem):
if self._rstate is None:
self._rstate = num.random.RandomState(self.seed)
self._rstate = problem.get_rstate_manager().get_rstate(
self.__class__.__name__, self.seed)
return self._rstate
@ -154,7 +155,8 @@ class UniformSamplerPhase(SamplerPhase):
def get_raw_sample(self, problem, iiter, chains):
xbounds = problem.get_parameter_bounds()
return Sample(model=problem.random_uniform(xbounds, self.get_rstate()))
return Sample(
model=problem.random_uniform(xbounds, self.get_rstate(problem)))
class DirectedSamplerPhase(SamplerPhase):
@ -200,7 +202,7 @@ class DirectedSamplerPhase(SamplerPhase):
return s or 1.0
def get_raw_sample(self, problem, iiter, chains):
rstate = self.get_rstate()
rstate = self.get_rstate(problem)
factor = self.get_scatter_scale_factor(iiter)
npar = problem.nparameters
pnames = problem.parameter_names
@ -476,10 +478,10 @@ class HighScoreOptimiser(Optimiser):
self._status_chains = None
self._rstate_bootstrap = None
def get_rstate_bootstrap(self):
def get_rstate_bootstrap(self, problem):
if self._rstate_bootstrap is None:
self._rstate_bootstrap = num.random.RandomState(
self.bootstrap_seed)
self._rstate_bootstrap = problem.get_rstate_manager().get_rstate(
'bootstraps', self.bootstrap_seed)
return self._rstate_bootstrap
@ -496,7 +498,7 @@ class HighScoreOptimiser(Optimiser):
ws = make_bayesian_weights(
self.nbootstrap,
nmisfits=nmisfits_w,
rstate=self.get_rstate_bootstrap())
rstate=self.get_rstate_bootstrap(problem))
imf = 0
for t in problem.targets:
@ -513,7 +515,7 @@ class HighScoreOptimiser(Optimiser):
for t in problem.targets:
if t.can_bootstrap_residuals:
t.init_bootstrap_residuals(
self.nbootstrap, rstate=self.get_rstate_bootstrap(),
self.nbootstrap, rstate=self.get_rstate_bootstrap(problem),
nthreads=self._nthreads)
else:
t.set_bootstrap_residuals(
@ -606,19 +608,30 @@ class HighScoreOptimiser(Optimiser):
self._tlog_last = t
def optimise(self, problem, rundir=None):
def optimise(self, problem, rundir=None, history=None):
if rundir is not None:
self.dump(filename=op.join(rundir, 'optimiser.yaml'))
history = ModelHistory(problem,
nchains=self.nchains,
path=rundir, mode='w')
if not history:
history = ModelHistory(
problem,
nchains=self.nchains,
path=rundir, mode='w')
chains = self.chains(problem, history)
if history.mode == 'r':
if history.nmodels == self.niterations:
return
logger.info('Continuing run at %d iterations...', history.nmodels)
history.mode = 'w'
chains.goto()
niter = self.niterations
isbad_mask = None
self._tlog_last = 0
for iiter in range(niter):
for iiter in range(history.nmodels, niter):
iphase, phase, iiter_phase = self.get_sampler_phase(iiter)
self.log_progress(problem, iiter, niter, phase, iiter_phase)

117
src/problems/base.py

@ -12,6 +12,7 @@ import logging
import os.path as op
import os
import time
import struct
from pyrocko import gf, util, guts
from pyrocko.guts import Object, String, List, Dict, Int
@ -34,7 +35,7 @@ g_rstate = num.random.RandomState()
def nextpow2(i):
return 2**int(math.ceil(math.log(i)/math.log(2.)))
return 2**int(math.ceil(math.log(i) / math.log(2.)))
def correlated_weights(values, weight_matrix):
@ -97,11 +98,20 @@ class Problem(Object):
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
@ -128,7 +138,7 @@ class Problem(Object):
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])
target.set_parameter_values(x[nprob:nprob + target.nparameters])
nprob += target.nparameters
def get_parameter_dict(self, model, group=None):
@ -145,6 +155,11 @@ class Problem(Object):
arr[ip] = d[p.name]
return arr
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)
@ -174,6 +189,9 @@ class Problem(Object):
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)
@ -292,12 +310,12 @@ class Problem(Object):
return xs[:, i]
else:
return self.make_dependant(
xs, self.dependants[i-self.nparameters].name)
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])
[target.get_combined_weight() for target in self.targets])
return self._target_weights
@ -385,7 +403,6 @@ class Problem(Object):
extra_residuals=None,
extra_correlated_weights=dict(),
get_contributions=False):
'''
Combine misfit contributions (residuals) to global or bootstrap misfits
@ -566,7 +583,7 @@ class Problem(Object):
result = target.finalize_modelling(
engine, source,
t2m_map[target],
modelling_results[imt:imt+nmt_this])
modelling_results[imt:imt + nmt_this])
imt += nmt_this
else:
@ -585,7 +602,7 @@ class Problem(Object):
imisfit = 0
for target, result in zip(self.targets, results):
if isinstance(result, MisfitResult):
misfits[imisfit:imisfit+target.nmisfits, :] = result.misfits
misfits[imisfit:imisfit + target.nmisfits, :] = result.misfits
imisfit += target.nmisfits
@ -798,20 +815,21 @@ class ModelHistory(object):
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_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, :, :]
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, :]
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, :] \
self._sample_contexts_buffer[nmodels:nmodels + n, :] \
= sampler_contexts
self.sampler_contexts = self._sample_contexts_buffer[:nmodels+n, :]
self.sampler_contexts = self._sample_contexts_buffer[
:nmodels + n, :]
if self.path and self.mode == 'w':
for i in range(n):
@ -823,7 +841,6 @@ class ModelHistory(object):
if sampler_contexts is not None else None)
self._sorted_misfit_idx.clear()
self.emit('extend', nmodels, n, models, misfits, sampler_contexts)
def append(
@ -1025,6 +1042,54 @@ class ModelHistory(object):
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 = {}
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 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 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:
@ -1075,8 +1140,8 @@ def load_problem_data(dirname, problem, nmodels_skip=0, nchains=None):
with open(fn, 'r') as f:
f.seek(nmodels_skip * problem.nparameters * 8)
models = num.fromfile(
f, dtype='<f8',
count=nmodels * problem.nparameters)\
f, dtype='<f8',
count=nmodels * problem.nparameters)\
.astype(num.float)
models = models.reshape((nmodels, problem.nparameters))
@ -1085,8 +1150,8 @@ def load_problem_data(dirname, problem, nmodels_skip=0, nchains=None):
with open(fn, 'r') as f:
f.seek(nmodels_skip * problem.nmisfits * 2 * 8)
misfits = num.fromfile(
f, dtype='<f8',
count=nmodels*problem.nmisfits*2)\
f, dtype='<f8',
count=nmodels * problem.nmisfits * 2)\
.astype(num.float)
misfits = misfits.reshape((nmodels, problem.nmisfits, 2))
@ -1096,8 +1161,8 @@ def load_problem_data(dirname, problem, nmodels_skip=0, nchains=None):
with open(fn, 'r') as f:
f.seek(nmodels_skip * nchains * 8)
chains = num.fromfile(
f, dtype='<f8',
count=nmodels*nchains)\
f, dtype='<f8',
count=nmodels * nchains)\
.astype(num.float)
chains = chains.reshape((nmodels, nchains))
@ -1108,11 +1173,15 @@ def load_problem_data(dirname, problem, nmodels_skip=0, nchains=None):
with open(fn, 'r') as f:
f.seek(nmodels_skip * 4 * 8)
sampler_contexts = num.fromfile(
f, dtype='<i8',
count=nmodels*4).astype(num.int)
f, dtype='<i8',
count=nmodels * 4).astype(num.int)
sampler_contexts = sampler_contexts.reshape((nmodels, 4))
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(

23
src/problems/rectangular/problem.py

@ -4,8 +4,7 @@ import logging
from pyrocko import gf, util
from pyrocko.guts import String, Float, Dict, Int
from grond.meta import expand_template, Parameter, has_get_plot_classes, \
GrondError
from grond.meta import expand_template, Parameter, has_get_plot_classes
from ..base import Problem, ProblemConfig
@ -52,27 +51,24 @@ class RectangularProblemConfig(ProblemConfig):
@has_get_plot_classes
class RectangularProblem(Problem):
# nucleation_x
# nucleation_y
# time
# stf
problem_parameters = [
Parameter('east_shift', 'm', label='Easting', **as_km),
Parameter('north_shift', 'm', label='Northing', **as_km),
Parameter('depth', 'm', label='Depth', **as_km),
Parameter('length', 'm', label='Length', **as_km),
Parameter('width', 'm', label='Width', **as_km),
Parameter('slip', 'm', label='Slip'),
Parameter('length', 'm', label='Length', optional=False, **as_km),
Parameter('width', 'm', label='Width', optional=False, **as_km),
Parameter('slip', 'm', label='Slip', optional=False),
Parameter('strike', 'deg', label='Strike'),
Parameter('dip', 'deg', label='Dip'),
Parameter('rake', 'deg', label='Rake')
]
]
problem_waveform_parameters = [
Parameter('nucleation_x', 'offset', label='Nucleation X'),
Parameter('nucleation_y', 'offset', label='Nucleation Y'),
Parameter('time', 's', label='Time'),
Parameter('velocity', 'm/s', label='Rupture Velocity')
]
dependants = []
@ -99,11 +95,6 @@ class RectangularProblem(Problem):
return source
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 = num.zeros(self.nparameters)
for i in range(self.nparameters):
x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
@ -114,7 +105,7 @@ class RectangularProblem(Problem):
# source = self.get_source(x)
# if any(self.distance_min > source.distance_to(t)
# for t in self.targets):
# raise Forbidden()
# raise Forbidden()
return x
@classmethod

57
test/test_rstate_manager.py

@ -0,0 +1,57 @@
import time
import tempfile
from grond.problems.base import RandomStateManager
def test_get_state():
rstate = RandomStateManager()
rs = rstate.get_rstate('test')
rs1 = rstate.get_rstate('test', seed=123123123)
assert rs is rs1
assert rstate.nstates == 1
rs2 = rstate.get_rstate('test2')
assert rstate.nstates == 2
def test_save_load_state():
rstate = RandomStateManager()
rstate.get_rstate('test')
rs = rstate.get_rstate('test2', seed=4414412)
rs.uniform()
with tempfile.NamedTemporaryFile() as f:
rstate.save_state(f.name)
unif = rs.uniform()
rstate_load = RandomStateManager()
rstate_load.load_state(f.name)
rs_loaded = rstate_load.get_rstate('test2')
assert unif == rs_loaded.uniform()
assert rstate_load.nstates == 2
def test_performance():
rstate_perf = RandomStateManager()
for r in range(1000):
rstate_perf.get_rstate('test-%d' % r)
with tempfile.NamedTemporaryFile() as f:
t0 = time.time()
rstate_perf.save_state(f.name)
print('saved %d states in %f s'
% (rstate_perf.nstates, time.time()-t0))
rstate_load = RandomStateManager()
t0 = time.time()
rstate_load.load_state(f.name)
print('loaded %d states in %f s'
% (rstate_load.nstates, time.time()-t0))
if __name__ == '__main__':
test_get_state()
test_save_load_state()
test_performance()
Loading…
Cancel
Save