Merge ab541fe into a8edcc1

.travis.yml

@@ -2,7 +2,7 @@ language: julia
 os:
     - linux
 julia:
-    - nightly
+    - release
 notifications:
     - email: false
 script:

README.md

@@ -8,6 +8,6 @@ Scattering coefficients are powerful features for signal-based machine learning
 while discarding non-informative variability, such as rigid motions, frequency transpositions, and elastic deformations.
 
 Its design is inspired from the [scattering.m](https://github.com/lostanlen/scattering.m) toolbox, while taking advantage of
-Julia's multiple dispatch and in-place operations. 
+Julia's multiple dispatch and in-place operations.
 
-The package is under ongoing development and is planned to be released in the fall 2015.
+The package is under ongoing development and is planned to be released in 2016.

REQUIRE

@@ -3,3 +3,5 @@ Clp 0.0.10
 DataStructures 0.3.11
 JuMP 0.10.1
 MathProgBase 0.3.17
+Mocha 0.0.9
+Wavelets 0.4.1

src/WaveletScattering.jl

@@ -1,16 +1,37 @@
 module WaveletScattering
 
-using Clp
-using DataStructures
+# using Clp
 using JuMP
 using MathProgBase
+using Mocha
 
-include("spec.jl")
+include("domain.jl")
+include("group.jl")
+include("waveletclass.jl")
+include("path.jl")
 include("meta.jl")
+include("spec.jl")
+include("spec1d.jl")
 include("filter.jl")
-include("fourierfilter.jl")
+include("fourier1dfilter.jl")
+include("weighting.jl")
+include("behavior1d.jl")
 include("bank.jl")
+include("bank1d.jl")
 include("morlet1d.jl")
-include("path.jl")
+include("node.jl")
+include("blob.jl")
+include("pointwise.jl")
+include("modulus.jl")
+include("symbols.jl")
+include("inputlayer.jl")
+include("pointwiselayer.jl")
+include("fourierlayer.jl")
+include("waveletlayer.jl")
+include("layerstate.jl")
+include("inputlayerstate.jl")
+include("pointwiselayerstate.jl")
+include("fourierlayerstate.jl")
+include("waveletlayerstate.jl")
 
 end

src/bank.jl

@@ -1,107 +1,8 @@
-"""A `Behavior` object contains the mutable information in a filter bank. The
-values of these fields may be changed *after* the construction of the filter
-bank, without having to recompute the underlying architecture. Fields:
+abstract AbstractBank{
+    T<:Number,
+    D<:AbstractDomain,
+    G<:AbstractPointGroup,
+    W<:RedundantWaveletClass}
 
-* `γ_range`: range of wavelet indices that are actually used in the
-convolutions. Default is all of them (see outer constructor).
-
-* `is_ϕ_applied`: true if and only if the lowpass filter `ϕ` (also known as
-scaling function) in addition to the bandpass filters `ψ`'s.
-
-* `log2_oversampling`: base-2 logarithm of the oversampling factor with respect
-to the critical sampling rate. Must be positive. Default is 0, i.e. no
-oversampling."""
-type Behavior
-    γ_range::UnitRange
-    is_ϕ_applied::Bool
-    log2_oversampling::Int
-    max_log2_stride::Int
-end
-
-"""An `AbstractBank` is a wavelet filter bank. Filter banks are of two kinds:
-* `AbstractNonOrientedBank`: no orientation variable `θ`, only a scale
-variable `γ`. One-dimensional filter banks with real inputs are non-oriented.
-* `AbstractOrientedBank`: wavelets have an orientation variable `θ` and a scale
-variable `γ`. Two-dimensional filter banks are oriented. One-dimensional filter
-banks with complex outputs are also oriented, in the sense that the sign
-of the center frequency is regarded as an orientation variable with two
-possible values."""
-abstract AbstractBank{T<:Number}
-abstract AbstractNonOrientedBank{T<:Number} <: AbstractBank{T}
-abstract AbstractOrientedBank{T<:Number} <: AbstractBank{T}
-
-"""A `FourierNonOriented1DBank` is a one-dimensional, non-oriented filter bank
-defined in the Fourier domain. It is not oriented in the sense that only the
-positive frequencies are guaranteed to be covered. Indeed, assuming that the
-input signal will be real — i.e. its Fourier will be symmetric — it can be
-recovered from the ""positive half"" of its Fourier spectrum. In summary, it is
-advisable to use this type of filter bank when handling real 1d data of
-moderate to large length."""
-immutable FourierNonOriented1DBank{T<:Number} <: AbstractNonOrientedBank{T}
-    ψs::Vector{AbstractFourier1DFilter{T}}
-    ϕ::Symmetric1DFilter{T}
-    behavior::Behavior
-    metas::Vector{NonOrientedMeta}
-    spec::Abstract1DSpec{T}
-    function call{T<:Number}(::Type{FourierNonOriented1DBank{T}},
-                             spec::Abstract1DSpec ;
-                             γ_range = nothing, is_ϕ_applied = false,
-                             log2_oversampling = 0, max_log2_stride = 0)
-        T == spec.signaltype || error("""Type parameter of
-        FourierNonOriented1DBankmust must be equal to spec.signaltype""")
-        γs, χs, js = gammas(spec), chromas(spec), octaves(spec)
-        ξs, qs = centerfrequencies(spec), qualityfactors(spec)
-        scs, bws = scales(spec), bandwidths(spec)
-        @inbounds metas = [ NonOrientedMeta(
-            γs[1+γ], χs[1+γ], bws[1+γ], ξs[1+γ], js[1+γ], qs[1+γ], scs[1+γ])
-            for γ in γs ]
-        ψs = pmap(fourierwavelet, metas, fill(spec, length(metas)))
-        ψs = convert(Array{AbstractFourier1DFilter{T}}, ψs)
-        ϕ = scalingfunction(spec)
-        renormalize!(ψs, ϕ, metas, spec)
-        (γ_range == nothing) && (γ_range = 0:(length(γs)-1))
-        behavior =
-            Behavior(γ_range, is_ϕ_applied, log2_oversampling, max_log2_stride)
-        new{T}(ψs, ϕ, behavior, metas, spec)
-    end
-end
-FourierNonOriented1DBank(spec::Abstract1DSpec ; args...) =
-    FourierNonOriented1DBank{spec.signaltype}(spec, args...)
-
-"""A `FourierOriented1DBank` is a one-dimensional, oriented filter bank defined
-in the Fourier domain. It is "oriented" insofar as its filters have negative
-center frequencies as well as positive center frequencies, as represented by
-the orientation parameter `θ`. It is advisable to use this type of filter bank
-when handling complex 1d data of moderate to large length."""
-immutable FourierOriented1DBank{T<:Number} <: AbstractOrientedBank{T}
-    ψs::Matrix{AbstractFourier1DFilter{T}}
-    ϕ::Symmetric1DFilter{T}
-    behavior::Behavior
-    metas::Matrix{OrientedMeta}
-    spec::Abstract1DSpec{T}
-    function call{T<:Number}(::Type{FourierOriented1DBank{T}},
-                             spec::Abstract1DSpec ;
-                             γ_range = nothing, is_ϕ_applied = false,
-                             log2_oversampling = 0, max_log2_stride = 0)
-        T == spec.signaltype || error("""Type parameter of
-        FourierNonOriented1DBankmust be equal to spec.signaltype""")
-        γs, χs, js = gammas(spec), chromas(spec), octaves(spec)
-        ξs, qs = centerfrequencies(spec), qualityfactors(spec)
-        scs, bws = scales(spec), bandwidths(spec)
-        θs = 0:1
-        @inbounds metas = [ OrientedMeta(
-            γs[γ], θs[θ], χs[γ], bws[γ], ξs[γ], js[γ], qs[γ], scs[γ])
-            for γ in eachindex(γs), θ in 1:2 ]
-        ψs = pmap(fourierwavelet, metas[:, 1], fill(spec, length(metas)))
-        ψs = convert(Array{AbstractFourier1DFilter{T}}, ψs)
-        ψs = hcat(ψs, map(spin, ψs))
-        ϕ = scalingfunction(spec)
-        renormalize!(ψs, ϕ, metas, spec)
-        (γ_range == nothing) && (γ_range = 0:(length(γs)-1))
-        behavior =
-            Behavior(γ_range, is_ϕ_applied, log2_oversampling, max_log2_stride)
-        new{T}(ψs, ϕ, behavior, metas, spec)
-    end
+immutable NullBank <: AbstractBank
 end
-FourierOriented1DBank(spec::Abstract1DSpec ; args...) =
-    FourierOriented1DBank{spec.signaltype}(spec, args...)

src/bank1d.jl

@@ -0,0 +1,74 @@
+"""A `Bank1D` is a one-dimensional wavelet filter bank, parametrized by
+* `T`: numeric type of input, e.g. Float32, Float64.
+* `D`: transform domain. Either `FourierDomain{1}` or `SpatialDomain{1}`.
+* `G`: point group. Either `TrivialGroup` or `ReflectionGroup`.
+* `W`: wavelet class, e.g. `Morlet` or `Gammatone`.
+Its fields are
+* `ϕ`: low-pass filter, also called scaling function.
+* `ψs`: 3d array of wavelets, indexed by spin, chroma, and octave.
+* `behavior`: mutable behavior, e.g. oversampling.
+* `spec`: immutable specifications, e.g. number of filters per octave.
+To create a `Bank1D`
+1. define a `Spec1D`,
+2. define a `PathKey` on which to apply the wavelet filterbank,
+3. if needed, provide behavior characteristics as keyword arguments.
+Example:
+spec = Spec1D(nFilters_per_octave = 12, nOctaves = 10)
+pathkey = PathKey(:time)
+bank = Bank1D(spec, pathkey, j_range = 2:9)
+"""
+immutable Bank1D{
+        T<:Number,
+        D<:LineDomains,
+        G<:LineGroups,
+        W<:RedundantWaveletClass} <: AbstractBank{T,D,G,W}
+    ϕ::AbstractFilter{T,D}
+    ψs::Array{AbstractFilter{T,D},3}
+    behavior::Behavior1D{T}
+    spec::Spec1D{T,D,G,W}
+    function call{T,D,G,W}(::Type{Bank1D},
+            spec::Spec1D{T,D,G,W},
+            pathkey::PathKey = PathKey(:time) ;
+            is_ϕ_applied::Bool = false,
+            j_range::UnitRange{Int} = 0:(spec.nOctaves-1),
+            log2_oversampling::Int = 0,
+            max_log2_stride::Int = spec.nOctaves-1,
+            weighting::AbstractWeighting = EqualWeighting())
+        (nΘs, nΧs, nJs) = size(spec.ψmetas)
+        ψs = Array(AbstractFilter{T,D}, (nΘs, nΧs, nJs))
+        ψs[1, :, :] =
+            pmap(AbstractFilter, spec.ψmetas[1, :, :], fill(spec, nΧs * nJs))
+        (nΘs > 1) && spin!(ψs)
+        ϕ = AbstractFilter(spec.ϕmeta, spec)
+        renormalize!(ϕ, ψs, spec)
+        behavior = Behavior1D(ϕ, ψs, spec, is_ϕ_applied, j_range,
+            log2_oversampling, max_log2_stride, pathkey, weighting)
+        new{T,D,G,W}(ϕ, ψs, behavior, spec)
+    end
+end
+
+Base.ndims(::Bank1D) = 1
+
+function call{T<:Real,DIM}(
+        bank::Bank1D{T,FourierDomain{1},TrivialGroup},
+        x::AbstractArray{T,DIM};
+        flags = FFTW.ESTIMATE,
+        timelimit = Inf,
+        verbose = 0)
+    syms = appendsymbols(fill(symbol(bank.behavior.pathkey), 1), DIM)
+    inputnode = Node(x, ntuple(k -> kthrange(syms, x, k), DIM))
+    fouriernode = AbstractFourierNode(inputnode, [1], flags, timelimit)
+    for j in bank.behavior.j_range
+        ψ_log2_sampling = bank.behavior.ψ_log2_samplings[1+j]
+        downsampled_length = size(x, 1) << (-ψ_log2_sampling)
+        octave_size = (downsampled_length,
+            size(x)[2:end]..., bank.spec.nFilters_per_octave)
+        octave_ft = zeros(Complex{T}, octave_size)
+        inds = [fill(Colon(), DIM) ; 0]
+        for χ in 0:(bank.spec.nFilters_per_octave-1)
+            ψ = bank.ψs[1, 1+χ, 1+j]
+            inds[end] = 1 + χ
+            transform!(sub(octave_ft, inds...), ψ, fouriernode, 1)
+        end
+    end
+end

src/behavior1d.jl

@@ -0,0 +1,77 @@
+"""A `Behavior` object contains the mutable information in a filter bank. The
+values of these fields may be changed *after* the construction of the filter
+bank, without having to recompute the underlying architecture. Fields:
+
+* `is_ϕ_applied`: true if and only if the lowpass filter `ϕ` (also known as
+scaling function) in addition to the bandpass filters `ψ`'s.
+
+* `j_range`: range of wavelet octaves that are actually used in the
+convolutions. Default is all of them.
+
+* `log2_oversampling`: base-2 logarithm of the oversampling factor with respect
+to the critical sampling rate. Must be positive. Default is 0, i.e. no
+oversampling.
+
+* `max_log2_stride`: base-2 logarithm of the maximum distance between
+neighboring coefficients after subsampling (also known as *hop size* or
+*stride*). Must be positive. Default is the number of octaves minus one, which
+imposes no oversampling per se.
+
+* `pathkey`: key of the variable over which the filter bank is applied.
+
+* `weighting`: scalar weighting of each wavelet frequency band. Default is
+`EqualWeighting()`, i.e. all weights are one. Can be set to
+`LoudnessWeighting(samplerate)` to model the relative loudness perceived by the
+human ear, as defined by the international standard 61672:2003.
+
+* `weights`: vector of floating-point weights corresponding to the
+frequencies of `ψ`'s. The weight assigned to `ϕ` is the same as the weight
+of the `ψ` with the lowest center frequency."""
+type Behavior1D{T<:Real}
+    ϕ_log2_sampling::Int
+    ψ_log2_samplings::Vector{Int}
+    is_ϕ_applied::Bool
+    j_range::UnitRange
+    log2_oversampling::Int
+    max_log2_stride::Int
+    pathkey::PathKey
+    weighting::AbstractWeighting
+    weights::Array{T, 3}
+end
+
+"""Given a lowpass filter `ϕ`, an array of wavelets `ψs`, and the corresponding
+filter bank specification `spec`, the `Behavior` outer constructor computes
+the critical sampling rates of all `ψs` and `ϕ`. Then, with the inputs
+`log2_oversampling` and `max_log2_stride`, these sampling rates are bounded:
+* from below, by the inverse of the maximum stride.
+* from above, by 1, to avoid unnecessary upsampling with respect to the original
+  sample rate.
+NB: if `log2_oversampling` and `max_log2_stride` are left as defaults, all
+sample rates are critical."""
+function Behavior1D{
+    T<:Number,
+    D<:AbstractDomain,
+    G<:AbstractPointGroup,
+    W<:RedundantWaveletClass}(
+        ϕ::AbstractFilter{T,D},
+        ψs::AbstractArray{AbstractFilter{T,D},3},
+        spec::AbstractSpec{T,D,G,W},
+        is_ϕ_applied::Bool,
+        j_range::UnitRange{Int},
+        log2_oversampling::Int,
+        max_log2_stride::Int,
+        pathkey::PathKey,
+        weighting::AbstractWeighting)
+    ϕ_critical_log2_sampling = critical_log2_sampling(ϕ, spec)
+    ϕ_log2_sampling =
+        clamp(ϕ_critical_log2_sampling + log2_oversampling, -max_log2_stride, 0)
+    ψ_critical_log2_samplings =
+        Int[ critical_log2_sampling(ψ, spec) for ψ in ψs[1, end, :] ]
+    ψ_log2_samplings = clamp(ψ_critical_log2_samplings + log2_oversampling,
+        -max_log2_stride, 0)
+    max_log2_stride = - min(ϕ_log2_sampling, minimum(ψ_log2_samplings))
+    ξs = map(get_centerfrequency, spec.ψmetas)
+    weights = map(T, weight_frequencies(weighting, ξs))
+    Behavior1D(ϕ_log2_sampling, ψ_log2_samplings, is_ϕ_applied, j_range,
+        log2_oversampling, max_log2_stride, pathkey, weighting, weights)
+end

src/blob.jl

@@ -0,0 +1,31 @@
+# ScatteredBlob
+immutable ScatteredBlob{NODE<:AbstractNode} <: Mocha.Blob
+    nodes::Dict{Path,NODE}
+end
+
+function Base.show(io::IO, blob::ScatteredBlob)
+    n_nodes = length(blob.nodes)
+    plural = n_nodes > 1 ? "s" : ""
+    print(io, "ScatteredBlob(", n_nodes, " node", plural, ")")
+end
+
+function forward!(
+        backend::Mocha.CPUBackend,
+        blob_out::ScatteredBlob,
+        bank::Bank1D,
+        blob_in::ScatteredBlob)
+    pathdepth(blob_in, bank.behavior.pathkey)
+
+    map(node -> pathdepth(bank.behavior.pathkey, keys(blob_in.nodes)))
+    γkey = cons(Literal(:γ, 1), bank.behavior.pathkey)
+    for j in bank.behavior.j_range
+        for χ in 0:(bank.spec.nFilters_per_octave-1)
+            ψ = bank.ψs[1 + θ, 1 + χ, 1 + j]
+            for (path_in, node_in) in input.nodes
+                path_out = copy(path_in)
+                path_out[γkey] = γ
+                transform!(blob[path_out], blob_in[path_in], ψ)
+            end
+        end
+    end
+end

src/domain.jl

@@ -0,0 +1,32 @@
+"""There are two subtypes of `AbstractDomain`:\n
+* `FourierDomain{K}`: filtering is implemented as a product in the Fourier
+domain. `K` is the dimension of the Fourier transform.\n
+* `SpatialDomain{K}`: filtering is implemented as a convolution in the spatial
+domain, or as a wavelet lifting scheme when it is possible. `K` is the number of
+subscripts over which the convolution operates.
+"""
+abstract AbstractDomain
+
+immutable FourierDomain{K} <: AbstractDomain
+    dim::Val{K}
+    function FourierDomain{K}(dim::Val{K})
+        @assert (K > 0)
+        @assert isa(K, Int)
+        new(dim)
+    end
+end
+FourierDomain(dim::Int) = FourierDomain{dim}(Val{dim}())
+
+
+immutable SpatialDomain{K} <: AbstractDomain
+    dim::Val{K}
+    function SpatialDomain{K}(dim::Val{K})
+        @assert (K > 0)
+        @assert isa(K, Int)
+        new(dim)
+    end
+end
+SpatialDomain(dim::Int) = SpatialDomain{dim}(Val{dim}())
+
+typealias LineDomains Union{FourierDomain{1},SpatialDomain{1}}
+typealias PlaneDomains Union{FourierDomain{2},SpatialDomain{2}}