Merge pull request #65 from GeoscienceAustralia/dev

Dev

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.1"
+version = "0.4.2"
 
 [deps]
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"

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)
+transD_GP.plotmodelfield!(tempest, mtest, nutest)

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

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,
@@ -414,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)
@@ -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
 

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

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

src/CommonToAll.jl

@@ -1096,24 +1096,47 @@ 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)   
+    zall = thicktodepth(thick; hasthick)
+    fstring = basename(fname)
+    dstring = dirname(fname)
+    fpath = joinpath(dstring, "LEI_Line_"*prefix*fstring)
+    colstovtk(X, Y, Z, σ, zall, fpath, islog10)
+end
+
+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]
-    fstring = basename(fname)
-    dstring = dirname(fname)
-    vtk_grid(joinpath(dstring, "LEI_Line_"*fstring), x, y, z) do vtk
-        vtk["cond_LEI"]  = log10.(σvtk)
+    T = islog10 ?  x->x : x->log10(x) # if log 10 leave alone
+    vtk_grid(fpath, x, y, z) do vtk
+        vtk["cond_LEI"]  = T.(σvtk)
     end
     nothing
 end
 
+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))
+    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_"*prefix*fstring)
+        colstovtk(X[idx], Y[idx], Z[idx], σ[idx,:], zall, fpath, islog10)
+    end
+end    
+
 function thicktodepth(thick; hasthick=true)
     if hasthick
         zb = [0; cumsum(thick)]