Line data Source code
1 15 : _as_int(i) = Int(i)
2 0 : _as_int(::Val{N}) where {N} = _as_int(N)
3 :
4 : """
5 : i = multiindices(d, q)
6 :
7 : Generate sequence of `d`-dimensional multi-indices which sum to
8 : `q`. `d` and `q` can be integers `n`n or `Val(n)`.
9 :
10 : Example: `multiindices(3, 2)` yields the sequence
11 :
12 : (2,0,0), (1,1,0), (0,2,0),
13 : (1,0,1), (0,1,1),
14 : (0,0,2)
15 :
16 : See also [`domainpoints`](@ref), [`mvbernstein`](@ref)
17 : """
18 8 : function multiindices(d, q)
19 8 : d = _as_int(d)
20 7 : q = _as_int(q)
21 :
22 159 : map(filter(i -> sum(i.I) == q, CartesianIndices(ntuple(_ -> 0:q, d)))) do i
23 30 : tuple(i.I...)
24 : end
25 : end
26 :
27 12 : monomial(q::Int, i::Int, x) = x ^ i
28 4 : monomial(::Val{Q}, i::Int, x) where {Q} = monomial(Q, i, x)
29 :
30 10 : monomial(q::Int, x) = ntuple(i -> monomial(q, i-1, x), q+1)
31 5 : monomial(q::Val{Q}, x) where {Q} = ntuple(i -> monomial(q, i-1, x), Val(Q+1))
32 :
33 : """
34 : y = monomial(q, i, x)
35 : y = monomial(Val(q), i, x)
36 :
37 : Evaluate `i`-th (univariate) monomials `x^i` of degree `q` at `x` for
38 : `i in 0:q`.
39 :
40 : ys = monomial(q, x)
41 : ys = monomial(Val(q), x)
42 :
43 : Evaluate all monomials of degree `q` at `x` to tuple.
44 :
45 : See also [`mvmonomial`](@ref), [`bernstein`](@ref)
46 : """ monomial
47 :
48 5 : function _exmvmq(d, q)
49 25 : function terms(idx)
50 20 : ts = map(1:length(idx)) do i
51 60 : quote
52 869 : x[$i] ^ $(idx[i])
53 : end
54 : end
55 :
56 20 : Expr(:call, :prod, Expr(:vect, ts...))
57 : end
58 :
59 5 : [ terms(idx) for idx in multiindices(d, q) ]
60 : end
61 :
62 2 : function _exmvm(d, q)
63 2 : foldl(vcat, [ _exmvmq(d, q) for q in 0:q ])
64 : end
65 :
66 11 : @generated function mvmonomials(::Val{Q}, x::NTuple{D, T}) where {Q, D, T}
67 2 : quote
68 11 : tuple($(_exmvm(D, Q)...))
69 : end
70 : end
71 :
72 1 : mvmonomials(q::Val{Q}, x::SVector) where {Q} = mvmonomials(q, tuple(x...))
73 6 : mvmonomials(q::Val{Q}, x...) where {Q} = mvmonomials(q, tuple(x...))
74 :
75 11 : mvmonomials(q::Int, args...) = mvmonomials(Val(q), args...)
76 :
77 : """
78 : values = mvmonomial(Val(Q), v, x...)
79 : values = mvmonomials(q, v, x...)
80 :
81 : Evaluate multivariate monomials of total degree `Q` or `q`. Evaluates
82 : - univariate polynomials `1, x, x^2, ...` if `x` is a number
83 : - bivariate polynomials `1, x, y, x^2, x*y, y^2, x^3, x^2*y, ...` if
84 : `x...` is `x, y` or a 2-tuple or `SVector{2}`
85 : - trivariate polynomials `1, x, y, z, x^2, x*y, y^2, y*z, z^2, ...`
86 : if `x...` is `x, y, z` or a 3-tuple or `SVector{3}`
87 : - etc.
88 :
89 : The polynomial basis is evaluated inthis canonical order and the
90 : result is returned as an n-tuple of `values`.
91 :
92 : !!! note
93 : Prefer "static dispatch" by type `Val(Q)`. This calls the
94 : appropriate method *without* branching over the polynomial degree.
95 :
96 : See also [`mvmonomials!`](@ref), [`mvpdim`](@ref)
97 : """ mvmonomials
98 :
99 :
100 3 : mvmonomials!(q, v::AbstractVector, x::NTuple{D, T}) where {D, T} = (v.= mvmonomials(q, x))
101 :
102 1 : mvmonomials!(q, v::AbstractVector, x::SVector) = mvmonomials!(q, v, tuple(x...))
103 1 : mvmonomials!(q, v::AbstractVector, x...) = mvmonomials!(q, v, tuple(x...))
104 :
105 : """
106 : mvmonomials!(Val(Q), v, x...)
107 : mvmonomials!(q, v, x...)
108 :
109 : Evaluate the multivariate monomials of total degree `Q` or `q` as for
110 : [`mvmonomials`](@ref) and store the computed values in `v`.
111 :
112 : !!! warning "precondition"
113 : `v` must have the appropriate size [`mvpdim`](@ref).
114 :
115 : See also [`mvmonomials`](@ref), [`bvmdim`](@ref)
116 : """ mvmonomials!
117 :
118 16 : _as_val(i::Int) = Val(i)
119 0 : _as_val(i::Val{N}) where {N} = i
120 0 : _as_val(::SVector{N, T}) where {N, T} = N
121 0 : _as_val(::NTuple{N, T}) where {N, T} = N
122 0 : _as_val(::Type{Val{N}}) where {N} = Val(N)
123 0 : _as_val(::Type{SVector{N, T}}) where {N, T} = N
124 0 : _as_val(::Type{NTuple{N, T}}) where {N, T} = N
125 :
126 3 : mvpdim(d, q) = mvpdim(_as_val(d), _as_val(q))
127 3 : mvpdim(::Val{D}, ::Val{Q}) where {D, Q} = binomial(Q + D, D)
128 :
129 : """
130 : n = mvpdim(Val(Q), Val(D))
131 : n = mvpdim(q, d)
132 :
133 : Get number of d-variate polynomials of total degree `q`.
134 :
135 : See also [`mvmonomials`](@ref), [`mvmonomials!`](@ref)
136 : """ mvpdim
137 :
138 12 : bernstein(q::Int, i::Int, x) = binomial(q, i) * x^i * (one(x)-x)^(q-i)
139 4 : bernstein(::Val{Q}, i::Int, x) where {Q} = bernstein(Q, i, x)
140 :
141 10 : bernstein(q::Int, x) = ntuple(i -> bernstein(q, i-1, x), q+1)
142 5 : bernstein(q::Val{Q}, x) where {Q} = ntuple(i -> bernstein(q, i-1, x), Val(Q+1))
143 :
144 : """
145 : y = bernstein(q, i, x)
146 : y = bernstein(Val(q), i, x)
147 :
148 : Evaluate `i`-th (univariate) Bernstein polynomial of degree `q` at `x` for
149 : `i in 0:q`.
150 :
151 : ys = bernstein(q, x)
152 : ys = bernstein(Val(q), x)
153 :
154 : @assert sum(ys) ≈ 1
155 :
156 : Evaluate all Bernstein polynomials of degree `q` at `x` to tuple.
157 :
158 : See also [`mvbernstein`](@ref), [`monomial`](@ref)
159 : """ bernstein
160 :
161 : """
162 : d = domainpoints(d, q)
163 :
164 : Generate sequence of domain points of `d`-variate Bernstein
165 : polynomials of total degree `q` in *barycentric coordinates* (w.r.t. a
166 : `d`-simplex).
167 :
168 : See also [`multiindices`](@ref), [`mvbernstein`](@ref)
169 : """
170 4 : domainpoints(d, q) = map(i -> SVector(i ./ q), multiindices(d, q))
171 :
172 2 : function _exmvb(d, q, typ)
173 9 : function terms(idx)
174 7 : ts = map(1:length(idx)) do i
175 17 : quote
176 187 : x[$i] ^ $(idx[i])
177 : end
178 : end
179 :
180 : # try to avoid overflow
181 14 : scale = factorial(q) // 1
182 7 : for i in idx
183 26 : scale *= 1 // factorial(i)
184 24 : end
185 :
186 7 : Expr(:call, :prod, Expr(:vect, typ(scale), ts...))
187 : end
188 :
189 2 : [ terms(idx) for idx in multiindices(d, q) ]
190 : end
191 :
192 13 : @generated function mvbernstein(::Val{Q}, x::NTuple{D, T}) where {Q, D, T}
193 2 : quote
194 13 : tuple($(_exmvb(D, Q, T)...))
195 : end
196 : end
197 :
198 1 : mvbernstein(q::Val{Q}, x::SVector) where {Q} = mvbernstein(q, tuple(x...))
199 0 : mvbernstein(q::Val{Q}, x...) where {Q} = mvbernstein(q, tuple(x...))
200 1 : mvbernstein(q::Val{Q}, x::Real) where {Q} = mvbernstein(q, tuple(1 - x, x))
201 :
202 13 : mvbernstein(q::Int, args...) = mvbernstein(Val(q), args...)
203 :
204 : """
205 : values = mvbernstein(Val(Q), λ₁, λ₂, ...)
206 : values = mvbernstein(Val(Q), λ::NTuple)
207 : values = mvbernstein(Val(Q), λ::SVector)
208 : values = mvbernstein(q, ...)
209 :
210 : Evaluate multivariate Bernstein polynomials of total degree `q` at
211 : *barycentric coordinates* `λ`. The sequence of basis functions
212 : correspond to [`multiindices`](@ref) or [`domainpoints`](@ref),
213 : likewise.
214 :
215 : !!! note
216 : In contrast to [`mvmonomials`](@ref), `mvbernstein` takes
217 : **barycentric coordinates** `λ` as input! Prefer passing a
218 : *tuple* or an *SVector* (with `sum(λ) = 1`).
219 :
220 : As an exception, univariate Bernstein polynomials can be evaluated
221 : at `x::Real` as
222 :
223 : values = mvbernstein(q, x)
224 :
225 : This call is equal to
226 :
227 : values = mvbernstein(q, 1-x, x) # pass barycentric coordinates
228 :
229 : and to
230 :
231 : values = bernstein(q, x)
232 :
233 : See also [`mvbernstein!`](@ref), [`domainpoints`](@ref),
234 : [`mvmonomials`](@ref), [`bernstein`](@ref), [`mvpdim`](@ref)
235 : """ mvbernstein
236 :
237 3 : mvbernstein!(q, v::AbstractVector, x::NTuple{D, T}) where {D, T} = (v.= mvbernstein(q, x))
238 :
239 1 : mvbernstein!(q, v::AbstractVector, x::SVector) = mvbernstein!(q, v, tuple(x...))
240 1 : mvbernstein!(q, v::AbstractVector, x...) = mvbernstein!(q, v, tuple(x...))
241 0 : mvbernstein!(q, v::AbstractVector, x::Real) = mvbernstein!(q, v, (1 - x, x))
242 :
243 : """
244 : mvbernstein!(Val(Q), v, λ...)
245 : mvbernstein!(q, v, λ...)
246 :
247 : Evaluate the multivariate Bernstein polynomials of total degree `Q` or
248 : `q` at *barycentric coordinates* `λ` as for [`mvbernstein`](@ref) and
249 : store the computed values in `v`.
250 :
251 : !!! warning "precondition"
252 : `v` must have the appropriate size [`mvpdim`](@ref).
253 :
254 : !!! note
255 : In contrast to [`mvmonomials`](@ref), `mvbernstein` takes
256 : **barycentric coordinates** `λ` as input! Prefer passing a
257 : *tuple* or an *SVector* (with `sum(λ) = 1`).
258 :
259 : As an exception, univariate Bernstein polynomials can be evaluated
260 : at `x::Real` as
261 :
262 : values = mvbernstein!(q, x)
263 :
264 : This call is equal to
265 :
266 : values = mvbernstein!(q, 1-x, x) # pass barycentric coordinates
267 :
268 : See also [`mvbernstein`](@ref), [`bvmdim`](@ref)
269 : """ mvbernstein!
270 :
271 : #------------------------------------------------------------------------------
272 :
273 : # NOTE: need to parameterize differently:
274 : # - polynomial space (defines, Q, D, N)
275 : # - basis
276 : # -
277 : # ...
278 :
279 : """
280 : Space of `D`-variate polynomials of total degree `Q`
281 : """
282 6 : struct PolynomialSpace{D, Q}
283 : end
284 :
285 5 : polynomialspace(d, q) = polynomialspace(_as_val(d), _as_val(q))
286 5 : polynomialspace(::Val{D}, ::Val{Q}) where {D, Q} = PolynomialSpace{D, Q}()
287 :
288 : """
289 : p = polynomialspace(d, q)
290 : p = polynomialspace(b::PolynomialBasis)
291 :
292 : Construct or query [`PolynomialSpace`](@ref)
293 :
294 : See also [`PolynomialBasis`](@ref)
295 : """ polynomialspace
296 :
297 : """
298 : q degree(p::PolynomialSpace)
299 :
300 : Get (total) degree of [`PolynomialSpace`](@ref)
301 : """
302 1 : degree(::PolynomialSpace{D, Q}) where {D, Q} = Q
303 :
304 : """
305 : d = nvariate(p::PolynomialSpace)
306 :
307 : Get number of variables in`p`.
308 : """
309 1 : nvariate(::PolynomialSpace{D, Q}) where {D, Q} = D
310 :
311 : """
312 : n = dimension(p::PolynomialSpace)
313 :
314 : Get dimensions (degrees of freedom) of [`PolynomialSpace`](@ref)
315 :
316 : See also [`mvpdim`](@ref)
317 : """
318 1 : dimension(::PolynomialSpace{D, Q}) where {D, Q} = mvpdim(D, Q)
319 :
320 : """
321 : Basis of a [`PolynomialSpace`](@ref)
322 : """
323 : struct PolynomialBasis{D, Q, B}
324 1 : basis::B
325 : end
326 :
327 : """
328 : values = (::PolynomialBasis)(xs...)
329 :
330 : Evaluate [`PolynomialBasis`](@ref).
331 :
332 : See also [`mvmonomials`](@ref), [`mvbernstein`](@ref)
333 : """
334 0 : (b::PolynomialBasis{D, Q})(args...) where {D, Q} = b.basis(Val(Q), args...)
335 :
336 1 : polynomialbasis(b::B, ::PolynomialSpace{D, Q}) where {B, D, Q} =
337 : PolynomialBasis{D, Q, B}(b)
338 :
339 0 : polynomialbasis(b::B, d, q) where {B} = polynomialbasis(b, polynomialspace(d, q))
340 :
341 : """
342 : b = polynomialbasis(f, p::PolynomialSpace)
343 : b = polynomialbasis(f, d, q)
344 :
345 : Construct basis for [`PolynomialSpace`](@ref), where `f` evaluates basis functions
346 : such as [`mvmonomials`](@ref) or [`mvbernstein`](@ref).
347 :
348 : See also [`PolynomialBasis`](@ref), [`PolynomialSpace`](@ref)
349 : """ polynomialbasis
350 :
351 1 : polynomialspace(::PolynomialBasis{D, Q, B}) where {D, Q, B} =
352 : PolynomialSpace{D, Q}()
353 :
354 0 : _scalartype(::Type{T}) where {T<:Real} = T
355 0 : _scalartype(::T) where {T<:Real} = T
356 :
357 0 : function _undef_ata(typ, p::PolynomialSpace{D, Q}) where {D, Q}
358 0 : N = dimension(p)
359 0 : T = _scalartype(typ)
360 :
361 0 : MMatrix{N, N, T, N * N}(undef)
362 : end
363 :
364 0 : function _undef_atb(typ, nrhs, p::PolynomialSpace{D, Q}) where {D, Q}
365 0 : N = dimension(p)
366 0 : T = _scalartype(typ)
367 0 : NRHS = _as_int(nrhs)
368 :
369 0 : MMatrix{N, NRHS, T, N * NRHS}(undef)
370 : end
371 :
372 0 : function _undef_atb(typ, ::AbstractVector, p::PolynomialSpace{D, Q}) where {D, Q}
373 0 : N = dimension(p)
374 0 : T = _scalartype(typ)
375 :
376 0 : MVector{N, T}(undef)
377 : end
378 :
379 0 : function _undef_atb(typ, B::AbstractMatrix, p::PolynomialSpace{D, Q}) where {D, Q}
380 0 : N = dimension(p)
381 0 : T = _scalartype(typ)
382 :
383 0 : Matrix{T}(undef, N, size(B, 2))
384 : end
385 :
386 : # more...?
387 :
388 0 : function fit(basis::PolynomialBasis{B, D, Q}, ts, bs) where {B, D, Q}
389 0 : sampletyp = eltype(ts) # expect Real, Tuple, SVector
390 0 : T = eltype(sampletyp)
391 :
392 0 : p = polynomialspace(basis)
393 0 : N = dimension(p)
394 0 : AtA = _undef_ata(T, p)
395 0 : Atb = _undef_atb(T, bs, p)
396 :
397 0 : AtA .= T(0)
398 0 : Atb .= T(0)
399 : # u = MMatrix{N, 1, T, N}(undef)
400 :
401 : # TODO: weights type
402 :
403 : # _weight(::Nothing, i) = T(1), T(1)
404 : # _weight(w, i) = w[i], sqrt(w[i])
405 : # _weight(i) = _weight(ws, i)
406 :
407 0 : for (ti, bi) in zip(ts, bs) # eachrow(bs)) # for a vectors it's not eachrow... # TODO: row must be an svector !!!
408 0 : w, sqrtw = 1, 1 # _weight(i)
409 :
410 0 : u = SVector(basis(ti...)...)
411 :
412 0 : AtA .+= (w * u) * u' # TODO: rankupdate
413 0 : Atb .+= u * (sqrtw * bi)
414 0 : end
415 :
416 0 : AtA, Atb
417 : end
418 :
419 : # TODO: how to B?
420 : # TODO: in-place
421 :
422 : # rename: assemble (A,b) (AtA, Atb), fit; given samples and weights
423 :
424 : # setup A'A (rank update) -> MMatrix and MVector and cholesky! and lowrankupdate!
425 :
426 : # TODO: tests, benchmarks !!!
427 :
428 :
429 : # @btime (A = BaseUtilities.assemblemvp(BaseUtilities.mvmonomials, Val(2), [2.0, 3.0, 5.0]); (A'A, A'b)
430 : # @btime BaseUtilities.assemblemvp(mvmonomials, Val(2), [2.0, 3.0, 5.0], [1.0, 2.0, 3.0])
431 :
432 : # TODO: wrap state in a struct; cholesky
433 :
434 : # struct MVPFit{B, D, Q, T, N, W, NRHS, NN, NNRHS}
435 : # AtA::MMatrix{N, N, T, NN}
436 : # Atb::MMatrix{N, NRHS, T, NNRHS}
437 : # weights::W
438 : # basis::B
439 :
440 : # function MVPFit{B, D, Q, T, N, W, NRHS, NN, NNRHS}(basis::B, d::Val{D}, q::Val{Q},
441 : # ::Type{T}, ::Val{N}, w::W, ::Val{NRHS},
442 : # ::Val{NN}, ::Val{NNRHS}) where {B, D, Q, T, N, W, NRHS, NN, NNRHS}
443 : # new(MMatrix{N, N, T, N*N}(undef), Atb::MMatrix{N, NRHS, T, N * NRHS}, weights, basis)
444 : # end
445 : # end
446 :
447 : # function mvpfit(basis::B, q, xs, rhs, weights = nothing) where {B}
448 : # vtyp = eltype(xs) # expect Real, Tuple, SVector
449 : # styp = eltype(vtyp)
450 :
451 : # _len(::NTuple{N, T}) where {N, T} = Val(N)
452 : # _len(::SVector{N, T}) where {N, T} = Val(N)
453 : # _len(::Type{T}) where {T<:Real} = Val(1)
454 :
455 : # N = mvpdim(q, _len(vtyp))
456 : # NRHS = _len(eltype(rhs
457 :
458 : # MVPFit{B, _as_val(q),
459 : # end
|
