text:
##
#
# Visualization and related functions
#

"""

    RootsImage(leftbottom, topright, width, height)

Essentially a two dimensional histogram covering a rectangle in the
complex plane specified by complex numbers `leftbottom` and `topright`
with `width` horizontal and `height` vertical bins.
"""
struct RootsImage{T <: Real}
  remin::T
  remax::T
  immin::T
  immax::T
  height::Int
  width::Int
  hbin::T
  wbin::T
  data::Array{T, 2}

  function RootsImage(leftbottom::Complex{T}, topright::Complex{T}, width::Int, height::Int) where {T <: Real}
    remin = real(leftbottom)
    immin = imag(leftbottom)
    remax = real(topright)
    immax = imag(topright)
    new{T}(remin, remax, immin, immax, height, width, (remax - remin) / width, (immax - immin) / height, zeros(T, width, height))
  end
end

"""

    add_root!(img, z)

Add a root `z` to the image `img`.
"""
function add_root!(img::RootsImage{T}, z::Complex{T}) where {T <: Real}
  r = real(z); i = imag(z)

  if r < img.remin || r >= img.remax || i < img.immin || i >= img.immax
    # out of range, skipping...
    return
  end

  x = min(max(1, 1 + floor(Int64,(r - img.remin) / img.hbin)), img.width)
  y = min(max(1, 1 + floor(Int64,(i - img.immin) / img.wbin)), img.height)

  img.data[x, y] += T(1)
end

"""

    plot(img[; mode])

Show a graphical representation of the image `img`.
The optional parameter `mode` controls how to scale number of roots in
histogram bins:

 * `:sharp`,
 * `:log`,
 * `:log_cutoff`,
 * `:linear`,
 * a function.

"""
function plot(img::RootsImage{T}; mode=:sharp) where {T <: Real}
  if isa(mode, Function)
    output = mode.(img.data)
  elseif mode == :sharp
    output = (x -> min(1, x)).(img.data)
  elseif mode == :log_cutoff
    m      = mean(img.data)
    output = (x -> log(1 + min(2m, x))).(img.data)
  elseif mode == :log
    output = (x -> log(1 + x)).(img.data)
  elseif mode == :linear
    output = img.data
  else
    error("Unknown mode!")
  end

  # normalization to [0, 1]
  output = output / maximum(output)

  return colorview(Gray, transpose(output)[end:-1:1,1:end])
end

"""

    save_image(img, filename[; mode])

Save the `img` to file `filename`.
Optional parameter `mode` has the same meaning as in the [`plot`](@ref) function.

"""
function save_image(img::RootsImage{T}, filename::AbstractString; mode=:log_cutoff, latex_filename=nothing) where {T <: Real}
  @info "Saving the image $filename..."

  if isa(mode, Function)
    output = mode.(img.data)
  elseif mode == :sharp
    output = (x -> min(1, x)).(img.data)
  elseif mode == :log_cutoff
    m      = mean(img.data)
    output = (x -> log(1 + min(2m, x))).(img.data)
  elseif mode == :log
    output = (x -> log(1 + x)).(img.data)
  elseif mode == :linear
    output = img.data
  else
    error("Unknown mode!")
  end

  output = output / maximum(output)

  Images.save(filename, colorview(Gray, transpose(output)[end:-1:1,1:end]))

  if latex_filename != nothing
    @info "Constructing the LaTeX file $latex_filename..."
    write_latex(img, filename, latex_filename)
  end
end

"""

    fill_images!(imgs, filename)

Read roots from file `filename` and draw them onto images contained
in the array `imgs`.

"""
function fill_images!(imgs::Array{RootsImage{T}}, filename::AbstractString) where {T <: Real}
  roots_count = div(filesize(filename), sizeof(Complex{T})) + 1
  progress = Progress(roots_count, 10)

  @info "Sifting through data (there are $roots_count roots)..."

  open(filename) do io
    while !eof(io)
      z = read(io, Complex{T})
      for img in imgs
        add_root!(img, z)
      end
      next!(progress)
    end
  end
end

"""

  animate_growth(poly_iter, img, filename[; output_filename])

Take roots of polynomials generated by `poly_iter` and saved in file `filename`
and draw them incrementaly (each step corresponds to change of degree) onto
image `img`.

"""
function animate_growth(poly_iter::PolynomialIterator{S}, img::RootsImage{T}, filename; output_filename=j->"animation_$j.png") where {S <: Number, T <: Real}
  progress = Progress(div(filesize(filename), sizeof(Complex{T})) + 1, 10)
  degree = 0

  open(filename) do io
    for p in poly_iter
      # save the current snapshot
      if degree != length(p) - 1
        save_image(img, output_filename(degree))
        degree += 1
      end

      # read roots
      for k in 1:degree
        z = read(io, Complex{T})
        add_root!(img, z)
        next!(progress)
      end
    end
  end
end

"""

    fill_image!(img, filename)

An analog of [`fill_images!`](@ref), only with one image.

"""
function fill_image!(img::RootsImage{T}, filename::AbstractString) where {T <: Real}
  fill_images!([img], filename)
end

"""

  write_latex(img, imagename, latex_filename)

Create a LaTeX PGFPlots version of the image with axes.

"""
function write_latex(img::RootsImage{T}, imagename::AbstractString, latex_filename::AbstractString) where {T <: Real}
  open(latex_filename, "w") do io
    write(io, "\\documentclass{standalone}\n")
    write(io, "\\usepackage[T1]{fontenc}\n")
    write(io, "\\usepackage[utf8]{inputenc}\n")
    write(io, "\\usepackage{pgfplots}\n")
    write(io, "\\usepackage{lmodern}\n")
    write(io, "\\begin{document}\n")
    write(io, "\\begin{tikzpicture}\n")
    write(io, "\\begin{axis}[enlargelimits=false, axis on top, xlabel={\$\\Re\$}, ylabel={\$\\Im\$}]\n")
    write(io, "\\addplot graphics [xmin=$(img.remin), xmax=$(img.remax), ymin=$(img.immin), ymax=$(img.immax)] {$imagename};\n")
    write(io, "\\end{axis}\n")
    write(io, "\\end{tikzpicture}\n")
    write(io, "\\end{document}\n")
  end
end


text:
using GPUifyLoops

kernel(A::Array) = kernel(Val(:CPU), A)
function kernel(::Val{Dev}, A) where Dev
    @setup Dev

    @loop for i in (1:size(A,1);
                    threadIdx().x)
        A[i] = 2*A[i]
    end
    @synchronize
end

data = Array{Float32}(undef, 1024)
kernel(data)

@static if Base.find_package("CuArrays") !== nothing
    using CuArrays
    using CUDAnative

    @eval function kernel(A::CuArray)
        @cuda threads=length(A) kernel(Val(:GPU), A)
    end

    data = CuArray{Float32}(undef, 1024)
    kernel(data)
end

ext: jl


text:
module WT
export
    DiscreteWavelet,
    ContinuousWavelet,
    FilterWavelet,
    LSWavelet,
    OrthoFilter,
    GLS,
    CFW,
    wavelet,
    boundaryType,
    eltypes,
    computeWavelets,
    Average,
    Father,
    Dirac,
    Daughter,
    findAveraging,
    getNWavelets,
    getScales,
    show

using ..Util
import Base.length
using SpecialFunctions
using LinearAlgebra

# TYPE HIERARCHY

abstract type DiscreteWavelet{T} end
abstract type ContinuousWavelet{Boundary,T} end
# discrete transforms via filtering
abstract type FilterWavelet{T} <: DiscreteWavelet{T} end
# discrete transforms via lifting
abstract type LSWavelet{T} <: DiscreteWavelet{T} end
# all wavelet types
#const WaveletTransformType = Union{DiscreteWavelet, ContinuousWavelet}

"""Get wavelet type name."""
function name(::DiscreteWavelet) end
"""Get wavelet filter length."""
function length(::FilterWavelet) end

struct FilterTransform end
struct LiftingTransform end
"""Transform by filtering."""
const Filter = FilterTransform()
"""Transform by lifting."""
const Lifting = LiftingTransform()

# BOUNDARY TYPES

abstract type WaveletBoundary end
# periodic (default)
struct PerBoundary <: WaveletBoundary end
# zero padding
struct ZPBoundary <: WaveletBoundary end
# constant padding
#struct CPBoundary <: WaveletBoundary end
struct NullBoundary <: WaveletBoundary end
# symmetric boundary (as in the DCTII)
struct SymBoundary <: WaveletBoundary end
# and so on...

const Periodic = PerBoundary()
const DEFAULT_BOUNDARY = PerBoundary()
const padded = ZPBoundary()
const NaivePer = NullBoundary()
const SymBound = SymBoundary()

# WAVELET CLASSES

"""
The `WaveletClass` type has subtypes `OrthoWaveletClass`,
 `BiOrthoWaveletClass`, and `ContinuousWaveletClass`

The `WT` module has for convenience constants defined named
as the class short name and optionally amended with a number
specifing the number of vanishing moments.

A class can also be explicitly constructed as e.g. `Daubechies{4}()`.

# Examples
`WT.db2`, `WT.haar`, `WT.cdf97`

**See also:** `WT.class`, `WT.name`, `WT.vanishingmoments`
"""
abstract type WaveletClass end
abstract type OrthoWaveletClass <: WaveletClass end
abstract type BiOrthoWaveletClass <: WaveletClass end
abstract type ContinuousWaveletClass <: WaveletClass end

# Single classes, orthogonal
for (TYPE, NAMEBASE, MOMENTS) in (
        (:Haar, "haar", 1),
        (:Beylkin, "beyl", -1), # TODO moments
        (:Vaidyanathan, "vaid",-1), # TODO moments
        )
    @eval begin
        struct $TYPE <: OrthoWaveletClass end
        class(::$TYPE) = $(string(TYPE))
        name(::$TYPE) = string($NAMEBASE)
        vanishingmoments(::$TYPE) = $MOMENTS
    end
    CONSTNAME = Symbol(NAMEBASE)
    @eval begin
        const $CONSTNAME = $TYPE()                  # type shortcut
    end
end

struct Morlet <: ContinuousWaveletClass
    σ::Float64 # \sigma is the time/space trade-off. as sigma->0, the spacial resolution increases; below 5, there is a danger of being non-analytic. Default is 5.8
    κσ::Float64
    cσ::Float64
end
"""
    morl = Morlet(σ::T) where T<: Real

    return the Morlet wavelet with first central frequency parameter σ, which controls the time-frequency trade-off. As σ goes to zero, all of the information becomes spatial. Default is `` 2π``, and it is recommended that you stay within 3-12. Above this range, the overlap between wavelets becomes impractically small. Below it, the mean subtraction term approaches the magnitude of the wavelet, so they become a sum of Gaussians rather than Gaussians.
"""
function Morlet(σ::T) where T<:Real
    κσ=exp(-σ^2/2)
    cσ=1. /sqrt(1+κσ^2-2*exp(-3*σ^2/4))
    Morlet(σ,κσ,cσ)
end
Morlet() = Morlet(2π)
class(::Morlet) = "Morlet"; name(::Morlet) = "morl"; vanishingmoments(::Morlet)=0
const morl = Morlet()
Base.show(io::IO, x::Morlet) = print(io,"Morlet mean $(x.σ)")

# TODO: include a sombrero wavelet, which is dog2.

# Parameterized classes

# abstract type Dog <: ContinuousWaveletClass end
# abstract type Paul <: ContinuousWaveletClass end

# continuous parameterized
for (TYPE, NAMEBASE, MOMENTS, RANGE) in (
        (:Paul, "paul", -1, 1:20), # moments? TODO: is this a good range of parameters?
        (:Dog, "dog",  -1, 0:6), # moments?
        )
    @eval begin
        struct $TYPE{N} <: ContinuousWaveletClass end#$TYPE end
        class(::$TYPE) = $(string(TYPE))
        name(::$TYPE{N}) where N = string($NAMEBASE,N)
        vanishingmoments(::$TYPE{N}) where N = -1
        order(::$TYPE{N}) where N = N # either order for Paul wavelets, or number of derivatives for DOGs
        Base.show(io::IO, x::$TYPE{N}) where N = print(io,$NAMEBASE * " order $(N)")
    end
    for NUM in RANGE
        CONSTNAME = Symbol(NAMEBASE,NUM)
        @eval begin
            const $CONSTNAME = $TYPE{$NUM}()                  # type shortcut
        end
    end
end

# discrete ortho classes
for (TYPE, NAMEBASE, RANGE) in (
        (:Daubechies, "db", 1:10),
        (:Coiflet, "coif", 2:2:8),
        (:Symlet, "sym", 4:10),
        (:Battle, "batt", 2:2:6),
        )
    @eval begin
        struct $TYPE{N} <: OrthoWaveletClass end
        class(::$TYPE) = $(string(TYPE))
        name(::$TYPE{N}) where N = string($NAMEBASE,N)
        vanishingmoments(::$TYPE{N}) where N = N
    end
    for NUM in RANGE
        CONSTNAME = Symbol(string(NAMEBASE, NUM))
        @eval begin
            const $CONSTNAME = $TYPE{$NUM}()        # type shortcut
        end
    end
end

# Parameterized BiOrtho classes
for (TYPE, NAMEBASE, RANGE1, RANGE2) in (
        (:CDF, "cdf", [9], [7]),
        )
    @eval begin
        struct $TYPE{N1, N2} <: BiOrthoWaveletClass end
        class(::$TYPE) = $(string(TYPE))
        name(::$TYPE{N1, N2}) where {N1, N2} = string($NAMEBASE,N1,"/",N2)
        vanishingmoments(::$TYPE{N1, N2}) where {N1, N2} = (N1, N2)
    end
    for i in length(RANGE1)
        CONSTNAME = Symbol(string(NAMEBASE,RANGE1[i],RANGE2[i]))
        @eval begin
            const $CONSTNAME = $TYPE{$RANGE1[$i],$RANGE2[$i]}()        # type shortcut
        end
    end
end

# IMPLEMENTATIONS OF FilterWavelet

"""
Wavelet type for discrete orthogonal transforms by filtering.

**See also:** `GLS`, `wavelet`
"""
struct OrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T}
    qmf     ::Vector{Float64}        # quadrature mirror filter
    name    ::String                 # filter short name
end

function OrthoFilter(w::WC, ::T=DEFAULT_BOUNDARY) where {WC<:OrthoWaveletClass, T<:WaveletBoundary}
    name = WT.name(w)
    if WC <: Daubechies
        qmf = daubechies(vanishingmoments(w))
    else
        qmf = get(FILTERS, name, nothing)
        qmf == nothing && throw(ArgumentError("filter not found"))
    end
    # make sure it is normalized in l2-norm
    return OrthoFilter{T}(qmf./norm(qmf), name)
end

length(f::OrthoFilter) = length(f.qmf)
qmf(f::OrthoFilter) = f.qmf
name(f::OrthoFilter) = f.name

"""Scale filter by scalar."""
function scale(f::OrthoFilter{T}, a::Number) where T<:WaveletBoundary
    return OrthoFilter{T}(f.qmf.*a, f.name)
end

"""Quadrature mirror filter pair."""
function makeqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f)))
    scfilter, dcfilter = makereverseqmfpair(f, fw, T)
    return reverse(scfilter), reverse(dcfilter)
end

"""Reversed quadrature mirror filter pair."""
function makereverseqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f)))
    h = convert(Vector{T}, qmf(f))
    if fw
        scfilter = reverse(h)
        dcfilter = mirror(h)
    else
        scfilter = h
        dcfilter = reverse(mirror(h))
    end
    return scfilter, dcfilter
end

#struct BiOrthoFilter{T<:WaveletBoundary} <: FilterWavelet{T}
#    qmf1::Vector{Float64}       # quadrature mirror filter 1
#    qmf2::Vector{Float64}       # quadrature mirror filter 2
#    name::String                # filter short name
#    BiOrthoFilter(qmf1, qmf2, name) = new(qmf1, qmf2, name)
#end

# IMPLEMENTATIONS OF LSWavelet

abstract type StepType end
struct PredictStep <: StepType end
struct UpdateStep <: StepType end
const Predict = PredictStep()
const Update = UpdateStep()

struct LSStepParam{T<:Number}
    coef    ::Vector{T}        # lifting coefficients
    shift   ::Int              # + left shift, - right shift
end

struct LSStep{T<:Number}
    param::LSStepParam{T}
    steptype::StepType
end

function LSStep(st::StepType, coef::Vector{T}, shift::Int) where T
    return LSStep{T}(LSStepParam{T}(coef, shift), st)
end

length(s::LSStep) = length(s.param)
length(s::LSStepParam) = length(s.coef)

"""
Wavelet type for discrete general (bi)orthogonal transforms
by using a lifting scheme.

**See also:** `OrthoFilter`, `wavelet`, `CWT`
"""
struct GLS{T<:WaveletBoundary} <: LSWavelet{T}
    step    ::Vector{LSStep{Float64}}    # steps to be taken
    norm1   ::Float64           # normalization of scaling coefs.
    norm2   ::Float64           # normalization of detail coefs.
    name    ::String            # name of scheme
end

function GLS(w::WC, ::T=DEFAULT_BOUNDARY) where {WC<:WaveletClass, T<:WaveletBoundary}
    name = WT.name(w)
    schemedef = get(SCHEMES, name, nothing)
    schemedef == nothing && throw(ArgumentError("scheme not found"))
    return GLS{T}(schemedef[1], schemedef[2], schemedef[3], name)
end

name(s::GLS) = s.name

######################################
######################################
######################################
######################################
# IMPLEMENTATIONS OF ContinuousWavelet

abstract type Average end
struct Dirac <: Average end
struct Father <: Average end
struct NoAve <: Average end

function Base.show(io::IO, cf::Dirac)
    print(io, "Dirac")
end
function Base.show(io::IO, cf::Father)
    print(io, "Father Wavelet")
end
function Base.show(io::IO, cf::NoAve)
    print(io, "No Averaging")
end

struct CFW{W, S, WT<:ContinuousWaveletClass, N} <: ContinuousWavelet{W, S}
    scalingFactor::S # the number of wavelets per octave, ie the scaling
                           # is s=2^(j/scalingfactor)
    decreasing::S # the amount that scalingFactor decreases per octave,
                        # to a minimum of 1
    fourierFactor::S
    coi          ::S
    α            ::Int   # the order for a Paul and the number of derivatives
                         # for a DOG
    σ            ::Array{S} # the morlet wavelet parameters
                            # (σ,κσ,cσ). NaN if not morlet.
    waveType     ::WT        # because multiple dispatch is good
    averagingLength::Int # the number of scales to override with averaging. If
                         # you want no averaging set it to zero
    averagingType  ::Average # either Dirac or mother; the first uniformly
                             # represents the lowest frequency information,
                             # while the second constructs
                             # a wavelet using Daughter that has mean frequency
                             # zero, and std equal to
                             # the first non-removed wavelet's mean
    frameBound     ::S       # if positive, set the frame bound of the
                             # transform to be frameBound. Otherwise leave it
                             # so that each wavelet has an L2 norm of 1
    normalization  ::N       # the normalization that is preserved for the
                             # wavelets as the scale changes. The conjugate
                             # p-norm is preserved for the  signal. Should be
                             # larger than 1
end

"""
    CFW(wave::WC, scalingFactor::S=8.0, averagingType::Symbol=:Father,
        boundary::T=WT.DEFAULT_BOUNDARY, averagingLength::Int =
        floor(Int,2*scalingFactor), frameBound::Float64=1.0,
        normalization::Float=Inf) where {WC<:WT.WaveletClass,
        T<:WT.WaveletBoundary, S<:Real}

The constructor for the CFW struct. w is a type of continuous wavelet,
scalingFactor is the number of wavelets between the octaves ``2^J`` and
``2^{J+1}`` (defaults to 8, which is most appropriate for music and other
audio). As this leads to excessively many high scale wavelets, `decreasing`
gives the amount that scalingFactor decreases per octave. The default boundary
condition is `periodic`, which is implemented by appending a flipped version of
the vector at the end (to eliminate edge discontinuities). Alternatives are
`ZPBoundary`, which pads with enough zeros to get to the nearest power of 2
(here the results returned by caveats are relevant, see Torrence and Compo
'97), and `NullBoundary`, which assumes the data is inherently periodic.

`averagingLength` and `averagingType` determine how wide scale information is
accounted for. `averagingLength` gives the number of wavelet octaves that are
covered by the averaging, while averaging type determines whether it is a
window `Dirac` or a wavelet `Father`. `frameBound` gives the total norm of the
whole collection, corresponding to the upper frame bound. `normalization`
refers to which p-norm is preserved as the scale changes. `normalization==2` is
the default scaling, while `normalization==Inf` gives all the same maximum
value, thus acting more like windows.
"""
function CFW(wave::WC, scalingFactor=8, boundary::T=WT.DEFAULT_BOUNDARY,
             averagingType::A = NoAve(),
             averagingLength::Int = 0,
             frameBound=-1, normalization::N=Inf,
             decreasing=1) where {WC<:WaveletClass, A <: Average,
                                  T <: WaveletBoundary, N <: Real}
    @assert scalingFactor > 0
    @assert normalization >= 1
    nameWavelet = name(wave)[1:3]
    tdef = get(CONT_DEFS, nameWavelet, nothing)
    tdef == nothing && error("transform definition not found; you gave $(nameWavelet)")
    # do some substitution of model parameters
    if nameWavelet=="mor"
        tdef = [eval(Meta.parse(replace(tdef[1],"σ" => wave.σ))), eval(Meta.parse(tdef[2])), -1, [wave.σ,wave.κσ,wave.cσ], wave]
    elseif nameWavelet[1:3]=="dog" || nameWavelet[1:3]=="pau"
        tdef = [eval(Meta.parse(replace(tdef[1], "α" => order(wave)))), eval(Meta.parse(replace(tdef[2], "α"=> order(wave)))), order(wave), [NaN], wave]
    else
        error("I'm not sure how you got here. Apparently the WaveletClass you gave doesn't have a name. Sorry about that")
    end
    if averagingLength <= 0 || typeof(averagingType) <: NoAve
        averagingLength = 0
        averagingType = NoAve()
    end
    # S is the most permissive type of the listed variables
    S = promote_type(typeof(scalingFactor), typeof(decreasing),
                     typeof(frameBound), typeof(normalization),
                     typeof(tdef[1]), typeof(tdef[2]))
    return CFW{T, S, WC, N}(S(scalingFactor), S(decreasing), tdef...,
                            averagingLength, averagingType, S(frameBound),
                            S(normalization))
end
name(s::CFW) = name(s.waveType)

function eltypes(::CFW{W, T, WT, N}) where {W, T, WT, N}
    T
end
function boundaryType(::CFW{W, T, WT, N}) where {W, T, WT, N}
    W
end
function waveletType(::CFW{W, T, WT, N}) where {W, T, WT, N}
    WT
end

function Base.show(io::IO, cf::CFW{W,S,WT,N}) where {W,S,WT,N}
    print("CFW[$(cf.waveType), $(cf.averagingType), decreasing rate = "*
          "$(cf.decreasing), aveLen = $(cf.averagingLength), frame = "*
          "$(cf.frameBound), norm=$(cf.normalization)]")
end

"""
    daughter = Daughter(this::CFW{W, T, <:ContinuousWaveletClass, N},
                        s::Real, nInOctave::Int, ω::AbstractArray{<:Real,1}) where {W, T, N}

given a CFW object, return a rescaled version of the mother wavelet, in the
fourier domain. ω is the frequency, which is fftshift-ed. s is the scale
variable.
"""
function Mother(this::CFW{W, T, Morlet, N}, s::Real, sWidth::Real,
                  ω::AbstractArray{<:Real,1}) where {W, T, N}
    constant = this.σ[3]*(π)^(1/4)
    gauss = exp.(-(this.σ[1].-ω/s).^2/(sWidth^2))
    shift = this.σ[2]*exp.(-1/2*(ω/s).^2)
    daughter = constant .* (gauss .- shift)
    return normalize(daughter, s, this.normalization)
end

function Mother(this::CFW{W, T, <:Paul, N}, s::Real, sWidth::Real,
                  ω::AbstractArray{<:Real,1}) where {W, T, N}
    daughter = zeros(length(ω))
    constant = (2^this.α) / sqrt((this.α) * gamma(2*(this.α)))
    polynomial = (ω[ω.>=0]/s).^(this.α)
    expDecay = exp.(-(ω[ω.>=0]/s))
    daughter[ω.>=0]= constant .* polynomial .* expDecay
    return normalize(daughter, s, this.normalization)
end

function Mother(this::CFW{W, T, <:Dog, N}, s::Real, sWidth::Real,
                  ω::AbstractArray{<:Real,1}) where {W, T, N}
    constant = im^(this.α)*sqrt(gamma((this.α)+1/2))
    polynomial = (ω/s).^(this.α)
    gauss = exp.(-(ω/s).^2/2)
    daughter =  constant .* polynomial .* gauss
    return normalize(daughter, s, this.normalization)
end

function normalize(daughter, s, p)
    if p == Inf
        normTerm = maximum(abs.(daughter))
    else
        normTerm = s^(1/p)
    end
    return daughter ./ normTerm
end

"""
    findAveraging(c::CFW, ω, averagingType::aT) where {aT<:Averaging}

this creates the averaging function, which covers the low frequency
information, and is emphatically not analytic. aT determines whether it has the
same form as the wavelets (`Father`), or just a bandpass `Dirac`.
`c.averagingLength` gives the number of octaves (base 2) that are covered by
the averaging function. The width is then derived so that it matches the next
wavelet at 1σ.

For the Morlet wavelet, the distribution is just a Gaussian, so it has variance
1/s^2 and mean σ[1]*s set the variance so that the averaging function has σ/2 at
the central frequency of the last scale.

For the Paul wavelets, it's a easy calculation to see that the mean of a paul
wavelet of order m is (m+1)/s, while σ=sqrt(m+1)/s. So we set the variance so
that the averaging function has 1σ at the central frequency of the last scale.

the derivative of a Gaussian (Dog) has a pretty nasty form for the mean and
variance; eventually, if you set σ_{averaging}=⟨ω⟩_{highest scale wavelet}, you
will get the scale of the averaging function to be
`s*gamma((c.α+2)/2)/sqrt(gamma((c.α+1)/2)*(gamma((c.α+3)/2)-gamma((c.α+2)/2)))`

"""
function findAveraging(c::CFW, ω, averagingType::Father, sWidth)
    s = 2^(c.averagingLength-1)
    s0, ω_shift = locationShift(c, s, ω, sWidth)
    averaging = Mother(c, s0, 1, ω_shift)
end

function locationShift(c::CFW{W, T, <:Morlet, N}, s, ω, sWidth) where {W,T,N}
        s0 = 3*c.σ[1]/4 *s*sWidth
        ω_shift = ω .+ c.σ[1] * s0
    return (s0, ω_shift)
end

function locationShift(c::CFW{W, T, <:Dog, N}, s, ω, sWidth) where {W,T,N}
    s0 = s*gamma((c.α+2)/2) / sqrt(gamma((c.α+1)/2)) * (gamma((c.α+3)/2) -
                                                        gamma((c.α+2)/2))
    μ = sqrt(2)*s0*gamma((c.α+2)/2)/gamma((c.α+1)/2)
    ω_shift = ω .+ μ
    return (s0, ω_shift)
end

function locationShift(c::CFW{W, T, <:Paul, N}, s, ω, sWidth) where {W,T,N}
        s0 = s*sqrt(c.α+1)
        ω_shift = ω .+ (c.α .+ 1) * s0
    return (s0, ω_shift)
end

# dirac version (that is, just a window around zero)
function findAveraging(c::CFW{<:WaveletBoundary, T}, ω,
                       averagingType::Dirac, sWidth) where {T}
    s = 2^(c.averagingLength) * sWidth
    averaging = zeros(T, size(ω))
    upperBound = getUpperBound(c, s)
    averaging[abs.(ω) .<= upperBound] .= 1
    return averaging
end

function getUpperBound(c::CFW{W, T, <:Morlet, N}, s) where {W,T,N}
    return c.σ[1] * s
end

function getUpperBound(c::CFW{W, T, <:Dog, N}, s) where {W,T,N}
    return sqrt(2)* s * gamma((c.α+2)/2) / gamma((c.α+1)/2)
end

function getUpperBound(c::CFW{W, T, <:Paul, N}, s) where {W,T,N}
    return (c.α + 1) * s
end

"""
    nOctaves, totalWavelets, sRanges, sWidth = getNWavelets(n1,c)

utility for understanding the spacing of the wavelets. `sRanges` is a list of
the s values used in each octave. sWidth is a list of the corresponding
variance adjustments
"""
function getNWavelets(n1,c)
    nOctaves = getNOctaves(n1,c)
    n = setn(n1,c)
    isAve = (c.averagingLength > 0 && !(typeof(c.averagingType) <: NoAve)) ? 1 : 0
    if round(nOctaves) < 0
        totalWavelets = 0
        sRanges = Array{Array{Float64,1},1}(undef,0)
        sWidth = [1.0]
        return nOctaves, totalWavelets, sRanges, sWidth
    end
    sRange = 2 .^ (polySpacing(nOctaves, c))
    totalWavelets = round(Int, length(sRange) + isAve)
    sWidth = varianceAdjust(c,totalWavelets)
    return nOctaves, totalWavelets, sRange, sWidth
end

function getNOctaves(n1,c::CFW{W,T, Morlet, N}) where {W, T, N}
    # choose the number of octaves so the last mean, which is at s*σ[1]
    # is 3 standard devations away from the end
    nOctaves = log2(max(n1*2π/(c.σ[1]+3), 2)) - c.averagingLength
end

function getNOctaves(n1,c::CFW{W,T, M, N}) where {W, T, N, M}
    nOctaves = log2(max(n1*2π, 2)) - c.averagingLength - 1
end

function varianceAdjust(this::CFW{W,T, M, N}, totalWavelets) where {W,T,N, M}
    # increases the width of the wavelets by σ[i] = (1+a(total-i)ᴾ)σₛ
    # which is a polynomial of order p going from 1 at the highest frequency to
    # sqrt(p) at the lowest
    p = this.decreasing
    a = (p^.8-1)/(totalWavelets-1)^p
    1 .+a .*(totalWavelets .- (1:totalWavelets)).^p
end

function polySpacing(nOct, c)
    a =c.averagingLength; O = nOct
    p = c.decreasing
    Q = c.scalingFactor
    # x is the index, y is the scale
    # y= aveLength + b*x^(1/p), solve for a and b with
    # (x₀,y₀)=(0,aveLength-1)
    # dy/dx = 1/s, so the Quality factor gives the slope at the last frequency
    b = (p/Q)^(1/p) * (O+a)^((p-1)/p)
    # the point x so that the second condition holds
    lastWavelet = Q * (O+a)/p
    # the point so that the first wavelet is at a
    firstWavelet = (a/b)^p
    # step size so that there are actually Q wavelets in the last octave
    startOfLastOctave = ((nOct+a-1)/b)^p
    stepSize = (lastWavelet - startOfLastOctave)/Q
    samplePoints = range(firstWavelet, stop=lastWavelet,
                         step = stepSize)#, length = round(Int, O * Q^(1/p)))
    return b .* (samplePoints).^(1/p)
end

function getNScales(n1, c)
    nOctaves = log2(max(n1, 2)) - c.averagingLength
    nWaveletsInOctave = reverse([max(1, round(Int,
                                              c.scalingFactor/x^(c.decreasing)))
                                 for x = 1:round(Int, nOctaves)])
    nScales = max(sum(nWaveletsInOctave), 0)
end

@doc """
    sj = getScales(n1::Integer, c::CFW)
since the spacing between wavelets isn't even, getting the actual scales is somewhat tricky
"""
function getScales(n1::Integer, c::CFW)
    nOctaves = log2(max(n1, 2)) - c.averagingLength
    nWaveletsInOctave = reverse([max(1, round(Int, c.scalingFactor /
                                              x^(c.decreasing))) for
                                 x=1:round(Int, nOctaves)])
    sj = zeros(sum(nWaveletsInOctave))
    for curOctave = 1:round(Int, nOctaves)
        sRange = (2 .^ (range(0, stop=1, length = nWaveletsInOctave[curOctave]+1) .+
                        curOctave .+ c.averagingLength .- 1))[1:end-1]
        ind = (sum(nWaveletsInOctave[1:curOctave-1])+1):sum(nWaveletsInOctave[1:curOctave])
        @debug "ind = $ind"
        sj[ind] = sRange
    end
    return sj
end

"""

adjust the length of the signal based on the boundary conditions
"""
function setn(n1, c)
    if boundaryType(c)() == WT.padded
        base2 = round(Int,log(n1 + 1)/log(2));   # power of 2 nearest to n1
        n = 2^(base2+1)
        n = n>>1 + 1
    elseif boundaryType(c)() == WT.DEFAULT_BOUNDARY
        # n1+1 rather than just n1 because this is going to be used in an rfft
        # for real data
        n = n1 + 1
    else
        n= n1>>1 + 1
    end
    return n
end

@doc """
      (daughters, ω) = computeWavelets(n1::Integer, c::CFW{W}; T=Float64, J1::Int64=-1, dt::S=NaN, s0::V=NaN) where {S<:Real,
                                                                   W<:WT.WaveletBoundary, V}
just precomputes the wavelets used by transform c::CFW{W}. For details, see cwt
"""
function computeWavelets(n1::Integer, c::CFW{W}; T=Float64, J1::Int64=-1, dt::S=NaN, s0::V=NaN) where {S<:Real,
                                                                   W<:WT.WaveletBoundary, V}
    # don't alter scaling with sampling information if it doesn't exists
    fλ = (4*π) / (c.σ[1] + sqrt(2 + c.σ[1]^2))
    if isnan(dt) || (dt<0)
        dt = 1
    end
    # smallest resolvable scale
    if isnan(s0) || (s0<0)
        s0 = 2 * dt / fλ
    end
    # J1 is the total number of scales
    # if J1<0
    #     J1 = Int(round(log2(n1 * dt / s0) * c.scalingFactor))
    # end

    nOctaves, totalWavelets, sRange, sWidth = getNWavelets(n1,c)

    # padding determines the actual number of elements
    n = setn(n1,c)
    # indicates whether we should keep a spot for the father wavelet
    isAve = (c.averagingLength > 0 && !(typeof(c.averagingType) <: NoAve)) ? 1 : 0
    # I guess matlab did occasionally do something useful

    ω, daughters = analyticOrNot(c, n, totalWavelets)

    # if the nOctaves is small enough there are none not covered by the
    # averaging, just use that
    if round(nOctaves) < 0
        father = findAveraging(c,ω, c.averagingType, sWidth[1])
        return father, ω
    end

    for (curWave, s) in enumerate(sRange)
        daughters[:,curWave+isAve] = Mother(c, s, sWidth[curWave], ω)
    end
    if isAve>0 # should we include the father?
        @debug "c = $(c), $(c.averagingType), size(ω)= $(size(ω))"
        daughters[:, 1] = findAveraging(c, ω, c.averagingType, sWidth[1])#[1:(n1+1)]
    end

    # adjust by the frame bound
    if c.frameBound > 0
        daughters = daughters.*(c.frameBound/norm(daughters, 2))
    end
    testFourierDomainProperties(daughters, isAve)
    return (daughters, ω)
end

# not
function analyticOrNot(c::CFW{W, T, <:Dog, N}, n, totalWavelets) where {W,T,N}
    ω = (0:(n-1))*2π
    # Odd derivatives introduce an imaginary term, so we need a complex representation
    if c.α % 2 == 1
        daughters = zeros(Complex{T}, n, totalWavelets)
    else
        daughters = zeros(T, n, totalWavelets)
    end
    return (ω, daughters)
end

function analyticOrNot(c::CFW{W, T, <:Union{Morlet, Paul}, N}, n, totalWavelets) where {W,T,N}
    ω = (0:(n-1))*2π
    daughters = zeros(T, n, totalWavelets)
    return (ω, daughters)
end

function testFourierDomainProperties(daughters, isAve)
    lowAprxAnalyt = daughters[1,2]./maximum(abs.(daughters[:,2]))
    if abs(lowAprxAnalyt) >= .01
        @warn "the lowest frequency wavelet has more than 1% its max at zero, so it may not be analytic. Think carefully" lowAprxAnalyt
    end
    highAprxAnalyt = abs.(daughters[end,end])./maximum(abs.(daughters[:,2]))
    if highAprxAnalyt >= .01
        @warn "the highest frequency wavelet has more than 1% its max at the end, so it may not be analytic. Think carefully" highAprxAnalyt
    end

    netWeight = sum(daughters, dims=2)
    centralFreqLast = argmax(abs.(daughters[:,end]))
    centralFreqFirst = argmax(abs.(daughters[:,1+isAve]))
    ratioOfCoverage = maximum(abs.(netWeight)) /
        minimum(abs.(netWeight[centralFreqFirst:centralFreqLast]))
    if ratioOfCoverage > 5
        @warn "there are frequencies which are significantly more covered than others, with a ratio of" ratioOfCoverage
    end
    minimalRegionComparedToLastPeak = maximum(abs.(daughters[:,end])) /
        minimum(abs.(netWeight[centralFreqFirst:centralFreqLast]))
    if minimalRegionComparedToLastPeak > 2
        @warn "there are wavelets whose peaks are far enough apart that the trough between them is less than half the height of the highest frequency wavelet" minimalRegionComparedToLastPeak
    end
end

# it's ok to just hand the total size, even if we're not transforming across
# all dimensions
function computeWavelets(Y::Tuple, c::CFW{W}; T=Float64) where {S<:Real,
                                                     W<:WT.WaveletBoundary}
    return computeWavelets(Y[1], c; T=T)
end
function computeWavelets(Y::AbstractArray{<:Integer}, c::CFW{W}; T=Float64) where {S<:Real,
                                                     W<:WT.WaveletBoundary}
    return computeWavelets(Y[1], c, T=T)
end

# also ok to just hand the whole thing being transformed
function computeWavelets(Y::AbstractArray{<:Number}, c::CFW{W}; T=Float64) where {S<:Real,
                                                     W<:WT.WaveletBoundary}
    return computeWavelets(size(Y)[1], c, T=T)
end

const CONT_DEFS = Dict{String,Tuple{String, String, String}}(
"mor" => ("(4*π)/(σ + sqrt(2 + σ.^2))", #FourierFactorFunction-- this will only work for σ≫1. Otherwise, you need a recursive algorithm to derive it.
      "1/sqrt(2)", #COIFunction
      "daughter = cσ*(π)^(1/4)*(exp.(-(σ-ω/s).^2/2)-κσ*exp.(-1/2*(ω/s).^2))", #This is based on the version in Wavelab, that corresponds to the Morlet wavelet on Wikipedia. It satisfies the admissibility condition exactly
      ),
"pau"=> ("4*π/(2*α+1)", #FourierFactorFunction
            "1*sqrt(2)", #COIFunction
            "daughter = ones(length(ω)); daughter[ω.<0]=0; daughter = daughter.*2α/sqrt(α*gamma(2*α)).*(ω/s).^α.*exp.(-(ω/s))", #MotherFunction
      ),
"dog" => ("2*π*sqrt(2 ./(2*α+1))", #FourierFactorFunction
      "1/sqrt(2)", #COIFunction
      "daughter = normalize(im^α*sqrt(gamma(α+1/2))*(ω/s).^α.*exp.(-(ω/s).^2/2))",
      )
)

# TRANSFORM TYPE CONSTRUCTORS

"""
    wavelet(c[, t=WT.Filter][, boundary=WT.Periodic][, s=Real])

Construct wavelet type where `c` is a wavelet class,
`t` is the transformation type (`WT.Filter` or `WT.Lifting`),
and `boundary` is the type of boundary treatment. In the continuous case, s>0
is the number of wavelets between the octaves ``2^J`` and ``2^{J+1}``
(defaults to 8, which is most appropriate for music and other audio).

# Examples
```julia
wavelet(WT.coif6)
wavelet(WT.db1, WT.Lifting)
wavelet(WT.morl,4)
```

**See also:** `WT.WaveletClass`
"""
function wavelet end

wavelet(c::K, boundary::WaveletBoundary=DEFAULT_BOUNDARY) where {K<:Union{BiOrthoWaveletClass, OrthoWaveletClass}}= wavelet(c, Filter, boundary)
wavelet(c::OrthoWaveletClass, t::FilterTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = OrthoFilter(c, boundary)
wavelet(c::WaveletClass, t::LiftingTransform, boundary::WaveletBoundary=DEFAULT_BOUNDARY) = GLS(c, boundary)

function wavelet(cw::T; s::S=8.0, boundary::WaveletBoundary=DEFAULT_BOUNDARY,
                 averagingType::A=Father(), averagingLength::Int = 4,
                 frameBound=1, normalization::N=Inf, decreasing=4) where {T<:ContinuousWaveletClass,
                                                                          A<:Average,
                                                                          S<:Real, N<:Real}
    if typeof(s) <: AbstractFloat
        decreasing = S(decreasing)
    end
    if typeof(decreasing) <: AbstractFloat
        s = typeof(decreasing)(s)
    end

    return CFW(cw,s, boundary, averagingType, averagingLength, frameBound,
               normalization, decreasing)
end

function wavelet(c::T, boundary::WaveletBoundary) where T<:WT.ContinuousWaveletClass
    CFW(c,8, boundary)
end
# ------------------------------------------------------------

# Compute filters from the Daubechies class
# N is the number of zeros at -1
function daubechies(N::Int)
    @assert N > 0
    # Create polynomial
    C = Vector{Int}(undef, N)
    @inbounds for n = 0:N-1
        C[N-n] = binomial(N-1+n, n)
    end

    # Find roots in y domain (truncated binomial series; (1 - y)^{-N})
    Y = roots(C)

    # Find roots in z domain:
    # z + z^{-1} = 2 - 4*y
    # where y is a root from above
    Z = zeros(ComplexF64, 2*N-2)
    @inbounds for i = 1:N-1
        Yi = Y[i]
        d = 2*sqrt( Yi*Yi - Yi )
        y2 = 1 - 2*Yi
        Z[i] = y2 + d
        Z[i+N-1] = y2 -d
    end

    # Retain roots inside unit circle
    nr = 0  # count roots
    @inbounds for i = eachindex(Z)
        if abs(Z[i]) <= 1 + eps()
            nr += 1
        end
    end

    # Find coefficients of the polynomial
    # (1 + z)^N * \prod_i (z - z_i)
    R = Vector{ComplexF64}(undef, N + nr)
    @inbounds for i = 1:N
        R[i] = -1
    end
    k = N
    @inbounds for i = eachindex(Z)
        if abs(Z[i]) <= 1 + eps()
            k += 1
            R[k] = Z[i]
        end
    end
    HH = vieta( R )

    # Normalize coefficients
    rmul!(HH, 1/norm(HH))
    return real(HH)
end

# Compute roots of polynomial
# Input is a coefficient vector with highest powers first
function roots(C::AbstractVector)
    A = compan(C)
    return eigvals(A)
end

# Create companion matrix for a polynomial
# Input is a coefficient vector with highest powers first
function compan(C::AbstractVector)
    n = length(C)
    A = zeros(n-1, n-1)

    if n > 1
        @inbounds A[1,:] .= -C[2:end] ./ C[1]
        @inbounds A[2:n:end] .= 1
    end
    return A
end

# Vieta-like formula for computing polynomial coefficients from roots
# See
# http://www.mathworks.se/help/matlab/ref/poly.html
function vieta(R::AbstractVector)
    n = length( R )
    C = zeros(ComplexF64, n+1)
    C[1] = 1
    Ci::ComplexF64 = 0
    Cig::ComplexF64 = 0

    @inbounds for k = 1:n
        Ci = C[1]
        for i = 1:k
            Cig = C[i+1]
            C[i+1] = Cig - R[k] * Ci
            Ci = Cig
        end
    end
    return C
end

# scaling filters h (low pass)
# the number at end of a filter name is the
# number of vanishing moments of the mother wavelet function
# sources:
# http://statweb.stanford.edu/~wavelab/ (Orthogonal/MakeONFilter.m)
# http://www.mathworks.com/matlabcentral/fileexchange/5502-filter-coefficients-to-popular-wavelets
### https://github.com/nigma/pywt/blob/master/src/wavelets_coeffs.template.h
# name => qmf
const FILTERS = Dict{String, Vector{Float64}}(
# Haar filter
"haar" =>
[0.7071067811865475,0.7071067811865475]
,
# Daubechies filters, see daubechies()

# Coiflet filters
"coif2" =>
[-0.072732619513,0.337897662458,0.852572020212,0.384864846864,-0.072732619513,-0.015655728135]
,
"coif4" =>
[0.0163873364635998,-0.0414649367819558,-0.0673725547222826,0.3861100668229939,0.8127236354493977,0.4170051844236707,-0.0764885990786692,-0.0594344186467388,0.0236801719464464,0.0056114348194211,-0.0018232088707116,-0.0007205494453679]
,
"coif6" =>
[-0.0037935128644910,0.0077825964273254,0.0234526961418362,-0.0657719112818552,-0.0611233900026726,0.4051769024096150,0.7937772226256169,0.4284834763776168,-0.0717998216193117,-0.0823019271068856,0.0345550275730615,0.0158805448636158,-0.0090079761366615,-0.0025745176887502,0.0011175187708906,0.0004662169601129,-0.0000709833031381,-0.0000345997728362]
,
"coif8" =>
[0.0008923136685824,-0.0016294920126020,-0.0073461663276432,0.0160689439647787,0.0266823001560570,-0.0812666996808907,-0.0560773133167630,0.4153084070304910,0.7822389309206135,0.4343860564915321,-0.0666274742634348,-0.0962204420340021,0.0393344271233433,0.0250822618448678,-0.0152117315279485,-0.0056582866866115,0.0037514361572790,0.0012665619292991,-0.0005890207562444,-0.0002599745524878,0.0000623390344610,0.0000312298758654,-0.0000032596802369,-0.0000017849850031]
,
"coif10" =>
[-0.0002120808398259,0.0003585896879330,0.0021782363583355,-0.0041593587818186,-0.0101311175209033,0.0234081567882734,0.0281680289738655,-0.0919200105692549,-0.0520431631816557,0.4215662067346898,0.7742896037334738,0.4379916262173834,-0.0620359639693546,-0.1055742087143175,0.0412892087544753,0.0326835742705106,-0.0197617789446276,-0.0091642311634348,0.0067641854487565,0.0024333732129107,-0.0016628637021860,-0.0006381313431115,0.0003022595818445,0.0001405411497166,-0.0000413404322768,-0.0000213150268122,0.0000037346551755,0.0000020637618516,-0.0000001674428858,-0.0000000951765727]
,
# Symmlet filter
"sym4" =>
[0.0455703458960000,-0.0178247014420000,-0.1403176241790000,0.4212345342040000,1.1366582434079999,0.7037390686560000,-0.0419109651250000,-0.1071489014180000]
,
"sym5" =>
[0.0276321529580000,-0.0298424998690000,-0.2479513626130000,0.0234789231360000,0.8965816483800000,1.0230529668940000,0.2819906968540000,-0.0553441861170000,0.0417468644220000,0.0386547959550000]
,
"sym6" =>
[-0.0110318675090000,0.0024999220930000,0.0632505626600000,-0.0297837512990000,-0.1027249698620000,0.4779043713330000,1.1138927839260000,0.6944579729580000,-0.0683231215870000,-0.1668632154120000,0.0049366123720000,0.0217847003270000]
,
"sym7" =>
[0.0145213947620000,0.0056713426860000,-0.1524638718960000,-0.1980567068070000,0.4081839397250000,1.0857827098140000,0.7581626019640000,0.0246656594890000,-0.0700782912220000,0.0960147679360000,0.0431554525820000,-0.0178704316510000,-0.0014812259150000,0.0037926585340000]
,
"sym8" =>
[-0.0047834585120000,-0.0007666908960000,0.0448236230420000,0.0107586117510000,-0.2026486552860000,-0.0866536154060000,0.6807453471900000,1.0991066305370001,0.5153986703740000,-0.0734625087610000,-0.0384935212630000,0.0694904659110000,0.0053863887540000,-0.0211456865280000,-0.0004283943000000,0.0026727933930000]
,
"sym9" =>
[0.0019811937360000,0.0008765025390000,-0.0187693968360000,-0.0163033512260000,0.0427444336020000,0.0008251409290000,-0.0771721610970000,0.3376589236020000,1.0152597908320000,0.8730484073490000,0.0498828309590000,-0.2708937835030000,-0.0257864459300000,0.0877912515540000,0.0125288962420000,-0.0145155785530000,-0.0006691415090000,0.0015124873090000]
,
"sym10" =>
[-0.0006495898960000,0.0000806612040000,0.0064957283750000,-0.0011375353140000,-0.0287862319260000,0.0081528167990000,0.0707035675500000,-0.0452407722180000,-0.0502565400920000,0.5428130112130000,1.0882515305000000,0.6670713381540000,-0.1002402150310000,-0.2255589722340000,0.0164188694260000,0.0649509245790000,-0.0020723639230000,-0.0122206426300000,0.0001352450200000,0.0010891704470000]
,
# Battle-Lemarie filter
"batt2" =>
[-0.0000867523000000,-0.0001586010000000,0.0003617810000000,0.0006529220000000,-0.0015570100000000,-0.0027458800000000,0.0070644200000000,0.0120030000000000,-0.0367309000000000,-0.0488618000000000,0.2809310000000000,0.5781630000000000,0.2809310000000000,-0.0488618000000000,-0.0367309000000000,0.0120030000000000,0.0070644200000000,-0.0027458800000000,-0.0015570100000000,0.0006529220000000,0.0003617810000000,-0.0001586010000000,-0.0000867523000000]
,
"batt4" =>
[0.0001033070000000,-0.0001642640000000,-0.0002018180000000,0.0003267490000000,0.0003959460000000,-0.0006556200000000,-0.0007804680000000,0.0013308600000000,0.0015462400000000,-0.0027452900000000,-0.0030786300000000,0.0057993200000000,0.0061414300000000,-0.0127154000000000,-0.0121455000000000,0.0297468000000000,0.0226846000000000,-0.0778079000000000,-0.0354980000000000,0.3068300000000000,0.5417360000000000,0.3068300000000000,-0.0354980000000000,-0.0778079000000000,0.0226846000000000,0.0297468000000000,-0.0121455000000000,-0.0127154000000000,0.0061414300000000,0.0057993200000000,-0.0030786300000000,-0.0027452900000000,0.0015462400000000,0.0013308600000000,-0.0007804680000000,-0.0006556200000000,0.0003959460000000,0.0003267490000000,-0.0002018180000000,-0.0001642640000000,0.0001033070000000]
,
"batt6" =>
[0.0001011130000000,0.0001107090000000,-0.0001591680000000,-0.0001726850000000,0.0002514190000000,0.0002698420000000,-0.0003987590000000,-0.0004224850000000,0.0006355630000000,0.0006628360000000,-0.0010191200000000,-0.0010420700000000,0.0016465900000000,0.0016413200000000,-0.0026864600000000,-0.0025881600000000,0.0044400200000000,0.0040788200000000,-0.0074684800000000,-0.0063988600000000,0.0128754000000000,0.0099063500000000,-0.0229951000000000,-0.0148537000000000,0.0433544000000000,0.0208414000000000,-0.0914068000000000,-0.0261771000000000,0.3128690000000000,0.5283740000000000,0.3128690000000000,-0.0261771000000000,-0.0914068000000000,0.0208414000000000,0.0433544000000000,-0.0148537000000000,-0.0229951000000000,0.0099063500000000,0.0128754000000000,-0.0063988600000000,-0.0074684800000000,0.0040788200000000,0.0044400200000000,-0.0025881600000000,-0.0026864600000000,0.0016413200000000,0.0016465900000000,-0.0010420700000000,-0.0010191200000000,0.0006628360000000,0.0006355630000000,-0.0004224850000000,-0.0003987590000000,0.0002698420000000,0.0002514190000000,-0.0001726850000000,-0.0001591680000000,0.0001107090000000,0.0001011130000000]
,
# Beylkin filter
"beyl" =>
[0.0993057653740000,0.4242153608130000,0.6998252140570000,0.4497182511490000,-0.1109275983480000,-0.2644972314460000,0.0269003088040000,0.1555387318770000,-0.0175207462670000,-0.0885436306230000,0.0196798660440000,0.0429163872740000,-0.0174604086960000,-0.0143658079690000,0.0100404118450000,0.0014842347820000,-0.0027360316260000,0.0006404853290000]
,
# Vaidyanathan filter
"vaid" =>
[-0.0000629061180000,0.0003436319050000,-0.0004539566200000,-0.0009448971360000,0.0028438345470000,0.0007081375040000,-0.0088391034090000,0.0031538470560000,0.0196872150100000,-0.0148534480050000,-0.0354703986070000,0.0387426192930000,0.0558925236910000,-0.0777097509020000,-0.0839288843660000,0.1319716614170000,0.1350842271290000,-0.1944504717660000,-0.2634948024880000,0.2016121617750000,0.6356010598720000,0.5727977932110000,0.2501841295050000,0.0457993341110000]

)

# biortho filters
# name => (qmf1,qmf2)
const BIFILTERS = Dict{String, NTuple{2, Vector{Float64}}}(
# test
"test" =>
([0.7071067811865475,0.7071067811865475],
 [0.7071067811865475,0.7071067811865475])

)

# in matlab (d,p) -> (predict, update)
const SCHEMES = Dict{String,NTuple{3, Any}}(
# cdf 5/3 -> bior 2.2, cdf 9/7 -> bior 4.4
# Cohen-Daubechies-Feauveau [Do Quan & Yo-Sung Ho. Optimized median lifting scheme for lossy image compression.]
"cdf9/7" => ([  LSStep(Update,  [1.0,1.0]*1.5861343420604, 0),
                LSStep(Predict, [1.0,1.0]*0.05298011857291494, 1),
                LSStep(Update,  [1.0,1.0]*-0.882911075531393, 0),
                LSStep(Predict, [1.0,1.0]*-0.44350685204384654, 1)],
                1.1496043988603355,
                0.8698644516247099)
,

# Haar, same as db1
"haar" => ([LSStep(Predict, [-1.0], 0),
            LSStep(Update,  [0.5], 0)],
            0.7071067811865475,
            1.4142135623730951)
,
# Daubechies
"db1" => ([ LSStep(Predict, [-1.0], 0),
            LSStep(Update,  [0.5], 0)],
            0.7071067811865475,
            1.4142135623730951)
,
"db2" => ([ LSStep(Predict, [-1.7320508075688772], 0),
            LSStep(Update,  [-0.0669872981077807,0.4330127018922193], 1),
            LSStep(Predict, [1.0], -1)],
            0.5176380902050414,
            1.9318516525781364)

)

end


text:
#!/usr/bin/env julia

using BinaryBuilder, LibGit2, SHA, Pkg
using BinaryBuilder.Wizard: yn_prompt, WizardState

function de_jll(name::String)
    if endswith(name, "_jll")
        name = name[1:end-4]
    end
    return name
end

function build_tarballs_path(dep_name::String, base::String = @__DIR__)
    return joinpath(base, uppercase(dep_name[1:1]), de_jll(dep_name), "build_tarballs.jl")
end

if !isfile(build_tarballs_path(get(Sys.ARGS, 1, "")))
    println(stderr, "usage: recursively_regenerate_jlls.jl <jll name>")
    exit(1)
end

json_objs = Dict()
function get_json_obj(dep_name::String)
    dep_name = de_jll(dep_name)
    if haskey(json_objs, dep_name)
        return json_objs[dep_name]
    end
    build_tarballs = build_tarballs_path(dep_name)
    if !isfile(build_tarballs)
        @warn("Unable to find build_tarballs.jl file for dependency \"$(dep["name"])\"")
    end
    @info("Parsing $(build_tarballs)")

    # We're going to `include()` the "build_tarballs.jl" function to avoid the overhead of `using BinaryBuilder`
    # every time we launch a new julia process.
    meta_json = tempname()
    m = Module(:__anon__)
    Core.eval(m, quote
        using BinaryBuilder

        # Our special overriding build_tarballs() function:
        function build_tarballs_meta_json(cli_args, args...; kwargs...)
            cli_args = vcat(string("--meta-json=", $(meta_json)), cli_args)
            BinaryBuilder.build_tarballs(cli_args, args...; kwargs...)
        end

        # Clear out ARGS from parent module.  :P
        empty!(ARGS)
    end)
    # Be super-sneaky and replace `build_tarballs()` funciton calls with `build_tarballs_meta_json()`
    Base.include(m, build_tarballs) do expr
        if hasproperty(expr, :head) && expr.head == :call && expr.args[1] == :build_tarballs
            expr.args[1] = :build_tarballs_meta_json
        end
        return expr
    end

    json_objs[dep_name] = BinaryBuilder.JSON.parse(String(read(meta_json)))
    BinaryBuilder.cleanup_merged_object!(json_objs[dep_name])
    rm(meta_json; force=true)
    return json_objs[dep_name]
end

function recursively_collect_dependencies(dep_name::String, dependencies = Set{String}())
    # Generate meta.json for this script, read it in
    meta_obj = get_json_obj(dep_name)
    for dep in meta_obj["dependencies"]
        name = BinaryBuilder.BinaryBuilderBase.getname(dep)
        if !(name ∈ dependencies)
            push!(dependencies, name)
            recursively_collect_dependencies(name, dependencies)
        end
    end
    return dependencies
end

function download_binaries_from_release(dep_name, code_dir, download_dir)
    Pkg.PlatformEngines.probe_platform_engines!()

    # Doownload the tarballs reading the information in the current `Artifacts.toml`.
    artifacts = Pkg.Artifacts.load_artifacts_toml(joinpath(code_dir, "Artifacts.toml"))
    artifact_listing = artifacts[dep_name]

    # Deal with non-platform-specific artifacts
    if isa(artifact_listing, Dict)
        artifact_listing = [artifact_listing]
    end

    for artifact in artifact_listing
        info = artifact["download"][1]
        url = info["url"]
        hash = info["sha256"]
        filename = basename(url)
        Pkg.PlatformEngines.download_verify(url, hash, joinpath(download_dir, filename); verbose=true)
    end
end

function regenerate_jll(dep_name::String)
    dep_name = de_jll(dep_name)
    json_obj = get_json_obj(dep_name)
    code_dir = joinpath(Pkg.devdir(), "$(dep_name)_jll")
    repo = "JuliaBinaryWrappers/$(dep_name)_jll.jl"
    BinaryBuilder.init_jll_package(
        dep_name,
        code_dir,
        repo,
    )

    download_dir = joinpath("/tmp/yggdrasil_download_temp", dep_name)
    mkpath(download_dir)
    download_binaries_from_release(dep_name, code_dir, download_dir)
    build_version = BinaryBuilder.get_next_wrapper_version(json_obj["name"], json_obj["version"])
    tag = "$(json_obj["name"])-v$(build_version)"
    upload_prefix = "https://github.com/$(repo)/releases/download/$(tag)"
    BinaryBuilder.rebuild_jll_package(
        json_obj;
        download_dir=download_dir,
        upload_prefix=upload_prefix,
        verbose=true,
        lazy_artifacts=json_obj["lazy_artifacts"],
        from_scratch=true,
    )
end

# This borrows heavily from BinaryBuilder.yggdrasil_deploy()
function open_jll_bump_pr(dep_name::String)
    @info("Opening JLL-bumping PR for $(dep_name)")
    gh_auth = BinaryBuilder.Wizard.github_auth(;allow_anonymous=false)
    fork = BinaryBuilder.GitHub.create_fork("JuliaPackaging/Yggdrasil"; auth=gh_auth)

    mktempdir() do tmp
        repo = LibGit2.clone(BinaryBuilder.Wizard.get_yggdrasil(), tmp)
        build_tarballs = build_tarballs_path(dep_name, tmp)
        if !isfile(build_tarballs)
            throw(ArgumentError("Invalid dep_name \"$(dep_name)\", no file $(build_tarballs)"))
        end

        # We're going to need to alter the build_tarballs.jl file.
        # We will just add a newline onto the end if it doesn't already
        # have two, otherwise we remove one.
        build_tarballs_content = String(read(build_tarballs))
        if build_tarballs_content[end-1:end] != "\n\n"
            build_tarballs_content = string(build_tarballs_content, "\n")
        else
            build_tarballs_content = build_tarballs_content[1:end-1]
        end

        recipe_hash = bytes2hex(sha256(build_tarballs_content))[end-3:end]
        branch_name = "jllbump/$(dep_name)_$(recipe_hash)"
        LibGit2.branch!(repo, branch_name)

        # Write out the modified build_tarballs.jl file
        open(build_tarballs, "w") do io
            write(build_tarballs, build_tarballs_content)
        end

        # Commit modified build_tarballs.jl file
        LibGit2.add!(repo, build_tarballs_path(dep_name, ""))
        LibGit2.commit(repo, "Bump JLL wrapper: $(dep_name)\n\n[skip build]")
        creds = LibGit2.UserPasswordCredential(
            dirname(fork.full_name),
            deepcopy(gh_auth.token),
        )
        try
            LibGit2.push(
                repo,
                refspecs=["+HEAD:refs/heads/$(branch_name)"],
                remoteurl="https://github.com/$(fork.full_name).git",
                credentials=creds,
                # This doesn't work :rage: instead we use `+HEAD:` at the beginning of our
                # refspec: https://github.com/JuliaLang/julia/issues/23057
                #force=true,
            )
        finally
            Base.shred!(creds)
        end

		params = Dict(
			"base" => "master",
			"head" => "$(dirname(fork.full_name)):$(branch_name)",
			"maintainer_can_modify" => true,
			"title" => "JLL bump: $(dep_name)",
			"body" => """
			This pull request bumps the JLL version of $(dep_name).
			It was generated via the `recursively_regenerate_jlls.jl` script.
			"""
		)
		pr = BinaryBuilder.create_or_update_pull_request("JuliaPackaging/Yggdrasil", params, auth=gh_auth)
        println("https://github.com/$(fork.full_name)/pull/new/$(BinaryBuilder.HTTP.escapeuri(branch_name))?expand=1")
    end
end

toplevel_dep_name = Sys.ARGS[1]
deps = recursively_collect_dependencies(toplevel_dep_name)
push!(deps, toplevel_dep_name)
println("Discovered dependencies: $(collect(deps))")

if yn_prompt(WizardState(), "Open JLL-bumping PRs?", :y) == :y
	for dep in deps
		open_jll_bump_pr(dep)
	end
elseif yn_prompt(WizardState(), "Generate new JLLs locally?", :y) == :y
    for dep in deps
        regenerate_jll(dep)
    end
end


text:
using MeshSlicer
using PyPlot

println("First mesh load")
@time mesh = PolygonMesh("cell_lamp.stl");
println("First rotation")
@time rotate!(mesh, 45.0, [1.0,0.0,0.0],[0.0,0.0,0.0])
println("Second mesh load")
@time mesh = PolygonMesh("cell_lamp.stl");
println("Second rotation")
@time rotate!(mesh, 45.0, [1.0,0.0,0.0],[0.0,0.0,0.0])
layers = 100;

# find slice heights
step = (mesh.bounds.zmax - mesh.bounds.zmin)/layers

heights = [mesh.bounds.zmin:step:mesh.bounds.zmax]
# create an array of slices
println("First Slice")
@time a = MeshSlice(mesh,heights)
println("Second Slice")
@time b = MeshSlice(mesh,heights)

for slice in b
    print(length(slice.polygons), " ")
    for poly in slice.polygons
        for line in poly.segments
            plot3D([line.start[1],line.finish[1]], [line.start[2],line.finish[2]], slice.layer)
        end
   end
end


text:
using StarWarsArrays
using Test

import StarWarsArrays: StarWarsError, order

@testset "OriginalOrder" begin
    v = StarWarsArray(collect(1:9))
    @test order(v) == OriginalOrder
    @test size(v) == size(parent(v))
    @test v[5] == 2
    @test v[3] == 6
    @test v[7] == 7
    @test_throws BoundsError v[0]
    v[4] = -42
    @test v[4] == -42

    m = StarWarsArray(reshape(collect(1:81), 9, 9), OriginalOrder)
    @test order(m) == OriginalOrder
    @test size(m) == size(parent(m))
    @test m[5] == 2
    @test m[3,4] == 6
    @test m[7,2] == 43
    @test m[4,4] == 1
    @test_throws BoundsError m[0, 12]
    m[4,4] = -42
    @test m[4] == -42
    m[1] = 0
    @test m[1,4] == 0
end

@testset "MacheteOrder" begin
    v = StarWarsArray(collect(1:9), MacheteOrder)
    @test order(v) == MacheteOrder
    @test size(v) == (8,)
    @test v[5] == 2
    @test v[3] == 4
    @test v[6] == 5
    @test v[7] == 6
    @test_throws StarWarsError v[1]
    @test_throws BoundsError v[0]
    v[4] = -42
    @test v[4] == -42

    m = StarWarsArray(reshape(collect(1:81), 9, 9), MacheteOrder)
    @test order(m) == MacheteOrder
    @test size(m) == (8, 8)
    @test m[5] == 2
    @test m[3,4] == 4
    @test m[7,2] == 24
    @test m[4,4] == 1
    @test m[9,4] == 8
    @test_throws BoundsError m[0, 12]
    @test_throws BoundsError m[82]
    @test_throws StarWarsError m[1]
    @test_throws StarWarsError m[3, 1]
    rpr = split(repr("text/plain", m), '\n')
    @test rpr[2] == " 21  30  3  12  39  48  57  66"
    @test rpr[6] == " 23  32  5  14  41  50  59  68"
    m[4,4] = -42
    @test m[4] == -42
    m[2] = 0
    @test m[2,4] == 0

    @test sprint(showerror, StarWarsArrays.StarWarsError(1, MacheteOrder)) ==
        "StarWarsError: there is no episode 1 in MacheteOrder"
end


text:
using OrdinaryDiffEq, Test, LinearAlgebra, Statistics

# TODO: clean up
@testset "Mass Matrix Accuracy Tests" begin
  function make_prob(mm_A; iip)
    mm_b = mm_A*ones(3)
    # iip
    function mm_f_iip(du,u,p,t)
      mul!(du,mm_A,u)
      tmp = t*mm_b
      du .+= tmp
    end
    mm_g_iip(du,u,p,t) = @. du = u + t
    # oop
    function mm_f_oop(u,p,t)
      du = u./t
      mm_f_iip(du,u,p,t)
      du
    end
    function mm_g_oop(u,p,t)
      du = u./t
      mm_g_iip(du,u,p,t)
      du
    end
    mm_analytic(u0,p,t) = @. 2u0*exp(t) - t - 1
    f = ((mm_f_oop, mm_g_oop), (mm_f_iip, mm_g_iip))[iip+1]
    prob = ODEProblem(ODEFunction(f[1],analytic=mm_analytic,mass_matrix=mm_A),ones(3),
                                     (0.0,1.0))
    prob2 = ODEProblem(ODEFunction(f[2],analytic=mm_analytic),ones(3),(0.0,1.0))
    return prob, prob2
  end
  for iip in (false, true)
    mm_A = [-2.0 1 4
             4 -2 1
             2 1 3]
    prob, prob2 = make_prob(mm_A; iip=iip)

    ######################################### Test each method for exactness

    if iip
      sol = solve(prob  ,Rosenbrock23())
      sol2 = solve(prob2,Rosenbrock23())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob, Rosenbrock32())
      sol2 = solve(prob2,Rosenbrock32())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob,  ROS3P())
      sol2 = solve(prob2,ROS3P())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob,  Rodas3())
      sol2 = solve(prob2,Rodas3())

      @test norm(sol .- sol2) ≈ 0 atol=1e-11

      sol = solve(prob,  RosShamp4())
      sol2 = solve(prob2,RosShamp4())

      @test norm(sol .- sol2) ≈ 0 atol=1e-10

      sol = solve(prob,  Rodas4())
      sol2 = solve(prob2,Rodas4())

      @test norm(sol .- sol2) ≈ 0 atol=1e-9

      sol = solve(prob,  Rodas5())
      sol2 = solve(prob2,Rodas5())

      @test norm(sol .- sol2) ≈ 0 atol=1e-7
    end

    sol = solve(prob,  ImplicitEuler(),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(),dt=1/10,adaptive=false)

    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob,  RadauIIA5(),dt=1/10,adaptive=false)
    sol2 = solve(prob2,RadauIIA5(),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-12

    sol = solve(prob,  ImplicitMidpoint(extrapolant = :constant),dt=1/10)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant),dt=1/10)

    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    # Anderson method
    sol = solve(prob,  ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)

    @test_broken norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob,  ImplicitMidpoint(nlsolve=NLAnderson(),extrapolant = :constant),dt=1/10)
    sol2 = solve(prob2,ImplicitMidpoint(nlsolve=NLAnderson(),extrapolant = :constant),dt=1/10)

    @test_broken norm(sol .- sol2) ≈ 0 atol=1e-7
    # Functional iteration
    prob, prob2 = make_prob(Matrix{Float64}(1.01I, 3, 3); iip=iip)
    sol = solve(prob,ImplicitEuler(
                          nlsolve=NLFunctional(κ=2000.,tol=1e-7,min_iter=10,max_iter=100)),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant,
                          nlsolve=NLFunctional(κ=2000.,tol=1e-7,min_iter=10,max_iter=100)),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant),dt=1/10)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7

    sol = solve(prob,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitEuler(nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7

    sol = solve(prob, ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    sol2 = solve(prob2,ImplicitMidpoint(extrapolant = :constant, nlsolve=NLAnderson()),dt=1/10,adaptive=false)
    @test norm(sol .- sol2) ≈ 0 atol=1e-7
    @test norm(sol[end] .- sol2[end]) ≈ 0 atol=1e-7
  end
end

#=

sol = solve(prob,  Trapezoid())
sol2 = solve(prob2,Trapezoid())

@test norm(sol .- sol2) ≈ 0 atol=1e-7

sol = solve(prob,  TRBDF2())
sol2 = solve(prob2,TRBDF2())

@test norm(sol .- sol2) ≈ 0 atol=1e-9

#sol = solve(prob, SDIRK2())
#sol2 = solve(prob2, SDIRK2())

#@test norm(sol .- sol2) ≈ 0 atol=1e-11

sol = solve(prob,   TRBDF2(),adaptive=false,dt=1/10)
sol2 = solve(prob2, TRBDF2(),adaptive=false,dt=1/10)

=#

# Singular mass matrices

@testset "Mass Matrix Tests with Singular Mass Matrices" begin
  function f!(du, u, p, t)
      du[1] = u[2]
      du[2] = u[2] - 1.
      return
  end

  u0 = [0.,1.]
  tspan = (0.0, 1.0)

  M = fill(0., 2,2)
  M[1,1] = 1.

  m_ode_prob = ODEProblem(ODEFunction(f!;mass_matrix=M), u0, tspan)
  @test_nowarn sol = solve(m_ode_prob, Rosenbrock23())

  M = [0.637947  0.637947
       0.637947  0.637947]

  inv(M) # not caught as singular

  function f2!(du, u, p, t)
      du[1] = u[2]
      du[2] = u[1]
      return
  end
  u0 = fill(0., 2)

  m_ode_prob = ODEProblem(ODEFunction(f2!;mass_matrix=M), u0, tspan)
  @test_nowarn sol = solve(m_ode_prob, Rosenbrock23())
end


text:
using NewtonsMethod
using Test
using LinearAlgebra

@testset "NewtonTests" begin

# Testing first function

g(x)=(x-1)^2
g′(x)=2*(x-1)
@test newtonroot(g,g′;x₀=2) ≈ 1.0 atol=1E-4
@test newtonroot(g;x₀=2)  ≈ 1.0 atol=1E-4

#Testing second function

h(x)=exp(x)-1
h′(x)=exp(x)
@test newtonroot(h,h′;x₀=3) ≈ 0.0 atol=1E-4
@test newtonroot(h;x₀=3) ≈ 0.0  atol=1E-4

#Testing function without convergence

j(x)=x^2+2
j′(x)=2*x

@test newtonroot(j,j′;x₀=1) == nothing
@test newtonroot(j;x₀=1) == nothing

#Testing maxiter

@test newtonroot(h,h′;x₀=3,maxiter=5) == nothing

#Testing tolerance
#From above
@test newtonroot(h,h′;x₀=3,tol=1E10) == 3

#From below
root1=newtonroot(g,g′;x₀=3,tol=1E-6)
root2=newtonroot(g,g′;x₀=3)
@test norm(root1-1)>norm(root2-1)

#Testing BigFloat

@test newtonroot(h,h′;x₀=BigFloat(3),tol=BigFloat(1E-7),maxiter=BigFloat(1000)) ≈ 0.0 atol=1E-4

l(x)=log(x)

@test_broken newtonroot(l;x₀=3) ≈ 1.0 atol=1E-4

end


text:
module TestGetObj

using Aqua: getobj
using Test

module ModuleA
function f end
end

module ModuleB
using ..ModuleA: ModuleA
ModuleA.f(::Int) = nothing
end

@testset begin
    m, = methods(ModuleA.f, Tuple{Int})
    @test getobj(m) === ModuleA.f
end

end  # module


text:
using FieldMetadata, Parameters, Test, Markdown, REPL

abstract type AbstractTest end

import FieldMetadata: @description, description, MetadataError

@description mutable struct Described{P}
   a::Int | "an Int"
   b      | "an untyped field"
   c::P   | "a parametric field"
   "An inner constructor should work, even with a docstring"
   Described(a, b, c) = begin
       new{typeof(c)}(a, b, c)
   end
   function Described(a, b)
       new{Nothing}(a, b, nothing)
   end
end

@test length(methods(Described).ms) == 2
@test length(methods(description).ms) == 12

d = Described(1, 1.0, nothing)
@test typeof(d) == typeof(Described(1, 1.0))

@test description(d, :a) == "an Int"
@test description(Described, Val{:a}) == "an Int"
@test description(Described, :a) == "an Int"
@test description(d, Val{:a}) == "an Int"
@test description(typeof(d), :b) == "an untyped field"
@test description(typeof(d), :c) == "a parametric field"
@test description(d) == ("an Int", "an untyped field", "a parametric field")

@test description(d, :d) == ""
@test description(d, Val{:d}) == ""
@test description(Described, Val{:d}) == ""
@test description(Described, :d) == ""

@inferred description(d, :a)
@inferred description(Described, Val{:a})
@inferred description(Described, :a)
@inferred description(d, Val{:a})
@inferred description(typeof(d), :b)
@inferred description(d, :c)
@inferred description(d, Val{:c})
@inferred description(Described, Val{:c})
@inferred description(Described, :c)
@inferred description(d)

ex = :(@some @arbitrary @macros struct TestMacros{T}
        a::T | u"1"
        b::T | u"2"
    end)
@test FieldMetadata.chained_macros(ex) == [Symbol("@some"), Symbol("@arbitrary"), Symbol("@macros")]

@metadata paramrange [0, 1] Array

# range array
@paramrange struct WithRange <: AbstractTest
    a::Int | [1, 4]
    b::Int | [4, 9]
end

w = WithRange(2,5)
@test paramrange(w, :a) == [1, 4]
@test paramrange(w, :b) == [4, 9]
@test paramrange(w) == ([1, 4], [4, 9])
@inferred paramrange(w, :a)

@paramrange @description struct WrongType <: AbstractTest
    a::Int | (1, 4) | :test
    b::Int | "4, 9" | AbstractTest
end

wt = WrongType(0, 0)
@test_throws MetadataError paramrange(wt, :a)
@test_throws MetadataError paramrange(WrongType, :b)
@test_throws MetadataError description(wt, :a)
@test_throws MetadataError description(WrongType, :b)

# combinations of metadata
@paramrange @description struct Combined{T} <: AbstractTest
    a::T | [1, 4] | "an Int with a range and a description"
    b::T | _      | "a Float with a range and a description"
end

c = Combined(3,5)
@test description(typeof(c), :a) == "an Int with a range and a description"
@test description(c, :b) == "a Float with a range and a description"
@test paramrange(typeof(c), :a) == [1, 4]
@test paramrange(c, :b) == [0, 1]
description(c, Val{:a})

# with Parameters.jl keywords
@paramrange @description @with_kw struct Keyword{T} <: AbstractTest
    a::T = 3 | [0, 100] | "an Int with a range and a description"
    b::T = 5 | [2, 9]   | "a Float with a range and a description"
end

k = Keyword()
@test paramrange(k, :a) == [0, 100]
@test description(k, :b) == "a Float with a range and a description"
@test k.a == 3
@test k.b == 5

# with missing keywords
@paramrange @description @with_kw struct MissingKeyword{T} <: AbstractTest
    a    = 3 | [0, 100] | "an Int with a range and a description"
    b::T     | [2, 9]   | "a Float with a range and a description"
    MissingKeyword{T}(a::T, b::T) where T = new{T}(a, b)
end

m = MissingKeyword(b = 99)
@test paramrange(m, :a) == [0, 100]
@test description(m, :b) == "a Float with a range and a description"
@test m.a == 3
@test m.b == 99
@test paramrange(m) == ([0, 100], [2, 9])

@metadata default nothing

@default @description @paramrange Described begin
    a | 1 | "a description updated in a begin block" | [22,33]
    b | 2 | "another updated description"            | [-8,-9]
end

@test paramrange(d, :a) == [22,33]
@test paramrange(d, :b) == [-8,-9]
@test description(d, :a) == "a description updated in a begin block"
@test description(d, :b) == "another updated description"

# Now eval in all the values
def_a = 3
def_b = 4
desc_a = "interpolated string a"
desc_b = "interpolated string b"
vec_a = [0, 0]
vec_b = [1, 1]

@eval   @default   @description   @paramrange :($(typeof(d))) begin
    a | $def_a   | $desc_a      | $vec_a
    b | $def_b   | $desc_b      | $vec_b
end

@test default(d, :a) == 3
@test default(d, :b) == 4
@test paramrange(d, :a) == [0, 0]
@test paramrange(d, :b) == [1, 1]
@test description(d, :a) == "interpolated string a"
@test description(d, :b) == "interpolated string b"

def_a = 44
@test default(d, :a) == 3

# docstrings
"The Docs"
@paramrange mutable struct Documented
    "Foo"
    a::Int     | [1,2]
    "Bar"
    b::Float64 | [3,4]
end

@test paramrange(Documented) == ([1,2], [3,4])

@test "The Docs\n" == Markdown.plain(REPL.doc(Documented))
@test "Foo\n" == Markdown.plain(REPL.fielddoc(Documented, :a))
@test "Bar\n" == Markdown.plain(REPL.fielddoc(Documented, :b))

# chaining macros

@chain columns @description @paramrange

@columns Described begin
    a::T | "a new Int description"     | [99,100]
    b::T | "a new Float64 description" | [-3,-4]
end

@test description(Described) == ("a new Int description", "a new Float64 description", "a parametric field")


text:
"""
This function adds the power limits of generators when there are no CommitmentVariables
"""
function activepower(m::JuMP.AbstractModel, devices::Array{R,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {R <: PowerSystems.RenewableGen, D <: AbstractRenewableDispatchForm, S <: PM.AbstractPowerFormulation}

    p_re = m[:p_re]
    time_index = m[:p_re].axes[2]
    name_index = m[:p_re].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    pmax_re = JuMP.JuMPArray(Array{ConstraintRef}(undef, length.(JuMP.axes(p_re))), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)

        if name == devices[ix].name
            pmax_re[name, t] = @constraint(m, p_re[name, t] <= devices[ix].tech.installedcapacity * values(devices[ix].scalingfactor)[t])
        else
            @error "Bus name in Array and variable do not match"
        end

    end

    JuMP.register_object(m, :pmax_re, pmax_re)

end

"""
This function adds the power limits of generators when there are no CommitmentVariables
"""
function reactivepower(m::JuMP.AbstractModel, devices::Array{R,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {R <: PowerSystems.RenewableGen, D <: AbstractRenewableDispatchForm,  S <: AbstractACPowerModel}

    q_re = m[:q_re]
    time_index = m[:q_re].axes[2]
    name_index = m[:q_re].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    qmax_re = JuMP.JuMPArray(Array{ConstraintRef}(undef,length.(JuMP.axes(q_re))), name_index, time_index)
    qmin_re = JuMP.JuMPArray(Array{ConstraintRef}(undef,length.(JuMP.axes(q_re))), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)

        if name == devices[ix].name

            qmax_re[name, t] = @constraint(m, q_re[name, t] <= devices[ix].tech.reactivepowerlimits.max)
            qmin_re[name, t] = @constraint(m, q_re[name, t] >= devices[ix].tech.reactivepowerlimits.min)

        else
            @error "Bus name in Array and variable do not match"
        end

    end

    JuMP.register_object(m, :qmax_re, qmax_re)
    JuMP.register_object(m, :qmin_re, qmin_re)

end


text:
##################################################
# Search index structures and associated methods #
##################################################
"""
HNSW index type for storing vectors. It is a wrapper around a
`HierarchicalNSW` (Hierarchical Navigable Small Worlds) NN graph
structure and performs a very efficient search using a distance-based
similarity between vectors.

# References
 * [Y. A. Malkov, D.A. Yashunin "Efficient and robust approximate nearest neighbor search using Hierarchical Navigable Small World graphs"](https://arxiv.org/abs/1603.09320)
"""
struct HNSWIndex{I,E,A,D} <: AbstractIndex
    tree::HierarchicalNSW{I,E,A,D}
end

HNSWIndex(data, args...; kwargs...) =
    HNSWIndex(
        add_to_graph!(
            HierarchicalNSW(
                [densify(data[:, i]) for i in 1:size(data, 2)];
                kwargs...)
           )
       )  # args are ignored

# Nearest neighbor search method
function knn_search(index::HNSWIndex{I,E,A,D},
                    point::AbstractVector,
                    k::Int,
                    keep::Vector{Int}=collect(1:length(index));
                    kwargs...
                   ) where {I<:Unsigned, E<:Real, A<:AbstractArray, D<:Metric}
    # Uses Euclidean distance by default
    _idxs, scores = knn_search(index.tree, Vector(point), k)
    idxs = Int.(_idxs)
    if length(keep) == length(index)
        # all data points are valid
        return idxs, scores
    else
        # this bit is slow if 'keep' is large
        mask = map(idx->in(idx, keep), idxs)
        return idxs[mask], scores[mask]
    end
end

# Length method
length(index::HNSWIndex) = length(index.tree.data)


text:
abstract type SCGALResults end
abstract type OptimizationLog end

struct SCGALSolution{T<:Real} <: SCGALResults
    UT::AbstractMatrix{T}
    ΛT::AbstractMatrix{T}
    ST::AbstractMatrix{T}
    Ω::AbstractMatrix{T}
    yT::AbstractVector{T}
    zT::AbstractVector{T}
    log::Union{OptimizationLog, Nothing}
end

struct SCGALLog{T<:Real} <: OptimizationLog
    dual_gap::AbstractVector{T}
    obj_val::AbstractVector{T}
    primal_infeas::AbstractVector{T}
    time_sec::AbstractVector{T}
    obj_val_Xhat::Union{AbstractVector{T}, Nothing}
    primal_infeas_Xhat::Union{AbstractVector{T}, Nothing}
    cut_value::Union{AbstractVector{T}, Nothing}
end


text:
using RealDot
using LinearAlgebra
using Test

# struct need to be defined outside of tests for julia 1.0 compat
# custom complex number (tests fallback definition)
struct CustomComplex{T} <: Number
    re::T
    im::T
end

Base.real(x::CustomComplex) = x.re
Base.imag(x::CustomComplex) = x.im

Base.conj(x::CustomComplex) = CustomComplex(x.re, -x.im)

function Base.:*(x::CustomComplex, y::Union{Real,Complex})
    return CustomComplex(reim(Complex(reim(x)...) * y)...)
end
Base.:*(x::Union{Real,Complex}, y::CustomComplex) = y * x
function Base.:*(x::CustomComplex, y::CustomComplex)
    return CustomComplex(reim(Complex(reim(x)...) * Complex(reim(y)...))...)
end

# custom quaternion to test definition for hypercomplex numbers
# adapted from Quaternions.jl
struct Quaternion{T<:Real} <: Number
    s::T
    v1::T
    v2::T
    v3::T
end

Base.real(q::Quaternion) = q.s
Base.conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3)

function Base.:*(q::Quaternion, w::Quaternion)
    return Quaternion(
        q.s * w.s - q.v1 * w.v1 - q.v2 * w.v2 - q.v3 * w.v3,
        q.s * w.v1 + q.v1 * w.s + q.v2 * w.v3 - q.v3 * w.v2,
        q.s * w.v2 - q.v1 * w.v3 + q.v2 * w.s + q.v3 * w.v1,
        q.s * w.v3 + q.v1 * w.v2 - q.v2 * w.v1 + q.v3 * w.s,
    )
end

function Base.:*(q::Quaternion, w::Union{Real,Complex,CustomComplex})
    a, b = reim(w)
    return q * Quaternion(a, b, zero(a), zero(a))
end
Base.:*(w::Union{Real,Complex,CustomComplex}, q::Quaternion) = conj(conj(q) * conj(w))

@testset "RealDot.jl" begin
    scalars = (
        randn(), randn(ComplexF64), CustomComplex(randn(2)...), Quaternion(randn(4)...)
    )
    arrays = (randn(10), randn(ComplexF64, 10))

    for inputs in (scalars, arrays)
        for x in inputs, y in inputs
            @test realdot(x, y) == real(dot(x, y))
        end
    end
end


text:
module PhasedArrayTracking
    using DocStringExtensions, StaticArrays, Tracking
    import Unitful: Hz, ms
    import Tracking: prompt, dump!

    export prompt,
        init_tracking,
        TrackingInitials

    include("discriminators.jl")
    include("tracking_loop.jl")
end


text:
using Base.Test

using MarketData
using TimeSeries

using MarketTechnicals

@testset "Volatility" begin

@testset "bollinger_bands" begin
    # TTR default uses sample=FALSE and value 117.3251
    @test isapprox(bollingerbands(cl)["up"].values[1]  , 117.3251, atol=.01)
    # TTR default uses sample=FALSE and value 89.39392
    @test isapprox(bollingerbands(cl)["down"].values[1], 89.39392, atol=.01)
    # TTR 103.3595
    @test isapprox(bollingerbands(cl)["mean"].values[1], 103.3595, atol=.01)
    @test bollingerbands(cl).timestamp[1]   == Date(2000,1,31)
    @test bollingerbands(cl).timestamp[end] == Date(2001,12,31)
end

@testset "truerange" begin
    @test isapprox(truerange(ohlc).values[end], 0.83, atol=.01) # TTR 0.83
    @test truerange(ohlc).timestamp[end] == Date(2001,12,31)
end

@testset "donchianchannels" begin
    ta = donchianchannels(ohlc)

    @test ta.meta      == ohlc.meta
    @test ta.timestamp == ohlc[20:end].timestamp
    @test isapprox(ta["up"].values[1]   , 121.5, atol=.01)
    @test isapprox(ta["down"].values[1] , 86.5, atol=.01)
    @test ta["mid"].values[1] == (ta["up"].values[1] + ta["down"].values[1]) / 2
end

@testset "atr" begin
    @test isapprox(atr(ohlc).values[1]  , 8.343571428571428)   # TTR value 8.343571
    @test isapprox(atr(ohlc).values[end], 0.9664561242651976)  # TTR value 0.9664561
    @test atr(ohlc).timestamp[1] == Date(2000,1,24)
end

@testset "keltner_bands" begin
    ta = keltnerbands(ohlc)

    @test ta["kup"].values[end] > ta["kma"].values[end]
    @test ta["kma"].values[end] > ta["kdn"].values[end]

    @test isapprox(ta["kup"].values[end], 23.3156, atol=.01)  # needs confirmation
    @test isapprox(ta["kma"].values[end], 21.3705, atol=.01)  # needs confirmation
    @test isapprox(ta["kdn"].values[end], 19.4254, atol=.01)  # needs confirmation
    @test ta.timestamp[1]   == Date(2000, 2, 1)
    @test ta.timestamp[end] == Date(2001, 12, 31)
end

@testset "chaikin_volatility" begin
    ta = chaikinvolatility(ohlc)

    @test ta.timestamp == ohlc[10+10:end].timestamp
    @test ta.colnames  == ["chaikinvolatility"]
    @test ta.meta      == ohlc.meta
    @test isapprox(ta.values[1], -2.2401, atol=.01)  # needs confirmation
    @test isapprox(ta.values[2], -3.2901, atol=.01)  # needs confirmation
end

end  # @testset "Volatility"


text:
struct Layer
    transformations::Tuple
    data::Any
    positional::Tuple
    named::NamedTuple
end

Layer(transformations::Tuple=()) = Layer(transformations, nothing, (), (;))

data(df) = Layer((), columns(df), (), (;))
mapping(args...; kwargs...) = Layer((), nothing, args, values(kwargs))

function Base.:*(l1::Layer, l2::Layer)
    t1, t2 = l1.transformations, l2.transformations
    d1, d2 = l1.data, l2.data
    p1, p2 = l1.positional, l2.positional
    n1, n2 = l1.named, l2.named
    transformations = (t1..., t2...)
    data = isnothing(d2) ? d1 : d2
    positional =(p1..., p2...)
    named = merge(n1, n2)
    return Layer(transformations, data, positional, named)
end


text:
# %% imports

include("../../Modules/save_param.jl")
include("../../Modules/Networks/net_arch.jl")

# %% usefull functions

function borne(x,a,b)
    y = a + 0.5*(b-a)*(1+cos(pi*x/10))
    return y
end

# %% function that creates the network

function set_weight!(net::Network, genes_w::Array{Float64})
    w = net.s.e
    for mat_idx in 1:length(w)
        for idx in eachindex(w[mat_idx])
            offset = sum([length(mat) for mat in w[1:mat_idx-1]])
            net.s.e[mat_idx][idx] = borne(genes_w[idx+offset], 0.0, 1.1)
        end
    end
end

function set_weight!(net::Network, genes_w::Array{Float64}, exc::Array{Int64}, inh::Array{Int64})
    w = net.s.e
    inh_count = 0
    for mat_idx in 1:length(w)
        for idx in eachindex(w[mat_idx])
            offset = sum([length(mat) for mat in w[1:mat_idx-1]]) - inh_count
            offset_neuron = sum(net.arch[1:mat_idx-1])
            (idx%size(w[mat_idx])[1]==0) ? line = size(w[mat_idx])[1] : line = idx%size(w[mat_idx])[1]
            current_neuron = line+offset_neuron
            if current_neuron in exc
                # net.s.e[mat_idx][idx] = abs(genes_w[idx+offset])
                net.s.e[mat_idx][idx] = borne(genes_w[idx+offset], 0.0, 1.1)
            elseif current_neuron in inh
                # net.s.e[mat_idx][idx] = - abs(genes_w[end])
                net.s.e[mat_idx][idx] = borne(genes_w[end], -2.0, 0.0)
                inh_count += 1
            end
        end
    end
end

function ei_idx(arch::Array{Int64})
    inh = Int64[]
    exc = Int64[]
    n_inh = [div(arch[i],4) for i in 1:length(arch)-1]
    push!(n_inh, 0)
    for l in 1:length(arch)
        for neuron in 1:arch[l]
            offset = sum(arch[1:l-1])
            if neuron <= n_inh[l]
                push!(inh, neuron+offset)
            else
                push!(exc, neuron+offset)
            end
        end
    end
    return exc, inh
end

function init_teacher_student(filepath::String, indiv::Int64)
    optim_param = YAML.load(open("Julia/CartPole/LIF_base/cfg_cmaes.yml"))
    arch = optim_param["arch"]
    results = load_param(filepath)
    genes = results["weights"][indiv]
    genes_alpha = map(x->borne(x,-1.0, 1.0),genes[1:arch[1]*4])
    matalpha = reshape(genes_alpha, (arch[1],4))
    genes_w = genes[arch[1]*4+1:end]
    params = YAML.load(open("Julia/Cartpole/Teacher/cfg_default.yml"))

    n_stud = Network(arch,params)
    n_stud.param["matalpha"] = matalpha
    n_teach = deepcopy(n_stud)

    # copy the weight of best indiv
    exc, inh = ei_idx(arch)
    set_weight!(n_teach,genes_w, exc, inh)
    return n_teach, n_stud
end

# # %% tests
#
#
# filepath = "Julia/CartPole/LIF_base/results/results_EI888.jld"
#
# n_teach, n_stud = init_teacher_student(filepath, 19)
#
# @show n_teach.param["matalpha"]

# # %% tests viz
#
# weight_viz(n_stud)
# weight_viz(n_teach)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using GSL, SpecialFunctions
using Test, LinearAlgebra

@testset "Special functions" begin

    @testset "bessel functions" begin

        x = rand() - 0.5 + (rand() - 0.5) * im

        @test diffbessel(besselj,2,x,2) ≈ diffbesselj(2,x,2)
        @test diffbessel(hankelh1,2,x,2) ≈ diffhankelh1(2,x,2)

        @test diffbessel(besselj,4,x,3) ≈ diffbesselj(4,x,3)
        @test diffbessel(hankelh1,4,x,3) ≈ diffhankelh1(4,x,3)

        # spherical bessel functions
        @test shankelh1(1, x) ≈ - exp(im*x) * (x + im) / (x^2)
        @test sbesselj(1, x) ≈ sin(x)/x^2 - cos(x)/x

        @test sbesselj(1, eps(zero(Float64))) ≈ 0.0
        @test sbesselj(0, eps(zero(Float64))) ≈ 1.0

        n = 1
        @test 2 * diffsbessel(shankelh1,n,x) ≈ shankelh1(n-1, x) - (shankelh1(n, x) + x*shankelh1(n+1, x))/x

        @test diffsbessel(shankelh1,n,x) ≈ diffshankelh1(n,x)
        @test diffsbessel(sbesselj,n,x) ≈ diffsbesselj(n,x)

        n = 2
        @test 2 * diffsbessel(shankelh1,n,x) ≈ shankelh1(n-1, x) - (shankelh1(n, x) + x*shankelh1(n+1, x))/x

        xs = 0.01:0.9:11.0;
        n = 15; ns = 0:n

        # Compare with spherical bessel from GSL, which can only accept real arguments..
        errs = map(xs) do x
            maximum(abs.(sf_bessel_jl_array(n,x) - [sbesselj(n,x) for n in ns]))
        end;
        @test maximum(errs) < 20*eps(Float64)

        errs = map(xs) do x
            maximum(
                abs.(
                    im .* sf_bessel_yl_array(n,x) -
                    [shankelh1(n,x) for n in ns] +
                    sf_bessel_jl_array(n,x)
                ) ./ abs.(sf_bessel_yl_array(n,x))
            )
        end;
        @test maximum(errs) < 1e4*eps(Float64)

    end

    @testset "Associated legendre functions" begin
        x = rand(1)[1] * 0.99
        l_max = 2

        # a few associated legendre functions without the Condon-Shortly factor
        Plm_arr = [1,x,sqrt(1-x^2),(3x^2-1)/2, 3x*sqrt(1-x^2),3*(1-x^2)]

        ls, ms = associated_legendre_indices(l_max)

        @test sf_legendre_array(GSL_SF_LEGENDRE_NONE, l_max, x)[1:length(ls)] ≈ Plm_arr

        sph_factors = (-1).^ms .* sqrt.((2 .* ls .+ 1) ./ (4pi) .* factorial.(ls-ms) ./ factorial.(ls+ms))

        condon_phase = (-1).^ms

        @test condon_phase .* sf_legendre_array(GSL_SF_LEGENDRE_SPHARM, l_max, x)[1:length(ls)] ≈ sph_factors .* Plm_arr

    end

    @testset "Spherical harmonics" begin
        θ = rand(1)[1] * 0.99
        φ = rand(1)[1] * 0.99

        l_max = 8 # small l_max due to factorial formula below

        ls, ms = spherical_harmonics_indices(l_max)
        sphs = spherical_harmonics(l_max, θ, φ)

        @test maximum(i - lm_to_spherical_harmonic_index(ls[i],ms[i]) for i in eachindex(ls)) == 0

        # check special case l == abs(m)
        is = findall(ls .== abs.(ms))

        for i in is
            @test sphs[i] ≈ (- sign(ms[i]))^ls[i] / (2^ls[i] * factorial(ls[i])) *
                sqrt(factorial(2*ls[i] + 1) / (4pi)) * sin(θ)^ls[i] * exp(im * ms[i] * φ)
        end

        l_max = 30
        ls, ms = spherical_harmonics_indices(l_max)
        sphs = spherical_harmonics(l_max, θ, φ)

        # special case m == 0, reduce to just Legendre polynomials
        Ps = sf_legendre_Pl_array(l_max, cos(θ))
        is = findall(ms .== 0)

        for l in 0:l_max
            @test sphs[is][l+1] ≈ sqrt((2l+1)/(4pi)) * Ps[l+1]
        end

        #sphs[inds] .≈ sqrt.((2 .* (0:l_max) .+ 1) ./ (4pi)) .* Ps

        # special case, north pole
        θ = 0.0
        sphs = spherical_harmonics(l_max, θ, φ)

        for i in eachindex(sphs)
            if ms[i] == 0
                @test sphs[i] ≈ sqrt((2*ls[i] + 1)/(4pi))
            else
                @test sphs[i] ≈ 0.0
            end
            i += 1
        end

    end

    @testset "Gaunt and Wigner symbols" begin
        l1 = rand(1:100)
        l2 = rand(1:100)
        l3 = rand(abs(l1-l2):2:(l1+l2)) # guarantees that iseven(l1+l2+l3)

        @test iseven(l1+l2+l3)

        m1 = rand(-l1:l1)
        m2 = rand(-l2:l2)
        m3 = m1-m2

        @test_throws(DomainError,gaunt_coefficient(l1,2*l1,l2,m2,l3,m3))
        @test_throws(MethodError,gaunt_coefficient(l1,m1,l2,m2,0.1,m3))

        # the spherical harmonics linearisation formula
        θ, φ = rand(2) .* pi

        l_small = 6
        l_max = 2*l_small # needs to be larger than l_small

        ls, ms = spherical_harmonics_indices(l_max);
        Ys = spherical_harmonics(l_max, θ, φ);

        cs = reshape(
            [
                gaunt_coefficient(l1,m1,l2,m2,l3,m3)

                    for l3 = 0:l_max for m3 = -l3:l3
                for l2 = 0:l_small for m2 = -l2:l2
            for l1 = 0:l_small for m1 = -l1:l1]
        , ((l_small + 1)^2,(l_small + 1)^2,(l_max+1)^2));

        for n2 in 1:(l_small + 1)^2, n3 in 1:(l_small + 1)^2
            @test 4pi * (-1)^(ms[n3]+ms[n2]) * (1.0im)^(ls[n3]-ls[n2]) * Ys[n3] * conj(Ys[n2]) ≈
                sum( (1.0im).^(-ls) .* (-1.0).^ms .* conj.(Ys) .* cs[n2,n3,:]) atol = 1e-12
        end

    end

    @testset "complex radial coordinates" begin

    # Test 3-dimensional transforms
        xs = [rand(-1.01:0.1:1.0,3) + rand(-1.01:0.1:1.0,3)*im for i = 1:100]
        rθφs = cartesian_to_radial_coordinates.(xs)
        @test maximum(norm.(xs - radial_to_cartesian_coordinates.(rθφs))) < 2e-14

        xs = [rand(-1.01:0.1:1.0,3) for i = 1:100]
        rθφs = cartesian_to_radial_coordinates.(xs)

        @test maximum(rθφ[2] for rθφ in rθφs) <= pi
        @test minimum(rθφ[2] for rθφ in rθφs) >= 0.0

        @test pi/2 < maximum(rθφ[3] for rθφ in rθφs) <= pi
        @test -pi <= minimum(rθφ[3] for rθφ in rθφs) < -pi/2

    # Test 2-dimensional transforms
        xs = [rand(-1.01:0.1:1.0,2) + rand(-1.01:0.1:1.0,2)*im for i = 1:100]
        rθs = cartesian_to_radial_coordinates.(xs)

        @test maximum(norm.(xs - radial_to_cartesian_coordinates.(rθs))) < 1e-14
    end
end


text:
struct ConvRBM{V,H,W,K,K2}
    visible::V
    hidden::H
    w::W
    stride::NTuple{K,Int}
    pad::NTuple{K2,Int}
    dilation::NTuple{K,Int}
    groups::Int
    pool::Bool
    function ConvRBM(
        visible::AbstractLayer, hidden::AbstractLayer, w::AbstractArray,
        stride::NTuple{K,Int}, pad::NTuple{K2,Int}, dilation::NTuple{K,Int}, groups::Int,
        pool::Bool = false
    ) where {K,K2}
        @assert K == ndims(w) - ndims(visible) - ndims(hidden)
        @assert K2 == 2K
        kernel_size = size(w)[(ndims(visible) + 1):(end - ndims(hidden))]
        @assert size(w) == (size(visible)..., kernel_size..., size(hidden)...)
        return new{typeof(visible), typeof(hidden), typeof(w), K, K2}(
            visible, hidden, w, stride, pad, dilation, groups, pool
        )
    end
end

"""
    ConvRBM(visible, hidden, weights; stride = 1, pad = 0, dilation = 1, groups = 1, pool = false)

Convolutional RBM.

`pad == :same` uses `samepad(...)` to infer a padding that preserves spatial dimensions.
"""
function ConvRBM(
    visible::AbstractLayer, hidden::AbstractLayer, w::AbstractArray;
    stride = 1, pad = 0, dilation = 1, groups::Int = 1, pool::Bool = false
)
    K = ndims(w) - ndims(visible) - ndims(hidden)
    kernel_size = size(w)[(ndims(visible) + 1):(ndims(visible) + K)]
    return ConvRBM(
        visible, hidden, w,
        expand_tuple(Val(K), stride),
        parsepad(kernel_size, dilation, pad),
        expand_tuple(Val(K), dilation),
        groups, pool
    )
end

"""
    ConvRBM(visible, hidden, kernel_size; kwargs...)

Convolutional RBM with given `kernel_size`.
"""
function ConvRBM(
    visible::AbstractLayer, hidden::AbstractLayer, kernel_size::TupleN{Int}, ::Type{T} = Float64;
    kwargs...
) where {T}
    w = zeros(T, size(visible)..., kernel_size..., size(hidden)...)
    return ConvRBM(visible, hidden, w; kwargs...)
end

RBMs.visible(rbm::ConvRBM) = rbm.visible
RBMs.hidden(rbm::ConvRBM) = rbm.hidden
RBMs.weights(rbm::ConvRBM) = rbm.w

channel_size(rbm::ConvRBM) = size(visible(rbm))
channel_length(rbm::ConvRBM) = length(visible(rbm))
channel_ndims(rbm::ConvRBM) = ndims(visible(rbm))

kernel_size(rbm::ConvRBM) = size(weights(rbm))[input_dims(rbm)]
kernel_ndims(rbm::ConvRBM) = length(kernel_size(rbm))

input_dims(rbm::ConvRBM) = (channel_ndims(rbm) + 1):(ndims(weights(rbm)) - ndims(hidden(rbm)))
output_dims(rbm::ConvRBM) = (ndims(hidden(rbm)) + 1):(ndims(hidden(rbm)) + kernel_ndims(rbm))

"""
    vsizes(convrbm, v) -> (channel_size, input_size, batch_size)

Returns a (named) tuple decomposition of the size of `v`, such that:

    size(v) == (channel_size..., input_size..., batch_size...)

Throws an error if `v` is not consistent with this size.
"""
function vsizes(rbm::ConvRBM, v::AbstractArray)
    input_size = size(v)[(channel_ndims(rbm) + 1):(channel_ndims(rbm) + kernel_ndims(rbm))]
    batch_size = size(v)[(channel_ndims(rbm) + kernel_ndims(rbm) + 1):end]
    @assert size(v) == (channel_size(rbm)..., input_size..., batch_size...)
    return (channel_size = channel_size(rbm), input_size = input_size, batch_size = batch_size)
end

"""
    hsizes(convrbm, h) -> (hidden_size, output_size, batch_size)

Returns a (named) tuple decomposition of the size of `h`, such that:

    size(h) == (hidden_size..., output_size..., batch_size...)

Throws an error if `h` is not consistent with this size.
"""
function hsizes(rbm::ConvRBM, h::AbstractArray)
    output_size = size(h)[(ndims(hidden(rbm)) + 1):(ndims(hidden(rbm)) + kernel_ndims(rbm))]
    batch_size = size(h)[(ndims(hidden(rbm)) + kernel_ndims(rbm) + 1):end]
    @assert size(h) == (size(hidden(rbm))..., output_size..., batch_size...)
    return (hidden_size = size(hidden(rbm)), output_size = output_size, batch_size = batch_size)
end

function RBMs.batch_size(rbm::ConvRBM, v::AbstractArray, h::AbstractArray)
    vsz = vsizes(rbm, v)
    hsz = hsizes(rbm, h)
    if !isempty(vsz.batch_size) && !isempty(hsz.batch_size)
        @assert vsz.batch_size == hsz.batch_size
        return vsz.batch_size
    elseif isempty(hsz.batch_size)
        return vsz.batch_size
    elseif isempty(vsz.batch_size)
        return hsz.batch_size
    end
end

RBMs.batch_size(rbm::ConvRBM, v::AbstractArray) = vsizes(rbm, v).batch_size

"""
    output_size(rbm, v)

Output size of the convolution.

If `v` has `batch_size` batches, then `output_size(rbm, v)` returns `output_size` such that:

    size(h) == (hidden_size..., output_size..., batch_size...)
"""
output_size(rbm::ConvRBM, v::AbstractArray) = output_size(rbm, vsizes(rbm, v).input_size)
output_size(rbm::ConvRBM, input_size::TupleN{Int}) = output_size(
    kernel_size(rbm), input_size, rbm.stride, rbm.pad, rbm.dilation
)

function RBMs.inputs_v_to_h(rbm::ConvRBM, v::AbstractArray)
    vsz = vsizes(rbm, v)
    wflat = reshape(rbm.w, channel_length(rbm), kernel_size(rbm)..., length(hidden(rbm)))
    vflat = reshape(v, channel_length(rbm), vsz.input_size..., prod(vsz.batch_size))
    Iflat = conv_v2h(
        wflat, activations_convert_maybe(wflat, vflat);
        rbm.stride, rbm.pad, rbm.dilation, rbm.groups
    )
    @assert size(Iflat) == (length(hidden(rbm)), size(Iflat)[2:end-1]..., prod(vsz.batch_size))
    return reshape(Iflat, size(hidden(rbm))..., size(Iflat)[2:end-1]..., vsz.batch_size...)
end

function RBMs.inputs_h_to_v(rbm::ConvRBM, h::AbstractArray)
    hsz = hsizes(rbm, h)
    wflat = reshape(rbm.w, length(visible(rbm)), kernel_size(rbm)..., length(hidden(rbm)))
    hflat = reshape(h, length(hidden(rbm)), hsz.output_size..., prod(hsz.batch_size))
    Iflat = conv_h2v(
        wflat, activations_convert_maybe(wflat, hflat);
        rbm.stride, rbm.pad, rbm.dilation, rbm.groups
    )
    @assert size(Iflat) == (channel_length(rbm), size(Iflat)[2:end-1]..., prod(hsz.batch_size))
    return reshape(Iflat, channel_size(rbm)..., size(Iflat)[2:end-1]..., hsz.batch_size...)
end

function RBMs.energy(rbm::ConvRBM, v::AbstractArray, h::AbstractArray)
    Ev = energy(visible(rbm), v)
    @assert size(Ev) == (vsizes(rbm, v).input_size..., vsizes(rbm, v).batch_size...)
    Eh = energy(hidden(rbm), h)
    @assert size(Eh) == (hsizes(rbm, h).output_size..., hsizes(rbm, h).batch_size...)
    Ew = interaction_energy(rbm, v, h)
    @assert size(Ew) == batch_size(rbm, v, h)
    Ev_reduced = reshape_maybe(sum(Ev; dims=1:kernel_ndims(rbm)), vsizes(rbm, v).batch_size)
    Eh_reduced = reshape_maybe(sum(Eh; dims=1:kernel_ndims(rbm)), hsizes(rbm, h).batch_size)
    E = Ev_reduced .+ Eh_reduced .+ Ew
    @assert size(E) == batch_size(rbm, v, h)
    return E
end

function RBMs.interaction_energy(rbm::ConvRBM, v::AbstractArray, h::AbstractArray)
    bsz = batch_size(rbm, v, h) # do size checks early
    I = inputs_v_to_h(rbm, v) # takes care of the convolution
    @assert size(I) == (size(hidden(rbm))..., hsizes(rbm, h).output_size..., vsizes(rbm, v).batch_size...)
    E = -sum(I .* h; dims=1:(ndims(hidden(rbm)) + kernel_ndims(rbm)))
    return reshape_maybe(E, bsz)
end

function RBMs.free_energy(rbm::ConvRBM, v::AbstractArray; β::Real = 1)
    if rbm.pool
        return free_energy_pooled(rbm, v; β)
    else
        return free_energy_nopool(rbm, v; β)
    end
end

function RBMs.sample_h_from_v(rbm::ConvRBM, v::AbstractArray; β::Real = true)
    if rbm.pool
        return sample_h_from_v_pooled(rbm, v; β)
    else
        return sample_h_from_v_nopool(rbm, v; β)
    end
end

function RBMs.sample_v_from_h(rbm::ConvRBM, h::AbstractArray; β::Real = true)
    inputs = RBMs.inputs_h_to_v(rbm, h)
    return RBMs.transfer_sample(visible(rbm), inputs; β)
end

function RBMs.sample_v_from_v(rbm::ConvRBM, v::AbstractArray; β::Real = true, steps::Int = 1)
    for _ in 1:steps
        v = oftype(v, RBMs.sample_v_from_v_once(rbm, v; β))
    end
    return v
end

function RBMs.sample_h_from_h(rbm::ConvRBM, h::AbstractArray; β::Real = true, steps::Int = 1)
    for _ in 1:steps
        h = oftype(h, RBMs.sample_h_from_h_once(rbm, h; β))
    end
    return h
end

function RBMs.sample_v_from_v_once(rbm::ConvRBM, v::AbstractArray; β::Real = true)
    h = RBMs.sample_h_from_v(rbm, v; β)
    v = RBMs.sample_v_from_h(rbm, h; β)
    return v
end

function RBMs.sample_h_from_h_once(rbm::ConvRBM, h::AbstractArray; β::Real = true)
    v = RBMs.sample_v_from_h(rbm, h; β)
    h = RBMs.sample_h_from_v(rbm, v; β)
    return h
end


text:
using Test, Random, LinearAlgebra
using Yao
using YaoExtensions

@testset "block extension" begin
    include("block_extension/blocks.jl")
end

@testset "easybuild" begin
    include("easybuild/easybuild.jl")
end

@testset "diff_extension" begin
    include("diff_extension/diffs.jl")
end

@testset "gatecount" begin
    include("gatecount.jl")
end