Report index
StatProfilerHTML.jl report
Generated on Thu, 26 Mar 2020 19:09:01
File source code
Line Exclusive Inclusive Code
1 Base.one(::Type{P}) where {C, T, P <: Polynomial{C, T}} = Polynomial(one(T))
2 Base.one(p::Polynomial) = one(typeof(p))
3
4 Base.zero(::Type{Polynomial{C, T, A}}) where {C, T, A} = Polynomial(A())
5 Base.zero(t::PolynomialLike) = zero(typeof(t))
6
7 combine(t1::Term, t2::Term) = combine(promote(t1, t2)...)
8 combine(t1::T, t2::T) where {T <: Term} = Term(t1.coefficient + t2.coefficient, t1.monomial)
9 compare(t1::Term, t2::Term) = monomial(t1) > monomial(t2)
10
11 # Graded Lexicographic order
12 # First compare total degree, then lexicographic order
13 function isless(m1::Monomial{V}, m2::Monomial{V}) where {V}
14 14 (2 %)
14 (2 %) samples spent in isless
14 (100 %) (incl.) when called from < line 268
14 (100 %) samples spent calling degree
d1 = degree(m1)
15 123 (15 %)
123 (15 %) samples spent in isless
123 (100 %) (incl.) when called from < line 268
123 (100 %) samples spent calling degree
d2 = degree(m2)
16 53 (6 %)
53 (6 %) samples spent in isless
53 (100 %) (incl.) when called from < line 268
53 (100 %) samples spent calling <
if d1 < d2
17 return true
18 9 (1 %)
9 (1 %) samples spent in isless
9 (100 %) (incl.) when called from < line 268
9 (100 %) samples spent calling >
elseif d1 > d2
19 return false
20 else
21 10 (1 %)
10 (1 %) samples spent in isless
10 (100 %) (incl.) when called from < line 268
1 (10 %) samples spent calling <
9 (90 %) samples spent calling <
return exponents(m1) < exponents(m2)
22 end
23 end
24 function grlex(m1::Monomial{V}, m2::Monomial{V}) where {V}
25 d1 = degree(m1)
26 d2 = degree(m2)
27 if d1 != d2
28 return d1 - d2
29 else
30 if exponents(m1) == exponents(m2)
31 return 0
32 elseif exponents(m1) < exponents(m2)
33 return -1
34 else
35 return 1
36 end
37 end
38 end
39 grlex(m1::Monomial, m2::Monomial) = grlex(promote(m1, m2)...)
40
41 jointerms(terms1::AbstractArray{<:Term}, terms2::AbstractArray{<:Term}) = mergesorted(terms1, terms2, compare, combine)
42 function jointerms!(output::AbstractArray{<:Term}, terms1::AbstractArray{<:Term}, terms2::AbstractArray{<:Term})
43 resize!(output, length(terms1) + length(terms2))
44 Sequences.mergesorted!(output, terms1, terms2, compare, combine)
45 end
46
47 (+)(p1::Polynomial, p2::Polynomial) = Polynomial(jointerms(terms(p1), terms(p2)))
48 (+)(p1::Polynomial, p2::Polynomial{<:UniformScaling}) = p1 + MP.mapcoefficientsnz(J -> J.λ, p2)
49 (+)(p1::Polynomial{<:UniformScaling}, p2::Polynomial) = MP.mapcoefficientsnz(J -> J.λ, p1) + p2
50 (+)(p1::Polynomial{<:UniformScaling}, p2::Polynomial{<:UniformScaling}) = MP.mapcoefficientsnz(J -> J.λ, p1) + p2
51 function MA.mutable_operate_to!(result::Polynomial, ::typeof(+), p1::Polynomial, p2::Polynomial)
52 if result === p1 || result === p2
53 error("Cannot call `mutable_operate_to!(output, +, p, q)` with `output` equal to `p` or `q`, call `mutable_operate!` instead.")
54 end
55 jointerms!(result.terms, terms(p1), terms(p2))
56 result
57 end
58 function MA.mutable_operate_to!(result::Polynomial, ::typeof(*), p::Polynomial, t::AbstractTermLike)
59 if iszero(t)
60 MA.mutable_operate!(zero, result)
61 else
62 resize!(result.terms, nterms(p))
63 for i in eachindex(p.terms)
64 # TODO could use MA.mul_to! for indices that were presents in `result` before the `resize!`.
65 result.terms[i] = p.terms[i] * t
66 end
67 return result
68 end
69 end
70 function MA.mutable_operate_to!(result::Polynomial, ::typeof(*), t::AbstractTermLike, p::Polynomial)
71 if iszero(t)
72 MA.mutable_operate!(zero, result)
73 else
74 resize!(result.terms, nterms(p))
75 for i in eachindex(p.terms)
76 # TODO could use MA.mul_to! for indices that were presents in `result` before the `resize!`.
77 result.terms[i] = t * p.terms[i]
78 end
79 return result
80 end
81 end
82 (-)(p1::Polynomial, p2::Polynomial) = Polynomial(jointerms(terms(p1), (-).(terms(p2))))
83 (-)(p1::Polynomial, p2::Polynomial{<:UniformScaling}) = p1 - MP.mapcoefficientsnz(J -> J.λ, p2)
84 (-)(p1::Polynomial{<:UniformScaling}, p2::Polynomial) = MP.mapcoefficientsnz(J -> J.λ, p1) - p2
85 (-)(p1::Polynomial{<:UniformScaling}, p2::Polynomial{<:UniformScaling}) = MP.mapcoefficientsnz(J -> J.λ, p1) - p2
86
87 function MA.mutable_operate!(op::Union{typeof(+), typeof(-)}, p::Polynomial{T}, q::Polynomial) where T
88 get1(i) = p.terms[i]
89 function get2(i)
90 t = q.terms[i]
91 Term(MA.scaling_convert(T, MA.operate(op, t.coefficient)), t.monomial)
92 end
93 set(i, t::Term) = p.terms[i] = t
94 push(t::Term) = push!(p.terms, t)
95 compare_monomials(t::Term, j::Int) = grlex(q.terms[j].monomial, t.monomial)
96 compare_monomials(i::Int, j::Int) = compare_monomials(get1(i), j)
97 combine(i::Int, j::Int) = p.terms[i] = Term(MA.operate!(op, p.terms[i].coefficient, q.terms[j].coefficient), p.terms[i].monomial)
98 combine(t::Term, j::Int) = Term(MA.operate!(op, t.coefficient, q.terms[j].coefficient), t.monomial)
99 resize(n) = resize!(p.terms, n)
100 # We can modify the coefficient since it's the result of `combine`.
101 keep(t::Term) = !MA.iszero!(t.coefficient)
102 keep(i::Int) = !MA.iszero!(p.terms[i].coefficient)
103 MP.polynomial_merge!(
104 nterms(p), nterms(q), get1, get2, set, push,
105 compare_monomials, combine, keep, resize
106 )
107 return p
108 end
109
110 (==)(::Variable{N}, ::Variable{N}) where {N} = true
111 (==)(::Variable, ::Variable) = false
112 57 (7 %)
57 (7 %) samples spent in ==
1 (2 %) (incl.) when called from uniqterms! line 108
56 (98 %) (incl.) when called from isless line 9
57 (100 %) samples spent calling ==
(==)(m1::Monomial{V}, m2::Monomial{V}) where {V} = exponents(m1) == exponents(m2)
113
114 # Multiplication is handled as a special case so that we can write these
115 # definitions without resorting to promotion:
116 MP.multconstant(α, v::Monomial) = Term(α, v)
117 MP.multconstant(α, v::Variable) = Term(α, Monomial(v))
118
119 (*)(v1::V, v2::V) where {V <: Variable} = Monomial{(V(),), 1}((2,))
120 (*)(v1::Variable, v2::Variable) = (*)(promote(v1, v2)...)
121
122 function MP.divides(m1::Monomial{V, N}, m2::Monomial{V, N}) where {V, N}
123 reduce((d, exp) -> d && (exp[1] <= exp[2]), zip(m1.exponents, m2.exponents), init=true)
124 end
125 MP.divides(m1::Monomial, m2::Monomial) = divides(promote(m1, m2)...)
126 function MP.mapexponents(op, m1::M, m2::M) where M<:Monomial
127 1 (0 %)
1 (0 %) samples spent in mapexponents
1 (100 %) (incl.) when called from * line 258
1 (100 %) samples spent calling map
M(map(op, m1.exponents, m2.exponents))
128 end
129 MP.mapexponents(op, m1::Monomial, m2::Monomial) = mapexponents(op, promote(m1, m2)...)
130 #function MP.mapexponents_to!(output::M, op, m1::M, m2::M) where M<:Monomial
131 # map!(op, output.exponents, m1.exponents, m2.exponents)
132 # return output
133 #end
134 #function MP.mapexponents!(op, m1::M, m2::M) where M<:Monomial
135 # map!(op, m1.exponents, m1.exponents, m2.exponents)
136 # return m1
137 #end
138
139 function MA.mutable_operate_to!(output::Polynomial, ::typeof(*), p::Polynomial, q::Polynomial)
140 empty!(output.terms)
141 MP.mul_to_terms!(output.terms, p, q)
142 sort!(output.terms, lt=(>))
143 MP.uniqterms!(output.terms)
144 return output
145 end
146 function MA.mutable_operate!(::typeof(*), p::Polynomial, q::Polynomial)
147 return MA.mutable_operate_to!(p, *, MA.mutable_copy(p), q)
148 end
149
150 function MA.mutable_operate!(::typeof(zero), p::Polynomial)
151 empty!(p.terms)
152 return p
153 end
154 function MA.mutable_operate!(::typeof(one), p::Polynomial{T}) where T
155 if isempty(p.terms)
156 push!(p.terms, constantterm(one(T), p))
157 else
158 t = p.terms[1]
159 p.terms[1] = Term(MA.one!(coefficient(t)), constantmonomial(t))
160 resize!(p.terms, 1)
161 end
162 return p
163 end
164
165 # The exponents are stored in a tuple, this is not mutable.
166 # We could remove these methods since it is the default.
167 MA.mutability(::Type{<:Monomial}) = MA.NotMutable()
168 MA.mutability(::Type{<:Term}) = MA.NotMutable()
169 # The polynomials can be mutated.
170 MA.mutability(::Type{<:Polynomial}) = MA.IsMutable()
171
172 ^(v::V, x::Integer) where {V <: Variable} = Monomial{(V(),), 1}((x,))
173
174 # dot(v1::AbstractVector{<:TermLike}, v2::AbstractVector) = dot(v1, v2)
175 # dot(v1::AbstractVector, v2::AbstractVector{<:TermLike}) = dot(v1, v2)
176 # dot(v1::AbstractVector{<:TermLike}, v2::AbstractVector{<:TermLike}) = dot(v1, v2)
177
178 # All of these types are immutable, so there's no need to copy anything to get
179 # a shallow copy.
180 Base.copy(x::TermLike) = x
181
182 # Terms are not mutable under the MutableArithmetics API
183 MA.mutable_copy(p::Polynomial) = Polynomial([Term(MA.copy_if_mutable(t.coefficient), t.monomial) for t in terms(p)])
184 Base.copy(p::Polynomial) = MA.mutable_copy(p)
185
186 adjoint(v::Variable) = v
187 adjoint(m::Monomial) = m
188 adjoint(t::Term) = Term(adjoint(coefficient(t)), monomial(t))
189 adjoint(x::Polynomial) = Polynomial(adjoint.(terms(x)))
190
191 function MP.mapcoefficientsnz_to!(output::Polynomial, f::Function, p::Polynomial)
192 resize!(output.terms, nterms(p))
193 for i in eachindex(p.terms)
194 t = p.terms[i]
195 output.terms[i] = Term(f(coefficient(t)), monomial(t))
196 end
197 return output
198 end
199 function MP.mapcoefficientsnz_to!(output::Polynomial, f::Function, p::AbstractPolynomialLike)
200 return MP.mapcoefficientsnz_to!(output, f, polynomial(p))
201 end