Merge pull request #59 from GeoscienceAustralia/convcheck

Convcheck to enable on ramp gates modeling

Project.toml

@@ -1,7 +1,7 @@
 name = "HiQGA"
 uuid = "01574e87-63b6-4466-925a-334cf7e21b6b"
 authors = ["richardt94 <[email protected]>"]
-version = "0.3.12"
+version = "0.3.13"
 
 [deps]
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"

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

@@ -15,6 +15,7 @@ d = [177,221]
 soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, 
 									fid, linenum, frame_height,
 									d, fname_dat, fname_specs_halt,
+									tx_rx_dz_pass_through = 0.01, # GA-AEM Z up, +ve is rx over tx
 									startfrom        = 1,
 									skipevery        = 1,
 									dotillsounding   = 2,

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

@@ -15,6 +15,7 @@ d = [177,221]
 soundings = transD_GP.VTEM1DInversion.read_survey_files(; X, Y, Z, 
 									fid, linenum, frame_height,
 									d, fname_dat, fname_specs_halt,
+									tx_rx_dz_pass_through = 0.01, # GA-AEM Z up, +ve is rx over tx
 									startfrom        = 1,
 									skipevery        = 50,
 									dotillsounding   = nothing,

examples/VTEM/DarlingParoo/gradientbased/electronics_halt.jl

@@ -50,5 +50,5 @@ times = vec(10 .^mean(log10.([
 			0.0100851	0.0113498
 ]), dims=2))
 include("hedgehog.jl")
-
+lowpassfcs = []
 σ_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/forwardtests/02_compare_ross.jl

@@ -13,11 +13,11 @@ z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
 rRx = 0.0
 zRx = -30.0
 zTx = -30.0
-freqlow = 1e-3
+freqlow = 1e-4
 freqhigh = 1e6
 nkᵣeval = 50
 ntimesperdecade = 10
-nfreqsperdecade = 6
+nfreqsperdecade = 5
 calcjacobian = true
 include("../waveletapprox/electronics_halt.jl")
 ## LM operator

examples/VTEM/synth/gradientbased/02_set_options.jl

@@ -7,4 +7,4 @@ regtype = :R2 # :R1 is first difference regularization, :R2 is two succesive fir
 lo, hi = -3, 1 # min, max log10 conductivity
 λ²min, λ²max=-0.5, 7 # min, max of regularization parameter (Tikhonov parameter)
 ntries = 8 # how many max tries between λ²min, λ²max
-β² = 0. # what fraction of roughness penalty to use to enforce a penalty in deviations from background
+β² = 0.001 # what fraction of roughness penalty to use to enforce a penalty in deviations from background

examples/VTEM/synth/waveletapprox/electronics_halt.jl

@@ -50,5 +50,5 @@ times = vec(10 .^mean(log10.([
 			0.0100851	0.0113498
 ]), dims=2))
 include("hedgehog.jl")
-
+lowpassfcs = []
 σ_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])

src/AEM_VMD_HMD.jl

@@ -21,6 +21,7 @@ mutable struct HFieldDHT <: HField
     ramp            :: Array{Float64, 2}
     log10ω          :: Array{Float64, 1}
     interptimes     :: Array{Float64, 1}
+    minresptime
     HFD_z           :: Array{ComplexF64, 1}
     HFD_r           :: Array{ComplexF64, 1}
     HFD_az          :: Array{ComplexF64, 1}
@@ -79,7 +80,7 @@ function HFieldDHT(;
       zTx       = -35.0,
       zRx       = -37.5,
       times     = 10 .^LinRange(-6, -1, 50),
-      ramp      = ones(10, 10),
+      ramp      = ones(10, 2),
       nfreqsperdecade = 5,
       ntimesperdecade = 10,
       glegintegorder = 5,
@@ -92,6 +93,7 @@ function HFieldDHT(;
       getazimH  = false,
       freqlow = 1e-4,
       freqhigh = 1e6,
+      minresptime = 1.e-6, # I think responses earlier than this are unstable
       calcjacobian = false
   )
     @assert all(freqs .> 0.)
@@ -103,20 +105,29 @@ function HFieldDHT(;
     ϵᵢ      = similar(pz)
     rTE       = zeros(length(pz)-1)
     rTM       = similar(rTE)
-    if freqhigh < 3/minimum(times)
-       freqhigh = 3/minimum(times)
+
+    ## all this below for pesky zero time at start instead of shutoff ...
+    mintime, maxtime = extrema(times)
+    @assert mintime > 0 # else interpolation in log10 trouble
+    mintime = 10^(log10(mintime) - 1) # go back a decade further than asked for
+    maxtime = 10^(log10(maxtime) + 1) # go ahead a decade further
+    if doconvramp
+         mintime, maxtime = checkrampformintime(times, ramp, minresptime, maxtime)
+    end
+    interptimes = 10 .^(log10(mintime) : 1/ntimesperdecade : log10(maxtime))
+    if freqhigh < 3/mintime
+        freqhigh = 3/mintime
     end
     if freqlow > 0.33/maximum(times)
        freqlow = 0.33/maximum(times)
     end
     if isempty(freqs)
-        freqs = 10 .^(log10(freqlow):1/nfreqsperdecade:log10(freqhigh))
+       freqs = 10 .^(log10(freqlow):1/nfreqsperdecade:log10(freqhigh))
     end
     interpkᵣ = 10 .^range(minimum(log10.(Filter_base))-0.5, maximum(log10.(Filter_base))+0.5, length = nkᵣeval)
     J0_kernel_h, J1_kernel_h, J0_kernel_v, J1_kernel_v = map(x->zeros(ComplexF64, length(interpkᵣ), length(freqs)), 1:4)
     J01kernelhold = zeros(ComplexF64, length(Filter_base))
-    log10ω = log10.(2*pi*freqs)
-    interptimes = 10 .^(minimum(log10.(times))-1:1/ntimesperdecade:maximum(log10.(times))+1)
+    log10ω = log10.(2*pi*freqs)    
     HFD_z       = zeros(ComplexF64, length(freqs)) # space domain fields in freq
     HFD_r       = zeros(ComplexF64, length(freqs)) # space domain fields in freq
     HFD_az       = zeros(ComplexF64, length(freqs)) # space domain fields in freq
@@ -161,7 +172,7 @@ function HFieldDHT(;
     #
     Jtemp = calcjacobian ? zeros(ComplexF64, nmax) : zeros(0)
     useprimary = modelprimary ? one(Float64) : zero(Float64)
-    HFieldDHT(thickness, pz, ϵᵢ, zintfc, rTE, rTM, zRx, zTx, rTx, rRx, freqs, times, ramp, log10ω, interptimes,
+    HFieldDHT(thickness, pz, ϵᵢ, zintfc, rTE, rTM, zRx, zTx, rTx, rRx, freqs, times, ramp, log10ω, interptimes, minresptime,
             HFD_z, HFD_r, HFD_az, HFD_z_interp, HFD_r_interp, HFD_az_interp,
             HTD_z_interp, HTD_r_interp, HTD_az_interp, dBzdt, dBrdt, dBazdt, J0_kernel_h, J1_kernel_h, J0_kernel_v, J1_kernel_v, J01kernelhold, lowpassfcs,
             quadnodes, quadweights, preallocate_ω_Hsc(interptimes, lowpassfcs)..., rxwithinloop, provideddt, doconvramp, useprimary,
@@ -170,6 +181,40 @@ function HFieldDHT(;
             HFD_r_J, HTD_r_J_interp, dBrdt_J)
 end
 
+function checkrampformintime(times, ramp, minresptime, maxtime)
+    # this checks we don't have ultra small tobs - ramp_time_a
+    minta = Inf
+    maxta = -Inf
+    for itime = 1:length(times)
+        for iramp = 1:size(ramp,1)-1
+            rta, rtb  = ramp[iramp,1], ramp[iramp+1,1]
+            if rta >= times[itime] # geq instead of gt as we could have an unlcky time
+                break
+            end
+            if rtb > times[itime] # end in this interval
+                rtb = times[itime]
+            end
+            # just so we know, rta < rtb and rta < t so rta < rtb <= t
+            ta = times[itime]-rta 
+            tb = max(times[itime]-rtb, minresptime) # rtb > rta, so make sure this is not zero because integ is in log10...
+            @assert ta>tb # else we're in trouble
+            if ta < minta
+                minta = ta
+            end
+            if ta > maxta
+                maxta = ta
+            end    
+        end
+    end
+        mintime = max(0.5minta, minresptime) 
+        # though I believe mintime = minta always is safe.
+    # if maxta > maxtime
+        maxtime = min(1.5maxta, maxtime)
+        # 1.5 to make sure we clear maxta in the interptimes
+    # end    
+    mintime, maxtime
+end
+
 #update geometry and dependent parameters - necessary for adjusting geometry
 #e.g. inverting nuisance parameters
 function update_ZR!(F::HFieldDHT, zTx, zRx, rTx, rRx; onlyprimary=false)
@@ -623,15 +668,15 @@ function convramp!(F::HFieldDHT, splz::CubicSpline, splr::CubicSpline, splaz::Cu
             dI   = F.ramp[iramp+1,2] - F.ramp[iramp,2]
             dIdt = dI/dt
 
-            if rta > F.times[itime]
+            if rta >= F.times[itime] # geq instead of eq as we could have an unlcky time
                 break
             end
             if rtb > F.times[itime] # end in this interval
                 rtb = F.times[itime]
             end
 
             ta = F.times[itime]-rta
-            tb = max(F.times[itime]-rtb, 1e-8) # rtb > rta, so make sure this is not zero...
+            tb = max(F.times[itime]-rtb, F.minresptime)# rtb > rta, so make sure this is not zero because integ is in log10...
             a, b = log10(ta), log10(tb)
             x, w = F.quadnodes, F.quadweights
             F.dBzdt[itime] += (b-a)/2*dot(getrampresponse((b-a)/2*x .+ (a+b)/2, splz), w)*dIdt

src/GradientInversion.jl

@@ -261,7 +261,12 @@ function gradientinv(   m::AbstractVector,
         istep += 1
         istep > nstepsmax && break
     end
-    isa(io, Nothing) || begin @info "Finished "*fname; close(io) end
+    if !isnothing(io) 
+        @info "Finished "*fname
+        close(io) 
+        nothing
+        return
+    end
     return map(x->x[1:istep], (mnew, χ², λ², oidx))
 end    
 
@@ -334,7 +339,12 @@ function gradientinv(   m::AbstractVector,
         istep += 1
         istep > nstepsmax && break
     end
-    isa(io, Nothing) || begin @info "Finished "*fname; close(io) end
+    if !isnothing(io) 
+        @info "Finished "*fname
+        close(io) 
+        nothing
+        return
+    end
     # nuisances mess up occam so this is a kludge as sometimes a smoothest root is not found
     (istep > nstepsmax) && (istep -= 1)
     return map(x->x[1:istep], (mnew, nunew, χ², χ²nu, λ², oidx))

src/TEMPEST1DInversion.jl

@@ -111,7 +111,7 @@ function Bfield(;
 	F = AEM_VMD_HMD.HFieldDHT(;nmax, calcjacobian,
 	                      ntimesperdecade = ntimesperdecade,
 	                      nfreqsperdecade = nfreqsperdecade,
-						  freqlow=1e-5,
+						#   freqlow =1e-5,
 						  nkᵣeval = nkᵣeval,
 	                      times  = times,
 	                      ramp   = ramp,

src/VTEM1DInversion.jl

@@ -34,6 +34,7 @@ function dBzdt(;times           = [1.],
         ramp            = [0. 1],
         rTx             = 10.,  
         zTx             = -30.,
+        zRx             = -30.01, # now Z down for use
         useML           = false,
         z               = [-1.],
         ρ               = [-1.],
@@ -46,14 +47,18 @@ function dBzdt(;times           = [1.],
         σ               = zeros(0),
         showgates       = false,
         modelprimary    = false,
+        lowpassfcs      = [],
+        freqlow         = 1e-4,
+        freqhigh        = 1e6, 
         )
 
     @assert size(σ)  == size(d)
     ndata, select = getndata(d)
 
     F = AEM_VMD_HMD.HFieldDHT(;
         times, ramp, nmax, zTx, rTx,
-        rRx = 0., zRx = zTx-0.01,
+        rRx = 0., zRx, 
+        lowpassfcs, freqlow, freqhigh,
         calcjacobian, nfreqsperdecade,
         ntimesperdecade, modelprimary
     )
@@ -100,6 +105,7 @@ mutable struct VTEMsoundingData <: Sounding
     fid             :: Float64
     linenum         :: Int
     zTx             :: Float64
+    zRx             :: Float64
     rTx             :: Float64
     lowpassfcs      :: Array{Float64, 1}
     times           :: Array{Float64, 1}
@@ -109,28 +115,29 @@ mutable struct VTEMsoundingData <: Sounding
 end
 
 returnforwrite(s::VTEMsoundingData) = [s.X, s.Y, s.Z, s.fid, 
-    s.linenum, s.zTx, s.rTx]
+    s.linenum, s.zTx, s.zRx, s.rTx]
 
 function getndata(S::VTEMsoundingData)
     getndata(S.data)[1]
 end    
 
 function VTEMsoundingData(;rRx=nothing, zRx=nothing, zTx=12.,
-                            rTx=-12.,lowpassfcs=[],
+                            rTx=-12.,lowpassfcs=[], 
+                            tx_rx_dz = -0.01, # not reading but using format Z down -ve is rx above tx
                             times=[1., 2.], ramp=[1 2; 3 4],
                             noise=[1.], data=[1.], 
                             sounding_string="sounding", X=nothing, Y=nothing, Z=nothing,
                             linenum=nothing, fid=nothing)
     @assert rTx > 0
     @assert zTx < 0 # my coordinate system z down 
     isnothing(rRx) && (rRx = 0.)
-    isnothing(zRx) && (zRx = zTx-0.1) # place receiver just above tx centre
+    isnothing(zRx) && (zRx = zTx + tx_rx_dz) 
     !isempty(lowpassfcs) && @assert all(lowpassfcs .> 0)
     @assert all(diff(times) .>0 )
     @assert all(diff(ramp[:,1]) .>0 )
     @assert all((noise .>0) .| isnan.(noise))
     @assert length(data) == length(noise)
-    VTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, zTx, rTx,
+    VTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, zTx, zRx, rTx,
         lowpassfcs, times, ramp, noise, data)
 end
 
@@ -152,6 +159,7 @@ function read_survey_files(;
     Z = -1,
     fid = -1,
     linenum = -1,
+    tx_rx_dz_pass_through = 0.01, # Z up GA-AEM reading convention +ve is rx above tx
     nanchar = "*")
 
     @assert frame_height > 0
@@ -171,24 +179,25 @@ function read_survey_files(;
     fiducial = soundings[:,fid]
     whichline = soundings[:,linenum]
     d = soundings[:,d[1]:d[2]]
-    zTx = -soundings[:,frame_height] # my coordinate system
-
+    zTx = soundings[:,frame_height] # read in Z up
+    zRx = -(zTx .+ tx_rx_dz_pass_through)  # my coordinate system Z down
+    zTx = -zTx # my coordinate system Z down
+
     @info "reading $fname_specs_halt"
     include(fname_specs_halt)
     @assert size(d, 2) == length(times)
     @assert size(d, 2) == length(σ_halt)
     σ_halt[:] .*= units
     d[:]      .*= units
     σ           = sqrt.((multnoise*d).^2 .+ (σ_halt').^2)
-
-    makeqcplots && plotsoundingdata(d, σ, times, zTx; figsize, fontsize)
+    makeqcplots && plotsoundingdata(d, σ, times, zTx, zRx; figsize, fontsize)
     nsoundings = size(d, 1)
     s_array = Array{VTEMsoundingData, 1}(undef, nsoundings)
     fracdone = 0 
     for is in 1:nsoundings
         l, f = Int(whichline[is]), fiducial[is]
-        s_array[is] = VTEMsoundingData(;zTx=zTx[is], rTx, 
-            times, ramp, noise=σ[is,:], data=d[is,:], 
+        s_array[is] = VTEMsoundingData(;zTx=zTx[is], zRx=zRx[is], rTx, 
+            times, ramp, noise=σ[is,:], data=d[is,:], lowpassfcs,
             sounding_string="sounding_$(l)_$f",
             X=easting[is], Y=northing[is], Z=topo[is], fid=f,
             linenum=l)
@@ -201,7 +210,7 @@ function read_survey_files(;
     return s_array
 end
 
-function plotsoundingdata(d, σ, times, zTx; figsize=(8,4), fontsize=1)
+function plotsoundingdata(d, σ, times, zTx, zRx; figsize=(8,4), fontsize=1)
     f, ax = plt.subplots(2, 2, figsize=figsize, gridspec_kw=Dict("width_ratios" => [1,0.01]))
     nsoundings = size(d, 1)
     plot_d = permutedims(d)
@@ -216,9 +225,11 @@ function plotsoundingdata(d, σ, times, zTx; figsize=(8,4), fontsize=1)
     axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "r")
     axx.semilogy(mean(σ./abs.(d), dims=1)[:], times, "--w")
     axx.set_xlabel("avg noise fraction")
+    ax[2].plot(1:nsoundings, zRx, label="zRx")
     ax[2].plot(1:nsoundings, zTx, label="zTx")
     ax[2].set_xlabel("sounding #")
     ax[2].set_ylabel("height m")
+    ax[2].invert_yaxis()
     ax[2].sharex(ax[1])
     ax[2].legend()
     ax[2,2].axis("off")
@@ -378,8 +389,8 @@ function makeoperator(sounding::VTEMsoundingData;
     z, ρ, = makezρ(zboundaries; zfixed, ρfixed)
     ρ[z.>=zstart] .= ρbg
     aem = dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary,
-        times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade,
-        rTx=sounding.rTx, zTx=sounding.zTx, z, ρ, calcjacobian, useML, showgates=plotfield)
+        times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade, lowpassfcs=sounding.lowpassfcs,
+        rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx, z, ρ, calcjacobian, useML, showgates=plotfield)
     plotfield && plotmodelfield!(aem, log10.(ρ[2:end]))
     aem, zall, znall, zboundaries
 end
@@ -388,9 +399,10 @@ function makeoperator(aem::dBzdt, sounding::VTEMsoundingData)
     ntimesperdecade = gettimesperdec(aem.F.interptimes)
     nfreqsperdecade = gettimesperdec(aem.F.freqs)
     modelprimary = aem.F.useprimary === 1. ? true : false
-    dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary,
+    dBzdt(;d=sounding.data/μ, σ=sounding.noise/μ, modelprimary, lowpassfcs=sounding.lowpassfcs,
         times=sounding.times, ramp=sounding.ramp, ntimesperdecade, nfreqsperdecade,
-        rTx=sounding.rTx, zTx=sounding.zTx, z=copy(aem.z), ρ=copy(aem.ρ), 
+        rTx=sounding.rTx, zTx=sounding.zTx, zRx=sounding.zRx,
+        z=copy(aem.z), ρ=copy(aem.ρ), 
         aem.F.calcjacobian, aem.useML, showgates=false)
 end