From d4d1aab89d8947a6847afa8ee299278a9173a6bf Mon Sep 17 00:00:00 2001
From: mloubout <mathias.louboutin@gmail.com>
Date: Tue, 12 Jul 2022 08:38:46 -0400
Subject: [PATCH] examples: switch to segy_scan for larger datasets

---
 examples/scripts/lsrtm_2D.jl   | 18 +++++++++++++-----
 examples/scripts/splsrtm_2D.jl |  9 ++++++---
 2 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/examples/scripts/lsrtm_2D.jl b/examples/scripts/lsrtm_2D.jl
index 79a555b91..3cca9ade0 100644
--- a/examples/scripts/lsrtm_2D.jl
+++ b/examples/scripts/lsrtm_2D.jl
@@ -19,7 +19,7 @@ model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)
 if ~isfile("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
     ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_2D.segy")
 end
-block = segy_read("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
+block = segy_scan(JUDI.JUDI_DATA, "marmousi_2D.segy", ["GroupX","GroupY","RecGroupElevation","SourceSurfaceElevation","dt"])
 d_lin = judiVector(block)   # linearized observed data
 
 # Set up wavelet
@@ -41,11 +41,17 @@ Mr = judiTopmute(model0.n, 52, 10)
 
 #' set up number of iterations
 niter = parse(Int, get(ENV, "NITER", "10"))
+# Default to 64, 5 for CI only with NITER=1
+nsrc = 5 * parse(Int, get(ENV, "NITER", "$(q.nsrc ÷ 5)"))
+indsrc = randperm(q.nsrc)[1:nsrc]
 lsqr_sol = zeros(Float32, prod(n))
 
 # LSQR
-Ml = judiMarineTopmute2D(30, d_lin.geometry)
-lsqr!(lsqr_sol, Ml*J*Mr, Ml*d_lin; maxiter=niter)
+dinv = d_lin[indsrc]
+Jinv = J[indsrc]
+
+Ml = judiMarineTopmute2D(30, dinv.geometry)
+lsqr!(lsqr_sol, Ml*Jinv*Mr, Ml*dinv; maxiter=niter)
 
 # Save final velocity model, function value and history
 h5open("lsrtm_marmousi_lsqr_result.h5", "w") do file
@@ -53,5 +59,7 @@ h5open("lsrtm_marmousi_lsqr_result.h5", "w") do file
 end
 
 # Plot final image
-figure(); imshow(copy(adjoint(reshape(lsqr_sol, model0.n))), extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
-title("LS-RTM with LSQR"); xlabel("Lateral position [km]"); ylabel("Depth [km]")
+figure(); imshow(reshape(lsqr_sol, model0.n)', extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
+title("LS-RTM with LSQR")
+xlabel("Lateral position [km]")
+ylabel("Depth [km]")
diff --git a/examples/scripts/splsrtm_2D.jl b/examples/scripts/splsrtm_2D.jl
index 8bafc7acf..36e7d942e 100644
--- a/examples/scripts/splsrtm_2D.jl
+++ b/examples/scripts/splsrtm_2D.jl
@@ -19,7 +19,7 @@ model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)
 if ~isfile("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
     ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/2DLSRTM/marmousi_2D.segy")
 end
-block = segy_read("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
+block = segy_scan(JUDI.JUDI_DATA, "marmousi_2D.segy", ["GroupX","GroupY","RecGroupElevation","SourceSurfaceElevation","dt"])
 d_lin = judiVector(block)   # linearized observed data
 
 # Set up wavelet
@@ -72,5 +72,8 @@ h5open("lsrtm_marmousi_breg_result.h5", "w") do file
 end
 
 # Plot final image
-figure(); imshow(copy(adjoint(reshape(solb.x, model0.n))), extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
-title("SPLS-RTM with Linearized Bregman"); xlabel("Lateral position [km]"); ylabel("Depth [km]")
\ No newline at end of file
+figure()
+imshow(reshape(solb.x, model0.n)', extent = (0, 7.99, 3.19, 0), cmap = "gray", vmin = -3e-2, vmax = 3e-2)
+title("SPLS-RTM with Linearized Bregman")
+xlabel("Lateral position [km]")
+ylabel("Depth [km]")
\ No newline at end of file