Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Precompute mismatch scaling factor #60

Open
briochemc opened this issue Feb 8, 2021 · 0 comments
Open

Precompute mismatch scaling factor #60

briochemc opened this issue Feb 8, 2021 · 0 comments

Comments

@briochemc
Copy link
Member

Specifically, refactor these to allow the user to provide transpose(o) * W * o instead of recalculating it every time:

## new functions for more generic obs packages
# TODO Add an optional function argument to transform the data before computingn the mismatch
# Example if for isotope tracers X where one ususally wants to minimize the mismatch in δ or ε.
function mismatch(x, grd::OceanGrid, obs; c=identity, W=I, M=interpolationmatrix(grd, obs.metadata), iwet=iswet(grd, obs))
o = view(obs, iwet)
δx = M * c(x) - o
return 0.5 * transpose(δx) * W * δx / (transpose(o) * W * o)
end
mismatch(x, grd::OceanGrid, ::Missing; kwargs...) = 0
function ∇mismatch(x, grd::OceanGrid, obs; c=identity, W=I, M=interpolationmatrix(grd, obs.metadata), iwet=iswet(grd, obs))
∇c = Diagonal(ForwardDiff.derivative-> c(x .+ λ), 0.0))
o = view(obs, iwet)
δx = M * c(x) - o
return transpose(W * δx) * M * ∇c / (transpose(o) * W * o)
end
∇mismatch(x, grd::OceanGrid, ::Missing; kwargs...) = transpose(zeros(length(x)))
# In case the mismatch is not based on the tracer but on some function of it
function indirectmismatch(xs::Tuple, grd::OceanGrid, modify::Function, obs, i, M=interpolationmatrix(grd, obs[i].metadata), iwet=iswet(grd, obs[i]))
x2 = modify(xs...)
out = 0.0
M = interpolationmatrix(grd, obs[i].metadata)
iwet = iswet(grd, obs[i])
o = view(obs[i], iwet)
δx = M * x2[i] - o
return 0.5 * transpose(δx) * δx / (transpose(o) * o)
end
function ∇indirectmismatch(xs::Tuple, grd::OceanGrid, modify::Function, obs, i, M=interpolationmatrix(grd, obs[i].metadata), iwet=iswet(grd, obs[i]))
nt, nb = length(xs), length(iswet(grd))
x2 = modify(xs...)
o = view(obs[i], iwet)
δx = M * x2[i] - o
∇modᵢ = ∇modify(modify, xs, i)
return transpose(δx) * M * ∇modᵢ / (transpose(o) * o)
end

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant