basis evaluation

src/dwt/evaluation.jl

@@ -132,18 +132,52 @@ function evaluate_in_dyadic_points{T}(side::DWT.Side, kind::DWT.Wvl, w::DWT.Disc
     end
 end
 
-function evaluate_periodic_in_dyadic_points(side::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d::Int=10;
-        points=false, options...) where {T}
-    f = cascade_algorithm_linear_combination(side, w, coeffs, d, perbound)
+evaluate_periodic_wavelet_basis_in_dyadic_points(s::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d::Int=10; options...) where {T} =
+    evaluate_periodic_scaling_basis_in_dyadic_points(s, w, DWT.idwt(coeffs, w, perbound), d; options...)
+
+function evaluate_periodic_scaling_basis_in_dyadic_points(s::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d::Int=10; options...) where {T}
+    f = zeros(T, evaluate_periodic_in_dyadic_points_scratch_length(s, scaling, w, Int(log2(length(coeffs))), 0, d))
+    scratch = zeros(T, evaluate_periodic_in_dyadic_points_scratch2_length(s, scaling, w, Int(log2(length(coeffs))), 0, d))
+    f_scaled = similar(f)
+    evaluate_periodic_wavelet_basis_in_dyadic_points!(zeros(T, 1<<d), s, w, coeffs, d, f, f_scaled, scratch; options...)
+end
+
+evaluate_periodic_wavelet_basis_in_dyadic_points!(y::Vector{T}, s::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d, f, f_scaled, scratch; options...) where T=
+    evaluate_periodic_scaling_basis_in_dyadic_points!(y, s, w, coeffs, d, f, f_scaled, scratch; options...)
+
+function evaluate_periodic_scaling_basis_in_dyadic_points!(y::Vector{T}, s::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d::Int, f, f_scaled, scratch;
+    points = false,
+    options...) where {T}
+    j = Int(log2(length(coeffs)))
+    DWT.evaluate_in_dyadic_points!(f, s, scaling, w, j, 0, d, scratch)
+    f_scaled = similar(f)
+    for k in 0:(1<<j)-1
+        f_scaled .= f .* coeffs[k+1]
+        DWT._periodize_add!(y, f_scaled, -ceil(Int,(1<<d)*DWT.support(s, scaling, w, j, k)[1])+1)
+    end
     if points
-        f, periodic_dyadic_points(d, T)
+        y, DWT.periodic_dyadic_points(d, T)
     else
-        f
+        y
     end
 end
 
-inv_evaluate_periodic_in_dyadic_points(side::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d::Int=10) where {T} =
-    inv_cascade_algorithm_linear_combination(side, w, coeffs, d, perbound)
+
+
+
+# Not correctly implented
+# function evaluate_periodic_in_dyadic_points(side::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d::Int=10;
+#         points=false, options...) where {T}
+#     f = cascade_algorithm_linear_combination(side, w, coeffs, d, perbound)
+#     if points
+#         f, periodic_dyadic_points(d, T)
+#     else
+#         f
+#     end
+# end
+#
+# inv_evaluate_periodic_in_dyadic_points(side::DWT.Side, w::DWT.DiscreteWavelet{T}, coeffs, d::Int=10) where {T} =
+#     inv_cascade_algorithm_linear_combination(side, w, coeffs, d, perbound)
 
 
 

src/dwt/scaling_coefficients.jl

@@ -95,42 +95,48 @@ _scalingcoefficient_filter(f::CompactSequence) =
     # scaling_coefficients_to_dyadic_grid(scaling_coefficients, Primal, w, bmd, d; options...)
 
 "Transforms scaling coeffients back to function evaluations on the dyadic grid."
-function scaling_coefficients_to_dyadic_grid{T}(scaling_coefficients::AbstractArray{T,1}, s::Side, w::DWT.DiscreteWavelet{T}, bnd::WaveletBoundary, d=ndyadicscales(scaling_coefficients); grid=false, options...)
+function scaling_coefficients_to_dyadic_grid{T}(scaling_coefficients::AbstractArray{T,1}, s::Side, w::DWT.DiscreteWavelet{T}, bnd::DWT.PeriodicBoundary, d=ndyadicscales(scaling_coefficients); grid=false, options...)
     function_evals = zeros(T,DWT.scaling_coefficients_to_dyadic_grid_length(s,w,d))
-    scratch = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch_length(s,w,d))
-    scratch2 = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(s,w,d))
-
-    scaling_coefficients_to_dyadic_grid!(function_evals, scaling_coefficients, w, bnd, scratch, scratch2; options...)
+    j = Int(log2(length(scaling_coefficients)))
+    f = zeros(T, evaluate_periodic_in_dyadic_points_scratch_length(s, scaling, w, j, 0, d))
+    scratch = zeros(T, evaluate_periodic_in_dyadic_points_scratch2_length(s, scaling, w, j, 0, d))
+    f_scaled = similar(f)
+    # scratch = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch_length(s,w,d))
+    # scratch2 = zeros(DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(s,w,d))
+
+    scaling_coefficients_to_dyadic_grid!(function_evals, scaling_coefficients, w, bnd, f, f_scaled, scratch; options...)
     grid ?
         (return function_evals, linspace(T(0), T(1), length(function_evals)+1)[1:end-1]) :
         (return function_evals)
 end
 
 # Scaling coefficients to function evaluations on dyadic grid (assumes periodic extension)
 "In place method of scaling_coefficients_to_dyadic_grid."
-scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray{T,1}, scaling_coefficients::AbstractArray{T,1}, w::DWT.DiscreteWavelet{T}, bnd::DWT.PeriodicBoundary, scratch, scratch2; options...) =
-    scaling_coefficients_to_dyadic_grid!(function_evals, scaling_coefficients, Primal, w, bnd, scratch, scratch2; options...)
+scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray{T,1}, scaling_coefficients::AbstractArray{T,1}, w::DWT.DiscreteWavelet{T}, bnd::DWT.PeriodicBoundary, f, f_scaled, scratch; options...) =
+    evaluate_periodic_scaling_basis_in_dyadic_points!(function_evals, Primal, w, scaling_coefficients, Int(log2(length(function_evals))), f, f_scaled, scratch; options...)
+
+    # scaling_coefficients_to_dyadic_grid!(function_evals, scaling_coefficients, Primal, w, bnd, scratch, scratch2; options...)
 
 
 "In place method of scaling_coefficients_to_dyadic_grid."
-function scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray{T,1}, scaling_coefficients::AbstractArray{T,1}, s::Side, w::DWT.DiscreteWavelet{T}, ::DWT.PeriodicBoundary, scratch, scratch2; grid=false, options...)
-    d = Int(log2(length(function_evals)))
-    @assert length(function_evals) == length(scratch)
-    @assert DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(s, w, d) == length(scratch2)
-    @assert isdyadic(scaling_coefficients)
-    j = ndyadicscales(scaling_coefficients)
-    function_evals[:] = T(0)
-    for (c_i,c) in enumerate(scaling_coefficients)
-        k = c_i - 1
-        DWT.evaluate_periodic_in_dyadic_points!(scratch, s, DWT.Scl(), w, j, k, d, scratch2, nothing)
-        scale!(scratch, c)
-        for i in 1:length(function_evals)
-            function_evals[i] += scratch[i]
-        end
-        # function_evals[:] += c*DWT.evaluate_periodic_in_dyadic_points(Prl(), Scl(), w, j, k, d)
-    end
-    nothing
-end
+# function scaling_coefficients_to_dyadic_grid!{T}(function_evals::AbstractArray{T,1}, scaling_coefficients::AbstractArray{T,1}, s::Side, w::DWT.DiscreteWavelet{T}, ::DWT.PeriodicBoundary, scratch, scratch2; grid=false, options...)
+#     d = Int(log2(length(function_evals)))
+#     @assert length(function_evals) == length(scratch)
+#     @assert DWT.scaling_coefficients_to_dyadic_grid_scratch2_length(s, w, d) == length(scratch2)
+#     @assert isdyadic(scaling_coefficients)
+#     j = ndyadicscales(scaling_coefficients)
+#     function_evals[:] = T(0)
+#     for (c_i,c) in enumerate(scaling_coefficients)
+#         k = c_i - 1
+#         DWT.evaluate_periodic_in_dyadic_points!(scratch, s, DWT.Scl(), w, j, k, d, scratch2, nothing)
+#         scale!(scratch, c)
+#         for i in 1:length(function_evals)
+#             function_evals[i] += scratch[i]
+#         end
+#         # function_evals[:] += c*DWT.evaluate_periodic_in_dyadic_points(Prl(), Scl(), w, j, k, d)
+#     end
+#     nothing
+# end
 
 scaling_coefficients_to_dyadic_grid_length(s::Side, w::DWT.DiscreteWavelet, d::Int) = 1<<d
 

src/dwt/util/cascade.jl

@@ -15,13 +15,17 @@
     function with period `1`.
 """
 function cascade_algorithm_linear_combination(side::Side, w::DiscreteWavelet{T}, coeffs, d::Int, bnd::WaveletBoundary) where {T}
+    # TODO implement for different lengths.
+    @assert d ≈ log2(length(coeffs))
     input = zeros(T,1<<d)
     copy!(input, coeffs)
     T(2)^(d//2)*idwt(input, side, w, bnd)
 end
 
 
 function inv_cascade_algorithm_linear_combination(side::Side, w::DiscreteWavelet{T}, coeffs, d::Int, bnd::WaveletBoundary) where {T}
+    # TODO implement for different lengths.
+    @assert d ≈ log2(length(coeffs))
     input = zeros(T,1<<d)
     copy!(input, coeffs)
     T(2)^(d//2)*dwt(input, side, w, bnd)

src/dwt/util/periodize.jl

@@ -44,3 +44,15 @@ function _periodize!{T}(dest::AbstractArray{T}, src::AbstractArray{T}, istart::I
         dest[i] = t
     end
 end
+
+function _periodize_add!{T}(dest::AbstractArray{T}, src::AbstractArray{T}, istart::Int, step::Int=1)
+    L = length(dest)
+    srclength = length(src)
+    for i in 1:L
+        t = T(0)
+        for m in mod(istart-1+step*(i-1),step*L)+1:step*L:srclength
+            t += src[m]
+        end
+        dest[i] += t
+    end
+end

src/dwt/util/scratchspace.jl

@@ -98,8 +98,23 @@ ScratchSpace(s::Vector{Vector{T}}) where {T} = ScratchSpace{T}(s, map(length, s)
 
 "Retrieve vectors of given lengths of the scratch space container."
 get_scratch_space(SS::ScratchSpace, lengths) =
-    [SS.scratch[index]  for index in [find(SS.lengths.==l)[1] for l in lengths]]
-
+    [SS.scratch[index]  for index in _get_indices(SS.lengths, lengths)]
+
+function _get_indices(lengths::Vector{Int}, input)
+    L = length(lengths)
+    LL = length(input)
+    r = Array{Int}(LL)
+    @inbounds for i in 1:LL
+        s = input[i]
+        for j in 1:L
+            if lengths[j] == s
+                r[i] = j
+                break
+            end
+        end
+    end
+    r
+end
 
 "Scratch space for a wavelet evaluation with l levels"
 EvalScratchSpace(side::Side, w::DiscreteWavelet, l::Int, d::Int) =

src/sequences/compactsequences.jl

@@ -43,13 +43,11 @@ function getindex(s::CompactSequence, c::Range)
 end
 
 function Base.getindex(s::CompactSequence, c::UnitRange)
-  e = zeros(eltype(s), length(c))
-  i1 = max(s.offset, c[1])
-  i2 = min(s.offset + s.n-1, c[end])
-  e1 = max(1, 1-c[1]+s.offset)
-  e2 = min(length(e),i2-i1+1-c[1]+s.offset)
-  e[e1:e2] = s.a[i1-s.offset+1:i2-s.offset+1]
-  e
+    e = zeros(eltype(s), length(c))
+    i1 = max(s.offset, c[1])
+    i2 = min(s.offset + s.n-1, c[end])
+    e[(1:i2-i1+1)+max(0,s.offset-c[1])] = s.a[i1-s.offset+1:i2-s.offset+1]
+    e
 end
 
 firstindex(s::CompactSequence) = imapindex(s, 1)
@@ -58,6 +56,27 @@ lastindex(s::CompactSequence) = imapindex(s, sublength(s))
 
 hascompactsupport{T}(::Type{CompactSequence{T}}) = True
 
+Base.isapprox(s1::CompactSequence, s2::CompactSequence) = s1.a≈s2.a
+
+function shifted_conv(c1::CompactSequence{ELT}, c2::CompactSequence{ELT}, shift::Int) where {ELT}
+    l1 = sublength(c1)
+    l2 = sublength(c2)
+    o1 = c1.offset
+    o2 = c2.offset
+    offset = o1+shift*o2
+    L = (l2-1)+shift*(l1-1)+1
+    a = zeros(ELT, L)
+    for ai in 0:L-1
+        t = ELT(0)
+        for k in max(0,floor(Int,(firstindex(c1)-l2+1)//shift)):max(0,floor(Int,(lastindex(c1)+l2)//shift))
+#         for k in firstindex(c1):lastindex(c1)
+            t += c1[k]*c2[ai-shift*k]
+        end
+        a[ai+1] = t
+    end
+    CompactSequence(a, offset)
+end
+
 """
 From a given filter h_i, compute a new filter satisfying the alternating flip relation,
 centered around the given pivot:
@@ -166,7 +185,6 @@ shift{L,OFS,T}(s::FixedSequence{L,OFS,T}, k::Int) = FixedSequence(s.a, Val{OFS+k
 
 hascompactsupport{L,OFS,T}(::Type{FixedSequence{L,OFS,T}}) = True
 
-
 for op in (:ctranspose, :evenpart, :oddpart, :alternating_flip, :reverse, :conj, :alternating)
     @eval $op{L,OFS,T}(s::FixedSequence{L,OFS,T}) = FixedSequence($op(CompactSequence{T}(s)))
 end

src/sequences/extensionsequences.jl

@@ -24,7 +24,7 @@ abstract type ExtensionSequence{A} <: Sequence end
 # imapindex(s::SomeSubType, i) -> the inverse map
 
 eltype{A}(::Type{ExtensionSequence{A}}) = eltype(A)
-eltype{E <: ExtensionSequence}(::Type{E}) = eltype(super(E))
+eltype{E <: ExtensionSequence}(::Type{E}) = eltype(supertype(E))
 
 "The subvector of the extension sequence."
 subvector(s::ExtensionSequence) = s.a
@@ -108,7 +108,7 @@ end
 ZeroPadding{A}(a::A) = ZeroPadding{A}(a)
 
 # We override getindex to return zero outside our embedded vector.
-getindex(s::ZeroPadding, k::Int) = (k < 0) || (k >= s.n) ? zero(eltype(s)) : getindex(s.a, k+1)
+Base.getindex(s::ZeroPadding, k::Int) = (k < 0) || (k >= s.n) ? eltype(s)(0) : getindex(s.a, k+1)
 
 firstindex(s::ZeroPadding) = 0
 

test/suite_evaluation.jl

@@ -372,7 +372,7 @@ function linear_combo_test()
         c = zeros(1<<l)
         for index in wavelet_indices(l)
             c[value(index)] = 1
-            t1 = evaluate_periodic_in_dyadic_points(Primal, cdf11, c, d)
+            t1 = DWT.evaluate_periodic_wavelet_basis_in_dyadic_points(Primal, cdf11, c, d)
             c[value(index)] = 0
             t2 = evaluate_periodic_in_dyadic_points(Primal, kind(index), cdf11, level(index), offset(index), d)
             @test t1≈t2