text:
using LabReports, Test
using CSV
using Unitful
using LabReports: sign_change, integration_domain, discharge_area,
    CVCapacitanceReport, CDCapacitanceReport

@testset "IO" begin
    df = CSV.read("fake_data/15/15_EIS", delim=';', datarow=2, copycols=true)
    @test df == results("fake_data", "EIS", "15", "15", processed=false)[2]

    df = CSV.read("fake_data/15/15_CV_3.dat", delim=';', datarow=3, copycols=true)
    @test df == results("fake_data", "CV", "15", "3", processed=true)[2]
end

@testset "CV Area" begin
    df = results("fake_data", "CV", "15", "14000", processed=true)[2]

    I = df[!, Symbol("Current (mA)")]
    V = df[!, Symbol("Potential (V)")]

    @testset "sign_change" begin
        @test sign_change(V, :- => :+) == 55
        @test sign_change(V, :+ => :-) == 33
        @test sign_change(I, :+ => :-) == 20
        @test sign_change(I, :- => :+) == 50
    end

    @testset "integration_domain" begin
        @test integration_domain(V, I, 1, :orar) == (55, 20)
        @test integration_domain(V, I, 2, :orar) == (50, 55)
        @test integration_domain(V, I, 3, :orar) == (33, 50)
        @test integration_domain(V, I, 4, :orar) == (20, 33)
    end

    @testset "4th quadrant" begin
        Δt, ΔV, area = discharge_area(df, 4, nothing)
        @test Δt ≈ 0.1299999996717u"s"
        @test ΔV ≈ -2.613535594192u"V"
        @test area ≈ -0.00190424207227u"V*A"
    end

    @testset "4th quadrant, fixed ΔV" begin
        Δt, ΔV, area = discharge_area(df, 4, 1.5u"V")
        @test Δt ≈ 0.1299999996717u"s"
        @test ΔV ≈ -1.5u"V"
        @test area ≈ -0.00190424207227u"V*A"
    end
end

@testset "CV capacitance" begin
    datafile, df = results("fake_data", "CV", "15", "14000", processed=true)

    setup = CSetup(a=0.1u"cm^2", A=10.0u"cm^2", fixed_ΔV=1.5u"V")
    cr = CVCapacitanceReport(datafile, df, 4, "fake_data", setup)
    @test cr.C ≈ 90.67819391u"μF"
    @test cr.C_specific ≈ 90678.19391788u"μF/g"
    @test cr.E ≈ 2.83369355e-8u"W*hr"
    @test cr.E_specific ≈ 0.0283369355u"W*hr/kg"
    @test cr.P ≈ 0.0007847151416u"W"
    @test cr.P_specific ≈ 784.71514165u"W/kg"
end

@testset "C&D capacitance" begin
    only_d(files) = begin
        for f in files
            if endswith(f.filename, "_D")
                return [f]
            end
        end
    end
    datafile, df = results("fake_data", "C&D", "15", "2e-3", processed=false, reduction=only_d)

    setup = CSetup(a=0.1u"cm^2", A=10.0u"cm^2", fixed_ΔV=1.4u"V")
    cr = CDCapacitanceReport(datafile, df, "fake_data", setup)
end


text:
# functions to describe a HybridNetwork object to avoid accessing the attributes directly
# Claudia August 2015

"""
    tipLabels(x)

Return a vector of taxon names at the leaves, for objects of various types:
`HybridNetwork`,
Vector of `HybridNetwork`s (in which case the union is taken then sorted),
Vector of `Quartet`s, `DataCF`,
`TraitSimulation`, `MatrixTopologicalOrder`.

For a network, the taxon names are coerced to strings.
"""
# see readData.jl and traits.jl
function tipLabels(net::HybridNetwork)
    return String[l.name for l in net.leaf] # AbstractString does not work for use by tree2Matrix
end

"""
`fittedQuartetCF(d::DataCF, format::Symbol)`

return a data frame with the observed and expected quartet concordance factors
after estimation of a network with snaq(T,d).
The format can be :wide (default) or :long.

- if wide, the output has one row per 4-taxon set, and each row has 10 columns: 4 columns
  for the taxon names, 3 columns for the observed CFs and 3 columns for the expected CF.
- if long, the output has one row per quartet, i.e. 3 rows per 4-taxon sets, and 7 columns:
  4 columns for the taxon names, one column to give the quartet resolution, one column for
  the observed CF and the last column for the expected CF.

see also: `topologyQPseudolik!` and `topologyMaxQPseudolik!` to update the fitted CF expected
under a specific network, inside the DataCF object `d`.
"""
function fittedQuartetCF(d::DataCF, format=:wide::Symbol)
    if format == :wide
        df=DataFrame(
                 tx1 = [q.taxon[1] for q in d.quartet],
                 tx2 = [q.taxon[2] for q in d.quartet],
                 tx3 = [q.taxon[3] for q in d.quartet],
                 tx4 = [q.taxon[4] for q in d.quartet],
                obsCF12=[q.obsCF[1] for q in d.quartet],
                obsCF13=[q.obsCF[2] for q in d.quartet],
                obsCF14=[q.obsCF[3] for q in d.quartet],
                expCF12=[q.qnet.expCF[1] for q in d.quartet],
                expCF13=[q.qnet.expCF[2] for q in d.quartet],
                expCF14=[q.qnet.expCF[3] for q in d.quartet])
    elseif format == :long
        nQ = length(d.quartet)
        df = DataFrame(
            tx1 = repeat([q.taxon[1] for q in d.quartet], inner=[3]),
            tx2 = repeat([q.taxon[2] for q in d.quartet], inner=[3]),
            tx3 = repeat([q.taxon[3] for q in d.quartet], inner=[3]),
            tx4 = repeat([q.taxon[4] for q in d.quartet], inner=[3]),
            quartet = repeat(["12_34","13_24","14_23"], outer=[nQ]),
            obsCF = Array{Float64}(3*nQ),
            expCF = Array{Float64}(3*nQ)  )
        row = 1
        for i in 1:nQ
            for j in 1:3
                df[row, 6] = d.quartet[i].obsCF[j]
                df[row, 7] = d.quartet[i].qnet.expCF[j]
                row += 1
            end
        end
    else
        error("format $(format) was not recognized. Should be :wide or :long.")
    end
    return df
end

"""
`setNonIdBL!(net)`

Set non-identifiable edge branch lengths to -1.0 (i.e. missing) for a level-1 network `net`,
except for edges in

- a good triangle: the edge below the hybrid is constrained to 0.
- a bad diamond II: the edge below the hybrid is constrained to 0
- a bad diamond I: the edges across from the hybrid node have non identifiable lengths
  but are kept, because the two γ*(1-exp(-t)) values are identifiable.

will break if `inCycle` attributes are not initialized (at -1) or giving a correct node number.
"""
# 'good' triangle in the paper = bad triangle in hyb node attribute: Not extremely or very bad
function setNonIdBL!(net::HybridNetwork)
    for e in net.edge
        if(!e.istIdentifiable)
            keeplength = any(n -> (n.isBadDiamondII || n.isBadTriangle), e.node)
            # if below 'bad' hybrid node, length=0 by constraint. If above, length estimated.
            # next: keep length if edge across from hybrid node in bad diamond I.
            if (!keeplength && e.inCycle != -1)
                hyb = net.node[getIndexNode(e.inCycle, net)]
                if (hyb.isBadDiamondI)
                    keeplength |= !any(n -> (n==hyb), e.node) # only if e across hyb: no touching it
                end
            end
            if !keeplength
                e.length = -1.0 # not setLength, which does not allow too negative BLs
            end
        end
    end
end

"""
setBLGammaParsimony!(net::HybridNetwork)

Maximum parsimony function does not provide estimates for branch lengths, or gamma.
But since the `maxParsimonyNet` function is using snaq move functions, branch lengths
and gamma values are set randomly (to initialize optimization).
We need to remove these random values before returning the maximum parsimony network.
"""
function setBLGammaParsimony!(net::HybridNetwork)
    for e in net.edge
        e.length = -1.0
        if e.hybrid
            e.gamma = -1.0
        end
    end
end

# function that we need to overwrite to avoid printing useless scary
# output for HybridNetworks
# PROBLEM: writeTopology changes the network and thus show changes the network
function Base.show(io::IO, obj::HybridNetwork)
    disp = "$(typeof(obj)), "
    if obj.isRooted
        disp = disp * "Rooted Network"
    else
        disp = disp * "Un-rooted Network"
    end
    disp = disp * "\n$(obj.numEdges) edges\n"
    disp = disp * "$(obj.numNodes) nodes: $(obj.numTaxa) tips, "
    disp = disp * "$(obj.numHybrids) hybrid nodes, "
    disp = disp * "$(obj.numNodes - obj.numTaxa - obj.numHybrids) internal tree nodes.\n"
    tipslabels = [n.name for n in obj.leaf]
    if length(tipslabels) > 1 || !all(tipslabels .== "")
        disptipslabels = "$(tipslabels[1])"
        for i in 2:min(obj.numTaxa, 4)
            disptipslabels = disptipslabels * ", $(tipslabels[i])"
        end
        if obj.numTaxa > 4 disptipslabels = disptipslabels * ", ..." end
        disp *= "tip labels: " * disptipslabels
    end
    par = ""
    try
        # par = writeTopology(obj,round=true) # but writeTopology changes the network, not good
        s = IOBuffer()
        writeSubTree!(s, obj, false,true, true,3)
        par = String(s)
    catch err
        println("ERROR with writeSubTree!:")
        showerror(STDOUT, err)
        println("Trying writeTopologyLevel1")
        par = writeTopologyLevel1(obj)
    end
    disp *= "\n$par"
    println(io, disp)
end

# and QuartetNetworks (which cannot be just written because they do not have root)
function Base.show(io::IO, net::QuartetNetwork)
    print(io,"taxa: $(net.quartetTaxon)\n")
    print(io,"number of hybrid nodes: $(net.numHybrids)\n")
    if(net.split != [-1,-1,-1,-1])
        print(io,"split: $(net.split)\n")
    end
end

function Base.show(io::IO,d::DataCF)
    print(io,"Object DataCF\n")
    print(io,"number of quartets: $(d.numQuartets)\n")
    if(d.numTrees != -1)
        print(io,"number of trees: $(d.numTrees)\n")
    end
end

function Base.show(io::IO,q::Quartet)
    print(io,"number: $(q.number)\n")
    print(io,"taxon names: $(q.taxon)\n")
    print(io,"observed CF: $(q.obsCF)\n")
    print(io,"pseudo-deviance under last used network: $(q.logPseudoLik) (meaningless before estimation)\n")
    print(io,"expected CF under last used network: $(q.qnet.expCF) (meaningless before estimation)\n")
    if(q.ngenes != -1)
        print(io,"number of genes used to compute observed CF: $(q.ngenes)\n")
    end
end

function Base.show(io::IO, obj::Node)
    disp = "$(typeof(obj)):"
    disp = disp * "\n number:$(obj.number)"
    if (obj.name != "") disp *= "\n name:$(obj.name)" end
    if (obj.hybrid)     disp *= "\n hybrid node" end
    if (obj.leaf)       disp *= "\n leaf node" end
    disp *= "\n attached to $(length(obj.edge)) edges, numbered:"
    for e in obj.edge disp *= " $(e.number)"; end
    println(io, disp)
end

function Base.show(io::IO, obj::Edge)
    disp = "$(typeof(obj)):"
    disp *= "\n number:$(obj.number)"
    disp *= "\n length:$(obj.length)"
    if (obj.hybrid)
        disp *= "\n " * (obj.isMajor ? "major" : "minor")
        disp *= " hybrid edge with gamma=$(obj.gamma)"
    elseif (!obj.isMajor)
        disp *= "\n minor tree edge"
    end
    disp *= "\n attached to $(length(obj.node)) node(s) (parent first):"
    if (length(obj.node)==1) disp *= " $(obj.node[1].number)";
    elseif (length(obj.node)==2)
        disp *= " $(obj.node[obj.isChild1 ? 2 : 1].number)"
        disp *= " $(obj.node[obj.isChild1 ? 1 : 2].number)"
    end
    println(io, disp)
end

ext: jl


text:
using ResourceContexts
using Test
using Logging
using Compat

# Use of @! to pass context to resource creation function
@! function foo(x, label)
    # Use of @defer inside a resource creation function
    @defer push!(x, label)
end

@! function bar(x; label=nothing)
    @defer push!(x, label)
end

@! function baz(x::T, label) where {T}
    @defer push!(x, label)
end

@testset "Cleanup ordering" begin
    cleanups = []
    @context begin
        @defer push!(cleanups, :A)
        @! foo(cleanups, :B)
        @! bar(cleanups; label=:C)
        @! baz(cleanups, :D)
        @test cleanups == []
    end
    @test cleanups == [:D, :C, :B, :A]
end

@testset "@context before function def" begin
    @context function func_with_context()
        cleanups = []
        @defer push!(cleanups, :A)
        @defer push!(cleanups, :B)
        @test isempty(cleanups)
        return cleanups
    end

    cleanups = func_with_context()
    @test cleanups == [:B, :A]
end

@testset "Exceptions during cleanup" begin
    try
        @context begin
            @defer error("A")
            @defer error("B")
        end
    catch exc
        stack = current_exceptions()
        @test stack[1][1] == ErrorException("B")
        @test stack[2][1] == ErrorException("A")
    end
end

@! function do_nothing()
    nothing
end

Main.eval(:(
# Eval in Main, in case test is eval'd elsewhere
function do_nothing_from_function()
    @! $(@__MODULE__).do_nothing()
end
))

@testset "Logging when using the global context" begin
    q = :(@! do_nothing())
    # Simulate call coming from the REPL in any module
    q.args[2] = LineNumberNode(1, Symbol("REPL[1]"))
    @test_logs (:debug, r"^Using global `ResourceContext`",
                Test.Ignored(), Test.Ignored(),
                Test.Ignored(), "REPL[1]", 1) #=
        =# min_level=Logging.Debug Main.eval(q)
    # Simulate call coming from Main module in global scope
    q.args[2] = LineNumberNode(1, Symbol("SomeOtherFile"))
    @test_logs (:debug,) min_level=Logging.Debug Main.eval(q)
    # Call coming from within a function
    @test_logs (:warn,) min_level=Logging.Debug Main.do_nothing_from_function()
end

@testset "enter_do — context management for `do` blocks" begin
    function invoke_user_func(f)
        fake_resource = 42
        f(fake_resource)
    end

    function fail_before_user_func(f)
        error("Oops1")
        f()
    end
    function fail_after_user_func(f)
        f()
        error("Oops2")
    end

    @context begin
        @test @!(enter_do(invoke_user_func)) == (42,)
    end

    @test try
        @context begin
            @! enter_do(fail_before_user_func)
        end
    catch e
        @test e isa TaskFailedException
        first(current_exceptions(e.task))[1]
    end == ErrorException("Oops1")

    @test try
        @context begin
            @! enter_do(fail_after_user_func)
        end
    catch e
        @test e isa TaskFailedException
        first(current_exceptions(e.task))[1]
    end == ErrorException("Oops2")
end

@testset "detach_context_cleanup" begin
    did_cleanup = false
    @noinline function use_context_and_detach()
        some_value = @context begin
            @defer did_cleanup = true
            @defer begin
                yield()  # Task switching during cleanup should be ok
            end
            @! ResourceContexts.detach_context_cleanup([1,2])
        end
        x = copy(some_value)  # prevent holding some_value live
        @test x == [1,2]
        @test !did_cleanup
    end
    use_context_and_detach()
    @test !did_cleanup
    GC.gc()
    # These yields are required because
    yield() # 1) The cleanup is run async within the finalizer
    yield() # 2) The cleanup calls yield() to assert task switching is ok
    @test did_cleanup
end

@testset "Cleanup robustness" begin
    cleanup_count = 0
    ctx = ResourceContexts.ResourceContext()
    ResourceContexts.defer(ctx) do
        cleanup_count += 1
    end
    ResourceContexts.cleanup!(ctx)
    ResourceContexts.cleanup!(ctx)
    @test cleanup_count == 1

    # Exceptions
    cleanup_count = 0
    ctx = ResourceContexts.ResourceContext()
    ResourceContexts.defer(ctx) do
        cleanup_count += 1
    end
    ResourceContexts.defer(ctx) do
        error("X")
    end
    try
        ResourceContexts.cleanup!(ctx)
    catch
    end
    # The second scheduled cleanup should still run even if the first throws
    @test cleanup_count == 1
    # Calling cleanup should not throw again
    @test isnothing(ResourceContexts.cleanup!(ctx))
    @test cleanup_count == 1
end

@testset "Manually managed contexts" begin
    @testset "Finalization" begin
        did_cleanup = false
        function forgotten_cleanup()
            ctx = ResourceContexts.ResourceContext()
            ResourceContexts.defer(ctx) do
                did_cleanup = true
            end
        end
        forgotten_cleanup()
        GC.gc()
        # Check that the finalizer of the manually managed `ctx` ran, cleaning up
        # the resources used.
        @test did_cleanup
    end
end

include("base_interop.jl")


text:
module DWT
export
    # WPD
    wpd,
    wpd!,
    iwpd,
    iwpd!,
    # Transform on all signals
    dwtall,
    wptall,
    wpdall,
    idwtall,
    iwptall,
    iwpdall

using
    Wavelets

using
    ..Utils

# ========== Wavelet Packet Decomposition ==========
# 1D Wavelet Packet Decomposition without allocated output array
"""
    wpd(x, wt[, L; standard])

Returns the wavelet packet decomposition (WPD) for L levels for input signal x.

!!! note
    Nonstandard transform is not yet available for 2D-WPD. Setting `standard = true`
    will currently result in an error message.

# Arguments
- `x::AbstractVector{T} where T<:Number` or `x::AbstractArray{T,2} where T<:Number`: Input
  vector/matrix. A vector input undergoes 1D wavelet decomposition whereas a matrix input
  undergoes 2D wavelet decomposition.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels for wavelet
  decomposition.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `::Array{T,2}` or `::Array{T,3}`: Decomposed signal. For an input vector `x`, output is a
  2D matrix where each column corresponds to a level of decomposition. For an input matrix
  `x`, output is a 3D array where each slice of dimension 3 corresponds to a level of
  decomposition.

# Examples:
```julia
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
wt = wavelet(WT.haar)
xw = wpd(x, wt)

# 2D wavelet decomposition
x = randn(8,8)
wt = wavelet(WT.haar)
xw = wpd(x, wt)
```

**See also:** [`wpd!`](@ref), [`iwpd`](@ref)
"""
function wpd(x::AbstractVector{T},
             wt::OrthoFilter,
             L::Integer=maxtransformlevels(x)) where T<:Number
    # Sanity check
    @assert isdyadic(x)
    @assert 0 ≤ L ≤ maxtransformlevels(x)

    # Allocate variable for result
    n = length(x)
    xw = Array{T,2}(undef, (n,L+1))
    # Wavelet packet decomposition
    wpd!(xw, x, wt, L)
    return xw
end

# 2D Wavelet Packet Decomposition without allocated output array
function wpd(x::AbstractArray{T,2},
             wt::OrthoFilter,
             L::Integer = maxtransformlevels(x);
             standard::Bool = true) where T<:Number
    # Sanity check
    @assert 0 ≤ L ≤ maxtransformlevels(x)

    # Allocate variable for result
    sz = size(x)
    xw = Array{T,3}(undef, (sz...,L+1))
    # Wavelet packet decomposition
    wpd!(xw, x, wt, L, standard=standard)
    return xw
end

# 1D Wavelet Packet Decomposition with allocated output array
"""
    wpd!(y, x, wt[, L; standard])

Same as `wpd` but without array allocation.

!!! note
    Nonstandard transform is not yet available for 2D-WPD. Setting `standard = true` will
    currently result in an error message.

# Arguments
- `y::AbstractArray{T,2} where T<:Number` or `y::AbstractArray{T,3} where T<:Number`: An
  allocated array to write the outputs of `x` onto.
- `x::AbstractVector{T} where T<:Number` or `x::AbstractArray{T,2} where T<:Number`: Input
  vector/matrix. A vector input undergoes 1D wavelet decomposition whereas a matrix input
  undergoes 2D wavelet decomposition.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels for wavelet
  decomposition.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `y::AbstractArray{T,2} where T<:Number` or `y::AbstractArray{T,3} where T<:Number`:
  Decomposed signal. For an input vector, output is a 2D matrix where each column
  corresponds to a level of decomposition. For an input matrix, output is a 3D array where
  each slice of dimension 3 corresponds to a level of decomposition.

# Examples:
```
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
xw = Array{eltype(x)}(undef, (8,3))
wt = wavelet(WT.haar)
wpd!(xw, x, wt)

# 2D wavelet decomposition
x = randn(8,8)
xw = Array{eltype(x)}(undef, (8,8,3))
wt = wavelet(WT.haar)
wpd!(xw, x, wt)
```

**See also:** [`wpd`](@ref), [`iwpd`](@ref)
"""
function wpd!(y::AbstractArray{T,2},
              x::AbstractArray{T,1},
              wt::DiscreteWavelet,
              L::Integer=maxtransformlevels(x)) where T<:Number
    # Sanity check
    n = length(x)
    @assert 0 ≤ L ≤ maxtransformlevels(x)
    @assert size(y) == (n,L+1)

    # Setup
    g, h = WT.makereverseqmfpair(wt, true)
    y[:,1] = x

    # Wavelet Decomposition
    for i in 0:(L-1)
        nₚ = nodelength(n, i)                   # Parent node length
        for j in 0:((1<<i)-1)
            colₚ = i+1                          # Parent column
            rngₚ = (j*nₚ+1):((j+1)*nₚ)           # Parent range
            @inbounds v = @view y[rngₚ, colₚ]    # Parent node
            colᵣ = colₚ+1                       # Child column
            nᵣ = nₚ÷2                           # Child node length
            rng₁ = (2*j*nᵣ+1):((2*j+1)*nᵣ)      # Approx coef range
            rng₂ = ((2*j+1)*nᵣ+1):(2*(j+1)*nᵣ)  # Detail coef range
            @inbounds w₁ = @view y[rng₁, colᵣ]  # Approx coef node
            @inbounds w₂ = @view y[rng₂, colᵣ]  # Detail coef node
            dwt_step!(w₁, w₂, v, h, g)          # Decompose
        end
    end
    return y
end

# 2D Wavelet Packet Decomposition with allocated output array
function wpd!(y::AbstractArray{T,3},
              x::AbstractArray{T,2},
              wt::OrthoFilter,
              L::Integer = maxtransformlevels(x);
              standard::Bool = true) where T<:Number
    # Sanity check
    m, n = size(x)
    @assert 0 ≤ L ≤ maxtransformlevels(x)
    @assert size(y) == (m, n,L+1)

    # Setup
    g, h = WT.makereverseqmfpair(wt, true)   # low & high pass filters
    temp = Array{T,2}(undef, (m,n))

    # First slice of y is level 0, ie. original signal
    @inbounds y[:,:,1] = x
    # ----- Compute L levels of decomposition -----
    for i in 0:(L-1)
        # Parent node width and height
        mₚ = nodelength(m, i)
        nₚ = nodelength(n, i)
        # Iterate over each nodes at current level
        lrange = 0:((1<<i)-1)
        for j in lrange
            for k in lrange
                # Extract parent node
                sliceₚ = i+1
                rng_row = (j*mₚ+1):((j+1)*mₚ)
                rng_col = (k*nₚ+1):((k+1)*nₚ)
                @inbounds v = @view y[rng_row, rng_col, sliceₚ]
                # Extract children nodes of current parent
                sliceᵣ = sliceₚ+1
                mᵣ = mₚ÷2
                nᵣ = nₚ÷2
                @inbounds w₁ = @view y[(2*j*mᵣ+1):((2*j+1)*mᵣ), (2*k*nᵣ+1):((2*k+1)*nᵣ), sliceᵣ]
                @inbounds w₂ = @view y[(2*j*mᵣ+1):((2*j+1)*mᵣ), ((2*k+1)*nᵣ+1):(2*(k+1)*nᵣ), sliceᵣ]
                @inbounds w₃ = @view y[((2*j+1)*mᵣ+1):(2*(j+1)*mᵣ), (2*k*nᵣ+1):((2*k+1)*nᵣ), sliceᵣ]
                @inbounds w₄ = @view y[((2*j+1)*mᵣ+1):(2*(j+1)*mᵣ), ((2*k+1)*nᵣ+1):(2*(k+1)*nᵣ), sliceᵣ]
                # Extract temp subarray (Same size as parent node)
                @inbounds tempₖ = @view temp[rng_row, rng_col]
                # Perform 1 level of wavelet decomposition
                @inbounds dwt_step!(w₁, w₂, w₃, w₄, v, h, g, tempₖ, standard=standard)
            end
        end
    end
    return y
end

# ========== Inverse Wavelet Packet Decomposition ==========
# IWPD by level without allocation
"""
    iwpd(xw, wt[, L; standard])
    iwpd(xw, wt, tree[; standard])

Computes the inverse wavelet packet decomposition (IWPD) for `L` levels or by given `tree.`

!!! note
    Nonstandard transform is not yet available. Setting `standard = true` will
    currently result in an error message.

# Arguments
- `x̂::AbstractVector{T} where T<:Number` or `x̂::AbstractArray{T,2} where T<:Number`:
  Allocated output vector/matrix for reconstructed signal.
- `xw::AbstractArray{T,2} where T<:Number` or `xw::AbstractArray{T,3} where T<:Number`:
  Input array. A 2D input undergoes 1D wavelet reconstruction whereas a 3D input undergoes
  2D wavelet reconstruction.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels of wavelet
  decomposition.
- `tree::BitVector`: Binary tree or Quadtree for inverse transform to be computed
  accordingly.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `x̂::Array{T,1}` or `x̂::Array{T,2}`: Reconstructed signal.

# Examples:
```
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
wt = wavelet(WT.haar)
xw = wpd(x, wt)

# 1D wavelet reconstruction
y = iwpd(xw, wt)

# 2D wavelet decomposition
x = randn(8,8)
wt = wavelet(WT.haar)
xw = wpd(x, wt)

# 2D wavelet reconstruction
y = iwpd(xw, wt)
```

**See also:** [`iwpd!`](@ref), [`iwpt`](@ref)
"""
function iwpd(xw::AbstractArray{T},
              wt::OrthoFilter,
              L::Integer = maxtransformlevels(xw,1);
              kwargs...) where T<:Number
    @assert ndims(xw) ≥ 2
    dim = size(xw)[1:(end-1)]
    x̂ = Array{T}(undef, dim)
    iwpd!(x̂, xw, wt, L, kwargs...)
    return x̂
end

# IWPD by tree without allocation
function iwpd(xw::AbstractArray{T},
              wt::OrthoFilter,
              tree::BitVector;
              kwargs...) where T<:Number
    @assert ndims(xw) ≥ 2
    dim = size(xw)[1:(end-1)]
    x̂ = Array{T}(undef, dim)
    iwpd!(x̂, xw, wt, tree, kwargs...)
    return x̂
end

# 1D IWPD by level with allocation
"""
    iwpd(x̂, xw, wt[, L; standard])
    iwpd(x̂, xw, wt, tree[; standard])

Same as `iwpd` but without array allocation.

!!! note
    Nonstandard transform is not yet available. Setting `standard = true` will
    currently result in an error message.

# Arguments
- `x̂::AbstractVector{T} where T<:Number` or `x̂::AbstractArray{T,2} where T<:Number`:
  Allocated output vector/matrix for reconstructed signal.
- `xw::AbstractArray{T,2} where T<:Number` or `xw::AbstractArray{T,3} where T<:Number`:
  Input array. A 2D input undergoes 1D wavelet reconstruction whereas a 3D input undergoes
  2D wavelet reconstruction.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels for wavelet
  decomposition.
- `tree::BitVector`: Binary tree or Quadtree to transform to be computed accordingly.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `x̂::Array{T,1}` or `x̂::Array{T,2}`: Reconstructed signal.

# Examples:
```
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
wt = wavelet(WT.haar)
xw = wpd(x, wt)
y = similar(x)

# 1D wavelet reconstruction
iwpd!(y, xw, wt)

# 2D wavelet decomposition
x = randn(8,8)
wt = wavelet(WT.haar)
xw = wpd(x, wt)
y = similar(x)

# 2D wavelet reconstruction
iwpd!(y, xw, wt)
```

**See also:** [`iwpd`](@ref), [`iwpt`](@ref)
"""
function iwpd!(x̂::AbstractVector{T},
               xw::AbstractArray{T,2},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x̂)) where T<:Number
    iwpd!(x̂, xw, wt, maketree(length(x̂), L, :full))
    return x̂
end

# 2D IWPD by level with allocation
function iwpd!(x̂::AbstractArray{T,2},
               xw::AbstractArray{T,3},
               wt::OrthoFilter,
               L::Integer = maxtransformlevels(x̂);
               standard::Bool = true) where T<:Number
    iwpd!(x̂, xw, wt, maketree(size(x̂)..., L, :full), standard=standard)
    return x̂
end

# 1D IWPD by tree with allocation
function iwpd!(x̂::AbstractVector{T},
               xw::AbstractArray{T,2},
               wt::OrthoFilter,
               tree::BitVector) where T<:Number
    # Sanity check
    @assert size(x̂,1) == size(xw,1)
    @assert isvalidtree(x̂, tree)
    # Get basis coefficients then compute iwpt
    w = getbasiscoef(xw, tree)
    iwpt!(x̂, w, wt, tree)
    return x̂
end

# 2D IWPD by tree with allocation
function iwpd!(x̂::AbstractArray{T,2},
               xw::AbstractArray{T,3},
               wt::OrthoFilter,
               tree::BitVector;
               standard::Bool = true) where T<:Number
    # Sanity check
    @assert size(x̂,1) == size(xw,1)
    @assert size(x̂,2) == size(xw,2)
    @assert isvalidtree(x̂, tree)

    # Setup
    m, n, _ = size(xw)
    g, h = WT.makereverseqmfpair(wt, true)
    xwₜ = copy(xw)
    temp = Array{T,2}(undef, (m,n))

    # Inverse transform
    for i in reverse(eachindex(tree))
        # Reconstruct node i if it has children
        if tree[i]
            # Extract parent node
            d = getdepth(i,:quad)
            rows = getrowrange(m, i)
            cols = getcolrange(n, i)
            @inbounds v = d == 0 ? x̂ : @view xwₜ[rows, cols, d+1]
            # Extract children nodes of current parent
            dᵣ = d+1
            rows₁ = getrowrange(m, getchildindex(i,:topleft))
            cols₁ = getcolrange(n, getchildindex(i,:topleft))
            rows₂ = getrowrange(m, getchildindex(i,:topright))
            cols₂ = getcolrange(n, getchildindex(i,:topright))
            rows₃ = getrowrange(m, getchildindex(i,:bottomleft))
            cols₃ = getcolrange(n, getchildindex(i,:bottomleft))
            rows₄ = getrowrange(m, getchildindex(i,:bottomright))
            cols₄ = getcolrange(n, getchildindex(i,:bottomright))
            @inbounds w₁ = @view xwₜ[rows₁, cols₁, dᵣ+1]
            @inbounds w₂ = @view xwₜ[rows₂, cols₂, dᵣ+1]
            @inbounds w₃ = @view xwₜ[rows₃, cols₃, dᵣ+1]
            @inbounds w₄ = @view xwₜ[rows₄, cols₄, dᵣ+1]
            # Extract temp subarray (Same size as parent node)
            @inbounds tempᵢ = @view temp[rows, cols]
            # Perform 1 level of wavelet reconstruction
            @inbounds idwt_step!(v, w₁, w₂, w₃, w₄, h, g, tempᵢ, standard=standard)
        else
            continue
        end
    end
    return x̂
end

# ========== Wavelet Packet Transform ==========
# 2D Wavelet Packet Transform without allocated output array
"""
    wpt(x, wt[, L; standard])
    wpt(x, wt, tree[; standard])

Returns the wavelet packet transform (WPT) for `L` levels or by given quadratic `tree`.

!!! note
    Nonstandard transform is not yet available. Setting `standard = true` will currently
    result in an error message.

# Arguments
- `x::AbstractVector{T} where T<:Number` or `x::AbstractArray{T,2} where T<:Number`: Input
  vector/matrix. A vector input undergoes 1D wavelet decomposition whereas a matrix input
  undergoes 2D wavelet decomposition.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels for wavelet
  decomposition.
- `tree::BitVector`: Quadtree to transform to be computed accordingly.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `::Array{T,1}` or `::Array{T,2}`: Transformed signal.

# Examples:
```
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
wt = wavelet(WT.haar)
xw = wpt(x, wt)

# 2D wavelet decomposition
x = randn(8,8)
wt = wavelet(WT.haar)
xw = wpt(x, wt)
```

**See also:** [`wpt!`](@ref), [`maketree`](@ref)
"""
function Wavelets.Transforms.wpt(x::AbstractArray{T,2},
                                 wt::OrthoFilter,
                                 L::Integer = maxtransformlevels(x);
                                 standard::Bool = true) where T<:Number
    return wpt(x, wt, maketree(size(x)..., L, :full), standard=standard)
end

function Wavelets.Transforms.wpt(x::AbstractArray{T,2},
                                 wt::OrthoFilter,
                                 tree::BitVector;
                                 standard::Bool = true) where T<:Number
    y = Array{T}(undef, size(x))
    return wpt!(y, x, wt, tree, standard=standard)
end

# 2D Wavelet Packet Transform with allocated output array
"""
    wpt!(y, x, wt[, L; standard])
    wpt!(y, x, wt, tree[; standard])

Same as `wpt` but without array allocation.

!!! note
    Nonstandard transform is not yet available. Setting `standard = true` will currently
    result in an error message.

# Arguments
- `y::AbstractVector{T} where T<:Number` or `y::AbstractArray{T,2} where T<:Number`:
  Allocated output vector/matrix.
- `x::AbstractVector{T} where T<:Number` or `x::AbstractArray{T,2} where T<:Number`: Input
  vector/matrix. A vector input undergoes 1D wavelet packet transform whereas a matrix input
  undergoes 2D wavelet packet transform.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels for wavelet
  decomposition.
- `tree::BitVector`: Quadtree to transform to be computed accordingly.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `::Array{T,1}` or `::Array{T,2}`: Transformed signal.

# Examples:
```
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
xw = similar(x)
wt = wavelet(WT.haar)
wpt!(xw, x, wt)

# 2D wavelet decomposition
x = randn(8,8)
xw = similar(x)
wt = wavelet(WT.haar)
wpt!(xw, x, wt)
```

**See also:** [`wpt`](@ref), [`maketree`](@ref)
"""
function Wavelets.Transforms.wpt!(y::AbstractArray{T,2},
                                  x::AbstractArray{T,2},
                                  wt::OrthoFilter,
                                  L::Integer = maxtransformlevels(x),
                                  standard::Bool = true) where T<:Number
    return wpt!(y, x, wt, maketree(size(x)..., L, :full), standard=standard)
end

function Wavelets.Transforms.wpt!(y::AbstractArray{T,2},
                                  x::AbstractArray{T,2},
                                  wt::OrthoFilter,
                                  tree::BitVector;
                                  standard::Bool = true) where T<:Number
    # Sanity check
    @assert isvalidtree(x, tree)
    @assert size(y) == size(x)

    # ----- Allocation and setup to match Wavelets.jl's function requirements -----
    m, n = size(x)
    g, h = WT.makereverseqmfpair(wt, true)      # low & high pass filters
    yₜ = copy(x)                                 # Placeholder of y
    temp = Array{T,2}(undef, (m,n))             # Temp array
    # ----- Compute transforms based on tree -----
    for i in eachindex(tree)
        # Decompose if node i has children
        if tree[i]
            # Extract parent node
            rows = getrowrange(m, i)
            cols = getcolrange(n, i)
            @inbounds v = @view yₜ[rows, cols]
            # Extract children nodes of current parent
            rows₁ = getrowrange(m, getchildindex(i,:topleft))
            cols₁ = getcolrange(n, getchildindex(i,:topleft))
            rows₂ = getrowrange(m, getchildindex(i,:topright))
            cols₂ = getcolrange(n, getchildindex(i,:topright))
            rows₃ = getrowrange(m, getchildindex(i,:bottomleft))
            cols₃ = getcolrange(n, getchildindex(i,:bottomleft))
            rows₄ = getrowrange(m, getchildindex(i,:bottomright))
            cols₄ = getcolrange(n, getchildindex(i,:bottomright))
            @inbounds w₁ = @view y[rows₁, cols₁]
            @inbounds w₂ = @view y[rows₂, cols₂]
            @inbounds w₃ = @view y[rows₃, cols₃]
            @inbounds w₄ = @view y[rows₄, cols₄]
            # Extract temp subarray (Same size as parent node)
            @inbounds tempᵢ = @view temp[rows, cols]
            # Perform 1 level of wavelet decomposition
            @inbounds dwt_step!(w₁, w₂, w₃, w₄, v, h, g, tempᵢ, standard=standard)
            # If current range not at final iteration, copy to yₜ
            if 4*i<length(tree)
                @inbounds yₜ[rows, cols] = @view y[rows, cols]
            end
        else
            continue
        end
    end
    return y
end

# ========== Inverse Wavelet Packet Transform ==========
# 2D Inverse Wavelet Packet Transform without allocated output array
"""
    iwpt(xw, wt[, L; standard])
    iwpt(xw, wt, tree[; standard])

Returns the inverse wavelet packet transform (IWPT) for `L` levels or by given quadratic
`tree`.

!!! note
    Nonstandard transform is not yet available. Setting `standard = true` will
    currently result in an error message.

# Arguments
- `xw::AbstractVector{T} where T<:Number` or `xw::AbstractArray{T,2} where T<:Number`: Input
  vector/matrix. A vector input undergoes 1D wavelet reconstruction whereas a matrix input
  undergoes 2D wavelet reconstruction.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels for wavelet
  decomposition.
- `tree::BitVector`: Quadtree to transform to be computed accordingly.

# Keyword Arguments
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `::Array{T,1}` or `::Array{T,2}`: Reconstructed signal.

# Examples:
```
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
wt = wavelet(WT.haar)
xw = wpt(x, wt)

# 1D wavelet reconstruction
y = iwpt(xw, wt)

# 2D wavelet decomposition
x = randn(8,8)
wt = wavelet(WT.haar)
xw = wpt(x, wt)

# 2D wavelet reconstruction
y = iwpt(xw, wt)
```

**See also:** [`iwpt!`](@ref), [`wpt!`](@ref), [`maketree`](@ref)
"""
function Wavelets.Transforms.iwpt(xw::AbstractArray{T,2},
                                  wt::OrthoFilter,
                                  L::Integer = maxtransformlevels(xw);
                                  standard::Bool = true) where T<:Number
    return iwpt(xw, wt, maketree(size(xw)..., L, :full), standard=standard)
end

function Wavelets.Transforms.iwpt(xw::AbstractArray{T,2},
                                  wt::OrthoFilter,
                                  tree::BitVector;
                                  standard::Bool = true) where T<:Number
    x̂ = similar(xw)
    return iwpt!(x̂, xw, wt, tree, standard=standard)
end

# 2D Inverse Wavelet Packet Transform with allocated output array
"""
    iwpt!(x̂, xw, wt[, L; standard])
    iwpt!(x̂, xw, wt, tree[; standard])

Same as `iwpt` but without array allocation.

!!! note
    Nonstandard transform is not yet available. Setting `standard = true` will currently
    result in an error message.

# Arguments
- `x̂::AbstractVector{T} where T<:Number` or `x̂::AbstractArray{T,2} where T<:Number`:
  Allocated output vector/matrix.
- `xw::AbstractVector{T} where T<:Number` or `xw::AbstractArray{T,2} where T<:Number`: Input
  vector/matrix. A vector input undergoes 1D wavelet packet reconstruction whereas a matrix
  input undergoes 2D wavelet packet reconstruction.
- `wt::OrthoFilter`: Wavelet filter.
- `L::Integer`: (Default: `maxtransformlevels(x)`) Number of levels for wavelet
  decomposition.
- `tree::BitVector`: Quadtree to transform to be computed accordingly.
- `standard::Bool`: (Default: `true`) Whether to compute the standard or non-standard
  wavelet transform. Only applicable for 2D signals.

# Returns
- `::Array{T,1}` or `::Array{T,2}`: Reconstructed signal.

# Examples:
```
using Wavelets, WaveletsExt

# 1D wavelet decomposition
x = randn(8)
xw = similar(x)
wt = wavelet(WT.haar)
wpt!(xw, x, wt)

# 1D wavelet reconstruction
y = similar(x)
iwpt!(y, xw, wt)

# 2D wavelet decomposition
x = randn(8,8)
xw = similar(x)
wt = wavelet(WT.haar)
wpt!(xw, x, wt)

# 2D wavelet reconstruction
y = similar(x)
iwpt!(y, xw, wt)
```

**See also:** [`iwpt`](@ref), [`wpt!`](@ref) [`maketree`](@ref)
"""
function Wavelets.Transforms.iwpt!(x̂::AbstractArray{T,2},
                                   xw::AbstractArray{T,2},
                                   wt::OrthoFilter,
                                   L::Integer = maxtransformlevels(xw);
                                   standard::Bool = true) where T<:Number
    return iwpt!(x̂, xw, wt, maketree(size(xw)..., L, :full), standard=standard)
end

function Wavelets.Transforms.iwpt!(x̂::AbstractArray{T,2},
                                   xw::AbstractArray{T,2},
                                   wt::OrthoFilter,
                                   tree::BitVector;
                                   standard::Bool = true) where T<:Number
    # Sanity check
    @assert isvalidtree(xw, tree)
    @assert size(x̂) == size(xw)

    # ----- Setup -----
    m, n = size(xw)
    g, h = WT.makereverseqmfpair(wt, true)      # low & high pass filters
    xwₜ = copy(xw)                               # Placeholder for xw
    temp = Array{T,2}(undef, (m,n))             # temp array
    # ----- Compute transforms based on tree -----
    for i in reverse(eachindex(tree))
        # Reconstruct to node i if it has children
        if tree[i]
            # Extract parent node
            rows = getrowrange(m, i)
            cols = getcolrange(n, i)
            @inbounds v = @view x̂[rows, cols]
            # Extract children nodes of current parent
            rows₁ = getrowrange(m, getchildindex(i,:topleft))
            cols₁ = getcolrange(n, getchildindex(i,:topleft))
            rows₂ = getrowrange(m, getchildindex(i,:topright))
            cols₂ = getcolrange(n, getchildindex(i,:topright))
            rows₃ = getrowrange(m, getchildindex(i,:bottomleft))
            cols₃ = getcolrange(n, getchildindex(i,:bottomleft))
            rows₄ = getrowrange(m, getchildindex(i,:bottomright))
            cols₄ = getcolrange(n, getchildindex(i,:bottomright))
            @inbounds w₁ = @view xwₜ[rows₁, cols₁]
            @inbounds w₂ = @view xwₜ[rows₂, cols₂]
            @inbounds w₃ = @view xwₜ[rows₃, cols₃]
            @inbounds w₄ = @view xwₜ[rows₄, cols₄]
            # Extract temp subarray (Same size as parent node)
            @inbounds tempᵢ = @view temp[rows, cols]
            # Perform 1 level of wavelet reconstruction
            @inbounds idwt_step!(v, w₁, w₂, w₃, w₄, h, g, tempᵢ, standard=standard)
            # If current range not at final iteration, copy to temp
            if i > 1
                @inbounds xwₜ[rows, cols] = v
            end
        else
            continue
        end
    end
    return x̂
end

include("dwt_one_level.jl")
include("dwt_all.jl")

end # end module


text:
using CausalInference
using LightGraphs
using Test

include("plotdag.jl")

g1 = g = DiGraph(7)
d = nv(g)
for (i,j) in [(1,2), (2,3), (2,4), (4,5), (3,5), (5,6), (7,5)]
    add_edge!(g, i, j)
end

tp = plot_dag(g)
save(PDF("dag1"), tp)

@testset "dsep g1" begin
    @test !dsep(g, 1, 2, [])
    @test !dsep(g, 2, 5, [])

    @test !dsep(g, 1, 5, [3])
    @test !dsep(g, 1, 5, [4])
    @test dsep(g, 1, 5, [3, 4])

    @test dsep(g, 3, 4, [2])
    @test dsep(g, 3, 4, [2, 7])

    @test !dsep(g, 3, 4, [5])
    @test !dsep(g, 3, 4, [2, 6])

    @test !dsep(g, 3, 4, [2, 5])
    @test !dsep(g, 3, 4, [2, 6])

    @test !dsep(g, 3, 5, [6])
    @test !dsep(g, 3, 5, [7])

end

g2 = g = DiGraph(7)
d = nv(g)
for (i,j) in [(1,3), (2,3), (3,4),(3,5), (4,6), (6, 7)]
    add_edge!(g, i, j)
end

tp = plot_dag(g)
save(PDF("dag2"), tp)

@testset "dsep g2" begin
    @test !dsep(g, 1, 2, [4, 5])
    @test dsep(g, 1, 2, [])
    @test !dsep(g, 1, 2, [3])
    @test dsep(g, 4, 5, [3])
    @test !dsep(g, 4, 5, [])
    @test !dsep(g, 4, 5, [1, 2])
end


text:
using MultivariateFunctions

tol = 10*eps()

# Testing constructors.
before = PE_Function(0.0,0.0,0.0,1)
first_ = PE_Function(1.0,0.0,0.0,1)
secc = PE_Function(20.0,0.0,0.0,2)
last_  = Sum_Of_Functions([first_, secc])

f1 = Piecewise_Function([before, first_, secc, last_], [-Inf, -1.0,3.0, 10.0])
abs(evaluate(f1, -10.0) - 0.0 ) < tol
abs(evaluate(f1,   1.0) - 1.0 ) < tol
abs(evaluate(f1,   5.0) - 500.0) < tol
abs(evaluate(f1,  20.0) - 20.0 - 20*20*20) < tol
f2 = Piecewise_Function([before, first_, secc, last_], [-0.1, 0.0,2.0, 40.0])

f3 =  Piecewise_Function([before, f1, secc, last_], [-0.1, 0.0,2.0, 40.0])
f4 =  Piecewise_Function([secc, f1, secc, last_], [-2.1, -1.1,4.0, 40.0])

rebadge_test = rebadge(f3, Dict{Symbol,Symbol}(:default => :tester))
pop!(underlying_dimensions(rebadge_test)) == :tester
function test_result(func, eval0, eval5, len = 1)
    val_test0 = abs(evaluate(func, 0.0) - eval0) < 1e-09
    if (!val_test0)
        print("Failed Val Test at 0.0")
    end
    val_test5 = abs(evaluate(func, 5.0) - eval5) < 1e-09
    if (!val_test5)
        print("Failed Val Test at 5.0")
    end
    if typeof(func) == MultivariateFunctions.Piecewise_Function
        len_test = length(func.functions_) == len
        if (!len_test)
            print("Failed Length Test")
        end
        return all([val_test0,val_test5,len_test])
    else
        return all([val_test0,val_test5])
    end
end

# Testing with Ints
test_result(f1 + 5, 5.0, 505.0, 4)
test_result(5 + f1, 5.0, 505.0, 4)
test_result(f1 - 5, -5.0, 495.0, 4)
test_result(5 - f1, 5.0, -495.0, 4)
test_result(f1 * 5, 0.0, 2500.0, 4)
test_result(5 * f1, 0.0, 2500.0, 4)
test_result(f1 / 5, 0.0, 100.0, 4)
# And with Float64s
test_result(f1 + 5.0, 5.0, 505.0, 4)
test_result(5.0 + f1, 5.0, 505.0, 4)
test_result(f1 - 5.0, -5.0, 495.0, 4)
test_result(5.0 - f1, 5.0, -495.0, 4)
test_result(f1 * 5.0, 0.0, 2500.0, 4)
test_result(5.0 * f1, 0.0, 2500.0, 4)
test_result(f1 / 5.0, 0.0, 100.0, 4)

typeof(f4 + Missing()) == Missing
typeof(Missing() + f4) == Missing
typeof(f4 - Missing()) == Missing
typeof(Missing() - f4) == Missing
typeof(f4 * Missing()) == Missing
typeof(Missing() * f4) == Missing

typeof(-1*f4) == MultivariateFunctions.Piecewise_Function

test_result(f1 + first_, 0.0, 505.0, 4)
test_result(first_ + f1, 0.0, 505.0, 4)
test_result(f1 - first_, 0.0, 495.0, 4)
test_result(first_ - f1, 0.0, -495.0, 4)
test_result(f1 * first_, 0.0, 2500.0, 4)
test_result(first_ * f1, 0.0, 2500.0, 4)

test_result(f1 + last_, 0.0, 1005.0, 4)
test_result(last_ + f1, 0.0, 1005.0, 4)
test_result(f1 - last_, 0.0, -5.0, 4)
test_result(last_ - f1, 0.0, 5.0, 4)
test_result(f1 * last_, 0.0, 252500.0, 4)
test_result(last_ * f1, 0.0, 252500.0, 4)

test_result(f1 + f4, 0.0, 1000.0, 8)
test_result(f1 - f4, 0.0, 0.0, 8)
test_result(f1 * f4, 0.0, 250000.0, 8)

# Testing linear rescaling
#test function will be f4
alpha = 0.95
beta = 3.0
X = [-0.3, 0.0, 1.3, 1.5, 4.0, 7.0, 20.0, 40.0, 70.0, 100.0]
X_converted = (alpha .* X) .+ beta
f4_converted = convert_to_linearly_rescale_inputs(f4, alpha, beta)
y = evaluate.(Ref(f4),X)
converted_y = evaluate.(Ref(f4_converted), X_converted)
sum(abs.(y .- converted_y)) < 1e-10


text:
using Test
using Random: seed!
using StaticArrays: SVector
using QuantumLattices.Essentials.Spatials
using QuantumLattices.Interfaces: decompose, rank, dimension, kind, expand, reset!
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Prerequisites.SimpleTrees: leaves

@testset "distance" begin
    @test distance([0.0, 0.0], [1.0, 1.0]) ≈ sqrt(2.0)
end

@testset "azimuthd" begin
    @test azimuthd([+1.0]) ≈ 0.0
    @test azimuthd([-1.0]) ≈ 180.0
    @test azimuthd([+1.0, +1.0]) ≈ 45.0
    @test azimuthd([-1.0, -1.0]) ≈ 225.0
    @test azimuthd([+1.0, +1.0, -2.0]) ≈ 45.0
    @test azimuthd([+1.0, -1.0, +2.0]) ≈ 315.0
end

@testset "azimuth" begin
    @test azimuth([+1.0]) ≈ 0.0
    @test azimuth([-1.0]) ≈ pi
    @test azimuth([+1.0, +1.0]) ≈ 1//4*pi
    @test azimuth([-1.0, -1.0]) ≈ 5//4*pi
    @test azimuth([+1.0, +1.0, -2.0]) ≈ 1//4*pi
    @test azimuth([+1.0, -1.0, +2.0]) ≈ 7//4*pi
end

@testset "polard" begin
    @test polard([1.0, 0.0, 1.0]) ≈ 45.0
end

@testset "polard" begin
    @test polar([1.0, 0.0, 1.0]) ≈ 1//4*pi
end

@testset "volume" begin
    @test volume([1, 0], [0, 1], [1, 1]) == 0
    @test volume([1, 0, 0], [0, 1, 0], [0, 0, 1]) == 1
end

@testset "isparallel" begin
    @test isparallel([0.0, 0.0], [1.0, 1.0]) == 1
    @test isparallel([2.0, 2.0], [1.0, 1.0]) == 1
    @test isparallel([-2.0, -2.0], [1.0, 1.0]) == -1
    @test isparallel([1.0, -1.0], [1.0, 1.0]) == 0
end

@testset "isonline" begin
    seed!()
    p1, p2 = [0.0, 0.0], [1.0, 1.0]
    @test isonline([0.0, 0.0], p1, p2, ends = (true, rand(Bool))) == true
    @test isonline([0.0, 0.0], p1, p2, ends = (false, rand(Bool))) == false
    @test isonline([1.0, 1.0], p1, p2, ends = (rand(Bool), true)) == true
    @test isonline([1.0, 1.0], p1, p2, ends = (rand(Bool), false)) == false
    @test isonline([0.5, 0.5], p1, p2, ends = (rand(Bool), rand(Bool))) == true
    @test isonline([1.1, 1.1], p1, p2, ends = (rand(Bool), rand(Bool))) == false
end

@testset "decompose" begin
    a, c = randn(3), randn()
    @test all(decompose(c*a, a) .≈ (c,))

    a1, a2 = randn(2), randn(2)
    c1, c2 = randn(2)
    @test all(decompose(c1*a1+c2*a2, a1, a2) .≈ (c1, c2))

    a1, a2 = randn(3), randn(3)
    c1, c2 = randn(2)
    @test all(decompose(c1*a1+c2*a2, a1, a2) .≈ (c1, c2))

    a1, a2, a3 = randn(3), randn(3), randn(3)
    c1, c2, c3 = randn(3)
    @test all(decompose(c1*a1+c2*a2+c3*a3, a1, a2, a3) .≈ (c1, c2, c3))
end

@testset "isintratriangle" begin
    seed!()
    p1, p2, p3 = [-1.0, 0.0], [0.0, 1.0], [1.0, 0.0]
    for (i, p) in enumerate((p1, p2, p3))
        @test isintratriangle(p, p1, p2, p3, vertexes=Tuple(j == i ? true : rand(Bool) for j = 1:3), edges=(rand(Bool), rand(Bool), rand(Bool))) == true
        @test isintratriangle(p, p1, p2, p3, vertexes=Tuple(j == i ? false : rand(Bool) for j = 1:3), edges=(rand(Bool), rand(Bool), rand(Bool))) == false
    end
    for (i, p) in enumerate(((p1+p2)/2, (p2+p3)/2, (p3+p1)/2))
        @test isintratriangle(p, p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=Tuple(j == i ? true : rand(Bool) for j = 1:3)) == true
        @test isintratriangle(p, p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=Tuple(j == i ? false : rand(Bool) for j = 1:3)) == false
    end
    @assert isintratriangle([0.0, 0.5], p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=(rand(Bool), rand(Bool), rand(Bool))) == true
    @assert isintratriangle([0.0, 1.5], p1, p2, p3, vertexes=(rand(Bool), rand(Bool), rand(Bool)), edges=(rand(Bool), rand(Bool), rand(Bool))) == false
end

@testset "issubordinate" begin
    @test issubordinate([2.0], [[1.0]]) == true
    @test issubordinate([1.5], [[1.0]]) == false
    @test issubordinate([0.5, 3.5], [[0.5, 0.5], [-0.5, 0.5]]) == true
    @test issubordinate([0.5, 3.0], [[0.5, 0.5], [-0.5, 0.5]]) == false
    @test issubordinate([4.0, 4.0, 4.0], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) == true
    @test issubordinate([4.0, 4.0, 2.5], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) == false
end

@testset "reciprocals" begin
    @test reciprocals([[1.0]]) ≈ [[1.0]]*2pi
    @test reciprocals([[1.0, 0.0], [0.0, 1.0]]) ≈ [[1.0, 0.0], [0.0, 1.0]]*2pi
    @test reciprocals([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]) ≈ [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]*2pi
end

@testset "translate" begin
    @test translate([1.0 2.0; 3.0 4.0], [1.0, 2.0]) == [2.0 3.0; 5.0 6.0]
end

@testset "rotate" begin
    @test rotate([1.0 2.0; 1.0 2.0], pi/4, axis=([1.0, 1.0], (0.0, 0.0))) ≈ [1.0 1.0; 1.0 1.0+√2.0]
    @test rotate(reshape([2.0, 2.0, 2.0], 3, 1), pi/12, axis=([0.0, 1.0, 1.0], (pi/2, 0.0))) ≈ reshape([2.0, 1.0+√2/2, 1.0+√6/2], 3, 1)
    @test rotate(reshape([2.0, 2.0, 2.0], 3, 1), pi/12, axis=([1.0, 0.0, 1.0], (pi/2, pi/2))) ≈ reshape([1.0+√6/2, 2.0, 1.0+√2/2], 3, 1)
    @test rotate(reshape([2.0, 2.0, 2.0], 3, 1), pi/12, axis=([1.0, 1.0, 0.0], (0.0, 0.0))) ≈ reshape([1.0+√2/2, 1.0+√6/2, 2.0], 3, 1)
end

@testset "tile" begin
    @test tile(reshape([0.0, 0.0], 2, 1), [[1.0, 0.0], [0.0, 1.0]], ((0.5, 0.5), (-0.5, -0.5))) == [0.5 -0.5; 0.5 -0.5]
end

@testset "minimumlengths" begin
    @test minimumlengths(reshape([0.0, 0.0], 2, 1), [[1.0, 0.0], [0.0, 1.0]], 7) ≈ [1.0, √2, 2.0, √5, 2*√2, 3.0, √10]
end

@testset "intralinks" begin
    ps, a1, a2 = [0.0 0.5; 0.0 0.5], [1.0, 0.0], [0.0, 1.0]
    neighbors = Dict(1=>1.0)
    links = [(1, 1, 1, SVector(0.0, -1.0)), (1, 1, 1, SVector(-1.0, 0.0)), (1, 2, 2, SVector(0.0, -1.0)), (1, 2, 2, SVector(-1.0, 0.0))]
    @test intralinks(ps, [a1, a2], neighbors) == links
end

@testset "interlinks" begin
    ps1 = [0.0 1.0; 0.0 0.0]
    ps2 = [0.0 1.0; 1.0 1.0]
    neighbors = Dict(1=>1.0)
    links = [(1, 1, 1), (1, 2, 2)]
    @test interlinks(ps1, ps2, neighbors) == links
end

@testset "PID" begin
    @test PID(1) == PID('T', 1) == PID(scope='T', site=1)
    @test PID(scope="tz", site=1)|>string == "PID(\"tz\", 1)"
    @test fieldnames(PID{String}) == (:scope, :site)
end

@testset "Point" begin
    point = Point(PID(0, 1), (0.0, 0.0), (0.0, 0.0))
    @test point == Point(PID(0, 1), [0.0, 0.0], [0.0, 0.0])
    @test isequal(point, Point{2}(PID(0, 1), [0.0, 0.0], [0.0, 0.0]))
    @test point|>length == point|>typeof|>length == 1
    @test point|>eltype == point|>typeof|>eltype == Point{2, PID{Int}}
    @test point|>rank == point|>typeof|>rank == 1
    @test point|>pidtype == point|>typeof|>pidtype == PID{Int}
    @test point|>dimension == point|>typeof|>dimension == 2
    @test point|>kind == point|>typeof|>kind == 0
    @test point|>string == "Point(PID(0, 1), [0.0, 0.0], [0.0, 0.0])"
    @test point|>collect == [point]
end

@testset "Bond" begin
    bond = Bond(1, Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0)), Point(PID(1, 2), (0.0, 1.0), (0.0, 1.0)))
    @test bond|>deepcopy == bond
    @test isequal(bond|>deepcopy, bond)
    @test bond|>string == "Bond(1, Point(PID(1, 1), [0.0, 0.0], [0.0, 0.0]), Point(PID(1, 2), [0.0, 1.0], [0.0, 1.0]))"
    @test bond|>reverse == Bond(1, Point(PID(1, 2), (0.0, 1.0), (0.0, 1.0)), Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0)))
    @test bond|>length == bond|>typeof|>length == 2
    @test bond|>eltype == bond|>typeof|>eltype == Point{2, PID{Int}}
    @test bond|>rank == bond|>typeof|>rank == 2
    @test bond|>pidtype == bond|>typeof|>pidtype == PID{Int}
    @test bond|>dimension == bond|>typeof|>dimension == 2
    @test bond|>kind == 1
    @test bond|>rcoord == [0.0, 1.0]
    @test bond|>icoord == [0.0, 1.0]
    @test bond|>isintracell == false
    @test bond|>collect == [Point(PID(1, 2), (0.0, 1.0), (0.0, 1.0)), Point(PID(1, 1), (0.0, 0.0), (0.0, 0.0))]
end

@testset "Lattice" begin
    lattice = Lattice("Tuanzi", [Point(PID(1, 1), (0.5, 0.5), (0.0, 0.0))], vectors=[[1.0, 0.0], [0.0, 1.0]], neighbors=1)
    @test lattice|>deepcopy == lattice
    @test isequal(lattice|>deepcopy, lattice)
    @test lattice|>string == "Lattice(Tuanzi)\n  with 1 point:\n    Point(PID(1, 1), [0.5, 0.5], [0.0, 0.0])\n  with 2 translation vectors:\n    [1.0, 0.0]\n    [0.0, 1.0]\n  with 1 order of nearest neighbors:\n    1 => 1.0\n"
    @test lattice|>length == 1
    @test lattice|>dimension == lattice|>typeof|>dimension == 2
    @test lattice|>keytype == lattice|>typeof|>keytype == PID{Int}
    @test lattice|>valtype == lattice|>typeof|>valtype == Point{2, PID{Int}}
    @test lattice[LatticeIndex{'R'}(1)] == lattice[LatticeIndex{'R'}(PID(1, 1))] == SVector(0.5, 0.5)
    @test lattice[LatticeIndex{'I'}(1)] == lattice[LatticeIndex{'I'}(PID(1, 1))] == SVector(0.0, 0.0)
    @test lattice[LatticeIndex{'P'}(1)] == lattice[LatticeIndex{'P'}(PID(1, 1))] == Point(PID(1, 1), (0.5, 0.5), (0.0, 0.0))
    @test lattice|>nneighbor == 1

    zerothbs = [Point(PID(1, 1), (0.5, 0.5), (0.0, 0.0))]
    insidebs = []
    acrossbs = [Bond(1, Point(PID(1, 1), [0.5, 0.5], [0.0, 0.0]), Point(PID(1, 1), [0.5, -0.5], [0.0, -1.0])),
                Bond(1, Point(PID(1, 1), [0.5, 0.5], [0.0, 0.0]), Point(PID(1, 1), [-0.5, 0.5], [-1.0, 0.0]))
    ]
    @test bonds(lattice, zerothbonds) == zerothbs
    @test bonds(lattice, insidebonds) == insidebs
    @test bonds(lattice, acrossbonds) == acrossbs
    @test setdiff(bonds(lattice), [zerothbs; insidebs; acrossbs])|>length == 0

    lattice = Lattice(  "SuperTuanzi",
                        [   Lattice("Tuanzi1", [Point(PID(1, 1), (0.0, 0.0))], neighbors=0),
                            Lattice("Tuanzi2", [Point(PID(2, 1), (0.5, 0.5))], neighbors=0),
                            ],
                        vectors=    [[1.0, 0.0], [0.0, 1.0]],
                        neighbors=  2,
                        )
    @test setdiff(bonds(lattice), bonds(lattice, zerothbonds, insidebonds, acrossbonds))|>length == 0
end

@testset "SuperLattice" begin
    lattice = SuperLattice( "SuperTuanzi",
                            [   Lattice("TuanziSys", [Point(PID(1, 1), (0.0, 0.0)), Point(PID(1, 2), (0.5, 0.5))], neighbors=Dict(1=>√2/2)),
                                Lattice("TuanziEnv", [Point(PID(2, 1), (-0.05, -0.05)), Point(PID(2, 2), (0.55, 0.55))], neighbors=Dict{Int, Float}()),
                                ],
                            vectors=    [[1.0, 0.0], [0.0, 1.0]],
                            neighbors=  Dict(1=>√2/2, -1=>√2/20)
                            )
    @test latticetype(lattice) == latticetype(typeof(lattice)) == Lattice{2, PID{Int}}
    @test setdiff(bonds(lattice), bonds(lattice, zerothbonds, insidebonds, acrossbonds))|>length == 0
    @test setdiff(bonds(lattice, insidebonds), bonds(lattice, intrabonds, interbonds))|>length == 0
end

@testset "Cylinder" begin
    cylinder = Cylinder{PID{String}}("Tuanzi", [0.0 0.0; 0.0 1.0], SVector(1.0, 0.0), vector=[0.0, 2.0], neighbors=1)
    insert!(cylinder, "A", "B")
    insert!(cylinder, "C3", cut=2, scopes=["C1", "C1", "C2", "C2"])
    @test cylinder.pids == [PID("C1", 1), PID("C3", 2), PID("C3", 1), PID("C1", 2), PID("C2", 1), PID("C2", 2)]
    @test cylinder.rcoords == [-1.0 -1.0 0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0 0.0 1.0]
    @test cylinder.icoords == [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0]
    lattice = cylinder("C1", "C2", "C3")
    @test lattice.pids == [PID("C1", 1), PID("C1", 2), PID("C2", 1), PID("C2", 2), PID("C3", 1), PID("C3", 2)]
    @test lattice.rcoords == [-1.0 -1.0 0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0 0.0 1.0]
    @test lattice.icoords == [0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0]
end

@testset "LatticeBonds" begin
    @test eltype(typeof(allbonds)) === AbstractBond
    @test eltype(typeof(zerothbonds)) === Point
    @test eltype(typeof(insidebonds)) === eltype(typeof(acrossbonds)) === Bond
    @test eltype(typeof(intrabonds)) === eltype(typeof(interbonds)) === Bond
    @test eltype(AbstractLattice{2, PID{Int}}, allbonds) == eltype(AbstractLattice{2, PID{Int}}, Val(allbonds)) == AbstractBond{2, PID{Int}}
    @test eltype(Lattice{3, PID{Char}}, zerothbonds) == eltype(Lattice{3, PID{Char}}, Val(zerothbonds)) == Point{3, PID{Char}}
    @test eltype(Lattice{1, PID{String}}, insidebonds) == eltype(Lattice{1, PID{String}}, Val(insidebonds)) == Bond{1, PID{String}}

    @test expand(Lattice, Val((allbonds,))) == Lattice|>latticebondsstructure|>leaves|>Tuple
    @test expand(Lattice, Val((zerothbonds,))) == (zerothbonds,)
    @test expand(Lattice, Val((insidebonds,))) == (insidebonds,)
    @test expand(Lattice, Val((acrossbonds,))) == (acrossbonds,)

    @test expand(SuperLattice, Val((allbonds,))) == SuperLattice|>latticebondsstructure|>leaves|>Tuple
    @test expand(SuperLattice, Val((zerothbonds,))) == (zerothbonds,)
    @test expand(SuperLattice, Val((insidebonds,))) == (intrabonds, interbonds)
    @test expand(SuperLattice, Val((acrossbonds,))) == (acrossbonds,)
    @test expand(SuperLattice, Val((intrabonds,))) == (intrabonds,)
    @test expand(SuperLattice, Val((interbonds,))) == (interbonds,)
end

@testset "Bonds" begin
    lattice = Lattice("Tuanzi", [Point(PID(1, 1), (0.5, 0.5), (0.0, 0.0))], vectors=[[1.0, 0.0], [0.0, 1.0]], neighbors=1)
    bs = Bonds(lattice)
    @test bs|>eltype == bs|>typeof|>eltype == AbstractBond{2, PID{Int}}
    @test bs==Bonds(lattice, allbonds) && isequal(bs, Bonds(lattice, allbonds))
    @test (bs|>length == 3) && (bs|>size == (3,))
    @test summary(bs) == "3-element Bonds"
    @test bs|>bondtypes == bs|>typeof|>bondtypes == (zerothbonds, insidebonds, acrossbonds)
    @test bs|>latticetype == bs|>typeof|>latticetype == Lattice{2, PID{Int}}
    @test bs|>rank == bs|>typeof|>rank == 3
    @test bs|>collect == bonds(lattice)
    @test (bs[1] == bs.bonds[1][1]) && (bs[2] == bs.bonds[3][1]) && (bs[3] == bs.bonds[3][2])
    @test filter(allbonds, bs, :include|>Val) == filter((zerothbonds, insidebonds, acrossbonds), bs, :include|>Val) == bs
    @test filter(zerothbonds, bs, :include|>Val) == filter((insidebonds, acrossbonds), bs, :exclude|>Val) == Bonds(lattice, zerothbonds)
    @test filter(insidebonds, bs, :include|>Val) == filter((zerothbonds, acrossbonds), bs, :exclude|>Val) == Bonds(lattice, insidebonds)
    @test filter(acrossbonds, bs, :include|>Val) == filter((zerothbonds, insidebonds), bs, :exclude|>Val) == Bonds(lattice, acrossbonds)
    emptybs = Bonds{(zerothbonds, insidebonds, acrossbonds), Lattice{2, PID{Int}}}((Point{2, PID{Int}}[], Bond{2, PID{Int}}[], Bond{2, PID{Int}}[]))
    @test filter(bond->(dimension(bond) == 3), bs) == empty!(deepcopy(bs)) == empty(bs) == emptybs
    @test reset!(emptybs, lattice) == bs
end


text:
module RocketTapOnUnsubscribeOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: tap_on_unsubscribe()" begin

    println("Testing: operator tap_on_unsubscribe()")

    run_proxyshowcheck("TapOnUnsubscribe", tap_on_unsubscribe(() -> begin end), args = (check_subscription = true, ))

    @testset begin
        values = Int[]
        source = from(1:5) |> tap_on_unsubscribe(() -> push!(values, -1))
        actor  = keep(Int)

        @test values == []

        subscription = subscribe!(source, actor)

        @test actor.values == [ 1, 2, 3, 4, 5 ]
        @test values == [ ]

        unsubscribe!(subscription)
        @test values == [ -1 ]
        unsubscribe!(subscription)
        @test values == [ -1 ]

        subscription = subscribe!(source, actor)

        @test actor.values == [ 1, 2, 3, 4, 5, 1, 2, 3, 4, 5 ]
        @test values       == [ -1 ]

        unsubscribe!(subscription)
        @test values       == [ -1, -1 ]
        unsubscribe!(subscription)
        @test values       == [ -1, -1 ]
    end

    struct DummySource <: Subscribable{Any}
        values :: Vector{Any}
    end

    function Rocket.on_subscribe!(source::DummySource, actor)
        next!(actor, 1)
        return DummySubscription(source.values)
    end

    struct DummySubscription <: Teardown
        values :: Vector{Any}
    end

    Rocket.as_teardown(::Type{ <: DummySubscription }) = UnsubscribableTeardownLogic()

    function Rocket.on_unsubscribe!(subscription::DummySubscription)
        push!(subscription.values, 0)
        return nothing
    end

    @testset begin
        values = []
        source = DummySource(values) |> tap_on_unsubscribe(() -> push!(values, -1), TapBeforeUnsubscription())

        @test values == []

        subscription = subscribe!(source, (e) -> push!(values, e))

        @test values == [ 1 ]

        unsubscribe!(subscription)

        @test values == [ 1, -1, 0 ]

        unsubscribe!(subscription)

        @test values == [ 1, -1, 0, 0 ]

        subscription = subscribe!(source, (e) -> push!(values, e))

        @test values == [ 1, -1, 0, 0, 1 ]

        unsubscribe!(subscription)

        @test values == [ 1, -1, 0, 0, 1, -1, 0 ]

        unsubscribe!(subscription)

        @test values == [ 1, -1, 0, 0, 1, -1, 0, 0 ]
    end

    @testset begin
        values = []
        source = DummySource(values) |> tap_on_unsubscribe(() -> push!(values, -1), TapAfterUnsubscription())

        @test values == []

        subscription = subscribe!(source, (e) -> push!(values, e))

        @test values == [ 1 ]

        unsubscribe!(subscription)

        @test values == [ 1, 0, -1 ]

        unsubscribe!(subscription)

        @test values == [ 1, 0, -1, 0 ]

        subscription = subscribe!(source, (e) -> push!(values, e))

        @test values == [ 1, 0, -1, 0, 1 ]

        unsubscribe!(subscription)

        @test values == [ 1, 0, -1, 0, 1, 0, -1 ]

        unsubscribe!(subscription)

        @test values == [ 1, 0, -1, 0, 1, 0, -1, 0 ]
    end

    @testset begin
        values = Int[]
        source = completed(Int) |> tap_on_unsubscribe(() -> push!(values, -1))
        actor  = keep(Int)

        @test values == []

        subscription = subscribe!(source, actor)

        @test actor.values == [ ]
        @test values       == [ ]

        unsubscribe!(subscription)

        @test actor.values == [ ]
        @test values       == [ -1 ]

        unsubscribe!(subscription)

        @test actor.values == [ ]
        @test values       == [ -1 ]
    end

end

end


text:
# include -
include("Include.jl")

# load my experimental_data file -
path_to_experimental_data_file = joinpath(_PATH_TO_DATA, "Data-3-5-AMP.txt")
experimental_data_array = load_experimental_data_set(path_to_experimental_data_file)

# load my parameters -> configuration_dictionary contains, two sub-dictionaries
path_to_configuration_file = joinpath(_PATH_TO_CONFIG, "Configuration.toml")
configuration_dictionary = load_configuration_dictionary(path_to_configuration_file)
v_variable_parameter_dictionary = configuration_dictionary["vvariableparameters"]
kinetic_limit_parameter_dictionary = configuration_dictionary["kineticparameters"]

# lets compute the kinetic limit -
kinetic_limit_value = compute_kinetic_limit(kinetic_limit_parameter_dictionary)

# ok, so let's compute a data driven v-variable -
data_driven_v_variable = Array{Float64,1}()
(number_of_effector_values, number_of_cols) = size(experimental_data_array)
for row_index = 1:number_of_effector_values
    value = experimental_data_array[row_index,2]
    scaled_value = value / kinetic_limit_value
    push!(data_driven_v_variable, scaled_value)
end

# make a plot -
plot(experimental_data_array[:,1],data_driven_v_variable[:,1],legend=false)
xlabel!("Effector 3-5-AMP (mM)",fontstze=16)
ylabel!("v-variable (AU)", fontsize=16)


text:
using Test

using FluidSynth

settings = Settings()
setstr(settings, "audio.file.name", "test_result.wav")
setstr(settings, "player.timing-source", "sample")
setint(settings, "synth.lock-memory", Int32(0));

synth = Synth(settings)
sfont_id = sfload(synth, "violin_sample.sf2")
program_select(synth, Int32(0), sfont_id, Int32(0), Int32(0))

player = Player(synth)
add(player, "test.mid")
play(player)

file_renderer = FileRenderer(synth)
while get_status(player) == Playing
    if process_block(file_renderer) != FLUID_OK
        break
    end
end

stop(player)
FluidSynth.join(player)


text:
module HatreeFock

if haskey(ENV,"LEARNHATREEFOCK_PATH")
    const SRC_PATH=ENV["LEARNHATREEFOCK_PATH"];
else
    println("! Environmental variable LEARNHARTREEFOCK_PATH not set, defaulting to pwd() !")
    const SRC_PATH=pwd();
end

#using CodeInfo
include("ModuleList.jl")

"""
                        HatreeFock(filename=user_filename)

Description: This is the main function responsible for executing the Hatree-Fock
routines for calculating the total energy of a given system.

Program Outline
1.  Parse atomic system info:
             a.) parse(filename) returns a datatype with atomic number and coordinates.
2.  Caluclate nuclear-nuclear repulsion
3.  Setup basis functions:
            a.) returns the basis-set functions and number of electrons
            b.) build orbital overlap given to handle non-ortho. Gaussian basis set
4.   Calculate kinetic energy of each electron
5.   Build/calculate electron-nuclear attraction
6.   Build/calculate electron-electron repulsion (mean-field and exchange)
7.   Self-consistent field solution for total energy

""" function main(filename="HatreeFock.in")

    #printcodeinfo()

    atomicsystem, basis = getcalcsetup(filename);
    natoms = atomicsystem.natoms;

    ZZ = calcnuclrepul(atomicsystem);

    basisfunc, numelec = buildbasisfunc(atomicsystem,basis);

    S = buildelecoverlap(natoms,basisfunc,basis);

    KE = buildkineticenergy(natoms,basisfunc,basis);

    Zq = buildnuclrattract(atomicsystem,basisfunc,basis);

    qq = buildelecelecrepulsion(natoms,basisfunc,basis);

    size(KE) == size(Zq) ? nothing : throw(AssertionError("Kinetic energy matrix size does not equal nuclear attraction matrix!"))

    Ho = KE + Zq ; #TODO: Need unit test

    energyscf = runscf(numelec,S,Ho,qq);

    E = energyscf + ZZ;

    println("\n")
    println("---------------------------------------------")
    println("Total Minimum Energy: $E [Ha]")
    println("---------------------------------------------")
    println("\n")

    return println("!!LearnHatreeFock.jl Succesfully Ran!!\n")
end #main

#Run Program
if !isempty(ARGS)
    main(ARGS[1])
else
    println("LearnHatreeFock.jl Module loaded!")
end
#End Program

end


text:
using Test
using GeometricIntegrators
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Utils
using GeometricProblems.PointVortices

const Δt = 0.01
const nt = 1000

const ref = [0.18722529318641928, 0.38967432450068706, 0.38125332930294187, 0.4258020604293123]

@testset "$(rpad("Point Vortices",80))" begin
    ode  = point_vortices_ode()
    iode = point_vortices_iode()
    idae = point_vortices_idae()

    int = Integrator(ode, TableauGauss(1), Δt)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, ref) < 3E-2

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(1), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, ref) < 3E-2

    int = Integrator(idae, TableauVSPARKGLRKpSymmetric(1), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, ref) < 3E-2
end


text:
# This file contains routine for bounds checking and the interface
# for the evaluation of dictionary elements.

####################
## Bounds checking
####################

# We hook into Julia's bounds checking system. See the Julia documentation.
# This is based on the set of indices, as returned by `indices(dict)`.
# One thing to take into account in our setting is that the map from linear indices
# to native indices and vice-versa requires knowledge of the dictionary. Hence,
# we do some conversions before `indices(dict)` is called and passed on.
checkbounds(dict::Dictionary, I...) = checkbounds(Bool, dict, I...) || throw(BoundsError())

# We make a special case for a linear index
checkbounds(::Type{Bool}, dict::Dictionary, i::Int) = checkindex(Bool, Base.OneTo(length(dict)), i)

# We also convert some native indices to linear indices, before moving on.
# (This is more difficult to do later on, e.g. in checkindex, because that routine
#  does not have access to the dict anymore)
checkbounds(::Type{Bool}, dict::Dictionary, i::NativeIndex) =
    checkbounds(Bool, dict, linear_index(dict, i))
checkbounds(::Type{Bool}, dict::Dictionary, i::AbstractShiftedIndex) =
    checkbounds(Bool, dict, linear_index(dict, i))
checkbounds(::Type{Bool}, dict::Dictionary, i::MultilinearIndex) =
    checkbounds(Bool, dict, linear_index(dict, i))
checkbounds(::Type{Bool}, dict::Dictionary, idx::NTuple{N,Int}) where {N} =
    checkbounds(Bool, d, CartesianIndex(idx))

# And here we call checkbounds_indices with indices(dict)
@inline checkbounds(::Type{Bool}, dict::Dictionary, I...) = checkbounds_indices(Bool, axes(dict), I)

"Return the support of the idx-th basis function. Default is support of the dictionary."
function support(dict::Dictionary, idx)
    checkbounds(dict, idx)
    support(dict)
end
# Warning: the functions above and below may be wrong for certain concrete
# dictionaries, for example for univariate functions with non-connected support.
# Make sure to override, and make sure that the overridden version is called.

"Does the given point lie inside the support of the given function or dictionary?"
in_support(dict::Dictionary, x) = dict_in_support(dict, x)
in_support(dict::Dictionary, idx, x) = dict_in_support(dict, idx, x)

# in_support(dict::Dictionary, x) =
#     dict_in_support(dict, element_type_check(x, domaintype(dict)))
# in_support(dict::Dictionary, idx, x) =
#     dict_in_support(dict, idx, element_type_check(x, domaintype(dict)))
#
# element_type_check(x, T) = _element_type_check(x, T, typeof(x))
# _element_type_check(x, T, S) = _element_type_check(x, T, S, promote_type(S,T))
# _element_type_check(x, T, S, V) = x
# _element_type_check(x, T, S, V) = x
# function _element_type_check(x, T, S, ::Type{Any})
#     @warn "Dictionary with domain type $(T) may not support evaluation with type $(S)."
#     x
# end
#
# # Some special cases
# _element_type_check(x, ::NTuple{N,T}, ::SVector{N,T}) where {T,N} = x
# _element_type_check(x, ::SVector{N,T}, ::NTuple{N,T}) where {T,N} = x

# The mechanism is as follows:
# - in_support(dict::Dictionary, ...) calls dict_in_support
# - any linear index is converted to a native index
# - concrete dictionary should implement dict_in_support
# The reasoning is that concrete dictionaries need not all worry about handling
# linear indices. Yet, they are free to implement other types of indices.
# If a more efficient algorithm is available for linear indices, then the concrete
# dictionary can still intercept the call to in_support.
# The delegation to a method with a different name (dict_in_support) makes it
# substantially easier to deal with ambiguity errors.

# This is the standard conversion to a native_index for an index of type
# Int. This calls a different function, hence it is fine if the native
# index happens to be a linear index (Int).
in_support(dict::Dictionary, idx::Int, x) =
    dict_in_support(dict, native_index(dict, idx), x)

# The default fallback is implemented below in terms of the support of the dictionary:
dict_in_support(dict::Dictionary, idx, x) = default_in_support(dict, idx, x)
dict_in_support(dict::Dictionary, x) = default_in_support(dict, x)

default_in_support(dict::Dictionary, idx, x) = approx_in(x, support(dict, idx), tolerance(dict))
default_in_support(dict::Dictionary, x) = approx_in(x, support(dict), tolerance(dict))

#################################################
## Evaluating dictionary elements and expansions
#################################################

"""
A member function of a dictionary is evaluated using the `eval_element` routine.
It takes as arguments the dictionary, the index of the member function and
the point in which to evaluate.

This function performs bounds checking on the index and also checks whether the
point x lies inside the support of the function. A `BoundsError` is thrown for
an index out of bounds. The value `0` is returned when x is outside the support.

After the check on the index, the function calls `unsafe_eval_element1.` This
function checks whether `x` lies in the support, and then calls
`unsafe_eval_element`. The latter function should be implemented by a concrete
dictionary. Any user who wants to avoid the bounds check or the support check
can intercept `eval_element` or `unsafe_eval_element1` respectively.
"""
function eval_element(dict::Dictionary, idx, x)
    # We convert to a native index before bounds checking
    idxn = native_index(dict, idx)
    @boundscheck checkbounds(dict, idxn)
    unsafe_eval_element1(dict, idxn, x)
end

# For linear indices, bounds checking is very efficient, so we intercept this case
# and only convert to a native index after the bounds check.
function eval_element(dict::Dictionary, idx::Int, x)
    @boundscheck checkbounds(dict, idx)
    unsafe_eval_element1(dict, native_index(dict, idx), x)
end

unsafe_eval_element1(dict::Dictionary, idx, x) =
    in_support(dict, idx, x) ? unsafe_eval_element(dict, idx, x) : zero(codomaintype(dict))

# Catch any index and convert to native index, in case it got through to here
unsafe_eval_element(dict::Dictionary, idx, x) =
    unsafe_eval_element(dict, native_index(dict, idx), x)

"""
Evaluate a member function with a boundscheck on the index, but without checking
the support of the function.
"""
function eval_element_extension(dict::Dictionary, idx, x)
    @boundscheck checkbounds(dict, idx)
    # We skip unsafe_evaluate_element1 and jump to unsafe_eval_element
    unsafe_eval_element(dict, native_index(dict, idx), x)
end

"""
This function is exactly like `eval_element`, but it evaluates the derivative
of the element instead.
"""
eval_element_derivative(dict::Dictionary, idx, x) =
    eval_element_derivative(dict, idx, x, difforder(dict))

function eval_element_derivative(dict::Dictionary, idx, x, order)
    if orderiszero(order)
        eval_element(dict, idx, x)
    else
        idxn = native_index(dict, idx)
        @boundscheck checkbounds(dict, idxn)
        unsafe_eval_element_derivative1(dict, idxn, x, order)
    end
end

function eval_element_derivative(dict::Dictionary, idx::Int, x, order)
    if orderiszero(order)
        eval_element(dict, idx, x)
    else
        @boundscheck checkbounds(dict, idx)
        unsafe_eval_element_derivative1(dict, native_index(dict, idx), x, order)
    end
end

function eval_element_extension_derivative(dict::Dictionary, idx, x, order)
    if orderiszero(order)
        eval_element_extension(dict, idx, x)
    else
        @boundscheck checkbounds(dict, idx)
        unsafe_eval_element_derivative(dict, native_index(dict, idx), x, order)
    end
end

function unsafe_eval_element_derivative1(dict::Dictionary, idx, x, order)
    in_support(dict, idx, x) ? unsafe_eval_element_derivative(dict, idx, x, order) : zero(codomaintype(dict))
end

# Fallback, may throw StackOverFlowError if the method is not defined for native indices
unsafe_eval_element_derivative(dict::Dictionary, idx, x, order) =
    unsafe_eval_element_derivative(dict, native_index(dict, idx), x, order)

eval_gradient(dict::Dictionary, idx, x) = _eval_gradient(dict, idx, x, domaintype(dict))

_eval_gradient(dict, idx, x, ::Type{<:Number}) = eval_element_derivative(dict, idx, x, 1)

function _eval_gradient(dict, idx, x, ::Type{SVector{N,T}}) where {N,T}
    @boundscheck checkbounds(dict, idx)
    if in_support(dict, x)
        SVector{N,codomaintype(dict)}((unsafe_eval_element_derivative(dict, idx, x, dimension_tuple(Val{N}(), dim)) for dim in 1:N)...)
    else
        zero(SVector{N,codomaintype(dict)})
    end
end

"""
Evaluate an expansion given by the set of coefficients in the point x.
"""
function eval_expansion(dict::Dictionary, coefficients, x; options...)
    @assert size(coefficients) == size(dict)
    in_support(dict, x) ? unsafe_eval_expansion(dict, coefficients, x) : zero(span_codomaintype(dict, coefficients))
end

unsafe_eval_expansion(dict::Dictionary, coefficients, x) =
    default_eval_expansion(dict, coefficients, x)

default_eval_expansion(dict::Dictionary, coefficients, x) =
    sum(coefficients[idx]*val for (idx,val) in pointvalues(dict, x))

function eval_expansion(dict::Dictionary, coefficients, grid::AbstractGrid; options...)
    @assert dimension(dict) == GridArrays.dimension(grid)
    @assert size(coefficients) == size(dict)
    # TODO: reenable test once product grids and product sets have compatible types again
    # @assert eltype(grid) == domaintype(dict)

    T = coefficienttype(dict)
    E = evaluation(T, dict, grid; options...)
    E * coefficients
end

# Evaluation of a dictionary means evaluation of all elements.
(dict::Dictionary)(x) = dict_eval(dict, x)

function dict_eval(dict::Dictionary, x)
    result = zeros(dict)
    dict_eval!(result, dict, x)
end

function dict_eval!(result, dict, x)
    for (idx,idxn) in enumerate(ordering(dict))
        result[idx] = unsafe_eval_element1(dict, idxn, x)
    end
    result
end

# evaluation of the dictionary is "unsafe" because the routine can assume that the
# support of x has already been checked
function unsafe_dict_eval(dict::Dictionary, x)
    # We allocate the zeros ourselves, in order to make sure that the output has
    # the same array structure as the dictionary itself. For example, it will be
    # a BlockArray for a composite dictionary, rather than a flat array.
    result = zeros(dict)
    unsafe_dict_eval!(result, dict, x)
end

function unsafe_dict_eval!(result, dict, x)
    for (idx, z) in pointvalues(dict, x)
        result[idx] = z
    end
    result
end

"""
Supertype of iterators over the pointwise evaluation of a dictionary.

The value iterator depends on a point `x` in the support of the dictionary, and
iterates over all the function values of the elements of the dictionary in that
point. Value iterators can be useful if the computation of all elements is required
and computing them in one go is more efficient than computing them one by one
independently.

For example, orthogonal polynomial iterators may implement the three-term recurrence
relation. The command
'[val for val in pointvalues(ops, x)]'
may be more efficient than
'[eval_element(dict, i, x) for i in eachindex(dict)]'
"""
abstract type DictionaryValueIterator{T} end

length(iter::DictionaryValueIterator) = length(dictionary(iter))
size(iter::DictionaryValueIterator) = size(dictionary(iter))
eltype(iter::DictionaryValueIterator{T}) where {T} = T

dictionary(iter::DictionaryValueIterator) = iter.dict
point(iter::DictionaryValueIterator) = iter.x

"Iterate over the values of a dictionary at a point."
struct GenericDictValueIterator{S,D<:Dictionary,I,T} <: DictionaryValueIterator{T}
    dict    ::  D
    x       ::  S
    idxiter ::  I
end

GenericDictValueIterator(Φ::Dictionary, x) = GenericDictValueIterator(Φ, x, eachindex(Φ))
GenericDictValueIterator(Φ::Dictionary, x, I) =
    GenericDictValueIterator{typeof(x),typeof(Φ),typeof(I),codomaintype(Φ)}(Φ, x, I)

function pointvalues(Φ::Dictionary, x)
    # Let's discourage the use of the iterator for points outside the support.
    @assert in_support(Φ, x)
    GenericDictValueIterator(Φ, x)
end

indexiterator(iter::GenericDictValueIterator) = iter.idxiter

function iterate(iter::GenericDictValueIterator)
    state = iterate(indexiterator(iter))
    if state != nothing
        idx, idx_state = state
        (idx, BasisFunctions.unsafe_eval_element(dictionary(iter), idx, point(iter))), idx_state
    else
        nothing
    end
end

function iterate(iter::GenericDictValueIterator, state)
    newstate = iterate(indexiterator(iter), state)
    if newstate != nothing
        idx, idx_state = newstate
        (idx, BasisFunctions.unsafe_eval_element(dictionary(iter), idx, point(iter))), idx_state
    end
end


text:
"""
    interpret_d(d::Number; rules::Union{String, Rules}="cohen1988")

Interpret a standardized difference (*Cohen's d*) with rules of thumb.

Cohen's *d* is a widely-used standardized effect size index. It's equivalent in many ways to a standardized regression coefficient, as both divide the size of the effect by the relevant standard deviations. Thus, instead of being in terms of the original units of X and Y, Cohen's *d* and standardized regression coefficients are expressed in terms of standard deviations.

# Arguments
- `d`: The *d* value.
- `rules`: Can be 'cohen1988', 'sawilowsky2009', or custom set of [Rules](@ref).

# Rules of thumb
- **Cohen (1988)**: The author suggests that *d* = 0.2 should be considered as small, *d* = 0.5 as medium, and *d* = 0.8 as large.
- **Sawilowsky (2009)**: The author suggests that values of 0.1, 0.2, 0.5, 0.8, 1.2 and 2.0 should be considered as very small, small, medium, large, very large and huge, respectively.

# Examples
```jldoctest
interpret_d(0.4, rules="cohen1988")

# output

"small"

# References
- Sawilowsky, S. S. (2009). New effect size rules of thumb.
- Cohen, J. (1988). Statistical power analysis for the behavioral sciences. 2nd.
```
"""
function interpret_d(d::Number; rules::Union{String, Rules}="cohen1988")

    if isa(rules, String)
        if rules == "cohen1988"
            rules = Rules(
                [0.2, 0.5, 0.8],
                ["very small", "small", "medium", "large"],
                true
                )
        elseif rules == "sawilowsky2009"
            rules = Rules(
                [0.1, 0.2, 0.5, 0.8, 1.2, 2.0],
                ["tiny", "very small", "small", "medium", "large", "very large", "huge"],
                true
                )
        else
            throw(ArgumentError("`rules` must be a String ('cohen1988', 'sawilowsky2009') or custom Rules."))
        end
    end

    interpretation = interpret(abs(d), rules)
    return interpretation
end


text:
#==========================================================
This routine initializes an optimization problem to recover
a CNOT gate on a single qudit with 4 energy levels (and 2
guard states). The drift Hamiltonian in the rotating frame
is
              H_0 = - 0.5*ξ_a(a^†a^†aa),
where a is the annihilation operator for the qudit.
Here the control Hamiltonian includes the usual symmetric
and anti-symmetric terms
     H_{sym} = p(t)(a + a^†),    H_{asym} = q(t)(a - a^†)
which come from the rotating frame approximation and hence
we refer to these as "coupled" controls.
The problem parameters are:
                ω_a =  2π × 4.10595   Grad/s,
                ξ_a =  2π × 2(0.1099) Grad/s.
We use Bsplines with carrier waves with frequencies
0, ξ_a, 2ξ_a Grad/s.
==========================================================#
using LinearAlgebra
using Plots
pyplot()
using FFTW
using DelimitedFiles
using Printf
using Ipopt
using Random

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

using Juqbox # quantum control module

Nosc = 1 # Number of coupled sub-systems = oscillators
N = 4 # Number of essential energy levels

Nguard = 2 # Number of extra guard levels
Ntot = N + Nguard # Total number of energy levels

T = 100.0 # Duration of gate.

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
fa = 4.10336
xa = 0.2198
rot_freq = [fa] # Used to calculate the lab frame ctrl function

# setup drift Hamiltonian
number = Diagonal(collect(0:Ntot-1))

H0 = -0.5*(2*pi)*xa* (number*number - number) # xa is in GHz

# lowering matrix
amat = Bidiagonal(zeros(Ntot),sqrt.(collect(1:Ntot-1)),:U) # standard lowering matrix
adag = Array(transpose(amat));
Hsym_ops=[Array(amat + adag)]
Hanti_ops=[Array(amat - adag)]
H0 = Array(H0)

# Estimate time step
maxctrl = 0.001*2*pi * 8.5 #  9, 10.5, 12, 15 MHz

# K1 =  H0 + maxctrl.*( amat +  amat') + 1im* maxctrl.*(amat - amat')
# lamb = eigvals(Array(K1))
# maxeig = maximum(abs.(lamb))

# Pmin = 40
# samplerate1 = maxeig*Pmin/(2*pi)
# nsteps = ceil(Int64, T*samplerate1)

nsteps = calculate_timestep(T, H0, Hsym_ops, Hanti_ops, [maxctrl])
println("# time steps: ", nsteps)

use_bcarrier = true # Use carrier waves in the control pulses?

if use_bcarrier
    Nfreq = 3 # number of carrier frequencies
else
    Nfreq = 1
end

Ncoupled = length(Hsym_ops) # Here, Ncoupled = 1
om = zeros(Ncoupled,Nfreq)

if use_bcarrier
    om[1:Ncoupled,2] .= -2.0*pi *xa       # Note negative sign
    om[1:Ncoupled,3] .= -2.0*pi* 2.0*xa
end
println("Carrier frequencies [GHz]: ", om[1,:]./(2*pi))

maxamp = zeros(Nfreq)

if Nfreq >= 3
    const_fact = 0.45
    maxamp[1] = maxctrl*const_fact
    maxamp[2:Nfreq] .= maxctrl*(1.0-const_fact)/(Nfreq-1) # split the remainder equally
else
    # same threshold for all frequencies
    maxamp .= maxctrl/Nfreq
end

maxpar = maximum(maxamp)

Ident = Matrix{Float64}(I, Ntot, Ntot)

# CNOT trarget
utarget = Ident[1:Ntot,1:N]
utarget[:,3] = Ident[:,4]
utarget[:,4] = Ident[:,3]

omega1 = Juqbox.setup_rotmatrices([N], [Nguard], [fa])

# Compute Ra*utarget
rot1 = Diagonal(exp.(im*omega1*T))

# target in the rotating frame
vtarget = rot1*utarget

# Initial conditions
U0 = Ident[1:Ntot,1:N]

params = Juqbox.objparams([N], [Nguard], T, nsteps, Uinit=U0, Utarget=vtarget, Cfreq=om, Rfreq=rot_freq,
                          Hconst=H0, Hsym_ops=Hsym_ops, Hanti_ops=Hanti_ops)

# initial parameter guess

# D1 smaller than 5 does not work
D1 = 10 # Number of B-spline coefficients per segment
nCoeff = 2*Ncoupled*Nfreq*D1 # factor '2' is for sin/cos

# Random.seed!(12456)

startFromScratch = true # false
startFile="cnot-pcof-opt.dat"

# initial parameter guess
if startFromScratch
    pcof0 = maxpar*0.01 * rand(nCoeff)
    println("*** Starting from random pcof with amplitude ", maxpar*0.01)
else
    # use if you want to read the initial coefficients from file
    pcof0 = vec(readdlm(startFile))
    println("*** Starting from B-spline coefficients in file: ", startFile)
    nCoeff = length(pcof0)
    D1 = div(nCoeff, 2*Ncoupled*Nfreq)  # number of B-spline coeff per control function
end

# min and max coefficient values
useBarrier = true
minCoeff, maxCoeff = Juqbox.assign_thresholds_freq(maxamp, Ncoupled, Nfreq, D1)

samplerate = 32 # only used for plotting
casename = "cnot1" # base file name (used in optimize-once.jl)

maxIter = 200 # 0  # optional argument
lbfgsMax = 250 # optional argument
ipTol = 1e-5   # optional argument
acceptTol = ipTol # 1e-4 # acceptable tolerance
acceptIter = 15

println("*** Settings ***")
println("System Hamiltonian coefficients [GHz]: (fa, xa) =  ", fa, xa)
println("Total number of states, Ntot = ", Ntot, " Total number of guard states, Nguard = ", Nguard)
if use_bcarrier
  println("Using B-spline basis functions with carrier wave, # freq = ", Nfreq)
else
  println("Using regular B-spline basis functions")
end
println("Number of coefficients per spline = ", D1, " Total number of parameters = ", nCoeff)
for q=1:Nfreq
    println("Carrier frequency: ", om[q]/(2*pi), " GHz, max parameter amplitude = ", 1000*maxamp[q]/(2*pi), " MHz")
end
println("Tikhonov coefficients: tik0 = ", params.tik0)

# Estimate number of terms in Neumann series for time stepping (Default 3)
tol = eps(1.0); # machine precision
Juqbox.estimate_Neumann!(tol, params, [maxpar])

wa = Juqbox.Working_Arrays(params,nCoeff)
prob = Juqbox.setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff, maxIter, lbfgsMax, startFromScratch)

println("Initial coefficient vector stored in 'pcof0'")