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.

Sign up for GitHub

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

Refactor Implicit Function Routine and Fix Bug #52

Merged
merged 2 commits into from
May 20, 2022
Merged

Refactor Implicit Function Routine and Fix Bug #52

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ae68d8b
Fix bug in implicit relaxation routine and add documentation
mewilhel May 20, 2022
1e9cc7f
Merge remote-tracking branch 'mewilhel/master'
mewilhel May 20, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "McCormick"
uuid = "53c679d3-6890-5091-8386-c291e8c8aaa1"
authors = ["Matthew Wilhelm <[email protected]>"]
version = "0.13.0"
version = "0.13.1"

[deps]
DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
Expand All @@ -15,6 +15,7 @@ NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

[compat]
DiffRules = "1.5"
Expand All @@ -27,6 +28,7 @@ NNlib = "~0.7"
Reexport = "~0.2, ~1"
StaticArrays = "~1.2"
SpecialFunctions = "~1, ~2"
UnPack = "~1"
julia = "~1.6, ~1.7"

[extras]
Expand Down
10 changes: 9 additions & 1 deletion src/McCormick.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ __precompile__()

module McCormick

using DocStringExtensions, LinearAlgebra
using DocStringExtensions, LinearAlgebra, UnPack
using DiffRules: diffrule
using StaticArrays: @SVector, SVector, zeros, ones
using ForwardDiff: Dual, Partials
Expand Down Expand Up @@ -618,6 +618,14 @@ function check_relaxation_error!(x::MC{N,T}) where {N, T<:RelaxTag}
return nothing
end

function create_vec_of_mutable(::Type{T}, f, k::Int) where T
z = T[]
for i = 1:k
push!(z, f())
end
z
end

include("forward_operators/forward.jl")
include("implicit_routines/implicit.jl")

Expand Down
78 changes: 38 additions & 40 deletions src/implicit_routines/contract.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,44 +71,43 @@ $(TYPEDSIGNATURES)

Applies the Gauss-Siedel variant of the Newton type contractor.
"""
function contract!(t::NewtonGS, d::MCCallback{FH,FJ,C,PRE,N,T}, k::Int, b::Bool) where {FH <: Function,
FJ <: Function,
C, PRE, N,
T<:RelaxTag}
function contract!(t::NewtonGS, d::MCCallback{FH,FJ,C,PRE,N,T}, k::Int, b::Bool) where {FH, FJ, C, PRE, N, T<:RelaxTag}
@unpack H, J0, J, p_mc, pref_mc, nx, x0_mc, x_mc, xz0, z_mc, use_apriori = d

S = zero(MC{N,T})
@. d.x0_mc = d.x_mc
for i = 1:d.nx
@. x0_mc = x_mc
for i = 1:nx
S = zero(MC{N,T})
for j = 1:d.nx
if (i !== j)
xv = @inbounds d.J[i,j]
yv = @inbounds d.x_mc[j] - d.z_mc[j]
for j = 1:nx
if i != j
xv = @inbounds J[i,j]
yv = @inbounds x_mc[j] - z_mc[j]
zv = xv*yv
if d.use_apriori
if use_apriori
if b
d.J0[k][i,j] = xv
d.xz0[k][i,j] = yv
J0[k][i,j] = xv
xz0[k][i,j] = yv
end
xr = d.J0[k][i,j]
yr = d.xz0[k][i,j]
u1max, u2max, v1nmax, v2nmax = estimator_extrema(xr, yr, d.p_ref)
xr = J0[k][i,j]
yr = xz0[k][i,j]
u1max, u2max, v1nmax, v2nmax = estimator_extrema(xr, yr, p_ref)
wIntv = zv.Intv
if (u1max < xv.Intv.hi) || (u2max < yv.Intv.hi)
u1cv, u2cv, u1cvg, u2cvg = estimator_under(xr, yr, d.p_mc, d.p_ref)
u1cv, u2cv, u1cvg, u2cvg = estimator_under(xr, yr, p_mc, p_ref)
za_l = mult_apriori_kernel(xv, yv, wIntv, u1cv, u2cv, u1max, u2max, u1cvg, u2cvg)
zv = zv ∩ za_l
end
if (v1nmax > -xv.Intv.lo) || (v2nmax > -yv.Intv.lo)
v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_over(xr, yr, d.p_mc, d.p_ref)
v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_over(xr, yr, p_mc, p_ref)
za_u = mult_apriori_kernel(-xv, -yv, wIntv, v1ccn, v2ccn, v1nmax, v2nmax, v1ccgn, v2ccgn)
zv = zv ∩ za_u
end
end
S += zv
end
end
@inbounds d.x_mc[i] = d.z_mc[i] - (d.H[i] + S)*McCormick.inv1(d.J[i,i], 1.0/d.J[i,i].Intv)
@inbounds d.x_mc[i] = final_cut(d.x_mc[i], d.x0_mc[i])
@inbounds x_mc[i] = z_mc[i] - (H[i] + S)*McCormick.inv1(J[i,i], 1.0/J[i,i].Intv)
@inbounds x_mc[i] = final_cut(x_mc[i], x0_mc[i])
end
return
end
Expand All @@ -118,42 +117,41 @@ $(TYPEDSIGNATURES)

Applies the componentwise variant of the Krawczyk type contractor.
"""
function contract!(t::KrawczykCW, d::MCCallback{FH,FJ,C,PRE,N,T}, k::Int, b::Bool) where {FH <: Function,
FJ <: Function,
C, PRE, N,
T<:RelaxTag}
S::MC{N,T} = zero(MC{N,T})
@. d.x0_mc = d.x_mc
for i=1:d.nx
function contract!(t::KrawczykCW, d::MCCallback{FH,FJ,C,PRE,N,T}, k::Int, b::Bool) where {FH, FJ, C, PRE, N, T<:RelaxTag}
@unpack H, J0, J, p_mc, pref_mc, nx, x0_mc, x_mc, xz0, z_mc, use_apriori = d

S = zero(MC{N,T})
@. x0_mc = x_mc
for i = 1:nx
S = zero(MC{N,T})
for j=1:d.nx
xv = (i !== j) ? -d.J[i,j] : (one(MC{N,T}) - d.J[i,j])
yv = @inbounds d.x_mc[j] - d.z_mc[j]
for j = 1:d.nx
xv = (i != j) ? -J[i,j] : (one(MC{N,T}) - J[i,j])
yv = @inbounds x_mc[j] - z_mc[j]
zv = xv*yv
if d.use_apriori
if use_apriori
if b
d.J0[k][i,j] = xv
d.xz0[k][i,j] = yv
J0[k][i,j] = xv
xz0[k][i,j] = yv
end
xr = d.J0[k][i,j]
yr = d.xz0[k][i,j]
u1max, u2max, v1nmax, v2nmax = estimator_extrema(xr, yr, d.p_ref)
xr = J0[k][i,j]
yr = xz0[k][i,j]
u1max, u2max, v1nmax, v2nmax = estimator_extrema(xr, yr, p_ref)
wIntv = zv.Intv
if (u1max < xv.Intv.hi) || (u2max < yv.Intv.hi)
u1cv, u2cv, u1cvg, u2cvg = estimator_under(xr, yr, d.p_mc, d.p_ref)
u1cv, u2cv, u1cvg, u2cvg = estimator_under(xr, yr, p_mc, p_ref)
za_l = mult_apriori_kernel(xv, yv, wIntv, u1cv, u2cv, u1max, u2max, u1cvg, u2cvg)
zv = zv ∩ za_l
end
if (v1nmax > -xv.Intv.lo) || (v2nmax > -yv.Intv.lo)
v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_over(xr, yr, d.p_mc, d.p_ref)
v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_over(xr, yr, p_mc, p_ref)
za_u = mult_apriori_kernel(-xv, -yv, wIntv, v1ccn, v2ccn, v1nmax, v2nmax, v1ccgn, v2ccgn)
zv = zv ∩ za_u
end
end
S += zv
end
@inbounds d.x_mc[i] = d.z_mc[i] - d.H[i] + S
@inbounds d.x_mc[i] = final_cut(d.x_mc[i], d.x0_mc[i])
@inbounds x_mc[i] = z_mc[i] - H[i] + S
@inbounds x_mc[i] = final_cut(x_mc[i], x0_mc[i])
end
return
end
Loading