WIP ended? McMC works on VTEM

.gitignore

@@ -14,3 +14,4 @@ Manifest.toml
 *err*
 docs/build/
 docs/src/generated/
+*summary*.txt

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

@@ -3,7 +3,7 @@ idx = 1:100:length(soundings)
 burninfrac = 0.5
 vmin, vmax = -2.5, 0.5
 transD_GP.summaryAEMimages(soundings, opt;
-    zall, dr, dz, cmap="turbo",
+    zall, dr, dz, cmap="turbo", burninfrac, 
     vmin=vmin, vmax=vmax, idx, useML,
     preferEright=true, saveplot=true)
 ##

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

@@ -1,3 +1,4 @@
+# MPI Init
 using MPIClusterManagers, Distributed
 import MPI
 MPI.Init()

examples/VTEM/DarlingParoo/McMC/01_read_data.jl

@@ -0,0 +1,21 @@
+using HiQGA.transD_GP
+## info to read data
+# datafile
+fname_dat = "../gradientbased/200.asc"
+zipsaveprefix = splitpath(fname_dat)[end]
+# electronics file
+fname_specs_halt = "../gradientbased/electronics_halt.jl"
+# column numbers from hdr file
+X, Y, Z = 28, 29, 634
+fid = 5
+linenum = 4
+frame_height = 30
+d = [177,221]
+##
+soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, 
+									fid, linenum, frame_height,
+									d, fname_dat, fname_specs_halt,
+									startfrom        = 1,
+									skipevery        = 1,
+									dotillsounding   = 2,
+									makesounding     = true)

examples/VTEM/DarlingParoo/McMC/02_set_options.jl

@@ -0,0 +1,58 @@
+## same for all soundingss
+nsamples = 2001
+nchainspersounding = 5
+ppn = 48
+@assert mod(ppn, nchainspersounding+1) == 0
+zfixed   = [-1e5]
+ρfixed   = [1e12]
+zstart             = 0.0
+extendfrac         = 1.06
+dz                 = 1.5
+ρbg                = 10
+nlayers            = 50
+showgeomplot = true
+plotfield = true
+ntimesperdecade = 10
+nfreqsperdecade = 5
+useML = false
+## make transD options
+nmin, nmax = 2, 40
+K = transD_GP.GP.OrstUhn()
+demean = false
+sampledc = true
+fbounds = [-1 3.0]
+sdpos = 0.05
+sdprop = 0.05
+sddc = 0.01
+λ, δ = [2], 0.1
+save_freq = 100
+Tmax = 2.50
+nchainsatone = 1
+## plot 1 random sounding and a background response
+aem, opt, zall = transD_GP.makeAEMoperatorandoptions(
+                            soundings[rand(1:length(soundings))];
+                            zfixed = zfixed,
+                            ρfixed = ρfixed,
+                            zstart = zstart,
+                            extendfrac = extendfrac,
+                            dz = dz,
+                            ρbg = ρbg,
+                            nlayers = nlayers,
+                            ntimesperdecade = ntimesperdecade,
+                            nfreqsperdecade = nfreqsperdecade,
+                            useML = useML,
+                            nmin = nmin,
+                            nmax = nmax,
+                            K = K,
+                            demean = demean,
+                            sdpos = sdpos,
+                            sdprop = sdprop,
+                            sddc = sddc,
+                            sampledc = sampledc,
+                            fbounds = fbounds,
+                            save_freq = save_freq,
+                            showgeomplot = showgeomplot,
+                            plotfield,
+                            λ = λ,
+                            δ = δ
+)

examples/VTEM/DarlingParoo/McMC/03_parallel_run.jl

@@ -0,0 +1,35 @@
+## MPI Init
+using MPIClusterManagers, Distributed
+import MPI
+MPI.Init()
+rank = MPI.Comm_rank(MPI.COMM_WORLD)
+sz = MPI.Comm_size(MPI.COMM_WORLD)
+if rank == 0
+    @info "size is $sz"
+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()
+@assert mod(ncores+1,nchainspersounding+1) == 0
+@assert mod(ppn, nchainspersounding+1) == 0
+nparallelsoundings = Int((ncores+1)/(nchainspersounding+1))
+nsequentialiters = ceil(Int, nsoundings/nparallelsoundings)
+@info "will require $nsequentialiters iterations of $nparallelsoundings soundings in parallel"
+## set up McMC
+@everywhere using Distributed
+@everywhere using HiQGA.transD_GP
+## do the parallel soundings
+@info "starting"
+transD_GP.loopacrossAEMsoundings(soundings, aem, opt;
+                    nsequentialiters, nparallelsoundings,
+                    Tmax, nsamples, nchainsatone, nchainspersounding )
+
+MPIClusterManagers.stop_main_loop(manager)
+rmprocs(workers())
+exit()

examples/VTEM/DarlingParoo/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", burninfrac,
+    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/VTEM/DarlingParoo/gradientbased/01_read_data.jl

@@ -0,0 +1,21 @@
+using HiQGA.transD_GP
+## info to read data
+# datafile
+fname_dat = "200.asc"
+zipsaveprefix = splitpath(fname_dat)[end]
+# electronics file
+fname_specs_halt = "electronics_halt.jl"
+# column numbers from hdr file
+X, Y, Z = 28, 29, 634
+fid = 5
+linenum = 4
+frame_height = 30
+d = [177,221]
+##
+soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, 
+									fid, linenum, frame_height,
+									d, fname_dat, fname_specs_halt,
+									startfrom        = 1,
+									skipevery        = 50,
+									dotillsounding   = nothing,
+									makesounding     = true)

examples/VTEM/DarlingParoo/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         = 1.06
+dz                 = 1.5
+ρbg                = 10
+nlayers            = 50
+ntimesperdecade    = 10
+nfreqsperdecade    = 5
+regtype            = :R1
+nstepsmax          = 40
+ntries             = 6
+lo                 = -3.
+hi                 = 1.
+λ²min              = -0.5
+λ²max              = 8
+β²                 = 0.01
+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/VTEM/DarlingParoo/gradientbased/03_parallel_run.jl

@@ -0,0 +1,33 @@
+## MPI Init
+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/VTEM/DarlingParoo/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"

examples/VTEM/synth/waveletapprox/electronics_halt.jl

@@ -0,0 +1,54 @@
+# Tx dimensions
+rTx = sqrt(530.92/pi)
+# halt_noise
+# gate times
+times = vec(10 .^mean(log10.([
+    0.0000180	0.0000230
+			0.0000230	0.0000290
+			0.0000290	0.0000340
+			0.0000340	0.0000390
+			0.0000390	0.0000450
+			0.0000450	0.0000510
+			0.0000510	0.0000590
+			0.0000590	0.0000680
+			0.0000680	0.0000780
+			0.0000780	0.0000900
+			0.0000900	0.0001030
+			0.0001030	0.0001180
+			0.0001180	0.0001360
+			0.0001360	0.0001560
+			0.0001560	0.0001790
+			0.0001790	0.0002060
+			0.0002060	0.0002360
+			0.0002360	0.0002710
+			0.0002710	0.0003120
+			0.0003120	0.0003580
+			0.0003580	0.0004110
+			0.0004110	0.0004720
+			0.0004720	0.0005430
+			0.0005430	0.0006230
+			0.0006230	0.0007160
+			0.0007160	0.0008230
+			0.0008230	0.0009450
+			0.0009450	0.0010860
+			0.0010860	0.0012470
+			0.0012470	0.0014320
+			0.0014320	0.0016460
+			0.0016460	0.0018910
+			0.0018910	0.0021720
+			0.0021720	0.0024950
+			0.0024950	0.0028650
+			0.0028650	0.0032920
+			0.0032920	0.0037810
+			0.0037810	0.0043410
+			0.0043410	0.0049870
+			0.0049870	0.0057290
+			0.0057290	0.0065810
+			0.0065810	0.0075600
+			0.0075600	0.0086850
+			0.0086850	0.0099770
+			0.0100851	0.0113498
+]), dims=2))
+include("hedgehog.jl")
+
+σ_halt = vec([0.238388  0.175366  0.126135  0.091645  0.071083  0.043692  0.030069  0.019933  0.015175  0.010694  0.008489  0.006127  0.005354  0.005058  0.003871  0.003524  0.002759  0.002185  0.001792  0.001505  0.001298  0.001112  0.000946  0.000805  0.000689  0.000620  0.000546  0.000465  0.000428  0.000399  0.000352  0.000333  0.000317  0.000293  0.000302  0.000334  0.000406  0.000465  0.000458  0.000356  0.000292  0.000193  0.000285  0.000195  0.000177])

examples/VTEM/synth/waveletapprox/hedgehog.jl

@@ -0,0 +1,25 @@
+ramp = [
+-0.007328119307	0.00281467
+-0.006828119307	0.39957362
+-0.006510410974	0.38746451
+-0.005999994307	0.6922533
+-0.005666660974	0.67181256
+-0.005182285974	0.83379579
+-0.004822910974	0.80791275
+-0.004348952641	0.88168368
+-0.003989577641	0.85425795
+-0.003515619307	0.92503455
+-0.003156244307	0.89625896
+-0.002682285974	0.96429834
+-0.002317702641	0.93442465
+-0.001848952641	0.9998768
+-0.001484369307	0.96887822
+-0.001322910974	0.9887351
+-0.000999994307	0.8390506
+-0.000901035974	0.77591249
+-0.000781244307	0.69091416
+-0.000614577641	0.55932957
+-0.000322910974	0.30213728
+5.693e-9	0.0
+5.214026e-6	0.0
+]

src/AEMnoNuisanceMcMCInversionTools.jl

@@ -3,6 +3,7 @@ using ..AbstractOperator, ..CommonToAll
 import ..AbstractOperator.makeoperator
 import ..AbstractOperator.loopacrossAEMsoundings
 import ..AbstractOperator.plotmodelfield!
+import ..AbstractOperator.getndata
 import ..Options, ..OptionsStat
 export makeAEMoperatorandoptions, loopacrossAEMsoundings, summaryAEMimages, plotindividualAEMsoundings
 import ..main # McMC function
@@ -114,8 +115,7 @@ function summarypost(soundings::Vector{S}, opt::Options;
                                                     qp1=qp1, qp2=qp2,
                                                     doplot=false)
             χ² = 2*CommonToAll.assembleTat1(opt, :U, temperaturenum=1, burninfrac=burninfrac)
-            ndata = sum(.!isnan.(soundings[idx].LM_data)) +
-                    sum(.!isnan.(soundings[idx].HM_data))
+            ndata = getndata(soundings[idx])
             χ²mean[idx] = mean(χ²)/ndata
             χ²sd[idx]   = std(χ²)/ndata
             if useML
@@ -137,6 +137,12 @@ function summarypost(soundings::Vector{S}, opt::Options;
     pl, pm, ph, ρmean, vdmean, vddev, χ²mean, χ²sd
 end
 
+function getndata(d)
+    select = .!isnan.(d)
+    ndata  = sum(select)
+    ndata, select
+end
+
 function plotindividualAEMsoundings(soundings::Vector{S}, aem_in::Operator1D, opt::Options, idxplot::Vector{Int};
     burninfrac=0.5,
     nbins = 50,

src/AbstractOperator.jl

@@ -24,6 +24,8 @@ function make_tdgp_opt end
 function loopacrossAEMsoundings end
 # to plot forwards from a sounding
 function plotmodelfield! end
+# to get number of data from sounding
+function getndata end
 
 export Operator, Operator1D, Operator2D, Sounding
 end

src/SkyTEM1DInversion.jl

@@ -4,6 +4,7 @@ import ..AbstractOperator.Sounding
 import ..AbstractOperator.makeoperator
 import ..AbstractOperator.getresidual
 import ..AbstractOperator.returnforwrite
+import ..AbstractOperator.getndata
 import ..AbstractOperator.plotmodelfield!
 using ..AbstractOperator, ..AEM_VMD_HMD, Statistics, Distributed, Printf, Dates, StatsBase,
       PyPlot, LinearAlgebra, ..CommonToAll, Random, DelimitedFiles, LinearMaps, SparseArrays
@@ -118,12 +119,6 @@ function getresidual(aem::dBzdt, log10σ::Vector{Float64}; computeJ=false)
     nothing    
 end
 
-function getndata(d)
-    select = .!isnan.(d)
-    ndata  = sum(select)
-    ndata, select
-end
-
 mutable struct SkyTEMsoundingData <: Sounding
     sounding_string :: String
     X :: Float64
@@ -152,6 +147,10 @@ returnforwrite(s::SkyTEMsoundingData) = [s.X, s.Y, s.Z, s.fid,
     s.linenum, s.rRx, s.zRxLM, s.zTxLM, s.zRxHM, 
     s.zTxHM, s.rTx]
 
+function getndata(S::SkyTEMsoundingData)
+    getndata(S.LM_data)[2] + getndata(S.HM_data)[2]
+end   
+
 function SkyTEMsoundingData(;rRx=-12., zRxLM=12., zTxLM=12.,
                             zRxHM=12., zTxHM=12., rTx=-12.,lowpassfcs=[-1, -2.],
                             LM_times=[1., 2.], LM_ramp=[1 2; 3 4],