WIP: SkyTEM all seems to work

examples/SkyTEM1D/Menindee/02_set_options.jl → examples/SkyTEM1D/Menindee/McMC/02_set_options.jl

@@ -28,7 +28,7 @@ save_freq = 100
 Tmax = 2.50
 nchainsatone = 1
 ## plot 1 random sounding and a background response
-aem, opt = transD_GP.makeAEMoperatorandoptions(
+aem, opt, zall = transD_GP.makeAEMoperatorandoptions(
                             soundings[rand(1:length(soundings))];
                             zfixed = zfixed,
                             ρfixed = ρfixed,

examples/SkyTEM1D/Menindee/03_parallel_run.jl → examples/SkyTEM1D/Menindee/McMC/03_parallel_run.jl

@@ -27,12 +27,8 @@ nsequentialiters = ceil(Int, nsoundings/nparallelsoundings)
 ## do the parallel soundings
 @info "starting"
 transD_GP.loopacrossAEMsoundings(soundings, aem, opt;
-                    nsequentialiters   = nsequentialiters,
-                    nparallelsoundings = nparallelsoundings,
-                    Tmax               = Tmax,
-                    nsamples           = nsamples,
-                    nchainsatone       = nchainsatone,
-                    nchainspersounding = nchainspersounding)
+                    nsequentialiters, nparallelsoundings,
+                    Tmax, nsamples, nchainsatone, nchainspersounding )
 
 MPIClusterManagers.stop_main_loop(manager)
 rmprocs(workers())

examples/SkyTEM1D/Menindee/McMC/04_plot.jl

@@ -0,0 +1,16 @@
+dr = 20
+idx = 1:100:length(soundings)
+burninfrac = 0.5
+vmin, vmax = -2.5, 0.5
+transD_GP.summaryAEMimages(soundings, opt;
+    zall, dr, dz, cmap="turbo",
+    vmin=vmin, vmax=vmax, idx, useML,
+    preferEright=true, saveplot=true)
+##
+computeforwards = true
+nforwards = 5
+nbins = 50
+transD_GP.plotindividualAEMsoundings(soundings, aem, opt, idx;
+    burninfrac, nbins, computeforwards,
+    nforwards)
+

examples/SkyTEM1D/Menindee/em_with_noise_calibration_lines.hdr → examples/SkyTEM1D/Menindee/McMC/em_with_noise_calibration_lines.hdr

@@ -1,31 +1,31 @@
-1	fiducial
-2	line
-3	flight
-4	datetime
-5	date
-6	utc_time
-7	angleX
-8	angleY
-9	frame_height
-10	longitude
-11	latitude
-12	easting
-13	northing
-14	ground_elevation_aem
-15	ground_elevation_lidar_aem_merged
-16	gps_height
-17	groundspeed
-18	frame_roll
-19	frame_pitch
-20	frame_yaw
-21	frame_dx
-22	frame_dy
-23	frame_dz
-24	HM_current
-25	LM_current
-26-42	LM_Z
-43-67	HM_Z
-68-84	LM_Z_rerr
-85-109	HM_Z_rerr
-110-126	LM_Z_noise
-127-151	HM_Z_noise
+1	fiducial
+2	line
+3	flight
+4	datetime
+5	date
+6	utc_time
+7	angleX
+8	angleY
+9	frame_height
+10	longitude
+11	latitude
+12	easting
+13	northing
+14	ground_elevation_aem
+15	ground_elevation_lidar_aem_merged
+16	gps_height
+17	groundspeed
+18	frame_roll
+19	frame_pitch
+20	frame_yaw
+21	frame_dx
+22	frame_dy
+23	frame_dz
+24	HM_current
+25	LM_current
+26-42	LM_Z
+43-67	HM_Z
+68-84	LM_Z_rerr
+85-109	HM_Z_rerr
+110-126	LM_Z_noise
+127-151	HM_Z_noise

examples/SkyTEM1D/Menindee/gradientbased/01_read_data.jl

@@ -0,0 +1,39 @@
+## info to read data
+# datafile
+fname_dat = "../McMC/coincidewithtempest.dat"
+zipsaveprefix = splitpath(fname_dat)[end]
+# electronics file
+fname_specs_halt = "../McMC/electronics_halt.jl"
+# column numbers from hdr file
+X, Y, Z = 12, 13, 15
+fid = 1
+linenum = 2
+frame_height = 9
+frame_dx = 21
+frame_dy = 22
+frame_dz = 23
+LM_Z = [26,42]
+HM_Z = [43,67]
+relerror = false
+units = 1e-12
+##
+using HiQGA.transD_GP
+soundings = transD_GP.SkyTEM1DInversion.read_survey_files(fname_dat = fname_dat,
+						             fname_specs_halt = fname_specs_halt,
+						             LM_Z             = LM_Z,
+									 HM_Z             = HM_Z,
+									 frame_height     = frame_height,
+									 frame_dz         = frame_dz,
+									 frame_dy         = frame_dy,
+									 frame_dx         = frame_dx,
+									 X                = X,
+									 Y                = Y,
+									 Z                = Z,     
+									 fid              = fid,
+									 linenum          = linenum,
+                                     startfrom        = 1,
+									 skipevery        = 50,
+									 dotillsounding   = nothing,
+									 relerror         = relerror,
+									 units            = units,     
+									 makesounding     = true)

examples/SkyTEM1D/Menindee/gradientbased/02_set_options.jl

@@ -0,0 +1,27 @@
+## make gradientinv options same for all soundings
+σstart, σ0         = -2, -2
+zfixed             = [-1e5]
+ρfixed             = [1e12]
+zstart             = 0.0
+extendfrac, dz     = 1.06, 1.25
+nlayers            = 52
+ρbg                = 10
+ntimesperdecade    = 10
+nfreqsperdecade    = 5
+regtype            = :R1
+nstepsmax          = 40
+ntries             = 6
+lo                 = -2.5
+hi                 = 0.5
+λ²min              = -0.5
+λ²max              = 8
+β²                 = 0.1
+knownvalue         = 0.7
+breakonknown       = true
+calcjacobian       = true
+## plot a random sounding and a background response
+aem, zall, = transD_GP.VTEM1DInversion.makeoperator(
+    soundings[rand(1:length(soundings))];
+    zfixed, ρfixed, zstart, extendfrac, calcjacobian,
+    dz, ρbg, nlayers, ntimesperdecade, nfreqsperdecade, plotfield=true)
+

examples/SkyTEM1D/Menindee/gradientbased/03_parallel_run.jl

@@ -0,0 +1,32 @@
+using MPIClusterManagers, Distributed
+import MPI
+MPI.Init()
+rank = MPI.Comm_rank(MPI.COMM_WORLD)
+size = MPI.Comm_size(MPI.COMM_WORLD)
+if rank == 0
+    @info "size is $size"
+end
+manager = MPIClusterManagers.start_main_loop(MPI_TRANSPORT_ALL)
+@info "there are $(nworkers()) workers"
+@everywhere @info gethostname()
+include("01_read_data.jl")
+include("02_set_options.jl")
+## MPI checks
+# split into sequential iterations of parallel soundings
+nsoundings = length(soundings)
+ncores = nworkers()
+nsequentialiters = ceil(Int, nsoundings/ncores)
+@info "will require $nsequentialiters iterations of $ncores soundings"
+## set up transD_GP
+@everywhere using Distributed
+@everywhere using HiQGA.transD_GP
+## do the parallel soundings
+@info "starting"
+transD_GP.loopacrossAEMsoundings(soundings, aem, σstart, σ0;
+                            nsequentialiters, regtype, nstepsmax, ntries,           
+                            lo, hi, λ²min, λ²max, β², knownvalue, breakonknown,              
+                            zipsaveprefix)                  
+
+MPIClusterManagers.stop_main_loop(manager)
+rmprocs(workers())
+exit()

examples/SkyTEM1D/Menindee/gradientbased/04_plot.jl

@@ -0,0 +1,6 @@
+transD_GP.plotconvandlast(soundings, 50, 2; zall, 
+    figsize=(15,7), vmin=lo, vmax=hi, postfix="_β²_$(β²)_$(regtype)_bg_$(round(10. ^σ0, sigdigits=4))Spm", prefix=zipsaveprefix,
+    preferEright=true, showplot=true, logscale=true)
+@info "writing PDF ..."
+run(`./pdfscript.sh`)
+@info "done"

src/AEMnoNuisanceGradientInversionTools.jl

@@ -1,6 +1,6 @@
 module AEMnoNuisanceGradientInversionTools
 using Distributed, Dates, Printf, PyPlot, DelimitedFiles, StatsBase
-using ..AbstractOperator, ..CommonToAll
+using ..AbstractOperator, ..CommonToAll, ..GP
 import ..AbstractOperator.makeoperator
 import ..AbstractOperator.Sounding
 import ..AbstractOperator.returnforwrite
@@ -36,7 +36,7 @@ end
 
 # plot the convergence and the result
 function plotconvandlast(soundings, delr, delz; 
-        zstart=-1, extendfrac=-1, dz=-1, nlayers=-1, cmapσ="turbo", vmin=-2.5, vmax=0.5, fontsize=12,
+        zall=[-1.], cmapσ="turbo", vmin=-2.5, vmax=0.5, fontsize=12,
         figsize=(20,5),
         topowidth=1,
         preferEright = false,
@@ -53,8 +53,8 @@ function plotconvandlast(soundings, delr, delz;
         dpi=400)
     linestartidx = splitsoundingsbyline(soundings)                    
     nlines = length(linestartidx)
+    nlayers = length(zall)
     fnamecheck = soundings[1].sounding_string*"_gradientinv.dat"
-    zall, = setupz(zstart, extendfrac, dz=dz, n=nlayers)
     isfile(fnamecheck) && compress(soundings, zall) # write everything in one file if not done yet
     fzipped = prefix == "" ? "zipped.dat" : prefix*"_zipped.dat"
     A = readdlm(fzipped)
@@ -161,18 +161,12 @@ function plotconvandlasteachline(soundings, σ, ϕd, delr, delz;
 end    
 
 # driver for gradient based AEM inversion
-function loopacrossAEMsoundings(soundings::Array{S, 1}, σstart, σ0; 
+function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0; 
                             nsequentialiters   =-1,
-                            zfixed             = [-1e5],
-                            ρfixed             = [1e12],
                             zstart             = 0.0,
                             extendfrac         = 1.06,
                             dz                 = 2.,
-                            ρbg                = 10,
                             nlayers            = 50,
-                            ntimesperdecade    = 10,
-                            nfreqsperdecade    = 5,
-                            modelprimary       = false,
                             regtype            = :R1,
                             nstepsmax          = 10,
                             ntries             = 6,
@@ -211,19 +205,7 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, σstart, σ0;
         @info "soundings in loop $iter of $nsequentialiters", ss
         pids = workers()
         @sync for (i, s) in enumerate(ss)
-            aem, = makeoperator(    soundings[s],
-                                    zfixed = zfixed,
-                                    ρfixed = ρfixed,
-                                    zstart = zstart,
-                                    extendfrac = extendfrac,
-                                    dz = dz,
-                                    ρbg = ρbg,
-                                    nlayers = nlayers,
-                                    modelprimary = modelprimary,
-                                    ntimesperdecade = ntimesperdecade,
-                                    nfreqsperdecade = nfreqsperdecade,
-                                    calcjacobian = true)
-
+            aem = makeoperator(aem_in, soundings[s])
             fname = soundings[s].sounding_string*"_gradientinv.dat"
             σstart_, σ0_ = map(x->x*ones(length(aem.ρ)-1), [σstart, σ0])
             @async remotecall_wait(gradientinv, pids[i], σstart_, σ0_, aem,