Skip to content

Commit 99865f7

Browse files
RXGottliebRXGottlieb
committed
Added documentation and unit tests, fixed bugs
Fixed some errors in negative-valued multiplication, and added unit testing to cover the different cases. Also added some unit testing for addition/subtraction, though it could be simplified.
1 parent e0a0424 commit 99865f7

File tree

6 files changed

+389
-221
lines changed

6 files changed

+389
-221
lines changed

‎README.md

Lines changed: 115 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,106 @@
11
# SourceCodeMcCormick.jl
2+
3+
This package is an experimental approach to use source-code transformation to apply McCormick relaxations
4+
to symbolic functions for use in deterministic global optimization. While packages like `McCormick.jl`
5+
take set-valued McCormick objects and utilize McCormick relaxation rules to overload standard math operations,
6+
`SourceCodeMcCormick.jl` aims to interpret symbolic expressions, apply McCormick transformations, and
7+
return a new symbolic expression representing those relaxations. This functionality is designed to be
8+
used for both algebraic and dynamic systems.
9+
210
Experimental Approach to McCormick Relaxation Source-Code Transformation for Differential Inequalities
311

4-
The purpose of this package is to transform a `ModelingToolkit` ODE system with factorable equations into a new `ModelingToolkit` ODE system with interval or McCormick relaxations applied. E.g.:
12+
## Algebraic Systems
13+
14+
For a given algebraic equation or system of equations, `SourceCodeMcCormick` is designed to provide
15+
symbolic transformations that represent the lower/upper bounds and convex/concave relaxations of the
16+
provided equation(s). Most notably, `SourceCodeMcCormick` uses this symbolic transformation to
17+
generate "evaluation functions" which, for a given expression, return the lower/upper bounds and
18+
convex/concave relaxations of an expression. E.g.:
19+
20+
```
21+
using SourceCodeMcCormick, Symbolics
22+
23+
@variables x
24+
to_compute = x*(15+x)^2
25+
x_lo_eval, x_hi_eval, x_cv_eval, x_cc_eval, order = all_evaluators(to_compute)
26+
```
27+
28+
Here, the outputs marked `_eval` are the evaluation functions for the lower bound (`lo`), upper
29+
bound (`hi`), convex underestimator (`cv`), and concave overestimator (`cc`) of the `to_compute`
30+
expression. The inputs to each of these functions are described by the `order` vector, which
31+
in this case is `[x_cc, x_cv, x_hi, x_lo]`.
32+
33+
There are several important benefits of using these functions. First, they are fast. Normally,
34+
applying McCormick relaxations takes some time as the relaxation depends on the provided bounds
35+
and relaxation values of the input(s), but also on the form of the expression itself. Packages
36+
like `McCormick.jl` use control flow to quickly determine the form of the relaxations and provide
37+
McCormick objects of the results, but because the forms of the expressions are not known *a priori*,
38+
this control flow takes time to process. In contrast, `SourceCodeMcCormick` effectively hard-codes
39+
the relaxations for a specific, pre-defined function into the evaluation functions, making it
40+
inherently faster. For example:
41+
42+
```
43+
using BenchmarkTools, McCormick
44+
45+
xMC = MC{1, NS}(2.5, Interval(-1.0, 4.0), 1)
46+
47+
@btime xMC*(15+xMC)^2
48+
# 228.381 ns (15 allocations: 384 bytes)
49+
# MC{1, NS}(683.5, 952.0, [-361, 1444], [501.0], [328.0], false)
50+
51+
@btime x_cv_eval(2.5, 2.5, 4.0, -1.0)
52+
# 30.382 ns (1 allocation: 16 bytes)
53+
# 683.5
54+
```
55+
56+
This is not an *entirely* fair example because the operation with xMC is simultaneously calculating
57+
lower/upper bounds, both relaxations, and (sub)gradients of the convex/concave relaxations, but the
58+
`SourceCodeMcCormick` result is just under an order of magnitude faster. As the expression inceases
59+
in complexity, this speedup becomes even more apparent.
60+
61+
Second, these evaluation functions are compatible with `CUDA.jl` in that they are broadcastable
62+
over `CuArray`s. Depending on the GPU, number of evaluations, and complexity of the function,
63+
this can dramatically decrease the time to compute the numerical values of function bounds and
64+
relaxations. E.g.:
65+
66+
```
67+
using CUDA
68+
69+
# Using McCormick.jl
70+
xMC_array = MC{1,NS}.(rand(10000), Interval.(zeros(10000), ones(10000)), ones(Int, 10000))
71+
@btime xMC_array.*(15 .+ xMC_array).^2
72+
# 1.616 ms (120012 allocations: 2.37 MiB)
73+
74+
# Using SourceCodeMcCormick.jl, broadcast using CPU
75+
xcc = rand(10000)
76+
xcv = copy(xcc)
77+
xhi = ones(10000)
78+
xlo = zeros(10000)
79+
@btime x_cv_eval.(xcc, xcv, xhi, xlo)
80+
# 100.100 μs (4 allocations: 78.27 KiB)
81+
82+
# Using SourceCodeMcCormick.jl and CUDA.jl, broadcast using GPU
83+
xcc_GPU = cu(xcc)
84+
xcv_GPU = cu(xcv)
85+
xhi_GPU = cu(xhi)
86+
xlo_GPU = cu(xlo)
87+
@btime x_cv_eval.(xcc_GPU, xcv_GPU, xhi_GPU, xlo_GPU)
88+
# 6.575 μs (33 allocations: 2.34 KiB)
89+
```
90+
91+
92+
## Dynamic Systems
93+
94+
For dynamic systems, `SourceCodeMcCormick` was built with a differential inequalities
95+
approach in mind where the relaxations of derivatives are calculated in advance and
96+
the resulting (larger) differential equation system, with explicit definitions of
97+
the relaxations of derivatives, can be solved. For algebraic systems, the main
98+
product of this package is the broadcastable evaluation functions. For dynamic
99+
systems, this package follows the same idea as in algebraic systems but stops at
100+
the symbolic representations of relaxations. This functionality is designed to work
101+
with a `ModelingToolkit`-type `ODESystem` with factorable equations--`SourceCodeMcCormick`
102+
will take such a system and return a new `ODESystem` with expanded equations to
103+
provide interval extensions and (if desired) McCormick relaxations. E.g.:
5104

6105
```
7106
using SourceCodeMcCormick, ModelingToolkit
@@ -40,4 +139,18 @@ Differential(t)(x_2_cv(t)) ~ p_2_cv + x_2_cv(t)
40139
Differential(t)(x_2_cc(t)) ~ p_2_cc + x_2_cc(t)
41140
```
42141

43-
where `x_lo < x_cv < x < x_cc < x_hi`. This new system of ODEs is generated using GPU-compatible language--i.e., any decision points in the form of the resulting equation based on some terms being positive or negative are handled by IfElse.ifelse statements and/or min/max evaluations. By using only these types of expressions, multiple trajectories of the resulting ODE system can be solved simultaneously on a GPU, such as by using `DiffEqGPU` in the SciML ecosystem.
142+
where `x_lo < x_cv < x < x_cc < x_hi`. Only addition is shown in this example, as other operations
143+
can appear very expansive, but the same operations available for algebraic systems are available
144+
for dynamic systems as well. As with the algebraic evaluation functions, equations created by
145+
`SourceCodeMcCormick` are GPU-ready--multiple trajectories of the resulting ODE system at
146+
different points and with different state/parameter bounds can be solved simultaneously using
147+
an `EnsembleProblem` in the SciML ecosystem, and GPU hardware can be applied for these solves
148+
using `DiffEqGPU.jl`.
149+
150+
## Limitations
151+
152+
Currently, as proof-of-concept, `SourceCodeMcCormick` can only handle functions with
153+
addition (+), subtraction (-), multiplication (\*), powers of 2 (^2), natural base
154+
exponentials (exp), and minimum/maximum (min/max) expressions. Future work will include
155+
adding other operations found in `McCormick.jl`.
156+

‎src/SourceCodeMcCormick.jl

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -15,14 +15,15 @@ abstract type AbstractTransform end
1515
function transform_rule end
1616

1717

18+
include(joinpath(@__DIR__, "interval", "interval.jl"))
19+
include(joinpath(@__DIR__, "relaxation", "relaxation.jl"))
20+
include(joinpath(@__DIR__, "transform", "transform.jl"))
21+
1822
export McCormickIntervalTransform
1923

2024
export apply_transform, extract_terms, genvar, genparam, get_name,
2125
factor!, binarize!, pull_vars, shrink_eqs, convex_evaluator,
2226
all_evaluators
2327

24-
include(joinpath(@__DIR__, "interval", "interval.jl"))
25-
include(joinpath(@__DIR__, "relaxation", "relaxation.jl"))
26-
include(joinpath(@__DIR__, "transform", "transform.jl"))
27-
28+
2829
end

‎src/interval/rules.jl

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -36,22 +36,22 @@ function transform_rule(::IntervalTransform, ::typeof(-), zL, zU, xL, xU, yL, yU
3636
return rl, ru
3737
end
3838
function transform_rule(::IntervalTransform, ::typeof(*), zL, zU, xL, xU, yL, yU)
39-
rl = Equation(zL, IfElse.ifelse(yL >= 0.0,
40-
IfElse.ifelse(xL >= 0.0, xL*yL,
41-
IfElse.ifelse(xU <= 0.0, xL*yU, xL*yU)),
42-
IfElse.ifelse(yU <= 0.0,
43-
IfElse.ifelse(xL >= 0.0, xU*yL,
44-
IfElse.ifelse(xU <= 0.0, xU*yU, xU*yL)),
45-
IfElse.ifelse(xL > 0.0, xU*yL,
46-
IfElse.ifelse(xU < 0.0, xL*yU, min(xU*yL, xU*yU))))))
47-
ru = Equation(zU, IfElse.ifelse(yL >= 0.0,
48-
IfElse.ifelse(xL >= 0.0, xU*yU,
49-
IfElse.ifelse(xU <= 0.0, xU*yL, xU*yU)),
50-
IfElse.ifelse(yU <= 0.0,
51-
IfElse.ifelse(xL >= 0.0, xL*yU,
52-
IfElse.ifelse(xU <= 0.0, xL*yL, xL*yL)),
53-
IfElse.ifelse(xL > 0.0, xU*yU,
54-
IfElse.ifelse(xU < 0.0, xL*yL, max(xL*yL, xU*yU))))))
39+
rl = Equation(zL, IfElse.ifelse(yL >= 0.0, #x*pos
40+
IfElse.ifelse(xL >= 0.0, xL*yL, #pos*pos
41+
IfElse.ifelse(xU <= 0.0, xL*yU, xL*yU)), #neg*pos, mix*pos
42+
IfElse.ifelse(yU <= 0.0, #x*neg
43+
IfElse.ifelse(xL >= 0.0, xU*yL, #pos*neg
44+
IfElse.ifelse(xU <= 0.0, xU*yU, xU*yL)), #neg*neg, mix*neg
45+
IfElse.ifelse(xL > 0.0, xU*yL, #(pos)*mix
46+
IfElse.ifelse(xU < 0.0, xL*yU, min(xL*yU, xU*yL)))))) #(neg)*mix, mix*mix
47+
ru = Equation(zU, IfElse.ifelse(yL >= 0.0, #x*pos
48+
IfElse.ifelse(xL >= 0.0, xU*yU, #pos*pos
49+
IfElse.ifelse(xU <= 0.0, xU*yL, xU*yU)), #neg*pos, mix*pos
50+
IfElse.ifelse(yU <= 0.0, #x*neg
51+
IfElse.ifelse(xL >= 0.0, xL*yU, #pos*neg
52+
IfElse.ifelse(xU <= 0.0, xL*yL, xL*yL)), #neg*neg, mix*neg
53+
IfElse.ifelse(xL > 0.0, xU*yU, #pos*mix
54+
IfElse.ifelse(xU < 0.0, xL*yL, max(xL*yL, xU*yU)))))) #neg*mix, mix*mix
5555
return rl, ru
5656
end
5757

‎src/relaxation/rules.jl

Lines changed: 52 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -50,28 +50,27 @@ end
5050
function transform_rule(::McCormickTransform, ::typeof(*), zL, zU, zcv, zcc, xL, xU, xcv, xcc, yL, yU, ycv, ycc)
5151
rcv = Equation(zcv, IfElse.ifelse(xL >= 0.0,
5252
IfElse.ifelse(yL >= 0.0, max(yU*xcv + xU*ycv - xU*yU, yL*xcv + xL*ycv - xL*yL),
53-
IfElse.ifelse(yU <= 0.0, -max((-yU)*xcv + xU*(-ycv) - xU*(-yU), (-yL)*xcv + xL*(-ycv) - xL*(-yL)),
53+
IfElse.ifelse(yU <= 0.0, -min((-yU)*xcc + xU*(-ycv) - xU*(-yU), (-yL)*xcc + xU*(-ycv) - xL*(-yL)),
5454
max(yU*xcv + xU*ycv - xU*yU, yL*xcc + xL*ycv - xL*yL))),
5555
IfElse.ifelse(xU <= 0.0,
56-
IfElse.ifelse(yL >= 0.0, -max(yU*(-xcv) + (-xU)*ycv - (-xU)*yU, yL*(-xcv) + (-xL)*ycv - (-xL)*yL),
57-
IfElse.ifelse(yU <= 0.0, max(yU*xcv + xU*ycv - xU*yU, yL*xcv + xL*ycv - xL*yL),
58-
-max(yU*(-xcv) + (-xU)*ycv - (-xU)*yU, yL*(-xcc) + (-xL)*ycv - (-xL)*yL))),
56+
IfElse.ifelse(yL >= 0.0, -min(yL*(-xcv) + (-xL)*ycc - (-xL)*yL, yU*(-xcv) + (-xU)*ycc - (-xU)*yU),
57+
IfElse.ifelse(yU <= 0.0, max(yL*xcc + xL*ycc - xL*yL, yU*xcc + xU*ycc - xU*yU),
58+
-min(yL*(-xcc) + (-xL)*ycc - (-xL)*yL, yU*(-xcv) + (-xU)*ycc - (-xU)*yU))),
5959
IfElse.ifelse(yL >= 0.0, max(xU*ycv + yU*xcv - yU*xU, xL*ycc + yL*xcv - yL*xL),
60-
IfElse.ifelse(yU <= 0.0, -max(xU*(-ycv) + (-yU)*xcv - (-yU)*xU, xL*(-ycc) + (-yL)*xcv - (-yL)*xL),
60+
IfElse.ifelse(yU <= 0.0, -min(xL*(-ycc) + (-yL)*xcc - (-yL)*xL, xU*(-ycv) + (-yU)*xcc - (-yU)*xU),
6161
max(yU*xcv + xU*ycv - xU*yU, yL*xcc + xL*ycc - xL*yL))))))
6262
rcc = Equation(zcc, IfElse.ifelse(xL >= 0.0,
63-
IfElse.ifelse(yL >= 0.0, min(yL*xcc + xU*ycc - xU*yL, yU*xcc + xU*ycc - xL*yU),
64-
IfElse.ifelse(yU <= 0.0, -min((-yL)*xcc + xU*(-ycc) - xU*(-yL), (-yU)*xcc + xU*(-ycc) - xL*(-yU)),
63+
IfElse.ifelse(yL >= 0.0, min(yL*xcc + xU*ycc - xU*yL, yU*xcc + xL*ycc - xL*yU),
64+
IfElse.ifelse(yU <= 0.0, -max((-yL)*xcv + xU*(-ycc) - xU*(-yL), (-yU)*xcv + xL*(-ycc) - xL*(-yU)),
6565
min(yL*xcv + xU*ycc - xU*yL, yU*xcc + xL*ycc - xL*yU))),
6666
IfElse.ifelse(xU <= 0.0,
67-
IfElse.ifelse(yL >= 0.0, -min(yL*(-xcc) + (-xU)*ycc - (-xU)*yL, yU*(-xcc) + (-xU)*ycc - (-xL)*yU),
68-
IfElse.ifelse(yU <= 0.0, min(yL*xcc + xU*ycc - xU*yL, yU*xcc + xU*ycc - xL*yU),
69-
-min(yL*(-xcv) + (-xU)*ycc - (-xU)*yL, yU*(-xcc) + (-xL)*ycc - (-xL)*yU))),
67+
IfElse.ifelse(yL >= 0.0, -max(yU*(-xcc) + (-xL)*ycv - (-xL)*yU, yL*(-xcc) + (-xU)*ycv - (-xU)*yL),
68+
IfElse.ifelse(yU <= 0.0, min(yU*xcv + xL*ycv - xL*yU, yL*xcv + xL*ycv - xU*yL),
69+
-max(yU*(-xcc) + (-xL)*ycv - (-xL)*yU, yL*(-xcv) + (-xU)*ycv - (-xU)*yL))),
7070
IfElse.ifelse(yL >= 0.0, min(xL*ycv + yU*xcc - yU*xL, xU*ycc + yL*xcc - yL*xU),
71-
IfElse.ifelse(yU <= 0.0, -min(xL*(-ycv) + (-yU)*xcc - (-yU)*xL, xU*(-ycc) + (-yL)*xcc - (-yL)*xU),
71+
IfElse.ifelse(yU <= 0.0, -max(xU*(-ycc) + (-yL)*xcv - (-yL)*xU, xL*(-ycv) + (-yU)*xcv - (-yU)*xL),
7272
min(yL*xcv + xU*ycc - xU*yL, yU*xcc + xL*ycv - xL*yU))))))
7373
return rcv, rcc
74-
7574

7675
# # This is the base if-else tree for the 9 multiplication cases. I think for the transform rule,
7776
# # the way this code needs to work, the entire tree needs to be available to the solver since
@@ -90,13 +89,13 @@ function transform_rule(::McCormickTransform, ::typeof(*), zL, zU, zcv, zcc, xL,
9089
# # [x+,y+]
9190
# # Normal case of multiply_STD_NS
9291
# cv = max(yU*xcv + xU*ycv - xU*yU, yL*xcv + xL*ycv - xL*yL)
93-
# cc = min(yL*xcc + xU*ycc - xU*yL, yU*xcc + xU*ycc - xL*yU)
92+
# cc = min(yL*xcc + xU*ycc - xU*yL, yU*xcc + xL*ycc - xL*yU) #typo fixed
9493

9594
# # [x+,y-]
96-
# # This returns -mult_kernel(x, -y), so put a negative in front, and then
97-
# # replace do [x+,y+] but replace all y's with -y because y's are secretly negative?
98-
# cv = -max((-yU)*xcv + xU*(-ycv) - xU*(-yU), (-yL)*xcv + xL*(-ycv) - xL*(-yL))
99-
# cc = -min((-yL)*xcc + xU*(-ycc) - xU*(-yL), (-yU)*xcc + xU*(-ycc) - xL*(-yU))
95+
# # This returns -mult_kernel(x, -y), so it's [x+, y+] but with y's all flipped,
96+
# # and then swap and negate cv/cc
97+
# cv = -min((-yU)*xcc + xU*(-ycv) - xU*(-yU), (-yL)*xcc + xU*(-ycv) - xL*(-yL))
98+
# cc = -max((-yL)*xcv + xU*(-ycc) - xU*(-yL), (-yU)*xcv + xL*(-ycc) - xL*(-yU))
10099

101100

102101
# # [x+,ym]
@@ -106,23 +105,37 @@ function transform_rule(::McCormickTransform, ::typeof(*), zL, zU, zcv, zcc, xL,
106105

107106

108107
# # [x-,y+]
109-
# # This returns -mult_kernel(-x, y). Negative in front, replace x's with -x's
110-
# cv = -max(yU*(-xcv) + (-xU)*ycv - (-xU)*yU, yL*(-xcv) + (-xL)*ycv - (-xL)*yL)
111-
# cc = -min(yL*(-xcc) + (-xU)*ycc - (-xU)*yL, yU*(-xcc) + (-xU)*ycc - (-xL)*yU)
108+
# # This returns -mult_kernel(-x, y). Replace x with flipped x, then swap/negate cv/cc
109+
# [x+, y+]
110+
# cv = max(yU*xcv + xU*ycv - xU*yU, yL*xcv + xL*ycv - xL*yL)
111+
# cc = min(yL*xcc + xU*ycc - xU*yL, yU*xcc + xL*ycc - xL*yU)
112+
# Then flip x's
113+
# cv = max(yU*(-xcc) + (-xL)*ycv - (-xL)*yU, yL*(-xcc) + (-xU)*ycv - (-xU)*yL)
114+
# cc = min(yL*(-xcv) + (-xL)*ycc - (-xL)*yL, yU*(-xcv) + (-xU)*ycc - (-xU)*yU)
115+
# Then swap/negate cv/cc
116+
# cv = -min(yL*(-xcv) + (-xL)*ycc - (-xL)*yL, yU*(-xcv) + (-xU)*ycc - (-xU)*yU)
117+
# cc = -max(yU*(-xcc) + (-xL)*ycv - (-xL)*yU, yL*(-xcc) + (-xU)*ycv - (-xU)*yL)
112118

113119

114120
# # [x-,y-]
115-
# # This returns mult_kernel(-x,-y), but since each term is of the form x_*y_,
116-
# # (-x)*(-y)=xy, so no change from the "normal case"
117-
# cv = max(yU*xcv + xU*ycv - xU*yU, yL*xcv + xL*ycv - xL*yL)
118-
# cc = min(yL*xcc + xU*ycc - xU*yL, yU*xcc + xU*ycc - xL*yU)
121+
# # This returns mult_kernel(-x,-y). Swap x with flipped x and y with flipped
122+
# # y, and since both are flipped, no need to worry about negative signs
123+
# cv = max(yL*xcc + xL*ycc - xL*yL, yU*xcc + xU*ycc - xU*yU)
124+
# cc = min(yU*xcv + xL*ycv - xL*yU, yL*xcv + xL*ycv - xU*yL)
119125

120126

121127
# # [x-,ym]
122-
# # This returns -mult_kernel(y,-x). Negative in front, then we're doing [xm, y+],
123-
# # but swap x's and y's and replace x's with (-x)'s
124-
# cv = -max(yU*(-xcv) + (-xU)*ycv - (-xU)*yU, yL*(-xcc) + (-xL)*ycv - (-xL)*yL)
125-
# cc = -min(yL*(-xcv) + (-xU)*ycc - (-xU)*yL, yU*(-xcc) + (-xL)*ycc - (-xL)*yU)
128+
# # This returns -mult_kernel(y,-x). So we do [xm, y+] but replace x's with y's,
129+
# # and y's with flipped x's. Then swap and negate cv/cc
130+
# This is [xm, y+]
131+
# cv = max(xU*ycv + yU*xcv - yU*xU, xL*ycc + yL*xcv - yL*xL)
132+
# cc = min(xL*ycv + yU*xcc - yU*xL, xU*ycc + yL*xcc - yL*xU)
133+
# Do the swaps:
134+
# cv = max(yU*(-xcc) + (-xL)*ycv - (-xL)*yU, yL*(-xcv) + (-xU)*ycv - (-xU)*yL)
135+
# cc = min(yL*(-xcc) + (-xL)*ycc - (-xL)*yL, yU*(-xcv) + (-xU)*ycc - (-xU)*yU)
136+
# Now swap definitions
137+
# cv = -min(yL*(-xcc) + (-xL)*ycc - (-xL)*yL, yU*(-xcv) + (-xU)*ycc - (-xU)*yU)
138+
# cc = -max(yU*(-xcc) + (-xL)*ycv - (-xL)*yU, yL*(-xcv) + (-xU)*ycv - (-xU)*yL)
126139

127140

128141
# # [xm,y+]
@@ -133,10 +146,17 @@ function transform_rule(::McCormickTransform, ::typeof(*), zL, zU, zcv, zcc, xL,
133146

134147

135148
# # [xm,y-]
136-
# # This one returns -mult_kernel(-y, x). So copy [x+,ym], but swap x's
137-
# # and y's, then replace y with (-y), and negative in front
138-
# cv = -max(xU*(-ycv) + (-yU)*xcv - (-yU)*xU, xL*(-ycc) + (-yL)*xcv - (-yL)*xL)
139-
# cc = -min(xL*(-ycv) + (-yU)*xcc - (-yU)*xL, xU*(-ycc) + (-yL)*xcc - (-yL)*xU)
149+
# # This one returns -mult_kernel(-y, x). So copy [x+,ym], but replace y with x,
150+
# and x with (-y), and then swap/negate cv/cc
151+
# [x+, ym]
152+
# cv = max(yU*xcv + xU*ycv - xU*yU, yL*xcc + xL*ycv - xL*yL)
153+
# cc = min(yL*xcv + xU*ycc - xU*yL, yU*xcc + xL*ycc - xL*yU)
154+
# swap y with x, x with (-y)
155+
# cv = max(xU*(-ycc) + (-yL)*xcv - (-yL)*xU, xL*(-ycv) + (-yU)*xcv - (-yU)*xL)
156+
# cc = min(xL*(-ycc) + (-yL)*xcc - (-yL)*xL, xU*(-ycv) + (-yU)*xcc - (-yU)*xU)
157+
# And now swap/negate
158+
# cv = -min(xL*(-ycc) + (-yL)*xcc - (-yL)*xL, xU*(-ycv) + (-yU)*xcc - (-yU)*xU)
159+
# cc = -max(xU*(-ycc) + (-yL)*xcv - (-yL)*xU, xL*(-ycv) + (-yU)*xcv - (-yU)*xL)
140160

141161

142162
# # [xm,ym]

0 commit comments

Comments
 (0)