From 01fe56e0a8ec217a350b37e6a0aef48aedda3085 Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Mon, 6 Dec 2021 11:34:16 +0100 Subject: [PATCH] new error map --- src/DIVAnd_diagapp.jl | 6 +- src/DIVAnd_errormap.jl | 190 +++++++++++++++++++++++++++++++++++---- src/DIVAnd_scalecpme!.jl | 20 ++++- 3 files changed, 195 insertions(+), 21 deletions(-) diff --git a/src/DIVAnd_diagapp.jl b/src/DIVAnd_diagapp.jl index c2a36be1..b611bfc4 100644 --- a/src/DIVAnd_diagapp.jl +++ b/src/DIVAnd_diagapp.jl @@ -101,6 +101,7 @@ function DIVAnd_diagapp( # in 2D needs to be adapted to cover the domain: pack unpack and sum on real domain ? # Allocate arrays once eij = zeros(Int, size(pmn[1])) + #@show size(eij) diagerror = zeros(Float64, size(pmn[1])) .* NaN tutuu = zeros(Float64, size(pmn[1])) tutu = statevector_pack(sv, (eij,)) @@ -169,10 +170,13 @@ function DIVAnd_diagapp( tutuub[:], = statevector_unpack(sv, zy) end # Now each point of eij sums up the contribution aournd its box center - for IG in findall(x -> x == 1, eij) + for IGtemp in findall(x -> x == 1, eij) + IG=CartesianIndex(IGtemp) #################################### # NEED TO TAKE INTO ACCOUNT MODDIM #################################### + #@show IG, IFI, mystep, ILA + #@show typeof(IG), typeof(IFI), typeof(mystep), typeof(ILA) diagerror[IG] = sum(tutuu[max(IFI, IG - mystep):min(IG + mystep, ILA)]) if Binv diagB[IG] = sum(tutuub[max(IFI, IG - mystep):min(IG + mystep, ILA)]) diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index c3e93912..3e9a75c8 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -13,13 +13,12 @@ function DIVAnd_errormap( ) # Criteria to define which fraction of the domain size L can be to be called small - LoverLlimit = 0.12 + LoverLlimit = 0.1 # Criteria for lot of data means lot of data in hypersphere of correlation length diameters. Need to think about L=0 case ... pointsperbubblelimit = 10 pointsperbubblelimitlow = 1 - errmethod = method ScalebyB = Bscale noP = s.P == () @@ -28,6 +27,7 @@ function DIVAnd_errormap( Bigdata = false Lowdata = false + if method == "auto" Lpmnrange = DIVAnd_Lpmnrange(pmn, len) @@ -40,7 +40,7 @@ function DIVAnd_errormap( LoverLdomain[i] = Lpmnrange[i][2] / size(mask)[i] end - if sum(Loverdomain .< LoverLlimit) == ndims(masl) + if sum(LoverLdomain .< LoverLlimit) == ndims(mask) smallL = true end @@ -49,22 +49,22 @@ function DIVAnd_errormap( realdims = ndims(mask) for i = 1:ndims(mask) LoverLdomain[i] = Lpmnrange[i][1] / size(mask)[i] + + if Lpmnrange[i][1] == 0 LoverLdomain[i] = 1.0 / size(mask)[i] realdims = realdims - 1 end end - - if prod(Loverdomain) * ndata > pointsperbubblelimit^realdims + #nbdonnee size of f a revoir en fonction depsilon2 + if prod(LoverLdomain)*(pi^realdims)/gamma((realdims/2)+1) * size(f)[1] > pointsperbubblelimit Bigdata = true end - if prod(Loverdomain) * ndata < pointsperbubblelimitlow^realdims + if prod(LoverLdomain)*(pi^realdims)/gamma((realdims/2)+1) * size(f)[1] < pointsperbubblelimitlow Lowdata = true end - - # try to guess # small L @@ -83,30 +83,128 @@ function DIVAnd_errormap( else errmethod = "diagapp" end - end - if Bigdata + else + if Bigdata errmethod = "scpme" else errmethod = "aexerr" + end + + end + + # So best guess up to now + @show errmethod + end + + + + if method == "cheap" + Lpmnrange = DIVAnd_Lpmnrange(pmn, len) + # L compared to domain size + + LoverLdomain = zeros(Float64, ndims(mask)) + + + for i = 1:ndims(mask) + LoverLdomain[i] = Lpmnrange[i][2] / size(mask)[i] + end + + if sum(LoverLdomain .< LoverLlimit) == ndims(mask) + smallL = true + end + + + # Now look at lower values to check for data coverage + realdims = ndims(mask) + for i = 1:ndims(mask) + LoverLdomain[i] = Lpmnrange[i][1] / size(mask)[i] + if Lpmnrange[i][1] == 0 + LoverLdomain[i] = 1.0 / size(mask)[i] + realdims = realdims - 1 end + end + #nbdonnee size of f a revoir en fonction depsilon2 + if prod(LoverLdomain) * size(f)[1] > pointsperbubblelimit^realdims + Bigdata = true + end + if prod(LoverLdomain) * size(f)[1] < pointsperbubblelimitlow^realdims + Lowdata = true + end + if smallL + if Lowdata + errmethod = "cpme" + else + if Bigdata + errmethod = "scpme" + else + errmethod = "cpme" + end + end else + if Bigdata + errmethod = "scpme" + else + errmethod = "cpme" + end + end + end + if method == "precise" + Lpmnrange = DIVAnd_Lpmnrange(pmn, len) + # L compared to domain size - end + LoverLdomain = zeros(Float64, ndims(mask)) + for i = 1:ndims(mask) + LoverLdomain[i] = Lpmnrange[i][2] / size(mask)[i] + end + if sum(LoverLdomain .< LoverLlimit) == ndims(mask) + smallL = true + end + # Now look at lower values to check for data coverage + realdims = ndims(mask) + for i = 1:ndims(mask) + LoverLdomain[i] = Lpmnrange[i][1] / size(mask)[i] + if Lpmnrange[i][1] == 0 + LoverLdomain[i] = 1.0 / size(mask)[i] + realdims = realdims - 1 + end + end + #nbdonnee size of f a revoir en fonction depsilon2 + if prod(LoverLdomain) * size(f)[1] > pointsperbubblelimit^realdims + Bigdata = true + end - # So best guess up to now - @show errmethod - end + if prod(LoverLdomain) * size(f)[1] < pointsperbubblelimitlow^realdims + Lowdata = true + end + if smallL + if Lowdata + errmethod = "aexerr" + else + if Bigdata + errmethod = "diagapp" + else + errmethod = "diagapp" + end + end + else + if Bigdata + errmethod = "diagapp" + else + errmethod = "aexerr" + end + end + end if errmethod == "cpme" && Bscale @@ -141,13 +239,71 @@ function DIVAnd_errormap( end + # Now calculate error depening on the method + @show errmethod, ScalebyB, pointsperbubblelimit + + if errmethod == "cpme" + errormap = DIVAnd_cpme( + mask, + pmn, + xi, + x, + f, + len, + epsilon2; + otherargs... + ) + + return errormap, errmethod + end - # Now calculate error depening on the method + if errmethod == "scpme" + errormap = DIVAnd_cpme( + mask, + pmn, + xi, + x, + f, + len, + epsilon2; + otherargs... + ) + + scpme=deepcopy(errormap) + @time DIVAnd_scalecpme!(scpme,s.P) + + return scpme, errmethod + end - @show errmethod, ScalebyB + if errmethod == "exact" + errormap =statevector_unpack(s.sv,diag(s.P)) + return errormap, errmethod + end - return errormap, errmethod + if errmethod == "diagapp" + errormap = DIVAnd_diagapp( + s.P, + pmn, + len, + s.sv, + ) + return errormap, errmethod + end + if errmethod == "aexerr" + errormap,bi,c,d =DIVAnd_aexerr( + mask, + pmn, + xi, + x, + f, + len, + epsilon2; + otherargs... + ) + return errormap, errmethod + end + @show "you should not be here" end diff --git a/src/DIVAnd_scalecpme!.jl b/src/DIVAnd_scalecpme!.jl index 3d729f09..92c25226 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -2,7 +2,10 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7) # IN PLACE rescaling of the clever poor mans estimate using a randomized estimate of the analysis error covariance # P returned in structure s (so s.P) from a previous run # nsamples is the number of random arrays used to estimate the value - z = randn((size(P)[1], nsamples)) + + fractionshift=0.5 + + z = randn((size(P)[1], nsamples)) errscale = 1 if P.factors != nothing ZZ = P.factors.PtL \ z @@ -11,7 +14,18 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7) ZZ = P * z errscale = mean(diag(z' * ZZ) ./ diag(z' * z)) end - errscale = errscale / mean(cpme[.!isnan.(cpme)]) - cpme[:] .= cpme[:] * errscale + # errscale here is the target + + # oldvalue is what we had + oldvalue=mean(cpme[.!isnan.(cpme)]) + + #try a shift of fraction fractionshift of the mismatch + + cpme[:] .= cpme[:] .+ fractionshift*(errscale-oldvalue) + + # the rest is corrected by the multiplicative scaling + + + cpme[:] .= cpme[:] * errscale / mean(cpme[.!isnan.(cpme)]) return cpme end