Skip to content

Commit 4710daa

Browse files
inverseprobleminverseproblem
committed
add test for Gauss-Newton
1 parent 63820e2 commit 4710daa

File tree

6 files changed

+238
-19
lines changed

6 files changed

+238
-19
lines changed

‎.github/workflows/Documentation.yml renamed to ‎.github/workflows/documentation.yml

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,18 @@
1-
name: Documentation
1+
name: Build documentation
22

33
on:
44
push:
55
branches:
66
- main # update to match your development branch (master, main, dev, trunk, ...)
7-
tags: '*'
7+
- develop
8+
tags:
9+
- 'v*'
810
pull_request:
11+
branches:
12+
- develop
913

1014
jobs:
11-
build:
15+
build_docs:
1216
permissions:
1317
actions: write
1418
contents: write
@@ -20,7 +24,9 @@ jobs:
2024
- uses: julia-actions/setup-julia@v2
2125
with:
2226
version: '1'
23-
- uses: julia-actions/cache@v2 # reduce GitHub Actions running time
27+
- name: Julia cache
28+
if: ${{ !env.ACT }}
29+
uses: julia-actions/cache@v2
2430
- name: Install dependencies
2531
run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
2632
- name: Build docs and deploy

‎.github/workflows/test.yml

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
name: Package tests
2+
3+
on:
4+
push:
5+
branches:
6+
- main
7+
- develop
8+
tags: '*'
9+
pull_request:
10+
branches:
11+
- develop
12+
types:
13+
- ready_for_review
14+
- synchronize
15+
16+
# needed to allow julia-actions/cache to delete old caches that it has created
17+
permissions:
18+
actions: write
19+
contents: read
20+
21+
jobs:
22+
run_tests:
23+
name: Julia ${{ matrix.julia-version }} - ${{ matrix.os }} - ${{ matrix.arch }}
24+
runs-on: ${{ matrix.os }}
25+
if : github.event.pull_request.draft == false
26+
strategy:
27+
matrix:
28+
julia-version: ['1'] #['lts', '1', 'pre']
29+
arch: [x64] #[x64, x86]
30+
os: [ubuntu-latest, macOS-latest] #[ubuntu-latest, windows-latest, macOS-latest]
31+
include:
32+
- os: macOS-latest
33+
julia-version: '1'
34+
arch: aarch64
35+
36+
steps:
37+
- uses: actions/checkout@v4
38+
- uses: julia-actions/setup-julia@v2
39+
with:
40+
version: ${{ matrix.julia-version }}
41+
arch: ${{ matrix.arch }}
42+
- name: Julia cache
43+
if: ${{ !env.ACT }}
44+
uses: julia-actions/cache@v2
45+
- uses: julia-actions/julia-buildpkg@v1
46+
- uses: julia-actions/julia-runtest@v1
47+
# with:
48+
# annotate: true

‎src/Optimizers/GaussNewton/gauss_newton.jl

Lines changed: 16 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,8 @@ function gaussnewton(calcfwd!::Function,
5757
end
5858
@assert length(xprior)==length(x0)
5959

60+
@info "\n*** Gauss-Newton optimization *** "
61+
6062
## init
6163
N = size(invCd,1)
6264
M = length(x0)
@@ -92,8 +94,8 @@ function gaussnewton(calcfwd!::Function,
9294
res_m .= xcur - xprior
9395
objval += 0.5 * dot(res_m,invCm,res_m)
9496
##=======================
95-
# Gradient
96-
# Jtres = J*res, i.e., the gradient
97+
# Gradient
98+
# grad = J^T * Cd^-1 * res_d + Cm^-1 * res_m
9799
mul!(tmpgr_d,invCd,res_d)
98100
mul!(grad,transpose(jac),tmpgr_d)
99101
mul!(tmpgr_m,invCm,res_m)
@@ -116,14 +118,15 @@ function gaussnewton(calcfwd!::Function,
116118
# and get the in-place jacobian as an optional argument
117119
misf[1] = fh!(grad,x[1])
118120
g0norm = norm(grad)
119-
121+
122+
@info "Initial misfit $(misf[1]) "
120123
## Loop
121124
for k=1:maxiter
122125

123126
##===============================
124127
# Jacobian
125128
calcjac!(jac,x[k])
126-
# Now (J^T*invCd*J + invCm) p_gn = - (J^T*invCm) *res
129+
# Now (J^T*invCd*J + invCm) p_gn = - (J^T*invCd) *res
127130
# invCd * J
128131
mul!(invCd_J,invCd,jac)
129132
# temporary store invCm in the Hessian array
@@ -132,13 +135,15 @@ function gaussnewton(calcfwd!::Function,
132135
# J^T * invCd * J + invCm
133136
mul!(H,transpose(jac),invCd_J,1.0,1.0)
134137
# solve linear system
135-
faJtJ = factorize(Symmetric(H))
136-
# J^t*J*p_gn = -J^t*res
137-
# pkgn is the descent direction
138+
#faJtJ = factorize(Symmetric(H)) # this seems to get singular...
139+
faJtJ = lu(H)
140+
## (J^T*invCd*J + invCm) p_gn = - (J^T*invCd) *res
141+
## pkgn is the descent direction
138142
ldiv!(pkgn,faJtJ,-grad)
143+
##===============================
139144

140145
##===============================
141-
# Line search
146+
## Line search
142147
if k==1
143148
α0 = target_update
144149
else
@@ -150,10 +155,10 @@ function gaussnewton(calcfwd!::Function,
150155
α0=α0,maxiterwolfe=maxiterwolfe,
151156
maxiterzoom=maxiterzoom,
152157
c1=c1,c2=c2)
153-
154-
##===============================
155158
# Update the solution
156-
x[k+1] .= x0αp #x[k] .+ α.*pkgn
159+
x[k+1] .= x0αp ###x[k] .+ α.*pkgn
160+
##===============================
161+
157162

158163
if bounds!=nothing
159164
# project x[k+1]

‎test/runtests.jl

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,17 +27,29 @@ with_logger(logger) do
2727
@test test_bfgs2()
2828

2929
end
30+
31+
@testset "Test Gauss-Newton algo" begin
32+
33+
printstyled("Gauss-Newton: Test 1\n", bold=true,color=:cyan)
34+
@test test_gaussnewton1()
35+
36+
printstyled("Gauss-Newton: Test 2\n", bold=true,color=:cyan)
37+
@test test_gaussnewton2()
38+
39+
printstyled("Gauss-Newton: Test 3\n", bold=true,color=:cyan)
40+
@test test_gaussnewton3()
41+
end
3042
end
3143

3244

3345
# KronLinInv
3446
@testset "Test KronLinInv" begin
3547

3648
printstyled("KronLinInv: Testing 2D example \n", bold=true,color=:cyan)
37-
@test test2D()
49+
@test test_KLI2D()
3850

3951
printstyled("KronLinInv: Testing 3D example \n", bold=true,color=:cyan)
40-
@test test3D()
52+
@test test_KLI3D()
4153

4254
println()
4355
end

‎test/test_gaussnewton.jl

Lines changed: 148 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,148 @@
1+
2+
3+
using LinearAlgebra
4+
## using GLMakie
5+
6+
7+
function test_gaussnewton1()
8+
9+
tval = −2.903534
10+
tvec = [tval,tval] # location of the minimum
11+
stmin = -39.16617*length(tvec) # value at the minimum
12+
13+
function calcfwd!(u,x)
14+
u .= styblinskitang(x)
15+
return
16+
end
17+
18+
function calcjac!(g,x)
19+
g[1,:] .= ∇styblinskitang(x)
20+
return
21+
end
22+
23+
xprior = [-0.5,-4.5]
24+
25+
obsdata = [stmin]
26+
invCd = [1.0;;]
27+
invCm = [0.0 0;
28+
0 0.0]
29+
30+
maxiter = 30
31+
32+
xout,misfout = gaussnewton(calcfwd!,calcjac!,
33+
obsdata=obsdata,
34+
invCd=invCd,
35+
invCm=invCm,
36+
xprior=xprior,
37+
maxiter=maxiter)
38+
39+
ce1 =isapprox(xout[end],tvec,rtol=1e-2)
40+
return ce1
41+
end
42+
43+
function test_gaussnewton2()
44+
45+
# generate function
46+
function fwd(a::Real,b::Real,x::Vector{Float64})
47+
return a.*x./(b.+x)
48+
end
49+
50+
function calcfwd!(u,m)
51+
u .= fwd(m[1],m[2],x)
52+
return nothing
53+
end
54+
55+
function calcjac!(jac,m)
56+
aj = m[1]
57+
bj = m[2]
58+
grad_a = x ./ (bj.+x)
59+
grad_b = - ((aj.*x) ./ ((bj .+ x).^2))
60+
jac .= [grad_a grad_b]
61+
return
62+
end
63+
64+
65+
a,b = 2,3
66+
x = collect(range(0,5,length=25))
67+
y = fwd(a,b,x)
68+
69+
ab0 = [1.0,5.0]
70+
xprior = copy(ab0)
71+
obsdata = y
72+
invCd = diagm(1.0 .* ones(length(obsdata)))
73+
invCm = diagm(0.0 .* ones(length(ab0)))
74+
75+
maxiter = 30
76+
77+
xout,misfout = gaussnewton(calcfwd!,calcjac!,
78+
obsdata=obsdata,
79+
invCd=invCd,invCm=invCm,
80+
xprior=xprior,maxiter=maxiter,
81+
target_update=3.0)
82+
83+
ce1 = xout[end] ≈ [a,b]
84+
return ce1
85+
end
86+
87+
88+
function test_gaussnewton3()
89+
90+
function calcfwd!(u,x)
91+
u .= sum(x.^2)
92+
return
93+
end
94+
95+
function calcjac!(J,x)
96+
J[1,:] .= 2 .* x
97+
return
98+
end
99+
100+
xprior = [10.0,-13.5]
101+
102+
obsdata = [0.0]
103+
invCd = [1.0;;]
104+
invCm = [0.0 0;
105+
0 0.0]
106+
107+
maxiter = 50
108+
109+
xout,misfout = gaussnewton(calcfwd!,calcjac!,
110+
obsdata=obsdata,
111+
invCd=invCd,
112+
invCm=invCm,
113+
xprior=xprior,
114+
maxiter=maxiter)
115+
116+
# N = size(invCd,1)
117+
# M = length(xprior)
118+
# u_calc = zeros(eltype(xprior),N)
119+
# grad = zeros(eltype(xprior),M)
120+
# jac = zeros(eltype(xprior),N,M)
121+
122+
# function fh_bfgs!(grad,xcur)
123+
# # Calculated data
124+
# calcfwd!(u_calc,xcur)
125+
# # Jacobian
126+
# calcjac!(jac,xcur)
127+
# ##=======================
128+
# # Value of objective function
129+
# # misfit
130+
# res_d = u_calc - obsdata
131+
# objval = 0.5 * dot(res_d,invCd,res_d)
132+
# # prior term
133+
# res_m = xcur - xprior
134+
# objval += 0.5 * dot(res_m,invCm,res_m)
135+
# grad .= transpose(jac) * invCd * res_d + invCm * res_m
136+
# return objval
137+
# end
138+
139+
# xout_bfgs,misfout_bfgs = lmbfgs(fh_bfgs!,
140+
# xprior,
141+
# mem=20,
142+
# maxiter=maxiter)
143+
# @show xout[end]
144+
# @show xout_bfgs[end]
145+
146+
ce1 = isapprox(xout[end],[0.0,0.0],atol=1e-1)
147+
return ce1
148+
end

‎test/test_kronlininv.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33

44
###############################################
55

6-
function test2D( )
6+
function test_KLI2D( )
77

88
## 2D problem, so set nx = 1
99
nx = 1
@@ -26,7 +26,7 @@ end
2626

2727
#########################################
2828

29-
function test3D( )
29+
function test_KLI3D( )
3030

3131
## 2D problem, so set nx = 1
3232
nx = 8

0 commit comments

Comments
 (0)