From 27519da05aa0e38d423a2dbdb86f7bb311628f62 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 12 Jan 2024 15:35:45 +1100 Subject: [PATCH 1/5] removed extraneous exports from AEM nuisance McMC helper --- .../gradientbased/03_run_inversion_nuisance.jl | 2 +- src/AEMnoNuisanceMcMCInversionTools.jl | 3 +++ src/AEMwithNuisanceMcMCInversionTools.jl | 15 ++++++++++----- src/AbstractOperator.jl | 4 ++++ src/CSEM1DInversion.jl | 2 +- 5 files changed, 19 insertions(+), 7 deletions(-) diff --git a/examples/tempest/synth/gradientbased/03_run_inversion_nuisance.jl b/examples/tempest/synth/gradientbased/03_run_inversion_nuisance.jl index 7042b6db..45390693 100644 --- a/examples/tempest/synth/gradientbased/03_run_inversion_nuisance.jl +++ b/examples/tempest/synth/gradientbased/03_run_inversion_nuisance.jl @@ -36,4 +36,4 @@ else @info transD_GP.TEMPEST1DInversion.get_misfit(-mtest, [zRx, x_rx, rx_pitch], tempest, nubounds) end @info transD_GP.TEMPEST1DInversion.get_misfit(-mtest, nu[end], tempest, nubounds) -transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest, mtest, nutest) \ No newline at end of file +transD_GP.plotmodelfield!(tempest, mtest, nutest) \ No newline at end of file diff --git a/src/AEMnoNuisanceMcMCInversionTools.jl b/src/AEMnoNuisanceMcMCInversionTools.jl index b73905f1..98385686 100644 --- a/src/AEMnoNuisanceMcMCInversionTools.jl +++ b/src/AEMnoNuisanceMcMCInversionTools.jl @@ -5,6 +5,9 @@ import ..AbstractOperator.loopacrossAEMsoundings import ..AbstractOperator.plotmodelfield! import ..AbstractOperator.getndata import ..Options, ..OptionsStat +import ..AbstractOperator.summaryAEMimages +import ..AbstractOperator.plotindividualAEMsoundings + export makeAEMoperatorandoptions, loopacrossAEMsoundings, summaryAEMimages, plotindividualAEMsoundings import ..main # McMC function import ..DEBUGLEVEL_TDGP diff --git a/src/AEMwithNuisanceMcMCInversionTools.jl b/src/AEMwithNuisanceMcMCInversionTools.jl index 278369eb..4d83c551 100644 --- a/src/AEMwithNuisanceMcMCInversionTools.jl +++ b/src/AEMwithNuisanceMcMCInversionTools.jl @@ -7,9 +7,12 @@ import ..AbstractOperator.getndata import ..AbstractOperator.makebounds import ..AbstractOperator.getoptnfromexisting import ..AbstractOperator.getnufromsounding +import ..AbstractOperator.summaryAEMimages +import ..AbstractOperator.plotindividualAEMsoundings import ..Options, ..OptionsStat, ..OptionsNuisance -export makeAEMoperatorandnuisanceoptions, loopacrossAEMsoundings, summaryAEMwithnuisanceimages, plotindividualAEMsoundingswithnuisance +export makeAEMoperatorandnuisanceoptions + import ..main # McMC function using ..SoundingDistributor import ..DEBUGLEVEL_TDGP @@ -194,7 +197,7 @@ end # plotting stuff -function summaryAEMwithnuisanceimages(soundings::Array{S, 1}, opt_in::Options, optn_in::OptionsNuisance; +function summaryAEMimages(soundings::Array{S, 1}, opt_in::Options, optn_in::OptionsNuisance; zall=[-1.], qp1=0.05, qp2=0.95, @@ -425,7 +428,7 @@ function plotnuquant(nqlow, nqmid, nqhigh, nunominal, s, gridx, icol, nrows, ms= icol end -function plotindividualsoundings(soundings::Vector{S}, +function plotindividualAEMsoundings(soundings::Vector{S}, aem_in::Operator1D, opt_in::Options, optn_in::OptionsNuisance, idxplot; zall = [-1.], @@ -492,7 +495,9 @@ function plotindividualsoundings(soundings::Vector{S}, end end -plotindividualAEMsoundingswithnuisance = plotindividualsoundings - +# legacy backwards compat +plotindividualAEMsoundingswithnuisance = plotindividualAEMsoundings +plotindividualsoundings = plotindividualAEMsoundings +summaryAEMwithnuisanceimages = summaryAEMimages end # module diff --git a/src/AbstractOperator.jl b/src/AbstractOperator.jl index 1533f330..9c1f5631 100644 --- a/src/AbstractOperator.jl +++ b/src/AbstractOperator.jl @@ -40,5 +40,9 @@ function plotconvandlast end function setnuforinvtype end # for the Line type from LineRegression which is useful function getsmoothline end +# summary stats from McMC AEM inversions +function summaryAEMimages end +# plot predictions from posterior for McMC AEM inversions +function plotindividualAEMsoundings end export Operator, Operator1D, Operator2D, Sounding end diff --git a/src/CSEM1DInversion.jl b/src/CSEM1DInversion.jl index f345cd9e..5e798f2b 100644 --- a/src/CSEM1DInversion.jl +++ b/src/CSEM1DInversion.jl @@ -4,7 +4,7 @@ import ..AbstractOperator.get_misfit import ..CSEM1DEr.getfield! using ..AbstractOperator, ..CSEM1DEr using PyPlot, LinearAlgebra, ..CommonToAll - +import ..AbstractOperator.plotmodelfield! import ..Model, ..Options export CSEMRadialEr, plotmodelfield!, addnoise From c6084a268cee7f49d7a131b1deba30092f34e382 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 12 Jan 2024 18:00:07 +1100 Subject: [PATCH 2/5] fixed silly subplot bug for nuisances --- src/AEMwithNuisanceMcMCInversionTools.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AEMwithNuisanceMcMCInversionTools.jl b/src/AEMwithNuisanceMcMCInversionTools.jl index 4d83c551..e165712d 100644 --- a/src/AEMwithNuisanceMcMCInversionTools.jl +++ b/src/AEMwithNuisanceMcMCInversionTools.jl @@ -417,7 +417,7 @@ function plotnuquant(nqlow, nqmid, nqhigh, nunominal, s, gridx, icol, nrows, ms= nnu = min(size(nqlow, 1), size(nunominal, 1)) # in case we've inverted a zero bounds nuisance by mistech... labelnu[1] == "" || @assert length(labelnu) == nnu for inu = 1:nnu - s[icol] = subplot(nrows, 1, icol, sharex=s[icol-1]) + s[icol].sharex(s[icol-1]) s[icol].fill_between(gridx, nqlow[inu,:], nqhigh[inu,:], alpha=0.5) s[icol].plot(gridx, nqmid[inu,:]) s[icol].plot(gridx, nunominal[inu,:], "o", markersize=ms) From 5a2303cdf675d8cd3a5c80233bf0ee71ba58f4f6 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Thu, 18 Jan 2024 16:18:09 +1100 Subject: [PATCH 3/5] can now make vtk from multiline ASEG-GDF --- src/CommonToAll.jl | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index dfbbeeeb..c3db49f3 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1097,23 +1097,45 @@ function readcols(cols::Vector, fname::String; decfactor=1, startfrom=1, dotill= end function colstovtk(cols::Vector, fname::String; decfactor=1, hasthick=true) + # for one file in a direcoty X, Y, Z, σ, thick = readcols(cols, fname; decfactor) thick = thick[1,:] - zall = thicktodepth(thick; hasthick) + zall = thicktodepth(thick; hasthick) + fstring = basename(fname) + dstring = dirname(fname) + fpath = joinpath(dstring, "LEI_Line_"*fstring) + colstovtk(X, Y, Z, σ, zall, fpath) +end + +function colstovtk(X, Y, Z, σ, zall, fpath::String) Ni = length(X) Nk = length(zall) x = [X[i] for i = 1:Ni, j = 1:1, k = 1:Nk] y = [Y[i] for i = 1:Ni, j = 1:1, k = 1:Nk] z = [Z[i] - zall[k] for i = 1:Ni, j = 1:1, k = 1:Nk] σvtk = [σ[i, k] for i = 1:Ni, j = 1:1, k = 1:Nk] - fstring = basename(fname) - dstring = dirname(fname) - vtk_grid(joinpath(dstring, "LEI_Line_"*fstring), x, y, z) do vtk + vtk_grid(fpath, x, y, z) do vtk vtk["cond_LEI"] = log10.(σvtk) end nothing end +function colstovtk(cols::Dict, fname::String; decfactor=1, hasthick=true) + Xc, Yc, Zc, σc, thickc, linec = map(x->get(cols, x, 0), ["X", "Y", "Z", "cond", "thick", "line"]) + X, Y, Z, σ, thick, lines = readcols([Xc, Yc, Zc, σc, thickc, linec], fname; decfactor) + linenos = unique(Int.(lines)) + thick = thick[1,:] # assumes all thicknesses are same + zall = thicktodepth(thick; hasthick) + dstring = dirname(fname) + for l in linenos + idx = lines .== l + fstring = "$l" + @info "doing Line "*fstring + fpath = joinpath(dstring, "LEI_Line_"*fstring) + colstovtk(X[idx], Y[idx], Z[idx], σ[idx,:], zall, fpath) + end +end + function thicktodepth(thick; hasthick=true) if hasthick zb = [0; cumsum(thick)] From 66f57b2572006dddc6ed43894abfb7ba527e76c4 Mon Sep 17 00:00:00 2001 From: Anandaroop Ray Date: Fri, 19 Jan 2024 18:28:33 +1100 Subject: [PATCH 4/5] bugfix for log10 in new colstovtk --- src/CommonToAll.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/CommonToAll.jl b/src/CommonToAll.jl index c3db49f3..ab75f650 100644 --- a/src/CommonToAll.jl +++ b/src/CommonToAll.jl @@ -1096,31 +1096,32 @@ function readcols(cols::Vector, fname::String; decfactor=1, startfrom=1, dotill= end end -function colstovtk(cols::Vector, fname::String; decfactor=1, hasthick=true) +function colstovtk(cols::Vector, fname::String; decfactor=1, hasthick=true, islog10=false, prefix="") # for one file in a direcoty X, Y, Z, σ, thick = readcols(cols, fname; decfactor) thick = thick[1,:] zall = thicktodepth(thick; hasthick) fstring = basename(fname) dstring = dirname(fname) - fpath = joinpath(dstring, "LEI_Line_"*fstring) - colstovtk(X, Y, Z, σ, zall, fpath) + fpath = joinpath(dstring, "LEI_Line_"*prefix*fstring) + colstovtk(X, Y, Z, σ, zall, fpath, islog10) end -function colstovtk(X, Y, Z, σ, zall, fpath::String) +function colstovtk(X, Y, Z, σ, zall, fpath::String, islog10::Bool) Ni = length(X) Nk = length(zall) x = [X[i] for i = 1:Ni, j = 1:1, k = 1:Nk] y = [Y[i] for i = 1:Ni, j = 1:1, k = 1:Nk] z = [Z[i] - zall[k] for i = 1:Ni, j = 1:1, k = 1:Nk] σvtk = [σ[i, k] for i = 1:Ni, j = 1:1, k = 1:Nk] + T = islog10 ? x->x : x->log10(x) # if log 10 leave alone vtk_grid(fpath, x, y, z) do vtk - vtk["cond_LEI"] = log10.(σvtk) + vtk["cond_LEI"] = T.(σvtk) end nothing end -function colstovtk(cols::Dict, fname::String; decfactor=1, hasthick=true) +function colstovtk(cols::Dict, fname::String; decfactor=1, hasthick=true, islog10=false, prefix="") Xc, Yc, Zc, σc, thickc, linec = map(x->get(cols, x, 0), ["X", "Y", "Z", "cond", "thick", "line"]) X, Y, Z, σ, thick, lines = readcols([Xc, Yc, Zc, σc, thickc, linec], fname; decfactor) linenos = unique(Int.(lines)) @@ -1131,8 +1132,8 @@ function colstovtk(cols::Dict, fname::String; decfactor=1, hasthick=true) idx = lines .== l fstring = "$l" @info "doing Line "*fstring - fpath = joinpath(dstring, "LEI_Line_"*fstring) - colstovtk(X[idx], Y[idx], Z[idx], σ[idx,:], zall, fpath) + fpath = joinpath(dstring, "LEI_Line_"*prefix*fstring) + colstovtk(X[idx], Y[idx], Z[idx], σ[idx,:], zall, fpath, islog10) end end From 5f77f17705bc41d899d2a409307a2ed185de34ff Mon Sep 17 00:00:00 2001 From: a2ray <52647259+a2ray@users.noreply.github.com> Date: Mon, 22 Jan 2024 10:24:39 +1100 Subject: [PATCH 5/5] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3d6aa02d..73eabf2a 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.0" +version = "0.4.2" [deps] Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"