gram schmidt orthogonalization

Author: daanhb

src/BasisFunctions.jl

@@ -322,6 +322,8 @@ include("sampling/quadweights.jl")
 
 include("bases/dictionary/indexing.jl")
 include("bases/dictionary/dictionary.jl")
+include("bases/dictionary/dict_evaluation.jl")
+include("bases/dictionary/dict_moments.jl")
 include("bases/dictionary/promotion.jl")
 include("bases/dictionary/span.jl")
 include("bases/dictionary/orthog.jl")

src/bases/dictionary/dictionary.jl

@@ -327,8 +327,3 @@ function iterate(d::Dictionary, state)
         (d[next_item], (iter,next_tuple))
     end
 end
-
-
-
-include("dict_evaluation.jl")
-include("dict_moments.jl")

src/bases/dictionary/expansions.jl

@@ -217,3 +217,6 @@ expansion_of_one(dict::Dictionary) = Expansion(dict, coefficients_of_one(dict))
 expansion_of_x(dict::Dictionary) = Expansion(dict, coefficients_of_x(dict))
 
 apply(op::DictionaryOperator, e::Expansion) = Expansion(dest(op), op * coefficients(e))
+
+norm(f::Expansion) = norm(f, measure(f))
+norm(f::Expansion, μ) = sqrt(innerproduct(f, f, μ))

src/bases/dictionary/orthog.jl

@@ -278,3 +278,28 @@ function mixedgrammatrix!(G, Φ1::Dictionary, Φ2::Dictionary, μ; options...)
     end
     G
 end
+
+
+#################################
+# Gram-Schmidt orthogonalization
+#################################
+
+
+function gram_schmidt(U::Dictionary, μ = measure(U))
+    N = length(U)
+    V = [Expansion(U, zeros(U)) for k in 1:N]
+    V[1].coefficients[1] = 1 / norm(U[1], μ)
+    for k in 2:N
+        V[k].coefficients[k] = 1
+        for j in 1:k-1
+            V[k] = V[k] - innerproduct(V[k], V[j], μ) * V[j]
+        end
+        V[k] = V[k] / norm(V[k], μ)
+    end
+
+    A = UpperTriangular(zeros(coefficienttype(U), N, N))
+    for k in 1:N
+        A[1:k,k] = V[k].coefficients[1:k]
+    end
+    A*U
+end

src/bases/poly/polynomials.jl

@@ -57,3 +57,20 @@ degree(p::Polynomial) = p.degree
 
 support(p::Polynomial) = support(dictionary(p))
 measure(p::Polynomial) = measure(dictionary(p))
+
+########################
+# Polynomial expansions
+########################
+
+"""
+    polynomialdegree(f::Expansion)
+
+If `f` is an expansion in a polynomial basis, return the degree of the polynomial function (not taking
+into account possible leading zero coefficients).
+"""
+polynomialdegree(f::Expansion) = maxpolynomialdegree(dictionary(f))
+
+"The maximal polynomial degree of any expansion in the given basis."
+maxpolynomialdegree(b::PolynomialBasis) = PolynomialDegree(length(b)-1)
+
+maxpolynomialdegree(b::OperatedDict) = maxpolynomialdegree(superdict(b))

src/computations/approximation.jl

@@ -3,6 +3,8 @@
 # Generic approximation
 ########################
 
+const ANYFUNCTION = Union{Function,Expansion,TypedFunction}
+
 """
 The `approximation` function returns an operator that can be used to approximate
 a function in the function set. This operator maps a grid to a set of coefficients.
@@ -23,27 +25,14 @@ end
 continuous_approximation(Φ::Dictionary; options...) = inv(gram(Φ; options...))
 
 # Automatically sample a function if an operator is applied to it
-(*)(op::DictionaryOperator, f::Function) = op * project(src(op), f)
-Base.:\(op::DictionaryOperator, f::Function) = op \ project(dest(op), f)
+(*)(op::DictionaryOperator, f::ANYFUNCTION) = op * project(src(op), f)
+Base.:\(op::DictionaryOperator, f::ANYFUNCTION) = op \ project(dest(op), f)
 
-project(Φ::GridBasis, f::Function; options...) = sample(Φ, f)
+project(Φ::GridBasis, f::ANYFUNCTION; options...) = sample(Φ, f)
 
 approximate(Φ::Dictionary, f::Function; options...) =
     Expansion(Φ, approximation(Φ; options...) * f)
 
-"""
-The 2-argument approximation exists to allow you to transform any
-Dictionary coefficients to any other Dictionary coefficients, including
-Discrete Grid Spaces.
-If a transform exists, the transform is used.
-"""
-function approximation(A::Dictionary, B::Dictionary, options...)
-    if hastransform(A, B)
-        return transform(A, B; options...)
-    else
-        error("Don't know a transform from ", A, " to ", B)
-    end
-end
 
 
 #################
@@ -78,7 +67,7 @@ default_interpolation(Φ::Dictionary, gb::GridBasis; options...) =
 default_interpolation(::Type{T}, Φ::Dictionary, gb::GridBasis; options...) where {T} =
     QR_solver(evaluation(T, Φ, grid(gb)))
 
-interpolate(Φ::Dictionary, pts, f, T=coefficienttype(s)) = Expansion(interpolation(T, Φ, pts) * f)
+interpolate(Φ::Dictionary, pts, f, T=coefficienttype(Φ)) = Expansion(Φ,interpolation(T, Φ, pts) * f)
 
 """
 Compute the weights of an interpolatory quadrature rule.