From f33f93491ff6babdc40380c4b79d7a70495de393 Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Mon, 10 May 2021 17:08:10 -0400 Subject: [PATCH 1/7] Added xabsx and new try of odd powers --- src/forward_operators/other.jl | 160 +++++++++++++++++++++++++++++++++ 1 file changed, 160 insertions(+) diff --git a/src/forward_operators/other.jl b/src/forward_operators/other.jl index 38d3d0e..17930ce 100644 --- a/src/forward_operators/other.jl +++ b/src/forward_operators/other.jl @@ -565,3 +565,163 @@ function xexpax(x::MC{N,T}, a::Float64) where {N, T <: RelaxTag} yMC, tp1, tp2 = xexpax_kernel(x, a, intv_xexpax, Inf, Inf) return yMC end + +""" +xabsx + +The function `xabsx` is defined as `xabsx(x) = x*abs(x)`. +""" +xabsx(x::Float64) = x*abs(x) +xabsx_deriv(x::Float64) = 2*abs(x) +xabsx_deriv2(x::Float64) = 2*abs(x)/x #BUT DOESN'T EXIST AT 0 +@inline function cc_xabsx(x::Float64, xL::Float64, xU::Float64) + (xU <= 0.0) && (return xabsx(x), xabsx_deriv(x)) + (0.0 <= xL) && (return dline_seg(^, pow_deriv, x, xL, xU, 2)) + + root_f(x::Float64, xL::Float64, xU::Float64) = (xabsx(xU) - xabsx(x)) - 2*(xU - x)*abs(x) + root_df(x::Float64, xL::Float64, xU::Float64) = xabsx_deriv(x) - xU*xabsx_deriv2(x) + inflection, flag = newton(xL, xL, xU, root_f, root_df, xL, xU) + flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) + if x <= inflection + return dline_seg(^, pow_deriv, x, xL, xU, 2) + else + return xabsx(x), xabsx_deriv(x) + end +end +@inline function cv_xabsx(x::Float64, xL::Float64, xU::Float64) + (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, 2)) + (0.0 <= xL) && (return xabsx(x), xabsx_deriv(x)) + + root_f(x::Float64, xL::Float64, xU::Float64) = (xabsx(x) - xabsx(xL)) - 2*(x - xL)*abs(x) + root_df(x::Float64, xL::Float64, xU::Float64) = -xabsx_deriv(x) + xL*xabsx_deriv2(x) + inflection, flag = newton(xL, xL, xU, root_f, root_df, xL, xU) + flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) + if x <= inflection + return xabsx(x), xabsx_deriv(x) + else + return dline_seg(^, pow_deriv, x, xL, xU, 2) + end +end +function xabsx(x::Interval{Float64}) + xabsx_xL = Interval(x.lo)*abs(Interval(x.lo)) + xabsx_xU = Interval(x.hi)*abs(Interval(x.hi)) + Interval{Float64}(xabsx_xL, xabsx_xU) +end + + +# REWRITTEN ODD POWER FUNCTIONS USING NEWTON'S METHOD +@inline function cv_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} + (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) + (0.0 <= xL) && (return x^n, n*x^(n - 1)) + + root_f(x::Float64, xL::Float64, xU::Float64) = (x^n - xL^n) - (x-xL)*(n)*(x^(n-1)); + root_df(x::Float64, xL::Float64, xU::Float64) = n*x^(n-1) + xL*n*(n-1)*x - n*x^(n-1) - n*(n-1)*x^(n-1); + inflection, flag = newton(xU, xL, xU, root_f, root_df, xL, xU) + flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) + if x <= inflection + return dline_seg(^, pow_deriv, x, xL, xU, n) + else + return x^n, n*x^(n-1) + end + + #val = (xL^n)*(xU - x)/(xU - xL) + (max(0.0, x))^n + #dval = -(xL^n)/(xU - xL) + n*(max(0.0, x))^(n-1) + #return val, dval +end +@inline function cc_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} + (xU <= 0.0) && (return x^n, n*x^(n - 1)) + (0.0 <= xL) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) + + root_f(x::Float64, xL::Float64, xU::Float64) = (xU^n-x^n) - (xU-x)*(n)*(x^(n-1)); + root_df(x::Float64, xL::Float64, xU::Float64) = -n*x^(n-1) - xU*n*(n-1)*x + n*x^(n-1) + n*(n-1)*x^(n-1); + inflection, flag = newton(xL, xL, xU, root_f, root_df, xL, xU) + flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) + if x <= inflection + return x^n, n*x^(n-1) + else + return dline_seg(^, pow_deriv, x, xL, xU, n) + end + + #val = (xU^n)*(x - xL)/(xU - xL) + (min(0.0, x))^n + #dval = (xU^n)/(xU - xL) + n*(min(0.0, x))^(n-1) + #return val, dval +end + +#= +# QUICK REWRITE OF NEWTON'S METHOD FOR FIXED ITERATIONS AND TOLERANCE +function newton(x0::Float64, xL::Float64, xU::Float64, f::Function, df::Function, + envp1::Float64, envp2::Float64) + + dfk = 0.0 + xk = max(xL, min(x0, xU)) + fk::Float64 = f(xk, envp1, envp2) + + for i = 1:5000 + dfk = df(xk, envp1, envp2) + if (abs(fk) < 1e-6) + return (xk, false) + end + (dfk == 0.0) && return (0.0, true) + if (xk == xL && fk/dfk > 0.0) + return (xk, false) + elseif (xk == xU && fk/dfk < 0.0) + return (xk, false) + end + xk = max(xL, min(xU, xk - fk/dfk)) + fk = f(xk, envp1, envp2) + end + + (0.0, true) +end + + + +function xlogx(x::Interval{Float64}) + min_pnt = one(Interval{Float64})/exp(one(Interval{Float64})) + xlogx_xL = Interval(x.lo)*log(Interval(x.lo)) + xlogx_xU = Interval(x.hi)*log(Interval(x.hi)) + xlogx_min_bool = isdisjoint(min_pnt, x) + if isdisjoint(min_pnt, x) + min_val = min(xlogx_xL.lo, xlogx_xU.hi) + else + min_val = -min_pnt.hi + end + Interval{Float64}(min_val, max(xlogx_xL.lo, xlogx_xU.hi)) +end +@inline function xlogx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}} + xL = x.Intv.lo + xU = x.Intv.hi + eps_max = xlogx(xU) > xlogx(xL) ? xU : xL + eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) + midcc, cc_id = mid3(x.cc, x.cv, eps_max) + midcv, cv_id = mid3(x.cc, x.cv, eps_min) + cc, dcc = cc_xlogx(midcc, xL, xU) + cv, dcv = cv_xlogx(midcv, xL, xU) + cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc + cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv + cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad) + return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst) +end +@inline function xlogx_kernel(x::MC{N,Diff}, y::Interval{Float64}) where N + xL = x.Intv.lo + xU = x.Intv.hi + eps_max = xlogx(xU) > xlogx(xL) ? xU : xL + eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) + midcc = mid3v(x.cv, x.cc, eps_max) + midcv = mid3v(x.cv, x.cc, eps_min) + cc, dcc = cc_xlogx(midcc, xL, xU) + cv, dcv = cv_xlogx(midcv, xL, xU) + gcc1, gdcc1 = cc_xlogx(x.cv, xL, xU) + gcv1, gdcv1 = cv_xlogx(x.cv, xL, xU) + gcc2, gdcc2 = cc_xlogx(x.cc, xL, xU) + gcv2, gdcv2 = cv_xlogx(x.cc, xL, xU) + cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad + cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad + return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst) +end +function xlogx(x::MC{N,T}) where {N, T <: RelaxTag} + xlogx_kernel(x, xlogx(x.Intv)) +end + + +=# \ No newline at end of file From d40d86e1246f7a67f7ce7a922b125ca76479530a Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Mon, 10 May 2021 17:35:06 -0400 Subject: [PATCH 2/7] Changes to odd powers --- src/McCormick.jl | 3 ++- src/forward_operators/other.jl | 43 ++++++++++++++++++++++++++++++---- src/forward_operators/power.jl | 2 ++ 3 files changed, 43 insertions(+), 5 deletions(-) diff --git a/src/McCormick.jl b/src/McCormick.jl index 03f329c..9977b39 100644 --- a/src/McCormick.jl +++ b/src/McCormick.jl @@ -71,7 +71,8 @@ export MC, cc, cv, Intv, lo, hi, cc_grad, cv_grad, cnst, +, -, *, /, convert, relu, param_relu, leaky_relu, maxsig, maxtanh, softplus, pentanh, sigmoid, bisigmoid, softsign, gelu, elu, selu, swish1, positive, negative, lower_bnd, upper_bnd, bnd, xlogx, - <, <=, ==, fma, cbrt, abs2, sinpi, cospi, arh, xexpax, trilinear + <, <=, ==, fma, cbrt, abs2, sinpi, cospi, arh, xexpax, trilinear, + xabsx # Export kernel operators export plus_kernel, minus_kernel, mult_kernel, div_kernel, max_kernel, diff --git a/src/forward_operators/other.jl b/src/forward_operators/other.jl index 17930ce..a573c4c 100644 --- a/src/forward_operators/other.jl +++ b/src/forward_operators/other.jl @@ -609,6 +609,41 @@ function xabsx(x::Interval{Float64}) end + +@inline function xlogx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}} + xL = x.Intv.lo + xU = x.Intv.hi + eps_max = xlogx(xU) > xlogx(xL) ? xU : xL + eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) + midcc, cc_id = mid3(x.cc, x.cv, eps_max) + midcv, cv_id = mid3(x.cc, x.cv, eps_min) + cc, dcc = cc_xlogx(midcc, xL, xU) + cv, dcv = cv_xlogx(midcv, xL, xU) + cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc + cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv + cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad) + return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst) +end +@inline function xlogx_kernel(x::MC{N,Diff}, y::Interval{Float64}) where N + xL = x.Intv.lo + xU = x.Intv.hi + eps_max = xlogx(xU) > xlogx(xL) ? xU : xL + eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) + midcc = mid3v(x.cv, x.cc, eps_max) + midcv = mid3v(x.cv, x.cc, eps_min) + cc, dcc = cc_xlogx(midcc, xL, xU) + cv, dcv = cv_xlogx(midcv, xL, xU) + gcc1, gdcc1 = cc_xlogx(x.cv, xL, xU) + gcv1, gdcv1 = cv_xlogx(x.cv, xL, xU) + gcc2, gdcc2 = cc_xlogx(x.cc, xL, xU) + gcv2, gdcv2 = cv_xlogx(x.cc, xL, xU) + cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad + cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad + return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst) +end + + + # REWRITTEN ODD POWER FUNCTIONS USING NEWTON'S METHOD @inline function cv_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) @@ -616,10 +651,10 @@ end root_f(x::Float64, xL::Float64, xU::Float64) = (x^n - xL^n) - (x-xL)*(n)*(x^(n-1)); root_df(x::Float64, xL::Float64, xU::Float64) = n*x^(n-1) + xL*n*(n-1)*x - n*x^(n-1) - n*(n-1)*x^(n-1); - inflection, flag = newton(xU, xL, xU, root_f, root_df, xL, xU) + inflection, flag = newton(xU, 0.0, xU, root_f, root_df, xL, xU) flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) if x <= inflection - return dline_seg(^, pow_deriv, x, xL, xU, n) + return dline_seg(^, pow_deriv, x, xL, inflection, n) else return x^n, n*x^(n-1) end @@ -634,12 +669,12 @@ end root_f(x::Float64, xL::Float64, xU::Float64) = (xU^n-x^n) - (xU-x)*(n)*(x^(n-1)); root_df(x::Float64, xL::Float64, xU::Float64) = -n*x^(n-1) - xU*n*(n-1)*x + n*x^(n-1) + n*(n-1)*x^(n-1); - inflection, flag = newton(xL, xL, xU, root_f, root_df, xL, xU) + inflection, flag = newton(xL, xL, 0.0, root_f, root_df, xL, xU) flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) if x <= inflection return x^n, n*x^(n-1) else - return dline_seg(^, pow_deriv, x, xL, xU, n) + return dline_seg(^, pow_deriv, x, inflection, xU, n) end #val = (xU^n)*(x - xL)/(xU - xL) + (min(0.0, x))^n diff --git a/src/forward_operators/power.jl b/src/forward_operators/power.jl index 0290665..db5f5b9 100644 --- a/src/forward_operators/power.jl +++ b/src/forward_operators/power.jl @@ -99,6 +99,7 @@ end return dline_seg(^, pow_deriv, x, xL, xU, n) end # convex/concave relaxation of odd powers +#= @inline function cv_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) (0.0 <= xL) && (return x^n, n*x^(n - 1)) @@ -113,6 +114,7 @@ end dval = (xU^n)/(xU - xL) + n*(min(0.0, x))^(n-1) return val, dval end +=# @inline function npp_or_pow4(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}} if (x.Intv.hi < 0.0) From 1afd93d975a8a1e18cab042cb35c4de92e4fe3a2 Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Mon, 10 May 2021 17:38:48 -0400 Subject: [PATCH 3/7] Added nonsmooth odd powers --- src/forward_operators/power.jl | 35 ++++++++++++++++++++++++++++++---- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/src/forward_operators/power.jl b/src/forward_operators/power.jl index db5f5b9..ac035fc 100644 --- a/src/forward_operators/power.jl +++ b/src/forward_operators/power.jl @@ -99,7 +99,6 @@ end return dline_seg(^, pow_deriv, x, xL, xU, n) end # convex/concave relaxation of odd powers -#= @inline function cv_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) (0.0 <= xL) && (return x^n, n*x^(n - 1)) @@ -114,7 +113,7 @@ end dval = (xU^n)/(xU - xL) + n*(min(0.0, x))^(n-1) return val, dval end -=# + @inline function npp_or_pow4(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}} if (x.Intv.hi < 0.0) @@ -163,8 +162,8 @@ end @inline function pos_odd(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}} midcc, cc_id = mid3(x.cc, x.cv, x.Intv.hi) midcv, cv_id = mid3(x.cc, x.cv, x.Intv.lo) - cc, dcc = cc_powodd(midcc, x.Intv.lo, x.Intv.hi, c) - cv, dcv = cv_powodd(midcv, x.Intv.lo, x.Intv.hi, c) + cc, dcc = cc_powodd_ns(midcc, x.Intv.lo, x.Intv.hi, c) + cv, dcv = cv_powodd_ns(midcv, x.Intv.lo, x.Intv.hi, c) cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad) @@ -183,6 +182,34 @@ end cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst) end +@inline function cv_powodd_ns(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} + (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) + (0.0 <= xL) && (return x^n, n*x^(n - 1)) + + root_f(x::Float64, xL::Float64, xU::Float64) = (x^n - xL^n) - (x-xL)*(n)*(x^(n-1)); + root_df(x::Float64, xL::Float64, xU::Float64) = n*x^(n-1) + xL*n*(n-1)*x - n*x^(n-1) - n*(n-1)*x^(n-1); + inflection, flag = newton(xU, 0.0, xU, root_f, root_df, xL, xU) + flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) + if x <= inflection + return dline_seg(^, pow_deriv, x, xL, inflection, n) + else + return x^n, n*x^(n-1) + end +end +@inline function cc_powodd_ns(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} + (xU <= 0.0) && (return x^n, n*x^(n - 1)) + (0.0 <= xL) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) + + root_f(x::Float64, xL::Float64, xU::Float64) = (xU^n-x^n) - (xU-x)*(n)*(x^(n-1)); + root_df(x::Float64, xL::Float64, xU::Float64) = -n*x^(n-1) - xU*n*(n-1)*x + n*x^(n-1) + n*(n-1)*x^(n-1); + inflection, flag = newton(xL, xL, 0.0, root_f, root_df, xL, xU) + flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) + if x <= inflection + return x^n, n*x^(n-1) + else + return dline_seg(^, pow_deriv, x, inflection, xU, n) + end +end # neg_powneg_odd computes the McComrick relaxation of x^c where x < 0.0 and c is odd @inline function cv_neg_powneg_odd(x::Float64, xL::Float64, xU::Float64, c::Z) where Z<:Integer From 21e07aa696048e356c1a5e776db0bbe743bba7da Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Wed, 12 May 2021 15:08:19 -0400 Subject: [PATCH 4/7] Fixes to xabsx and addition of tests --- src/forward_operators/other.jl | 205 +++++++-------------------------- test/forward_mccormick.jl | 28 +++++ 2 files changed, 69 insertions(+), 164 deletions(-) diff --git a/src/forward_operators/other.jl b/src/forward_operators/other.jl index a573c4c..859bfb9 100644 --- a/src/forward_operators/other.jl +++ b/src/forward_operators/other.jl @@ -576,187 +576,64 @@ xabsx_deriv(x::Float64) = 2*abs(x) xabsx_deriv2(x::Float64) = 2*abs(x)/x #BUT DOESN'T EXIST AT 0 @inline function cc_xabsx(x::Float64, xL::Float64, xU::Float64) (xU <= 0.0) && (return xabsx(x), xabsx_deriv(x)) - (0.0 <= xL) && (return dline_seg(^, pow_deriv, x, xL, xU, 2)) + (0.0 <= xL) && (return dline_seg(xabsx, xabsx_deriv, x, xL, xU)) - root_f(x::Float64, xL::Float64, xU::Float64) = (xabsx(xU) - xabsx(x)) - 2*(xU - x)*abs(x) - root_df(x::Float64, xL::Float64, xU::Float64) = xabsx_deriv(x) - xU*xabsx_deriv2(x) - inflection, flag = newton(xL, xL, xU, root_f, root_df, xL, xU) + root_f(x, xL, xU) = (xabsx(xU) - xabsx(x)) - 2*(xU - x)*abs(x) + root_df(x, xL, xU) = xabsx_deriv(x) - xU*xabsx_deriv2(x) + inflection, flag = newton(xL, xL, 0.0, root_f, root_df, xL, xU) flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) if x <= inflection - return dline_seg(^, pow_deriv, x, xL, xU, 2) - else return xabsx(x), xabsx_deriv(x) + else + return dline_seg(xabsx, xabsx_deriv, x, inflection, xU) end end @inline function cv_xabsx(x::Float64, xL::Float64, xU::Float64) - (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, 2)) + (xU <= 0.0) && (return dline_seg(xabsx, xabsx_deriv, x, xL, xU)) (0.0 <= xL) && (return xabsx(x), xabsx_deriv(x)) - root_f(x::Float64, xL::Float64, xU::Float64) = (xabsx(x) - xabsx(xL)) - 2*(x - xL)*abs(x) - root_df(x::Float64, xL::Float64, xU::Float64) = -xabsx_deriv(x) + xL*xabsx_deriv2(x) - inflection, flag = newton(xL, xL, xU, root_f, root_df, xL, xU) + root_f(x, xL, xU) = (xabsx(x) - xabsx(xL)) - 2*(x - xL)*abs(x) + root_df(x, xL, xU) = -xabsx_deriv(x) + xL*xabsx_deriv2(x) + inflection, flag = newton(xU, 0.0, xU, root_f, root_df, xL, xU) flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) if x <= inflection - return xabsx(x), xabsx_deriv(x) + return dline_seg(xabsx, xabsx_deriv, x, xL, inflection) else - return dline_seg(^, pow_deriv, x, xL, xU, 2) + return xabsx(x), xabsx_deriv(x) end end function xabsx(x::Interval{Float64}) xabsx_xL = Interval(x.lo)*abs(Interval(x.lo)) xabsx_xU = Interval(x.hi)*abs(Interval(x.hi)) - Interval{Float64}(xabsx_xL, xabsx_xU) -end - - - -@inline function xlogx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}} - xL = x.Intv.lo - xU = x.Intv.hi - eps_max = xlogx(xU) > xlogx(xL) ? xU : xL - eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) - midcc, cc_id = mid3(x.cc, x.cv, eps_max) - midcv, cv_id = mid3(x.cc, x.cv, eps_min) - cc, dcc = cc_xlogx(midcc, xL, xU) - cv, dcv = cv_xlogx(midcv, xL, xU) - cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc - cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv - cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad) - return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst) -end -@inline function xlogx_kernel(x::MC{N,Diff}, y::Interval{Float64}) where N - xL = x.Intv.lo - xU = x.Intv.hi - eps_max = xlogx(xU) > xlogx(xL) ? xU : xL - eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) - midcc = mid3v(x.cv, x.cc, eps_max) - midcv = mid3v(x.cv, x.cc, eps_min) - cc, dcc = cc_xlogx(midcc, xL, xU) - cv, dcv = cv_xlogx(midcv, xL, xU) - gcc1, gdcc1 = cc_xlogx(x.cv, xL, xU) - gcv1, gdcv1 = cv_xlogx(x.cv, xL, xU) - gcc2, gdcc2 = cc_xlogx(x.cc, xL, xU) - gcv2, gdcv2 = cv_xlogx(x.cc, xL, xU) - cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad - cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad - return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst) + Interval{Float64}(xabsx_xL.lo, xabsx_xU.hi) end - - - -# REWRITTEN ODD POWER FUNCTIONS USING NEWTON'S METHOD -@inline function cv_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} - (xU <= 0.0) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) - (0.0 <= xL) && (return x^n, n*x^(n - 1)) - - root_f(x::Float64, xL::Float64, xU::Float64) = (x^n - xL^n) - (x-xL)*(n)*(x^(n-1)); - root_df(x::Float64, xL::Float64, xU::Float64) = n*x^(n-1) + xL*n*(n-1)*x - n*x^(n-1) - n*(n-1)*x^(n-1); - inflection, flag = newton(xU, 0.0, xU, root_f, root_df, xL, xU) - flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) - if x <= inflection - return dline_seg(^, pow_deriv, x, xL, inflection, n) - else - return x^n, n*x^(n-1) - end - - #val = (xL^n)*(xU - x)/(xU - xL) + (max(0.0, x))^n - #dval = -(xL^n)/(xU - xL) + n*(max(0.0, x))^(n-1) - #return val, dval -end -@inline function cc_powodd(x::Float64, xL::Float64, xU::Float64, n::Z) where {Z <: Integer} - (xU <= 0.0) && (return x^n, n*x^(n - 1)) - (0.0 <= xL) && (return dline_seg(^, pow_deriv, x, xL, xU, n)) - - root_f(x::Float64, xL::Float64, xU::Float64) = (xU^n-x^n) - (xU-x)*(n)*(x^(n-1)); - root_df(x::Float64, xL::Float64, xU::Float64) = -n*x^(n-1) - xU*n*(n-1)*x + n*x^(n-1) + n*(n-1)*x^(n-1); - inflection, flag = newton(xL, xL, 0.0, root_f, root_df, xL, xU) - flag && (inflection = golden_section(xL, xU, root_f, xL, xU)) - if x <= inflection - return x^n, n*x^(n-1) - else - return dline_seg(^, pow_deriv, x, inflection, xU, n) - end - - #val = (xU^n)*(x - xL)/(xU - xL) + (min(0.0, x))^n - #dval = (xU^n)/(xU - xL) + n*(min(0.0, x))^(n-1) - #return val, dval -end - -#= -# QUICK REWRITE OF NEWTON'S METHOD FOR FIXED ITERATIONS AND TOLERANCE -function newton(x0::Float64, xL::Float64, xU::Float64, f::Function, df::Function, - envp1::Float64, envp2::Float64) - - dfk = 0.0 - xk = max(xL, min(x0, xU)) - fk::Float64 = f(xk, envp1, envp2) - - for i = 1:5000 - dfk = df(xk, envp1, envp2) - if (abs(fk) < 1e-6) - return (xk, false) - end - (dfk == 0.0) && return (0.0, true) - if (xk == xL && fk/dfk > 0.0) - return (xk, false) - elseif (xk == xU && fk/dfk < 0.0) - return (xk, false) - end - xk = max(xL, min(xU, xk - fk/dfk)) - fk = f(xk, envp1, envp2) - end - - (0.0, true) -end - - - -function xlogx(x::Interval{Float64}) - min_pnt = one(Interval{Float64})/exp(one(Interval{Float64})) - xlogx_xL = Interval(x.lo)*log(Interval(x.lo)) - xlogx_xU = Interval(x.hi)*log(Interval(x.hi)) - xlogx_min_bool = isdisjoint(min_pnt, x) - if isdisjoint(min_pnt, x) - min_val = min(xlogx_xL.lo, xlogx_xU.hi) - else - min_val = -min_pnt.hi - end - Interval{Float64}(min_val, max(xlogx_xL.lo, xlogx_xU.hi)) -end -@inline function xlogx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}} - xL = x.Intv.lo - xU = x.Intv.hi - eps_max = xlogx(xU) > xlogx(xL) ? xU : xL - eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) - midcc, cc_id = mid3(x.cc, x.cv, eps_max) - midcv, cv_id = mid3(x.cc, x.cv, eps_min) - cc, dcc = cc_xlogx(midcc, xL, xU) - cv, dcv = cv_xlogx(midcv, xL, xU) - cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc - cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv - cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad) - return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst) -end -@inline function xlogx_kernel(x::MC{N,Diff}, y::Interval{Float64}) where N - xL = x.Intv.lo - xU = x.Intv.hi - eps_max = xlogx(xU) > xlogx(xL) ? xU : xL - eps_min = in(-1.0/exp(1.0), x) ? -1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ? xL : xU) - midcc = mid3v(x.cv, x.cc, eps_max) - midcv = mid3v(x.cv, x.cc, eps_min) - cc, dcc = cc_xlogx(midcc, xL, xU) - cv, dcv = cv_xlogx(midcv, xL, xU) - gcc1, gdcc1 = cc_xlogx(x.cv, xL, xU) - gcv1, gdcv1 = cv_xlogx(x.cv, xL, xU) - gcc2, gdcc2 = cc_xlogx(x.cc, xL, xU) - gcv2, gdcv2 = cv_xlogx(x.cc, xL, xU) - cv_grad = max(0.0, gdcv1)*x.cv_grad + min(0.0, gdcv2)*x.cc_grad - cc_grad = min(0.0, gdcc1)*x.cv_grad + max(0.0, gdcc2)*x.cc_grad - return MC{N,Diff}(cv, cc, y, cv_grad, cc_grad, x.cnst) +@inline function xabsx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}} + xL = x.Intv.lo + xU = x.Intv.hi + eps_max = xU + eps_min = xL + midcc, cc_id = mid3(x.cc, x.cv, eps_max) + midcv, cv_id = mid3(x.cc, x.cv, eps_min) + cc, dcc = cc_xabsx(midcc, xL, xU) + cv, dcv = cv_xabsx(midcv, xL, xU) + cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc + cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv + cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad) + return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst) +end +@inline function xabsx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Diff} + xL = x.Intv.lo + xU = x.Intv.hi + eps_max = xU + eps_min = xL + midcc, cc_id = mid3(x.cc, x.cv, eps_max) + midcv, cv_id = mid3(x.cc, x.cv, eps_min) + cc, dcc = cc_xabsx(midcc, xL, xU) + cv, dcv = cv_xabsx(midcv, xL, xU) + cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc + cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv + return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst) end -function xlogx(x::MC{N,T}) where {N, T <: RelaxTag} - xlogx_kernel(x, xlogx(x.Intv)) +function xabsx(x::MC{N,T}) where {N, T <: RelaxTag} + xabsx_kernel(x, xabsx(x.Intv)) end - - -=# \ No newline at end of file diff --git a/test/forward_mccormick.jl b/test/forward_mccormick.jl index 00e1a68..a71d0bc 100644 --- a/test/forward_mccormick.jl +++ b/test/forward_mccormick.jl @@ -554,6 +554,34 @@ end @test check_vs_ref1(atanh, x_atanh_z2, yref_atanh_d1_z2, mctol) @test check_vs_ref1(atanh, x_atanh_z1_ns, yref_atanh_ns, mctol) + ##### x Times Abs(x) ##### + x_xabsx_p = MC{2,Diff}(0.3, 0.3, Interval{Float64}(0.1,0.7), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) + x_xabsx_n = MC{2,Diff}(-0.3, -0.3, Interval{Float64}(-0.7,-0.1), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) + x_xabsx_z1 = MC{2,Diff}(0.2, 0.2, Interval{Float64}(-0.7,0.5), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) + x_xabsx_z2 = MC{2,Diff}(0.1, 0.1, Interval{Float64}(-0.5,0.7), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) + + yref_xabsx_d1_p = MC{2, Diff}(0.09, 0.16999999999999998, Interval{Float64}(0.01, 0.49), @SVector[0.6, 0.0], @SVector[0.7999999999999999, 0.0], false) + yref_xabsx_d1_n = MC{2, Diff}(-0.16999999999999998, -0.09, Interval{Float64}(-0.49, -0.01), @SVector[0.7999999999999999, 0.0], @SVector[0.6, 0.0], false) + yref_xabsx_d1_z1 = MC{2, Diff}(0.03190908859009977, 0.1257359312880715, Interval{Float64}(-0.49, 0.25), @SVector[0.579898987322333, 0.0], @SVector[0.4142135623730951, 0.0], false) + yref_xabsx_d1_z2 = MC{2, Diff}(-0.001471862576142969, 0.14206060760660016, Interval{Float64}(-0.25, 0.49), @SVector[0.4142135623730951, 0.0], @SVector[0.579898987322333, 0.0], false) + + @test check_vs_ref1(xabsx, x_xabsx_p, yref_xabsx_d1_p, mctol) + @test check_vs_ref1(xabsx, x_xabsx_n, yref_xabsx_d1_n, mctol) + @test check_vs_ref1(xabsx, x_xabsx_z1, yref_xabsx_d1_z1, mctol) + @test check_vs_ref1(xabsx, x_xabsx_z2, yref_xabsx_d1_z2, mctol) + + #= + "Convex relaxation" + cv::Float64 + "Concave relaxation" + cc::Float64 + "Interval bounds" + Intv::Interval{Float64} + "(Sub)gradient of convex relaxation" + cv_grad::SVector{N,Float64} + "(Sub)gradient of concave relaxation" + cc_grad::SVector{N,Float64} +=# # CONVERSION X = MC{2,NS}(4.5, 4.5, Interval{Float64}(-3.0,8.0), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) X1 = convert(MC{2,NS}, 1) From 6714b2dbe832cd754d6f7a3b3f2138df24640999 Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Wed, 12 May 2021 15:11:13 -0400 Subject: [PATCH 5/7] Minor fix --- src/forward_operators/power.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/forward_operators/power.jl b/src/forward_operators/power.jl index ac035fc..9f210ba 100644 --- a/src/forward_operators/power.jl +++ b/src/forward_operators/power.jl @@ -114,7 +114,6 @@ end return val, dval end - @inline function npp_or_pow4(x::MC{N,T}, c::Z, y::Interval{Float64}) where {N, Z<:Integer, T<:Union{NS,MV}} if (x.Intv.hi < 0.0) eps_min = x.Intv.hi From 301ade922305d2b2b532018264a36500bafa283f Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Wed, 12 May 2021 15:13:10 -0400 Subject: [PATCH 6/7] Removing extraneous comment --- test/forward_mccormick.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/test/forward_mccormick.jl b/test/forward_mccormick.jl index a71d0bc..7861dbb 100644 --- a/test/forward_mccormick.jl +++ b/test/forward_mccormick.jl @@ -570,18 +570,6 @@ end @test check_vs_ref1(xabsx, x_xabsx_z1, yref_xabsx_d1_z1, mctol) @test check_vs_ref1(xabsx, x_xabsx_z2, yref_xabsx_d1_z2, mctol) - #= - "Convex relaxation" - cv::Float64 - "Concave relaxation" - cc::Float64 - "Interval bounds" - Intv::Interval{Float64} - "(Sub)gradient of convex relaxation" - cv_grad::SVector{N,Float64} - "(Sub)gradient of concave relaxation" - cc_grad::SVector{N,Float64} -=# # CONVERSION X = MC{2,NS}(4.5, 4.5, Interval{Float64}(-3.0,8.0), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) X1 = convert(MC{2,NS}, 1) From 81cd02c711d4660fe7ac6d26c9fba664eeefa00a Mon Sep 17 00:00:00 2001 From: Robert Gottlieb Date: Wed, 12 May 2021 15:23:12 -0400 Subject: [PATCH 7/7] Change to xabsx test --- test/forward_mccormick.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/forward_mccormick.jl b/test/forward_mccormick.jl index 7861dbb..a497088 100644 --- a/test/forward_mccormick.jl +++ b/test/forward_mccormick.jl @@ -558,12 +558,12 @@ end x_xabsx_p = MC{2,Diff}(0.3, 0.3, Interval{Float64}(0.1,0.7), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) x_xabsx_n = MC{2,Diff}(-0.3, -0.3, Interval{Float64}(-0.7,-0.1), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) x_xabsx_z1 = MC{2,Diff}(0.2, 0.2, Interval{Float64}(-0.7,0.5), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) - x_xabsx_z2 = MC{2,Diff}(0.1, 0.1, Interval{Float64}(-0.5,0.7), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) + x_xabsx_z2 = MC{2,Diff}(0.4, 0.4, Interval{Float64}(-0.7,0.5), seed_gradient(1,Val(2)), seed_gradient(1,Val(2)), false) yref_xabsx_d1_p = MC{2, Diff}(0.09, 0.16999999999999998, Interval{Float64}(0.01, 0.49), @SVector[0.6, 0.0], @SVector[0.7999999999999999, 0.0], false) yref_xabsx_d1_n = MC{2, Diff}(-0.16999999999999998, -0.09, Interval{Float64}(-0.49, -0.01), @SVector[0.7999999999999999, 0.0], @SVector[0.6, 0.0], false) yref_xabsx_d1_z1 = MC{2, Diff}(0.03190908859009977, 0.1257359312880715, Interval{Float64}(-0.49, 0.25), @SVector[0.579898987322333, 0.0], @SVector[0.4142135623730951, 0.0], false) - yref_xabsx_d1_z2 = MC{2, Diff}(-0.001471862576142969, 0.14206060760660016, Interval{Float64}(-0.25, 0.49), @SVector[0.4142135623730951, 0.0], @SVector[0.579898987322333, 0.0], false) + yref_xabsx_d1_z2 = MC{2, Diff}(0.16000000000000003, 0.20857864376269047, Interval{Float64}(-0.49, 0.25), @SVector[0.8, 0.0], @SVector[0.4142135623730951, 0.0], false) @test check_vs_ref1(xabsx, x_xabsx_p, yref_xabsx_d1_p, mctol) @test check_vs_ref1(xabsx, x_xabsx_n, yref_xabsx_d1_n, mctol)