Skip to content

Commit

Permalink
speed up
Browse files Browse the repository at this point in the history
  • Loading branch information
@cspencerjones
cspencerjones committed Jan 16, 2024
1 parent 2e10e84 commit f04fa98
Show file tree
Hide file tree
Showing 2 changed files with 88 additions and 43 deletions.
85 changes: 61 additions & 24 deletions scripts/mit_inv.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,13 @@ flush(stdout)
# The first guess for the tracer concentration should be close to the actual tracer concentration
# take first guess as θtrue
cvec=(vec(θtrue))
θguess = unvec(θtrue,cvec)

#θguess = unvec(θtrue,cvec)

ucnew = load("ucnew.jld")["ucnew"]
#first guess tracer control vector is near zero, and we want this to remain relatively small
u = Field(-0.01.*ones(size(γ.wet)),γ,θtrue.name,θtrue.longname,θtrue.units)
uvec = vec(u)
u1 = Field(0.0.*ones(size(γ.wet)),γ,θtrue.name,θtrue.longname,θtrue.units)
uvec = ucnew
u = unvec(u1,ucnew)
q = zeros(size(uvec))
q[surfind] = cvec[surfind]

Expand All @@ -63,8 +64,10 @@ Q⁻ = Diagonal(Qdiag)

ind_1 = load("ind_1.jld")["ind_1"]
ind_2 = load("ind_2.jld")["ind_2"]
ufnew = load("ufnew.jld")["ufnew"]


A0 = sparse(ind_1, ind_2, ones(length(ind_1))*0.1)
Adummy = sparse(ind_1, ind_2, ones(length(ind_1))*0.1)


#test: apply A0 a bunch of times
Expand All @@ -91,9 +94,11 @@ A0 = sparse(ind_1, ind_2, ones(length(ind_1))*0.1)
#readline()


non_zero_indices1, non_zero_indices2, non_zero_values = findnz(A0)
non_zero_indices1, non_zero_indices2, non_zero_values = findnz(Adummy)

Afguess = sparse(non_zero_indices1, non_zero_indices2, ufnew)

non_zero_indices1, non_zero_indices2, non_zero_values = findnz(Afguess)

non_zero_indices = hcat(non_zero_indices1, non_zero_indices2)

Expand All @@ -105,18 +110,19 @@ print(ulength)
print(sum(γ.wet))

# get sample J value
F = costfunction_gridded_model(convec,non_zero_indices,u,A0,ctrue,cvec,q,W⁻,Q⁻,γ)
fg!(F,G,x) = costfunction_gridded_model!(F,G,x,non_zero_indices,u,A0,ctrue,cvec,q,W⁻,Q⁻,γ)
fg(x) = costfunction_gridded_model(x,non_zero_indices,u,A0,ctrue,cvec,q,W⁻,Q⁻,γ)
F = costfunction_gridded_model(convec,non_zero_indices,u,Adummy,ctrue,cvec,q,W⁻,Q⁻,γ)
fg!(F,G,x) = costfunction_gridded_model!(F,G,x,non_zero_indices,u,Adummy,ctrue,cvec,q,W⁻,Q⁻,γ)
fg(x) = costfunction_gridded_model(x,non_zero_indices,u,Adummy,ctrue,cvec,q,W⁻,Q⁻,γ)
f(x) = fg(x)[1]
J₀,gJ₀ = fg(convec)

#### gradient check ###################
# check with forward differences
ϵ = 1e-3
println(size(length(convec)))
ii = rand(1:length(convec))
ii=748529+10000
#ii = rand(1:length(convec))
#ii=748529+10000
ii=4000
println("Location for test =",ii)
flush(stdout)
δu = copy(convec); δu[ii] += ϵ
Expand All @@ -132,47 +138,63 @@ println("Percent error ",100*abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite))
flush(stdout)
#### end gradient check #################

#print(length(convec))
print(length(convec))
# filter the data with an Optim.jl method
iterations = 2
iterations = 20000
out = steadyclimatology_optim(convec,fg!,iterations)

# reconstruct by hand to double-check.
ũ = unvec(θguess,Float32.(out.minimizer[begin:ulength]))
#ũ = unvec(θguess,Float32.(out.minimizer[begin:ulength]))

#println(typeof(ũ))

# reconstruct tracer map
c₀ = θguess
c̃ = θguess+ũ
#c₀ = θguess
#c̃ = θguess+ũ

Δc̃ = c̃ - θtrue
Δc₀ = θguess - θtrue
#Δc̃ = c̃ - θtrue
#Δc₀ = θguess - θtrue

Anew = A0 + sparse(non_zero_indices[:, 1], non_zero_indices[:, 2], out.minimizer[ulength+1:end])
Anew = Adummy + sparse(non_zero_indices[:, 1], non_zero_indices[:, 2], out.minimizer[ulength+1:end])

maxA = maximum(Anew)
minA = minimum(Anew)

println("Max A:$maxA")
println("Min A:$minA")

ucout = out.minimizer[begin:ulength]

maxuc = maximum(ucout)
minuc = minimum(ucout)

println("Max uc:$maxuc")
println("Min uc:$minuc")


flush(stdout)
onesvec = ones(size(q))

#oldf = sum((non_zero_values).^2)
#newf = sum((out.minimizer[ulength+1:end]).^2)
tracer_cons1 = sum((Q⁻*A0*cvec - q).^2)
tracer_cons2 = sum((Q⁻*Anew*(cvec+out.minimizer[begin:ulength]) - q).^2)
mass_cons1 = sum((A0*onesvec).^2)
tracer_cons1 = sum(((Adummy+Afguess)*(cvec+u) - q).^2)
tracer_cons2 = sum((Anew*(cvec+out.minimizer[begin:ulength]) - q).^2)
mass_cons1 = sum(((Adummy+ Afguess)*onesvec).^2)
mass_cons2 = sum((Anew*onesvec).^2)

println("old tracer cons:$tracer_cons1")
println("new tracer cons:$tracer_cons2")
println("old mass cons:$mass_cons1")
println("new mass cons:$mass_cons2")

#b = getsurfaceboundary(θtrue)
b = getsurfaceboundary(θtrue)

save("ufnew.jld","ufnew", out.minimizer[ulength+1:end])
save("ucnew.jld","ucnew", out.minimizer[begin:ulength])





#A0

Expand All @@ -181,7 +203,22 @@ println("new mass cons:$mass_cons2")
#Alu2 = lu(Anew)


#final_θ = steadyinversion(Anew,b,γ)
final_θ = steadyinversion(Anew,b,γ)
#

c1 = (vec(final_θ))

test = sum(Anew * c1 - q)


max_θ = maximum(filter(!isnan,final_θ.tracer))
min_θ = minimum(filter(!isnan,final_θ.tracer))

println("final_θ max:$max_θ")
println("final_θ min:$min_θ")
println("test:$test")

#writefield("final_theta.nc",final_θ)

# plot the difference
#level = 15 # your choice 1-33
Expand Down
46 changes: 27 additions & 19 deletions src/TMI.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1672,22 +1672,26 @@ function costfunction_gridded_model(convec,non_zero_indices,u₀::Field{T},A0,y:
#dummy,dummy,fguess = findnz(A0)
#dummy,dummy,fnow = findnz(A)
onesvec = ones(size(uvec))
#csum = Wⁱ * uvec+c
csum = uvec+c


# find lagrange multipliers
#muk = transpose(A) * Qⁱ * (A * c - q)
dAcdf = spzeros(length(ufvec),length(c))
dA1df = spzeros(length(ufvec),length(c))
muk = transpose(A) * Qⁱ * (A * csum - q)
dAcdf = spzeros(length(ufvec))#,length(c))
dA1df = spzeros(length(ufvec))#,length(c))
A1 = (A * onesvec)
Acq = Qⁱ*(A * csum - q)
for ii in eachindex(ufvec)
dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
dA1df[ii,non_zero_indices[ii, 1]] = 1
#dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
#dA1df[ii,non_zero_indices[ii, 1]] = 1
dAcdf[ii] = sum(csum[non_zero_indices[ii, 2]].*Acq[non_zero_indices[ii, 1]])
dA1df[ii] = sum(A1[non_zero_indices[ii, 1]])
end


dAdf_terms = dA1df * (A * onesvec) + dAcdf * Qⁱ * (A * c - q)# + dA1df * (A * onesvec)
dAdf_terms = dA1df + dAcdf #* Qⁱ * (A * c - q)# + dA1df * (A * onesvec)

J = 10^3 * uvec ⋅ uvec + transpose(A * c - q) * Qⁱ * (A*c - q) +
J = 10 *uvec ⋅ uvec + transpose(A * csum - q) * Qⁱ * (A*csum - q) +
#2 * transpose(muk)*( Wⁱ * uvec+c-y)# + sum(ufvec)#sum(A*onesvec)
transpose(A * onesvec)* (A* onesvec)# + 10 * ufvec ⋅ ufvec.^3

Expand All @@ -1698,7 +1702,7 @@ function costfunction_gridded_model(convec,non_zero_indices,u₀::Field{T},A0,y:
if ii <= ulength
#this is the derivative of the cost function wrt the part of the control vector
# associated with the tracer concentration
guvec[ii] = 2 * 10^3 * uvec[ii] #- (2 * transpose(muk) * Wⁱ)[ii]
guvec[ii] = 20 * uvec[ii] + 2 * transpose(muk)[ii]
else
#this is the derivative of the cost function wrt the part of the control vector
# associated with the transport vector
Expand All @@ -1723,33 +1727,37 @@ function costfunction_gridded_model!(J,guvec,convec::Vector{T},non_zero_indices,
#dummy,dummy,fguess = findnz(A0)
#dummy,dummy,fnow = findnz(A)
onesvec = ones(size(uvec))
#csum = Wⁱ * uvec+c
csum = uvec+c

# find lagrange multipliers
#muk = transpose(A) * Qⁱ * (A * c - q)
dAcdf = spzeros(length(ufvec),length(c))
dA1df = spzeros(length(ufvec),length(c))
muk = transpose(A) * Qⁱ * (A * csum - q)
dAcdf = spzeros(length(ufvec))#,length(c))
dA1df = spzeros(length(ufvec))#,length(c))
A1 = (A * onesvec)
Acq = Qⁱ*(A * csum - q)
for ii in eachindex(ufvec)
dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
dA1df[ii,non_zero_indices[ii, 1]] = 1
#dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
#dA1df[ii,non_zero_indices[ii, 1]] = 1
dAcdf[ii] = sum(csum[non_zero_indices[ii, 2]].*Acq[non_zero_indices[ii, 1]])
dA1df[ii] = sum(A1[non_zero_indices[ii, 1]])
end
dAdf_terms = dA1df * (A * onesvec) + dAcdf * Qⁱ * (A * c - q)# + dA1df * (A * onesvec)
dAdf_terms = dA1df + dAcdf #* Qⁱ * (A * c - q)# + dA1df * (A * onesvec)

if guvec != nothing
tmp = guvec
for (ii,vv) in enumerate(tmp)
if ii <= ulength
guvec[ii] = 2 * 10^3 * uvec[ii]#-(2 * transpose(muk) * Wⁱ)[ii]
guvec[ii] = 20 * uvec[ii]+2 * muk[ii]
else
guvec[ii]= 2 * dAdf_terms[ii-ulength]# + 4 * 10 * convec[ii].^3
end
end
end

if J !=nothing
return 10^3 * uvec ⋅ uvec + transpose(A * c - q) * Qⁱ * (A*c - q)+
return 10 * uvec ⋅ uvec + transpose(A * csum - q) * Qⁱ * (A*csum - q)+
#2 * transpose(muk)*( Wⁱ * uvec+c-y)+ sum(ufvec)#sum(A*onesvec)
transpose(A * onesvec) * (A* onesvec)# + 10 * ufvec ⋅ ufvec.^3
transpose(A * onesvec) * (A* onesvec)# + 10 * ufvec ⋅ ufvec.^3
end
end

Expand Down

0 comments on commit f04fa98

Please sign in to comment.