Pointwise dot() for powerlaw Legendre (#34)

* move over legendredot code

* generic types

* remove todo

* minor tidying

* minor tidying

* Update test_fourier.jl

* Fix broadcasting notation

* Try avoiding Clenshaw

* big stability overhaul

* minor adjustments

* add dot()

* adjust dot

* Fix inconsistent resizedata!

* generalised weighted OP re-expansion

* Fix \(::Weighted, ::Weighted)

* Tests pass

* v0.4.1

* trigger build

* increase coverage

Co-authored-by: Sheehan Olver <[email protected]>

Author: TSGut

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.4"
+version = "0.4.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,7 +29,7 @@ ArrayLayouts = "0.7"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.15"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.8"
+ContinuumArrays = "0.8.1"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
@@ -42,7 +42,7 @@ IntervalSets = "0.5"
 LazyArrays = "0.21"
 LazyBandedMatrices = "0.5.5"
 QuasiArrays = "0.6"
-SpecialFunctions = "1"
+SpecialFunctions = "1.0"
 julia = "1.6"
 
 [extras]

src/ClassicalOrthogonalPolynomials.jl

@@ -1,4 +1,5 @@
 module ClassicalOrthogonalPolynomials
+using IntervalSets: UnitRange
 using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, BlockArrays,
     IntervalSets, DomainSets, ArrayLayouts, SpecialFunctions,
     InfiniteLinearAlgebra, InfiniteArrays, LinearAlgebra, FastGaussQuadrature, FastTransforms, FFTW,
@@ -11,8 +12,8 @@ import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout,
                 sub_materialize, arguments, sub_paddeddata, paddeddata, PaddedLayout, resizedata!, LazyVector, ApplyLayout, call,
-                _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout, 
-                AbstractCachedVector, AbstractCachedMatrix
+                _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout,
+                AbstractCachedVector, AbstractCachedMatrix, paddeddata
 import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
                         subdiagonaldata, diagonaldata, supdiagonaldata
 import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, AbstractLazyBandedLayout
@@ -39,20 +40,17 @@ import FastGaussQuadrature: jacobimoment
 import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray
 import BandedMatrices: bandwidths
 
-export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial, 
+export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial,
             Hermite, Jacobi, Legendre, Chebyshev, ChebyshevT, ChebyshevU, ChebyshevInterval, Ultraspherical, Fourier, Laguerre,
             HermiteWeight, JacobiWeight, ChebyshevWeight, ChebyshevGrid, ChebyshevTWeight, ChebyshevUWeight, UltrasphericalWeight, LegendreWeight, LaguerreWeight,
             WeightedUltraspherical, WeightedChebyshev, WeightedChebyshevT, WeightedChebyshevU, WeightedJacobi,
-            ∞, Derivative, .., Inclusion, 
+            ∞, Derivative, .., Inclusion,
             chebyshevt, chebyshevu, legendre, jacobi,
             legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
             jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted
 
-if VERSION < v"1.6-"
-    oneto(n) = Base.OneTo(n)
-else
-    import Base: oneto
-end
+
+import Base: oneto
 
 
 include("interlace.jl")
@@ -278,10 +276,10 @@ end
 
 _tritrunc(X, n) = _tritrunc(MemoryLayout(X), X, n)
 
-jacobimatrix(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Inclusion,OneTo}}) = 
+jacobimatrix(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Inclusion,OneTo}}) =
     _tritrunc(jacobimatrix(parent(V)), maximum(parentindices(V)[2]))
 
-grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,AbstractUnitRange}}) = 
+grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,AbstractUnitRange}}) =
     eigvals(symtridiagonalize(jacobimatrix(P)))
 
 function golubwelsch(X)
@@ -326,11 +324,11 @@ function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}},
     parent(A) \ parent(B)
 end
 
+# assume we can expand w_B in wA to reduce to polynomial multiplication
 function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
-    w_A,A = arguments(wA)
+    _,A = arguments(wA)
     w_B,B = arguments(wB)
-    w_A == w_B || error("Not implemented")
-    A\B
+    A \ ((A * (wA \ w_B)) .* B)
 end
 
 

src/classical/jacobi.jl

@@ -177,11 +177,13 @@ broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::HalfWeighted
 \(w_A::HalfWeighted, B::AbstractQuasiArray) = convert(WeightedOrthogonalPolynomial, w_A) \ B
 \(A::AbstractQuasiArray, w_B::HalfWeighted) = A \ convert(WeightedOrthogonalPolynomial, w_B)
 
-function \(A::AbstractQuasiArray, w_B::WeightedOrthogonalPolynomial{<:Any,<:Weight,<:Normalized})
+function _norm_expand_ldiv(A, w_B)
     w,B = w_B.args
     B̃,D = arguments(ApplyLayout{typeof(*)}(), B)
     (A \ (w .* B̃)) * D
 end
+\(A::AbstractQuasiArray, w_B::WeightedOrthogonalPolynomial{<:Any,<:Weight,<:Normalized}) = _norm_expand_ldiv(A, w_B)
+\(A::WeightedOrthogonalPolynomial, w_B::WeightedOrthogonalPolynomial{<:Any,<:Weight,<:Normalized}) = _norm_expand_ldiv(A, w_B)
 
 axes(::AbstractJacobi{T}) where T = (Inclusion{T}(ChebyshevInterval{real(T)}()), oneto(∞))
 ==(P::Jacobi, Q::Jacobi) = P.a == Q.a && P.b == Q.b
@@ -292,6 +294,17 @@ function recurrencecoefficients(P::Jacobi)
     (A,B,C)
 end
 
+# explicit special case for normalized Legendre
+# todo: do we want these explicit constructors for normalized Legendre?
+# function jacobimatrix(::Normalized{<:Any,<:Legendre{T}}) where T
+#     b = (one(T):∞) ./sqrt.(4 .*(one(T):∞).^2 .-1)
+#     Symmetric(_BandedMatrix(Vcat(zeros(∞)', (b)'), ∞, 1, 0), :L)
+# end
+# function recurrencecoefficients(::Normalized{<:Any,<:Legendre{T}}) where T
+#     n = zero(T):∞
+#     nn = one(T):∞
+#     ((2n .+ 1) ./ (n .+ 1) ./ sqrt.(1 .-2 ./(3 .+2n)), Zeros{T}(∞), Vcat(zero(T),nn ./ (nn .+ 1) ./ sqrt.(1 .-4 ./(3 .+2nn))))
+# end
 
 @simplify *(X::Identity, P::Legendre) = ApplyQuasiMatrix(*, P, P\(X*P))
 

src/lanczos.jl

@@ -49,7 +49,7 @@ function LanczosData(w::AbstractQuasiVector, P::AbstractQuasiMatrix)
     x = axes(P,1)
     wP = weighted(P)
     X = jacobimatrix(P)
-    W = Clenshaw(P * (wP \ w), P)
+    W = wP \ (w .* P)
     LanczosData(X, W)
 end
 
@@ -58,6 +58,7 @@ function resizedata!(L::LanczosData, N)
     resizedata!(L.R, N, N)
     resizedata!(L.γ, N)
     resizedata!(L.β, N)
+    resizedata!(L.W, N, N)
     lanczos!(L.ncols+1:N, L.X, L.W, L.γ, L.β, L.R)
     L.ncols = N
     L

src/normalized.jl

@@ -166,12 +166,27 @@ broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::OrthonormalW
 
 
 abstract type AbstractWeighted{T} <: Basis{T} end
+struct WeightedOPLayout <: AbstractBasisLayout end
 
-MemoryLayout(::Type{<:AbstractWeighted}) = WeightedBasisLayout()
+# make act like WeightedBasisLayout
+ContinuumArrays._grid(::WeightedOPLayout, P) = ContinuumArrays._grid(WeightedBasisLayout(), P)
+ContinuumArrays._factorize(::WeightedOPLayout, P) = ContinuumArrays._factorize(WeightedBasisLayout(), P)
+ContinuumArrays.sublayout(::WeightedOPLayout, inds::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractVector}}) = sublayout(WeightedBasisLayout(), inds)
+ContinuumArrays.sublayout(::WeightedOPLayout, inds::Type{<:Tuple{<:Inclusion,<:AbstractVector}}) = sublayout(WeightedBasisLayout(), inds)
+
+MemoryLayout(::Type{<:AbstractWeighted}) = WeightedOPLayout()
 ContinuumArrays.unweightedbasis(wP::AbstractWeighted) = wP.P
-\(w_A::AbstractWeighted, w_B::AbstractWeighted) = convert(WeightedOrthogonalPolynomial, w_A) \ convert(WeightedOrthogonalPolynomial, w_B)
-\(w_A::AbstractWeighted, B::AbstractQuasiArray) = convert(WeightedOrthogonalPolynomial, w_A) \ B
-\(A::AbstractQuasiArray, w_B::AbstractWeighted) = A \ convert(WeightedOrthogonalPolynomial, w_B)
+function copy(L::Ldiv{WeightedOPLayout,WeightedOPLayout})
+    L.A.P == L.B.P && return Eye{eltype(L)}(∞)
+    convert(WeightedOrthogonalPolynomial, L.A) \ convert(WeightedOrthogonalPolynomial, L.B)
+end
+
+copy(L::Ldiv{WeightedOPLayout,<:AbstractLazyLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
+copy(L::Ldiv{WeightedOPLayout,<:AbstractBasisLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
+copy(L::Ldiv{<:AbstractLazyLayout,WeightedOPLayout}) = L.A \ convert(WeightedOrthogonalPolynomial, L.B)
+copy(L::Ldiv{<:AbstractBasisLayout,WeightedOPLayout}) = L.A \ convert(WeightedOrthogonalPolynomial, L.B)
+copy(L::Ldiv{WeightedOPLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
+copy(L::Ldiv{WeightedOPLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
 
 """
     Weighted(P)
@@ -186,7 +201,6 @@ axes(Q::Weighted) = axes(Q.P)
 copy(Q::Weighted) = Q
 
 ==(A::Weighted, B::Weighted) = A.P == B.P
-
 weight(wP::Weighted) = orthogonalityweight(wP.P)
 
 convert(::Type{WeightedOrthogonalPolynomial}, P::Weighted) = weight(P) .* unweightedbasis(P)

src/stieltjes.jl

@@ -175,8 +175,6 @@ end
     T[kr,:] * Diagonal(Vcat(2*convert(V,π)*c*log(c), 2c*D.diag.args[2]))
 end
 
-
-
 ### generic fallback
 for Op in (:Hilbert, :StieltjesPoint, :LogKernel, :PowKernel)
     @eval @simplify function *(H::$Op, wP::WeightedBasis{<:Any,<:Weight,<:Any}) 
@@ -185,3 +183,200 @@ for Op in (:Hilbert, :StieltjesPoint, :LogKernel, :PowKernel)
         (H * Weighted(Q)) * (Q \ P)
     end
 end
+
+
+#################################################
+# ∫f(x)g(x)(t-x)^a dx evaluation where f and g given in coefficients
+#################################################
+# recognize structure of W = ((t .- x).^a
+const PowKernelPoint{T,V,D,F} =  BroadcastQuasiVector{T, typeof(^), Tuple{ContinuumArrays.AffineQuasiVector{T, V, Inclusion{V, D}, T}, F}}
+
+####
+# cached operator implementation
+####
+# Constructors support BigFloat and it's recommended to use them for high orders.
+mutable struct PowerLawMatrix{T, PP<:Normalized{<:Any,<:Legendre{<:Any}}} <: AbstractCachedMatrix{T}
+    P::PP # OPs - only normalized Legendre supported for now
+    a::T  # naming scheme follows (t-x)^a
+    t::T
+    data::Matrix{T}
+    datasize::Tuple{Int,Int}
+    function PowerLawMatrix{T, PP}(P::PP, a::T, t::T) where {T, PP<:AbstractQuasiMatrix}
+        new{T, PP}(P,a,t, gennormalizedpower(a,t,50),(50,50))
+    end
+end
+PowerLawMatrix(P::AbstractQuasiMatrix, a::T, t::T) where T = PowerLawMatrix{T,typeof(P)}(P,a,t)
+size(K::PowerLawMatrix) = (ℵ₀,ℵ₀) # potential to add maximum size of operator
+copy(K::PowerLawMatrix{T,PP}) where {T,PP} = K # Immutable entries
+
+# data filling
+#TODO: fix the weird inds
+cache_filldata!(K::PowerLawMatrix, inds, _) = fillcoeffmatrix!(K, inds)
+
+# because it really only makes sense to compute this symmetric operator in square blocks, we have to slightly rework some of LazyArrays caching and resizing
+function getindex(K::PowerLawMatrix, k::Int, j::Int)
+    resizedata!(K, k, j)
+    K.data[k, j]
+end
+function getindex(K::PowerLawMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange)
+    resizedata!(K, maximum(kr),maximum(jr))
+    K.data[kr, jr]
+end
+function resizedata!(K::PowerLawMatrix, n::Integer, m::Integer)
+    olddata = K.data
+    νμ = size(olddata)
+    nm = max.(νμ,max(n,m))
+    if νμ ≠ nm
+        K.data = similar(K.data, nm...)
+        K.data[axes(olddata)...] = olddata
+    end
+    if maximum(nm) > maximum(νμ)
+        inds = maximum(νμ):maximum(nm)
+        cache_filldata!(K, inds, inds)
+        K.datasize = nm
+    end
+    K
+end
+
+####
+# methods
+####
+function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), K::PowKernelPoint, Q::Normalized{<:Any,<:Legendre})
+    tx,a = K.args
+    t = tx[zero(typeof(a))]
+    return Q*PowerLawMatrix(Q,a,t)
+end
+
+####
+# operator generation
+####
+
+# evaluates ∫(t-x)^a Pn(x)P_m(x) dx for the case m=0, i.e. first row of operator, via backwards recursion
+function powerleg_backwardsfirstrow(a::T, t::T, ℓ::Int) where T <: Real
+    ℓ = ℓ+200
+    coeff = zeros(BigFloat,ℓ)
+    coeff[end-1] =  one(BigFloat)
+    for m = reverse(1:ℓ-2)
+        coeff[m] = (coeff[m+2]-t/((a+m+2)/(2*m+1))*normconst_Pnadd1(m,t)*coeff[m+1])/(((a-m+1)/(2*m+1))/((a+m+2)/(2*m+1))*normconst_Pnsub1(m,t))
+    end
+    coeff = PLnorminitial00(t,a)/coeff[1].*coeff
+    return T.(coeff[1:ℓ-200])
+end
+# modify recurrence coefficients to work for normalized Legendre
+normconst_Pnadd1(m::Int, settype::T) where T<:Real = sqrt(2*m+3*one(T))/sqrt(2*m+one(T))
+normconst_Pnsub1(m::Int, settype::T) where T<:Real = sqrt(2*m+3*one(T))/sqrt(2*m-one(T))
+normconst_Pmnmix(n::Int, m::Int, settype::T) where T<:Real = sqrt(2*m+3*one(T))*sqrt(2*n+one(T))/(sqrt(2*m+one(T))*sqrt(2*n-one(T)))
+# useful explicit initial case
+function PLnorminitial00(t::Real, a::Real)
+    return ((t+1)^(a+1)-(t-1)^(a+1))/(2*(a+1))
+end
+
+# compute r-th coefficient of product expansion of order p and order q normalized Legendre polynomials
+productseriescfs(p::T, q::T, r::T) where T = sqrt((2*p+1)*(2*q+1)/((2*(p+q-2*r)+1)))*(2*(p+q-2*r)+1)/(2*(p+q-r)+1)*exp(loggamma(r+one(T)/2)+loggamma(p-r+one(T)/2)+loggamma(q-r+one(T)/2)-loggamma(q-r+one(T))-loggamma(p-r+one(T))-loggamma(r+one(T))-loggamma(q+p-r+one(T)/2)+loggamma(q+p-r+one(T)))/π
+
+# # this generates the entire operator via normalized product Legendre decomposition
+# # this is very stable but scales rather poorly with high orders, so we only use it for testing
+# function productoperator(a::T, t::T, ℓ::Int) where T
+#     op::Matrix{T} = zeros(T,ℓ,ℓ)
+#     # first row where n arbitrary and m==0
+#     first = powerleg_backwardsfirstrow(a,t,2*ℓ+1)
+#     op[1,:] = first[1:ℓ]
+#     # generate remaining rows
+#     for p = 1:ℓ-1
+#         for q = p:ℓ-1
+#             productcfs = zeros(T,2*ℓ+1)
+#             for i = 0:min(p,q)
+#                 productcfs[1+q+p-2*i] = productseriescfs(p,q,i)
+#             end
+#             op[p+1,q+1] = dot(first,productcfs)
+#         end
+#     end
+#     # matrix is symmetric
+#     for m = 1:ℓ
+#         for n = m+1:ℓ
+#             op[n,m] = op[m,n]
+#         end
+#     end
+#     return op
+# end
+
+# This function returns the full ℓ×ℓ dot product operator, relying on several different methods for first row, second row, diagonal and remaining elements. We don't use this outside of the initial block.
+function gennormalizedpower(a::T, t::T, ℓ::Int) where T <: Real
+    # initialization
+    ℓ = ℓ+3
+    coeff = zeros(T,ℓ,ℓ)
+    # construct first row via stable backwards recurrence
+    first = powerleg_backwardsfirstrow(a,t,2*ℓ+1)
+    coeff[1,:] = first[1:ℓ]
+    # contruct second row via normalized product Legendre decomposition
+    @inbounds for q = 1:ℓ-1
+        productcfs = zeros(T,2*ℓ+1)
+        productcfs[q+2] = productseriescfs(1,q,0)
+        productcfs[q] = productseriescfs(1,q,1)
+        coeff[2,q+1] = dot(first,productcfs)
+    end
+    # contruct the diagonal via normalized product Legendre decomposition
+    @inbounds for q = 2:ℓ-1
+        productcfs = zeros(T,2*ℓ+1)
+        @inbounds for i = 0:q
+            productcfs[1+2*q-2*i] = productseriescfs(q,q,i)
+        end
+        coeff[q+1,q+1] = dot(first,productcfs)
+    end
+    #the remaining cases can be constructed iteratively by means of a T-shaped recurrence
+    @inbounds for m = 2:ℓ-2
+        # build remaining row elements
+        @inbounds for j = 1:m-2
+            coeff[j+2,m+1] = (t/((j+1)*(a+m+j+2)/((2*j+1)*(m+j+1)))*normconst_Pnadd1(j,t)*coeff[j+1,m+1]+((a+1)*m/(m*(m+1)-j*(j+1)))/((j+1)*(a+m+j+2)/((2*j+1)*(m+j+1)))*normconst_Pmnmix(m,j,t)*coeff[j+1,m]-(j*(a+m-j+1)/((2*j+1)*(m-j)))*1/((j+1)*(a+m+j+2)/((2*j+1)*(m+j+1)))*normconst_Pnsub1(j,t)*coeff[j,m+1])
+        end
+    end
+    #matrix is symmetric
+    @inbounds for m = 1:ℓ
+        @inbounds for n = m+1:ℓ
+            coeff[n,m] = coeff[m,n]
+        end
+    end
+    return coeff[1:ℓ-3,1:ℓ-3]
+end
+
+# the following version takes a previously computed block that has been resized and fills in the missing data guided by indices in inds
+function fillcoeffmatrix!(K::PowerLawMatrix, inds::AbstractUnitRange)
+    # the remaining cases can be constructed iteratively
+    a = K.a; t = K.t; T = eltype(promote(a,t));
+    ℓ = maximum(inds)
+    # fill in first row via stable backwards recurrence
+    first = powerleg_backwardsfirstrow(a,t,2*ℓ+1)
+    K.data[1,inds] = first[inds]
+    # fill in second row via normalized product Legendre decomposition
+    @inbounds for q = minimum(inds):ℓ-1
+        productcfs = zeros(T,2*ℓ+1)
+        productcfs[q+2] = productseriescfs(1,q,0)
+        productcfs[q] = productseriescfs(1,q,1)
+        K.data[2,q+1] = dot(first,productcfs)
+    end
+    # fill in the diagonal via normalized product Legendre decomposition
+    @inbounds for q = minimum(inds):ℓ-1
+        productcfs = zeros(T,2*ℓ+1)
+        @inbounds for i = 0:q
+            productcfs[1+2*q-2*i] = productseriescfs(q,q,i)
+        end
+        K.data[q+1,q+1] = dot(first,productcfs)
+    end
+    @inbounds for m in inds
+        m = m-2
+        # build remaining row elements
+        @inbounds for j = 1:m-1
+            K.data[j+2,m+2] = (t/((j+1)*(a+m+j+3)/((2*j+1)*(m+j+2)))*normconst_Pnadd1(j,t)*K.data[j+1,m+2]+((a+1)*(m+1)/((m+1)*(m+2)-j*(j+1)))/((j+1)*(a+m+j+3)/((2*j+1)*(m+j+2)))*normconst_Pmnmix(m+1,j,t)*K.data[j+1,m+1]-(j*(a+m-j+2)/((2*j+1)*(m+1-j)))*1/((j+1)*(a+m+j+3)/((2*j+1)*(m+j+2)))*normconst_Pnsub1(j,t)*K.data[j,m+2])
+        end
+    end
+    # matrix is symmetric
+    @inbounds for m in reverse(inds)
+        K.data[m,1:end] = K.data[1:end,m]
+    end
+end
+
+function dot(v::AbstractVector{T}, W::PowerLawMatrix, q::AbstractVector{T}) where T
+    vpad, qpad = paddeddata(v), paddeddata(q)
+    vl, ql = length(vpad), length(qpad)
+    return dot(vpad,W[1:vl,1:ql]*qpad)
+end