From 2016b3f8900c42e816f6c56f19328d1e0ee07f60 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Wed, 10 Jul 2024 11:12:04 +0200
Subject: [PATCH] up

---
 src/MultivariateBases.jl |  1 +
 src/chebyshev.jl         | 21 +++++++++-----
 src/lagrange.jl          | 61 ++++++++++++++++++++++++++++++++++++++--
 src/orthogonal.jl        | 15 ++++++++++
 4 files changed, 88 insertions(+), 10 deletions(-)

diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl
index c03ac5a..fcef0f7 100644
--- a/src/MultivariateBases.jl
+++ b/src/MultivariateBases.jl
@@ -64,6 +64,7 @@ include("hermite.jl")
 include("laguerre.jl")
 include("legendre.jl")
 include("chebyshev.jl")
+include("lagrange.jl")
 include("quotient.jl")
 
 function algebra(
diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index f94ab74..24aa106 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -89,6 +89,17 @@ function SA.coeffs!(
     return res
 end
 
+function transformation_to(
+    source::SubBasis{Chebyshev},
+    target::SubBasis{Monomial},
+)
+    A = zeros(Float64, length(target), length(source))
+    for (i, cheby) in enumerate(source)
+        A[:, i] = SA.coeffs(algebra_element(cheby), target)
+    end
+    return A
+end
+
 function SA.coeffs(
     cfs,
     source::SubBasis{Monomial},
@@ -97,17 +108,13 @@ function SA.coeffs(
     sub = explicit_basis_covering(target, source)
     # Need to make A square so that it's UpperTriangular
     extended = SubBasis{Monomial}(sub.monomials)
-    A = zeros(Float64, length(extended), length(sub))
-    for (i, cheby) in enumerate(sub)
-        A[:, i] = SA.coeffs(algebra_element(cheby), extended)
-    end
     ext = SA.coeffs(algebra_element(cfs, source), extended)
     return SA.SparseCoefficients(
         sub.monomials,
-        #LinearAlgebra.UpperTriangular(A) \ ext, # Julia v1.6 converts `A` to the eltype of the `result` which is bad for JuMP
+        #transformation_to(sub, extended) \ ext, # Julia v1.6 converts the matrix to the eltype of the `result` which is bad for JuMP
         LinearAlgebra.ldiv!(
-            zeros(_promote_coef(eltype(ext), Chebyshev), size(A, 2)),
-            LinearAlgebra.UpperTriangular(A),
+            zeros(_promote_coef(eltype(ext), Chebyshev), length(sub)),
+            transformation_to(sub, extended),
             ext,
         ),
     )
diff --git a/src/lagrange.jl b/src/lagrange.jl
index 744c61f..4788177 100644
--- a/src/lagrange.jl
+++ b/src/lagrange.jl
@@ -1,10 +1,65 @@
-struct ImplicitLagrangeBasis
+# Inspired by `Hypatia.jl/src/PolyUtils/realinterp.jl`
+
+export BoxSampling
+
+abstract type AbstractNodes{T,V} end
+
+mutable struct BoxSampling{T,V} <: AbstractNodes{T,V}
+    lower::V
+    upper::V
+    sample_factor::Int
+    orthogonalize::Bool
+    function BoxSampling(lower, upper; sample_factor = 0, orthogonalize = false)
+        @assert length(lower) == length(upper)
+        l = float.(lower)
+        u = float.(upper)
+        V = promote_type(typeof(l), typeof(u))
+        return new{eltype(V),V}(l, u, sample_factor, orthogonalize)
+    end
 end
 
-struct LagrangePolynomial{T,V<:AbstractVector{T}}
+function sample(s::BoxSampling, n::Integer) where {T}
+    samples = rand(T, length(s.lower), n) .- inv(T(2))
+    shift = (dom.u .+ dom.l) .* inv(T(2))
+    for i in 1:n
+        for j in eachindex(s.lower)
+            samples[j, i] = samples[j, i] * (dom.u[j] - dom.l[j]) + shift[j]
+        end
+    end
+    return samples
+end
+
+struct LagrangePolynomial{T,V}
     point::V
 end
 
-struct LagrangeBasis{T,P<:AbstractVector{T},V<:AbstractVector{P}}
+struct ImplicitLagrangeBasis{T,V,N<:AbstractNodes{T,V}} <:
+       SA.ImplicitBasis{LagrangePolynomial{T,V},V}
+    sampling::AbstractNodes{T,V}
+    function ImplicitLagrangeBasis(nodes::AbstractNodes{T,V}) where {T,V}
+        return new{T,V,typeof(nodes)}(nodes)
+    end
+end
+struct LagrangeBasis{T,P,V<:AbstractVector{P}} <:
+       SA.ExplicitBasis{LagrangePolynomial{T,V},Int}
     points::V
 end
+
+Base.length(basis::LagrangeBasis) = length(basis.points)
+
+function transformation_to(basis::SubBasis, lag::LagrangeBasis{T}) where {T}
+    V = Matrix{T}(undef, length(lag), length(basis))
+    for 
+end
+
+function sample(s::AbstractNodes, basis::SubBasis)
+    samples = sample(s, length(basis) * s.sample_factor)
+    return eachcol(samples)
+end
+
+function explicit_basis_covering(
+    implicit::ImplicitLagrangeBasis,
+    basis::SubBasis,
+)
+    return LagrangeBasis(sample(implicit.sampling, length(basis)))
+end
diff --git a/src/orthogonal.jl b/src/orthogonal.jl
index 9f5e0c6..093b416 100644
--- a/src/orthogonal.jl
+++ b/src/orthogonal.jl
@@ -62,6 +62,21 @@ d_k p_k(x_i) = (a_k x_i + b_k) p_{k-1}(x_i) + c_k p_{k-2}(x_i)
 """
 function reccurence_deno_coef end
 
+function univariate_eval!(
+    values::Vector,
+    value,
+    ::Type{B},
+) where {B}
+    if 1 in eachindex(values)
+        values[1] = one(value)
+    end
+    if 2 in eachindex(values)
+        values[1] = one(value)
+    end
+    for d in 2:(length(values) - 1)
+    end
+end
+
 """
     univariate_orthogonal_basis(
         B::Type{<:AbstractMultipleOrthogonal},