text:
@testset "Accuracy" begin
    # test against data made by data/make_filter_data.m in Octave
    wname = Any["Daubechies","Coiflet","Haar","Symmlet","Battle","Vaidyanathan","Beylkin"];
    wtype = Any[WT.Daubechies, WT.Coiflet, WT.Haar, WT.Symlet, WT.Battle, WT.Vaidyanathan, WT.Beylkin];
    wnum = Any[collect(4:2:20),collect(2:5),[0],collect(4:10),[1,3,5],[0],[0]];
    wvm = Any[collect(2:1:10),collect(4:2:10),[0],collect(4:10),[2,4,6],[0],[0]];

    #print("transforms: reading 1d and 2d data ...\n")
    name = "filter"
    data = vec(readdlm(joinpath(dirname(@__FILE__), "data", "filter1d_data.txt"),'\t'))
    data2 = readdlm(joinpath(dirname(@__FILE__), "data", "filter2d_data.txt"),'\t')

    #print("transforms: testing 1d and 2d data ...\n")
    stderr1 = 1e-9*sqrt(length(data))
    stderr2 = 1e-9*sqrt(length(data2))

    for i in eachindex(wname)
        for num in eachindex(wnum[i])
            # expected transforms
            ye = vec(readdlm(joinpath(dirname(@__FILE__), "data", string(name,"1d_",wname[i],wnum[i][num],".txt")),'\t'))
            ye2 = readdlm(joinpath(dirname(@__FILE__), "data", string(name,"2d_",wname[i],wnum[i][num],".txt")),'\t')

            wtc = wtype[i]
            if wnum[i][num] != 0
                class = wtc{wvm[i][num]}()
            else
                class = wtc()
            end
            wt = wavelet(class, WT.Filter)

            # transform data
            y = dwt(data, wt)
            y2 = dwt(data2, wt)

            @test (@vecnorm_eq_eps y ye stderr1)
            @test (@vecnorm_eq_eps y2 ye2 stderr2)
            # a few bad cases have quite large norms
            if wname[i] != "Battle" && (wname[i] != "Coiflet" && wvm[i][num]==10)
                @test abs(norm(data)-norm(y)) < 1e-9
                @test abs(norm(data2)-norm(y2)) < 1e-9
                @test (@vecnorm_eq_eps idwt(y,wt) data stderr1*100)
                @test (@vecnorm_eq_eps idwt(y2,wt) data2 stderr2*100)
            end
        end
    end
end

@testset "Accuracy non-square" begin
    data2 = readdlm(joinpath(dirname(@__FILE__), "data", "filter2d_nonsquare_data.txt"))
    y2  = dwt(data2, wavelet(WT.haar), 1)
    ye2 = readdlm(joinpath(dirname(@__FILE__), "data", "filter2d_nonsquare_Haar0.txt"))
    stderr2 = 1e-9*sqrt(length(ye2))
    @test (@vecnorm_eq_eps y2 ye2 stderr2)
end

@testset "Lifting vs filter" begin
    n = 32
    stderr1 = 1e-10*sqrt(n)
    stderr2 = 1e-10*sqrt(n*n)
    # 1-D, 2-D, 3-D lifting vs filtering / inverse vs original tests
    @testset "wclass $wclass" for wclass in (WT.db1, WT.db2)
        wf = wavelet(wclass, WT.Filter)
        wls = wavelet(wclass, WT.Lifting)
        x = randn(n)
        x2 = copy(x)

        for L in (ndyadicscales(n),0,1,2)

            yf = dwt(x, wf, L)
            yls = dwt(x, wls, L)

            # filter vs lifting
            @test (@vecnorm_eq_eps yf yls stderr1)

            ytf = idwt(yf, wf, L)
            ytls = idwt(yls, wls, L)

            # inverse vs original
            @test (@vecnorm_eq_eps ytf x stderr1)
            @test (@vecnorm_eq_eps ytls x stderr1)
        end

        x = randn(n,n)
        x2 = copy(x)
        for L in (ndyadicscales(n),0,1,2)

            yf = dwt(x, wf, L)
            yls = dwt(x, wls, L)

            # filter vs lifting
            @test (@vecnorm_eq_eps yf yls stderr2)

            ytf = idwt(yf, wf, L)
            ytls = idwt(yls, wls, L)

            # inverse vs original
            @test (@vecnorm_eq_eps ytf x stderr2)
            @test (@vecnorm_eq_eps ytls x stderr2)
        end

        x = randn(n,n,n)
        x2 = copy(x)
        for L in (ndyadicscales(n),0,1,2)

            yf = dwt(x, wf, L)
            #if L == 5
            #    @profile dwt(x, wls, L)
            #    Profile.print()
            #end
            yls = dwt(x, wls, L)

            # filter vs lifting
            #@test (@vecnorm_eq_eps yf yls stderr2)
            @test yf ≈ yls rtol=stderr2

            ytf = idwt(yf, wf, L)
            ytls = idwt(yls, wls, L)

            # inverse vs original
            #@test (@vecnorm_eq_eps ytf x stderr2)
            #@test (@vecnorm_eq_eps ytls x stderr2)
            @test ytf ≈ x rtol=stderr2
            @test ytls ≈ x rtol=stderr2

        end
    end
end

@testset "Types and sizes" begin
    function makedwt(ft::Type, n, wf, L)
    	x0 = rand(-5:5, n)
    	x = zeros(ft, n)
    	copyto!(x,x0)
    	return (x, dwt(x, wf, L))
    end

    # 1-d
    n = 8
    wf = wavelet(WT.db2)
    L = 2

    x, y = makedwt(Float64, n, wf, L)
    @test Vector{Float64} == typeof(y) && length(y) == n
    x, y = makedwt(Float32, n, wf, L)
    @test Vector{Float32} == typeof(y) && length(y) == n
    x, y = makedwt(Int64, n, wf, L)
    @test Vector{typeof(float(x[1]))} == typeof(y) && length(y) == n
    x, y = makedwt(Int32, n, wf, L)
    @test Vector{typeof(float(x[1]))} == typeof(y) && length(y) == n
    @test dwt(x,wf) ≈ dwt(float(x),wf)

    # Complex types
    for wfc in (wavelet(WT.db2), wavelet(WT.db2, WT.Lifting))
        xc = zeros(ComplexF64, n)
        map!(i->rand(ComplexF64), xc, xc)
        yc = dwt(xc, wfc, L)
        @test Vector{ComplexF64} == typeof(yc) && length(yc) == n
        @test xc ≈ idwt(dwt(xc, wfc), wfc)
    end

    # 2-d
    sett = ((n,n),wf,L)
    ft = Float64; x, y = makedwt(ft, sett...)
    @test Array{ft,2} == typeof(y) && length(y) == n*n
    ft = Float32; x, y = makedwt(ft, sett...)
    @test Array{ft,2} == typeof(y) && length(y) == n*n
    ft = Int64; x, y = makedwt(ft, sett...)
    @test Array{typeof(float(x[1])),2} == typeof(y) && length(y) == n*n
    ft = Int32; x, y = makedwt(ft, sett...)
    @test Array{typeof(float(x[1])),2} == typeof(y) && length(y) == n*n
    @test dwt(x,wf) ≈ dwt(float(x),wf)

    # non-Array type
    wt = wavelet(WT.db2, WT.Lifting)
    x = randn(16, 16)
    xs = view(copy(x), 1:16, 1:16)
    @test dwt(x,wt,2) ≈ dwt(xs,wt,2)

    # util functions
    for class in (WT.haar, WT.db2, WT.cdf97)
        WT.class(class)
        WT.name(class)
        WT.vanishingmoments(class)
    end

    wt = wavelet(WT.db1, WT.Filter)
    @test length(wt) == 2
    @test wt.qmf*0.7 ≈ WT.scale(wt, 0.7).qmf

    # inplace methods
    wt = wavelet(WT.db1, WT.Filter)
    x = randn(8)
    y = similar(x)
    @test dwt!(y, x, wt) ≈ dwt(x, wt)

    wt = wavelet(WT.db1, WT.Lifting)
    x = randn(8)
    y = copy(x)
    @test dwt!(x, wt) ≈ dwt(y, wt)
end

struct wunknownt <: DiscreteWavelet{Float64} end
@testset "Errors" begin
    uwt = wunknownt()
    EE = Exception
    @test_throws EE dwt(randn(4),uwt)
    @test_throws EE dwt(randn(4,4),uwt)
    @test_throws EE wavelet(WT.Coiflet{33}())
    @test_throws EE wavelet("db2asdsad")
    @test_throws EE wavelet("db2", "ppppp")
end

# 2-d
function makedwt(ft::Type, n, wf, L)
    x0 = rand(-5:5, n)
    x = zeros(ft, n)
    copyto!(x,x0)
    return (x, dwt(x, wf, L))
end

n=8
wf = wavelet(WT.db2)
L = 2
sett = ((n,n),wf,L)
ft = Float64; x, y = makedwt(ft, sett...)
@test Array{ft,2} == typeof(y) && length(y) == n*n
ft = Float32; x, y = makedwt(ft, sett...)
@test Array{ft,2} == typeof(y) && length(y) == n*n
ft = Int64; x, y = makedwt(ft, sett...)
@test Array{typeof(float(x[1])),2} == typeof(y) && length(y) == n*n
ft = Int32; x, y = makedwt(ft, sett...)
@test Array{typeof(float(x[1])),2} == typeof(y) && length(y) == n*n
@test dwt(x,wf) ≈ dwt(float(x),wf)

# non-Array type
wt = wavelet(WT.db2, WT.Lifting)
x = randn(16, 16)
xs = view(copy(x), 1:16, 1:16)
@test dwt(x,wt,2) ≈ dwt(xs,wt,2)

# util functions
for class in (WT.haar, WT.db2, WT.cdf97)
    WT.class(class)
    WT.name(class)
    WT.vanishingmoments(class)
end
class = WT.db1
wt = wavelet(class, WT.Filter)
@test length(wt) == 2
@test wt.qmf*0.7 ≈ WT.scale(wt, 0.7).qmf

# inplace methods
class = WT.db1
wt = wavelet(class, WT.Filter)
x = randn(8)
y = similar(x)
@test dwt!(y, x, wt) ≈ dwt(x, wt)

wt = wavelet(class, WT.Lifting)
x = randn(8)
y = copy(x)
@test dwt!(x, wt) ≈ dwt(y, wt)

# continuous 1-d; different scalings should lead to different sizes, different boundary condtions shouldn't
@testset "Continuous Wavelet Transform" begin
    for xSize = (33, 67)
        for boundary = (WT.DEFAULT_BOUNDARY, WT.padded, WT.NaivePer)
            for s=1:2:8
                for wfc in (wavelet(WT.morl,s=s,boundary=boundary), wavelet(WT.dog0,s=s,boundary=boundary), wavelet(WT.paul4,s=s,boundary=boundary))
                    xc = rand(Float64,xSize)
                    yc = cwt(xc,wfc)
                    if typeof(wfc.waveType) <: Union{WT.Morlet, WT.Paul}
                        @test Array{ComplexF64,2}==typeof(yc)
                    else
                        @test Array{Float64,2}==typeof(yc)
                    end
                    nOctaves= log2(xSize) - wfc.averagingLength;
                    nWaveletsInOctave = reverse([max(1, round(Int, s / x^(1))) for
                                                 x=1:round(Int, nOctaves)])
                    totalWavelets = round(Int, sum(nWaveletsInOctave))
                    @test size(yc) == (xSize, totalWavelets+1)
                end
            end
        end
    end
end
# TODO: test actual values
#       test averaging types
#            various extra dimensions

@testset "WPT" begin
    wf = wavelet(WT.db2, WT.Filter)
    x = randn(16)

    L = 1
    wp = wpt(x,wf,L)
    dw = dwt(x,wf,L)
    @test wp ≈ dw
    @test iwpt(wp,wf,L) ≈ x

    L = 2
    wp = wpt(x,wf,L)
    dw = dwt(x,wf,L)
    dw2 = copy(dw)
    dw2[9:end] = dwt(dw[9:end],wf,1)
    @test dw[1:8] ≈ wp[1:8]
    @test dw2 ≈ wp
    @test iwpt(wp,wf,L) ≈ x

    L = 3
    wp = wpt(x,wf,L)
    dw = dwt(x,wf,L)
    @test dw[1:4] ≈ wp[1:4]
    @test dwt(dw2[5:8],wf,1) ≈ wp[5:8]
    @test dwt(dw2[9:12],wf,1) ≈ wp[9:12]
    @test dwt(dw2[13:16],wf,1) ≈ wp[13:16]
    @test iwpt(wp,wf,L) ≈ x

    wl = wavelet(WT.db2, WT.Lifting)
    x = randn(128)
    @test iwpt(wpt(x,wf),wf) ≈ x
    @test iwpt(wpt(x,wl),wl) ≈ x
    @test wpt(x,wl,1) ≈ wpt(x,wf,1)
    @test wpt(x,wl,2) ≈ wpt(x,wf,2)
    @test wpt(x,wl,4) ≈ wpt(x,wf,4)
    @test wpt(x,wl) ≈ wpt(x,wf)

    wl = wavelet(WT.db2, WT.Lifting)
    n = 128
    x = randn(n)
    for L = 0:ndyadicscales(n)
        @test wpt(x, wl, maketree(n, L, :dwt)) ≈ dwt(x, wl, L)
        @test iwpt(x, wl, maketree(n, L, :dwt)) ≈ idwt(x, wl, L)
        @test wpt(x, wf, maketree(n, L, :dwt)) ≈dwt(x, wf, L)
        @test iwpt(x, wf, maketree(n, L, :dwt)) ≈ idwt(x, wf, L)
    end

    # non-dyadic tests
    wt = wavelet(WT.db2, WT.Filter)
    n = 5*8
    x = randn(n)
    for L in 0:maxtransformlevels(n)
        for wt in (wavelet(WT.db2, WT.Filter), wavelet(WT.db2, WT.Lifting))
            @test wpt(x, wt, maketree(n, L, :dwt)) ≈ dwt(x, wt, L)
            @test iwpt(x, wt, maketree(n, L, :dwt)) ≈ idwt(x, wt, L)
        end
    end
end

@testset "MODWT" begin
    wf = wavelet(WT.db4)
    # power-of-two
    x = randn(128)
    W = modwt(x, wf)
    x_back = imodwt(W, wf)
    @test x ≈ x_back

    # non power-of-two
    x = cumsum(randn(129))
    W = modwt(x, wf)
    x_back = imodwt(W, wf)
    @test x ≈ x_back
    @test size(W) == (length(x), maxmodwttransformlevels(x)+1)

    # less than max number of levels
    L = 4
    Wl = modwt(x, wf, L)
    @test W[:, 1:L-1] ≈ Wl[:, 1:L-1]
end

@testset "Wavelet types" begin
    wt = wavelet(WT.db2, WT.Filter)
    @test length(wt) == 4
    same = WT.name(wt) == "db2"
    @test same

    wt = wavelet(WT.db2, WT.Lifting)
    @test WT.name(wt) == "db2"
end


text:
@testset "gmres" begin
	@testset "real matrix" begin
		A  = sprandn(100,100,.1) + 10*I
		n  = size(A,2)
		D  = Array(diag(A))
		M2 = x -> D.\x
		rhs = randn(100)
		tol = 1e-6;

		# test printing and behaviour for early stopping
		xtt = gmres(A,rhs ,3,tol=1e-10,maxIter=3,out=2,storeInterm=true)
		@test xtt[2]==-1

		# test behaviour for zero rhs
		xtt = gmres(A,0*rhs,5,tol=tol,maxIter=10,out=2)
		@test xtt[2]==-9
		@test all(xtt[1].==0)
		@test length(xtt[1])==100
		@test eltype(xtt[1])==eltype(rhs)

		x1 = gmres(A,rhs ,5,tol=tol,maxIter=100,out=1)
		x3 = gmres(A,rhs,5,tol=tol,maxIter=100,x=randn(size(rhs)))
		x4 = gmres(A,rhs,5,tol=tol,maxIter=100,M=M2)

		@test norm(A*x1[1]-rhs)/norm(rhs) < tol
		@test norm(A*x3[1]-rhs)/norm(rhs) < tol
		# @test norm(A*x4[1]-rhs)/norm(rhs) < tol
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5
	end

	@testset "complex matrix" begin
		A  = sprandn(100,100,.1) + 10*I + im*(sprandn(100,100,.1) + 10*I )
		D  = Array(diag(A))
		M3 = x -> D.\x
		rhs = complex(randn(100))
		tol = 1e-6;

		# test behaviour for zero rhs
		xtt = gmres(A,0*rhs,5,tol=tol,maxIter=10,out=2)
		@test xtt[2]==-9
		@test all(xtt[1].==0)
		@test length(xtt[1])==100
		@test eltype(xtt[1])==eltype(rhs)

		x1 = gmres(A,rhs ,5,tol=tol,maxIter=100)
		x3 = gmres(A,rhs,5,tol=tol,maxIter=100,x=randn(size(rhs)))
		x4 = gmres(A,rhs,5,tol=tol,maxIter=100,M=M3)

		@test norm(A*x1[1]-rhs)/norm(rhs) < tol
		@test norm(A*x3[1]-rhs)/norm(rhs) < tol
		@test norm(x3[1]-x1[1])/norm(x1[1]) < 1e-5
	end
end

ext: jl


text:
__precompile__()
module ProgressLoggers

export ProgressLogger,
       update!,
       flush!

@inline isatty(stream::IO) = isa(stream, Base.TTY)

immutable ProgressLogger
    sink::IO
    format::AbstractString
    interval::Int
    printer::Function
    noun::AbstractString
    starttime::UInt64

    function ProgressLogger(sink::IO, format::AbstractString, interval::Int,
                            noun::AbstractString)
        printer = isatty(sink) ? print : println
        new(sink, format, interval, printer, noun, time_ns())
    end
end

function ProgressLogger(sink::IO, interval::Int, noun::AbstractString)
    if isatty(sink)
        fmt = "\x1b[1G\x1b[2K- "
    else
        fmt = "    ... "
    end
    ProgressLogger(sink, fmt, interval, noun)
end

function ProgressLogger(interval::Int, noun::AbstractString)
    ProgressLogger(STDERR, interval, noun)
end

function ProgressLogger(interval::Int)
    ProgressLogger(STDERR, interval, "item")
end

function humanreadable(count::Real)
    if count < 1
        logc = 0
    else
        logc = ceil(log10(count))
    end
    if logc <= 3
        # Special case: Just return the
        return @sprintf("%d", count)
    elseif logc <= 6
        unit = "K"
        div = 1000
    elseif logc <= 9
        unit = "M"
        div = 1000000
    elseif logc <= 12
        unit = "G"
        div = 1000000000
    else
        unit = "T"
        div =  1000000000000
    end
    return @sprintf("%0.2f%s", count/div, unit)
end

function flush!(log::ProgressLogger, item::Int, etc...)
    if log.interval > 0
        elapsed = (time_ns() - log.starttime) / 1e9
        persec = humanreadable(item / elapsed)
        item = humanreadable(item)
        println(log.sink, log.format, item,
                " $(log.noun) ($(persec) / sec)", etc...)
    end
end

function update!(log::ProgressLogger, item::Int, etc...)
    if log.interval > 0 && item % log.interval == 0
        item = humanreadable(item)
        log.printer(log.sink, log.format, item, " $(log.noun)", etc...)
    end
end

end # module


text:
export Var, asvar, data, gradient, initgrad!, cleargrad!, isgraddefined, addname!, name
export Param, isparam

mutable struct Var{T <: AbstractArray,F <: Func}
    data::T
    creator::F
    grad::Union{Var,Nothing}
end

function Var(data::T, creator::F=nullfunc) where {T <: AbstractArray,F <: Func}
    Var{T,F}(data, creator, nothing)
end
Var(data::Number, creator=nullfunc) = Var([data], creator)

asvar(x) = x isa Var ? x : Var(x)
asvar(x, creator) = Var(x, creator)

data(x::Var) = x.data

gradient(x::Var) = x.grad

Base.size(x::Var) = size(x.data)

Base.size(x::Var, i::Int) = size(x.data, i)

Base.length(x::Var) = length(x.data)

Base.ndims(x::Var) = ndims(x.data)

Base.eltype(x::Var) = eltype(x.data)

function Base.show(io::IO, v::Var)
    print(io, "Var($(v.data))")
end

function cleargrad!(x::Var)
    x.grad = nothing
end

cleargrad!(x::Var...) = cleargrad!.(x)

function initgrad!(x::Var, v=1.0)
    x.grad = Var(fill(v, size(x.data)))
end

isgraddefined(x::Var) = !isnothing(x.grad)

const varname = WeakKeyDict{Var, String}()

function addname!(x::Var, name)
    varname[x] = name
end

getname(x::Var) = haskey(varname, x) ? varname[name] : ""

Param(data) = Var(data, paramcreator)

isparam(x) = x isa Var && x.creator isa ParamCreator


text:
# Susceptibility.jl
#  Frequency-dependent susceptibility (and thus mobility) of the Feynman polaron state
#    A WORK IN PROGRESS.

# Data structure to store the results
struct susceptibility
    nu
    ImX
    μ
end
Susceptibility()=susceptibility([],[],[])

"""
function ImX(nurange,v,w,βred,α,ω,mb)

Impedance in (47a) from Feynman1962, directly solving freq dep without taking
Hellwarth1999 limit of v->0 .

Calculates a frequency dependent (over range of nu) susceptibility which can be linked back to mobility.

HERE BE DRAGONS!
Not well tested or validated code; the central numeric integration is highly
oscillatory and becomes intractable for large nu.
"""
function ImX(nurange,v,w,βred,α,ω,mb)
        @printf("\nAin't no Imaginary Impedance like a Feynman 1962 ImX...\n")
        println("ImX Call signature: v: $v w: $w βred: $βred α: $α ω: $ω mb: $mb")
    # Feynman, I love you - but using Nu, v; Omega, w in the same paper + formulas, for similar objects?!
    s=Susceptibility()
    for nu in nurange
        R=(v^2-w^2)/(w^2*v)     # FHIP1962, page 1011, eqn (47c). Note this is wrong in some textbooks / 1990s PRB.
        b=R*βred/sinh(βred*v/2) # FHIP1962, page 1010, eqn (47b)
        a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2)) # FHIP1962, page 1010, eqn (47b)
        k(u,a,b,v,nu) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u)*cos(nu*u) # integrand with cos(vu) term, as (47a)

        @printf("Numerical integration of FHIP1962(42a): nu=%.2f ",nu)
        println("   a: $(a) b:$(b) ")

        # This is a simple approximation of part of the
        # https://www.gnu.org/software/gsl/doc/html/integration.html#qawf-adaptive-integration-for-fourier-integrals
        # GSL adaptive method for Fourier integrals
        # We split the integral into a load of integrals, balanced at some of
        # the roots.
        c=(2*floor(nu)+1)*π/nu
        println("Fourier integral c: $(c)")
        fourier_range = [ c*i for i in 0:2501 ]

        # Catch Inf c range
        if nu==0
            fourier_range = [0, Inf]
        end

        # These params tweaked to get the best behaviour at reproducing
        # Mishchenko-Fig4 in a timely manner.
        @time n=quadgk(u->k(u,a,b,v,nu),fourier_range... ,
                       maxevals=10^6,rtol=0.0001, atol=1e-15, order=7) # numerical quadrature integration of (2)
        K=n[1]
        err=n[2]
        @printf(" quadgk: K=%g err=%g\n",K,err)

        if K<err # we've lost control of our errors, due to losing the oscillatory war with nu
            break # --> so give up.
        end

        # Full 47a constructed here
        ImX= 2*α/(3*sqrt(π)) * βred^(3/2) * (sinh(βred*nu/2))/sinh(βred/2) * (v^3/w^3) * K

        μ=ImX^-1 * (q)/(ω*mb)

        @printf(" %.3f %g %g\n",nu,ImX,μ)

        append!(s.nu,nu)
        append!(s.ImX,ImX)
        append!(s.μ,μ)
    end

    @printf("\n\n")
    return(s)
end