Silentlog

Project.toml

@@ -1,7 +1,7 @@
 name = "HiQGA"
 uuid = "01574e87-63b6-4466-925a-334cf7e21b6b"
 authors = ["richardt94 <[email protected]>"]
-version = "0.4.7"
+version = "0.4.8"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -25,6 +25,7 @@ KernelDensitySJ = "49dc5b4e-6806-4504-b2cc-a81764e7a38b"
 LazyGrids = "7031d0ef-c40d-4431-b2f8-61a8d2f650db"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -50,7 +51,7 @@ Dates = "1.7"
 DelimitedFiles = "1.7"
 Distances = "0.10.7"
 Distributed = "1.7"
-DistributedArrays = "0.6.6"
+DistributedArrays = "0.6.7"
 Distributions = "0.25.58"
 FastGaussQuadrature = "0.4.9"
 FileIO = "1.10"

examples/VTEM/DarlingParoo/McMC/04_plot.jl

@@ -1,4 +1,4 @@
-dr = 20
+dr = 10
 idx = [1,2]
 burninfrac = 0.5
 vmin, vmax = -2.5, 0.5

examples/tempest/Menindee/McMC/01_read_data.jl

@@ -25,5 +25,5 @@ soundings = transD_GP.TEMPEST1DInversion.read_survey_files(
 	startfrom = 1,
 	skipevery = 1,
 	linenum = 1,
-	dotillsounding = 2,
+	dotillsounding = 4,
 	makeqcplots = true)

examples/tempest/Menindee/McMC/02_set_options.jl

@@ -1,5 +1,5 @@
 ## same for all soundingss
-nsamples = 201
+nsamples = 3001
 nchainspersounding = 3
 ppn = 4
 @assert mod(ppn, nchainspersounding+1) == 0

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

@@ -3,7 +3,7 @@ idxplot = [1,2]
 burninfrac = 0.5
 nbins = 50
 computeforwards, nforwards = true, 2
-interpolatedr = 10
+interpolatedr = 1
 vmax, vmin = .- extrema(fbounds)
 ## get the summary images of resistivity and nuisances over the profile
 transD_GP.AEMwithNuisanceMcMCInversionTools.summaryAEMwithnuisanceimages(soundings, opt, optn; zall,

src/AEM_VMD_HMD.jl

@@ -233,8 +233,10 @@ function checkrampformintime(times, ramp, minresptime, maxtime, rampchoice, isdI
             end    
         end
     end
-        mintime = max(0.5minta, minresptime) 
-        # though I believe mintime = minta always is safe.
+        mintime = max(0.25minta, minresptime) 
+        # I believe mintime = 0.25minta is safe.
+        # later times can cause errors in tb, unless extrapolating earlier
+        # than F.interptimes[1]. Extrapolation I'm hesitant to do.
     # if maxta > maxtime
         maxtime = min(1.5maxta, maxtime)
         # 1.5 to make sure we clear maxta in the interptimes

src/AEMnoNuisanceGradientInversionTools.jl

@@ -1,6 +1,6 @@
 module AEMnoNuisanceGradientInversionTools
 using Distributed, Dates, Printf, PyPlot, DelimitedFiles, StatsBase
-using ..AbstractOperator, ..CommonToAll, ..GP
+using ..AbstractOperator, ..CommonToAll, ..GP, ..SoundingDistributor
 import ..AbstractOperator.makeoperator
 import ..AbstractOperator.Sounding
 import ..AbstractOperator.returnforwrite
@@ -15,26 +15,26 @@ function compress(soundings, zall; prefix="", rmfile=true, isfirstparalleliterat
     fname = soundings[1].sounding_string*"_gradientinv.dat"    
     !isfile(fname) && throw(AssertionError("file does not exist perhaps soundings already zipped?"))
     fout = prefix == "" ? "zipped.dat" : prefix*"_zipped.dat"
-    iomode = "w"
-    if isfile(fout)
-        isfirstparalleliteration && throw(AssertionError("Zipped file "*fout*" exists, will not overwrite!"))
-        iomode = "a"
-    end    
-    io = open(fout, iomode)
+    isfirstparalleliteration && isfile(fout) && throw(AssertionError("Zipped file "*fout*" exists, will not overwrite!"))
     for (i, s) in enumerate(soundings)
         fname = s.sounding_string*"_gradientinv.dat"
         A = readdlm(fname)
         ϕd = A[end,2]
         σgrid = vec(A[end,3:end])
-        for el in [returnforwrite(s)...; vec(zall); σgrid; ϕd]
-            msg = @sprintf("%.4f\t", el)
-            write(io, msg)
-        end
-        write(io, "\n")                
-        flush(io) # slower but ensures write is complete
+        elinonerow = [returnforwrite(s)..., vec(zall), σgrid, ϕd]
+        nelinonerow = length(elinonerow)
+        writenames = [string.(s.writefields)..., "zcenter", "log10_cond", "phid_err"]
+        sfmt = fill("%15.3f", nelinonerow)
+        ϕd > 1e4 && (ϕd = 1e4 )
+        iomode = (isfirstparalleliteration && i==1) ? "w" : "a"
+        writeasegdat(elinonerow, sfmt, fout[1:end-4], iomode)
         rmfile && rm(fname)
+        if isfirstparalleliteration && i == 1
+            channel_names = [writenames, fill("", nelinonerow), writenames]
+            writeasegdfnfromonerow(elinonerow, channel_names, sfmt, fout[1:end-4])
+            dfn2hdr(fout[1:end-4]*".dfn")
+        end
     end
-    close(io)    
 end
 
 # plot the convergence and the result
@@ -237,21 +237,22 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0;
                             compresssoundings  = true,
                             zipsaveprefix      = "",
                             minimprovfrac      = nothing,
+                            verbose            = false,
                             minimprovkickinstep = round(Int, nstepsmax/2),        
                             ) where S<:Sounding
 
     @assert nsequentialiters  != -1
     nparallelsoundings = nworkers()
     nsoundings = length(soundings)
     zall = zboundarytocenter(aem_in.z[aem_in.nfixed+1:end]) # needed for sounding compression in write
+
+    writetogloballog("will require $nsequentialiters iterations of $nparallelsoundings soundings", iomode="w")
+    writetogloballog("starting sequential parallel iterations at $(Dates.now())")
     for iter = 1:nsequentialiters
-        if iter<nsequentialiters
-            ss = (iter-1)*nparallelsoundings+1:iter*nparallelsoundings
-        else
-            ss = (iter-1)*nparallelsoundings+1:nsoundings
-        end
-        @info "soundings in loop $iter of $nsequentialiters", ss
+        ss = getss_deterministic(iter, nsequentialiters, nparallelsoundings, nsoundings)
+        writetogloballog("soundings in loop $iter of $nsequentialiters $ss")
         pids = workers()
+        t2 = time()
         @sync for (i, s) in enumerate(ss)
             aem = makeoperator(aem_in, soundings[s])
             fname = soundings[s].sounding_string*"_gradientinv.dat"
@@ -279,15 +280,16 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0;
                                                 breakonknown       = breakonknown    ,              
                                                 dobo               = dobo            ,
                                                 fname              = fname           ,
-                                                minimprovfrac,
+                                                minimprovfrac, verbose,
                                                 minimprovkickinstep) 
 
 
         end # @sync
         isfirstparalleliteration = iter == 1 ? true : false
         compresssoundings && compress(soundings[ss[1]:ss[end]], zall, 
             isfirstparalleliteration = isfirstparalleliteration, prefix=zipsaveprefix)
-        @info "done $iter out of $nsequentialiters at $(Dates.now())"
+        dt = time() - t2 # seconds
+        writetogloballog("done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec")
     end
 end
 

src/AEMnoNuisanceMcMCInversionTools.jl

@@ -37,6 +37,7 @@ function summaryAEMimages(soundings::Array{S, 1}, opt::Options;
                         yl = nothing,
                         showplot = true,
                         showmean = false,
+                        Rmax = nothing,
                         logscale = true,
                         dpi = 300) where S<:Sounding
     compatidxwarn(idx, lnames)
@@ -48,7 +49,7 @@ function summaryAEMimages(soundings::Array{S, 1}, opt::Options;
         continueflag && continue
         summaryimages(soundings[a:b], opt; qp1, qp2, burninfrac, zall,dz, dr, 
             fontsize, vmin, vmax, cmap, figsize, topowidth, idx=idspec, omitconvergence, useML, 
-            preferEright, showplot, preferNright, saveplot, yl, dpi, showmean, logscale)
+            preferEright, showplot, preferNright, saveplot, yl, dpi, showmean, logscale, Rmax)
     end
     nothing    
 end
@@ -76,7 +77,8 @@ function summaryimages(soundings::Array{S, 1}, opt::Options;
                         logscale = true,
                         showplot = true,
                         showmean = false,
-                        dpi = 300) where S<:Sounding
+                        dpi = 300,
+                        Rmax=nothing) where S<:Sounding
     @assert !(preferNright && preferEright) # can't prefer both labels to the right
     pl, pm, ph, ρmean, χ²mean, χ²sd  = summarypost(soundings, opt; zall, qp1, qp2, burninfrac, useML)
 
@@ -85,9 +87,10 @@ function summaryimages(soundings::Array{S, 1}, opt::Options;
                                                         zall, dz; dr)
 
     lname = "Line_$(soundings[1].linenum)"
+
     plotsummarygrids1(soundings, σmeangrid, phgrid, plgrid, pmgrid, gridx, gridz, topofine, R, Z, χ²mean, χ²sd, lname; qp1, qp2,
                         figsize, fontsize, cmap, vmin, vmax, 
-                        topowidth, idx, omitconvergence, useML,
+                        topowidth, idx, omitconvergence, useML, Rmax,
                         preferEright, preferNright, saveplot, showplot, dpi,
                         yl, showmean, logscale)                  
 end
@@ -145,7 +148,7 @@ function processonesounding(opt_in::Options, sounding::Sounding, zall, burninfra
     ndata = getndata(sounding)
     χ²mean = mean(χ²)/ndata
     χ²sd   = std(χ²)/ndata
-    if hasproperty(sounding, :force_ML) 
+    if hasproperty(sounding, :forceML) 
         if sounding.forceML
             χ²mean = 1.
             χ²sd   = 0.
@@ -334,9 +337,10 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_
     nsoundings = length(soundings)
     nsequentialiters, nparallelsoundings = splittasks(soundings; nchainspersounding, ppn)
 
+    writetogloballog("starting sequential parallel iterations at $(Dates.now())")
     for iter = 1:nsequentialiters
         ss = getss(iter, nsequentialiters, nparallelsoundings, nsoundings)
-        @info "soundings in loop $iter of $nsequentialiters", ss
+        writetogloballog("soundings in loop $iter of $nsequentialiters $ss")
         t2 = time()
         @sync for (i, s) in Iterators.reverse(enumerate(ss))
             pids = getpids(i, nchainspersounding)
@@ -350,8 +354,8 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_
 
         end # @sync
         dt = time() - t2 #seconds
-        t2 = time()
-        @info "done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec"
+        catlocallogs(nparallelsoundings, nchainspersounding)
+        writetogloballog("done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec")
     end
 end
 

src/AEMwithNuisanceGradientInversionTools.jl

@@ -1,6 +1,6 @@
 module AEMwithNuisanceGradientInversionTools
 using Distributed, Dates, Printf, PyPlot, DelimitedFiles, StatsBase
-using ..AbstractOperator, ..CommonToAll, ..GP
+using ..AbstractOperator, ..CommonToAll, ..GP, ..SoundingDistributor
 import ..AbstractOperator.makeoperator
 import ..AbstractOperator.setnuboundsandstartforgradinv
 import ..AbstractOperator.Sounding
@@ -17,27 +17,27 @@ function compress(soundings, zall, nnu; prefix="", rmfile=true, isfirstparalleli
     fname = soundings[1].sounding_string*"_gradientinv.dat"    
     !isfile(fname) && throw(AssertionError("file does not exist perhaps soundings already zipped?"))
     fout = prefix == "" ? "zipped.dat" : prefix*"_zipped.dat"
-    iomode = "w"
-    if isfile(fout)
-        isfirstparalleliteration && throw(AssertionError("Zipped file "*fout*" exists, will not overwrite!"))
-        iomode = "a"
-    end    
-    io = open(fout, iomode)
+    isfirstparalleliteration && isfile(fout) && throw(AssertionError("Zipped file "*fout*" exists, will not overwrite!")) 
     for (i, s) in enumerate(soundings)
         fname = s.sounding_string*"_gradientinv.dat"
         A = readdlm(fname)
         ϕd = A[end,2]
         σgrid = vec(A[end,3:end-nnu])
         nu = vec(A[end,end-nnu+1:end])
-        for el in [returnforwrite(s)...; vec(zall); σgrid; nu; ϕd]
-            msg = @sprintf("%.4f\t", el)
-            write(io, msg)
-        end
-        write(io, "\n")                
-        flush(io) # slower but ensures write is complete
+        elinonerow = [returnforwrite(s)..., vec(zall), σgrid, nu, ϕd]
+        nelinonerow = length(elinonerow)
+        writenames = [string.(s.writefields)..., "zcenter", "log10_cond", "nuisance_param_in_order", "phid_err"]
+        sfmt = fill("%15.3f", nelinonerow)
+        ϕd > 1e4 && (ϕd = 1e4 )
+        iomode = (isfirstparalleliteration && i==1) ? "w" : "a"
+        writeasegdat(elinonerow, sfmt, fout[1:end-4], iomode)
         rmfile && rm(fname)
+        if isfirstparalleliteration && i == 1
+            channel_names = [writenames, fill("", nelinonerow), writenames]
+            writeasegdfnfromonerow(elinonerow, channel_names, sfmt, fout[1:end-4])
+            dfn2hdr(fout[1:end-4]*".dfn")
+        end
     end
-    close(io)    
 end
 
 # plot the convergence and the result
@@ -234,6 +234,9 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0, nu
                             knownvalue         = 0.7,
                             compresssoundings  = true,
                             zipsaveprefix      = "",
+                            minimprovfrac      = nothing,
+                            verbose            = false,
+                            minimprovkickinstep = round(Int, nstepsmax/2),
                             # optim stuff
                             ntriesnu           = 5,
                             boxiters           = 3, 
@@ -247,14 +250,14 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0, nu
     nparallelsoundings = nworkers()
     nsoundings = length(soundings)
     zall = zboundarytocenter(aem_in.z[aem_in.nfixed+1:end]) # needed for sounding compression in write
+
+    writetogloballog("will require $nsequentialiters iterations of $nparallelsoundings soundings", iomode="w")
+    writetogloballog("starting sequential parallel iterations at $(Dates.now())")
     for iter = 1:nsequentialiters
-        if iter<nsequentialiters
-            ss = (iter-1)*nparallelsoundings+1:iter*nparallelsoundings
-        else
-            ss = (iter-1)*nparallelsoundings+1:nsoundings
-        end
-        @info "soundings in loop $iter of $nsequentialiters", ss
+        ss = getss_deterministic(iter, nsequentialiters, nparallelsoundings, nsoundings)
+        writetogloballog("soundings in loop $iter of $nsequentialiters $ss")
         pids = workers()
+        t2 = time()
         @sync for (i, s) in enumerate(ss)
             aem = makeoperator(aem_in, soundings[s])
             nubounds, nustart = setnuboundsandstartforgradinv(aem, nuboundsΔ)
@@ -273,14 +276,16 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0, nu
                                                 regularizeupdate   = regularizeupdate,              
                                                 knownvalue         = knownvalue      ,              
                                                 fname              = fname           ,
-                                                ntriesnu, nubounds, boxiters, usebox, reducenuto, debuglevel, breaknuonknown) 
+                                                ntriesnu, nubounds, boxiters, usebox, reducenuto, debuglevel, breaknuonknown,
+                                                minimprovfrac, verbose, minimprovkickinstep) 
 
 
         end # @sync
         isfirstparalleliteration = iter == 1 ? true : false
         compresssoundings && compress(soundings[ss[1]:ss[end]], zall, size(nuboundsΔ, 1),
             isfirstparalleliteration = isfirstparalleliteration, prefix=zipsaveprefix)
-        @info "done $iter out of $nsequentialiters at $(Dates.now())"
+        dt = time() - t2 # seconds
+        writetogloballog("done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec")
     end
 end
 

src/AEMwithNuisanceMcMCInversionTools.jl

@@ -173,9 +173,10 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_
     nsoundings = length(soundings)
     nsequentialiters, nparallelsoundings = splittasks(soundings; nchainspersounding, ppn)
 
+    writetogloballog("starting sequential parallel iterations at $(Dates.now())")
     for iter = 1:nsequentialiters
         ss = getss(iter, nsequentialiters, nparallelsoundings, nsoundings)
-        @info "soundings in loop $iter of $nsequentialiters", ss
+        writetogloballog("soundings in loop $iter of $nsequentialiters $ss")
         t2 = time()
         @sync for (i, s) in Iterators.reverse(enumerate(ss))
             pids = getpids(i, nchainspersounding)
@@ -190,8 +191,8 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_
 
         end # @sync
         dt = time() - t2 #seconds
-        t2 = time()
-        @info "done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec"
+        catlocallogs(nparallelsoundings, nchainspersounding)
+        writetogloballog("done $iter out of $nsequentialiters at $(Dates.now()) in $dt sec")
     end
 end