Report index
StatProfilerHTML.jl report
Generated on Thu, 26 Mar 2020 19:09:01
File source code
Line Exclusive Inclusive Code
1 # We reverse the order of comparisons here so that the result
2 # of x < y is equal to the result of Monomial(x) < Monomial(y)
3 Base.@pure Base.isless(v1::AbstractVariable, v2::AbstractVariable) = name(v1) > name(v2)
4 Base.isless(m1::AbstractTermLike, m2::AbstractTermLike) = isless(promote(m1, m2)...)
5
6 function Base.isless(t1::AbstractTerm, t2::AbstractTerm)
7 15 (2 %) 256 (31 %)
256 (31 %) samples spent in isless
15 (6 %) (ex.), 256 (100 %) (incl.) when called from < line 268
241 (100 %) samples spent calling <
if monomial(t1) < monomial(t2)
8 true
9 56 (7 %)
56 (7 %) samples spent in isless
56 (100 %) (incl.) when called from < line 268
56 (100 %) samples spent calling ==
elseif monomial(t1) == monomial(t2)
10 36 (4 %)
36 (4 %) samples spent in isless
36 (100 %) (incl.) when called from < line 268
26 (72 %) samples spent calling <
10 (28 %) samples spent calling abs
abs(coefficient(t1)) < abs(coefficient(t2))
11 else
12 false
13 end
14 end
15
16 # promoting multiplication is not a good idea
17 # For example a polynomial of Float64 * a polynomial of JuMP affine expression
18 # is a polynomial of JuMP affine expression but if we promote it would be a
19 # polynomial of quadratic expression
20 for op in [:+, :-, :(==)]
21 @eval Base.$op(p1::APL, p2::APL) = $op(promote(p1, p2)...)
22 end
23 # Promotion between `I` and `1` is `Any`.
24 # Promotion between `I` and `2I` is `UniformScaling`.
25 for op in [:+, :-]
26 @eval Base.$op(p1::APL, p2::APL{<:LinearAlgebra.UniformScaling}) = $op(p1, mapcoefficientsnz(J -> J.λ, p2))
27 @eval Base.$op(p1::APL{<:LinearAlgebra.UniformScaling}, p2::APL) = $op(mapcoefficientsnz(J -> J.λ, p1), p2)
28 @eval Base.$op(p1::APL{<:LinearAlgebra.UniformScaling}, p2::APL{<:LinearAlgebra.UniformScaling}) = $op(mapcoefficientsnz(J -> J.λ, p1), p2)
29 end
30 Base.isapprox(p1::APL, p2::APL; kwargs...) = isapprox(promote(p1, p2)...; kwargs...)
31
32 # @eval $op(p::APL, α) = $op(promote(p, α)...) would be less efficient
33 for (op, fun) in [(:+, :plusconstant), (:-, :minusconstant), (:*, :multconstant), (:(==), :eqconstant)]
34 @eval Base.$op(p::APL, α) = $fun(p, α)
35 @eval Base.$op(α, p::APL) = $fun(α, p)
36 end
37 # Fix ambiguity between above methods and methods in MA
38 Base.:+(::MA.Zero, p::APL) = MA.copy_if_mutable(p)
39 Base.:+(p::APL, ::MA.Zero) = MA.copy_if_mutable(p)
40
41 # Special case AbstractArrays of APLs
42 # We add these instead of relying on the broadcasting API since the above method definitions are very wide.
43 # In particular, there is support for Matrices as coefficents. In order to avoid issues like #104 we therefore
44 # explicitly define this (instead of implictly getting unexpected results).
45 for op in [:+, :-]
46 @eval Base.$op(p::APL, A::AbstractArray{<:APL}) = map(f -> $op(p, f), A)
47 @eval Base.$op(A::AbstractArray{<:APL}, p::APL) = map(f -> $op(f, p), A)
48 end
49 Base.:*(p::APL, A::AbstractArray) = map(f -> p * f, A)
50 Base.:*(A::AbstractArray, p::APL) = map(f -> f * p, A)
51 Base.:/(A::AbstractArray, p::APL) = map(f -> f / p, A)
52
53 constant_function(::typeof(+)) = plusconstant
54 constant_function(::typeof(-)) = minusconstant
55 MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::APL, α) = MA.mutable_operate!(constant_function(op), p, α)
56 MA.mutable_operate_to!(output::AbstractPolynomial, op::typeof(*), α, p::APL) = MA.mutable_operate_to!(output, multconstant, α, p)
57 MA.mutable_operate_to!(output::AbstractPolynomial, op::typeof(*), p::APL, α) = MA.mutable_operate_to!(output, multconstant, p, α)
58
59 function polynomial_merge!(
60 n1::Int, n2::Int, get1::Function, get2::Function,
61 set::Function, push::Function, compare_monomials::Function,
62 combine::Function, keep::Function, resize::Function)
63 buffer = nothing
64 i = j = k = 1
65 # Invariant:
66 # The terms p[0] -> p[k-1] are sorted and are smaller than the remaining terms.
67 # The terms p[k] -> p[i-1] are garbage.
68 # The terms p[i] -> p[end] are sorted and still need to be added.
69 # The terms q[j] -> p[end] are sorted and still need to be added.
70 # If `buffer` is not empty:
71 # The terms in `buffer` are sorted and still need to be added.
72 # Moreover, they are smaller than the terms p[i] -> p[end].
73 while i <= n1 && j <= n2
74 @assert buffer === nothing || isempty(buffer)
75 comp = compare_monomials(i, j)
76 if comp > 0
77 if k == i
78 t0 = get1(i)
79 if buffer === nothing
80 buffer = DataStructures.Queue{typeof(t0)}()
81 end
82 DataStructures.enqueue!(buffer, t0)
83 i += 1
84 end
85 set(k, get2(j))
86 j += 1
87 k += 1
88 elseif iszero(comp)
89 combine(i, j)
90 if keep(i)
91 if k != i
92 @assert k < i
93 set(k, get1(i))
94 end
95 k += 1
96 end
97 i += 1
98 j += 1
99 else
100 if k != i
101 set(k, get1(i))
102 end
103 i += 1
104 k += 1
105 end
106 while buffer !== nothing && !isempty(buffer) && j <= n2
107 @assert i == k
108 t = first(buffer)
109 comp = compare_monomials(t, j)
110 if comp >= 0
111 if comp > 0
112 t = get2(j)
113 else
114 t = combine(t, j)
115 end
116 j += 1
117 end
118 if comp <= 0
119 DataStructures.dequeue!(buffer)
120 end
121 # if `comp` is zero, we called `combine` so `t`
122 # might not be kept. If `comp` is not zero, we
123 # skip the `keep` call that might be costly.
124 if iszero(comp) && !keep(t)
125 continue
126 end
127 if k <= n1
128 DataStructures.enqueue!(buffer, get1(i))
129 set(k, t)
130 else
131 push(t)
132 n1 += 1
133 end
134 i += 1
135 k += 1
136 end
137 end
138 if buffer !== nothing && !isempty(buffer)
139 @assert j == n2 + 1
140 @assert i == k
141 n = length(buffer)
142 resize(n1 + n)
143 for k in n1:-1:i
144 set(k + n, get1(k))
145 end
146 for k in i:(i + n - 1)
147 set(k, DataStructures.dequeue!(buffer))
148 end
149 n1 += n
150 else
151 len = (k - 1) + (n2 - (j - 1)) + (n1 - (i - 1))
152 if n1 < len
153 resize(len)
154 end
155 while j <= n2
156 set(k, get2(j))
157 j += 1
158 k += 1
159 end
160 @assert j == n2 + 1
161 while i <= n1
162 set(k, get1(i))
163 i += 1
164 k += 1
165 end
166 if len < n1
167 resize(len)
168 end
169 @assert i == n1 + 1
170 @assert k == len + 1
171 end
172 return
173 end
174
175
176 #MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomial) = MA.mutable_operate_to!(p, op, p, q)
177 MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomialLike) = MA.mutable_operate!(op, p, polynomial(q))
178
179
135 (16 %) samples spent in mul_to_terms!
135 (100 %) (incl.) when called from * line 188
function mul_to_terms!(ts::Vector{<:AbstractTerm}, p1::APL, p2::APL)
180 for t1 in terms(p1)
181 for t2 in terms(p2)
182 135 (16 %)
2 (1 %) samples spent calling iterate
1 (1 %) samples spent calling *
86 (64 %) samples spent calling push!
46 (34 %) samples spent calling push!
push!(ts, t1 * t2)
183 end
184 end
185 return ts
186 end
187
833 (100 %) samples spent in *
3 (0 %) (ex.), 783 (94 %) (incl.) when called from power_by_squaring line 228
50 (6 %) (incl.) when called from power_by_squaring line 226
function Base.:*(p::AbstractPolynomial, q::AbstractPolynomial)
188 3 (0 %) 833 (100 %)
694 (84 %) samples spent calling polynomial!
1 (0 %) samples spent calling mul_to_terms!
135 (16 %) samples spent calling mul_to_terms!
polynomial!(mul_to_terms!(MA.promote_operation(*, termtype(p), termtype(q))[], p, q))
189 end
190
191 Base.isapprox(p::APL, α; kwargs...) = isapprox(promote(p, α)...; kwargs...)
192 Base.isapprox(α, p::APL; kwargs...) = isapprox(promote(p, α)...; kwargs...)
193
194 # `MA.operate(-, p)` redirects to `-p` as it assumes that `-p` can be modified
195 # through the MA API without modifying `p`. We should either copy the monomial
196 # here or implement a `MA.operate(-, p)` that copies it. We choose the first
197 # option.
198 Base.:-(m::AbstractMonomialLike) = (-1) * MA.copy_if_mutable(m)
199 Base.:-(t::AbstractTermLike) = MA.operate(-, coefficient(t)) * MA.copy_if_mutable(monomial(t))
200 Base.:-(p::APL) = polynomial!((-).(terms(p)))
201 Base.:+(p::Union{APL, RationalPoly}) = p
202 Base.:*(p::Union{APL, RationalPoly}) = p
203
204 # Avoid adding a zero constant that might artificially increase the Newton polytope
205 # Need to add polynomial conversion for type stability
206 plusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(+, S, T)) : p + constantterm(α, p)
207 plusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(+, S, T)) : constantterm(α, p) + p
208 function MA.mutable_operate!(::typeof(plusconstant), p::APL, α)
209 if !iszero(α)
210 MA.mutable_operate!(+, p, constantterm(α, p))
211 end
212 return p
213 end
214 minusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(-, S, T)) : p - constantterm(α, p)
215 minusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial(-p, MA.promote_operation(-, S, T)) : constantterm(α, p) - p
216 function MA.mutable_operate!(::typeof(minusconstant), p::APL, α)
217 if !iszero(α)
218 MA.mutable_operate!(-, p, constantterm(α, p))
219 end
220 return p
221 end
222
223 # Coefficients and variables commute
224 multconstant(α, v::AbstractVariable) = multconstant(α, monomial(v)) # TODO linear term
225 multconstant(m::AbstractMonomialLike, α) = multconstant(α, m)
226
227 _multconstant(α, f, t::AbstractTermLike) = mapcoefficientsnz(f, t)
228 function _multconstant(α::T, f, p::AbstractPolynomial{S}) where {S, T}
229 if iszero(α)
230 zero(polynomialtype(p, MA.promote_operation(*, T, S)))
231 else
232 mapcoefficientsnz(f, p)
233 end
234 end
235 _multconstant(α, f, p::AbstractPolynomialLike) = _multconstant(α, f, polynomial(p))
236
237 multconstant(α, p::AbstractPolynomialLike) = _multconstant(α, β -> α*β, p)
238 multconstant(p::AbstractPolynomialLike, α) = _multconstant(α, β -> β*α, p)
239
240 function mapcoefficientsnz_to! end
241
242 function _multconstant_to!(output, α, f, p)
243 if iszero(α)
244 MA.mutable_operate!(zero, output)
245 else
246 mapcoefficientsnz_to!(output, f, p)
247 end
248 end
249 function MA.mutable_operate_to!(output, ::typeof(multconstant), p::APL, α)
250 _multconstant_to!(output, α, β -> β*α, p)
251 end
252 function MA.mutable_operate_to!(output, ::typeof(multconstant), α, p::APL)
253 _multconstant_to!(output, α, β -> α*β, p)
254 end
255
256 MA.mutable_operate_to!(output::AbstractMonomial, ::typeof(*), m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents_to!(output, +, m1, m2)
257 MA.mutable_operate!(::typeof(*), m1::AbstractMonomial, m2::AbstractMonomialLike) = mapexponents!(+, m1, m2)
258 1 (0 %)
1 (0 %) samples spent in *
1 (100 %) (incl.) when called from * line 263
1 (100 %) samples spent calling mapexponents
Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(+, m1, m2)
259 #Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = *(monomial(m1), monomial(m2))
260
261 Base.:*(m::AbstractMonomialLike, t::AbstractTermLike) = coefficient(t) * (m * monomial(t))
262 Base.:*(t::AbstractTermLike, m::AbstractMonomialLike) = coefficient(t) * (monomial(t) * m)
263 1 (0 %)
1 (0 %) samples spent in *
1 (100 %) (incl.) when called from mul_to_terms! line 182
1 (100 %) samples spent calling *
Base.:*(t1::AbstractTermLike, t2::AbstractTermLike) = (coefficient(t1) * coefficient(t2)) * (monomial(t1) * monomial(t2))
264
265 Base.:*(t::AbstractTermLike, p::APL) = polynomial!(map(te -> t * te, terms(p)))
266 Base.:*(p::APL, t::AbstractTermLike) = polynomial!(map(te -> te * t, terms(p)))
267 Base.:*(p::APL, q::APL) = polynomial(p) * polynomial(q)
268
269 # guaranteed that monomial(t1) > monomial(t2)
270 function _polynomial_2terms(t1::TT, t2::TT, ::Type{T}) where {TT<:AbstractTerm, T}
271 if iszero(t1)
272 polynomial(t2, T)
273 elseif iszero(t2)
274 polynomial(t1, T)
275 else
276 # not `polynomial!` because we `t1` and `t2` cannot be modified
277 polynomial(termtype(TT, T)[t1, t2], SortedUniqState())
278 end
279 end
280 for op in [:+, :-]
281 @eval begin
282 Base.$op(t1::AbstractTermLike, t2::AbstractTermLike) = $op(term(t1), term(t2))
283 Base.$op(t1::AbstractTerm, t2::AbstractTerm) = $op(_promote_terms(t1, t2)...)
284 function Base.$op(t1::TT, t2::TT) where {T, TT <: AbstractTerm{T}}
285 S = MA.promote_operation($op, T, T)
286 # t1 > t2 would compare the coefficient in case the monomials are equal
287 # and it will throw a MethodError in case the coefficients are not comparable
288 if monomial(t1) == monomial(t2)
289 polynomial($op(coefficient(t1), coefficient(t2)) * monomial(t1), S)
290 elseif monomial(t1) > monomial(t2)
291 ts = _polynomial_2terms(t1, $op(t2), S)
292 else
293 ts = _polynomial_2terms($op(t2), t1, S)
294 end
295 end
296 end
297 end
298 _promote_terms(t1, t2) = promote(t1, t2)
299 # Promotion between `I` and `1` is `Any`.
300 _promote_terms(t1::AbstractTerm, t2::AbstractTerm{<:LinearAlgebra.UniformScaling}) = _promote_terms(t1, coefficient(t2).λ * monomial(t2))
301 _promote_terms(t1::AbstractTerm{<:LinearAlgebra.UniformScaling}, t2::AbstractTerm) = _promote_terms(coefficient(t1).λ * monomial(t1), t2)
302 # Promotion between `I` and `2I` is `UniformScaling`, not `UniformScaling{Int}`.
303 function _promote_terms(t1::AbstractTerm{LinearAlgebra.UniformScaling{S}}, t2::AbstractTerm{LinearAlgebra.UniformScaling{T}}) where {S<:Number, T<:Number}
304 U = LinearAlgebra.UniformScaling{promote_type(S, T)}
305 _promote_terms(MA.scaling_convert(U, coefficient(t1)) * monomial(t1), MA.scaling_convert(U, coefficient(t2)) * monomial(t2))
306 end
307 function _promote_terms(t1::AbstractTerm{LinearAlgebra.UniformScaling{T}}, t2::AbstractTerm{LinearAlgebra.UniformScaling{T}}) where T<:Number
308 return promote(t1, t2)
309 end
310
311 LinearAlgebra.adjoint(v::AbstractVariable) = v
312 LinearAlgebra.adjoint(m::AbstractMonomial) = m
313 LinearAlgebra.adjoint(t::AbstractTerm) = LinearAlgebra.adjoint(coefficient(t)) * monomial(t)
314 LinearAlgebra.adjoint(p::AbstractPolynomialLike) = polynomial(map(LinearAlgebra.adjoint, terms(p)))
315 LinearAlgebra.adjoint(r::RationalPoly) = adjoint(numerator(r)) / adjoint(denominator(r))
316
317 LinearAlgebra.transpose(v::AbstractVariable) = v
318 LinearAlgebra.transpose(m::AbstractMonomial) = m
319 LinearAlgebra.transpose(t::AbstractTerm) = LinearAlgebra.transpose(coefficient(t)) * monomial(t)
320 LinearAlgebra.transpose(p::AbstractPolynomialLike) = polynomial(map(LinearAlgebra.transpose, terms(p)))
321 LinearAlgebra.transpose(r::RationalPoly) = transpose(numerator(r)) / transpose(denominator(r))
322
323 LinearAlgebra.dot(p1::AbstractPolynomialLike, p2::AbstractPolynomialLike) = p1' * p2
324 LinearAlgebra.dot(x, p::AbstractPolynomialLike) = x' * p
325 LinearAlgebra.dot(p::AbstractPolynomialLike, x) = p' * x
326
327 LinearAlgebra.symmetric_type(PT::Type{<:APL}) = PT
328 LinearAlgebra.symmetric(p::APL, ::Symbol) = p
329
330 # Amazingly, this works! Thanks, StaticArrays.jl!
331 """
332 Convert a tuple of variables into a static vector to allow array-like usage.
333 The element type of the vector will be Monomial{vars, length(vars)}.
334 """
335 Base.vec(vars::Tuple{Vararg{AbstractVariable}}) = [vars...]
336 # vec(vars::Tuple{Vararg{<:TypedVariable}}) = SVector(vars)
337
338 # https://github.com/JuliaLang/julia/pull/23332
339 833 (100 %)
833 (100 %) samples spent in ^
833 (100 %) (incl.) when called from macro expansion line 0
833 (100 %) samples spent calling power_by_squaring
Base.:^(x::AbstractPolynomialLike, p::Integer) = Base.power_by_squaring(x, p)
340
341 function MA.mutable_operate_to!(output::AbstractPolynomial, op::MA.AddSubMul, x, args::Vararg{Any, N}) where N
342 return MA.mutable_operate_to!(output, MA.add_sub_op(op), x, *(args...))
343 end
344 function MA.mutable_operate!(op::MA.AddSubMul, x, y, z, args::Vararg{Any, N}) where N
345 return MA.mutable_operate!(MA.add_sub_op(op), x, *(y, z, args...))
346 end
347 MA.buffer_for(::MA.AddSubMul, ::Type{<:AbstractPolynomial}, args::Vararg{Type, N}) where {N} = zero(MA.promote_operation(*, args...))
348 function MA.mutable_buffered_operate_to!(buffer::AbstractPolynomial, output::AbstractPolynomial, op::MA.AddSubMul, x, y, z, args::Vararg{Any, N}) where N
349 product = MA.operate_to!(buffer, *, y, z, args...)
350 return MA.mutable_operate_to!(output, MA.add_sub_op(op), x, product)
351 end
352 function MA.mutable_buffered_operate!(buffer::AbstractPolynomial, op::MA.AddSubMul, x, y, z, args::Vararg{Any, N}) where N
353 product = MA.operate_to!(buffer, *, y, z, args...)
354 return MA.mutable_operate!(MA.add_sub_op(op), x, product)
355 end