example.cfg

[model]

# Set custom output label (e.g. label_YYYYMMDDTHHMMSSS).
# Default is "output"
output_label = test

# Output path (default current working directory)
# Converts to $output_path/$output_label_YYYYMMDDTHHMMSSS
#output_dir = 

# Scratch directory (where traveltimes, raypaths, etc stored. Recommend RAM partition (e.g. /dev/shm)
# Default is output_dir/scratch
scratch_dir = /dev/shm/pyvoro_test

# Log file name (default pyvorotomo.log)
#log_file =

# Define event catalog loadable via pd.read_hdf(path, key="events") (and key="arrivals")
events_path = /path/to/catalog.h5

# Define network/stations loadable via pd.read_hdf(path, key="stations")
stations_path = /path/to/stations.h5

# Velocity model loadable via pykonal.fields.read_hdf(path)
initial_pwave_path = /path/to/P.vel.h5
initial_swave_path = /path/to/S.vel.h5

# Generate 1D velocity models (as csv) during the final iteration
# These are stacked and QC'd the same as 3D calculation with uniform
# dz depths (same as input 3D model)
# Default is "True" (trival extra processing time)
output_1d_model = True

# Perform a "checkerboard" resolution test at the end?
perform_res_test = True
# Checkerboard horizontal size (km), and strength of anomaly (decimal ratio)
res_test_size_mag = 75,0.08
# Depths (km) where checkerboard will flip polarity (default is 10,25,50,70,120,170,230)
res_test_layers = 6,18,50

# If pyvorotomo is ran with the -t flag, then ONLY do the checkerboard test assuming
#   the latest model and catalogs found in the rerun_restest directory
rerun_restest = /path/to/pyvorotomo_data/output_SWAN06

# IF SET, filter stations (only! events will be subsequently filtered via max_dist)
# format is: minlat,minlon,maxlat,maxlon (degrees)
#map_filter = -31,115,-27,120



[algorithm]
# Number of iterations to perform.
# Recommend 5+, monitor mean velocity variance in log for guidance.
niter    = 9

# Number of random realizations per iteration.
# At least 50, 100 is safer? They get faster over time.
nreal    = 100

# Number of Voronoi cells per realization.
# Keep an eye on avg cell width and rays per cell to decide on something reasonable.
# Only cells with >= min_rays_per_cell will be used.
# Bigger is not always better!
nvoronoi = 600

# How many raypaths have to travel through a cell for it to be counted.
# Will depend on your experiment (default is 3)
min_rays_per_cell = 5

# Number of Voronoi cells distributed using k-medians clustering.
#  Increase this to better handle sharper velocity gradients.
#  This is in PERCENT of NVORONOI (default 15%)
#  Increases increases meshing time by ~1.5x
kvoronoi = 10

# Number of raypath points used in k-medians clustering.
#  This controls how well the raypath distribution is sampled
#  This is in PERCENT OF ARRIVALS (default 15%)
#  Increases increases meshing time linearly ~ 1x
#  Note this only applies for non-adaptive cells
k_medians_percent = 10

# Horizontal-to-vertical ratio for stretching Voronoi cells.
# A ratio of max distance / max raypath depth seems sensible (e.g., for shallow events 150/25 = 6?).
# If your data is ringy or noisy, increasing this value can help.
#  It is also a good means to increase vertical resolution (along with paretos_alpha)
hvr = 6

# Scale factor alpha in Pareto's distribution from which cell depths are sampled.
# Higher value yields greater concentration of cells near the surface.
# Set this to 0 if you want to draw from a uniform distribution.
# Roughly scales linear, meaning a value of 3 will have about 3x more cells at the surface than bottom
paretos_alpha = 2

# Should Voronoi cells be generated adaptively and to what degree?
# This controls how far (%) cells will be adjusted based on raypath density and velocity gradients
# Recommended! (Default is 0.6, disabled/purely random if set to 0)
adaptive_data_weight = 0.6

# What should the adaptive algorithm prioritize? This sets the relative weight between
#   raypath density and HORIZONTAL velocity model gradients. Values > 0.5 will place more adaptive
#   cells in areas of higher raypath density, < 0.5 prefer higher velocity gradients.
# Default is 0.5, an even mix. Only applies if adaptive_data_weight > 0.
density_to_gradient_weight = 0.5

# Model iterations are computed via stacking all realizations and taking
#  the MEAN value of the TRIMMED stack by default.
# This can cause issues for heavily skewed distributions, or for highly variable/insufficient realizations
# Users can choose "mean" (default) or "median"
stack_type = mean
# Users may also optionally trim stack_trim_percent (%) of the lowest and highest values in
#  the stack (both 3D and 1D), mostly benefiting the "mean" calculation. By default it is 20%
#  i.e. 40% of the stack comprising the highest and lowest values are cut prior to taking mean or median
stack_trim_percent = 20

# Number of events & arrivals to sample per realization.
# Confirm good numbers with average rays per cell. 
# Ballpark: events 40% of total and arrivals 20% of total? (per phase!)
# Note that many arrivals may be QC'd out via min residual AND tukey fencing, 
#   so best to juice narrival a bit higher than your desired amount
nevent = 3000
narrival = 20000

# alternatively, can express via % of total available (per run!)
#   if set, over-rides nevent and narrival
#     recommended!
nevent_percent = 50
# note that narrivals will be a ratio of the remaining arrivals pre-decimated
#   so from your total you divide /2 for P/S then multiply by nevent_percent
narrival_percent = 40

# Threshold for # arrivals for each event.
#    Events with fewer than min_arrival arrivals are removed at the start.
# If using a "cleaned" catalog, may need to lower this
min_narrival = 13

# Remove arrivals by distance and depth (km)
# Does not apply to events with an depth greater than cutoff_depth (km, default 50)
#    this allows, for example, deep events to pass regardless of distance
cutoff_depth=50
# Recommend min_dist at 2x the model resolution and max_dist to be about 150 km or less
#    to avoid refracted waves and reflections
min_dist=7
max_dist=150
#min_depth=0
#max_depth=15

# Remove arrivals & events if residuals (seconds) too large-- HARD LIMIT!
#    Don't set TOO tightly.. things may improve!
#    Remember, Tukey fencing will remove relative arrival residuals per realization
#     and you can also downweigh by residual (solver_weight params below)
max_arrival_residual = 1.1
max_event_residual = 1.3

# Multiplicative factor for arrivals outlier removal using Tukey fences
# Applied per realization to the selected samples
#   Smaller values = stricter QC
#   ex IQR = 1.349 x std. 
#     so, IQR 0.741 = 1.0 std, IQR 1.0 = 2.0 std, IQR 1.5 = 2.7 std, IQR 3.0 = 4.7 std
outlier_removal_factor = 1.3

# Solver weight relative to arrival residuals
# Applied to the observations in the G matrix calc in compute_model_update
# Options are: huber (default) or linear 
solver_weight_method = huber
# Optionally define the weight tuning. If < 0, use default from method
#  "huber": down-weigh beyond MAD threshold (default 1.5), weights approach but never reach 0
#  "linear": down-weigh beyond threshold percentile (default 80%), then linear decay to 0
solver_weight_tuning = -1
# The weight scheme increases "blends" linearly with each iteration so that early iterations are not/less affected. Scale 0-1
#  Set both to 0 to skip weighing entirely, or both to the same value for constant weighting (not advised!)
#   default start = 0.0, default end = 0.8 (ie maximum 80% weight scaling at the last iteration)
solver_weight_start=0.0
solver_weight_end=0.8


# Remove events that have migrated beyond these limits. Only checks first 3 iterations.
#   lat/lon in degrees, depth in km, time in seconds
# May be good to restart with these relocated events pre-removed from your starting catalog
max_dlat = 0.15
max_dlon = 0.15
max_ddepth = 25
max_dtime = 2.0

# Maximum allowed depth in km.
# If set, assign hard boundaries that events cannot migrate past
#    to avoid "runaway relocations" in poorly resolved areas
# min_depth = 0
max_depth = 20

# Set which phases to run, and in what order
#   Default is P and then S, but world's your oyster
phase_order = P,S

# The following parameters (atol, btol, maxiter, conlim, and damp)
#   are passed through directly to scipy.sparse.linalg.lsmr().
# You may need to adjust based on the number of LSMR iterations you are seeing (hopefully < 300)
atol = 1e-4
btol = 0.01

# You shouldn't see anything close to 500 hopefully...
maxiter = 500
# Stops iterating if ||G^-1||||G|| (condition number estimate, conda) exceeds this
conlim = 500

# If damp is < 0, will compute a reasonable 3D estimate per iteration (recommended!) 
# Otherwise you can eyeball this and set a consistent factor for all iterations
damp = -1



[relocate]
# Set the "method" parameter here to "linear" for linearized relocation
#   or "DE" for differential evolution relocation.
# STRONGLY recommend using DE unless you have a specific reason not to
method = DE


[de_relocation]
# This section is only needed if the "method" parameter of the "relocate" section is set to "DE".
# These set the absolute maximum Differential Evolution search range for a new event.
# Keeping these values smaller increases both speed and stability, but will require more iterations 
depth_min = -1.5
dlat = 0.03
dlon = 0.03
ddepth = 3
dtime = 0.25


[linearized_relocation]
# This section is only needed if the "method" parameter of the "relocate" section is set to "linear".
atol = 1e-4
btol = 0.005
conlim = 200
damp = 0.02
maxiter = 300