Open
Description
Problem description
I am trying to implement geostatistical regularization in my TravelTime inversion, but I encounter a segmentation fault at the GeostatisticConstraintsMatrix creation step.
Interestingly, the code runs fine on my personal machine but not on the remote cluster I want to use. (I use the same code and I installed pygimli in the same way on both).
Moreover, everything works correctly on both systems if I call the invert() method without the correlationLengths argument.
I checked issues #659 and #790 and tried the suggested solutions, but they did not resolve the problem. (The segmentation fault occurs on the remote machine as soon as I call the pg.matrix.GeostatisticConstraintsMatrix() method.)
Your environment
Local machine
Date: Fri Jan 31 13:51:24 2025 CET
OS : Linux (Ubuntu 22.04)
CPU(s) : 14
Machine : x86_64
Architecture : 64bit
Environment : Python
Python 3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:36:13) [GCC
12.3.0]
pygimli : 1.5.3
pgcore : 1.5.0
numpy : 1.26.4
matplotlib : 3.9.1
scipy : 1.14.1
tqdm : 4.67.1
IPython : 8.31.0
pyvista : 0.44.2
Remote machine (the problematic one)
Date: Fri Jan 31 13:50:20 2025 CET
OS : Linux (Debian GNU/Linux 11)
CPU(s) : 32
Machine : x86_64
Architecture : 64bit
Environment : Python
Python 3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:36:13) [GCC
12.3.0]
pygimli : 1.5.3
pgcore : 1.5.0
numpy : 1.26.4
matplotlib : 3.9.1
scipy : 1.15.1
tqdm : 4.67.1
IPython : 8.31.0
pyvista : 0.44.2
Steps to reproduce
Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:
import matplotlib.pyplot as plt
import pygimli as pg
import pygimli.physics.traveltime as tt
import pygimli.meshtools as mt
from pygimli.viewer import pv
import pyvista
import numpy as np
import ipdb
import glob
import os
import time
import random
###############################################################################
# # Set data parameters
data_file = "../picking_corrected_002ms/pygimly_data_reciprocity_error_based.sgt"
sensors_to_remove = [65,66,67,68,69,75,76,77,78,79,81,82,83,84]
max_validity = 10
max_error = 0.5
###############################################################################
# # Set Mesh parameters
paraDX = 0.5
paraDepth = 20
paraBoundary = None
paraMaxCellSize = 400
boundary = None
boundaryMaxCellSize = 0
surfaceMeshQuality = 20
surfaceMeshArea = 20
addTopo = None
isClosed = False
###############################################################################
# # Set Inversion parameters
secNodes = 0
vTop = 300
vBottom = 3000
zWeight = 0.2
correlationLengths = [5,5,5]
resultFolder = "../TT_tomography/run_test/
###############################################################################
# # Load data
data = pg.load(data_file, verbose=True)
###############################################################################
# We remove sensors that we don't want, which automatically removes the measurements associated with them.
data.removeSensorIdx(np.array(sensors_to_remove)-1) # -1 to move from sensor name to sensor index
###############################################################################
# We center the sensors positions
sensors = np.array(data.sensors())
sensors[:,0]=sensors[:,0] - np.median(sensors[:,0])
sensors[:,1]=sensors[:,1] - np.median(sensors[:,1])
sensors[:,2]=(sensors[:,2] - np.median(sensors[:,2]))
data.setSensors(sensors)
###############################################################################
# We remove the lowest-quality measures and the largest errors.
data.markInvalid(data("validity") >= max_validity)
if data.haveData("relative_err"):
data.markInvalid(data("relative_err") >= max_error)
data.removeInvalid()
###############################################################################
# We initialize the refraction manager.
mgr = tt.TravelTimeManager(data)
###############################################################################
# # PLC mesh with topography
plc = mt.createParaMeshPLC3D(data.sensors(), paraDX=paraDX, paraDepth=paraDepth, paraBoundary=paraBoundary, paraMaxCellSize=paraMaxCellSize, boundary=boundary, boundaryMaxCellSize=boundaryMaxCellSize, surfaceMeshQuality=surfaceMeshQuality, surfaceMeshArea=surfaceMeshArea, addTopo=addTopo, isClosed=isClosed)
mesh = mt.createMesh(plc)
mesh = mesh.createSubMesh(mesh.cells(mesh.cellMarkers()==2))
# ###############################################################################
# # Finally, we call the `invert` method and plot the result.
mgr.invert(mesh = mesh, secNodes=secNodes, vTop=vTop, vBottom=vBottom, verbose=2, zWeight=zWeight, correlationLengths = correlationLengths)
# ###############################################################################
# # Save the results
models_history = mgr.inv.modelHistory
mgr.saveResult(folder=resultFolder)
np.save(resultFolder+"/TravelTimeManager/models_history.npy", models_history, allow_pickle=True)Expected behavior
To work properly on both systems
Actual behavior
The output on my personal machine
Opening /tmp/tmp0gsyjkvi.poly.
Segments are connected properly.
31/01/25 - 12:02:34 - pyGIMLi - INFO - Found 1 regions.
31/01/25 - 12:02:34 - pyGIMLi - INFO - Creating forward mesh from region infos.
31/01/25 - 12:02:34 - pyGIMLi - INFO - Creating refined mesh (secnodes: 2) to solve forward task.
min/max(dweight) = 164.026/3728.56
31/01/25 - 12:02:35 - pyGIMLi - INFO - Create gradient starting model. 300.0: 5000.0
31/01/25 - 12:02:35 - pyGIMLi - INFO - Created startmodel from forward operator:3360, min/max=0.000200/0.003333
31/01/25 - 12:02:35 - pyGIMLi - INFO - Starting inversion.
fop: <pygimli.physics.traveltime.modelling.TravelTimeDijkstraModelling object at 0x73da8d945d00>
Data transformation: Identity transform
Model transformation (cumulative):
0 Logarithmic LU transform, lower bound 0.0, upper bound 0.0
min/max (data): 0.0089/0.1
min/max (error): 3%/9%
min/max (start model): 2.0e-04/0.0033
--------------------------------------------------------------------------------
use model trans from RegionManager
min/max(dweight) = 164.026/3728.56
Building constraints matrix
31/01/25 - 12:03:03 - pyGIMLi - INFO - Creating GeostatisticConstraintsMatrix for region 2 with: I=[5, 5, 5], dip=0, strike=0
31/01/25 - 12:03:03 - pyGIMLi - INFO - Cache /home/herbelom/anaconda3/envs/pg/lib/python3.11/site-packages/pygimli/math/matrix.py:createCm05 restored (27.9s x 3): /home/herbelom/.cache/pygimli/4315247502913639306
constraint matrix of size(nBounds x nModel) 3360 x 3360
31/01/25 - 12:03:03 - Core - ERROR - no cWeights defined. You should create constraints matrix first.
check Jacobian: wrong dimensions: (0x0) should be (3496x3360) force: 1
jacobian size invalid, forced recalc: 1
Calculating Jacobian matrix (forced=1)...... 27.6059 s
min data = 0.00894 max data = 0.10161 (3496)
min error = 0.03 max error = 0.09 (3496)
min response = 0.0024178 max response = 0.0377571 (3496)
calc without reference model
0: rms/rrms(data, response) = 0.0159771/47.8529%
0: chi^2(data, response, error, log) = 188.18
0: Phi = py + 6529.49 * 20 = 788466
inv.iter 0 ... chi² = 188.18
--------------------------------------------------------------------------------
solve CGLSCDWWtrans with lambda = 20
1: LS newModel: min = 4.08786e-05; max = 0.0503592
1: LS newResponse: min = 0.00150161; max = 0.10306
1: rms/rrms(data, LS newResponse) = 0.0123024/39.7984%
1: chi^2(data, LS newResponse, error, log) = 130.411
1: Phi = 455918+1059.4*20=477106
Performing line search with tau = 0.83
1: Model: min = 6.12102e-05; max = 0.023899
1: Response: min = 0.00216225; max = 0.0858382
1: rms/rrms(data, Response) = 0.0109195/35.7095%
1: chi^2(data, Response, error, log) = 109.904
1: Phi = 384224+1097.68*20=406177
inv.iter 1 ... chi² = 109.90 (dPhi = 48.49%) lam: 20.0
The segmentation fault (on the remote cluster)
Opening /tmp/tmp3nylsgmc.poly.
Segments are connected properly.
31/01/25 - 13:51:29 - pyGIMLi - INFO - Found 1 regions.
31/01/25 - 13:51:29 - pyGIMLi - INFO - Creating forward mesh from region infos.
31/01/25 - 13:51:29 - pyGIMLi - INFO - Creating refined mesh (secnodes: 2) to solve forward task.
min/max(dweight) = 164.026/3728.56
31/01/25 - 13:51:30 - pyGIMLi - INFO - Create gradient starting model. 300.0: 5000.0
31/01/25 - 13:51:30 - pyGIMLi - INFO - Created startmodel from forward operator:3360, min/max=0.000200/0.003333
31/01/25 - 13:51:30 - pyGIMLi - INFO - Starting inversion.
fop: <pygimli.physics.traveltime.modelling.TravelTimeDijkstraModelling object at 0x7f171d75fb50>
Data transformation: Identity transform
Model transformation (cumulative):
0 Logarithmic LU transform, lower bound 0.0, upper bound 0.0
min/max (data): 0.0089/0.1
min/max (error): 3%/9%
min/max (start model): 2.0e-04/0.0033
-------------------------------------
use model trans from RegionManager
min/max(dweight) = 164.026/3728.56
Building constraints matrix
31/01/25 - 13:55:44 - pyGIMLi - INFO - Creating GeostatisticConstraintsMatrix for region 2 with: I=[5, 5, 5], dip=0, strike=0
Segmentation fault
If possible, please add one or more labels to your issue, e.g. if you expect that your issue is rather a question than a problem with the code, please add the label "question".