search

hvasbath / beat

Bayesian Earthquake Analysis Tool
GNU General Public License v3.0
133 stars 43 forks source link

Examples in document out of date? #129

Closed sangchengfang closed 1 month ago

sangchengfang commented 3 months ago

I tried to learn BEAT following the examples in the document (https://pyrocko.org/beat/docs/current/examples/index.html) today, however, I ran into several issues.

For example, when I tried to run

beat build_gfs FullMT --datatypes='seismic' --force --execute

for example 1, I got the following output

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
beat         - INFO     Creating one global Green's Function store, which is being used by all stations!
beat         - INFO     Store name: AqabaMT
heart        - INFO     Using custom model from config file
heart        - INFO     Station AqabaMT
heart        - INFO     ---------------------
heart        - WARNING  Minimum grid distance is below zero. Setting it to zero!
heart        - INFO     Creating Store at /home/scf/master/research/code/Inversion/beat/GFs/AqabaMT_ak135_2.000Hz_0
heart        - INFO     Using custom model from config file
heart        - INFO     Receiver and source site structures have to be identical as distance and ray depth not high enough for common receiver depth!
heart        - INFO     Using modelling code: qseis with version: 2006a
heart        - INFO     Filling store ...
pyrocko.gf.store - INFO     making "phase" table for phasegroup "any_P"
pyrocko.spit - INFO     at level  0:    0.0% covered,      1 cell
pyrocko.spit - INFO     at level  1:   50.0% covered,      3 cells
pyrocko.spit - INFO     at level  2:   75.0% covered,      5 cells
pyrocko.spit - INFO     at level  3:   87.5% covered,      7 cells
pyrocko.spit - INFO     at level  4:   87.5% covered,      9 cells
pyrocko.spit - INFO     at level  5:   96.9% covered,     13 cells
pyrocko.spit - INFO     at level  6:   98.4% covered,     15 cells
pyrocko.spit - INFO     at level  7:  100.0% covered,     17 cells
pyrocko.fomosto.qseis - INFO     Starting block 1 / 11
pyrocko.fomosto.qseis - INFO     Starting block 2 / 11
pyrocko.fomosto.qseis - INFO     Starting block 3 / 11
pyrocko.fomosto.qseis - INFO     Starting block 4 / 11
pyrocko.fomosto.qseis - INFO     Starting block 5 / 11
pyrocko.fomosto.qseis - WARNING  qseis emitted something via stderr:

At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7f9f4f70e179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7f9f4f70e844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x55bae70d7649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x55bae70dd4e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x55bae70d328e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78

pyrocko.fomosto.qseis - WARNING  not removing temporary directory: /tmp/qseisrun-mpubqa_8
pyrocko.fomosto.qseis - INFO     Starting block 6 / 11
pyrocko.fomosto.qseis - WARNING  qseis emitted something via stderr:

At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7f234355d179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7f234355d844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x563128ae5649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x563128aeb4e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x563128ae128e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78

pyrocko.fomosto.qseis - INFO     Starting block 7 / 11
pyrocko.fomosto.qseis - WARNING  qseis emitted something via stderr:

At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7f46fc358179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7f46fc358844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x55df45013649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x55df450194e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x55df4500f28e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78

pyrocko.fomosto.qseis - WARNING  qseis emitted something via stderr:

At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7f54f9598179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7f54f9598844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x556bde08d649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x556bde0934e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x556bde08928e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78

pyrocko.fomosto.qseis - WARNING  qseis emitted something via stderr:

At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7fe72fdf8179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7fe72fdf8844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x5582d5000649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x5582d50064e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x5582d4ffc28e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78

pyrocko.fomosto.qseis - WARNING  qseis emitted something via stderr:

At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7f79eade1179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7f79eade1844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x5616912f7649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x5616912fd4e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x5616912f328e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
pyrocko.fomosto.qseis - WARNING  qseis emitted something via stderr:

At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7ff634a4a179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7ff634a4a844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x56197aaec649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x56197aaf24e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x56197aae828e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78

    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1475, in command_build_gfs
    heart.seis_construct_gf(
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/heart.py", line 2318, in seis_construct_gf
    backend_builders[sf.code](store_dir, nworkers=sf.nworkers, force=force)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/fomosto/qseis.py", line 1139, in build
    return QSeisGFBuilder.build(
           ^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/gf/builder.py", line 183, in build
    for x in parimap(
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/parimap.py", line 171, in parimap
    raise exc
pyrocko.fomosto.qseis.QSeisError: ===== begin qseis input =====
# autogenerated QSEIS input by qseis.py
#
# This is the input file of FORTRAN77 program "qseis06" for calculation of
# synthetic seismograms based on a layered halfspace earth model.
#
# by
# Rongjiang  Wang <wang@gfz-potsdam.de>
# GeoForschungsZentrum Potsdam
# Telegrafenberg, D-14473 Potsdam, Germany
#
# Last modified: Potsdam, Nov., 2006
#
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# If not specified, SI Unit System is used overall!
#
# Coordinate systems:
# cylindrical (z,r,t) with z = downward,
#                          r = from source outward,
#                          t = azmuth angle from north to east;
# cartesian (x,y,z) with   x = north,
#                          y = east,
#                          z = downward;
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
#
#       SOURCE PARAMETERS
#       =================
# 1. source depth [km]
#------------------------------------------------------------------------------
 8.000000e+00                    |dble: source_depth;
#------------------------------------------------------------------------------
#
#       RECEIVER PARAMETERS
#       ===================
# 1. receiver depth [km]
# 2. switch for distance sampling role (1/0 = equidistant/irregular); switch
#    for unit used (1/0 = km/deg)
# 3. number of distance samples
# 4. if equidistant, then start and end trace distance (> 0); else distance
#    list (please order the receiver distances from small to large)
# 5. (reduced) time begin [sec] & length of time window [sec], number of time
#    samples (<= 2*nfmax in qsglobal.h)
# 6. switch for unit of the following time reduction parameter: 1 = velocity
#    [km/sec], 0 = slowness [sec/deg]; time reduction parameter
#------------------------------------------------------------------------------
 0.000000e+00                         |dble: receiver_depth;
 0  1   |int: sw_equidistant, sw_d_unit;
 101                            |int: no_distances;
 0.000000e+00 1.000000e+00 2.000000e+00 3.000000e+00 4.000000e+00 5.000000e+00 6.000000e+00 7.000000e+00 8.000000e+00 9.000000e+00 1.000000e+01 1.100000e+01 1.200000e+01 1.300000e+01 1.400000e+01 1.500000e+01 1.600000e+01 1.700000e+01 1.800000e+01 1.900000e+01 2.000000e+01 2.100000e+01 2.200000e+01 2.300000e+01 2.400000e+01 2.500000e+01 2.600000e+01 2.700000e+01 2.800000e+01 2.900000e+01 3.000000e+01 3.100000e+01 3.200000e+01 3.300000e+01 3.400000e+01 3.500000e+01 3.600000e+01 3.700000e+01 3.800000e+01 3.900000e+01 4.000000e+01 4.100000e+01 4.200000e+01 4.300000e+01 4.400000e+01 4.500000e+01 4.600000e+01 4.700000e+01 4.800000e+01 4.900000e+01 5.000000e+01 5.100000e+01 5.200000e+01 5.300000e+01 5.400000e+01 5.500000e+01 5.600000e+01 5.700000e+01 5.800000e+01 5.900000e+01 6.000000e+01 6.100000e+01 6.200000e+01 6.300000e+01 6.400000e+01 6.500000e+01 6.600000e+01 6.700000e+01 6.800000e+01 6.900000e+01 7.000000e+01 7.100000e+01 7.200000e+01 7.300000e+01 7.400000e+01 7.500000e+01 7.600000e+01 7.700000e+01 7.800000e+01 7.900000e+01 8.000000e+01 8.100000e+01 8.200000e+01 8.300000e+01 8.400000e+01 8.500000e+01 8.600000e+01 8.700000e+01 8.800000e+01 8.900000e+01 9.000000e+01 9.100000e+01 9.200000e+01 9.300000e+01 9.400000e+01 9.500000e+01 9.600000e+01 9.700000e+01 9.800000e+01 9.900000e+01 1.000000e+03                          |dble: d_1,d_n; or d_1,d_2, ...(no comments in between!);
 -6.500000e+01 5.115000e+02 1024  |dble: t_start,t_window; int: no_t_samples;
 1 0.000000e+00  |int: sw_t_reduce; dble: t_reduce;
#------------------------------------------------------------------------------
#
#       WAVENUMBER INTEGRATION PARAMETERS
#       =================================
# 1. select slowness integration algorithm (0 = suggested for full wave-field
#    modelling; 1 or 2 = suggested when using a slowness window with narrow
#    taper range - a technique for suppressing space-domain aliasing);
# 2. 4 parameters for low and high slowness (Note 1) cut-offs [s/km] with
#    tapering: 0 < slw1 < slw2 defining cosine taper at the lower end, and 0 <
#    slw3 < slw4 defining the cosine taper at the higher end. default values
#    will be used in case of inconsistent input of the cut-offs (possibly with
#    much more computational effort);
# 3. parameter for sampling rate of the wavenumber integration (1 = sampled
#    with the spatial Nyquist frequency, 2 = sampled with twice higher than
#    the Nyquist, and so on: the larger this parameter, the smaller the space-
#    domain aliasing effect, but also the more computation effort);
# 4. the factor for suppressing time domain aliasing (> 0 and <= 1) (Note 2).
#------------------------------------------------------------------------------
 0                    |int: sw_algorithm;
 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00             |dble: slw(1-4);
 2.500000e+00             |dble: sample_rate;
 1.000000e-01     |dble: supp_factor;
#------------------------------------------------------------------------------
#
#               OPTIONS FOR PARTIAL SOLUTIONS
#       (only applied to the source-site structure)
#           ===========================================
#
# 1. switch for filtering free surface effects (0 = with free surface, i.e.,
#    do not select this filter; 1 = without free surface; 2 = without free
#    surface but with correction on amplitude and wave form. Note switch 2
#    can only be used for receivers at the surface)
# 2. switch for filtering waves with a shallow penetration depth (concerning
#    their whole trace from source to receiver), penetration depth limit [km]
#
#    if this option is selected, waves whose travel path never exceeds the
#    given depth limit will be filtered ("seismic nuting"). the condition for
#    selecting this filter is that the given shallow path depth limit should
#    be larger than both source and receiver depth.
#
# 3. number of depth ranges where the following selected up/down-sp2oing P or
#    SV waves should be filtered
# 4. the 1. depth range: upper and lower depth [km], switch for filtering P
#    or SV wave in this depth range:
#
#    switch no:              1      2        3       4         other
#    filtered phase:         P(up)  P(down)  SV(up)  SV(down)  Error
#
# 5. the 2. ...
#
#    The partial solution options are useful tools to increase the numerical
#    significance of desired wave phases. Especially when the desired phases
#    are smaller than the undesired phases, these options should be selected
#    and carefully combined.
#------------------------------------------------------------------------------
 0                  |int: isurf;
 0 0.000000e+00  |int: sw_path_filter; dble:shallow_depth_limit;
 0 
#------------------------------------------------------------------------------
#
#       SOURCE TIME FUNCTION (WAVELET) PARAMETERS (Note 3)
#       ==================================================
# 1. wavelet duration [unit = time sample rather than sec!], that is about
#    equal to the half-amplitude cut-off period of the wavelet (> 0. if <= 0,
#    then default value = 2 time samples will be used), and switch for the
#    wavelet form (0 = user's own wavelet; 1 = default wavelet: normalized
#    square half-sinusoid for simulating a physical delta impulse; 2 = tapered
#    Heaviside wavelet, i.e. integral of wavelet 1)
# 2. IF user's own wavelet is selected, then number of the wavelet time samples
#    (<= 1024), and followed by
# 3. equidistant wavelet time samples
# 4  ...(continue) (! no comment lines allowed between the time sample list!)
#    IF default, delete line 2, 3, 4 ... or comment them out!
#------------------------------------------------------------------------------
 1.000000e-03 2
#------------------------------------------------------------------------------
#
#        FILTER PARAMETERS OF RECEIVERS (SEISMOMETERS OR HYDROPHONES)
#        ============================================================
# 1. constant coefficient (normalization factor)
# 2. number of roots (<= nrootmax in qsglobal.h)
# 3. list of the root positions in the complex format (Re,Im). If no roots,
#    comment out this line
# 4. number of poles (<= npolemax in qsglobal.h)
# 5. list of the pole positions in the complex format (Re,Im). If no poles,
#    comment out this line
#------------------------------------------------------------------------------
 (1.0,0.0)
0
#
0
#------------------------------------------------------------------------------
#
#       OUTPUT FILES FOR GREEN'S FUNCTIONS (Note 4)
#       ===========================================
# 1. selections of source types (yes/no = 1/0)
# 2. file names of Green's functions (please give the names without extensions,
#    which will be appended by the program automatically: *.tz, *.tr, *.tt
#    and *.tv are for the vertical, radial, tangential, and volume change (for
#    hydrophones) components, respectively)
#------------------------------------------------------------------------------
#  explosion   strike-slip dip-slip   clvd       single_f_v  single_f_h
#------------------------------------------------------------------------------
 1 1 1 1 0 0
 'ex' 'ss' 'ds' 'cl' 'fz' 'fh'
#------------------------------------------------------------------------------
#       OUTPUT FILES FOR AN ARBITRARY POINT DISLOCATION SOURCE
#               (for applications to earthquakes)
#       ======================================================
# 1. selection (0 = not selected; 1 or 2 = selected), if (selection = 1), then
#    the 6 moment tensor elements [N*m]: Mxx, Myy, Mzz, Mxy, Myz, Mzx (x is
#    northward, y is eastward and z is downard); else if (selection = 2), then
#    Mis [N*m] = isotropic moment part = (MT+MN+MP)/3, Mcl = CLVD moment part
#    = (2/3)(MT+MP-2*MN), Mdc = double-couple moment part = MT-MN, Strike [deg],
#    Dip [deg] and Rake [deg].
#
#    Note: to use this option, the Green's functions above should be computed
#          (selection = 1) if they do not exist already.
#
#                 north(x)
#                  /
#                 /\ strike
#                *----------------------->  east(y)
#                |\                       #                |-\                       #                |  \     fault plane       #                |90 \                       #                |-dip\                       #                |     \                       #                |      \                       #           downward(z)  \-----------------------\
#
# 2. switch for azimuth distribution of the stations (0 = uniform azimuth,
#    else = irregular azimuth angles)
# 3. list of the azimuth angles [deg] for all stations given above (if the
#    uniform azimuth is selected, then only one azimuth angle is required)
#
#------------------------------------------------------------------------------
#     Mis        Mcl        Mdc        Strike     Dip        Rake      File
#------------------------------------------------------------------------------
#  2   0.00       1.00       6.0E+19    120.0      30.0       25.0      'seis'
#------------------------------------------------------------------------------
#     Mxx        Myy        Mzz        Mxy        Myz        Mzx       File
#------------------------------------------------------------------------------
1    1.000000e+00    0.000000e+00    0.000000e+00    1.000000e+00    0.000000e+00    0.000000e+00  'seis'
1
0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#------------------------------------------------------------------------------
#
#       GLOBAL MODEL PARAMETERS (Note 5)
#       ================================
# 1. switch for flat-earth-transform
# 2. gradient resolution [%] of vp, vs, and ro (density), if <= 0, then default
#    values (depending on wave length at cut-off frequency) will be used
#------------------------------------------------------------------------------
 0     |int: sw_flat_earth_transform;
 0.000000e+00 0.000000e+00 0.000000e+00   |dble: vp_res, vs_res, ro_res;
#------------------------------------------------------------------------------
#
#                       LAYERED EARTH MODEL
#       (SHALLOW SOURCE + UNIFORM DEEP SOURCE/RECEIVER STRUCTURE)
#       =========================================================
# 1. number of data lines of the layered model (source site)
#------------------------------------------------------------------------------
 43                   |int: no_model_lines;
#------------------------------------------------------------------------------
#
#       MULTILAYERED MODEL PARAMETERS (source site)
#       ===========================================
# no  depth[km]  vp[km/s]  vs[km/s]  ro[g/cm^3] qp      qs
#------------------------------------------------------------------------------
1 0.000000e+00 5.510000e+00 3.100000e+00 2.600000e+00 1.264000e+03 6.000000e+02
2 7.200000e+00 5.510000e+00 3.100000e+00 2.600000e+00 1.264000e+03 6.000000e+02
3 7.200000e+00 6.230000e+00 3.600000e+00 2.800000e+00 1.283000e+03 6.000000e+02
4 2.164000e+01 6.230000e+00 3.600000e+00 2.800000e+00 1.283000e+03 6.000000e+02
5 2.164000e+01 8.036000e+00 4.486000e+00 3.630000e+00 9.531000e+02 3.960000e+02
6 8.000000e+01 8.040000e+00 4.481000e+00 3.500000e+00 1.008000e+03 4.176000e+02
7 8.000000e+01 8.045000e+00 4.490000e+00 3.502000e+00 1.820000e+02 7.560000e+01
8 1.200000e+02 8.050000e+00 4.500000e+00 3.427000e+00 1.826000e+02 7.606000e+01
9 1.200000e+02 8.050000e+00 4.500000e+00 3.427000e+00 1.826000e+02 7.606000e+01
10 1.650000e+02 8.175000e+00 4.509000e+00 3.371000e+00 1.887000e+02 7.655000e+01
11 2.100000e+02 8.301000e+00 4.518000e+00 3.324000e+00 2.010000e+02 7.940000e+01
12 2.100000e+02 8.300000e+00 4.519000e+00 3.323000e+00 3.382000e+02 1.337000e+02
13 3.000000e+02 8.628000e+00 4.679000e+00 3.401000e+00 3.536000e+02 1.387000e+02
14 4.100000e+02 9.030000e+00 4.870000e+00 3.506000e+00 3.775000e+02 1.465000e+02
15 4.100000e+02 9.360000e+00 5.080000e+00 3.929000e+00 4.141000e+02 1.627000e+02
16 6.600000e+02 1.020000e+01 5.611000e+00 3.918000e+00 4.285000e+02 1.729000e+02
17 6.600000e+02 1.079000e+01 5.965000e+00 4.240000e+00 1.349000e+03 5.495000e+02
18 7.640000e+02 1.107000e+01 6.215000e+00 4.359000e+00 1.276000e+03 5.371000e+02
19 8.491000e+02 1.121000e+01 6.272000e+00 4.463000e+00 1.263000e+03 5.273000e+02
20 1.038000e+03 1.152000e+01 6.407000e+00 4.616000e+00 1.230000e+03 5.069000e+02
21 1.227000e+03 1.181000e+01 6.527000e+00 4.710000e+00 1.198000e+03 4.880000e+02
22 1.416000e+03 1.208000e+01 6.636000e+00 4.812000e+00 1.168000e+03 4.704000e+02
23 1.605000e+03 1.233000e+01 6.735000e+00 4.909000e+00 1.141000e+03 4.544000e+02
24 1.795000e+03 1.255000e+01 6.827000e+00 5.004000e+00 1.111000e+03 4.377000e+02
25 1.984000e+03 1.278000e+01 6.913000e+00 5.096000e+00 1.065000e+03 4.156000e+02
26 2.173000e+03 1.300000e+01 6.997000e+00 5.186000e+00 1.037000e+03 4.006000e+02
27 2.362000e+03 1.321000e+01 7.079000e+00 5.273000e+00 9.954000e+02 3.812000e+02
28 2.551000e+03 1.343000e+01 7.162000e+00 5.357000e+00 9.681000e+02 3.673000e+02
29 2.740000e+03 1.365000e+01 7.248000e+00 5.439000e+00 9.325000e+02 3.506000e+02
30 2.740000e+03 1.365000e+01 7.248000e+00 5.693000e+00 7.227000e+02 2.717000e+02
31 2.790000e+03 1.365000e+01 7.259000e+00 5.720000e+00 7.269000e+02 2.740000e+02
32 2.839000e+03 1.366000e+01 7.270000e+00 5.746000e+00 7.251000e+02 2.740000e+02
33 2.892000e+03 1.366000e+01 7.282000e+00 5.772000e+00 7.231000e+02 2.740000e+02
34 2.892000e+03 7.972000e+00 0.000000e+00 9.928000e+00 5.782000e+04 0.000000e+00
35 3.344000e+03 8.744000e+00 0.000000e+00 1.059000e+01 5.782000e+04 0.000000e+00
36 3.797000e+03 9.338000e+00 0.000000e+00 1.112000e+01 5.782000e+04 0.000000e+00
37 4.249000e+03 9.759000e+00 0.000000e+00 1.155000e+01 5.782000e+04 0.000000e+00
38 4.702000e+03 1.009000e+01 0.000000e+00 1.188000e+01 5.782000e+04 0.000000e+00
39 5.154000e+03 1.032000e+01 0.000000e+00 1.214000e+01 5.782000e+04 0.000000e+00
40 5.154000e+03 1.104000e+01 3.505000e+00 1.270000e+01 6.328000e+02 8.503000e+01
41 5.397000e+03 1.112000e+01 3.567000e+00 1.282000e+01 6.192000e+02 8.503000e+01
42 5.884000e+03 1.123000e+01 3.647000e+00 1.297000e+01 6.050000e+02 8.503000e+01
43 6.370999e+03 1.127000e+01 3.672000e+00 1.302000e+01 6.006000e+02 8.503000e+01
#------------------------------------------------------------------------------
#
#                 LAYERED EARTH MODEL
#       (ONLY THE SHALLOW RECEIVER STRUCTURE)
#       =====================================
# 1. number of data lines of the layered model
#
#    Note: if the number = 0, then the receiver site is the same as the
#          source site, else different receiver-site structure is considered.
#          please be sure that the lowest interface of the receiver-site
#          structure given given below can be found within the source-site
#          structure, too.
#
#------------------------------------------------------------------------------
 0                               |int: no_model_lines;
#------------------------------------------------------------------------------
#
#       MULTILAYERED MODEL PARAMETERS (shallow receiver-site structure)
#       ===============================================================
# no  depth[km]    vp[km/s]    vs[km/s]   ro[g/cm^3]   qp      qs
#------------------------------------------------------------------------------
# no receiver side model
#---------------------------------end of all inputs----------------------------

Note 1:

The slowness is defined by inverse value of apparent wave velocity = sin(i)/v
with i = incident angle and v = true wave velocity.

Note 2:

The suppression of the time domain aliasing is achieved by using the complex
frequency technique. The suppression factor should be a value between 0 and 1.
If this factor is set to 0.1, for example, the aliasing phase at the reduced
time begin is suppressed to 10%.

Note 3:

The default basic wavelet function (option 1) is (2/tau)*sin^2(pi*t/tau),
for 0 < t < tau, simulating physical delta impuls. Its half-amplitude cut-off
frequency is 1/tau. To avoid high-frequency noise, tau should not be smaller
than 4-5 time samples.

Note 4:

  Double-Couple   m11/ m22/ m33/ m12/ m23/ m31  Azimuth_Factor_(tz,tr,tv)/(tt)
  ============================================================================
  explosion       1.0/ 1.0/ 1.0/ -- / -- / --       1.0         /   0.0
  strike-slip     -- / -- / -- / 1.0/ -- / --       sin(2*azi)  /   cos(2*azi)
                  1.0/-1.0/ -- / -- / -- / --       cos(2*azi)  /  -sin(2*azi)
  dip-slip        -- / -- / -- / -- / -- / 1.0      cos(azi)    /   sin(azi)
                  -- / -- / -- / -- / 1.0/ --       sin(azi)    /  -cos(azi)
  clvd           -0.5/-0.5/ 1.0/ -- / -- / --       1.0         /   0.0
  ============================================================================
  Single-Force    fx / fy / fz                  Azimuth_Factor_(tz,tr,tv)/(tt)
  ============================================================================
  fz              -- / -- / 1.0                        1.0      /   0.0
  fx              1.0/ -- / --                         cos(azi) /   sin(azi)
  fy              -- / 1.0/ --                         sin(azi) /  -cos(azi)
  ============================================================================

Note 5:

Layers with a constant gradient will be discretized with a number of homogeneous
sublayers. The gradient resolutions are then used to determine the maximum
allowed thickness of the sublayers. If the resolutions of Vp, Vs and Rho
(density) require different thicknesses, the smallest is first chosen. If this
is even smaller than 1% of the characteristic wavelength, then the latter is
taken finally for the sublayer thickness.
===== end qseis input =====
===== begin qseis output =====
 ######################################################
 #                                                    #
 #               Welcome to the program               #
 #                                                    #
 #                                                    #
 #        QQQ     SSSS    EEEEE    III     SSSS       #
 #       Q   Q   S        E         I     S           #
 #       Q Q Q    SSS     EEEE      I      SSS        #
 #       Q  QQ       S    E         I         S       #
 #        QQQQ   SSSS     EEEEE    III    SSSS        #
 #                                                    #
 #                  (Version 2006a)                   #
 #                                                    #
 #                                                    #
 #                      by                            #
 #                 Rongjiang Wang                     #
 #              (wang@gfz-potsdam.de)                 #
 #                                                    #
 #           GeoForschungsZentrum Potsdam             #
 #           Last modified: October, 2006             #
 ######################################################

 the input data file is ===== end qseis output =====
===== begin qseis error =====
At line 133 of file qsgetinp.f (unit = 10, file = 'input')
Fortran runtime error: Bad real number in item 2 of list input

Error termination. Backtrace:
#0  0x7f9f4f70e179 in read_real
        at ../../../libgfortran/io/list_read.c:2080
#1  0x7f9f4f70e844 in list_formatted_read_scalar
        at ../../../libgfortran/io/list_read.c:2236
#2  0x55bae70d7649 in qsgetinp_
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsgetinp.f:133
#3  0x55bae70dd4e3 in qseis
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:44
#4  0x55bae70d328e in main
        at /home/scf/master/research/code/Inversion/beat/fomosto-qseis/src/qsmain.f:78
===== end qseis error =====
qseis had a non-zero exit state: 2
qseis has been invoked as "fomosto_qseis2006a"
in the directory /tmp/qseisrun-mpubqa_8

When I tried to run

beat sample Laquila_dc --hypers

for example 2, I got the following output

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using DCSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila_dc/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

models       - INFO     ... Building Hyper model ...

models       - INFO     Initialised hyperparameter h_any_P_0_Z with size 35 
models       - INFO     Initialized 35 hyperparameters in total!
seismic      - INFO     Retrieving seismic data-covariances with structure "variance" for any_P_0 ...
covariance   - INFO     Variance estimate of IU, TSUM, 0, Z = 7.21748e-15
covariance   - INFO     Variance estimate of IU, RCBR, 0, Z = 1.1586e-14
covariance   - INFO     Variance estimate of CN, RES, 0, Z = 3.08076e-14
covariance   - INFO     Variance estimate of IC, KMI, 0, Z = 1.66608e-15
covariance   - INFO     Variance estimate of IU, SSPA, 0, Z = 2.74838e-14
covariance   - INFO     Variance estimate of MS, BTDF, 0, Z = 1.27172e-15
covariance   - INFO     Variance estimate of CU, ANWB, 0, Z = 7.10948e-15
covariance   - INFO     Variance estimate of G, MBO, 0, Z = 3.15833e-14
covariance   - INFO     Variance estimate of IU, TEIG, 0, Z = 3.89182e-15
covariance   - INFO     Variance estimate of IU, KBL, 0, Z = 1.22667e-15
covariance   - INFO     Variance estimate of II, MSEY, 0, Z = 1.05611e-14
covariance   - INFO     Variance estimate of GT, DBIC, 0, Z = 4.86832e-15
covariance   - INFO     Variance estimate of IC, HIA, 0, Z = 2.51146e-14
covariance   - INFO     Variance estimate of IU, YAK, 0, Z = 2.31697e-15
covariance   - INFO     Variance estimate of IC, QIZ, 0, Z = 2.81997e-15
covariance   - INFO     Variance estimate of GE, KBS, 0, Z = 8.1559e-15
covariance   - INFO     Variance estimate of IU, LSZ, 0, Z = 6.80444e-15
covariance   - INFO     Variance estimate of II, TLY, 0, Z = 8.67367e-16
covariance   - INFO     Variance estimate of IU, INCN, 0, Z = 2.95215e-15
covariance   - INFO     Variance estimate of IU, DWPF, 0, Z = 4.46816e-15
covariance   - INFO     Variance estimate of IC, WMQ, 0, Z = 8.53728e-14
covariance   - INFO     Variance estimate of II, AAK, 0, Z = 1.02798e-15
covariance   - INFO     Variance estimate of IC, XAN, 0, Z = 5.60621e-15
covariance   - INFO     Variance estimate of IU, KMBO, 0, Z = 2.74861e-15
covariance   - INFO     Variance estimate of II, KURK, 0, Z = 1.94683e-15
covariance   - INFO     Variance estimate of ZF, ADYE, 0, Z = 5.08499e-15
covariance   - INFO     Variance estimate of IU, SDV, 0, Z = 3.43033e-15
covariance   - INFO     Variance estimate of II, ARU, 0, Z = 3.25339e-15
covariance   - INFO     Variance estimate of IU, SAML, 0, Z = 8.99711e-15
covariance   - INFO     Variance estimate of II, ALE, 0, Z = 5.47305e-15
covariance   - INFO     Variance estimate of II, ABKT, 0, Z = 1.12205e-15
covariance   - INFO     Variance estimate of KZ, AKTO, 0, Z = 2.29098e-15
covariance   - INFO     Variance estimate of IC, BJT, 0, Z = 1.18455e-15
covariance   - INFO     Variance estimate of IC, SSE, 0, Z = 1.56772e-14
covariance   - INFO     Variance estimate of IU, CHTO, 0, Z = 9.11924e-15
seismic      - INFO     Initialising weights ...
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/heart.py", line 3056, in prepare_data
    arrival_times[i] = marker_arrival_times[target.codes]
                       ~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^
KeyError: ('IU', 'TSUM', '0', 'Z')

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1061, in command_sample
    problem = load_model(project_dir, options.mode, options.hypers)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 920, in load_model
    problem.built_hyper_model()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 291, in built_hyper_model
    composite.update_llks(point)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/seismic.py", line 520, in update_llks
    results = self.assemble_results(point, chop_bounds=["b", "c"])
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/seismic.py", line 470, in assemble_results
    syn_proc_traces, obs_proc_traces = self.get_synthetics(
                                       ^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/seismic.py", line 872, in get_synthetics
    wmap.prepare_data(
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/heart.py", line 3059, in prepare_data
    arrival_times[i] = get_phase_arrival_time(
                       ^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/heart.py", line 2571, in get_phase_arrival_time
    atime = store.t(wavename, (source.depth, dist)) + source.time
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/gf/store.py", line 1753, in t
    return timing.evaluate(self.get_phase, args, attributes=attributes)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/gf/meta.py", line 580, in evaluate
    phases = [
             ^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/gf/meta.py", line 581, in <listcomp>
    get_phase(phase_def, *extra_args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/gf/store.py", line 1622, in get_phase
    return self.get_stored_phase(phase_def)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pyrocko/gf/store.py", line 1600, in get_stored_phase
    raise NoSuchPhase(phase_id)
pyrocko.gf.store.NoSuchPhase: phase for key "any_P.phase" not found. Running "fomosto ttt" or "fomosto sat" may be needed.

when I tried to run

beat sample Laquila --hypers

for example 3, I got the following output

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

models       - INFO     ... Building Hyper model ...

geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 137, in get_context
    candidate: T | None = cls.get_contexts()[-1]
                          ~~~~~~~~~~~~~~~~~~^^^^
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 504, in __new__
    model = Model.get_context()
            ^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 142, in get_context
    raise TypeError(f"No {cls} on context stack")
TypeError: No <class 'pymc.model.core.Model'> on context stack

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 41, in init_uniform_random
    dist = Uniform(**kwargs)
           ^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 506, in __new__
    raise TypeError(
TypeError: No model on context stack, which is needed to instantiate distributions. Add variable inside a 'with model:' block, or use the '.dist' syntax for a standalone distribution.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1061, in command_sample
    problem = load_model(project_dir, options.mode, options.hypers)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 920, in load_model
    problem.built_hyper_model()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 273, in built_hyper_model
    self.init_hierarchicals()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 599, in init_hierarchicals
    composite.init_hierarchicals(self.config.problem_config)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/geodetic.py", line 395, in init_hierarchicals
    ] = init_uniform_random(kwargs)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 45, in init_uniform_random
    kwargs.pop("transform")
KeyError: 'transform'

I am trying to figure out what is wrong, but I noticed that there is a brand new release of BEAT. Maybe the examples in the document have not been updated yet? Or maybe there are some issues with my installation of BEAT and compilation of the backend of BEAT? Could you help me solve these problems? Thanks.

hvasbath commented 3 months ago

Thank you for reporting! For the first example it seems there is a compilation issue with fomosto_qseis_2006a. Can you please try to recompile that? Also there is a new version out in the meanwhile fomosto_qseis_2006b that you may want to try. You would need to set in the seismic_config.gf_config.version = "2006b" to use that.

The second example I will need to investigate- it can indeed be due to the update- I havent run the tutorial in a while- thank you for your patience. I will get back to you.

sangchengfang commented 3 months ago

Thanks for your reply and suggestions! For example 1, I did not try fomosto_qseis_2006, but set the version to 2006b and it works. I was able to run this example. However, I ran into new problems related to plotting. I can run this command

beat plot FullMT_nont waveform_fits --nensemble=100

The result is waveforms_-1_max_100_0.pdf.

When I ran

beat plot FullMT_nont  waveform_fits  --reference

I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using MTSource for 1 sources for event --- !pf.Event
lat: 29.07
lon: 34.73
time: '1995-11-22 04:15:26.200000048'
depth: 8000.0
name: '112295A'
magnitude: 7.20583885303153
region: 'ARAB REPUBLIC OF EGYPT'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: -0.43283071
  mee: 0.65741974
  mdd: -0.22458903
  mne: 0.63839719
  mnd: 0.50698292
  med: 0.02063122
  strike1: 294.09064689577235
  dip1: 77.25911269444555
  rake1: -148.53160423717637
  strike2: 196.4032557270442
  dip2: 59.39114178630978
  rake2: -14.847482440839393
  moment: 0.999999992919433
  magnitude: -6.033333335383367
duration: 22.0

seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/FullMT_nont/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 20 

seismic      - INFO     Initialising seismic wavemap for "slowest" ...
seismic      - INFO     The waveform defined in "slowest 1" config is not included in the optimization!
plotting.seismic - INFO     Drawing Waveform fits ...
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
plotting.common - INFO     Non-Toeplitz noise structure: Using TestPoint for Covariance!
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
seismic      - INFO     Retrieving seismic data-covariances with structure "non-toeplitz" for any_P_0 ...
seismic      - INFO     Initialising weights ...
seismic      - INFO     Updating data-covariances ...
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
seismic      - INFO     Retrieving seismic data-covariances with structure "non-toeplitz" for any_P_0 ...
seismic      - INFO     Not updating seismic velocity model-covariances because number of model variations is too low! < 5
seismic      - INFO     Updating weights of wavemap any_P_0
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1957, in command_plot
    plotting.plots_catalog[plot](problem, po)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/seismic.py", line 1125, in draw_seismic_fits
    event_figs = seismic_fits(problem, stage, po)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/seismic.py", line 828, in seismic_fits
    bvar_reductions = composite.get_variance_reductions(
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/seismic.py", line 609, in get_variance_reductions
    hp = get_hypervalue_from_point(point, tr, counter, hp_specific=hp_specific)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 31, in get_hypervalue_from_point
    hp = point[hp_name][counter(hp_name)]
         ~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^
IndexError: index 1 is out of bounds for axis 0 with size 1

When I ran

beat plot FullMT_nont stage_posteriors --reference --stage_number=-1 --format='png'

I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using MTSource for 1 sources for event --- !pf.Event
lat: 29.07
lon: 34.73
time: '1995-11-22 04:15:26.200000048'
depth: 8000.0
name: '112295A'
magnitude: 7.20583885303153
region: 'ARAB REPUBLIC OF EGYPT'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: -0.43283071
  mee: 0.65741974
  mdd: -0.22458903
  mne: 0.63839719
  mnd: 0.50698292
  med: 0.02063122
  strike1: 294.09064689577235
  dip1: 77.25911269444555
  rake1: -148.53160423717637
  strike2: 196.4032557270442
  dip2: 59.39114178630978
  rake2: -14.847482440839393
  moment: 0.999999992919433
  magnitude: -6.033333335383367
duration: 22.0

seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/FullMT_nont/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 20 

seismic      - INFO     Initialising seismic wavemap for "slowest" ...
seismic      - INFO     The waveform defined in "slowest 1" config is not included in the optimization!
plotting.marginals - INFO     Plotting "pdf"
config       - INFO     not solving for depth, got fixed at 8.0
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 137, in get_context
    candidate: T | None = cls.get_contexts()[-1]
                          ~~~~~~~~~~~~~~~~~~^^^^
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 504, in __new__
    model = Model.get_context()
            ^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 142, in get_context
    raise TypeError(f"No {cls} on context stack")
TypeError: No <class 'pymc.model.core.Model'> on context stack

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 41, in init_uniform_random
    dist = Uniform(**kwargs)
           ^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 506, in __new__
    raise TypeError(
TypeError: No model on context stack, which is needed to instantiate distributions. Add variable inside a 'with model:' block, or use the '.dist' syntax for a standalone distribution.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1957, in command_plot
    plotting.plots_catalog[plot](problem, po)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/marginals.py", line 880, in draw_posteriors
    problem.varnames + problem.hypernames + problem.hierarchicalnames + ["like"]
                       ^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 355, in hypernames
    self.init_hyperparams()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 421, in init_hyperparams
    hyperparams[hp_name] = init_uniform_random(kwargs)
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 45, in init_uniform_random
    kwargs.pop("transform")
KeyError: 'transform'

The command

beat plot FullMT_nont correlation_hist --reference --stage_number=-1 --format='png' --varnames='mee, med, mdd, mnn, mnd, mne, north_shift, east_shift, magnitude'

ran successfully. The result is corr_hist_-1_ref_0 The code may be running but there are some problems with the plot functions of BEAT. This is the summary of this example. summary.txt I hope this is helpful to you.

sangchengfang commented 3 months ago

I also tried example 3 and it works. This is the summary of this example. summary.txt But I was only able to run this command to plot the result

beat plot Laquila correlation_hist --format=png --varnames=east_shift,north_shift,depth,length,width,strike,dip,rake,slip

This is the result and this figure seems plausible. corr_hist_-1_max_0

The other three plotting commands give me many error messages, which look similar to the failed plot commands of example 1. So I gusse the major problems may relate to plotting of BEAT. By the way, I have installed Kite in my system.

The commands and corresponding error messages are listed here for you.

beat plot Laquila scene_fits
To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

plotting.geodetic - INFO     Drawing SAR misfits ...
backend      - INFO     Loading multitrace from /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/geometry/stage_-1
geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 137, in get_context
    candidate: T | None = cls.get_contexts()[-1]
                          ~~~~~~~~~~~~~~~~~~^^^^
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 504, in __new__
    model = Model.get_context()
            ^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 142, in get_context
    raise TypeError(f"No {cls} on context stack")
TypeError: No <class 'pymc.model.core.Model'> on context stack

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 41, in init_uniform_random
    dist = Uniform(**kwargs)
           ^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 506, in __new__
    raise TypeError(
TypeError: No model on context stack, which is needed to instantiate distributions. Add variable inside a 'with model:' block, or use the '.dist' syntax for a standalone distribution.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1957, in command_plot
    plotting.plots_catalog[plot](problem, po)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/geodetic.py", line 1147, in draw_scene_fits
    figs = scene_fits(problem, stage, po)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/geodetic.py", line 503, in scene_fits
    problem.init_hierarchicals()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 599, in init_hierarchicals
    composite.init_hierarchicals(self.config.problem_config)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/geodetic.py", line 395, in init_hierarchicals
    ] = init_uniform_random(kwargs)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 45, in init_uniform_random
    kwargs.pop("transform")
KeyError: 'transform'
beat plot Laquila scene_fits --plot_projection=latlon
To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

plotting.geodetic - INFO     Drawing SAR misfits ...
backend      - INFO     Loading multitrace from /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/geometry/stage_-1
geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 137, in get_context
    candidate: T | None = cls.get_contexts()[-1]
                          ~~~~~~~~~~~~~~~~~~^^^^
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 504, in __new__
    model = Model.get_context()
            ^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 142, in get_context
    raise TypeError(f"No {cls} on context stack")
TypeError: No <class 'pymc.model.core.Model'> on context stack

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 41, in init_uniform_random
    dist = Uniform(**kwargs)
           ^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 506, in __new__
    raise TypeError(
TypeError: No model on context stack, which is needed to instantiate distributions. Add variable inside a 'with model:' block, or use the '.dist' syntax for a standalone distribution.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1957, in command_plot
    plotting.plots_catalog[plot](problem, po)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/geodetic.py", line 1147, in draw_scene_fits
    figs = scene_fits(problem, stage, po)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/geodetic.py", line 503, in scene_fits
    problem.init_hierarchicals()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 599, in init_hierarchicals
    composite.init_hierarchicals(self.config.problem_config)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/geodetic.py", line 395, in init_hierarchicals
    ] = init_uniform_random(kwargs)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 45, in init_uniform_random
    kwargs.pop("transform")
KeyError: 'transform'
beat plot Laquila waveform_fits --nensemble=100
To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

plotting.seismic - INFO     Drawing Waveform fits ...
backend      - INFO     Loading multitrace from /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/geometry/stage_-1
geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 137, in get_context
    candidate: T | None = cls.get_contexts()[-1]
                          ~~~~~~~~~~~~~~~~~~^^^^
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 504, in __new__
    model = Model.get_context()
            ^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 142, in get_context
    raise TypeError(f"No {cls} on context stack")
TypeError: No <class 'pymc.model.core.Model'> on context stack

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 41, in init_uniform_random
    dist = Uniform(**kwargs)
           ^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 506, in __new__
    raise TypeError(
TypeError: No model on context stack, which is needed to instantiate distributions. Add variable inside a 'with model:' block, or use the '.dist' syntax for a standalone distribution.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1957, in command_plot
    plotting.plots_catalog[plot](problem, po)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/seismic.py", line 1125, in draw_seismic_fits
    event_figs = seismic_fits(problem, stage, po)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/seismic.py", line 786, in seismic_fits
    problem.init_hierarchicals()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 599, in init_hierarchicals
    composite.init_hierarchicals(self.config.problem_config)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/geodetic.py", line 395, in init_hierarchicals
    ] = init_uniform_random(kwargs)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 45, in init_uniform_random
    kwargs.pop("transform")
KeyError: 'transform'
~/master/research/code/Inversion/beat/exp_scf                                                             Py beat 22:28:55
❯ beat plot Laquila stage_posteriors --stage_number=-1
To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: /home/scf/master/research/code/Inversion/beat/exp_scf/Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

plotting.marginals - INFO     Plotting "pdf"
Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 137, in get_context
    candidate: T | None = cls.get_contexts()[-1]
                          ~~~~~~~~~~~~~~~~~~^^^^
IndexError: list index out of range

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 504, in __new__
    model = Model.get_context()
            ^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/model/core.py", line 142, in get_context
    raise TypeError(f"No {cls} on context stack")
TypeError: No <class 'pymc.model.core.Model'> on context stack

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 41, in init_uniform_random
    dist = Uniform(**kwargs)
           ^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/pymc/distributions/distribution.py", line 506, in __new__
    raise TypeError(
TypeError: No model on context stack, which is needed to instantiate distributions. Add variable inside a 'with model:' block, or use the '.dist' syntax for a standalone distribution.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/scf/miniforge3/envs/beat/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/apps/beat.py", line 1957, in command_plot
    plotting.plots_catalog[plot](problem, po)
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/plotting/marginals.py", line 880, in draw_posteriors
    problem.varnames + problem.hypernames + problem.hierarchicalnames + ["like"]
                       ^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 355, in hypernames
    self.init_hyperparams()
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/problems.py", line 421, in init_hyperparams
    hyperparams[hp_name] = init_uniform_random(kwargs)
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/scf/miniforge3/envs/beat/lib/python3.11/site-packages/beat/models/base.py", line 45, in init_uniform_random
    kwargs.pop("transform")
KeyError: 'transform'

I hope what I present here is clear and helpful to you. Thanks again for your reply.

hvasbath commented 3 months ago

The first error is related to the fact that you can plot --reference only with dataset_specific_residual_noise=False` the other stuff is related to me introducing a bug with the latest maintenance release. Will fix that and update the release. Will let you know.

hvasbath commented 3 months ago

Hi @sangchengfang can you please install the current master branch from github and try if that resolves your problems? Locally I could reproduce your errors and could resolve them. Still its good to double check,... Thanks!

sangchengfang commented 2 months ago

Hi, sorry for the reply. Thanks for your efforts to solve my problems. I have been testing the latest version of BEAT for the last few days. The previous problems have been solved. But I got new problems with the last step of example 1, i.e. the repeat of example 1 with the noise structure set to non-toeplitz.

When I ran

beat sample FullMT_nont

I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using MTSource for 1 sources for event --- !pf.Event
lat: 29.07
lon: 34.73
time: '1995-11-22 04:15:26.200000048'
depth: 8000.0
name: '112295A'
magnitude: 7.20583885303153
region: 'ARAB REPUBLIC OF EGYPT'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: -0.43283071
  mee: 0.65741974
  mdd: -0.22458903
  mne: 0.63839719
  mnd: 0.50698292
  med: 0.02063122
  strike1: 294.09064689577235
  dip1: 77.25911269444555
  rake1: -148.53160423717637
  strike2: 196.4032557270442
  dip2: 59.39114178630978
  rake2: -14.847482440839393
  moment: 0.999999992919433
  magnitude: -6.033333335383367
duration: 22.0

seismic      - INFO     Loading seismic data for event 0 from: /work/csang/exp/FullMT_nont/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 20 

seismic      - INFO     Initialising seismic wavemap for "slowest" ...
seismic      - INFO     The waveform defined in "slowest 1" config is not included in the optimization!
models       - INFO     ... Building model ...

config       - INFO     not solving for depth, got fixed at 8.0
models       - INFO     Initialised hyperparameter h_any_P_0_Z with size 10 
models       - INFO     Initialised hyperparameter h_any_P_0_T with size 10 
models       - INFO     Initialized 20 hyperparameters in total!
seismic      - INFO     Seismic optimization on: 
  duration, east_shift, magnitude, mdd, med, mee, mnd, mne, mnn, north_shift, time
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
seismic      - INFO     Retrieving seismic data-covariances with structure "non-toeplitz" for any_P_0 ...
seismic      - INFO     Initialising weights ...
seismic      - INFO     Preparing data of "any_P_0" for optimization
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
seismic      - INFO     Initializing synthesizer for "any_P_0"
seismic      - INFO     Using all sources for wavemap any_P_0 !
models       - INFO     Model building was successful! 

models       - INFO     Using "bin" backend to store samples!
models       - INFO     ... Initiate Sequential Monte Carlo ... 
 n_chains=2000, tune_interval=10, n_jobs=100, proposal_distribution: MultivariateNormal, 

metropolis   - INFO     Creating initial population for 2000 chains ...
Sampling: [duration, east_shift, h_any_P_0_T, h_any_P_0_Z, magnitude, mdd, med, mee, mnd, mne, mnn, north_shift, time]
pymc.sampling.forward - INFO     Sampling: [duration, east_shift, h_any_P_0_T, h_any_P_0_Z, magnitude, mdd, med, mee, mnd, mne, mnn, north_shift, time]
metropolis   - INFO     Compiling model graph ...
metropolis   - INFO     Initializing proposal distribution ...MultivariateNormal
metropolis   - INFO     Time for proposal covariance init: 0.022191
models       - INFO     Compilation time: 5.248945 

models.base  - INFO     ... Starting SMC ...

backend      - INFO     Removing previous sampling results ... /work/csang/exp/FullMT_nont/geometry/stage_0
backend      - INFO     Found no sampling results under /work/csang/exp/FullMT_nont/geometry/stage_0 
backend      - INFO     Init new trace!
smc          - INFO     Sample initial stage: ...
smc          - INFO     Beta: 0.000000 Stage: 0
sampler      - INFO     Initialising 2000 chain traces ...
sampler      - INFO     Serial time per sample: 0.008421
sampler      - INFO     Chunksize per worker is 20
sampler      - INFO     Data to be memory shared: 
sampler      - INFO     No data to be memshared!
sampler      - INFO     Sampling ...
sampler      - INFO     Worker 11: Finished 4 / 20 chains
...... (redundancy lines are omitted here)
sampler      - INFO     Worker 45: Finished 20 / 20 chains
parallel     - INFO     
 Feierabend! Done with the work!
backend      - INFO     Loading multitrace from /work/csang/exp/FullMT_nont/geometry/stage_0
smc          - INFO     Updating Covariances ...
Traceback (most recent call last):
  File "/home/csang/miniforge3/envs/beat_dev/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/apps/beat.py", line 1068, in command_sample
    sample(step, problem)
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/models/base.py", line 262, in sample
    smc_sample(
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/sampler/smc.py", line 494, in smc_sample
    map_pt = step.get_map_end_points()
             ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/sampler/smc.py", line 288, in get_map_end_points
    return self.bij.rmap(self.array_population[idx, :])
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/pymc/blocking.py", line 105, in rmap
    raise TypeError("`array` must be a `RaveledVars` type")
TypeError: `array` must be a `RaveledVars` type

In short, I cannot reproduce the last figure of example 1, but I can reproduce other steps of example 1. Attached is my config file for example 1. FullMT_config.zip

sangchengfang commented 2 months ago

For example 2, my results do not seem right. The beach ball and probability distributions are not good, but I did not find the problems. So I need some help to reproduce the example. Attached is my config file, the log file, and the results for example 2.

example2.zip

sangchengfang commented 2 months ago

for example 3 All steps run successfully. But when I want to plot the InSAR data fitting with

beat plot Laquila scene_fits

and 

beat plot Laquila scene_fits --plot_projection=latlon

I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: ./Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

plotting.geodetic - INFO     Drawing SAR misfits ...
backend      - INFO     Loading multitrace from ./Laquila/geometry/stage_-1
geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
geodetic     - INFO     Evaluating config for Euler Pole corrections for datasets...
geodetic     - INFO     No Euler Pole correction!
geodetic     - INFO     Evaluating config for Strain Rate corrections for datasets...
geodetic     - INFO     No Strain Rate correction!
geodetic     - INFO     Initialized 6 hierarchical parameters.
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for Laquila_dscxn ...
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for Laquila_ascxn ...
geodetic     - INFO     Not updating geodetic velocity model-covariances because number of model variations is too low! < 5
plotting.geodetic - INFO     Not plotting shaded relief for nensemble < 301.
plotting.geodetic - INFO     Loading full resolution kite scene: ./Laquila_dscxn
plotting.geodetic - WARNING  Full resolution data could not be loaded! Skipping ...
plotting.geodetic - INFO     Loading full resolution kite scene: ./Laquila_ascxn
plotting.geodetic - WARNING  Full resolution data could not be loaded! Skipping ...
plotting.common - INFO     saving figures to ./Laquila/geometry/figures/scenes_-1_max_local_0.pdf

However, the fig is empty. This is the same with InSAR data fitting plots with example 4. It is likely that there are some problems with the plotting function. Attached are the plots and my config file. example3_geometry.zip

sangchengfang commented 2 months ago

For example 4, I have a question. I can run all the steps of example 4 but the results are not as good as that in the tutorial. I tried several times but I cannot reproduce similar results. I know for nonlinear problems, the results of every execution are not the same, but I think my results are not good. Maybe there is something wrong with my config file. Are there any tricks to improve the results of the finite fault inversion? Attached is my config file and the results for example 4. example4_ffi.zip

hvasbath commented 2 months ago

for example 3 All steps run successfully. But when I want to plot the InSAR data fitting with

beat plot Laquila scene_fits

and

beat plot Laquila scene_fits --plot_projection=latlon

I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: ./Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

plotting.geodetic - INFO     Drawing SAR misfits ...
backend      - INFO     Loading multitrace from ./Laquila/geometry/stage_-1
geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
geodetic     - INFO     Evaluating config for Euler Pole corrections for datasets...
geodetic     - INFO     No Euler Pole correction!
geodetic     - INFO     Evaluating config for Strain Rate corrections for datasets...
geodetic     - INFO     No Strain Rate correction!
geodetic     - INFO     Initialized 6 hierarchical parameters.
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for Laquila_dscxn ...
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for Laquila_ascxn ...
geodetic     - INFO     Not updating geodetic velocity model-covariances because number of model variations is too low! < 5
plotting.geodetic - INFO     Not plotting shaded relief for nensemble < 301.
plotting.geodetic - INFO     Loading full resolution kite scene: ./Laquila_dscxn
plotting.geodetic - WARNING  Full resolution data could not be loaded! Skipping ...
plotting.geodetic - INFO     Loading full resolution kite scene: ./Laquila_ascxn
plotting.geodetic - WARNING  Full resolution data could not be loaded! Skipping ...
plotting.common - INFO     saving figures to ./Laquila/geometry/figures/scenes_-1_max_local_0.pdf

However, the fig is empty. This is the same with InSAR data fitting plots with example 4. It is likely that there are some problems with the plotting function. Attached are the plots and my config file. example3_geometry.zip

for example 3 All steps run successfully. But when I want to plot the InSAR data fitting with

beat plot Laquila scene_fits

and

beat plot Laquila scene_fits --plot_projection=latlon

I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using RectangularSource for 1 sources for event --- !pf.Event
lat: 42.29
lon: 13.35
time: '2009-04-06 01:32:49.190000057'
depth: 12000.0
name: '200904060132A'
magnitude: 6.343080192483292
region: 'CENTRAL ITALY'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: 1.43e+18
  mee: 1.87e+18
  mdd: -3.3e+18
  mne: 1.77e+18
  mnd: -1.43e+18
  med: 2.6900000000000003e+17
  strike1: 120.23408298515041
  dip1: 54.240869089580485
  rake1: -112.81739742081386
  strike2: 335.98575923255856
  dip2: 41.58440373860804
  rake2: -61.69749587601104
  moment: 3.6696131948749036e+18
  magnitude: 6.343080192483292
duration: 7.0

/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/pytensor/tensor/sharedvar.py:14: FutureWarning: The class `ScalarSharedVariable` has been deprecated. Use `TensorSharedVariable` instead and check for `ndim==0`.
  warnings.warn(
geodetic     - INFO     Number of geodetic datasets: 2 
geodetic     - INFO     Initialising corrections ...
heart        - INFO     Setting up Ramps correction for Laquila_dscxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_dscxn for Euler Pole
heart        - INFO     Not correcting Laquila_dscxn for Strain Rate
heart        - INFO     Setting up Ramps correction for Laquila_ascxn
heart        - INFO     Masking data for Ramps estimation!
heart        - INFO     Not correcting Laquila_ascxn for Euler Pole
heart        - INFO     Not correcting Laquila_ascxn for Strain Rate
geodetic     - INFO     Number of geodetic data points: 419 
geodetic     - INFO     Initialising geometry geodetic composite ...
seismic      - INFO     Loading seismic data for event 0 from: ./Laquila/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 35 

plotting.geodetic - INFO     Drawing SAR misfits ...
backend      - INFO     Loading multitrace from ./Laquila/geometry/stage_-1
geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
geodetic     - INFO     Evaluating config for Euler Pole corrections for datasets...
geodetic     - INFO     No Euler Pole correction!
geodetic     - INFO     Evaluating config for Strain Rate corrections for datasets...
geodetic     - INFO     No Strain Rate correction!
geodetic     - INFO     Initialized 6 hierarchical parameters.
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for Laquila_dscxn ...
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for Laquila_ascxn ...
geodetic     - INFO     Not updating geodetic velocity model-covariances because number of model variations is too low! < 5
plotting.geodetic - INFO     Not plotting shaded relief for nensemble < 301.
plotting.geodetic - INFO     Loading full resolution kite scene: ./Laquila_dscxn
plotting.geodetic - WARNING  Full resolution data could not be loaded! Skipping ...
plotting.geodetic - INFO     Loading full resolution kite scene: ./Laquila_ascxn
plotting.geodetic - WARNING  Full resolution data could not be loaded! Skipping ...
plotting.common - INFO     saving figures to ./Laquila/geometry/figures/scenes_-1_max_local_0.pdf

However, the fig is empty. This is the same with InSAR data fitting plots with example 4. It is likely that there are some problems with the plotting function. Attached are the plots and my config file. example3_geometry.zip

It tries to load the data, but was not successful- so it couldnt plot it. So please check the path and data format etc. It needs kite scenes. Did you download them? https://github.com/braunfuss/laquila_kite_container image

hvasbath commented 2 months ago

For example 2, my results do not seem right. The beach ball and probability distributions are not good, but I did not find the problems. So I need some help to reproduce the example. Attached is my config file, the log file, and the results for example 2.

example2.zip

Do you have a 72core computer? To me it looks like a problem of convergence. Meaning that some of the chains sampling in the posterior simply do not sample at the correct position yet. This can happen if the master chain is often swapping only the same workers, which can happen when the other workers do not sample at all, because you oversubscribed too much. May I ask for your mpiversion info? Which os are you using? Also the PT sampler wasnt in use much lately. So if you want to do research- I do not recommend using it- since it is not well tested in the last years.

hvasbath commented 2 months ago

Hi, sorry for the reply. Thanks for your efforts to solve my problems. I have been testing the latest version of BEAT for the last few days. The previous problems have been solved. But I got new problems with the last step of example 1, i.e. the repeat of example 1 with the noise structure set to non-toeplitz.

When I ran

beat sample FullMT_nont

I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Geometry Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

models       - INFO     Using MTSource for 1 sources for event --- !pf.Event
lat: 29.07
lon: 34.73
time: '1995-11-22 04:15:26.200000048'
depth: 8000.0
name: '112295A'
magnitude: 7.20583885303153
region: 'ARAB REPUBLIC OF EGYPT'
catalog: 'gCMT'
moment_tensor: !pf.MomentTensor
  mnn: -0.43283071
  mee: 0.65741974
  mdd: -0.22458903
  mne: 0.63839719
  mnd: 0.50698292
  med: 0.02063122
  strike1: 294.09064689577235
  dip1: 77.25911269444555
  rake1: -148.53160423717637
  strike2: 196.4032557270442
  dip2: 59.39114178630978
  rake2: -14.847482440839393
  moment: 0.999999992919433
  magnitude: -6.033333335383367
duration: 22.0

seismic      - INFO     Loading seismic data for event 0 from: /work/csang/exp/FullMT_nont/seismic_data.pkl 
seismic      - INFO     Initialising seismic wavemap for "any_P" ...
heart        - INFO     Consistent number of datasets and targets in any_P_0 wavemap!
heart        - INFO     Number of seismic datasets for wavemap: any_P_0: 20 

seismic      - INFO     Initialising seismic wavemap for "slowest" ...
seismic      - INFO     The waveform defined in "slowest 1" config is not included in the optimization!
models       - INFO     ... Building model ...

config       - INFO     not solving for depth, got fixed at 8.0
models       - INFO     Initialised hyperparameter h_any_P_0_Z with size 10 
models       - INFO     Initialised hyperparameter h_any_P_0_T with size 10 
models       - INFO     Initialized 20 hyperparameters in total!
seismic      - INFO     Seismic optimization on: 
  duration, east_shift, magnitude, mdd, med, mee, mnd, mne, mnn, north_shift, time
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
seismic      - INFO     Retrieving seismic data-covariances with structure "non-toeplitz" for any_P_0 ...
seismic      - INFO     Initialising weights ...
seismic      - INFO     Preparing data of "any_P_0" for optimization
heart        - INFO     Did not find custom arrival times.
heart        - INFO     Using theoretical arrival times for "any_P_0"
seismic      - INFO     Initializing synthesizer for "any_P_0"
seismic      - INFO     Using all sources for wavemap any_P_0 !
models       - INFO     Model building was successful! 

models       - INFO     Using "bin" backend to store samples!
models       - INFO     ... Initiate Sequential Monte Carlo ... 
 n_chains=2000, tune_interval=10, n_jobs=100, proposal_distribution: MultivariateNormal, 

metropolis   - INFO     Creating initial population for 2000 chains ...
Sampling: [duration, east_shift, h_any_P_0_T, h_any_P_0_Z, magnitude, mdd, med, mee, mnd, mne, mnn, north_shift, time]
pymc.sampling.forward - INFO     Sampling: [duration, east_shift, h_any_P_0_T, h_any_P_0_Z, magnitude, mdd, med, mee, mnd, mne, mnn, north_shift, time]
metropolis   - INFO     Compiling model graph ...
metropolis   - INFO     Initializing proposal distribution ...MultivariateNormal
metropolis   - INFO     Time for proposal covariance init: 0.022191
models       - INFO     Compilation time: 5.248945 

models.base  - INFO     ... Starting SMC ...

backend      - INFO     Removing previous sampling results ... /work/csang/exp/FullMT_nont/geometry/stage_0
backend      - INFO     Found no sampling results under /work/csang/exp/FullMT_nont/geometry/stage_0 
backend      - INFO     Init new trace!
smc          - INFO     Sample initial stage: ...
smc          - INFO     Beta: 0.000000 Stage: 0
sampler      - INFO     Initialising 2000 chain traces ...
sampler      - INFO     Serial time per sample: 0.008421
sampler      - INFO     Chunksize per worker is 20
sampler      - INFO     Data to be memory shared: 
sampler      - INFO     No data to be memshared!
sampler      - INFO     Sampling ...
sampler      - INFO     Worker 11: Finished 4 / 20 chains
...... (redundancy lines are omitted here)
sampler      - INFO     Worker 45: Finished 20 / 20 chains
parallel     - INFO     
 Feierabend! Done with the work!
backend      - INFO     Loading multitrace from /work/csang/exp/FullMT_nont/geometry/stage_0
smc          - INFO     Updating Covariances ...
Traceback (most recent call last):
  File "/home/csang/miniforge3/envs/beat_dev/bin/beat", line 8, in <module>
    sys.exit(main())
             ^^^^^^
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/apps/beat.py", line 2441, in main
    globals()["command_" + command](args)
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/apps/beat.py", line 1068, in command_sample
    sample(step, problem)
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/models/base.py", line 262, in sample
    smc_sample(
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/sampler/smc.py", line 494, in smc_sample
    map_pt = step.get_map_end_points()
             ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/beat/sampler/smc.py", line 288, in get_map_end_points
    return self.bij.rmap(self.array_population[idx, :])
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/csang/miniforge3/envs/beat_dev/lib/python3.11/site-packages/pymc/blocking.py", line 105, in rmap
    raise TypeError("`array` must be a `RaveledVars` type")
TypeError: `array` must be a `RaveledVars` type

In short, I cannot reproduce the last figure of example 1, but I can reproduce other steps of example 1. Attached is my config file for example 1. FullMT_config.zip

Please disable update_covariances and rerun. Non-toeplitz may not be used safely with that option- it is very unstable and you need to know what you are doing.

hvasbath commented 2 months ago

For example 4, I have a question. I can run all the steps of example 4 but the results are not as good as that in the tutorial. I tried several times but I cannot reproduce similar results. I know for nonlinear problems, the results of every execution are not the same, but I think my results are not good. Maybe there is something wrong with my config file. Are there any tricks to improve the results of the finite fault inversion? Attached is my config file and the results for example 4. example4_ffi.zip

Your posteriors for h_laplacian and h_SAR hit the lower prior bounds. So please reduce them both further, maybe -5 and -1 respectively?

sangchengfang commented 2 months ago

Following your suggestions, I ran example 1 and example 2 successfully today. For example 1, I followed your instructions and the results seem alright. For example 2, I changed “PT” sampler to “SMC” sampler and got good results. I have 104 cores on the cluster of my university and my MPI version is 5.0.5, I think my installation is correct since I can reproduce the example 2 successfully.

For the plotting issue, I downloaded the InSAR data from the GitHub repository and got the plots successfully. I am sorry I didn't notice this detail. I thought the InSAR data to be plotted is already in “geodetic_data.pkl” under the project folder.

Thanks for your help!

sangchengfang commented 2 months ago

I adjust the range of hyperparameters h_SAR and h_laplacian to [-5, 5] and [0, 7], respectively. The result seems much better. Attached is the new results.

example4_ffi_adjust_hyper_range.zip

Are those results reasonable?

I also wonder what the black arrows and grey arrows mean in the slip-distribution figure. Is the black arrows the inverted fault slip? I read both the document and the source code but found no description of them.

hvasbath commented 2 months ago

Yes, I think those results wont get much better than that. In the paper I scaled the vectors and ellipses not to be as big as in the figure now, which is why the results appear different. The black arrows are the MAP model and the grey arrows are the mean of the distributions. Sorry yes, the figure caption only describes the black. My bad. image

Please use the --nensemble=300 arguemnt in the scene_fit plot to get the variance reduction value histogram in the top right. Thats more informative and tells you the percantage of data that is well explained.

sangchengfang commented 2 months ago

Thanks for the suggestions! The “ffi” inversion and InSAR fit plots seem alright now.

I have another question about "gnss_fits". I finished my “ffi” inversion and now I want to plot the data fitting to check the performance of the inverted fault slip model, so I run

beat plot 202404022358B_e gnss_fits --mode=ffi

and I got

To enable 'bem' mode packages 'pygmsh' and 'cutde' need to be installed.
beat         - INFO     Loading problem ...
config       - INFO     All hierarchicals ok!
config       - INFO     All hyperparameters ok!
config       - INFO     All priors ok!
models       - INFO     ... Initialising Distribution Optimizer ... 

models       - INFO     Analysing problem ...
models       - INFO     ---------------------

geodetic     - INFO     Number of geodetic datasets: 3 
geodetic     - INFO     Number of geodetic data points: 174 
config       - INFO     not solving for utens, got fixed at 0.0
plotting.geodetic - INFO     Drawing GNSS misfits ...
backend      - INFO     Loading multitrace from ./202404022358B_e/ffi/stage_-1
plotting.geodetic - INFO     Trying to load GNSS data from: ./202404022358B_e
config       - INFO     Loading file NGL_TW_1d_20240604_inv.txt ...
pyrocko.model.gnss - WARNING  Station 0000 does not exist in campaign, do nothing.
inputf       - INFO     Loaded data of 58 GNSS stations
config       - INFO     Successfully loaded GNSS data from file NGL_TW_1d_20240604_inv.txt
geodetic     - INFO     Evaluating config for Ramps corrections for datasets...
geodetic     - INFO     No Ramps correction!
geodetic     - INFO     Evaluating config for Euler Pole corrections for datasets...
geodetic     - INFO     No Euler Pole correction!
geodetic     - INFO     Evaluating config for Strain Rate corrections for datasets...
geodetic     - INFO     No Strain Rate correction!
geodetic     - INFO     Initialized 0 hierarchical parameters.
plotting.geodetic - INFO     FFI gnss fit, using reference source ...
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for NGL_TW_1d_20240604_inv.txt_north ...
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for NGL_TW_1d_20240604_inv.txt_east ...
geodetic     - INFO     Retrieving geodetic data-covariances with structure "import" for NGL_TW_1d_20240604_inv.txt_up ...
geodetic     - INFO     Not updating geodetic velocity model-covariances because number of model variations is too low! < 5
plotting.geodetic - WARNING  More than 40 stations disabling station labels ..
mapproject [INFORMATION]: Processing input table data
mapproject [INFORMATION]: Transform 121.024/122.133/23.564/24.5768 -> 0/452.614/0/450.220152817 [point]
mapproject [INFORMATION]: Reading Data Table from Standard Input stream
mapproject [INFORMATION]: Writing Data Table to Standard Output stream
mapproject [INFORMATION]: Projected 2 points
mapproject [INFORMATION]: Input extreme values:  Xmin: 121.023574767 Xmax: 122.13299521 Ymin: 23.563954162 Ymax 24.5768390897
mapproject [INFORMATION]: Output extreme values: Xmin: -0.174309686616 Xmax: 452.611982794 Ymin: -0.0196969664487 Ymax 450.237522235
mapproject [INFORMATION]: Mapped 2 lon-lat pairs to x-y pairs[point]
pscoast [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
pscoast [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
pscoast [INFORMATION]: Selected ice front line as Antarctica coastline
pscoast [INFORMATION]: GSHHG version 2.3.7
pscoast [INFORMATION]: Derived from World Vector Shoreline, CIA WDB-II, and Atlas of the Cryosphere
pscoast [INFORMATION]: Processed by Paul Wessel and Walter H. F. Smith, 1994-2017
pscoast [INFORMATION]: Painting entire map with ocean color first, then draw land on top later
pscoast [INFORMATION]: Adding Rivers...psvelo [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psvelo [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psvelo [INFORMATION]: Reading Data Table from Standard Input stream
psvelo [INFORMATION]: psvelo: 2-D confidence interval and scaling factor 0.950000 2.447747
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 27. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 40. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 44. Consider adding +n<norm> to -A
psvelo [INFORMATION]: 3 vector heads had length exceeding the vector length and were skipped. Consider the +n<norm> modifier to -A
psvelo [INFORMATION]: Number of records read: 58
psvelo [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psvelo [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psvelo [INFORMATION]: Reading Data Table from Standard Input stream
psvelo [INFORMATION]: psvelo: 2-D confidence interval and scaling factor 0.950000 2.447747
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 27. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 40. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 44. Consider adding +n<norm> to -A
psvelo [INFORMATION]: 3 vector heads had length exceeding the vector length and were skipped. Consider the +n<norm> modifier to -A
psvelo [INFORMATION]: Number of records read: 58
psxy [INFORMATION]: Processing input table data
psxy [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psxy [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psxy [INFORMATION]: Reading Data Table from Stream 7f140361aa00
psxy [INFORMATION]: Plotting segment 0
psxy [INFORMATION]: Processing input table data
psxy [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psxy [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psxy [INFORMATION]: Reading Data Table from Stream 7fd9fe18fa00
psxy [INFORMATION]: Plotting segment 0
psbasemap [INFORMATION]: Constructing the basemap
psbasemap [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psbasemap [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psbasemap [INFORMATION]: Save current gridline information to gmt.history
plotting.geodetic - WARNING  More than 40 stations disabling station labels ..
mapproject [INFORMATION]: Processing input table data
mapproject [INFORMATION]: Transform 121.024/122.133/23.564/24.5768 -> 0/452.614/0/450.220152817 [point]
mapproject [INFORMATION]: Reading Data Table from Standard Input stream
mapproject [INFORMATION]: Writing Data Table to Standard Output stream
mapproject [INFORMATION]: Projected 2 points
mapproject [INFORMATION]: Input extreme values:  Xmin: 121.023574767 Xmax: 122.13299521 Ymin: 23.563954162 Ymax 24.5768390897
mapproject [INFORMATION]: Output extreme values: Xmin: -0.174309686616 Xmax: 452.611982794 Ymin: -0.0196969664487 Ymax 450.237522235
mapproject [INFORMATION]: Mapped 2 lon-lat pairs to x-y pairs[point]
pscoast [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
pscoast [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
pscoast [INFORMATION]: Selected ice front line as Antarctica coastline
pscoast [INFORMATION]: GSHHG version 2.3.7
pscoast [INFORMATION]: Derived from World Vector Shoreline, CIA WDB-II, and Atlas of the Cryosphere
pscoast [INFORMATION]: Processed by Paul Wessel and Walter H. F. Smith, 1994-2017
pscoast [INFORMATION]: Painting entire map with ocean color first, then draw land on top later
pscoast [INFORMATION]: Adding Rivers...psvelo [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psvelo [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psvelo [INFORMATION]: Reading Data Table from Standard Input stream
psvelo [INFORMATION]: psvelo: 2-D confidence interval and scaling factor 0.950000 2.447747
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 13. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 16. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 18. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 20. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 23. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 25. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 27. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 29. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 44. Consider adding +n<norm> to -A
psvelo [INFORMATION]: 9 vector heads had length exceeding the vector length and were skipped. Consider the +n<norm> modifier to -A
psvelo [INFORMATION]: Number of records read: 58
psvelo [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psvelo [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psvelo [INFORMATION]: Reading Data Table from Standard Input stream
psvelo [INFORMATION]: psvelo: 2-D confidence interval and scaling factor 0.950000 2.447747
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 13. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 16. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 18. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 20. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 23. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 25. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 27. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 29. Consider adding +n<norm> to -A
psvelo [INFORMATION]: Vector head length exceeds overall vector length near line 44. Consider adding +n<norm> to -A
psvelo [INFORMATION]: 9 vector heads had length exceeding the vector length and were skipped. Consider the +n<norm> modifier to -A
psvelo [INFORMATION]: Number of records read: 58
psxy [INFORMATION]: Processing input table data
psxy [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psxy [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psxy [INFORMATION]: Reading Data Table from Stream 7f8c6dbc2a00
psxy [INFORMATION]: Plotting segment 0
psxy [INFORMATION]: Processing input table data
psxy [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psxy [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psxy [INFORMATION]: Reading Data Table from Stream 7f5b9676ba00
psxy [INFORMATION]: Plotting segment 0
psbasemap [INFORMATION]: Constructing the basemap
psbasemap [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psbasemap [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psbasemap [INFORMATION]: Save current gridline information to gmt.history
plotting.geodetic - INFO     saving figures to ./202404022358B_e/ffi/figures/gnss_-1_max_0_local
mapproject [INFORMATION]: Processing input table data
mapproject [INFORMATION]: Transform 121.024/122.133/23.564/24.5768 -> 0/452.614/0/450.220152817 [point]
mapproject [INFORMATION]: Reading Data Table from Standard Input stream
mapproject [INFORMATION]: Writing Data Table to Standard Output stream
mapproject [INFORMATION]: Projected 2 points
mapproject [INFORMATION]: Input extreme values:  Xmin: 121.023574767 Xmax: 122.13299521 Ymin: 23.563954162 Ymax 24.5768390897
mapproject [INFORMATION]: Output extreme values: Xmin: -0.174309686616 Xmax: 452.611982794 Ymin: -0.0196969664487 Ymax 450.237522235
mapproject [INFORMATION]: Mapped 2 lon-lat pairs to x-y pairs[point]
psxy [INFORMATION]: Processing input table data
psxy [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psxy [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psxy [INFORMATION]: Reading Data Table from Stream 7f770740ca00
psxy [WARNING]: File <stdin> is empty!
mapproject [INFORMATION]: Processing input table data
mapproject [INFORMATION]: Transform 121.024/122.133/23.564/24.5768 -> 0/452.614/0/450.220152817 [point]
mapproject [INFORMATION]: Reading Data Table from Standard Input stream
mapproject [INFORMATION]: Writing Data Table to Standard Output stream
mapproject [INFORMATION]: Projected 2 points
mapproject [INFORMATION]: Input extreme values:  Xmin: 121.023574767 Xmax: 122.13299521 Ymin: 23.563954162 Ymax 24.5768390897
mapproject [INFORMATION]: Output extreme values: Xmin: -0.174309686616 Xmax: 452.611982794 Ymin: -0.0196969664487 Ymax 450.237522235
mapproject [INFORMATION]: Mapped 2 lon-lat pairs to x-y pairs[point]
psxy [INFORMATION]: Processing input table data
psxy [INFORMATION]: gmt_map_setup perimeter search region: 120.9196564616369/122.233/23.46399765223102/24.67780919276663.
psxy [INFORMATION]: Map scale is 7.06293 km per cm or 1:706293.
psxy [INFORMATION]: Reading Data Table from Stream 7fd4e0072a00
psxy [WARNING]: File <stdin> is empty!

However, the figures are not correct. As shown in the figures in the attached zip file below, there are only 10 stations in the horizontal and vertical plots, but I have 58 stations included in the “ffi” inversion. The GPS data file is an ASCII file and its format follows the standard format shown in the document (https://pyrocko.org/beat/docs/current/getting_started/import_data.html#geodetic-data). How can I fix that?

gnss_plot.zip

I also want to export the surface displacement due to the inverted fault slip model (i.e., results of forward modeling) and discretized fault model, so I can plot the data fitting with my code and calculate the Coulomb stress change due to fault slip. Are there any options to do this or do I need to write some code? Could you give me some hints If I need to write some code to export the synthetic surface displacements and fault model cause I do not know where to start? Thanks a lot.

hvasbath commented 2 months ago

To ensure all stations are included in the plotting area please run:

beat plot 202404022358B_e gnss_fits --mode=ffi --plot_projection=latlon

The default is a zoom in on the fault area.

To export the MAP model and most important parameters of interest under $projectdir/$mode/results please run:

beat export 202404022358B_e --mode=ffi --stage_number=-1

Again the hint to please use: beat --help

once in a while to try to find out some of that yourself. Each option again supports again --help, like: beat plot --help

Regards, Hannes

hvasbath commented 1 month ago

Released version 2.0.3 with fix.