OSCR

Three-dimensional inversion of gravity data using implicit neural representations and scientific machine learning.

Code ↔ Paper

15 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 15 matches
  1. [1] § Methods › Implicit neural representation and spatial encodings ↔ 001-EncodingComparisons_blocky.py, lines 120–245 · score 0.76 · learnable feature vectors, trilinear interpolation, hash encoding, concatenated, corner, embeddings
  2. [2] § Methods › Implicit neural representation and spatial encodings ↔ 001-EncodingComparisons_blocky.py, lines 118–243 · score 0.76 · learnable feature vectors, trilinear interpolation, hash encoding, concatenated, corner, embeddings
  3. [3] § Results › Spectral bias and spatial encoding comparison ↔ v01/05-TestingNoiseSensitivity.py, lines 1–54 · score 0.74 · rectangular prism kernel, stationary Gaussian, Gaussian noise, gravity inversion, whitened, surface
  4. [4] § Methods › Forward modelling of gravity data ↔ v01/04-BlockModel.py, lines 109–130 · score 0.69 · vertical component, rectangular prism, sensitivity matrix, cells, gravity, model
  5. [5] § Results › Network capacity and grid refinement ↔ 002a-BlockNetworkSizeComparison.py, lines 268–399 · score 0.66 · model error, training loss, voxel model, ratio, INR parameter, largest
  6. [6] § Methods › Forward modelling of gravity data ↔ v01/05-TestingNoiseSensitivity.py, lines 76–90 · score 0.64 · vertical component, rectangular prism, arctan, sensitivity, gravity
  7. [7] § Results › Implicit versus explicit regularisation ↔ v01/04-BlockModel.py, lines 232–281 · score 0.64 · smallness term, linear system, depth weighted, gradient, smoothness, noise
  8. [8] § Results › Spectral bias and spatial encoding comparison ↔ 001-EncodingComparisons_blocky.py, lines 120–245 · score 0.63 · finest resolution, noise floor, hash encoding, encoding comparison, geometric, blocky
  9. [9] § Results › Spectral bias and spatial encoding comparison ↔ 001-EncodingComparisons_blocky.py, lines 118–243 · score 0.63 · finest resolution, noise floor, hash encoding, encoding comparison, geometric, blocky
  10. [10] § Results › Implicit versus explicit regularisation ↔ 003-ImplictvsExplicitRegularisation.py, lines 1–21 · score 0.60 · explicit smoothness, hidden layers, frequency bands, positional encoding, implicit
  11. [11] § Results › Spectral bias and spatial encoding comparison ↔ 001-EncodingComparisons_blocky.py, lines 73–117 · score 0.60 · overly smooth, sinusoidal positional encoding, lowest, blocky, MLP, bias
  12. [12] § Results › Spectral bias and spatial encoding comparison ↔ 001-EncodingComparisons_blocky.py, lines 71–115 · score 0.60 · overly smooth, sinusoidal positional encoding, lowest, blocky, MLP, bias
  13. [13] § Results › Network capacity and grid refinement ↔ 002b-BlockGridSizeComparison.py, lines 1–20 · score 0.53 · voxel model parameter, INR architectures, INR parameter, gravity inversion, block model, discretisation
  14. [14] § Results › Network capacity and grid refinement ↔ 004-ModelEnsambles.py, lines 1–23 · score 0.51 · kept fixed, INR architecture, block model, medium, positional encoding, Network
  15. [15] § Results › Network capacity and grid refinement ↔ 002b-BlockGridSizeComparison.py, lines 1–20 · score 0.50 · voxel model parameter, INR parameter, refined, architecture, discretisation, inversion

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 759 lines · 30 KB · MIT · 3 matches

  1. import os
  2. import sys
  3. import time
  4. import random
  5. import numpy as np
  6. import torch
  7. import torch.nn as nn
  8. import matplotlib
  9. import matplotlib.pyplot as plt
  13. # --- Random seed -------------------------------------------------------
  14. SEED = 42 #10,42, 40, 41
  15. DATA_SEED = SEED
  17. # --- Grid / domain -----------------------------------------------------
  18. DX = 50.0 # cell size in x (m)
  19. DY = 50.0 # cell size in y (m)
  20. DZ = 50.0 # cell size in z (m)
  21. X_MAX = 1000.0 # domain extent in x (m)
  22. Y_MAX = 1000.0 # domain extent in y (m)
  23. Z_MAX = 500.0 # domain extent in z (m)
  25. # --- Block model -------------------------------------------------------
  26. RHO_BG = 0.0 # background density contrast (kg/m³)
  27. RHO_BLK = 400.0 # block density contrast (kg/m³)
  29. # --- Noise level -------------------------------------------------------
  30. NOISE_LEVEL = 0.01 # fraction of gz_true std
  32. # --- Training / optimisation -------------------------------------------
  33. GAMMA = 1.0 # data-term weight
  34. EPOCHS = 500 # number of training epochs
  35. LR = 1e-2 # Adam learning rate
  37. # --- Early stopping ----------------------------------------------------
  38. # Stop near the expected noise floor to reduce overfitting to 1% noise.
  39. USE_EARLY_STOPPING = True
  40. EARLY_STOP_MIN_EPOCHS = 100
  41. EARLY_STOP_PATIENCE = 25
  42. EARLY_STOP_TARGET = 1.0
  43. EARLY_STOP_TOL = 0.05
  44. EARLY_STOP_OVERFIT_PATIENCE = 5
  46. # --- INR network -------------------------------------------------------
  47. HIDDEN = 256 # hidden-layer width
  48. DEPTH = 4 # number of hidden layers
  49. RHO_ABS_MAX = 600.0 # tanh output scaling (kg/m³)
  51. # --- Encoding strategy -------------------------------------------------
  52. # Compare only the two encoded models in this script.
  53. ENCODING_ORDER = ('positional', 'hash')
  55. ENCODING_CONFIGS = {
  56. 'positional': dict(num_freqs=2),
  57. 'hash': dict(n_levels=2, n_features_per_level=2,
  58. log2_hashmap_size=17, base_resolution=4,
  59. finest_resolution=128),
  60. }
  62. # --- Plotting -----------------------------------------------------------
  63. CMAP = 'Spectral_r' # colormap for all plots
  64. INV_VMAX = 250 # fixed colorbar max for inverted model
  65. FIG_DPI = 300
  66. TITLE_FONTSIZE = 14
  67. LABEL_FONTSIZE = 12
  68. TICK_FONTSIZE = 11
  69. LEGEND_FONTSIZE = 10
  73. class PositionalEncoding(nn.Module):
  74. """Sinusoidal positional encoding (NeRF-style Fourier features).
  76. Intuition
  77. ---------
  78. Plain MLPs are biased toward learning low-frequency functions — they
  79. tend to produce overly smooth outputs and struggle to represent sharp
  80. edges or fine detail ("spectral bias"). Positional encoding fights
  81. this by lifting each input coordinate into a set of sine and cosine
  82. waves at *exponentially increasing* frequencies:
  84. γ(p) = [sin(2⁰πp), cos(2⁰πp), sin(2¹πp), cos(2¹πp), …]
  86. Think of it as giving the network a set of "rulers" at different
  87. scales: the lowest frequency captures the overall trend, while
  88. higher frequencies let the network resolve increasingly fine spatial
  89. variation. The resulting feature vector is fixed (no learnable
  90. parameters in the encoding itself).
  92. Trade-offs
  93. ----------
  94. • Simple and deterministic — no extra learnable parameters.
  95. • The frequency ladder is fixed (powers of 2), so the spectrum can
  96. have gaps and may not cover all scales equally well.
  97. • Too few frequencies → overly smooth model; too many → potential
  98. noise fitting and slower convergence.
  99. • Works best with a sufficiently deep/wide MLP.
  101. Args:
  102. num_freqs: Number of frequency octaves (L). Output grows
  103. as input_dim × (1 + 2L) with include_input.
  104. include_input: Whether to prepend the raw (x, y, z) values.
  105. input_dim: Spatial dimensionality (default 3).
  106. """
  107. def __init__(self, num_freqs=8, include_input=True, input_dim=3):
  108. super().__init__()
  109. self.include_input = include_input
  110. self.register_buffer('freqs', 2.0 ** torch.arange(0, num_freqs))
  111. self.out_dim = input_dim * (1 + 2 * num_freqs) if include_input else input_dim * 2 * num_freqs
  113. def forward(self, x):
  114. parts = [x] if self.include_input else []
  115. for f in self.freqs:
  116. parts += [torch.sin(f * x), torch.cos(f * x)]
  117. return torch.cat(parts, dim=-1)
  120. class HashEncoding(nn.Module):
  121. """Multi-resolution hash encoding (Müller et al., 2022 / Instant-NGP).
  123. Intuition
  124. ---------
  125. Imagine overlaying your 3-D domain with a stack of voxel grids, from
  126. very coarse (e.g. 4³) to very fine (e.g. 128³). At each resolution
  127. level, every voxel stores a small learnable feature vector. To
  128. encode a point, you look up its 8 surrounding voxel corners at every
  129. level, trilinearly interpolate, and concatenate across levels.
  131. The "hash" trick makes this memory-efficient: instead of allocating a
  132. full 3-D grid (which grows as O(R³)), each level maps voxel corners
  133. to a fixed-size hash table. Collisions (two corners sharing one
  134. slot) are resolved implicitly by the gradient-based optimisation —
  135. the network learns to disentangle them.
  137. The result is an encoding that is:
  138. • **Adaptive**: features are *learned*, so detail concentrates where
  139. the data demands it (unlike fixed Fourier features).
  140. • **Multi-scale**: coarse levels capture large-scale trends; fine
  141. levels capture sharp boundaries and small anomalies.
  142. • **Memory-bounded**: hash table size is constant regardless of
  143. how fine the resolution is — crucial for large 3-D domains.
  145. Trade-offs
  146. ----------
  147. • Most flexible encoding — can represent both smooth and sharp models.
  148. • Hash collisions at fine levels can introduce small artifacts or a
  149. noise floor if the hash table is too small (increase
  150. log2_hashmap_size to mitigate).
  151. • Many hyperparameters to tune (n_levels, base/finest resolution,
  152. hash table size).
  153. • Learnable parameters mean more total parameters to optimise, and
  154. they may overfit noisy data without appropriate regularisation.
  156. Args:
  157. n_levels: Number of resolution levels. More levels
  158. give smoother multi-scale interpolation.
  159. n_features_per_level: Feature dimensions stored per hash entry.
  160. Typically 2; larger values add capacity.
  161. log2_hashmap_size: Log₂ of hash-table size per level.
  162. 2^19 ≈ 500 k entries is a common default.
  163. base_resolution: Coarsest voxel-grid resolution (e.g. 4).
  164. finest_resolution: Finest voxel-grid resolution (e.g. 128).
  165. Intermediate levels are spaced
  166. geometrically between base and finest.
  167. input_dim: Spatial dimensionality (default 3).
  168. """
  169. def __init__(self, n_levels=16, n_features_per_level=2,
  170. log2_hashmap_size=19, base_resolution=16,
  171. finest_resolution=512, input_dim=3):
  172. super().__init__()
  173. self.n_levels = n_levels
  174. self.n_features_per_level = n_features_per_level
  175. self.input_dim = input_dim
  176. self.out_dim = n_levels * n_features_per_level
  177. self.hashmap_size = 2 ** log2_hashmap_size
  179. if n_levels > 1:
  180. self.growth_factor = np.exp(
  181. (np.log(finest_resolution) - np.log(base_resolution))
  182. / (n_levels - 1))
  183. else:
  184. self.growth_factor = 1.0
  185. self.base_resolution = base_resolution
  187. # Learnable hash tables (one per level)
  188. self.hash_tables = nn.ModuleList([
  189. nn.Embedding(self.hashmap_size, n_features_per_level)
  190. for _ in range(n_levels)
  191. ])
  192. for tbl in self.hash_tables:
  193. nn.init.uniform_(tbl.weight, -1e-4, 1e-4)
  195. # Large primes for the spatial hash
  196. self.register_buffer(
  197. 'primes', torch.tensor([1, 2654435761, 805459861], dtype=torch.long))
  199. def _hash(self, coords_int):
  200. """Spatial hash: integer grid coords -> hash-table index."""
  201. result = torch.zeros(coords_int.shape[:-1],
  202. dtype=torch.long, device=coords_int.device)
  203. for d in range(self.input_dim):
  204. result ^= coords_int[..., d] * self.primes[d]
  205. return result % self.hashmap_size
  207. def forward(self, x):
  208. # Normalise to [0, 1] using per-batch bounds
  209. x_min = x.min(dim=0, keepdim=True).values
  210. x_max = x.max(dim=0, keepdim=True).values
  211. x_scaled = (x - x_min) / (x_max - x_min + 1e-8)
  213. outputs = []
  214. for level in range(self.n_levels):
  215. resolution = int(self.base_resolution * (self.growth_factor ** level))
  216. x_grid = x_scaled * resolution # (N, 3)
  217. x_floor = torch.floor(x_grid).long() # voxel origin
  218. x_frac = x_grid - x_floor.float() # interpolation weight
  220. # Eight voxel corners
  221. corners = []
  222. for dz in (0, 1):
  223. for dy in (0, 1):
  224. for dx in (0, 1):
  225. corners.append(
  226. x_floor + torch.tensor([dx, dy, dz], device=x.device))
  227. corners = torch.stack(corners, dim=1) # (N, 8, 3)
  229. indices = self._hash(corners) # (N, 8)
  230. features = self.hash_tables[level](indices) # (N, 8, F)
  232. # Trilinear interpolation weights
  233. wx, wy, wz = (x_frac[:, 0:1],
  234. x_frac[:, 1:2],
  235. x_frac[:, 2:3])
  236. weights = torch.stack([
  237. (1-wx)*(1-wy)*(1-wz), wx*(1-wy)*(1-wz),
  238. (1-wx)* wy *(1-wz), wx* wy *(1-wz),
  239. (1-wx)*(1-wy)* wz , wx*(1-wy)* wz ,
  240. (1-wx)* wy * wz , wx* wy * wz ,
  241. ], dim=1) # (N, 8, 1)
  243. outputs.append((weights * features).sum(dim=1)) # (N, F)
  245. return torch.cat(outputs, dim=-1) # (N, n_levels*F)
  248. def create_encoding(encoding_type, **kwargs):
  249. """Factory: build an encoding module by name.
  251. Supported types
  252. ---------------
  253. positional – Sinusoidal Fourier features (Mildenhall et al., 2020
  254. / NeRF). Fixed frequencies on a power-of-2 ladder.
  255. Simple, no learnable params; good general-purpose
  256. baseline. Tune `num_freqs` for resolution.
  258. hash – Multi-resolution hash tables (Müller et al., 2022 /
  259. Instant-NGP). Learnable feature grids at multiple
  260. resolutions compressed via spatial hashing. Most
  261. flexible; best for sharp boundaries. More params &
  262. hyperparameters.
  263. """
  264. if encoding_type == 'positional':
  265. enc = PositionalEncoding(
  266. num_freqs=kwargs.get('num_freqs', 8),
  267. include_input=kwargs.get('include_input', True))
  268. return enc
  269. if encoding_type == 'hash':
  270. return HashEncoding(
  271. n_levels=kwargs.get('n_levels', 16),
  272. n_features_per_level=kwargs.get('n_features_per_level', 2),
  273. log2_hashmap_size=kwargs.get('log2_hashmap_size', 19),
  274. base_resolution=kwargs.get('base_resolution', 16),
  275. finest_resolution=kwargs.get('finest_resolution', 512))
  276. raise ValueError(f"Unknown encoding type: {encoding_type}. Expected 'positional' or 'hash'.")
  279. class DensityContrastINR(nn.Module):
  280. """INR density-contrast model with pluggable spatial encoding."""
  281. def __init__(self, encoding_type='positional', hidden=256, depth=5,
  282. rho_abs_max=600.0, **encoding_kwargs):
  283. super().__init__()
  284. self.pe = create_encoding(encoding_type, **encoding_kwargs)
  285. in_dim = self.pe.out_dim
  286. layers = []
  287. h = hidden
  288. layers += [nn.Linear(in_dim, h), nn.LeakyReLU(0.01)]
  289. for _ in range(depth - 1):
  290. layers += [nn.Linear(h, h), nn.LeakyReLU(0.01)]
  291. layers += [nn.Linear(h, 1)]
  292. self.net = nn.Sequential(*layers)
  293. self.rho_abs_max = float(rho_abs_max)
  295. def forward(self, x):
  296. z = self.pe(x)
  297. out = self.net(z)
  298. return self.rho_abs_max * torch.tanh(out)
  300. # ──────────────────────────────────────────────────────────────────────
  301. # UTILITY FUNCTIONS
  302. # ──────────────────────────────────────────────────────────────────────
  304. def set_seed(seed: int = 42):
  305. random.seed(seed)
  306. np.random.seed(seed)
  307. torch.manual_seed(seed)
  308. if torch.cuda.is_available():
  309. torch.cuda.manual_seed(seed)
  310. torch.cuda.manual_seed_all(seed)
  311. torch.backends.cudnn.deterministic = True
  312. torch.backends.cudnn.benchmark = False
  313. os.environ['PYTHONHASHSEED'] = str(seed)
  314. print(f"Seed = {seed}")
  317. def capture_rng_state():
  318. state = {
  319. 'python': random.getstate(),
  320. 'numpy': np.random.get_state(),
  321. 'torch': torch.get_rng_state(),
  322. }
  323. if torch.cuda.is_available():
  324. state['cuda'] = torch.cuda.get_rng_state_all()
  325. return state
  328. def restore_rng_state(state):
  329. random.setstate(state['python'])
  330. np.random.set_state(state['numpy'])
  331. torch.set_rng_state(state['torch'])
  332. if torch.cuda.is_available() and 'cuda' in state:
  333. torch.cuda.set_rng_state_all(state['cuda'])
  335. def A_integral_torch(x, y, z):
  336. eps = 1e-20
  337. r = torch.sqrt(x**2 + y**2 + z**2).clamp_min(eps)
  338. return -(x * torch.log(torch.abs(y + r) + eps) +
  339. y * torch.log(torch.abs(x + r) + eps) -
  340. z * torch.atan2(x * y, z * r + eps))
  342. @torch.inference_mode()
  343. def construct_sensitivity_matrix_G_torch(cell_grid, data_points, d1, d2, device):
  344. Gamma = 6.67430e-11
  345. cx = cell_grid[:, 0].unsqueeze(0)
  346. cy = cell_grid[:, 1].unsqueeze(0)
  347. cz = cell_grid[:, 2].unsqueeze(0)
  348. czh = cell_grid[:, 3].unsqueeze(0)
  349. ox = data_points[:, 0].unsqueeze(1)
  350. oy = data_points[:, 1].unsqueeze(1)
  351. oz = data_points[:, 2].unsqueeze(1)
  352. x2, x1 = (cx + d1/2) - ox, (cx - d1/2) - ox
  353. y2, y1 = (cy + d2/2) - oy, (cy - d2/2) - oy
  354. z2, z1 = (cz + czh) - oz, (cz - czh) - oz
  355. A = (A_integral_torch(x2, y2, z2) - A_integral_torch(x2, y2, z1) -
  356. A_integral_torch(x2, y1, z2) + A_integral_torch(x2, y1, z1) -
  357. A_integral_torch(x1, y2, z2) + A_integral_torch(x1, y2, z1) +
  358. A_integral_torch(x1, y1, z2) - A_integral_torch(x1, y1, z1))
  359. return (Gamma * A).to(device)
  361. def generate_grf_torch(nx, ny, nz, dx, dy, dz, lam, nu, sigma, device):
  362. kx = torch.fft.fftfreq(nx, d=dx, device=device) * 2 * torch.pi
  363. ky = torch.fft.fftfreq(ny, d=dy, device=device) * 2 * torch.pi
  364. kz = torch.fft.fftfreq(nz, d=dz, device=device) * 2 * torch.pi
  365. Kx, Ky, Kz = torch.meshgrid(kx, ky, kz, indexing='ij')
  366. k2 = Kx**2 + Ky**2 + Kz**2
  367. P = (k2 + (1/lam**2))**(-nu - 1.5)
  368. P[0, 0, 0] = 0
  369. noise = torch.randn(nx, ny, nz, dtype=torch.complex64, device=device)
  370. f = noise * torch.sqrt(P)
  371. m = torch.real(torch.fft.ifftn(f))
  372. m = (m - m.mean()) / (m.std() + 1e-9)
  373. return sigma * m
  375. def train_inr(model, opt, coords_norm, G, gz_obs, Wd, Nx, Ny, Nz, dx, dy, dz, cfg):
  376. history = {"total": [], "gravity": []}
  377. use_early_stopping = cfg.get('use_early_stopping', False)
  378. min_epochs = cfg.get('early_stop_min_epochs', 0)
  379. patience = cfg.get('early_stop_patience', 0)
  380. target = cfg.get('early_stop_target', 1.0)
  381. tol = cfg.get('early_stop_tol', 0.0)
  382. overfit_patience = cfg.get('early_stop_overfit_patience', 0)
  383. best_gap = float('inf')
  384. best_epoch = -1
  385. best_weighted_mse = None
  386. best_state = None
  387. in_band_count = 0
  388. overfit_count = 0
  390. for ep in range(cfg['epochs']):
  391. opt.zero_grad()
  392. m_pred = model(coords_norm).view(-1)
  393. gz_pred = torch.matmul(G, m_pred.unsqueeze(1)).squeeze(1)
  394. residual = gz_pred - gz_obs
  395. data_term = cfg['gamma'] * torch.mean((Wd * residual) ** 2)
  396. loss = data_term
  397. loss.backward()
  398. opt.step()
  399. history['gravity'].append(float(data_term.item()))
  400. history['total'].append(float(loss.item()))
  402. if use_early_stopping:
  403. weighted_mse = float(data_term.item() / cfg['gamma'])
  404. gap = abs(weighted_mse - target)
  405. if gap < best_gap:
  406. best_gap = gap
  407. best_epoch = ep
  408. best_weighted_mse = weighted_mse
  409. best_state = {
  410. name: value.detach().cpu().clone()
  411. for name, value in model.state_dict().items()
  412. }
  414. if ep + 1 >= min_epochs:
  415. if target - tol <= weighted_mse <= target + tol:
  416. in_band_count += 1
  417. else:
  418. in_band_count = 0
  420. if weighted_mse < target - tol:
  421. overfit_count += 1
  422. else:
  423. overfit_count = 0
  425. if in_band_count >= patience:
  426. print(
  427. f"Early stopping at epoch {ep:4d} | "
  428. f"reason = near noise floor | "
  429. f"loss = {weighted_mse:.3f} | "
  430. f"best epoch = {best_epoch}"
  431. )
  432. break
  434. if overfit_count >= overfit_patience:
  435. print(
  436. f"Early stopping at epoch {ep:4d} | "
  437. f"reason = below noise floor | "
  438. f"loss = {weighted_mse:.3f} | "
  439. f"best epoch = {best_epoch}"
  440. )
  441. break
  443. if ep % 50 == 0 or ep == cfg['epochs'] - 1:
  444. print(f"Epoch {ep:4d} | loss {history['gravity'][-1]:.3e}")
  446. if use_early_stopping and best_state is not None:
  447. model.load_state_dict(best_state)
  448. history['best_epoch'] = best_epoch
  449. history['best_weighted_mse'] = best_weighted_mse
  450. else:
  451. history['best_epoch'] = len(history['gravity']) - 1
  452. history['best_weighted_mse'] = history['gravity'][-1] / cfg['gamma']
  454. return history
  456. def make_block_model(Nx, Ny, Nz, dx, dy, dz, rho_bg=0.0, rho_blk=400.0):
  457. m = torch.full((Nx, Ny, Nz), rho_bg)
  458. for i in range(7):
  459. z_idx = 1 + i
  460. y_start, y_end = 11 - i, 16 - i
  461. x_start, x_end = 7, 13
  462. if 0 <= z_idx < Nz:
  463. ys, ye = max(0, y_start), min(Ny, y_end)
  464. xs, xe = max(0, x_start), min(Nx, x_end)
  465. m[xs:xe, ys:ye, z_idx] = rho_blk
  466. return m.view(-1), m
  468. def get_block_boundaries(Nx, Ny, Nz):
  469. boundaries = []
  470. for i in range(7):
  471. z_idx = 1 + i
  472. y_start, y_end = 11 - i, 16 - i
  473. x_start, x_end = 7, 13
  474. if 0 <= z_idx < Nz:
  475. ys, ye = max(0, y_start), min(Ny, y_end)
  476. xs, xe = max(0, x_start), min(Nx, x_end)
  477. boundaries.append((xs, xe, ys, ye, z_idx))
  478. return boundaries
  481. def add_xy_outline(ax, boundaries, x, y, dx, dy, iz):
  482. boundary_for_z = next((b for b in boundaries if b[4] == iz), None)
  483. if boundary_for_z:
  484. xs, xe, ys, ye, _ = boundary_for_z
  485. rect = plt.Rectangle((x[xs] - dx / 2, y[ys] - dy / 2),
  486. (xe - xs) * dx,
  487. (ye - ys) * dy,
  488. edgecolor='white', facecolor='none', linewidth=2)
  489. ax.add_patch(rect)
  492. def add_xz_outline(ax, boundaries, x, z, dx, dz, iy):
  493. z_indices_in_slice = []
  494. x_range = None
  495. for xs, xe, ys, ye, z_idx in boundaries:
  496. if ys <= iy < ye:
  497. z_indices_in_slice.append(z_idx)
  498. if x_range is None:
  499. x_range = (xs, xe)
  500. if z_indices_in_slice and x_range:
  501. min_z_idx, max_z_idx = min(z_indices_in_slice), max(z_indices_in_slice)
  502. xs, xe = x_range
  503. rect = plt.Rectangle((x[xs] - dx / 2, z[min_z_idx] - dz / 2),
  504. (xe - xs) * dx,
  505. (max_z_idx - min_z_idx + 1) * dz,
  506. edgecolor='white', facecolor='none', linewidth=2)
  507. ax.add_patch(rect)
  510. def add_yz_outline(ax, boundaries, y, z, dy, dz, ix):
  511. for xs, xe, ys, ye, z_idx in boundaries:
  512. if xs <= ix < xe:
  513. rect = plt.Rectangle((y[ys] - dy / 2, z[z_idx] - dz / 2),
  514. (ye - ys) * dy,
  515. dz,
  516. edgecolor='white', facecolor='none', linewidth=2)
  517. ax.add_patch(rect)
  520. def style_axes(ax, xlabel, ylabel):
  521. ax.set_xlabel(xlabel, fontsize=LABEL_FONTSIZE)
  522. ax.set_ylabel(ylabel, fontsize=LABEL_FONTSIZE)
  523. ax.tick_params(labelsize=TICK_FONTSIZE)
  526. def encoding_label(name):
  527. return 'Positional Encoding' if name == 'positional' else 'Hash Encoding'
  530. def print_runtime_info(device):
  531. print('Runtime info:')
  532. print(f" Python = {sys.version.split()[0]}")
  533. print(f" NumPy = {np.__version__}")
  534. print(f" Matplotlib = {matplotlib.__version__}")
  535. print(f" PyTorch = {torch.__version__}")
  536. print(f" Device = {device}")
  537. print(f" CUDA = {torch.cuda.is_available()}")
  539. def run():
  540. set_seed(DATA_SEED)
  541. device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
  542. print_runtime_info(device)
  543. os.makedirs('plots', exist_ok=True)
  545. dx, dy, dz = DX, DY, DZ
  546. x = np.arange(0.0, X_MAX + dx, dx)
  547. y = np.arange(0.0, Y_MAX + dy, dy)
  548. z = np.arange(0.0, Z_MAX + dz, dz)
  549. Nx, Ny, Nz = len(x), len(y), len(z)
  551. Xc = x.astype(float)
  552. Yc = y.astype(float)
  553. Zc = z.astype(float)
  554. X3, Y3, Z3 = np.meshgrid(Xc, Yc, Zc, indexing='ij')
  555. grid_coords = np.stack([X3.ravel(), Y3.ravel(), Z3.ravel()], axis=1)
  557. c_mean = grid_coords.mean(axis=0, keepdims=True)
  558. c_std = grid_coords.std(axis=0, keepdims=True)
  559. coords_norm = (grid_coords - c_mean) / (c_std + 1e-12)
  560. coords_norm = torch.tensor(coords_norm, dtype=torch.float32, device=device, requires_grad=True)
  562. dz_half = dz / 2.0
  563. cell_grid = np.hstack([grid_coords, np.full((grid_coords.shape[0], 1), dz_half)])
  564. cell_grid = torch.tensor(cell_grid, dtype=torch.float32, device=device)
  566. XX, YY = np.meshgrid(x, y, indexing='ij')
  567. obs = np.column_stack([XX.ravel(), YY.ravel(), -np.ones(XX.size)])
  568. obs = torch.tensor(obs, dtype=torch.float32, device=device)
  570. print("Assembling sensitivity G ...")
  571. t0 = time.time()
  572. G = construct_sensitivity_matrix_G_torch(cell_grid, obs, dx, dy, device)
  573. G = G.clone().detach().requires_grad_(False)
  574. print(f"G shape = {tuple(G.shape)}, time = {time.time() - t0:.2f}s")
  576. rho_true_vec, rho_true_3d = make_block_model(Nx, Ny, Nz, dx, dy, dz, rho_bg=RHO_BG, rho_blk=RHO_BLK)
  577. rho_true_vec = rho_true_vec.to(device)
  579. with torch.no_grad():
  580. gz_true = (G @ rho_true_vec.unsqueeze(1)).squeeze(1)
  582. sigma = NOISE_LEVEL * gz_true.std()
  583. noise = sigma * torch.randn_like(gz_true)
  584. gz_obs = gz_true + noise
  585. Wd = 1.0 / sigma
  586. model_rng_state = capture_rng_state()
  588. cfg = dict(gamma=GAMMA, epochs=EPOCHS, lr=LR)
  589. cfg.update(
  590. use_early_stopping=USE_EARLY_STOPPING,
  591. early_stop_min_epochs=EARLY_STOP_MIN_EPOCHS,
  592. early_stop_patience=EARLY_STOP_PATIENCE,
  593. early_stop_target=EARLY_STOP_TARGET,
  594. early_stop_tol=EARLY_STOP_TOL,
  595. early_stop_overfit_patience=EARLY_STOP_OVERFIT_PATIENCE,
  596. )
  598. results = {}
  599. for encoding_type in ENCODING_ORDER:
  600. restore_rng_state(model_rng_state)
  601. print(f"\n▶ Encoding strategy: {encoding_type}")
  602. model = DensityContrastINR(
  603. encoding_type=encoding_type, hidden=HIDDEN, depth=DEPTH,
  604. rho_abs_max=RHO_ABS_MAX,
  605. **ENCODING_CONFIGS[encoding_type]
  606. ).to(device)
  607. opt = torch.optim.Adam(model.parameters(), lr=cfg['lr'])
  608. hist = train_inr(model, opt, coords_norm, G, gz_obs, Wd, Nx, Ny, Nz, dx, dy, dz, cfg)
  610. with torch.no_grad():
  611. m_inv = model(coords_norm).view(-1)
  612. gz_pred = (G @ m_inv.unsqueeze(1)).squeeze(1)
  614. results[encoding_type] = {
  615. 'm_inv': m_inv.detach().cpu().numpy().reshape(Nx, Ny, Nz),
  616. 'gz_pred': gz_pred.detach().cpu().numpy(),
  617. 'hist': hist,
  618. 'rms_rho': torch.sqrt(torch.mean((m_inv - rho_true_vec.to(device)) ** 2)).item(),
  619. 'rms_gz': torch.sqrt(torch.mean((gz_pred - gz_obs) ** 2)).item() * 1e5,
  620. }
  622. def get_axes_coords():
  623. x1d = grid_coords[:, 0].reshape(Nx, Ny, Nz)[:, 0, 0]
  624. y1d = grid_coords[:, 1].reshape(Nx, Ny, Nz)[0, :, 0]
  625. z1d = grid_coords[:, 2].reshape(Nx, Ny, Nz)[0, 0, :]
  626. return x1d, y1d, z1d
  628. x1d, y1d, z1d = get_axes_coords()
  629. block_boundaries = get_block_boundaries(Nx, Ny, Nz)
  631. ix, iy, iz = Nx // 2, Ny // 2, min(Nz - 1, 5)
  632. tru = rho_true_3d.cpu().numpy()
  633. inv_pos = results['positional']['m_inv']
  634. inv_hash = results['hash']['m_inv']
  636. tru_max = tru.max()
  637. inv_max = INV_VMAX
  639. fig1, axes = plt.subplots(4, 3, figsize=(16, 20))
  640. x1d, y1d, z1d = get_axes_coords()
  642. # cell-edge limits
  643. x_edge_min, x_edge_max = x1d[0] - dx/2, x1d[-1] + dx/2
  644. y_edge_min, y_edge_max = y1d[0] - dy/2, y1d[-1] + dy/2
  645. z_edge_min, z_edge_max = z1d[0] - dz/2, z1d[-1] + dz/2
  647. # use edges for all extents
  648. extent_xy = [x_edge_min, x_edge_max, y_edge_min, y_edge_max]
  649. # for depth plots keep depth increasing downward by reversing z limits
  650. extent_xz = [x_edge_min, x_edge_max, z_edge_max, z_edge_min]
  651. extent_yz = [y_edge_min, y_edge_max, z_edge_max, z_edge_min]
  654. im = axes[0, 0].imshow(tru[:, :, iz].T, origin='lower', extent=extent_xy,
  655. aspect='auto', vmin=0, vmax=tru_max, cmap=CMAP)
  656. axes[0, 0].set_title(f"True Model XY @ z≈{z1d[iz]:.0f} m", fontsize=TITLE_FONTSIZE)
  657. cbar = fig1.colorbar(im, ax=axes[0, 0], label='kg/m³', fraction=0.046, pad=0.04)
  658. cbar.ax.tick_params(labelsize=TICK_FONTSIZE)
  660. im = axes[0, 1].imshow(tru[:, iy, :].T, origin='upper', extent=extent_xz,
  661. aspect='auto', vmin=0, vmax=tru_max, cmap=CMAP)
  662. axes[0, 1].set_title(f"True Model XZ @ y≈{y1d[iy]:.0f} m", fontsize=TITLE_FONTSIZE)
  664. im = axes[0, 2].imshow(tru[ix, :, :].T, origin='upper', extent=extent_yz,
  665. aspect='auto', vmin=0, vmax=tru_max, cmap=CMAP)
  666. axes[0, 2].set_title(f"True Model YZ @ x≈{x1d[ix]:.0f} m", fontsize=TITLE_FONTSIZE)
  668. model_rows = [
  669. ('positional', inv_pos, 1),
  670. ('hash', inv_hash, 2),
  671. ]
  672. for encoding_type, inv, row in model_rows:
  673. title = encoding_label(encoding_type)
  675. im = axes[row, 0].imshow(inv[:, :, iz].T, origin='lower', extent=extent_xy,
  676. aspect='auto', vmin=0, vmax=inv_max, cmap=CMAP)
  677. axes[row, 0].set_title(f"Recovered Model XY, {title}", fontsize=TITLE_FONTSIZE)
  678. add_xy_outline(axes[row, 0], block_boundaries, x, y, dx, dy, iz)
  679. cbar = fig1.colorbar(im, ax=axes[row, 0], label='kg/m³', fraction=0.046, pad=0.04)
  680. cbar.ax.tick_params(labelsize=TICK_FONTSIZE)
  682. im = axes[row, 1].imshow(inv[:, iy, :].T, origin='upper', extent=extent_xz,
  683. aspect='auto', vmin=0, vmax=inv_max, cmap=CMAP)
  684. axes[row, 1].set_title(f"Recovered Model XZ, {title}", fontsize=TITLE_FONTSIZE)
  685. add_xz_outline(axes[row, 1], block_boundaries, x, z, dx, dz, iy)
  687. im = axes[row, 2].imshow(inv[ix, :, :].T, origin='upper', extent=extent_yz,
  688. aspect='auto', vmin=0, vmax=inv_max, cmap=CMAP)
  689. axes[row, 2].set_title(f"Recovered Model YZ, {title}", fontsize=TITLE_FONTSIZE)
  690. add_yz_outline(axes[row, 2], block_boundaries, y, z, dy, dz, ix)
  692. for ax in [axes[0, 1], axes[0, 2], axes[1, 1], axes[1, 2], axes[2, 1], axes[2, 2]]:
  693. ax.set_aspect(1.0)
  694. for ax in [axes[0, 0], axes[1, 0], axes[2, 0]]:
  695. ax.set_aspect(1.0)
  697. for row in range(3):
  698. style_axes(axes[row, 0], 'x (m)', 'y (m)')
  699. style_axes(axes[row, 1], 'x (m)', 'Depth (m)')
  700. style_axes(axes[row, 2], 'y (m)', 'Depth (m)')
  702. def to_mgal(g):
  703. return 1e5 * g.detach().cpu().numpy()
  704. obs_mgal = to_mgal(gz_obs)
  705. res_pos_mgal = obs_mgal - 1e5 * results['positional']['gz_pred']
  706. res_hash_mgal = obs_mgal - 1e5 * results['hash']['gz_pred']
  708. obs_x = obs[:, 0].cpu().numpy()
  709. obs_y = obs[:, 1].cpu().numpy()
  711. axes[3, 0].plot(results['positional']['hist']['gravity'], label='Positional', color='tab:red')
  712. axes[3, 0].plot(results['hash']['hist']['gravity'], label='Hash', color='black')
  713. axes[3, 0].set_title('Training Convergence by Encoding', fontsize=TITLE_FONTSIZE)
  714. axes[3, 0].set_yscale('log')
  715. axes[3, 0].grid(True, which='both', ls='--', alpha=0.3)
  716. axes[3, 0].legend(fontsize=LEGEND_FONTSIZE)
  717. style_axes(axes[3, 0], 'Epoch', 'Loss = mean((residual / sigma)^2)')
  719. vmax_res = max(np.abs(res_pos_mgal).max(), np.abs(res_hash_mgal).max())
  720. sc = axes[3, 1].scatter(obs_x, obs_y, c=res_pos_mgal, s=80, cmap=CMAP,
  721. vmin=-vmax_res, vmax=vmax_res, marker='o', edgecolors='none')
  722. axes[3, 1].set_title(
  723. f"Data Misfit, Positional Encoding (RMS = {np.sqrt(np.mean(res_pos_mgal**2)):.3f} mGal)",
  724. fontsize=TITLE_FONTSIZE,
  725. )
  726. cbar = fig1.colorbar(sc, ax=axes[3, 1], fraction=0.046, pad=0.04)
  727. cbar.ax.tick_params(labelsize=TICK_FONTSIZE)
  729. sc = axes[3, 2].scatter(obs_x, obs_y, c=res_hash_mgal, s=80, cmap=CMAP,
  730. vmin=-vmax_res, vmax=vmax_res, marker='o', edgecolors='none')
  731. axes[3, 2].set_title(
  732. f"Data Misfit, Hash Encoding (RMS = {np.sqrt(np.mean(res_hash_mgal**2)):.3f} mGal)",
  733. fontsize=TITLE_FONTSIZE,
  734. )
  735. cbar = fig1.colorbar(sc, ax=axes[3, 2], fraction=0.046, pad=0.04)
  736. cbar.ax.tick_params(labelsize=TICK_FONTSIZE)
  738. for ax in [axes[3, 1], axes[3, 2]]:
  739. style_axes(ax, 'x (m)', 'y (m)')
  740. ax.set_aspect('equal')
  742. fig1.tight_layout()
  743. fig1.savefig('plots/BlockModel_EncodingComparison.png', dpi=FIG_DPI)
  744. plt.close(fig1)
  746. for encoding_type in ENCODING_ORDER:
  747. label = encoding_label(encoding_type)
  748. hist = results[encoding_type]['hist']
  749. print(
  750. f"{label} RMS density-contrast error ≈ {results[encoding_type]['rms_rho']:.2f} kg/m^3"
  751. )
  752. print(
  753. f"{label} RMS data misfit ≈ {results[encoding_type]['rms_gz']:.3f} mGal | "
  754. f"best epoch = {hist['best_epoch']} | "
  755. f"best loss = {hist['best_weighted_mse']:.3f}"
  756. )
  758. if __name__ == '__main__':
  759. run()

001-EncodingComparisons_blocky.py at commit 3d5a05d, under MIT · at the source

Overview

Authors: Pankaj K Mishra1, Sanni Laaksonen1, Jochen Kamm1, Anand Singh2
ORCID iDs: Pankaj K Mishra
  1. Geological Survey of Finland,Vuorimiehentie 5, Espoo, Finland
  2. Indian Institute of Technology Bombay,Mumbai, India
Institutions: Geological Survey of Finland (Finland); Indian Institute of Technology Bombay (India)
Journal: Scientific reports, volume 16, issue 1, article 25630
Dates: received 14 October 2025; accepted 28 May 2026; published online 4 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41598-026-55960-5 · PMID 42243533 · PMCID PMC13478334 · OpenAlex W4416054661
Open access: gold, a free copy (OpenAlex)
Preprint: arxiv:2510.17876
Status: code verified
Categories: systems (subfield)
Methods: Machine learning
Keywords: Neural fields, Gravity, Physics-based deep learning, Scientific machine learning, Inversion, Mathematics and computing, Solid Earth sciences
Topic: Geophysical and Geoelectrical Methods (Geophysics, Earth and Planetary Sciences), according to OpenAlex
Funding: Research Council of Finland
Citations: not cited yet (Europe PMC); 39 references in the paper

Abstract

Inversion of gravity data is an important method for investigating subsurface density variations relevant to mineral exploration, geothermal assessment, carbon storage, natural hydrogen, groundwater resources, and tectonic evolution. Here we present a scientific machine-learning approach for three-dimensional gravity inversion that represents subsurface density as a continuous field using an implicit neural representation (INR). The method trains a deep neural network directly through a physics-based forward-model loss, mapping spatial coordinates to a continuous density field without predefined meshes or discretisation. Spatial encoding enhances the network’s capacity to capture sharp contrasts and short-wavelength features that conventional coordinate-based networks tend to oversmooth due to spectral bias. We demonstrate the approach on synthetic examples including smooth models, representing realistic geological complexity, and a dipping block model to assess recovery of structures at different depths. The INR framework reconstructs detailed structure and geologically plausible boundaries without explicit regularisation or depth weighting, while reducing the number of inversion parameters as the problem size grows bigger. These results highlight the potential of implicit representations to enable scalable, flexible, and interpretable large-scale geophysical inversion. This framework could generalise to other geophysical methods and for joint/multiphysics inversion.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 15 matches between paragraphs and lines of code.

Zenodo 19440024

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (12 files), NumPy (12 files), PyTorch (12 files)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
14 files
At the source: zenodo.org/records/19440024

pankajkmishra/inrgravity3dinv

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 3d5a05d1b1c59d909be7dc765f888ea5f2247849, 13 April 2026
Languages: Python (12)
Size: 16 files, 12 scripts
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, license file, environment (requirements.txt)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (12 files), NumPy (12 files), PyTorch (12 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
14 files
At the source: github.com/pankajkmishra/INRGravity3DInv/tree/v1.1

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 24 scripts, each with its path and the digest of its content;
  • 15 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

The codes for reproducing the results can be found at https://zenodo.org/records/19440024.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 7 keywords, 1 funder, 23 references.

Cite

This paper

Mishra, P. K., Laaksonen, S., Kamm, J., & Singh, A. (2026). Three-dimensional inversion of gravity data using implicit neural representations and scientific machine learning. Scientific reports, 16(1), 25630. https://doi.org/10.1038/s41598-026-55960-5

BibTeX

@article{mishra2026three,
author = {Mishra, Pankaj K and Laaksonen, Sanni and Kamm, Jochen and Singh, Anand},
title = {{Three-dimensional inversion of gravity data using implicit neural representations and scientific machine learning}},
journal = {Scientific reports},
year = {2026},
month = jun,
volume = {16},
number = {1},
pages = {25630},
publisher = {Nature Publishing Group},
issn = {2045-2322},
doi = {10.1038/s41598-026-55960-5},
url = {https://doi.org/10.1038/s41598-026-55960-5},
pmid = {42243533},
pmcid = {PMC13478334}
}

RIS

TY - JOUR
AU - Mishra, Pankaj K
AU - Laaksonen, Sanni
AU - Kamm, Jochen
AU - Singh, Anand
TI - Three-dimensional inversion of gravity data using implicit neural representations and scientific machine learning
T2 - Scientific reports
J2 - Sci Rep
PY - 2026
DA - 2026/06/04
VL - 16
IS - 1
SP - 25630
SN - 2045-2322
PB - Nature Publishing Group
DO - 10.1038/s41598-026-55960-5
UR - https://doi.org/10.1038/s41598-026-55960-5
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41598-026-55960-5",
"type": "article-journal",
"title": "Three-dimensional inversion of gravity data using implicit neural representations and scientific machine learning",
"container-title": "Scientific reports",
"author": [
{
"family": "Mishra",
"given": "Pankaj K"
},
{
"family": "Laaksonen",
"given": "Sanni"
},
{
"family": "Kamm",
"given": "Jochen"
},
{
"family": "Singh",
"given": "Anand"
}
],
"container-title-short": "Sci Rep",
"volume": "16",
"issue": "1",
"page": "25630",
"DOI": "10.1038/s41598-026-55960-5",
"PMID": "42243533",
"PMCID": "PMC13478334",
"ISSN": "2045-2322",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41598-026-55960-5",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
4
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-74092-y [code]
High-rate phase association with travel time neural fields.
Journal: Nature communications
In common: PyTorch, Matplotlib, NumPy, 1 reference
[2] doi:10.3390/jimaging12040170
ARS-GS: Anisotropic Reflective Spherical 3D Gaussian Splatting.
Journal: Journal of imaging
In common: 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.