From 41b281725628dc6a14d5c29f38b5c95520141d74 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Tue, 17 Sep 2024 16:41:54 +1000 Subject: [PATCH 01/22] WIP:MLnoise for SkyTEM soundings --- src/SkyTEM1DInversion.jl | 49 ++++++++++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 9 deletions(-) diff --git a/src/SkyTEM1DInversion.jl b/src/SkyTEM1DInversion.jl index 7d95a08..bf5948f 100644 --- a/src/SkyTEM1DInversion.jl +++ b/src/SkyTEM1DInversion.jl @@ -5,6 +5,7 @@ import ..AbstractOperator.makeoperator import ..AbstractOperator.getresidual import ..AbstractOperator.returnforwrite import ..AbstractOperator.getndata +import ..AbstractOperator.checkifdatalow import ..AbstractOperator.plotmodelfield! using ..AbstractOperator, ..AEM_VMD_HMD, Statistics, Distributed, Printf, Dates, StatsBase, PyPlot, LinearAlgebra, ..CommonToAll, Random, DelimitedFiles, LinearMaps, SparseArrays @@ -143,6 +144,7 @@ mutable struct SkyTEMsoundingData <: Sounding HM_noise :: Array{Float64, 1} LM_data :: Array{Float64, 1} HM_data :: Array{Float64, 1} + forceML end returnforwrite(s::SkyTEMsoundingData) = [s.X, s.Y, s.Z, s.fid, @@ -159,7 +161,7 @@ function SkyTEMsoundingData(;rRx=-12., zRxLM=12., zTxLM=12., HM_times=[1., 2.], HM_ramp=[1 2; 3 4], LM_noise=[1.], HM_noise=[1.], LM_data=[1.], HM_data=[1.], sounding_string="sounding", X=nothing, Y=nothing, Z=nothing, - linenum=nothing, fid=nothing) + linenum=nothing, fid=nothing, forceML=false) @assert rRx > 0 && rTx > 0 @assert zRxLM <0 && zTxLM <0 @assert zRxHM <0 && zTxHM <0 @@ -174,7 +176,7 @@ function SkyTEMsoundingData(;rRx=-12., zRxLM=12., zTxLM=12., @assert length(HM_data) == length(HM_noise) SkyTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, rRx, zRxLM, zTxLM, zRxHM, zTxHM, rTx, lowpassfcs, LM_times, LM_ramp, HM_times, HM_ramp, LM_noise, HM_noise, - LM_data, HM_data) + LM_data, HM_data, forceML) end function read_survey_files(; @@ -208,7 +210,10 @@ function read_survey_files(; Z = -1, fid = -1, linenum = -1, - nanchar = "*") + nanchar = "*", + lineslessthan = nothing, + forceML = false, # for very low amp, in conjunction with datacutoff_LM and datacutoff_HM + ) @assert frame_height > 0 @assert (frame_dz > 0) | !isnothing(tx_rx_dz_pass_through) @assert (frame_dx > 0) | !isnothing(tx_rx_dx_pass_through) @@ -224,6 +229,13 @@ function read_survey_files(; @assert Z > 0 @assert linenum > 0 @assert fid > 0 + if forceML + @assert !isnothing(datacutoff_LM) + @assert !isnothing(datacutoff_HM) + end + if !isnothing(lineslessthan) + lineslessthan::Int + end @info "reading $fname_dat" soundings = readlargetextmatrix(fname_dat, startfrom, skipevery, dotillsounding) soundings[soundings .== nanchar] .= NaN @@ -273,12 +285,12 @@ function read_survey_files(; σ_LM = σ_LM.*(noise_scalevec[1:length(LM_times)])' σ_HM = σ_HM.*(noise_scalevec[1:length(HM_times)])' end - if !isnothing(datacutoff_LM) + if !isnothing(datacutoff_LM) && !forceML # since my dBzdt is +ve idxbad = d_LM .< datacutoff_LM d_LM[idxbad] .= NaN end - if !isnothing(datacutoff_HM) + if !isnothing(datacutoff_HM) && !forceML # since my dBzdt is +ve idxbad = d_HM .< datacutoff_HM d_HM[idxbad] .= NaN @@ -293,10 +305,16 @@ function read_survey_files(; s_array = Array{SkyTEMsoundingData, 1}(undef, nsoundings) fracdone = 0 + countforceML = 0 for is in 1:nsoundings idxbadz[is] && continue l, f = Int(whichline[is]), fiducial[is] + if !isnothing(lineslessthan) + l > lineslessthan && continue # skips high_alt and repeat lines if specified + end dlow, dhigh = vec(d_LM[is,:]), vec(d_HM[is,:]) + lowampflag = checkifdatalow(dlow, dhigh, datacutoff_LM, datacutoff_HM) + countforceML += lowampflag ? 1 : 0 s_array[is] = SkyTEMsoundingData(rRx=rRx[is], zRxLM=zRx[is], zTxLM=zTx[is], zRxHM=zRx[is], zTxHM=zTx[is], rTx=rTx, lowpassfcs=lowpassfcs, LM_times=LM_times, LM_ramp=LM_ramp, @@ -304,14 +322,16 @@ function read_survey_files(; LM_noise=σ_LM[is,:], HM_noise=σ_HM[is,:], LM_data=dlow, HM_data=dhigh, sounding_string="sounding_$(l)_$f", X=easting[is], Y=northing[is], Z=topo[is], fid=f, - linenum=l) + linenum=l, forceML=lowampflag) fracnew = round(Int, is/nsoundings*100) if (fracnew-fracdone)>10 fracdone = fracnew @info "read $is out of $nsoundings" end end - return s_array[.!idxbadz] + idx = [isassigned(s_array, i) for i in 1:length(s_array)] + forceML && @info("low-amplitude forceML on $countforceML out of $(sum(idx)): $(round(100*countforceML/sum(idx)))%") + return s_array[idx] end function plotsoundingdata(nsoundings, LM_times, HM_times, d_LM, d_HM, @@ -359,7 +379,16 @@ function plotsoundingdata(nsoundings, LM_times, HM_times, d_LM, d_HM, ax[4].set_xlabel("sounding #") ax[4].set_ylabel("rRx m") plt.tight_layout() -end +end + +function checkifdatalow(d_LM, d_HM, datacutoff_LM, datacutoff_HM) + lowflag = false + if (mean(abs.(d_LM)) < datacutoff_LM) | (mean(abs.(d_HM)) < datacutoff_HM) + lowflag = true + end + lowflag +end + # all calling functions here for misfit, field, etc. assume model is in log10 resistivity # SANS the top. For lower level field calculation use AEM_VMD_HMD structs @@ -528,8 +557,10 @@ end function makeoperator(aem::dBzdt, sounding::SkyTEMsoundingData) ntimesperdecade = gettimesperdec(aem.Flow.interptimes) nfreqsperdecade = gettimesperdec(aem.Flow.freqs) + @assert !(aem.useML & sounding.forceML) "useML and forceML cannot both be true" + useML = (aem.useML | sounding.forceML), # OR logic for useML with true, true disqualified earlier modelprimary = aem.Flow.useprimary === 1. ? true : false - makeoperator(sounding, ntimesperdecade, nfreqsperdecade, modelprimary, aem.Flow.calcjacobian, aem.useML, copy(aem.z), copy(aem.ρ)) + makeoperator(sounding, ntimesperdecade, nfreqsperdecade, modelprimary, aem.Flow.calcjacobian, useML, copy(aem.z), copy(aem.ρ)) end # all plotting codes here assume that the model is in log10 resistivity, SANS From b230bb60f0c7897ec84fa4e65e4463b232ee9721 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Wed, 18 Sep 2024 14:42:54 +1000 Subject: [PATCH 02/22] modifications for forceML in SkyTEM low amp --- src/AEMnoNuisanceMcMCInversionTools.jl | 2 +- src/SkyTEM1DInversion.jl | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/AEMnoNuisanceMcMCInversionTools.jl b/src/AEMnoNuisanceMcMCInversionTools.jl index 9838568..a66cd82 100644 --- a/src/AEMnoNuisanceMcMCInversionTools.jl +++ b/src/AEMnoNuisanceMcMCInversionTools.jl @@ -145,7 +145,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. diff --git a/src/SkyTEM1DInversion.jl b/src/SkyTEM1DInversion.jl index bf5948f..d61a4ac 100644 --- a/src/SkyTEM1DInversion.jl +++ b/src/SkyTEM1DInversion.jl @@ -5,7 +5,6 @@ import ..AbstractOperator.makeoperator import ..AbstractOperator.getresidual import ..AbstractOperator.returnforwrite import ..AbstractOperator.getndata -import ..AbstractOperator.checkifdatalow import ..AbstractOperator.plotmodelfield! using ..AbstractOperator, ..AEM_VMD_HMD, Statistics, Distributed, Printf, Dates, StatsBase, PyPlot, LinearAlgebra, ..CommonToAll, Random, DelimitedFiles, LinearMaps, SparseArrays @@ -313,7 +312,7 @@ function read_survey_files(; l > lineslessthan && continue # skips high_alt and repeat lines if specified end dlow, dhigh = vec(d_LM[is,:]), vec(d_HM[is,:]) - lowampflag = checkifdatalow(dlow, dhigh, datacutoff_LM, datacutoff_HM) + lowampflag = checkifdatalow(dlow, dhigh, datacutoff_LM, datacutoff_HM, forceML) countforceML += lowampflag ? 1 : 0 s_array[is] = SkyTEMsoundingData(rRx=rRx[is], zRxLM=zRx[is], zTxLM=zTx[is], zRxHM=zRx[is], zTxHM=zTx[is], rTx=rTx, lowpassfcs=lowpassfcs, @@ -381,9 +380,11 @@ function plotsoundingdata(nsoundings, LM_times, HM_times, d_LM, d_HM, plt.tight_layout() end -function checkifdatalow(d_LM, d_HM, datacutoff_LM, datacutoff_HM) +function checkifdatalow(d_LM, d_HM, datacutoff_LM, datacutoff_HM, forceML) lowflag = false - if (mean(abs.(d_LM)) < datacutoff_LM) | (mean(abs.(d_HM)) < datacutoff_HM) + !forceML && return lowflag + mLM, mHM = map(xx->mean(log10.(abs.(xx))), [d_LM, d_HM]) + if (mLM < log10(datacutoff_LM)) || (mHM < log10(datacutoff_HM)) lowflag = true end lowflag @@ -558,7 +559,7 @@ function makeoperator(aem::dBzdt, sounding::SkyTEMsoundingData) ntimesperdecade = gettimesperdec(aem.Flow.interptimes) nfreqsperdecade = gettimesperdec(aem.Flow.freqs) @assert !(aem.useML & sounding.forceML) "useML and forceML cannot both be true" - useML = (aem.useML | sounding.forceML), # OR logic for useML with true, true disqualified earlier + useML = (aem.useML | sounding.forceML) # OR logic for useML with true, true disqualified earlier modelprimary = aem.Flow.useprimary === 1. ? true : false makeoperator(sounding, ntimesperdecade, nfreqsperdecade, modelprimary, aem.Flow.calcjacobian, useML, copy(aem.z), copy(aem.ρ)) end From dfdac1ca21d8485453daa03d1050b91d2374a17a Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Sat, 21 Sep 2024 17:26:48 +1000 Subject: [PATCH 03/22] MLnoise SkyTEM plot changes --- src/AEMnoNuisanceMcMCInversionTools.jl | 9 ++-- src/CommonToAll.jl | 70 ++++++++++++++------------ 2 files changed, 43 insertions(+), 36 deletions(-) diff --git a/src/AEMnoNuisanceMcMCInversionTools.jl b/src/AEMnoNuisanceMcMCInversionTools.jl index a66cd82..2575c87 100644 --- a/src/AEMnoNuisanceMcMCInversionTools.jl +++ b/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 diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index d867b30..1b9ac08 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1669,40 +1669,44 @@ function plotsummarygrids1(soundings, meangrid, phgrid, plgrid, pmgrid, gridx, g if isnothing(Rmax) Rmax = maximum(gridx) end - idx_split, thereisrmn = getRsplits(gridx, Rmax) - nimages = length(idx_split) - dr = gridx[2] - gridx[1] - if iszero(idx_split[1]) && thereisrmn # Rmax is larger than section - nx = length(range(gridx[1], Rmax, step=dr)) - elseif !iszero(idx_split[1]) && !thereisrmn # Rmax is exactly the section length - nx = length(gridx) - elseif iszero(idx_split[2]) # Rmax is smaller than the section - nx = idx_split[1] - else - nx = idx_split[2]-idx_split[1] # There are many Rmax length splits - end - ab, nimages = getinterpsplits(idx_split, nimages, gridx) - i_idx = 1:nimages - for i in i_idx - a, b = ab[i,1], ab[i,2] - a_uninterp = i == firstindex(i_idx) ? 1 : findlast(R.<=gridx[a]) - b_uninterp = i != lastindex(i_idx) ? findlast(R.<=gridx[b]) : lastindex(soundings) - if thereisrmn && i == lastindex(i_idx) - xrangelast = range(gridx[a], step=dr, length=nx) - else - xrangelast = nothing + while true + idx_split, thereisrmn = getRsplits(gridx, Rmax) + nimages = length(idx_split) + dr = gridx[2] - gridx[1] + if iszero(idx_split[1]) && thereisrmn # Rmax is larger than section + nx = length(range(gridx[1], Rmax, step=dr)) + elseif !iszero(idx_split[1]) && !thereisrmn # Rmax is exactly the section length + nx = length(gridx) + elseif iszero(idx_split[2]) # Rmax is smaller than the section + nx = idx_split[1] + else + nx = idx_split[2]-idx_split[1] # There are many Rmax length splits + end + ab, nimages = getinterpsplits(idx_split, nimages, gridx) + i_idx = 1:nimages + for i in i_idx + a, b = ab[i,1], ab[i,2] + a_uninterp = i == firstindex(i_idx) ? 1 : findlast(R.<=gridx[a]) + b_uninterp = i != lastindex(i_idx) ? findlast(R.<=gridx[b]) : lastindex(soundings) + if thereisrmn && i == lastindex(i_idx) + xrangelast = range(gridx[a], step=dr, length=nx) + else + xrangelast = nothing + end + f, s, icol = setupconductivityplot(gridx[a:b], omitconvergence, showmean, R[a_uninterp:b_uninterp], + figsize, fontsize, lname, χ²mean[a_uninterp:b_uninterp], χ²sd[a_uninterp:b_uninterp], useML, i, nimages, logscale) + + summaryconductivity(s, icol, f, soundings[a_uninterp:b_uninterp], + meangrid[:,a:b], phgrid[:,a:b], plgrid[:,a:b], pmgrid[:,a:b], + gridx[a:b], gridz, topofine[a:b], R[a_uninterp:b_uninterp], Z[a_uninterp:b_uninterp], ; qp1, qp2, fontsize, + cmap, vmin, vmax, topowidth, idx, omitconvergence, preferEright, preferNright, yl, showmean, xrangelast) + postfix = nimages == 1 ? "_whole" : "_split_$(i)_of_$(nimages)" + saveplot && savefig(lname*postfix*".png", dpi=dpi) + showplot || close(f) end - f, s, icol = setupconductivityplot(gridx[a:b], omitconvergence, showmean, R[a_uninterp:b_uninterp], - figsize, fontsize, lname, χ²mean[a_uninterp:b_uninterp], χ²sd[a_uninterp:b_uninterp], useML, i, nimages, logscale) - - summaryconductivity(s, icol, f, soundings[a_uninterp:b_uninterp], - meangrid[:,a:b], phgrid[:,a:b], plgrid[:,a:b], pmgrid[:,a:b], - gridx[a:b], gridz, topofine[a:b], R[a_uninterp:b_uninterp], Z[a_uninterp:b_uninterp], ; qp1, qp2, fontsize, - cmap, vmin, vmax, topowidth, idx, omitconvergence, preferEright, preferNright, yl, showmean, xrangelast) - - saveplot && savefig(lname*"_split_$(i)_of_$(nimages).png", dpi=dpi) - showplot || close(f) - end + nimages == 1 && break # don't need another pass if it was whole line to start with + Rmax = maximum(gridx) # this will ensure nimages=1 on next pass + end end function setupconductivityplot(gridx, omitconvergence, showmean, R, figsize, fontsize, lname, χ²mean, χ²sd, useML, iimage, nimages, logscale) From c27e9f031cbaabda923a68b0f1141c1f6173db52 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Wed, 2 Oct 2024 12:19:49 +1000 Subject: [PATCH 04/22] small fix for no improvement --- src/GradientInversion.jl | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/GradientInversion.jl b/src/GradientInversion.jl index 3b7a85e..bd8a466 100644 --- a/src/GradientInversion.jl +++ b/src/GradientInversion.jl @@ -78,6 +78,8 @@ function occamstep(m::AbstractVector, m0::AbstractVector, Δm::AbstractVector, m end if all(χ² .>= χ²₀) α = count < countmax - 1 ? α/2 : 0. # make sure we don't start going uphill again + else # if any χ² < what we started with + count = countmax end count += 1 count > countmax && break @@ -86,7 +88,7 @@ function occamstep(m::AbstractVector, m0::AbstractVector, Δm::AbstractVector, m if all(χ² .> target) idx = argmin(χ²) else - idx = findlast(χ² .<= target) + idx = findfirst(χ² .<= target) a = log10(λ²[idx][1]) b = λ²max # isa(bidx, Nothing) ? λ²max : log10(λ²[bidx][1]) # maybe just set default λ²max dophase2(m, m0, F, R, regularizeupdate, lo, hi, target, a, b, mnew, χ², λ², β², idx) @@ -275,7 +277,7 @@ function gradientinv( m::AbstractVector, end noimprovement = iszero(λ²[istep][idx][2]) ? true : false if (istep == nstepsmax - 1) || noimprovement || littleimprovement - target = χ²[istep][idx] # exit with smoothest + target = χ²[istep][idx] # exit with smoothest if possible end istep += 1 istep > nstepsmax && break @@ -312,6 +314,8 @@ function gradientinv( m::AbstractVector, debuglevel = 0, usebox = true, boxiters = 2, + minimprovfrac = nothing, + minimprovkickinstep = round(Int, nstepsmax/2), fname="") R = makereg(regtype, F) ndata = length(F.res) @@ -327,6 +331,7 @@ function gradientinv( m::AbstractVector, ndata = length(F.res) istep = 1 io = open_history(fname) + littleimprovement = false while true # Optim stuff for nuisances # f(x) = 2*get_misfit(m, x, F, nubounds) # why? always set usebox=true @@ -354,9 +359,15 @@ function gradientinv( m::AbstractVector, oidx[istep] = idx isa(io, Nothing) || write_history(io, [istep; χ²[istep][idx]/target₀; vec(m); nu]) foundroot && break + if !isnothing(minimprovfrac) && (istep > minimprovkickinstep) + prevχ²= χ²[istep-1][oidx[istep-1]] + if (χ²[istep][idx] - prevχ²)/prevχ² < minimprovfrac + littleimprovement = true + end + end noimprovement = iszero(λ²[istep][idx][2]) ? true : false - if (istep == nstepsmax - 1) || noimprovement - target = χ²[istep][idx] # exit with smoothest + if (istep == nstepsmax - 1) || noimprovement || littleimprovement + target = χ²[istep][idx] # exit with smoothest if possible end istep += 1 istep > nstepsmax && break From 5e23ae7445503f2b3c0c5a6a8d942c9f65863d5f Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Wed, 2 Oct 2024 15:24:55 +1000 Subject: [PATCH 05/22] WIP: dfn for dat files --- src/AEMnoNuisanceGradientInversionTools.jl | 8 ++++++++ src/AEMwithNuisanceGradientInversionTools.jl | 9 +++++++++ src/CommonToAll.jl | 2 +- src/VTEM1DInversion.jl | 9 +++++---- 4 files changed, 23 insertions(+), 5 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index 6385901..273398d 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -33,6 +33,14 @@ function compress(soundings, zall; prefix="", rmfile=true, isfirstparalleliterat write(io, "\n") flush(io) # slower but ensures write is complete rmfile && rm(fname) + if isfirstparalleliteration && i == 1 + elinonerow = [returnforwrite(s), vec(zall), σgrid, ϕd] + nelinonerow = length(elinonerow) + writenames = [[f for f in fieldnames(typeof(s))], "zcenter", "log10_cond", "ϕd_err"] + sfmt = fill("%15.3f", nelinonerow) + channel_names = [writenames, fill("", nelinonerow), writenames] + writeasegdfnfromonerow(elinonerow, channel_names, sfmt, fout[1:end-4]) + end end close(io) end diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl index f33c4d6..56dacf2 100644 --- a/src/AEMwithNuisanceGradientInversionTools.jl +++ b/src/AEMwithNuisanceGradientInversionTools.jl @@ -36,6 +36,15 @@ function compress(soundings, zall, nnu; prefix="", rmfile=true, isfirstparalleli write(io, "\n") flush(io) # slower but ensures write is complete rmfile && rm(fname) + if i == 1 + elinonerow = [returnforwrite(s), vec(zall), σgrid, nu, ϕd] + nelinonerow = length(elinonerow) + writenames = [[f for f in fieldnames(typeof(a))], "zcenter", "log10_cond", "nuisance_param", "ϕd_err"] + sfmt = fill("%15.3f", nelinonerow) + channel_names = [writenames, fill("", nelinonerow, writenames)] + writeasegdfnforgradinv(elinonerow, fname[1:end-4]) + writeasegdfnfromonerow(elinonerow, channel_names, sfmt, fname[1:end-4]) + end end close(io) end diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 1b9ac08..97801e7 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -18,7 +18,7 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo summaryconductivity, plotsummarygrids1, getVE, writevtkfromsounding, trapezoidrule, readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, writeijkfromsounding, nanmean, infmean, nanstd, infstd, infnanmean, infnanstd, - kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, + kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, writeasegdfnfromonerow, getprobabilisticoutputs,readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getRandgridr, getallxyinr, getXYlast, getprobabilisticlinesfromdirectory, readxyzrhoϕnu, plotgausshist, findMAE, getclosestidx, getpostatXY, getpostatwell diff --git a/src/VTEM1DInversion.jl b/src/VTEM1DInversion.jl index 54154ac..1c3939d 100644 --- a/src/VTEM1DInversion.jl +++ b/src/VTEM1DInversion.jl @@ -117,10 +117,10 @@ mutable struct VTEMsoundingData <: Sounding ramp :: Array{Float64, 2} noise :: Array{Float64, 1} data :: Array{Float64, 1} + writefields end -returnforwrite(s::VTEMsoundingData) = [s.X, s.Y, s.Z, s.fid, - s.linenum, s.zTx, s.zRx, s.rTx] +returnforwrite(s::VTEMsoundingData) = getfield.(Ref(s), s.writefields) function getndata(S::VTEMsoundingData) getndata(S.data)[1] @@ -132,7 +132,8 @@ function VTEMsoundingData(;rRx=nothing, zRx=nothing, zTx=12., 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) + linenum=nothing, fid=nothing, writefields=[:X, :Y, :Z, :fid, + :linenum, :zTx, :zRx, :rTx]) @assert rTx > 0 @assert zTx < 0 # my coordinate system z down isnothing(rRx) && (rRx = 0.) @@ -143,7 +144,7 @@ function VTEMsoundingData(;rRx=nothing, zRx=nothing, zTx=12., @assert all((noise .>0) .| isnan.(noise)) @assert length(data) == length(noise) VTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, zTx, zRx, rTx, - lowpassfcs, times, ramp, noise, data) + lowpassfcs, times, ramp, noise, data, writefields) end function read_survey_files(; From f7669bfb9dd3a26c133b924a9754653c9c29b177 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Wed, 2 Oct 2024 17:46:16 +1000 Subject: [PATCH 06/22] WIP:need to check if DFN writing works for gradinv --- src/AEMnoNuisanceGradientInversionTools.jl | 18 ++++-------- src/AEMwithNuisanceGradientInversionTools.jl | 29 +++++++++----------- src/CommonToAll.jl | 2 +- src/SkyTEM1DInversion.jl | 10 +++---- src/TEMPEST1DInversion.jl | 11 ++++---- 5 files changed, 30 insertions(+), 40 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index 273398d..659a0cb 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -20,29 +20,23 @@ function compress(soundings, zall; prefix="", rmfile=true, isfirstparalleliterat isfirstparalleliteration && throw(AssertionError("Zipped file "*fout*" exists, will not overwrite!")) iomode = "a" end - io = open(fout, iomode) 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.(soundings[1].writefields)..., "zcenter", "log10_cond", "ϕd_err"] + sfmt = fill("%15.3f", nelinonerow) + ϕd > 1e4 && (ϕd = 1e4 ) + writeasegdat(elinonerow, sfmt, fout[1:end-4], iomode) rmfile && rm(fname) if isfirstparalleliteration && i == 1 - elinonerow = [returnforwrite(s), vec(zall), σgrid, ϕd] - nelinonerow = length(elinonerow) - writenames = [[f for f in fieldnames(typeof(s))], "zcenter", "log10_cond", "ϕd_err"] - sfmt = fill("%15.3f", nelinonerow) channel_names = [writenames, fill("", nelinonerow), writenames] writeasegdfnfromonerow(elinonerow, channel_names, sfmt, fout[1:end-4]) end end - close(io) end # plot the convergence and the result diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl index 56dacf2..84b759e 100644 --- a/src/AEMwithNuisanceGradientInversionTools.jl +++ b/src/AEMwithNuisanceGradientInversionTools.jl @@ -22,31 +22,28 @@ function compress(soundings, zall, nnu; prefix="", rmfile=true, isfirstparalleli isfirstparalleliteration && throw(AssertionError("Zipped file "*fout*" exists, will not overwrite!")) iomode = "a" end - io = open(fout, iomode) 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) + # for el in [returnforwrite(s)...; vec(zall); σgrid; nu; ϕd] + # msg = @sprintf("%.4f\t", el) + # write(io, msg) + # end + writenames = [string.(soundings[1].writefields), "zcenter", "log10_cond", "nuisance_param_in_order", "ϕd_err"] + sfmt = fill("%15.3f", nelinonerow) + ϕd > 1e4 && (ϕd = 1e4 ) + writeasegdat(elinonerow, sfmt, fout[1:end-4], iomode) rmfile && rm(fname) - if i == 1 - elinonerow = [returnforwrite(s), vec(zall), σgrid, nu, ϕd] - nelinonerow = length(elinonerow) - writenames = [[f for f in fieldnames(typeof(a))], "zcenter", "log10_cond", "nuisance_param", "ϕd_err"] - sfmt = fill("%15.3f", nelinonerow) - channel_names = [writenames, fill("", nelinonerow, writenames)] - writeasegdfnforgradinv(elinonerow, fname[1:end-4]) - writeasegdfnfromonerow(elinonerow, channel_names, sfmt, fname[1:end-4]) + if isfirstparalleliteration && i == 1 + channel_names = [writenames, fill("", nelinonerow), writenames] + writeasegdfnfromonerow(elinonerow, channel_names, sfmt, fout[1:end-4]) end end - close(io) end # plot the convergence and the result diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 97801e7..652d412 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -19,7 +19,7 @@ export trimxft, assembleTat1, gettargtemps, checkns, getchi2forall, nicenup, plo readcols, colstovtk, findclosestidxincolfile, zcentertoboundary, zboundarytocenter, writeijkfromsounding, nanmean, infmean, nanstd, infstd, infnanmean, infnanstd, kde_sj, plotmanygrids, readwell, getlidarheight, plotblockedwellonimages, getdeterministicoutputs, writeasegdfnfromonerow, - getprobabilisticoutputs,readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getRandgridr, getallxyinr, getXYlast, + writeasegdat, getprobabilisticoutputs, readfzipped, readxyzrhoϕ, writevtkxmlforcurtain, getRandgridr, getallxyinr, getXYlast, getprobabilisticlinesfromdirectory, readxyzrhoϕnu, plotgausshist, findMAE, getclosestidx, getpostatXY, getpostatwell # Kernel Density stuff diff --git a/src/SkyTEM1DInversion.jl b/src/SkyTEM1DInversion.jl index d61a4ac..406cbe6 100644 --- a/src/SkyTEM1DInversion.jl +++ b/src/SkyTEM1DInversion.jl @@ -144,11 +144,10 @@ mutable struct SkyTEMsoundingData <: Sounding LM_data :: Array{Float64, 1} HM_data :: Array{Float64, 1} forceML + writefields end -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] +returnforwrite(s::SkyTEMsoundingData) = getfield.(Ref(s), s.writefields) function getndata(S::SkyTEMsoundingData) getndata(S.LM_data)[1] + getndata(S.HM_data)[1] @@ -160,7 +159,8 @@ function SkyTEMsoundingData(;rRx=-12., zRxLM=12., zTxLM=12., HM_times=[1., 2.], HM_ramp=[1 2; 3 4], LM_noise=[1.], HM_noise=[1.], LM_data=[1.], HM_data=[1.], sounding_string="sounding", X=nothing, Y=nothing, Z=nothing, - linenum=nothing, fid=nothing, forceML=false) + linenum=nothing, fid=nothing, forceML=false, writefields=[:X, :Y, :Z, :fid, + :linenum, :rRx, :zRxLM, :zTxLM, :zRxHM, :zTxHM, :rTx]) @assert rRx > 0 && rTx > 0 @assert zRxLM <0 && zTxLM <0 @assert zRxHM <0 && zTxHM <0 @@ -175,7 +175,7 @@ function SkyTEMsoundingData(;rRx=-12., zRxLM=12., zTxLM=12., @assert length(HM_data) == length(HM_noise) SkyTEMsoundingData(sounding_string, X, Y, Z, fid, linenum, rRx, zRxLM, zTxLM, zRxHM, zTxHM, rTx, lowpassfcs, LM_times, LM_ramp, HM_times, HM_ramp, LM_noise, HM_noise, - LM_data, HM_data, forceML) + LM_data, HM_data, forceML, writefields) end function read_survey_files(; diff --git a/src/TEMPEST1DInversion.jl b/src/TEMPEST1DInversion.jl index 786a106..127d2a1 100644 --- a/src/TEMPEST1DInversion.jl +++ b/src/TEMPEST1DInversion.jl @@ -746,7 +746,8 @@ mutable struct TempestSoundingData <: Sounding Hx_data :: Array{Float64, 1} Hz_data :: Array{Float64, 1} forceML :: Bool - peakcurrent + peakcurrent + writefields end function getnufromsounding(t::TempestSoundingData) @@ -764,10 +765,7 @@ function getnufromsounding(t::TempestSoundingData) t.roll_tx, t.pitch_tx, t.yaw_tx] end - -returnforwrite(s::TempestSoundingData) = [s.X, s.Y, s.Z, s.fid, - s.linenum, getnufromsounding(s)...] - +returnforwrite(s::TEMPESTsoundingData) = getfield.(Ref(s), s.writefields) function getndata(s::TempestSoundingData, vectorsum) n = getndata(s.Hx_data)[1] + getndata(s.Hz_data)[1] @@ -933,7 +931,8 @@ function read_survey_files(; z_tx[is], times, ramp, σ_Hx[is,:], σ_Hz[is,:], dHx, dHz, forceML, - peakcurrent + peakcurrent, [:X, :Y, :Z, :fid, :linenum, :z_tx, :z_rx, :x_rx, :y_rx, :roll_rx, :pitch_rx, :yaw_rx, + :roll_tx, :pitch_tx, :yaw_tx] ) fracnew = round(Int, is/nsoundings*100) if (fracnew-fracdone)>10 From 51b12512d65325dc774911d6a1dfda61f44fb04d Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Wed, 2 Oct 2024 18:09:25 +1000 Subject: [PATCH 07/22] added in writing DFN->HDR --- src/AEMnoNuisanceGradientInversionTools.jl | 1 + src/AEMwithNuisanceGradientInversionTools.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index 659a0cb..4c5678d 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -35,6 +35,7 @@ function compress(soundings, zall; prefix="", rmfile=true, isfirstparalleliterat 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 end diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl index 84b759e..c2a016a 100644 --- a/src/AEMwithNuisanceGradientInversionTools.jl +++ b/src/AEMwithNuisanceGradientInversionTools.jl @@ -42,6 +42,7 @@ function compress(soundings, zall, nnu; prefix="", rmfile=true, isfirstparalleli 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 end From 33fea83a4b860dcb8d2852068b6b554973fbe8a5 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 3 Oct 2024 10:42:25 +1000 Subject: [PATCH 08/22] case for TEMPEST --- src/TEMPEST1DInversion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TEMPEST1DInversion.jl b/src/TEMPEST1DInversion.jl index 127d2a1..b7e6ca5 100644 --- a/src/TEMPEST1DInversion.jl +++ b/src/TEMPEST1DInversion.jl @@ -765,7 +765,7 @@ function getnufromsounding(t::TempestSoundingData) t.roll_tx, t.pitch_tx, t.yaw_tx] end -returnforwrite(s::TEMPESTsoundingData) = getfield.(Ref(s), s.writefields) +returnforwrite(s::TempestSoundingData) = getfield.(Ref(s), s.writefields) function getndata(s::TempestSoundingData, vectorsum) n = getndata(s.Hx_data)[1] + getndata(s.Hz_data)[1] From d9b181fbeb780ec6b1f47003b64266980c8bc04e Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 3 Oct 2024 11:12:38 +1000 Subject: [PATCH 09/22] WIP:writing ASEG-GDF --- src/AEMwithNuisanceGradientInversionTools.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl index c2a016a..20fd0cb 100644 --- a/src/AEMwithNuisanceGradientInversionTools.jl +++ b/src/AEMwithNuisanceGradientInversionTools.jl @@ -28,13 +28,9 @@ function compress(soundings, zall, nnu; prefix="", rmfile=true, isfirstparalleli ϕd = A[end,2] σgrid = vec(A[end,3:end-nnu]) nu = vec(A[end,end-nnu+1:end]) - elinonerow = [returnforwrite(s), vec(zall), σgrid, nu, ϕd] + elinonerow = [returnforwrite(s)..., vec(zall), σgrid, nu, ϕd] nelinonerow = length(elinonerow) - # for el in [returnforwrite(s)...; vec(zall); σgrid; nu; ϕd] - # msg = @sprintf("%.4f\t", el) - # write(io, msg) - # end - writenames = [string.(soundings[1].writefields), "zcenter", "log10_cond", "nuisance_param_in_order", "ϕd_err"] + writenames = [string.(soundings[1].writefields)..., "zcenter", "log10_cond", "nuisance_param_in_order", "ϕd_err"] sfmt = fill("%15.3f", nelinonerow) ϕd > 1e4 && (ϕd = 1e4 ) writeasegdat(elinonerow, sfmt, fout[1:end-4], iomode) From dc4927dae4b8a4d69be6a1d9034e6374ed455a41 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 3 Oct 2024 11:58:45 +1000 Subject: [PATCH 10/22] WIP:ASEG-GDF stop stomping on existing file --- src/AEMnoNuisanceGradientInversionTools.jl | 9 +++------ src/AEMwithNuisanceGradientInversionTools.jl | 9 +++------ src/GradientInversion.jl | 6 ++++-- 3 files changed, 10 insertions(+), 14 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index 4c5678d..a650395 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -15,11 +15,7 @@ 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 + 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) @@ -27,9 +23,10 @@ function compress(soundings, zall; prefix="", rmfile=true, isfirstparalleliterat σgrid = vec(A[end,3:end]) elinonerow = [returnforwrite(s)..., vec(zall), σgrid, ϕd] nelinonerow = length(elinonerow) - writenames = [string.(soundings[1].writefields)..., "zcenter", "log10_cond", "ϕd_err"] + writenames = [string.(s.writefields)..., "zcenter", "log10_cond", "ϕd_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 diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl index 20fd0cb..5caa297 100644 --- a/src/AEMwithNuisanceGradientInversionTools.jl +++ b/src/AEMwithNuisanceGradientInversionTools.jl @@ -17,11 +17,7 @@ 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 + 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) @@ -30,9 +26,10 @@ function compress(soundings, zall, nnu; prefix="", rmfile=true, isfirstparalleli nu = vec(A[end,end-nnu+1:end]) elinonerow = [returnforwrite(s)..., vec(zall), σgrid, nu, ϕd] nelinonerow = length(elinonerow) - writenames = [string.(soundings[1].writefields)..., "zcenter", "log10_cond", "nuisance_param_in_order", "ϕd_err"] + writenames = [string.(s.writefields)..., "zcenter", "log10_cond", "nuisance_param_in_order", "ϕd_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 diff --git a/src/GradientInversion.jl b/src/GradientInversion.jl index bd8a466..892b3c2 100644 --- a/src/GradientInversion.jl +++ b/src/GradientInversion.jl @@ -231,6 +231,7 @@ function gradientinv( m::AbstractVector, knownvalue=0.7, firstvalue=:last, κ = GP.Mat52(), + verbose = false, breakonknown=false, minimprovfrac = nothing, minimprovkickinstep = round(Int, nstepsmax/2), @@ -264,7 +265,7 @@ function gradientinv( m::AbstractVector, lo, hi, λ²min, λ²max, β², ntries, knownvalue=knownvalue, regularizeupdate = regularizeupdate) end prefix = isempty(fname) ? fname : fname*" : " - @info prefix*"iteration: $istep χ²: $(χ²[istep][idx]) target: $target" + verbose && @info prefix*"iteration: $istep χ²: $(χ²[istep][idx]) target: $target" m = mnew[istep][idx] oidx[istep] = idx isa(io, Nothing) || write_history(io, [istep; χ²[istep][idx]/target₀; vec(m)]) @@ -314,6 +315,7 @@ function gradientinv( m::AbstractVector, debuglevel = 0, usebox = true, boxiters = 2, + verbose = false, minimprovfrac = nothing, minimprovkickinstep = round(Int, nstepsmax/2), fname="") @@ -354,7 +356,7 @@ function gradientinv( m::AbstractVector, idx, foundroot = occamstep(m, m0, Δm, mnew[istep], χ²[istep], λ²[istep], F, R, target, lo, hi, λ²min, λ²max, β², ntries, knownvalue=knownvalue, regularizeupdate = regularizeupdate) prefix = isempty(fname) ? fname : fname*" : " - @info prefix*"iteration: $istep χ²: $(χ²[istep][idx]) target: $target" + verbose && @info prefix*"iteration: $istep χ²: $(χ²[istep][idx]) target: $target" m = mnew[istep][idx] oidx[istep] = idx isa(io, Nothing) || write_history(io, [istep; χ²[istep][idx]/target₀; vec(m); nu]) From e8d555036b988675461ba97d94e48702b7bbef9e Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 3 Oct 2024 12:10:31 +1000 Subject: [PATCH 11/22] even less verbose grad by default --- src/GradientInversion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GradientInversion.jl b/src/GradientInversion.jl index 892b3c2..b68073a 100644 --- a/src/GradientInversion.jl +++ b/src/GradientInversion.jl @@ -284,7 +284,7 @@ function gradientinv( m::AbstractVector, istep > nstepsmax && break end if !isnothing(io) - @info "Finished "*fname + verbose && @info "Finished "*fname close(io) nothing return @@ -375,7 +375,7 @@ function gradientinv( m::AbstractVector, istep > nstepsmax && break end if !isnothing(io) - @info "Finished "*fname + verbose && @info "Finished "*fname close(io) nothing return From 9512f2fc41671dea7d767bde17f9320ce280dc04 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 3 Oct 2024 12:21:17 +1000 Subject: [PATCH 12/22] verbosity on for interactive, not loopacross --- src/AEMnoNuisanceGradientInversionTools.jl | 3 ++- src/AEMwithNuisanceGradientInversionTools.jl | 6 +++++- src/GradientInversion.jl | 4 ++-- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index a650395..835c713 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -237,6 +237,7 @@ 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 @@ -279,7 +280,7 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0; breakonknown = breakonknown , dobo = dobo , fname = fname , - minimprovfrac, + minimprovfrac, verbose, minimprovkickinstep) diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl index 5caa297..cf0f3a5 100644 --- a/src/AEMwithNuisanceGradientInversionTools.jl +++ b/src/AEMwithNuisanceGradientInversionTools.jl @@ -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, @@ -273,7 +276,8 @@ 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 diff --git a/src/GradientInversion.jl b/src/GradientInversion.jl index b68073a..091bc7b 100644 --- a/src/GradientInversion.jl +++ b/src/GradientInversion.jl @@ -231,7 +231,7 @@ function gradientinv( m::AbstractVector, knownvalue=0.7, firstvalue=:last, κ = GP.Mat52(), - verbose = false, + verbose = true, breakonknown=false, minimprovfrac = nothing, minimprovkickinstep = round(Int, nstepsmax/2), @@ -315,7 +315,7 @@ function gradientinv( m::AbstractVector, debuglevel = 0, usebox = true, boxiters = 2, - verbose = false, + verbose = true, minimprovfrac = nothing, minimprovkickinstep = round(Int, nstepsmax/2), fname="") From db76934e2cee235ce9cb3a9eb2d75eed268f97a7 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 3 Oct 2024 12:26:10 +1000 Subject: [PATCH 13/22] first parallel iter start time --- src/AEMnoNuisanceGradientInversionTools.jl | 1 + src/AEMwithNuisanceGradientInversionTools.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index 835c713..b12ca44 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -245,6 +245,7 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0; nparallelsoundings = nworkers() nsoundings = length(soundings) zall = zboundarytocenter(aem_in.z[aem_in.nfixed+1:end]) # needed for sounding compression in write + @info "starting sequential parralel iterations at $(Dates.now())" for iter = 1:nsequentialiters if iter Date: Thu, 3 Oct 2024 12:40:22 +1000 Subject: [PATCH 14/22] first parallel iter start time --- src/AEMnoNuisanceGradientInversionTools.jl | 2 +- src/AEMwithNuisanceGradientInversionTools.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index b12ca44..e16db40 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -245,7 +245,7 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0; nparallelsoundings = nworkers() nsoundings = length(soundings) zall = zboundarytocenter(aem_in.z[aem_in.nfixed+1:end]) # needed for sounding compression in write - @info "starting sequential parralel iterations at $(Dates.now())" + @info "starting sequential parallel iterations at $(Dates.now())" for iter = 1:nsequentialiters if iter Date: Thu, 3 Oct 2024 13:25:24 +1000 Subject: [PATCH 15/22] timestamps and dt cleanups --- src/AEMnoNuisanceGradientInversionTools.jl | 5 ++++- src/AEMnoNuisanceMcMCInversionTools.jl | 2 +- src/AEMwithNuisanceGradientInversionTools.jl | 5 ++++- src/AEMwithNuisanceMcMCInversionTools.jl | 2 +- 4 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index e16db40..c3463d9 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -245,6 +245,7 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0; nparallelsoundings = nworkers() nsoundings = length(soundings) zall = zboundarytocenter(aem_in.z[aem_in.nfixed+1:end]) # needed for sounding compression in write + @info "starting sequential parallel iterations at $(Dates.now())" for iter = 1:nsequentialiters if iter Date: Thu, 3 Oct 2024 15:33:31 +1000 Subject: [PATCH 16/22] phi no symbol in dfn --- src/AEMnoNuisanceGradientInversionTools.jl | 2 +- src/AEMwithNuisanceGradientInversionTools.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index c3463d9..305fb89 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/src/AEMnoNuisanceGradientInversionTools.jl @@ -23,7 +23,7 @@ function compress(soundings, zall; prefix="", rmfile=true, isfirstparalleliterat σgrid = vec(A[end,3:end]) elinonerow = [returnforwrite(s)..., vec(zall), σgrid, ϕd] nelinonerow = length(elinonerow) - writenames = [string.(s.writefields)..., "zcenter", "log10_cond", "ϕd_err"] + writenames = [string.(s.writefields)..., "zcenter", "log10_cond", "phid_err"] sfmt = fill("%15.3f", nelinonerow) ϕd > 1e4 && (ϕd = 1e4 ) iomode = (isfirstparalleliteration && i==1) ? "w" : "a" diff --git a/src/AEMwithNuisanceGradientInversionTools.jl b/src/AEMwithNuisanceGradientInversionTools.jl index 2f22668..caca7aa 100644 --- a/src/AEMwithNuisanceGradientInversionTools.jl +++ b/src/AEMwithNuisanceGradientInversionTools.jl @@ -26,7 +26,7 @@ function compress(soundings, zall, nnu; prefix="", rmfile=true, isfirstparalleli nu = vec(A[end,end-nnu+1:end]) elinonerow = [returnforwrite(s)..., vec(zall), σgrid, nu, ϕd] nelinonerow = length(elinonerow) - writenames = [string.(s.writefields)..., "zcenter", "log10_cond", "nuisance_param_in_order", "ϕd_err"] + 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" From ced067fe81c283b95cc25175387312812f9b0f49 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 3 Oct 2024 20:07:34 +1000 Subject: [PATCH 17/22] 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. --- src/AEM_VMD_HMD.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/AEM_VMD_HMD.jl b/src/AEM_VMD_HMD.jl index 37e5b46..e0e1476 100644 --- a/src/AEM_VMD_HMD.jl +++ b/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 From c50d63e5cb0df27820819fce59a9942b6446c856 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Sun, 6 Oct 2024 23:40:42 +1100 Subject: [PATCH 18/22] silent logging works for McMC --- Project.toml | 1 + examples/VTEM/DarlingParoo/McMC/04_plot.jl | 2 +- .../tempest/Menindee/McMC/01_read_data.jl | 2 +- .../tempest/Menindee/McMC/02_set_options.jl | 2 +- examples/tempest/Menindee/McMC/04_plot.jl | 2 +- src/AEMnoNuisanceMcMCInversionTools.jl | 7 +-- src/AEMwithNuisanceMcMCInversionTools.jl | 7 +-- src/CommonToAll.jl | 1 + src/MCMC_Driver.jl | 49 +++++++++++-------- src/SoundingDistributor.jl | 22 ++++++++- 10 files changed, 62 insertions(+), 33 deletions(-) diff --git a/Project.toml b/Project.toml index a355622..59b1684 100644 --- a/Project.toml +++ b/Project.toml @@ -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" diff --git a/examples/VTEM/DarlingParoo/McMC/04_plot.jl b/examples/VTEM/DarlingParoo/McMC/04_plot.jl index 48241d7..4c849c9 100755 --- a/examples/VTEM/DarlingParoo/McMC/04_plot.jl +++ b/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 diff --git a/examples/tempest/Menindee/McMC/01_read_data.jl b/examples/tempest/Menindee/McMC/01_read_data.jl index ac5c22b..011fc97 100644 --- a/examples/tempest/Menindee/McMC/01_read_data.jl +++ b/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) \ No newline at end of file diff --git a/examples/tempest/Menindee/McMC/02_set_options.jl b/examples/tempest/Menindee/McMC/02_set_options.jl index 14ae8c6..49f082e 100644 --- a/examples/tempest/Menindee/McMC/02_set_options.jl +++ b/examples/tempest/Menindee/McMC/02_set_options.jl @@ -1,5 +1,5 @@ ## same for all soundingss -nsamples = 201 +nsamples = 2001 nchainspersounding = 3 ppn = 4 @assert mod(ppn, nchainspersounding+1) == 0 diff --git a/examples/tempest/Menindee/McMC/04_plot.jl b/examples/tempest/Menindee/McMC/04_plot.jl index d804bbf..1fb30c8 100644 --- a/examples/tempest/Menindee/McMC/04_plot.jl +++ b/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, diff --git a/src/AEMnoNuisanceMcMCInversionTools.jl b/src/AEMnoNuisanceMcMCInversionTools.jl index be46711..c1bcc47 100644 --- a/src/AEMnoNuisanceMcMCInversionTools.jl +++ b/src/AEMnoNuisanceMcMCInversionTools.jl @@ -337,10 +337,10 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_ nsoundings = length(soundings) nsequentialiters, nparallelsoundings = splittasks(soundings; nchainspersounding, ppn) - @info "starting sequential parallel iterations at $(Dates.now())" + 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) @@ -354,7 +354,8 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_ end # @sync dt = time() - t2 #seconds - @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 diff --git a/src/AEMwithNuisanceMcMCInversionTools.jl b/src/AEMwithNuisanceMcMCInversionTools.jl index fb62716..b2a6899 100644 --- a/src/AEMwithNuisanceMcMCInversionTools.jl +++ b/src/AEMwithNuisanceMcMCInversionTools.jl @@ -173,10 +173,10 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_ nsoundings = length(soundings) nsequentialiters, nparallelsoundings = splittasks(soundings; nchainspersounding, ppn) - @info "starting sequential parallel iterations at $(Dates.now())" + 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) @@ -191,7 +191,8 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in::Operator1D, opt_ end # @sync dt = time() - t2 #seconds - @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 diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index 652d412..9823878 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1673,6 +1673,7 @@ function plotsummarygrids1(soundings, meangrid, phgrid, plgrid, pmgrid, gridx, g idx_split, thereisrmn = getRsplits(gridx, Rmax) nimages = length(idx_split) dr = gridx[2] - gridx[1] + dr >= Rmax && error("dr:$dr >= Rmax:$Rmax, decrease dr") if iszero(idx_split[1]) && thereisrmn # Rmax is larger than section nx = length(range(gridx[1], Rmax, step=dr)) elseif !iszero(idx_split[1]) && !thereisrmn # Rmax is exactly the section length diff --git a/src/MCMC_Driver.jl b/src/MCMC_Driver.jl index 750ec25..21043b5 100755 --- a/src/MCMC_Driver.jl +++ b/src/MCMC_Driver.jl @@ -1,5 +1,5 @@ using Distributed, DistributedArrays, - PyPlot, LinearAlgebra, Dates + PyPlot, LinearAlgebra, Dates, Logging import .AbstractOperator.Operator import .AbstractOperator.Sounding @@ -411,7 +411,7 @@ function domcmciters(iterlast, nsamples, chains, mns::DArray{ModelNonstat}, m::D current_misfit, F, chain.T, isample, wp, chain_idx, chain.master_pid) end - t, tlong, doquit = disptime(isample, t, tlong, iterlast, nsamples, nominaltime) + t, tlong, doquit = disptime(isample, t, tlong, nsamples, nominaltime) doquit && break end end @@ -525,7 +525,7 @@ function domcmciters(iterlast, nsamples, chains, mns::DArray{ModelNonstat}, m::D current_misfit, F, chain.T, isample, wp, chain_idx, chain.master_pid) end - t, tlong, doquit = disptime(isample, t, tlong, iterlast, nsamples, nominaltime) + t, tlong, doquit = disptime(isample, t, tlong, nsamples, nominaltime) doquit && break end end @@ -637,7 +637,7 @@ function domcmciters(iterlast, nsamples, chains, m::DArray{ModelStat}, mn::DArra current_misfit, F, chain.T, isample, wp, chain_idx, chain.master_pid) end - t, tlong, doquit = disptime(isample, t, tlong, iterlast, nsamples, nominaltime) + t, tlong, doquit = disptime(isample, t, tlong, nsamples, nominaltime) doquit && break end end @@ -729,7 +729,7 @@ function domcmciters(iterlast, nsamples, chains, m::DArray{ModelStat}, current_misfit, F, chain.T, isample, wp, chain_idx, chain.master_pid) end - t, tlong, doquit = disptime(isample, t, tlong, iterlast, nsamples, nominaltime) + t, tlong, doquit = disptime(isample, t, tlong, nsamples, nominaltime) doquit && break end end @@ -783,25 +783,32 @@ function swap_temps(chains::Array{Chain, 1}) end end -function disptime(isample, t, tlong, iterlast, nsamples, nominaltime) +function disptime(isample, t, tlong, nsamples, nominaltime) + iomode = isample == 1 ? "w" : "a" + io = open("$(myid()).log", iomode) + logger = ConsoleLogger(io) doquit = false - windowtime = 1000 - nw = 5 - if mod(isample-1, windowtime) == 0 - dt = time() - t #seconds - t = time() - @info("on pid $(myid()) **$(@sprintf("%.2f", dt))**sec** $isample out of $(nsamples)") - end - if !isnothing(nominaltime) - if mod(isample-1, nw*windowtime) == 0 - dt = time() - tlong #seconds - tlong = time() - if (dt/nw) > nominaltime - doquit = true - @info("QUIT: PID: $(myid()) TIME:**$(@sprintf("%.2f", dt/nw))**SEC AVG** SAMPLE: $isample") + with_logger(logger) do + windowtime = 1000 + nw = 5 + if mod(isample-1, windowtime) == 0 + dt = time() - t #seconds + t = time() + @info("on pid $(myid()) **$(@sprintf("%.2f", dt))**sec** $isample out of $(nsamples)") + end + if !isnothing(nominaltime) + if mod(isample-1, nw*windowtime) == 0 + dt = time() - tlong #seconds + tlong = time() + if (dt/nw) > nominaltime + doquit = true + @info("QUIT: PID: $(myid()) TIME:**$(@sprintf("%.2f", dt/nw))**SEC AVG** SAMPLE: $isample") + end end end - end + end # logger + flush(io) + close(io) t, tlong, doquit end diff --git a/src/SoundingDistributor.jl b/src/SoundingDistributor.jl index f487330..8867f11 100644 --- a/src/SoundingDistributor.jl +++ b/src/SoundingDistributor.jl @@ -1,6 +1,6 @@ module SoundingDistributor using Distributed -export splittasks, getss, getpids +export splittasks, getss, getpids, writetogloballog, catlocallogs function splittasks(soundings::AbstractVector; nchainspersounding=nothing, ppn=nothing) nsoundings = length(soundings) @@ -16,7 +16,7 @@ function splittasks(;nsoundings=nothing, ncores=nothing, nchainspersounding=noth subtractone = 1 nparallelsoundings = Int((ncores+1)/(nchainspersounding+1)) - subtractone nsequentialiters = ceil(Int, nsoundings/nparallelsoundings) - @info "will require $nsequentialiters iterations of $nparallelsoundings soundings in one iteration" + writetogloballog( "will require $nsequentialiters iterations of $nparallelsoundings soundings in one iteration", iomode="w") nsequentialiters, nparallelsoundings end @@ -33,4 +33,22 @@ function getpids(i, nchainspersounding) pids = reverse((nprocs()+2) .- (i*(nchainspersounding+1)+1:(i+1)*(nchainspersounding+1))) end +# Global logging function +function writetogloballog(str::String; fname="global.log", iomode="a", newline=true) + io = open(fname, iomode) + nwchar = newline ? "\n" : "" + write(io, str*nwchar) + flush(io) + close(io) +end + +function catlocallogs(nparallelsoundings, nchainspersounding) + managerpids = [p[1] for p in getpids.(1:nparallelsoundings, nchainspersounding)] + map(managerpids) do p + locallogfile = "$p.log" + writetogloballog(read(locallogfile, String), newline=false) + rm(locallogfile) + end +end + end From 9df7d17f1b8beef573c0114f180bf3b5a8833ab0 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Mon, 7 Oct 2024 00:11:09 +1100 Subject: [PATCH 19/22] seems to work with silent logging --- src/AEMnoNuisanceGradientInversionTools.jl | 17 +++++++---------- src/AEMwithNuisanceGradientInversionTools.jl | 17 +++++++---------- src/SoundingDistributor.jl | 10 +++++++++- 3 files changed, 23 insertions(+), 21 deletions(-) diff --git a/src/AEMnoNuisanceGradientInversionTools.jl b/src/AEMnoNuisanceGradientInversionTools.jl index 305fb89..4c8f894 100644 --- a/src/AEMnoNuisanceGradientInversionTools.jl +++ b/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 @@ -246,14 +246,11 @@ function loopacrossAEMsoundings(soundings::Array{S, 1}, aem_in, σstart, σ0; nsoundings = length(soundings) zall = zboundarytocenter(aem_in.z[aem_in.nfixed+1:end]) # needed for sounding compression in write - @info "starting sequential parallel iterations at $(Dates.now())" + writetogloballog("will require $nsequentialiters iterations of $nparallelsoundings soundings", iomode="w") + writetogloballog("starting sequential parallel iterations at $(Dates.now())") for iter = 1:nsequentialiters - if iter Date: Mon, 7 Oct 2024 01:19:50 +1100 Subject: [PATCH 20/22] only cat logs that exist --- src/SoundingDistributor.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/SoundingDistributor.jl b/src/SoundingDistributor.jl index 267b30b..271c15e 100644 --- a/src/SoundingDistributor.jl +++ b/src/SoundingDistributor.jl @@ -54,8 +54,10 @@ function catlocallogs(nparallelsoundings, nchainspersounding) managerpids = [p[1] for p in getpids.(1:nparallelsoundings, nchainspersounding)] map(managerpids) do p locallogfile = "$p.log" - writetogloballog(read(locallogfile, String), newline=false) - rm(locallogfile) + if isfile(locallogfile) + writetogloballog(read(locallogfile, String), newline=false) + rm(locallogfile) + end end end From 9fdca760609543d0408a317195d2eb1ba7ffacd5 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Mon, 7 Oct 2024 15:56:49 +1100 Subject: [PATCH 21/22] DArray changes --- Project.toml | 2 +- examples/tempest/Menindee/McMC/02_set_options.jl | 2 +- src/MCMC_Driver.jl | 4 ++++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 59b1684..5108d08 100644 --- a/Project.toml +++ b/Project.toml @@ -51,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" diff --git a/examples/tempest/Menindee/McMC/02_set_options.jl b/examples/tempest/Menindee/McMC/02_set_options.jl index 49f082e..7895182 100644 --- a/examples/tempest/Menindee/McMC/02_set_options.jl +++ b/examples/tempest/Menindee/McMC/02_set_options.jl @@ -1,5 +1,5 @@ ## same for all soundingss -nsamples = 2001 +nsamples = 3001 nchainspersounding = 3 ppn = 4 @assert mod(ppn, nchainspersounding+1) == 0 diff --git a/src/MCMC_Driver.jl b/src/MCMC_Driver.jl index 21043b5..2e5cdbf 100755 --- a/src/MCMC_Driver.jl +++ b/src/MCMC_Driver.jl @@ -437,6 +437,7 @@ function main(opt_in ::OptionsStat, optn, statns, stat, statn, current_misfit, F, wpns, wp, wpn, nominaltime) close_history.([wp, wpns, wpn]) + d_closeall() nothing end @@ -550,6 +551,7 @@ function main(opt_in ::OptionsStat, statns, stat, current_misfit, F, wpns, wp, nominaltime) close_history.([wp, wpns]) + d_closeall() nothing end @@ -662,6 +664,7 @@ function main(opt_in ::OptionsStat, optn, stat, statn, current_misfit, F, wp, wpn, nominaltime) close_history.([wp, wpn]) + d_closeall() nothing end @@ -753,6 +756,7 @@ function main(opt_in ::OptionsStat, stat, current_misfit, F, wp, nominaltime) close_history(wp) + d_closeall() nothing end From 9954a092e68a29dce71cc4ee07c98aa0eeb63123 Mon Sep 17 00:00:00 2001 From: a2ray <52647259+a2ray@users.noreply.github.com> Date: Mon, 7 Oct 2024 21:10:27 +1100 Subject: [PATCH 22/22] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5108d08..8ad52a1 100644 --- a/Project.toml +++ b/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.7" +version = "0.4.8" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"