Merge pull request #78 from GeoscienceAustralia/dev

Dev

Author: a2ray

Project.toml

@@ -1,7 +1,7 @@
 name = "HiQGA"
 uuid = "01574e87-63b6-4466-925a-334cf7e21b6b"
 authors = ["richardt94 <29562583+richardt94@users.noreply.github.com>"]
-version = "0.4.11"
+version = "0.4.12"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

src/AEMnoNuisanceMcMCInversionTools.jl

@@ -293,14 +293,23 @@ function makeAEMoperatorandoptions(sounding::Sounding;
                         showgeomplot = false,
                         plotfield = false,
                         useML = false,
-                        modelprimary = false
+                        modelprimary = false,
+                        isdIdt        = nothing,
+                        rampchoice    = nothing,
                         )
-    aem, zall, znall,  = makeoperator(sounding;
-                        zfixed, ρfixed, zstart, extendfrac,
-                        dz = dz, ρbg, nlayers, ntimesperdecade,
-                        nfreqsperdecade, useML, showgeomplot,
-                        modelprimary, plotfield)
-
+    if isnothing(isdIdt) # call VTEM or SkyTEM
+        aem, zall, znall,  = makeoperator(sounding;
+                            zfixed, ρfixed, zstart, extendfrac,
+                            dz = dz, ρbg, nlayers, ntimesperdecade,
+                            nfreqsperdecade, useML, showgeomplot,
+                            modelprimary, plotfield)
+    else # call VTEM for dIdt defined
+        aem, zall, znall,  = makeoperator(sounding; isdIdt, rampchoice,
+                            zfixed, ρfixed, zstart, extendfrac,
+                            dz = dz, ρbg, nlayers, ntimesperdecade,
+                            nfreqsperdecade, useML, showgeomplot,
+                            modelprimary, plotfield)
+    end
     opt = make_tdgp_opt(znall = znall,
                         fileprefix = sounding.sounding_string,
                         nmin = nmin,
@@ -357,7 +366,7 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_
         dt = time() - t2 #seconds
         catlocallogs(nparallelsoundings, nchainspersounding, batchstr)
         writetogloballog("done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec"; fname)
-        GC.gc()
+        # GC.gc()
     end
 end
 

src/AEMwithNuisanceMcMCInversionTools.jl

@@ -194,7 +194,7 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_
         dt = time() - t2 #seconds
         catlocallogs(nparallelsoundings, nchainspersounding, batchstr)
         writetogloballog("done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec"; fname)
-        GC.gc()
+        # GC.gc()
     end
 end
 

src/CommonToAll.jl

@@ -439,7 +439,7 @@ function main(opt_in     ::OptionsStat,
 
     close_history.([wp, wpns, wpn])
     d_closeall()
-    explicit_gc(master_pid, chainprocs)
+    # explicit_gc(master_pid, chainprocs)
     nothing
 end
 
@@ -557,7 +557,7 @@ function main(opt_in     ::OptionsStat,
 
     close_history.([wp, wpns])
     d_closeall()
-    explicit_gc(master_pid, chainprocs)
+    # explicit_gc(master_pid, chainprocs)
     nothing
 end
 
@@ -673,7 +673,7 @@ function main(opt_in       ::OptionsStat,
 
     close_history.([wp, wpn])
     d_closeall()
-    explicit_gc(master_pid, chainprocs)
+    # explicit_gc(master_pid, chainprocs)
     nothing
 end
 
@@ -768,7 +768,7 @@ function main(opt_in ::OptionsStat,
 
     close_history(wp)
     d_closeall()
-    explicit_gc(master_pid, chainprocs)
+    # explicit_gc(master_pid, chainprocs)
     nothing
 end
 

src/MCMC_Driver.jl

@@ -202,6 +202,8 @@ function read_survey_files(;
     startfrom = 1,
     skipevery = 1,
     multnoise = 0.03,
+    multnoise_LM = nothing,
+    multnoise_HM = nothing,
     datacutoff_LM = nothing,
     datacutoff_HM = nothing,
     X = -1,
@@ -228,6 +230,12 @@ function read_survey_files(;
     @assert Z > 0
     @assert linenum > 0
     @assert fid > 0
+    if !isnothing(multnoise_LM) | !isnothing(multnoise_HM)
+        @assert isnothing(multnoise) "set multnoise to nothing to be consistent!"
+    else
+        multnoise_LM = multnoise
+        multnoise_HM = multnoise
+    end        
     if forceML
         @assert !isnothing(datacutoff_LM)
         @assert !isnothing(datacutoff_HM)
@@ -277,8 +285,8 @@ function read_survey_files(;
     if !relerror
         @assert size(d_LM, 2) == length(LM_noise)
         @assert size(d_HM, 2) == length(HM_noise)
-        σ_LM = sqrt.((multnoise*d_LM).^2 .+ (LM_noise').^2)
-        σ_HM = sqrt.((multnoise*d_HM).^2 .+ (HM_noise').^2)
+        σ_LM = sqrt.((multnoise_LM*d_LM).^2 .+ (LM_noise').^2)
+        σ_HM = sqrt.((multnoise_HM*d_HM).^2 .+ (HM_noise').^2)
         if !isempty(noise_scalevec) 
             @assert length(noise_scalevec) == length(LM_times)+length(HM_times)
             σ_LM = σ_LM.*(noise_scalevec[1:length(LM_times)])'

src/SkyTEM1DInversion.jl

@@ -409,7 +409,7 @@ function makeoperator(sounding::VTEMsoundingData;
             calcjacobian  = false,
             plotfield     = false,
             isdIdt        = false,
-            rampchoice    = :mid,
+            rampchoice    = :next,
             )
 
     zall, znall, zboundaries = setupz(zstart, extendfrac, dz=dz, n=nlayers, showplot=showgeomplot)

src/VTEM1DInversion.jl

@@ -22,7 +22,7 @@ x, y = RDP.smoothline(XY1b, λ²=10000, finefactor=1000, regtype=:R2, fname="bas
 x, y = [[xy[i] for xy in RDP.readpoints(["base_menindee.txt", "base_1_menindee.txt", "base_2a_menindee.txt", "base_2b_menindee.txt", 
                     "base_3a_menindee.txt", "base_3b_menindee.txt"])] for i in 1:2]
 map(axp) do ax
-    for a in ax[1:6]
+    for a in ax[1:8]
         for (xx, yy) in zip(x, y)
             a.plot(xx, yy, "--k", linewidth=0.5)
         end