Skip to content

Commit

Permalink
Merge pull request #91 from gher-ulg/errormap
Browse files Browse the repository at this point in the history 
Error map (scpme)
  • Loading branch information
@Alexander-Barth
Alexander-Barth authored Dec 8, 2021
2 parents 4375a03 + 8777608 commit 68e6396
Show file tree
Hide file tree
Showing 3 changed files with 195 additions and 21 deletions.
6 changes: 5 additions & 1 deletion src/DIVAnd_diagapp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,))
Expand Down Expand Up @@ -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)])
Expand Down
190 changes: 173 additions & 17 deletions src/DIVAnd_errormap.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 == ()
Expand All @@ -28,6 +27,7 @@ function DIVAnd_errormap(
Bigdata = false
Lowdata = false


if method == "auto"

Lpmnrange = DIVAnd_Lpmnrange(pmn, len)
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
20 changes: 17 additions & 3 deletions src/DIVAnd_scalecpme!.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

0 comments on commit 68e6396

Please sign in to comment.