text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Distribute the computation over multiple patches
#--------------------------------------------------------------------

"""
    Set up an abstract computational product spaces given the grid parameters
    Input: Grid parameters
    Output: Array of computational elements with their individual product spaces.
"""

function setup(params::Parameters{T}) where {T}
    ustops = range(params.ubounds[1], stop=params.ubounds[2], length=params.npatchs[1]+1)
    vstops = range(params.vbounds[1], stop=params.vbounds[2], length=params.npatchs[2]+1)
    AoT = Array{Union{ProductSpace, NTuple{params.nfields, Field}}}(undef, params.npatchs[1], params.npatchs[2])
    for index in CartesianIndices(AoT)
        AoT[index] = ProductSpace(ChebyshevGL{U, params.npoints[1], T}(ustops[index.I[1]], ustops[index.I[1]+1]),
                                  ChebyshevGL{V, params.npoints[2], T}(vstops[index.I[2]], vstops[index.I[2]+1]))
    end
    return AoT
end

"""
    Main code to compute the solution over multiple patches
    Input: grid parameters, excision function, uboundary, vboundary
    Output: Array of ProductSpace or NTuple of solutions
"""
function distribute(params::Parameters{T}, computePatch::Function,
        computeUboundary::Function, computeVboundary::Function, debug::Int=0) where {T}
    AoT = setup(params)
    for index in CartesianIndices(AoT)
        if (debug >= 2)
            println("Computing patch ", index.I)
            ((umin, vmin), (umax, vmax)) = range(AoT[index])
            @printf("  umin = %.2f umax = %.2f\n", umin, umax)
            @printf("  vmin = %.2f vmax = %.2f\n", vmin, vmax)
        end
        uboundary = index.I[1] == 1 ?  computeUboundary(AoT[index]) : extractUboundary(AoT[index - CartesianIndex((1,0))], :outgoing)
        vboundary = index.I[2] == 1 ?  computeVboundary(AoT[index]) : extractVboundary(AoT[index - CartesianIndex((0,1))], :outgoing)
        AoT[index] = computePatch(AoT[index], uboundary, vboundary)
    end
    return AoT
end

function setup_(params::Parameters{T})::Array{Union{ProductSpace, Future, NTuple{3, Field}}} where {T}
    ustops = range(params.ubounds[1], stop=params.ubounds[2], length=params.npatchs[1]+1)
    vstops = range(params.vbounds[1], stop=params.vbounds[2], length=params.npatchs[2]+1)
    AoT = Array{Union{ProductSpace, Future, NTuple{3, Field}}}(undef, params.npatchs[1], params.npatchs[2])
    for index in CartesianIndices(AoT)
        AoT[index] = @spawnat :any ProductSpace(ChebyshevGL{U, params.npoints[1], T}(ustops[index.I[1]], ustops[index.I[1]+1]),
                                         ChebyshevGL{V, params.npoints[2], T}(vstops[index.I[2]], vstops[index.I[2]+1]))
    end
    return AoT
end

function distribute_(params::Parameters{T}, computePatch::Function,
                     computeUboundary::Function, computeVboundary::Function, debug::Int=0) where {T}
    function computeUboundary_(PS::Future)::Future where {S1, S2}
        @spawnat :any computeUboundary(fetch(PS))
    end

    function computeVboundary_(PS::Future)::Future where {S1, S2}
        @spawnat :any computeVboundary(fetch(PS))
    end

    function extractUboundary_(U::Future, boundarytype::Symbol)::Future
        @spawnat :any extractUboundary(fetch(U), boundarytype)
    end

    function extractVboundary_(U::Future, boundarytype::Symbol)::Future
        @spawnat :any extractVboundary(fetch(U), boundarytype)
    end

    function computePatch_(PS::Future, Uboundary::Future, Vboundary::Future)::Future where {S1, S2}
        @spawnat :any computePatch(fetch(PS), fetch(Uboundary), fetch(Vboundary))
    end

    AoT = setup_(params)
    for index in CartesianIndices(AoT)
        if (debug >= 2)
            println("Computing patch ", index.I)
            println("  (umin, vmin) = ", range(fetch(AoT[index]))[1])
            println("  (umax, vmax) = ", range(fetch(AoT[index]))[2])
        end
        uboundary = index.I[1] == 1 ?  computeUboundary_(AoT[index]) : extractUboundary_(AoT[index - CartesianIndex((1,0))], :outgoing)
        vboundary = index.I[2] == 1 ?  computeVboundary_(AoT[index]) : extractVboundary_(AoT[index - CartesianIndex((0,1))], :outgoing)
        AoT[index] = computePatch_(AoT[index], uboundary, vboundary)
    end
    return fetch.(AoT)
end


text:
include("macros.jl")
include("dynamic_dsl.jl")
include("static_dsl.jl")

ext: jl

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
using DTW
using Base.Test

# write your own tests here
@test 1 == 1

x=[1:5]
y=[1:2:5]

a = dtw(x,y)


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

n = 10
plot(x=fill(1, n), y=rand(n), Geom.bar)


text:
using Logging

"""
Unhandled event exception.

This gets thrown if an event is not handled by the state machine when
`handle_event!()` is called. This usually means the root state machine state
does not have an event handler for the event type, or is not returning `true`.
"""
struct HsmUnhandledEventError <: Exception
    msg::String
end

"""
State transition failed exception.

This gets thrown if the destination state is not a part of the state machine
when `transition_to_state!()` is called.
"""
struct HsmStateTransitionError <: Exception
    msg::String
end

"""
    AbstractHsmEvent

Abstract HSM event type.

Concrete events must extend `AbstractHsmEvent`.
"""
abstract type AbstractHsmEvent end

"""
    AbstractHsmState

Abstract HSM state type.

Concrete states must extend `AbstractHsmState`, and either contain an
`HsmStateInfo`  struct called `state_info`, or implement the `parent_state()`,
`parent_state!()`, `active_substate()`, and `active_substate!()` getters and
setters.

Pass the parent state to the HsmStateInfo constructor in the concrete state's
constructor if using `HsmStateInfo`.

Example:

```julia
struct MyState <: AbstractHsmState
    state_info::HsmStateInfo
    app_info::MyStatefulApplicationInfo

    MyState(parent_state, app_info) = new(HSM.HsmStateInfo(parent_state), app_info)
end
```

Otherwise, call `parent_state(obj, <parent state>)` and
`active_substate(obj, nothing)` when constructing your concrete state.
"""
abstract type AbstractHsmState end

"""
    parent_state(obj::AbstractHsmState)

parent_state getter.

Extend this for your concrete class if your concrete struct does not include
`state_info::HsmStateInfo`.
"""
function parent_state(obj::AbstractHsmState)
    obj.state_info.parent_state
end

"""
    parent_state!(obj::AbstractHsmState, value::Union{AbstractHsmState, Nothing})

parent_state setter.

Extend this for your concrete class if your concrete struct does not include
`state_info::HsmStateInfo`.
"""
function parent_state!(obj::AbstractHsmState, value::Union{AbstractHsmState, Nothing})
    obj.state_info.parent_state = value
end

"""
    active_substate(obj::AbstracHsmState)

active_substate getter.

Extend this for your concrete class if your concrete struct does not include
`state_info::HsmStateInfo`.
"""
function active_substate(obj::AbstractHsmState)
    obj.state_info.active_substate
end

"""
    active_substate!(obj::AbstractHsmState, value::Union{AbstractHsmState, Nothing})

active_substate setter.

Extend this for your concrete class if your concrete struct does not include
`state_info::HsmStateInfo`.
"""
function active_substate!(obj::AbstractHsmState, value::Union{AbstractHsmState, Nothing})
    obj.state_info.active_substate = value
end

"""
    HsmStateInfo

State information needed by HSM interfaces that act on `AbstractHsmState`s.

* `parent_state::AbstractHsmState` -- initialized to the parent state of the
    concrete state, or `nothing` if it is the root state machine state.
* `active_substate::AbstractHsmState` -- initialized to `nothing`.
"""
mutable struct HsmStateInfo
    parent_state::Union{AbstractHsmState, Nothing}
    active_substate::Union{AbstractHsmState, Nothing}

    """
        HsmStateInfo(parent_state::Union{AbstractHsmState, Nothing})

    Constructor for HsmStateInfo.

    `parent_state` argument is the state containing the current state, or
    `nothing`, if the state is the root state machine.
    """
    function HsmStateInfo(parent_state)
        new(parent_state, nothing)
    end
end

#####
# Global Handlers
#####

"""
    on_event!(state::AbstractHsmState, event::AbstractHsmEvent)

Global event handler.

Create event handlers for your derived states and events to handle those
specific events in those states.

```julia
struct MyState <: HSM.AbstractHsmState ; end
struct MyEvent <: HSM.AbstractHsmEvent ; end

function on_event!(state::MyState, event::MyEvent)
    # process event.
    return true
end
```

State-specific event handlers are expected to return `true` when they handle an
event, unless they specifically want their parent state to process the event as
well.

State-specific instances of `on_event!()` correspond to "event / action" labels
in state blocks or on state transitions in UML state charts.

```plantuml
@startuml
state ExampleMachine {
    state MyState
    state MyOtherState

    MyState --> MyOtherState : MyEvent / do something
    MyOtherState : MyEvent / do something
}
@enduml
```

The global event handler is called as a last resort, and is expected to return
`false` to generate an `HsmUnhandledEventError`.
"""
function on_event!(state::AbstractHsmState, event::AbstractHsmEvent)
    @debug "on_event!(AbstractHsmState, AbstractHsmEvent) for ($(string(typeof(state))), $(string(typeof(event))))"
    return false
end

"""
    on_entry!(state::AbstractHsmState)

Global state entry handler.

Create entry handlers for your derived states when you want that state to
execute something whenever it is entered. This includes transitions to itself.

```julia
struct MyState <: HSM.AbstractHsmState ; end

function on_entry!(state::MyState)
    # do something.
end
```

State-specific instances of `on_entry!()` correspond to "on entry / action"
labels in state blocks in UML state charts.

```plantuml
@startuml
state Machine {
    state MyState
    MyState : on entry / do something
}
@enduml
```

The global state entry handler is the default behavior of a state. It does not
do anything.
"""
function on_entry!(state::AbstractHsmState)
    @debug "on_entry!(AbstractHsmState) for $(string(typeof(state)))"
    # do nothing.
end

"""
    on_entry!(state::AbstractHsmState)

Global state exit handler.

Create exit handlers for your derived states when you want that state to execute
something whenever it is exited. This includes transitions to itself.

```julia
struct MyState <: HSM.AbstractHsmState ; end

function on_exit!(state::MyState)
    # do something.
end
```

State-specific instances of `on_exit!()` correspond to "on exit / action"
labels in state blocks in UML state charts.

```plantuml
@startuml
state Machine {
    state MyState
    MyState : on exit / do something
}
@enduml
```

The global state exit handler is the default behavior of a state. It does not do
anything.
"""
function on_exit!(state::AbstractHsmState)
    @debug "on_exit!(AbstractHsmState) for $(string(typeof(state)))"
    # do nothing.
end

"""
    on_initialize!(state::AbstractHsmState)

Global state initializer.

Create initializers for states that must be initialized upon entry. Generally,
this is used when a state has sub-states, where one must be transitioned to
when the state is entered. For example, the root state machine.

```julia
struct Machine <: HSM.AbstractHsmState ; end
struct MyState <: HSM.AbstractHsmState ; end

machine = Machine(nothing)
my_state = MyState(machine)

function on_initialize!(state::Machine)
    transition_to_state!(my_state)
end
```

State-specific instances of `on_initialize!()` correspond to the starting point
trasitions in UML state charts.

```plantuml
@startuml
state Machine {
    state MyState

    [*] --> MyState
}
@enduml
```

The global state initializer is the default behavior of a state. It does not do
anything.
"""
function on_initialize!(state::AbstractHsmState)
    @debug "on_initialize!(AbstractHsmState) for $(string(typeof(state)))"
    # Do nothing.
end

#####
# Internals
#####

"""
    root_state(current_state::AbstractHsmState)

Returns the root state of `current_state`.

Follows the state's parent state up, until there is no more parent state.

Example:

Given the following state machine:

```plantuml
@startuml
state A {
    state B
    state C {
        state D
        state E
    }
}
@enduml
```

The root state of any state in the state machine is `A`.
"""
function root_state(current_state::AbstractHsmState)
    s = current_state
    while !isnothing(parent_state(s))
        s = parent_state(s)
    end
    return s
end

"""
    active_state(current_state::AbstractHsmState)

Returns the active state of the state machine that `current_state` resides in.

Example:

Given the following state machine:

```plantuml
@startuml
state A {
    state B
    state C {
        state D
        state E
    }
}
@enduml

If state `B` is the last state transitioned to, then the active state for any
state in the state machine will be `B`.
```
"""
function active_state(current_state::AbstractHsmState)
    s = root_state(current_state)
    while !isnothing(active_substate(s))
        s = active_substate(s)
    end
    return s
end

"""
    common_parent(left_state::AbstractHsmState, right_state::AbstractHsmState)

Returns the common parent of `left_state` and `right_state`.

Example:

Given the following state machine:

```plantuml
@startuml
state A {
    state B
    state C {
        state D
        state E
    }
}
@enduml
```

The common parent of `A` and any other state is `A`.

The common parent of `B` and any other state is `A`.

The common parent of `C` and any of its substates is `C`.

The common parent of `D` and `E` is `C`.
"""
function common_parent(
    left_state::AbstractHsmState,
    right_state::AbstractHsmState
)
    l = left_state
    while !isnothing(l)
        r = right_state
        while !isnothing(r)
            if r == l
                # Common parent found.
                return r
            end

            r = parent_state(r)
        end

        l = parent_state(l)
    end

    # No common parent.
    return nothing
end

#####
# Public Interfaces
#####

"""
    handle_event!(state_machine::AbstractHsmState, event::AbstractHsmEvent)

Pass `event` to `state_machine` for processing.

This is the public interface to use to pass an event to a state machine for
processing.
"""
function handle_event!(state_machine::AbstractHsmState, event::AbstractHsmEvent)
    @debug "handle_event!($(string(typeof(state_machine))), $(string(typeof(event))))"

    handled = false

    s = active_state(state_machine)
    while !isnothing(s) && !(handled = on_event!(s, event))
        # Event not handled by current state. Try the parent.
        s = parent_state(s)
    end

    if !handled
        # Unhandled event. This is a really good indicator that, at the very
        # least, your root state machine does not have handlers for all possible
        # events.
        throw(
            HsmUnhandledEventError("Unhandled Event: $(string(typeof(event))). Does your root state machine have an event handler (`on_event!()`) for this event?")
        )
    end
end

"""
    transition_to_state!(machine::AbstractHsmState, state::AbstractHsmState)

Change the active state of `state_machine` to `state`.

This is the interface to use in a state's event handler (`on_event!()`,
`on_initialize()`) to transition states.

```julia
struct Machine <: HSM.AbstractHsmState ; end
struct MyState <: HSM.AbstractHsmState ; end
struct MyOtherState <: HSM.AbstractHsmState ; end
struct MyEvent <: HSM.AbstractHsmEvent ; end

machine = Machine(nothing)
my_state = MyState(machine)
my_other_state = MyOtherState(machine)

function on_event!(state::MyState, event::MyEvent)
    # Active state becomes my_other_state.
    transition_to_state!(machine, my_other_state)
end
```

This interface corresponds to arrows from one state to another in a UML state
chart.

```plantuml
@startuml
state Machine {
    state MyState
    state MyOtherState

    MyState --> MyOtherState
}
@enduml
```
"""
function transition_to_state!(
    state_machine::AbstractHsmState,
    state::AbstractHsmState
)
    @debug "transition_to_state!($(string(typeof(state_machine))), $(string(typeof(state))))"

    s = active_state(state_machine)
    cp = common_parent(s, state)

    if isnothing(cp)
        throw(
            HsmStateTransitionError(
                "Destination state " *
                string(typeof(state)) *
                " does not exist in state machine " *
                string(typeof(state_machine))
            )
        )
    end

    # Call `on_exit!()` from active state up to, but not including common
    # parent.
    while s != cp
        on_exit!(s)
        s = parent_state(s)
    end

    # Update active state pointers from common parent to `state`.
    active_substate!(state, nothing)
    s = state
    while s != cp
        active_substate!(parent_state(s), s)
        s = parent_state(s)
    end

    # Call `on_entry!()` for common parent's active state to `state`.
    s = active_substate(cp)
    while !isnothing(s)
        on_entry!(s)
        s = active_substate(s)
    end

    on_initialize!(state)
end

"""
    transition_to_shallow_history!(state_machine::AbstractHsmState, state::AbstractHsmState)

Change the active state of `state_machine`, following the active state of
`state` as much as one layer.

This is the interface to use in a state's event handler (`on_event!()`) to
transition to the history of a state.

```julia
struct Machine <: HSM.AbstractHsmState ; end
struct MyState <: HSM.AbstractHsmState ; end
struct MyOtherState <: HSM.AbstractHsmState ; end
struct MySubState1 <: HSM.AbstractHsmState ; end
struct MySubState2 <: HSM.AbstractHsmEvent ; end
struct MyEvent <: HSM.AbstractHsmEvent ; end

machine = Machine(nothing)
my_state = MyState(machine)
my_other_state = MyOtherState(machine)
my_sub_state_1 = MySubState1(my_other_state)
my_sub_state_2 = MySubState2(my_other_state)

function on_event!(state::MyState, event::MyEvent)
    # Will transition to last active state in my_other_state -- my_other_state,
    # my_sub_state_1, or my_sub_state_2.
    transition_to_shallow_history!(machine, my_other_state)
end
```

This interface corresponds to an arrow to a history marker in a state.

```plantuml
@startuml
state Machine {
    state MyState
    state MyOtherState {
        state MySubState1
        state MySubState2
    }

    MyState --> MyOtherState[H]
}
@enduml
```
"""
function transition_to_shallow_history!(
    state_machine::AbstractHsmState,
    state::AbstractHsmState
)
    @debug "transition_to_shallow_history!($(string(typeof(state_machine))), $(string(typeof(state))))"

    as = active_substate(state)
    s = isnothing(as) ? state : as
    transition_to_state!(state_machine, s)
end

"""
    transition_to_deep_history!(state_machine::AbstractHsmState, state::AbstractHsmState)

Change the active state of `state_machine`, following the active state of
`state` as far as it goes.

This is the interface to use in a state's event handler (`on_event!()`) to
transition to the deep history of a state.

```julia
struct Machine <: HSM.AbstractHsmState ; end
struct MyState <: HSM.AbstractHsmState ; end
struct MyOtherState <: HSM.AbstractHsmState ; end
struct MySubState1 <: HSM.AbstractHsmState ; end
struct MySubState2 <: HSM.AbstractHsmEvent ; end
struct MySubSubState1 <: HSM.AbstractHsmEvent ; end
struct MySubSubState2 <: HSM.AbstractHsmEvent ; end
struct MyEvent <: HSM.AbstractHsmEvent ; end

machine = Machine(nothing)
my_state = MyState(machine)
my_other_state = MyOtherState(machine)
my_sub_state_1 = MySubState1(my_other_state)
my_sub_state_2 = MySubState2(my_other_state)
my_sub_sub-state_1 = MySubSubState1(my_sub_state_1)
my_sub_sub-state_2 = MySubSubState2(my_sub_state_1)

function on_event!(state::MyState, event::MyEvent)
    # Will transition last active state in my_other_state -- my_other_state,
    # my_sub_state_1, my_sub-state_2, my_sub_sub_state_1, or my_sub_sub_state_2.
    transition_to_shallow_history!(machine, my_other_state)
end
```

This interface corresponds to an arrow to a deep history marker in a state.

```plantuml
@startuml
state Machine {
    state MyState
    state MyOtherState {
        state MySubState1 {
            state MySubSubState1
            state MySubSubState2
        }
        state MySubState2
    }

    MyState --> MyOtherState[H*]
}
@enduml
```
"""
function transition_to_deep_history!(
    state_machine::AbstractHsmState,
    state::AbstractHsmState
)
    @debug "transition_to_deep_history!($(string(typeof(state_machine))), $(string(typeof(state))))"

    s = state
    while !isnothing(active_substate(s))
        s = active_substate(s)
    end
    transition_to_state!(state_machine, s)
end


text:
#= As in Binary units, we prefer to encode binary states with
floats 0.0 / 1.0 instead of Booleans, because they hit faster
BLAS linear algebra routines. =#

"""
    Potts(θ)

Potts layer, with external fields `θ`.
Encodes categorical variables as one-hot vectors.
The number of classes is the size of the first dimension.
"""
struct Potts{N, T, A <: AbstractArray{T,N}} <: AbstractLayer{N}
    θ::A
end
Potts(::Type{T}, q::Int, n::Int...) where {T} = Potts(zeros(T, q, n...))
Potts(q::Int, n::Int...) = Potts(Float64, q, n...)

Base.repeat(l::Potts, n::Int...) = Potts(repeat(l.θ, n...))

free_energies(layer::Potts, inputs::Union{Real,AbstractArray} = 0) = -logsumexp(layer.θ .+ inputs; dims=1)
mean_from_inputs(layer::Potts, inputs::Union{Real,AbstractArray} = 0) = softmax(layer.θ .+ inputs; dims=1)
mean_abs_from_inputs(layer::Potts, inputs::Union{Real,AbstractArray} = 0) = mean_from_inputs(layer, inputs)
std_from_inputs(layer::Potts, inputs::Union{Real,AbstractArray} = 0) = sqrt.(var_from_inputs(layer, inputs))

function mode_from_inputs(layer::Potts, inputs::Union{Real,AbstractArray} = 0)
    return layer.θ .+ inputs .== maximum(layer.θ .+ inputs; dims=1)
end

function var_from_inputs(layer::Potts, inputs::Union{Real,AbstractArray} = 0)
    μ = mean_from_inputs(layer, inputs)
    return μ .* (1 .- μ)
end

function meanvar_from_inputs(layer::Potts, inputs::Union{Real,AbstractArray} = 0)
    μ = mean_from_inputs(layer, inputs)
    ν = μ .* (1 .- μ)
    return μ, ν
end

function sample_from_inputs(layer::Potts, inputs::Union{Real,AbstractArray} = 0)
    c = categorical_sample_from_logits(layer.θ .+ inputs)
    return onehot_encode(c, 1:size(layer, 1))
end


text:
function Base.display(plot::Plot)
    @assert length(plot.layers) > 0

    # compute ranges
    xmin, xmax = Inf, -Inf
    ymin, ymax = Inf, -Inf
    for layer in plot.layers
        xmin = min(xmin, x_minimum(layer))
        xmax = max(xmax, x_maximum(layer))
        ymin = min(ymin, y_minimum(layer))
        ymax = max(ymax, y_maximum(layer))
    end

    # adjust ranges
    xmin, xmax = adjustrange(xmin, xmax)
    ymin, ymax = adjustrange(ymin, ymax)

    # set window
    viewport = (0.13, 0.95, 0.10, 0.95)
    setviewport(viewport...)
    setcharheight(0.020)
    settextcolorind(82)
    clearws()
    setwindow(xmin, xmax, ymin, ymax)

    # grid and axes
    xmajor = ymajor = 5
    xtick = tick(xmin, xmax) / xmajor
    ytick = tick(ymin, ymax) / ymajor
    grid(xtick, ytick, xmax, ymax, xmajor, ymajor)
    ticksize = 0.0125 * (viewport[2] - viewport[1])
    axes(xtick, ytick, xmin, ymin,  xmajor,  ymajor,  ticksize)
    axes(xtick, ytick, xmax, ymax, -xmajor, -ymajor, -ticksize)
    text(0.50, 0.03, string(plot.aes.x))
    setcharup(-1, 0)
    text(0.03, 0.50, string(plot.aes.y))

    # plot layers
    for layer in plot.layers
        plotlayer(layer)
    end
    updatews()

    return nothing
end


text:
# Script to test text layout observable updating on toggle.
using GLMakie
let
    f = Figure(backgroundcolor = :gray70);
    passbox = f[1,2] = Axis(f[1,2], title = "Passes")
    toggle = Toggle(f)
    label = Label(f, lift(x -> x ? "ON TOGGLE" : "OFF TOGGLE", toggle.active), textsize = 48f0)
    f[1,2] = grid!(hcat([toggle], label), tellheight = false)
    textbox = f[1,1] = Axis(f[1,1], title = "Compiled code")
    text!(textbox, lift(x -> x ? "ON TOGGLE" : "OFF TOGGLE", toggle.active))
    f
end


text:
using EffectiveWaves, Test

@testset "match purely numerical solution" begin

## high attenuating material
    medium = Medium(ρ=1.0, c=1.0)
    # cs = [0.1,0.5]
    cs = [0.5]
    species = [
        Specie(ρ = c, c = 0.6-c, r = 2*c, volfrac=c)
    for c in cs]

    ω = 1.1
    k = ω/medium.c
    θin = 0.3
    tol = 1e-8
    hankel_order = 2
    radius_multiplier = 1.005

    wave_effs_arr = [
        effective_waves(ω, medium, [s];
            radius_multiplier=radius_multiplier,
            hankel_order=hankel_order,
            mesh_points=6,
            num_wavenumbers=28,
            tol = tol,  θin = θin,
            extinction_rescale = false)
    for s in species];

   # use only 20 least attenuating
   wave_effs_arr = [w[1:20] for w in wave_effs_arr]

    match_ws = [
        MatchWave(ω, medium, species[i];
            radius_multiplier=radius_multiplier,
            hankel_order=hankel_order,
            θin = θin, tol = tol,
            wave_effs = wave_effs_arr[i],
            max_size=80)
            # x = 0.0:(radius_multiplier*2*species[i].r/100):4)
    for i in eachindex(species)];

    @test maximum(match_error.(match_ws)) < 1e-4

    avgs = [
        AverageWave(ω, medium, species[i];
                radius_multiplier=radius_multiplier,
                hankel_order=hankel_order,
                tol = tol, θin = θin,
                wave_effs = wave_effs_arr[i], max_size=700)
    for i in eachindex(species)]

    R_ms = [reflection_coefficient(ω, match_ws[i], medium, species[i]) for i in eachindex(species)]
    R_ds = [reflection_coefficient(ω, avgs[i], medium, species[i]) for i in eachindex(species)]

    avg_eff = AverageWave(match_ws[1].x_match[end]:0.002:40, match_ws[1].effective_waves);
    R1 = reflection_coefficient(ω, avg_eff, medium, species[1])
    R2 = reflection_coefficient(ω, match_ws[1].effective_waves, medium, species; x=avg_eff.x[1])
    @test norm(R1 - R2) < 1e-7

    @test norm(R_ms - R_ds) < 5e-4

    map(eachindex(species)) do i
        j0 = findmin(abs.(avgs[i].x .- match_ws[i].x_match[1]))[2]
        x0 = avgs[i].x[j0+1:end]
        avg_m = AverageWave(x0, match_ws[i].effective_waves)
        maximum(abs.(avgs[i].amplitudes[j0+1:end,:,:][:] - avg_m.amplitudes[:]))
        @test norm(avgs[i].amplitudes[j0+1:end,:,:][:] - avg_m.amplitudes[:])/norm(avg_m.amplitudes[:]) < 4e-3
        @test maximum(abs.(avgs[i].amplitudes[j0+1:end,:,:][:] - avg_m.amplitudes[:])) < 2e-3
    end
end


text:
mutable struct fPEPS
    Nx::Int
    Ny::Int
    A_::AbstractMatrix{ITensor}

    fPEPS() = new(0, 0, Matrix{ITensor}(),0,0)

    fPEPS(Nx::Int, Ny::Int, A::Matrix{ITensor}) = new(Nx, Ny, A)
    function fPEPS(sites, lattice::Lattice, Nx::Int, Ny::Int; mindim::Int=1, is_gpu::Bool=false)
        p           = Matrix{ITensor}(undef, Ny, Nx)
        right_links = [ Index(mindim, "Link,c$j,r$i,r") for i in 1:Ny, j in 1:Nx ]
        up_links    = [ Index(mindim, "Link,c$j,r$i,u") for i in 1:Ny, j in 1:Nx ]
        T           = is_gpu ? cuITensor : ITensor
        @inbounds for ii in eachindex(sites)
            col = div(ii-1, Ny) + 1
            row = mod(ii-1, Ny) + 1
            old_s = sites[ii]
            s = Index(dim(old_s), "Site,r$row,c$col")
            if 1 < row < Ny && 1 < col < Nx
                p[row, col] = T(right_links[row, col], up_links[row, col], right_links[row, col-1], up_links[row-1, col], s)
            elseif row == 1 && 1 < col < Nx
                p[row, col] = T(right_links[row, col], up_links[row, col], right_links[row, col-1], s)
            elseif 1 < row < Ny && col == 1
                p[row, col] = T(right_links[row, col], up_links[row, col], up_links[row-1, col], s)
            elseif row == Ny && 1 < col < Nx
                p[row, col] = T(right_links[row, col], right_links[row, col-1], up_links[row-1, col], s)
            elseif 1 < row < Ny && col == Nx
                p[row, col] = T(up_links[row, col], right_links[row, col-1], up_links[row-1, col], s)
            elseif row == Ny && col == 1
                p[row, col] = T(right_links[row, col], up_links[row-1, col], s)
            elseif row == Ny && col == Nx
                p[row, col] = T(right_links[row, col-1], up_links[row-1, col], s)
            elseif row == 1 && col == 1
                p[row, col] = T(right_links[row, col], up_links[row, col], s)
            elseif row == 1 && col == Nx
                p[row, col] = T(up_links[row, col], right_links[row, col-1], s)
            end
            @assert p[row, col] == p[ii]
        end
        new(Nx, Ny, p)
    end
end
Base.eltype(A::fPEPS) = eltype(A.A_[1,1])

function cudelt(left::Index, right::Index)
    d_data   = CUDA.zeros(Float64, ITensors.dim(left), ITensors.dim(right))
    ddi = diagind(d_data, 0)
    d_data[ddi] = 1.0
    delt = cuITensor(vec(d_data), left, right)
    return delt
end

function mydelt(left::Index, right::Index)
    d_data   = zeros(Float64, ITensors.dim(left), ITensors.dim(right))
    ddi = diagind(d_data, 0)
    d_data[ddi] = ones(Float64, length(ddi))
    delt = ITensor(vec(d_data), left, right)
    return delt
end

function checkerboardfPEPS(sites, Nx::Int, Ny::Int; mindim::Int=1)
    lattice = square_lattice(Nx, Ny,yperiodic=false)
    A = fPEPS(sites, lattice, Nx, Ny, mindim=mindim)
    @inbounds for ii ∈ eachindex(sites)
        row = div(ii-1, Nx) + 1
        col = mod(ii-1, Nx) + 1
        spin_side = isodd(row - 1) ⊻ isodd(col - 1) ? 2 : 1
        si     = firstind(A[ii], "Site")
        lis = filter(inds(A[ii]), "Link")
        ivs = [li(1) for li in lis]
        ivs = vcat(ivs, si(spin_side))
        A[ii][ivs...] = 1.0
    end
    for row in 1:Ny, col in 1:Nx
        A[row, col] += randomITensor(inds(A[row, col]))/10.0
    end
    return A
end

function randomfPEPS(sites, Nx::Int, Ny::Int; mindim::Int=1)
    lattice = square_lattice(Nx,Ny,yperiodic=false)
    A = fPEPS(sites, lattice, Nx, Ny, mindim=mindim)
    @inbounds for ii ∈ eachindex(sites)
        randn!(A[ii])
        normalize!(A[ii])
    end
    return A
end

is_gpu(A::fPEPS)    = all(is_gpu.(A[:,:]))
is_gpu(A::Vector{ITensor})    = all(is_gpu.(A[:]))
is_gpu(A::ITensor) = (NDTensors.data(store(A)) isa CuArray)

include("environments.jl")
include("ancillaries.jl")
include("gauge.jl")
include("observables.jl")
include("hamiltonians.jl")

function randomCufPEPS(sites, Nx::Int, Ny::Int; mindim::Int=1)
    lattice = square_lattice(Nx,Ny,yperiodic=false)
    A = fPEPS(sites, lattice, Nx, Ny; mindim=mindim, is_gpu=true)
    @inbounds for ii ∈ eachindex(sites)
        randn!(A[ii])
        normalize!(A[ii])
    end
    return A
end

function cufPEPS(A::fPEPS)
    Ny, Nx = size(A)
    cA     = similar(A)
    @inbounds for i in 1:Ny, j in 1:Nx
        cA[i, j] = cuITensor(A[i, j])
    end
    return cA
end

function Base.collect(cA::fPEPS)
    Ny, Nx = size(A)
    A = similar(cA)
    @inbounds for i in 1:Ny, j in 1:Nx
        A[i, j] = collect(cA[i, j])
    end
    return A
end

store(A::fPEPS)   = A.A_
Base.size(A::fPEPS) = (A.Ny, A.Nx)

Base.getindex(A::fPEPS, i::Integer, j::Integer) = getindex(store(A), i, j)::ITensor
Base.getindex(A::fPEPS, ::Colon,    j::Integer) = getindex(store(A), :, j)::Vector{ITensor}
Base.getindex(A::fPEPS, i::Integer, ::Colon)    = getindex(store(A), i, :)::Vector{ITensor}
Base.getindex(A::fPEPS, ::Colon,    j::UnitRange{Int}) = getindex(store(A), :, j)::Matrix{ITensor}
Base.getindex(A::fPEPS, i::UnitRange{Int}, ::Colon)    = getindex(store(A), i, :)::Matrix{ITensor}
Base.getindex(A::fPEPS, ::Colon, ::Colon)       = getindex(store(A), :, :)::Matrix{ITensor}
Base.getindex(A::fPEPS, i::Integer)             = getindex(store(A), i)::ITensor

Base.setindex!(A::fPEPS, val::ITensor, i::Integer, j::Integer)       = setindex!(store(A), val, i, j)
Base.setindex!(A::fPEPS, vals::Vector{ITensor}, ::Colon, j::Integer) = setindex!(store(A), vals, :, j)
Base.setindex!(A::fPEPS, vals::Vector{ITensor}, i::Integer, ::Colon) = setindex!(store(A), vals, i, :)
Base.setindex!(A::fPEPS, vals::Matrix{ITensor}, ::Colon, j::UnitRange{Int}) = setindex!(store(A), vals, :, j)
Base.setindex!(A::fPEPS, vals::Matrix{ITensor}, i::UnitRange{Int}, ::Colon) = setindex!(store(A), vals, i, :)

Base.copy(A::fPEPS)    = fPEPS(A.Nx, A.Ny, copy(store(A)))
Base.similar(A::fPEPS) = fPEPS(A.Nx, A.Ny, similar(store(A)))

function Base.show(io::IO, A::fPEPS)
  print(io,"fPEPS")
  (size(A)[1] > 0 && size(A)[2] > 0) && print(io,"\n")
  @inbounds for i in 1:A.Nx, j in 1:A.Ny
      println(io,"$i $j $(A[i,j])")
  end
end

@enum Op_Type Field=0 Vertical=1 Horizontal=2 Diag=3
struct Operator
    sites::Vector{Pair{Int,Int}}
    ops::Vector{ITensor}
    site_ind::Index
    dir::Op_Type
end

Base.show(io::IO, o::Operator) = show(io, "Direction: $(o.dir)\nSites [row => col]: $(o.sites)\n")

getDirectional(ops::Vector{Operator}, dir::Op_Type) = collect(filter(x->x.dir==dir, ops))

include("single_site.jl")
include("two_site_hori.jl")
include("two_site_vert.jl")

function spinI(s::Index; is_gpu::Bool=false)::ITensor
    I_data      = is_gpu ? CUDA.zeros(Float64, ITensors.dim(s)*ITensors.dim(s)) : zeros(Float64, ITensors.dim(s), ITensors.dim(s))
    idi         = diagind(reshape(I_data, ITensors.dim(s), ITensors.dim(s)), 0)
    I_data[idi] = is_gpu ? CUDA.ones(Float64, ITensors.dim(s)) : ones(Float64, ITensors.dim(s))
    I           = is_gpu ? cuITensor( I_data, IndexSet(s, s') ) : itensor(I_data, IndexSet(s, s'))
    return I
end

function intraColumnGauge(A::fPEPS, col::Int; kwargs...)::fPEPS
    Ny, Nx = size(A)
    @inbounds for row in reverse(2:Ny)
        @debug "\tBeginning intraColumnGauge for col $col row $row"
        cmb_is     = IndexSet(firstind(A[row, col], "Site"))
        if col > 1
            cmb_is = IndexSet(cmb_is..., commonind(A[row, col], A[row, col - 1]))
        end
        if col < Nx
            cmb_is = IndexSet(cmb_is..., commonind(A[row, col], A[row, col + 1]))
        end
        cmb = combiner(cmb_is, tags="CMB")
        Lis = IndexSet(combinedind(cmb)) #cmb_is
        if row < Ny
            Lis = IndexSet(Lis..., commonind(A[row, col], A[row + 1, col]))
        end
        old_ci         = commonind(A[row, col], A[row-1, col])
        Ac             = A[row, col]*cmb
        U, S, V        = svd(Ac, Lis; kwargs...)
        A[row, col]    = U*cmb
        A[row-1, col] *= (S*V)
        new_ci         = commonind(A[row, col], A[row-1, col])
        A[row, col]    = replaceind!(A[row, col], new_ci, old_ci)
        A[row-1, col]  = replaceind!(A[row-1, col], new_ci, old_ci)
    end
    return A
end

function simpleUpdate(A::fPEPS, col::Int, next_col::Int, H; kwargs...)::fPEPS
    do_side::Bool = get(kwargs, :do_side, true)
    τ::Float64    = get(kwargs, :tau, -0.1)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    @inbounds for row in Iterators.reverse(1:Ny)
        if do_side
            hori_col   = next_col < col ? next_col : col
            nhori_col  = next_col < col ? col : next_col
            si_a       = firstind(A[row, col], "Site")
            si_b       = firstind(A[row, next_col], "Site")
            ci         = commonind(A[row, col], A[row, next_col])
            min_dim    = ITensors.dim(ci)
            Ua, Sa, Va = svd(A[row, col], si_a, ci; mindim=min_dim, kwargs...)
            Ub, Sb, Vb = svd(A[row, next_col], si_b, ci; mindim=min_dim, kwargs...)
            Hab_hori   = is_cu ? cuITensor() : ITensor()
            horiH      = getDirectional(vcat(H[:, hori_col]...), Horizontal)
            horiH      = filter(x->x.sites[1][1] == row, horiH)
            for hH in horiH
                op_a = replaceind!(copy(hH.ops[1]), hH.site_ind, hori_col == col ? si_a : si_b)
                op_a = replaceind!(op_a, hH.site_ind', hori_col == col ? si_a' : si_b')
                op_b = replaceind!(copy(hH.ops[2]), hH.site_ind, hori_col == col ? si_b : si_a)
                op_b = replaceind!(op_b, hH.site_ind', hori_col == col ? si_b' : si_a')
                Hab_hori = ITensors.dim(Hab_hori) < 2 ? op_a * op_b : Hab_hori + op_a * op_b
            end
            cmb       = combiner(findinds(Hab_hori, 0), tags="hab,Site")
            ci        = combinedind(cmb)
            Hab_hori *= cmb
            Hab_hori *= cmb'
            Hab_mat   = is_cu ? matrix(collect(Hab_hori)) : matrix(Hab_hori)
            expiH_mat = exp(τ*Hab_mat)
            expiH     = is_cu ? cuITensor(vec(expiH_mat), ci, ci') : itensor(expiH_mat, ci, ci')
            expiH *= cmb
            expiH *= cmb'

            bond  = noprime(expiH * Ua * Ub)
            Uf, Sf, Vf = svd(bond, si_a, commonind(Ua, Sa); vtags="r,Link,r$row,c$hori_col", mindim=min_dim, kwargs...)
            A[row, col] = Sa * Va * Uf * Sf
            A[row, next_col] = Sb * Vb * Vf
        end
        if row < Ny
            si_a       = firstind(A[row, col], "Site")
            si_b       = firstind(A[row+1, col], "Site")
            ci         = commonind(A[row, col], A[row+1, col])
            min_dim    = ITensors.dim(ci)
            Ua, Sa, Va = svd(A[row, col], si_a, ci; mindim=min_dim, kwargs...)
            Ub, Sb, Vb = svd(A[row+1, col], si_b, ci; mindim=min_dim, kwargs...)
            Hab_vert   = is_cu ? cuITensor() : ITensor()
            vertH      = getDirectional(vcat(H[:, col]...), Vertical)
            vertH      = filter(x->x.sites[1][1] == row, vertH)
            for vH in vertH
                op_a = replaceind!(copy(vH.ops[1]), vH.site_ind, si_a)
                op_a = replaceind!(op_a, vH.site_ind', si_a')
                op_b = replaceind!(copy(vH.ops[2]), vH.site_ind, si_b)
                op_b = replaceind!(op_b, vH.site_ind', si_b')
                Hab_vert = ITensors.dim(Hab_vert) < 2 ? op_a * op_b : Hab_vert + op_a * op_b
            end
            cmb       = combiner(findinds(Hab_vert, 0), tags="hab,Site")
            ci        = combinedind(cmb)
            Hab_vert *= cmb
            Hab_vert *= cmb'
            Hab_mat   = is_cu ? matrix(collect(Hab_vert)) : matrix(Hab_vert)
            expiH_mat = exp(τ*Hab_mat)
            expiH     = is_cu ? cuITensor(vec(expiH_mat), ci, ci') : itensor(expiH_mat, ci, ci')
            expiH *= cmb
            expiH *= cmb'
            bond  = noprime(expiH * Ua * Ub)
            Uf, Sf, Vf = svd(bond, si_a, commonind(Ua, Sa); vtags="u,Link,r$row,c$col", mindim=min_dim, kwargs...)
            A[row, col] = Sa * Va * Uf * Sf
            A[row+1, col] = Sb * Vb * Vf
        end
    end
    return A
end

function measureEnergy(A::fPEPS,
                       L::Environments, R::Environments,
                       AncEnvs, H,
                       row::Int, col::Int)::Tuple{Float64, Float64}
    Ny, Nx    = size(A)
    Hs, N     = buildLocalH(A, L, R, AncEnvs, H, row, col, A[row, col])
    initial_N = collect(N * dag(A[row, col])')
    localH    = sum(Hs)
    initial_E = collect(localH * dag(A[row, col])')
    return real(scalar(initial_N)), real(scalar(initial_E))/real(scalar(initial_N))
end

function rightwardSweep(A::fPEPS,
                        Ls::Vector{Environments},
                        Rs::Vector{Environments},
                        H;
                        kwargs...)
    simple_update_cutoff = get(kwargs, :simple_update_cutoff, 4)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    dummyI = is_cu ? MPS([cuITensor(1.0) for ii in 1:Ny], 0, Ny+1) : MPS([ITensor(1.0) for ii in 1:Ny], 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny), fill(ITensor(), 1, Ny))
    sweep::Int = get(kwargs, :sweep, 0)
    sweep_width::Int = get(kwargs, :sweep_width, Nx)
    offset    = mod(Nx, 2)
    midpoint  = div(Nx, 2)
    rightmost = sweep_width == Nx ? Nx - 1 : midpoint + div(sweep_width, 2) + offset
    leftmost  = sweep_width == Nx ? 1 : midpoint - div(sweep_width, 2)
    ts_hori::Bool = get(kwargs, :two_site_hori, false)
    if ts_hori
        @inbounds for col in leftmost:rightmost
            next_col = col + 1
            L = col == 1 ? dummyEnv : Ls[col - 1]
            R = next_col == Nx ? dummyEnv : Rs[col + 2]
            @debug "Sweeping col $col"
            if sweep >= simple_update_cutoff
                @timeit "sweep" begin
                    #A = sweepColumnHorizontal(A, L, R, H, col, next_col; kwargs...)
                    A, Aj = sweepColumnHorizontal(A, L, R, H, col, next_col; kwargs...)
                end
            end
            if sweep < simple_update_cutoff
                # Simple update...
                A = simpleUpdate(A, col, col+1, H; do_side=(col < Nx), kwargs...)
            end
            if sweep >= simple_update_cutoff
                # Gauge
                #A = gaugeColumn(A, col, :right; kwargs...)
                A = gaugeColumnHorizontal(A, Aj, col, next_col, :right; kwargs...)
            end
            if col == 1
                left_H_terms = getDirectional(vcat(H[:, 1]...), Horizontal)
                @timeit "left edge env" begin
                    Ls[col] = buildEdgeEnvironment(A, H, left_H_terms, :left, 1; kwargs...)
                end
            else
                @timeit "left next env" begin
                    Ls[col] = buildNextEnvironment(A, Ls[col-1], H, :left, col; kwargs...)
                end
            end
        end
        Ls[Nx-1] = buildNextEnvironment(A, Ls[Nx-2], H, :left, Nx-1; kwargs...)
        @timeit "intraColumnGauge" begin
            A = intraColumnGauge(A, Nx; kwargs...)
        end
        EAncEnvs = buildAncs(A, Ls[Nx - 1], dummyEnv, H, Nx)
        N, E = measureEnergy(A, Ls[Nx - 1], dummyEnv, EAncEnvs, H, 1, Nx)
        println("Energy at Nx: ", E/(Nx*Ny), " and norm: ", N)
        println()
        println()
    else
        @inbounds for col in leftmost:rightmost
            L = col == 1 ? dummyEnv : Ls[col - 1]
            @debug "Sweeping col $col"
            if sweep >= simple_update_cutoff
                @timeit "sweep" begin
                    A = sweepColumn(A, L, Rs[col+1], H, col; kwargs...)
                end
            end
            if sweep < simple_update_cutoff
                # Simple update...
                A = simpleUpdate(A, col, col+1, H; do_side=(col < Nx), kwargs...)
            end
            if sweep >= simple_update_cutoff
                # Gauge
                A = gaugeColumn(A, col, :right; kwargs...)
            end
            if col == 1
                left_H_terms = getDirectional(vcat(H[:, 1]...), Horizontal)
                @timeit "left edge env" begin
                    Ls[col] = buildEdgeEnvironment(A, H, left_H_terms, :left, 1; kwargs...)
                end
            else
                @timeit "left next env" begin
                    Ls[col] = buildNextEnvironment(A, Ls[col-1], H, :left, col; kwargs...)
                end
            end
        end
    end
    return A, Ls, Rs
end

function leftwardSweep(A::fPEPS,
                       Ls::Vector{Environments},
                       Rs::Vector{Environments},
                       H;
                       kwargs...)
    simple_update_cutoff = get(kwargs, :simple_update_cutoff, 4)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    dummyI = is_cu ? MPS([cuITensor(1.0) for ii in 1:Ny], 0, Ny+1) : MPS([ITensor(1.0) for ii in 1:Ny], 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny), fill(ITensor(), 1, Ny))
    sweep::Int = get(kwargs, :sweep, 0)
    sweep_width::Int = get(kwargs, :sweep_width, Nx)
    offset    = mod(Nx, 2)
    midpoint  = div(Nx, 2)
    rightmost = midpoint + div(sweep_width, 2) + offset
    leftmost  = sweep_width == Nx ? 2 : midpoint - div(sweep_width, 2)
    ts_hori::Bool = get(kwargs, :two_site_hori, false)
    if ts_hori
        @inbounds for col in reverse(leftmost:rightmost)
            next_col = col - 1
            R = col == Nx    ? dummyEnv    : Rs[col + 1]
            L = next_col > 1 ? Ls[col - 2] : dummyEnv
            @debug "Sweeping col $col and next col $(col - 1)"
            if sweep >= simple_update_cutoff
                @timeit "sweep" begin
                    A, Aj = sweepColumnHorizontal(A, L, R, H, col, next_col; kwargs...)
                    #A = sweepColumnHorizontal(A, L, R, H, col, next_col; kwargs...)
                end
            end
            if sweep < simple_update_cutoff
                # Simple update...
                A = simpleUpdate(A, col, col-1, H; do_side=(col>1), kwargs...)
            end
            if sweep >= simple_update_cutoff
                # Gauge
                A = gaugeColumnHorizontal(A, Aj, col, next_col, :left; kwargs...)
                #A = gaugeColumn(A, col, :left; kwargs...)
            end
            if col == Nx
                right_H_terms  = getDirectional(vcat(H[:, Nx - 1]...), Horizontal)
                @timeit "right edge env" begin
                    Rs[col] = buildEdgeEnvironment(A, H, right_H_terms, :right, Nx; kwargs...)
                end
            else
                @timeit "right next env" begin
                    Rs[col] = buildNextEnvironment(A, Rs[col+1], H, :right, col; kwargs...)
                end
            end
        end
        Rs[2] = buildNextEnvironment(A, Rs[3], H, :right, 2; kwargs...)
        @timeit "intraColumnGauge" begin
            A = intraColumnGauge(A, 1; kwargs...)
        end
        EAncEnvs = buildAncs(A, dummyEnv, Rs[2], H, 1)
        N, E = measureEnergy(A, dummyEnv, Rs[2], EAncEnvs, H, 1, 1)
        println("Energy at 1: ", E/(Nx*Ny), " and norm: ", N)
        println()
        println()
    else
        @inbounds for col in reverse(leftmost:rightmost)
            R = col == Nx ? dummyEnv : Rs[col + 1]
            @debug "Sweeping col $col"
            if sweep >= simple_update_cutoff
                @timeit "sweep" begin
                    A = sweepColumn(A, Ls[col - 1], R, H, col; kwargs...)
                end
            end
            if sweep < simple_update_cutoff
                # Simple update...
                A = simpleUpdate(A, col, col-1, H; do_side=(col>1), kwargs...)
            end
            if sweep >= simple_update_cutoff
                # Gauge
                A = gaugeColumn(A, col, :left; kwargs...)
            end
            if col == Nx
                right_H_terms  = getDirectional(vcat(H[:, Nx - 1]...), Horizontal)
                @timeit "right edge env" begin
                    Rs[col] = buildEdgeEnvironment(A, H, right_H_terms, :right, Nx; kwargs...)
                end
            else
                @timeit "right next env" begin
                    Rs[col] = buildNextEnvironment(A, Rs[col+1], H, :right, col; kwargs...)
                end
            end
        end
    end
    return A, Ls, Rs
end

function doSweeps(A::fPEPS,
                  Ls::Vector{Environments},
                  Rs::Vector{Environments},
                  H;
                  mindim::Int=1,
                  maxdim::Int=1,
                  simple_update_cutoff::Int=4,
                  sweep_start::Int=1,
                  sweep_count::Int=10,
                  cutoff::Float64=0.,
                  env_maxdim=2maxdim,
                  do_mag::Bool=false,
                  prefix="mag",
                  max_gauge_iter::Int=100,
                  model::Symbol=:XXZ,
                  two_site_vert::Bool=false,
                  two_site_hori::Bool=false,
                  rotation::Bool=false)
    Ny, Nx = size(A)
    for sweep in sweep_start:sweep_count
        overlap_cutoff = 0.999
        if iseven(sweep)
            (A, Ls, Rs), this_time, bytes, gctime, memallocs = @timed rightwardSweep(A, Ls, Rs, H; sweep=sweep, mindim=mindim, maxdim=maxdim, simple_update_cutoff=simple_update_cutoff, overlap_cutoff=overlap_cutoff, cutoff=cutoff, env_maxdim=env_maxdim, max_gauge_iter=max_gauge_iter, two_site_vert=two_site_vert, two_site_hori=two_site_hori)
            println("SWEEP RIGHT $sweep, time $this_time")
        else
            (A, Ls, Rs), this_time, bytes, gctime, memallocs = @timed leftwardSweep(A, Ls, Rs, H; sweep=sweep, mindim=mindim, maxdim=maxdim, simple_update_cutoff=simple_update_cutoff, overlap_cutoff=overlap_cutoff, cutoff=cutoff, env_maxdim=env_maxdim, max_gauge_iter=max_gauge_iter, two_site_vert=two_site_vert, two_site_hori=two_site_hori)
            println("SWEEP LEFT $sweep, time $this_time")
        end
        flush(stdout)
        if sweep == simple_update_cutoff - 1
            for col in reverse(2:Nx)
                A = gaugeColumn(A, col, :left; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim, max_gauge_iter=max_gauge_iter)
            end
            Ls = buildLs(A, H; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
            Rs = buildRs(A, H; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
        end
        # make this dispatch at some point
        if do_mag && model == :XXZ
            measure_correlators_heisenberg(copy(A), H, Ls, Rs, sweep;  mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim, prefix=prefix)
        elseif do_mag && model == :Ising
            measure_correlators_ising(copy(A), H, Ls, Rs, sweep; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim, prefix=prefix)
        end
        if rotation && mod(sweep, 10) == 0
            new_A = copy(A)
            for row in 1:Ny, col in 1:Nx
                new_A[col, row] = A[row, col]
            end
            Ls = buildLs(A, H; mindim=maxdim, maxdim=maxdim, cutoff=cutoff, env_maxdim=env_maxdim)
        end
    end
    return A
end


text:
lower(x::Dict{Symbol}) = BSONDict(x)

# Basic Types

ismutable(::Type{Symbol}) = false
lower(x::Symbol) = BSONDict(:tag => "symbol", :name => String(x))

lower(x::Tuple) = BSONDict(:tag => "tuple", :data => Any[x...])

ismutable(::Type{SimpleVector}) = false
lower(x::SimpleVector) = BSONDict(:tag => "svec", :data => Any[x...])

# References

ref(path::Symbol...) = BSONDict(:tag => "ref", :path => Base.string.([path...]))

resolve(fs) = reduce((m, f) -> getfield(m, Symbol(f)), fs; init = Main)

function modpath(x::Module)
  y = parentmodule(x)
  x == y ? [nameof(x)] : [modpath(y)..., nameof(x)]
end

ismutable(::Type{Module}) = false
lower(m::Module) = ref(modpath(m)...)

# Types

ismutable(::Type{<:Type}) = false

typepath(x::DataType) = [modpath(x.name.module)..., x.name.name]

function lower(v::DataType)
  BSONDict(:tag => "datatype",
           :name => Base.string.(typepath(v)),
           :params => [v.parameters...])
end

constructtype(T, Ts) = isempty(Ts) ? T : T{Ts...}
constructtype(T::Type{Tuple}, Ts) = T{Ts...}

lower(v::UnionAll) =
  BSONDict(:tag => "unionall",
           :body => v.body,
           :var => v.var)

# Arrays

lower(x::Vector{Any}) = copy(x)
lower(x::Vector{UInt8}) = x

function reinterpret_(::Type{T}, x) where T
  length(x) < 1 && return T[]

  len = Int(length(x) * (sizeof(eltype(x)) / sizeof(T)))

  GC.@preserve x begin
    p = Ptr{UInt8}(pointer(x))
    t = Ptr{UInt8}(Base.Libc.calloc(len, sizeof(T)))
    ccall(:memset, Cvoid, (Ptr{UInt8}, Cint, Csize_t), t, 0, len * sizeof(T))
    ccall(:memcpy, Cvoid, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), t, p, sizeof(x))
    Base.unsafe_wrap(Array, Ptr{T}(t), len; own=true)
  end
end

function lower(x::Array)
  ndims(x) == 1 && !isbitstype(eltype(x)) && return Any[x...]
  BSONDict(:tag => "array", :type => eltype(x), :size => Any[size(x)...],
           :data => isbitstype(eltype(x)) ? reinterpret_(UInt8, reshape(x, :)) : Any[x...])
end

# Structs

isprimitive(T::Type)::Bool = fieldcount(T) == 0 && T.size > 0

structdata(x) = isprimitive(typeof(x)) ? reinterpret_(UInt8, [x]) :
    Any[getfield(x, f) for f in fieldnames(typeof(x))]

function lower(x)
  BSONDict(:tag => "struct", :type => typeof(x), :data => structdata(x))
end

initstruct(T) = ccall(:jl_new_struct_uninit, Any, (Any,), T)

function newstruct!(x, fs...)
  #@debug "newstruct!" typeof(x) fs
  for (i, f) = enumerate(fs)
    f = convert(fieldtype(typeof(x),i), f)
    ccall(:jl_set_nth_field, Nothing, (Any, Csize_t, Any), x, i-1, f)
  end
  return x
end

function newstruct(T::Type, xs...)
  if !T.mutable
    flds = Any[convert(fieldtype(T, i), x) for (i,x) in enumerate(xs)]
    return ccall(:jl_new_structv, Any, (Any,Ptr{Cvoid},UInt32), T, flds, length(flds))
  else
    # Manual inline of newstruct! to work around bug
    # https://github.com/MikeInnes/BSON.jl/issues/2#issuecomment-452204339
    x = initstruct(T)

    for (i, f) = enumerate(xs)
      f = convert(fieldtype(typeof(x),i), f)
      ccall(:jl_set_nth_field, Nothing, (Any, Csize_t, Any), x, i-1, f)
    end
    x

  end
end

function newstruct_raw(cache::IdDict{Any, Any}, T::Type, d::TaggedStruct)
  @assert isstructtype(T) "$T is not struct type"
  #@debug "newstruct_raw" T d

  x = cache[d] = initstruct(T)
  fs = (_raise_recursive(x, cache) for x in d.data)

  newstruct!(x, fs...)
end

newprimitive(T, data) = reinterpret_(T, data)[1]

iscyclic(T) = ismutable(T)

lower(v::Type{Union{}}) = BSONDict(:tag=>"jl_bottom_type")

# Base data structures

structdata(d::Dict) = Any[collect(keys(d)), collect(values(d))]

initstruct(D::Type{<:Dict}) = D()

function newstruct!(d::Dict, ks, vs)
  for (k, v) in zip(ks, vs)
    d[k] = v
  end
  return d
end


text:
using Bamberg
using Base.Test, StaticArrays

for n in 1000:1000:10000
    gridpoints_fa = generateidealdata(n,"flat")
    gridpoints_na = generateidealdata(n,"nested")
    gridpoints_st = generateidealdata(n,"static-nested")
    @testset "Number of points: $n" begin
        @test findminnorm(gridpoints_fa, [n-1 1]')[2]  == n-1
        @test fnp(gridpoints_fa, [n-1 1]')[2]  == n-1
        @test pointsvecnorm(gridpoints_na, [n-1 1]')[2] == n-1
        @test pointsvecnorm(gridpoints_st, SVector{2, Int}(n-1, 1))[2] == n-1
        @test pointsvecnorm_opt(gridpoints_st, SVector{2, Int}(n-1, 1))[2] == n-1
    end
end

max_stars_repo_name: Affos/Bamberg

max_issues_repo_name: Affos/Bamberg

max_forks_repo_name: Affos/Bamberg


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Dimension independent code related to containers of the DG solver
# with the mesh type TreeMesh

function reinitialize_containers!(mesh::TreeMesh, equations, dg::DGSEM, cache)
  # Get new list of leaf cells
  leaf_cell_ids = local_leaf_cells(mesh.tree)

  # re-initialize elements container
  @unpack elements = cache
  resize!(elements, length(leaf_cell_ids))
  init_elements!(elements, leaf_cell_ids, mesh, dg.basis)

  # re-initialize interfaces container
  @unpack interfaces = cache
  resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids))
  init_interfaces!(interfaces, elements, mesh)

  # re-initialize boundaries container
  @unpack boundaries = cache
  resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids))
  init_boundaries!(boundaries, elements, mesh)

  # re-initialize mortars container
  @unpack mortars = cache
  resize!(mortars, count_required_mortars(mesh, leaf_cell_ids))
  init_mortars!(mortars, elements, mesh)

  if mpi_isparallel()
    # re-initialize mpi_interfaces container
    @unpack mpi_interfaces = cache
    resize!(mpi_interfaces, count_required_mpi_interfaces(mesh, leaf_cell_ids))
    init_mpi_interfaces!(mpi_interfaces, elements, mesh)

    # re-initialize mpi_mortars container
    @unpack mpi_mortars = cache
    resize!(mpi_mortars, count_required_mpi_mortars(mesh, leaf_cell_ids))
    init_mpi_mortars!(mpi_mortars, elements, mesh)

    # re-initialize mpi cache
    @unpack mpi_cache = cache
    init_mpi_cache!(mpi_cache, mesh, elements, mpi_interfaces, mpi_mortars,
                    nvariables(equations), nnodes(dg))
  end
end

# Dimension-specific implementations
include("containers_1d.jl")
include("containers_2d.jl")
include("containers_3d.jl")

end # @muladd


text:
import Laplacian
import p3

# Generate a clean Hamiltonian
function q3ham(N::Integer, q::Real)

    # Boundary conditions as matrices
    Ψb, B = p3.boundaryconditions3(N)

    # Reshape conditions into vectors
    ψb = reshape(Ψb, (N+1)^2)
    β  = reshape(B,  (N+1)^2)
    V  = reshape(p3.V(N,q), (N+1)^2)

    # Hamiltonian matrix
    H = -5Laplacian.lap(N) - V.*eye((N+1)^2)
    Hf, Hb = Laplacian.opsplit(H, β)
    He = Laplacian.deltaeff(Hf, β)
    Hm = p3.minor(He, β)
    return Hm, He, ψb, β

end

function q3groundeig(Hm::Matrix, steps::Integer, tolerance::Real)

    v = rand(size(Hm)[1])
    v /= norm(v)

    # using the minor matrix kill most zero-valued eigenvalues
    T, K = Lanczos.lanczos(Hm,v,steps)

    # find smallest (nonzero) eigenvector
    D, V = eig(T)
    k = findnext(x -> >(x,1.0), D, 1)

    # old eigenvalue
    λ = D[k]
    err = 100.0

    # best guess at small eigenvalue
    vold = v
    v = K*V[:,k]
    v /= norm(v)

    while err > tolerance
        T, K = Lanczos.lanczos(Hm,v,steps)
        D, V = eig(T)
        k = findnext(x -> >(x,1.0), D, 1)

        # how far off our eigenvalue guess is from the real thing
        err = norm((Hm - λ*eye(size(Hm)[1]))*v)
        λ = D[k]                    # smallest nonzero eigenvalue

        vold = v
        v = K*V[:,k]
        v /= norm(v)

    end

    # return the eigenvalue and eigenvector
    return λ, v

end

# Get the ground state energy and wavefunction, with the latter as a matrix
function q3groundpsi(N::Integer, q::Real, steps::Integer, tolerance::Real)
    Hm, He, ψb, β = q3ham(N, q)
    E0, ψ0 = q3groundeig(Hm, steps, tolerance)
    return E0, reshape(p3.addwater(ψ0,β) + ψb, N+1, N+1)
end

q3groundpsi(N::Integer, q::Real) = q3groundpsi(N, q, 10, 1e-3)

# Return ground state energies and wavefunctions for different q values
function q3groundpsi(N::Integer, qs::Vector)

    Es = Array(Float64, length(qs))
    ψs = Array(Float64, N+1, N+1, length(qs))

    @sync @parallel for i in 1:length(qs)
        Es[i], ψs[:,:,i] = q3groundpsi(N,qs[i])
    end

    return Es, ψs

end

# Find probability that we're in the lower half
function q3P(N::Integer, q::Real, steps::Integer, tolerance::Real)
    E0, ψ0 = q3groundpsi(N, q, steps, tolerance)
    r = ψ0[:,1:int(floor(N/2+1))]   # just bottom half
    return sum(r.*r)                # ψ*ψ
end

# Something we can feed to an equation solver
q3P(q::Real) = q3P(32, q, 20, 1e-3)


text:
# Romeo LTI frequency responses.
#
# Implementation of `freqresp` function of various systems.

freq_resp(sys::Diff{<:Number}, ω::AbstractVector{<:Number}) = (1im * ω) .^ sys.α

export freq_resp


text:
using Printf
include("util.jl")

lines = getlines("inputs/7.txt")

function parse_step(line::String)::Pair{Char,Char}
    matches = match(r"Step (.) must be finished before step (.) can begin.", line)
    x = matches[1][1]
    y = matches[2][1]
    return x => y
end

function get_steps(lines::Array{String})
    steps::Dict{Char,Array{Char}} = Dict()
    letters::Set{Char} = Set()
    for line in lines
        p = parse_step(line)
        push!(letters, p[1])
        push!(letters, p[2])
        step = get(steps, p[2], [])
        push!(step, p[1])
        steps[p[2]] = step
    end
    for l in letters
        if !(l in keys(steps))
            steps[l] = []
        end
    end
    return steps
end

function and(arr::Array{Bool})::Bool
    for a in arr
        if !a
            return false
        end
    end
    return true
end

function ready(step::Char, steps::Dict{Char,Array{Char}}, runsteps::Array{Char}, takensteps::Array{Char})::Bool
    if step in runsteps
        return false
    elseif step in takensteps
        return false
    else
        return and([(s in runsteps) for s in steps[step]])
    end
end

function run_steps(steps::Dict{Char,Array{Char}})
    runsteps::Array{Char} = []
    sl = sort(collect(keys(steps)))
    println(sl)
    while length(runsteps) < length(sl)
        ts::Array{Char} = []
        readysteps = filter(x->ready(x, steps, runsteps, ts), sl)
        next_step = sort(readysteps)[1]
        push!(runsteps, next_step)
    end
    return runsteps
end

steps = get_steps(lines)
part_1 = run_steps(steps)

mutable struct Worker
    timeleft::Int64
    letter::Char
end

function do_work(steps::Dict{Char,Array{Char}})
    w() = Worker(0,'.')
    workers = [w(),w(),w(),w(),w()]
    base_time = 61
    sl = sort(collect(keys(steps)))
    times::Dict{Char,Int64} = Dict()
    for l in sl
        times[l] = base_time
        base_time += 1
    end

    runsteps::Array{Char} = []
    takensteps::Array{Char} = []
    timetaken = 0
    while length(runsteps) < length(sl)
        for worker in workers
            if worker.timeleft > 1
                worker.timeleft -= 1
            elseif worker.timeleft <= 1
                push!(runsteps,worker.letter)
                readysteps = filter(x->ready(x, steps, runsteps, takensteps), sl)
                if length(readysteps) > 0
                    next_step = sort(readysteps)[1]
                    worker.letter = next_step
                    push!(takensteps, next_step)
                    worker.timeleft = times[next_step]
                else
                    worker.letter = '.'
                end
            end
        end
        s = @sprintf "%s %s %s %s %s %s\n" timetaken workers[1].letter workers[2].letter workers[3].letter workers[4].letter workers[5].letter;
        print(s)
        runsteps = filter(x->x != '.', runsteps)
        timetaken += 1
    end
    println(runsteps)
    println(timetaken-1)
end

do_work(steps)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Base64URL

# Base64EncodePipe is a pipe-like IO object, which converts into base64 data
# sent to a stream. (You must close the pipe to complete the encode, separate
# from closing the target stream).  We also have a function base64encode(f,
# args...) which works like sprint except that it produces base64-encoded data,
# along with base64encode(args...)  which is equivalent to base64encode(write,
# args...), to return base64 strings.  A Base64DecodePipe object can be used to
# decode base64-encoded data read from a stream , while function base64decode is
# useful for decoding strings

include("buffer.jl")
include("encode.jl")
include("decode.jl")

end


text:
using Yao, YaoExtensions
using Test

@testset "SqrtX" begin
    @test mat(SqrtX)^2 ≈ mat(X)
    @test mat(SqrtY)^2 ≈ mat(Y)
    @test cz_entangler(5, [1=>2, 4=>5]) == chain(control(5, 1, 2=>Z), control(5, 4, 5=>Z))
end

@testset "supremacy" begin
    c = rand_supremacy2d(4, 4, 8)
    @test length(c) == 8
    @test length(pair_supremacy(6,6)) == 8
    @test pair_supremacy(6,6)[1] == [7=>8, 19=>20, 31=>32, 3=>4, 15=>16, 27=>28, 11=>12, 23=>24, 35=>36]
end


text:
module PlinkPlumber
	include("plink.jl")
end # module


text:
export l2Norm_project

"""
Force the rows of L and R to have norm at most B
"""
function l2Norm_project(x::Tx, B, weights, params::Dict{String,<:Any}) where {ETx<:Number, Tx<:AbstractVector{ETx}}

if isapprox(norm(x), 0.0)
    println("WARNING: norm(x) cannot be 0 in l2Norm_project")
end

c = B/norm(x)
xout = min(1,c).*x

#Dummy
itn = 1

return xout, itn

end


text:
using Printf

mutable struct GameLayer <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    # Collection of nodes.
    children::Array{Ranger.AbstractNode,1}

    min::Geometry.Point{Float64}
    max::Geometry.Point{Float64}
    buc_min::Geometry.Point{Float64}
    buc_max::Geometry.Point{Float64}

    orbit_system::OrbitSystemNode
    yellow_rect::Custom.OutlinedRectangle
    circle::Custom.OutlinedCircle
    ship::Ship

    zoom::Custom.ZoomNode

    function GameLayer(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.children = Array{Ranger.AbstractNode,1}[]
        o.transform = Nodes.TransformProperties{Float64}()

        o.min = Geometry.Point{Float64}()
        o.max = Geometry.Point{Float64}()
        o.buc_min = Geometry.Point{Float64}()
        o.buc_max = Geometry.Point{Float64}()

        o
    end
end

function build(layer::GameLayer, world::Ranger.World)
    hw = Float64(world.view_width) / 2.0
    hh = Float64(world.view_height) / 2.0

    # top-left
    Geometry.set!(layer.min, -hw, -hh)

    # bottom-right
    Geometry.set!(layer.max, hw, hh)

    layer.orbit_system = OrbitSystemNode(world, "OrbitSystemNode", layer)
    build(layer.orbit_system, world)
    set!(layer.orbit_system, -0.5, -0.5, 0.5, 0.5)
    Nodes.set_rotation_in_degrees!(layer.orbit_system, 30.0)
    Nodes.set_scale!(layer.orbit_system, 100.0)
    Nodes.set_position!(layer.orbit_system, 100.0, -100.0)
    layer.orbit_system.color = RangerGame.olive
    push!(layer.children, layer.orbit_system);

    # Zoom Testing -----------------------------------------------------
    layer.zoom = Custom.ZoomNode(world, "ZoomNode", layer)

    rect = Custom.OutlinedRectangle(world, "YellowOutlinedRectangle", layer.zoom)
    layer.yellow_rect = rect
    Custom.set!(rect, -0.5, -0.5, 0.5, 0.5) # centered
    Nodes.set_scale!(rect, 200.0)
    Nodes.set_position!(rect, -300.0, 300.0)
    Nodes.set_rotation_in_degrees!(rect, 30.0)
    rect.color = RangerGame.yellow
    push!(layer.zoom.children, rect);

    circle = Custom.OutlinedCircle(world, "OutlinedCircle", layer.zoom)
    layer.circle = circle
    Custom.set!(circle, 18.0)
    Nodes.set_scale!(circle, 100.0)
    Nodes.set_position!(circle, 300.0, 300.0)
    circle.color = RangerGame.lightpurple
    push!(layer.zoom.children, circle);

    push!(layer.children, layer.zoom);
    # ------------------------------------------------------------

    layer.ship = Ship(world, "Ship", layer)
    Nodes.set_scale!(layer.ship, 30.0)
    set_position!(layer.ship, -100.0, -100.0)
    set_direction!(layer.ship, 30.0)
    layer.ship.color = RangerGame.peach
    push!(layer.children, layer.ship);
end

# --------------------------------------------------------
# Timing
# --------------------------------------------------------
function Nodes.update(layer::GameLayer, dt::Float64)
    Nodes.update(layer.yellow_rect, dt)
    Nodes.update(layer.circle, dt)
end

# --------------------------------------------------------
# Rendering
# --------------------------------------------------------
function Nodes.draw(layer::GameLayer, context::Rendering.RenderContext)
    if Nodes.is_dirty(layer)
        Rendering.transform!(context, layer.min, layer.max, layer.buc_min, layer.buc_max)
        Nodes.set_dirty!(layer, false)
    end

    # Render background (i.e. a solid gray color)
    Rendering.set_draw_color(context, RangerGame.darkgray)
    Rendering.render_aa_rectangle(context, layer.buc_min, layer.buc_max, Rendering.FILLED)

    # Draw debug layer name
    Rendering.set_draw_color(context, RangerGame.white)
    Rendering.draw_text(context, 10, 10, layer.base.name, 2, 2, false)
end

# --------------------------------------------------------
# Life cycle events
# --------------------------------------------------------
function Nodes.enter_node(layer::GameLayer, man::Nodes.NodeManager)
    println("enter ", layer)
    # Register node as a timing target in order to receive updates
    Nodes.register_target(man, layer)
    Nodes.register_target(man, layer.ship)

    # Register for io events such as keyboard or mouse
    Nodes.register_event_target(man, layer)
    Nodes.register_event_target(man, layer.ship)
    Nodes.register_event_target(man, layer.zoom);
end

function Nodes.exit_node(layer::GameLayer, man::Nodes.NodeManager)
    println("exit ", layer)
    Nodes.unregister_target(man, layer)
    Nodes.unregister_target(man, layer.ship)

    Nodes.unregister_event_target(man, layer)
    Nodes.unregister_event_target(man, layer.ship);
    Nodes.unregister_event_target(man, layer.zoom);
end

# --------------------------------------------------------
# Events
# --------------------------------------------------------
function Nodes.io_event(node::GameLayer, event::Events.MouseMotionEvent)
    Nodes.io_event(node.yellow_rect, event)
    Nodes.io_event(node.circle, event)
end

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::GameLayer)
    node.children
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Metatheory
using Metatheory.Classic
using Metatheory.EGraphs
using Metatheory.Library
using Metatheory.Util

using DataStructures
using Test

# Metatheory.options[:verbose] = true
# Metatheory.options[:printiter] = true

@metatheory_init

Test.FallbackTestSet(desc) = Test.FallbackTestSet()

ts = Test.FallbackTestSet

@testset ts "Metatheory Tests" begin
   @timev begin
      include("test_egraphs.jl")
      include("test_ematch.jl")
      include("test_ematch_assertions.jl")
      include("test_analysis.jl")
      include("test_extraction.jl")
      include("test_dynamic_ematch.jl")
      # FIXME include("test_mu.jl")
      include("test_boson.jl")
      include("test_while_interpreter.jl")
      include("test_theories.jl")
      include("test_reductions.jl")
      include("test_taylor.jl")
      include("test_while_superinterpreter.jl")
      include("test_fibonacci.jl")
      include("test_cat.jl")
      include("test_calculational_logic.jl")
      include("test_logic.jl")
      include("test_inequality.jl")

      # exported consistency test
      for m ∈ [Metatheory, Metatheory.Util, Metatheory.Classic, Metatheory.EGraphs, Metatheory.EGraphs.Schedulers]
         for i ∈ propertynames(m)
            xxx = getproperty(m, i)
         end
      end
   end
end


text:
draw(PNG("square-scale.png", 140, 100),
     Paper("lightgrey"),
     Axes(centred=false),
     Scale(0.5)) do
    move(0, 0)
    square(1, align=7)
end


text:
using Printf
using Plots
using JULES
using Oceananigans
using Oceananigans.Advection

using Oceananigans.Grids: Cell, xnodes

ENV["GKSwstype"] = "100"

N = 512
L = 1
T = 0.25
Δt = 1e-3
Nt = Int(T/Δt)

topo = (Bounded, Periodic, Periodic)
domain = (x=(0, L), y=(0, 1), z=(0, 1))
grid = RegularCartesianGrid(topology=topo, size=(N, 1, 1), halo=(4, 4, 4); domain...)

model = CompressibleModel(
                      grid = grid,
                 advection = WENO5(),
                     gases = DryEarth(),
    thermodynamic_variable = Entropy(),
                   closure = IsotropicDiffusivity(ν=0, κ=0)
)

gas = model.gases.ρ
R, cₚ, cᵥ = gas.R, gas.cₚ, gas.cᵥ
u₀₀, T₀₀, ρ₀₀, s₀₀ = gas.u₀, gas.T₀, gas.ρ₀, gas.s₀
g  = model.gravity

ρₗ, ρᵣ = 1.0, 0.125
pₗ, pᵣ = 1.0, 0.1
uₗ, uᵣ = 0.0, 0.0

ρ₀(x, y, z) = x < 0.5 ? ρₗ : ρᵣ
p₀(x, y, z) = x < 0.5 ? pₗ : pᵣ
u₀(x, y, z) = x < 0.5 ? uₗ : uᵣ

T₀(x, y, z) = p₀(x, y, z) / (R*ρ₀(x, y, z))
ρs₀(x, y, z) = ρ₀(x, y, z) * (s₀₀ + cᵥ * log(T₀(x, y, z)/T₀₀) - R * log(ρ₀(x, y, z)/ρ₀₀))

set!(model.tracers.ρ, ρ₀)
set!(model.tracers.ρs, ρs₀)
update_total_density!(model)

anim = @animate for n in 1:Nt
    @info "iteration $n/$Nt"
    time_step!(model, Δt)

    title = @sprintf("Shock tube t=%.3f", model.clock.time)

    x = xnodes(Cell, grid)
    ρ = interior(model.total_density)[:]
    plot(x, ρ, lw=2, label="", title=title, xlims=(0, 1), ylims=(0, 1), dpi=200)
end

mp4(anim, "shock_tube.mp4", fps=60)


text:
using DataStructures

import Base: eof, close, position
export ContextMap, close, value, eof, advance!

type ContextMap
    reader
    contextBefore::Int64
    contextAfter::Int64
    position::Int64
    posQueue::Deque{Tuple{Int64,Float64}}
    value::Float64
end

function ContextMap(reader, contextBefore::Int64, contextAfter::Int64)
    cm = ContextMap(reader, contextBefore, contextAfter, 0, Deque{Tuple{Int64,Float64}}(), 0.0)

    advance!(cm)
    cm
end
close(cm::ContextMap) = close(cm.reader)
value(cm::ContextMap) = cm.value
position(cm::ContextMap) = cm.position
eof(cm::ContextMap) = cm.position < 0

function advance!(cm::ContextMap)
    cm.position += 1

    # remove old positions that fell outside the sliding window
    while length(cm.posQueue) > 0 && front(cm.posQueue)[1] < cm.position-cm.contextBefore
        cm.value -= front(cm.posQueue)[2]
        shift!(cm.posQueue)
    end

    # skip ahead to the next spot we have data if our queue is empty
    if length(cm.posQueue) == 0
        cm.position = position(cm.reader)-cm.contextBefore
    end

    # add new positions in the sliding window
    while !eof(cm.reader) && position(cm.reader) <= cm.position+cm.contextAfter
        v = value(cm.reader)
        cm.value += v
        push!(cm.posQueue, (position(cm.reader),v))
        advance!(cm.reader)
    end
end


text:
@testset "Select" begin
    @testset "Dense" begin
        #dense first
        X = GBMatrix(rand(-10.0:0.05:10.0, 10, 10))
        test_frule(map, UnaryOps.SIN, X)
        test_rrule(map, UnaryOps.SIN, X)
        test_frule(map, UnaryOps.COS, X)
        test_rrule(map, UnaryOps.COS, X)
        test_frule(map, UnaryOps.COSH, X)
        test_rrule(map, UnaryOps.COSH, X)
        test_frule(map, UnaryOps.SINH, X)
        test_rrule(map, UnaryOps.SINH, X)
        test_frule(map, UnaryOps.TANH, X)
        test_rrule(map, UnaryOps.TANH, X)
        test_frule(map, UnaryOps.MINV, X)
        test_rrule(map, UnaryOps.MINV, X)
        test_frule(map, UnaryOps.ABS, X)
        test_rrule(map, UnaryOps.ABS, X)
        test_frule(map, UnaryOps.IDENTITY, X)
        test_rrule(map, UnaryOps.IDENTITY, X)
        test_frule(map, UnaryOps.EXP, X)
        test_rrule(map, UnaryOps.EXP, X)
    end
    @testset "Sparse" begin
        X = GBMatrix(sprand(50, 50, 0.15))
    end
end


text:
"""
Copyright Marcus Greiff 2020

Tests for the power distribution functionality
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

# Testing fixture
function fixture()
    J  = rand(3,3); J  = J + J'; J  = J + I*maximum(abs.(eigvals(J)));
    kR, kc, kw = 1.0, 1.0, 1.0;
    a, b, c, d = 0.35, 0.26, 1.0+rand(), 1.0+rand();
    C  = con_state_qw_fsf_t(J, kR, kc, kw; a=a, b=b, c=c, d=d);
    C.thrust = 0.5;
    t = rand(3,1);
    C.torques = ntuple(i->Cdouble(t[i,1]), 3)
    Cc = deepcopy(C);
    d, dc, dm = get_matrix_and_copy(rand(4,1))
    return C, Cc, t, d, dc, dm;
end

@verbose(1, "Testing the power distribution...")
@testset "Test the power distribution" begin
    @testset "Good inputs" begin
        # Fixture
        C, Cc, t, d, dc, dm = fixture();

        apar = C.param_a;
        bpar = C.param_b;
        cpar = C.param_c;
        dpar = C.param_d;
        T = C.thrust;

        # Modified configuration for the razor drone
        A = [ +1.00  +1.00  +1.00  +1.00;
              -dpar  -dpar  +dpar  +dpar;
              -dpar  +dpar  +dpar  -dpar;
              +cpar  -cpar  +cpar  -cpar]

        #A = [+1    +1    +1    +1   ;
        #     +0    -dpar +0    +dpar;
        #     +dpar +0    -dpar +0   ;
        #     -cpar +cpar -cpar +cpar]

        f_julia = A \ ([T; t]);

        f_julia[f_julia .<0.0] .= 0.0

        d_julia = -(apar/(2.0*bpar)) .+ sqrt.((apar/(2*bpar)).^2 .+ f_julia./bpar);

        # Test
        status = ccall((:compute_power_distribution, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t}, Ref{matrix_double_t}),
            C,
            dm)

        # Assert
        @test isapprox(status, 1);
        for i = 1:4
            @test isapprox(d_julia[i,1], d[i,1]);
        end
    end
    @testset "Bad inputs (negative thrust)" begin
        # Fixture
        C, Cc, t, d, dc, dm = fixture();
        C.thrust = -C.thrust;
        Cc = deepcopy(C);

        # Test
        status = ccall((:compute_power_distribution, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t}, Ref{matrix_double_t}),
            C,
            dm)

        # Assert
        @test isapprox(status, 0);
    end
    @testset "Bad inputs (negative parameter value)" begin
        for i = 1:4
            # Fixture
            C, Cc, t, d, dc, dm = fixture();
            C.param_a = -C.param_a;
            Cc = deepcopy(C);

            # Test
            status = ccall((:compute_power_distribution, AerialVehicleControl.CONT_LIB_PATH),
                Cint,
                (Ref{con_state_qw_fsf_t}, Ref{matrix_double_t}),
                C,
                dm)

            # Assert
            @test isapprox(status, 0);
        end
    end
    @testset "Bad inputs (wrong duty cycle row dimensions)" begin
        # Fixture
        C, Cc, t, _, _, _ = fixture();
        d, dc, dm = get_matrix_and_copy(rand(3,1));

        # Test
        status = ccall((:compute_power_distribution, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t}, Ref{matrix_double_t}),
            C,
            dm)
        # Assert
        @test isapprox(status, 0);
    end
    @testset "Bad inputs (wrong duty cycle column dimensions)" begin
        # Fixture
        C, Cc, t, _, _, _ = fixture();
        d, dc, dm = get_matrix_and_copy(rand(4,2));

        # Test
        status = ccall((:compute_power_distribution, AerialVehicleControl.CONT_LIB_PATH),
            Cint,
            (Ref{con_state_qw_fsf_t}, Ref{matrix_double_t}),
            C,
            dm)
        # Assert
        @test isapprox(status, 0);
    end
end


text:
config = MOI.Test.TestConfig()

@testset "Model" begin
    optimize!(mock) = MOIU.mock_optimize!(mock, [0.0, 0.0],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1.0]],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0]])
    for mock in mocks(optimize!)
        Tests.sos_term_test(mock, config)
        Tests.sdsos_term_test(mock, config)
        Tests.dsos_term_test(mock, config)
    end
end

# The VectorOfVariables-in-SOSPolynomialSet is bridged by VectorFunctionizeBridge
# since the free variable is bridged. This tests that the GramMatrixAttribute, ...
# are passed by the VectorFunctionizeBridge.
@testset "NoFreeVariable" begin
    optimize!(mock) = MOIU.mock_optimize!(mock, [0.0, 0.0, 0.0],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1.0]],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.0]])
    nofree_mock = bridged_mock(optimize!, model = NoFreeVariable{Float64}())
    for mock in mocks(optimize!)
        Tests.sos_term_test(nofree_mock, config)
        Tests.sdsos_term_test(nofree_mock, config)
        Tests.dsos_term_test(nofree_mock, config)
    end
end


text:
using Multilane
using MCTS
using POMCP

using ArgParse

s = ArgParseSettings()

@add_arg_table s begin
    "dir"
        help = "data directory - will look here for results_list.txt"
        nargs = 1
        required = false
end

args = parse_args(ARGS, s)

if length(args["dir"]) > 0
    dir = first(args["dir"])
else
    parent = get(ENV, "SCRATCH", tempdir())
    dirs = readdir(parent)
    latest = 0.0
    latest_dir = ""
    for d in dirs
        abspath = joinpath(parent, d)
        if ctime(abspath) > latest
            latest = ctime(abspath)
            latest_dir = abspath
        end
    end
    dir = latest_dir
end

println("Looking at $(joinpath(dir, "results_list.txt"))")
files = []
open(joinpath(dir, "results_list.txt")) do list
    for f in eachline(list)
        push!(files, chomp(f))
    end
end

filename = string("combined_results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
results = gather_results(files, save_file=filename)
println("saved to $filename")


text:
abstract type ParametricModel <: AbstractYieldCurve end
Base.Broadcast.broadcastable(x::T) where {T<:ParametricModel} = Ref(x)
__ratetype(::Type{T}) where {T<:ParametricModel}= Yields.Rate{Float64, typeof(DEFAULT_COMPOUNDING)}

"""
    NelsonSiegel(rates::AbstractVector, maturities::AbstractVector; τ_initial=1.0)

Return the NelsonSiegel fitted parameters. The rates should be zero spot rates. If `rates` are not `Rate`s, then they will be interpreted as `Continuous` `Rate`s.

    NelsonSiegel(β₀, β₁, β₂, τ₁)

Parameters of Nelson and Siegel (1987) parametric model:

- β₀ represents a long-term interest rate
- β₁ represents a time-decay component
- β₂ represents a hump
- τ₁ controls the location of the hump

# Examples

```julia-repl
julia> β₀, β₁, β₂, τ₁ = 0.6, -1.2, -1.9, 3.0
julia> nsm = Yields.NelsonSiegel.(β₀, β₁, β₂, τ₁)

# Extend Help

## References
- https://onriskandreturn.com/2019/12/01/nelson-siegel-yield-curve-model/
- https://www.bis.org/publ/bppdf/bispap25.pdf

```
"""
struct NelsonSiegelCurve{T} <: ParametricModel
    β₀::T
    β₁::T
    β₂::T
    τ₁::T

    function NelsonSiegelCurve(β₀::T, β₁::T, β₂::T, τ₁::T) where {T<:Real}
        (τ₁ <= 0) && throw(DomainError("Wrong parameter ranges"))
        return new{T}(β₀, β₁, β₂, τ₁)
    end
end

"""
    NelsonSiegel(τ_initial)
    NelsonSiegel() # defaults to τ_initial=1.0

This parameter set is used to fit the Nelson-Siegel parametric model to given rates. `τ_initial` should be a scalar and is used as the starting τ value in the optimization. The default value for `τ_initial` is 1.0.

When fitting rates using this `YieldCurveFitParameters` object, the Nelson-Siegel model is used. If constructing curves and the rates are not `Rate`s (ie you pass a `Vector{Float64}`), then they will be interpreted as `Continuous` `Rate`s.

See for more:

- [`Zero`](@ref)
- [`Forward`](@ref)
- [`Par`](@ref)
- [`CMT`](@ref)
- [`OIS`](@ref)
"""
struct NelsonSiegel{T} <: YieldCurveFitParameters
    τ_initial::T
end
NelsonSiegel() = NelsonSiegel(1.0)

function Base.zero(ns::NelsonSiegelCurve, t)
    if iszero(t)
        # zero rate is undefined for t = 0
        t += eps()
    end
    Continuous.(ns.β₀ .+ ns.β₁ .* (1.0 .- exp.(-t ./ ns.τ₁)) ./ (t ./ ns.τ₁) .+ ns.β₂ .* ((1.0 .- exp.(-t ./ ns.τ₁)) ./ (t ./ ns.τ₁) .- exp.(-t ./ ns.τ₁)))
end
discount(ns::NelsonSiegelCurve, t) = discount.(zero.(ns,t),t)

function fit_β(ns::NelsonSiegel,func,yields,maturities,τ)
    Δₘ = vcat([maturities[1]], diff(maturities))
    param₀ = [1.0, 0.0, 0.0]
    _rate(m, p) = rate.(func.(NelsonSiegelCurve(p[1], p[2], p[3],only(τ)), m))

    return LsqFit.curve_fit(_rate, maturities, yields, Δₘ,param₀)
end

function __fit_NS(ns::NelsonSiegel,func,yields,maturities,τ)
    f(τ) = β_sum_sq_resid(ns,func,yields,maturities,τ)
    r = Optim.optimize(f, [ns.τ_initial])

    τ = only(Optim.minimizer(r))

    return τ, fit_β(ns,func,yields,maturities,τ)
end

function β_sum_sq_resid(ns,func,yields,maturities,τ)
    result = fit_β(ns,func,yields,maturities,τ)
    return sum(r^2 for r in result.resid)
end

function Zero(ns::NelsonSiegel,yields, maturities=eachindex(yields))
    yields = rate.(Continuous().(yields))
    func = zero
    τ, result = __fit_NS(ns,func,yields,maturities,ns.τ_initial)
    return NelsonSiegelCurve(result.param[1], result.param[2], result.param[3], τ)
end

function Par(ns::NelsonSiegel,yields, maturities=eachindex(yields))
    yields = rate.(Continuous().(yields))
    func = par
    τ, result = __fit_NS(ns,func,yields,maturities,ns.τ_initial)
    return NelsonSiegelCurve(result.param[1], result.param[2], result.param[3],result.param[4],  first(τ), last(τ))
end

function Forward(ns::NelsonSiegel,yields, maturities=eachindex(yields))
    yields = rate.(Continuous().(yields))
    func = forward
    τ, result = __fit_NS(ns,func,yields,maturities,ns.τ_initial)
    return NelsonSiegelCurve(result.param[1], result.param[2], result.param[3],result.param[4],  first(τ), last(τ))
end

"""
    NelsonSiegelSvensson(yields::AbstractVector, maturities::AbstractVector; τ_initial=[1.0,1.0])

Return the NelsonSiegelSvensson fitted parameters. The rates should be continuous zero spot rates. If `rates` are not `Rate`s, then they will be interpreted as `Continuous` `Rate`s.

When fitting rates using this `YieldCurveFitParameters` object, the Nelson-Siegel model is used. If constructing curves and the rates are not `Rate`s (ie you pass a `Vector{Float64}`), then they will be interpreted as `Continuous` `Rate`s.

See for more:

    - [`Zero`](@ref)
    - [`Forward`](@ref)
    - [`Par`](@ref)
    - [`CMT`](@ref)
    - [`OIS`](@ref)

    NelsonSiegelSvensson(β₀, β₁, β₂, β₃, τ₁, τ₂)

Parameters of Svensson (1994) parametric model:

- β₀ represents a long-term interest rate
- β₁ represents a time-decay component
- β₂ represents a hump
- β₃ represents a second hum
- τ₁ controls the location of the hump
- τ₁ controls the location of the second hump

# Examples

```julia-repl
julia> β₀, β₁, β₂, β₃, τ₁, τ₂ = 0.6, -1.2, -2.1, 3.0, 1.5
julia> nssm = NelsonSiegelSvensson.NelsonSiegelSvensson.(β₀, β₁, β₂, β₃, τ₁, τ₂)

## References
- https://onriskandreturn.com/2019/12/01/nelson-siegel-yield-curve-model/
- https://www.bis.org/publ/bppdf/bispap25.pdf

```
"""
struct NelsonSiegelSvenssonCurve{T} <: ParametricModel
    β₀::T
    β₁::T
    β₂::T
    β₃::T
    τ₁::T
    τ₂::T

    function NelsonSiegelSvenssonCurve(β₀::T, β₁::T, β₂::T, β₃::T, τ₁::T, τ₂::T) where {T<:Real}
        (τ₁ <= 0 || τ₂ <= 0) && throw(DomainError("Wrong parameter ranges"))
        return new{T}(β₀, β₁, β₂, β₃, τ₁, τ₂)
    end
end

"""
    NelsonSiegelSvensson(τ_initial)
    NelsonSiegelSvensson() # defaults to τ_initial=[1.0,1.0]

This parameter set is used to fit the Nelson-Siegel parametric model to given rates. `τ_initial` should be a two element vector and is used as the starting τ value in the optimization. The default value for `τ_initial` is [1.0,1.0].

See for more:

- [`Zero`](@ref)
- [`Forward`](@ref)
- [`Par`](@ref)
- [`CMT`](@ref)
- [`OIS`](@ref)
"""
struct NelsonSiegelSvensson{T} <: YieldCurveFitParameters
    τ_initial::T
end
NelsonSiegelSvensson() = NelsonSiegelSvensson([1.0,1.0])

function fit_β(ns::NelsonSiegelSvensson,func,yields,maturities,τ)
    Δₘ = vcat([maturities[1]], diff(maturities))
    param₀ = [1.0, 0.0, 0.0, 0.0]
    _rate(m, p) = rate.(func.(NelsonSiegelSvenssonCurve(p[1], p[2], p[3],p[4],first(τ),last(τ)), m))
    return LsqFit.curve_fit(_rate, maturities, yields, Δₘ,param₀)
end

function __fit_NS(ns::NelsonSiegelSvensson,func,yields,maturities,τ)
    f(τ) = β_sum_sq_resid(ns,func,yields,maturities,τ)
    r = Optim.optimize(f, ns.τ_initial)

    τ = Optim.minimizer(r)[[1,2]]

    return τ, fit_β(ns,func,yields,maturities,τ)
end

function Zero(ns::NelsonSiegelSvensson,yields, maturities=eachindex(yields))
    yields = rate.(Continuous().(yields))
    func = zero
    τ, result = __fit_NS(ns,func,yields,maturities,ns.τ_initial)
    return NelsonSiegelSvenssonCurve(result.param[1], result.param[2], result.param[3],result.param[4],  first(τ), last(τ))
end

function Par(ns::NelsonSiegelSvensson,yields, maturities=eachindex(yields))
    yields = rate.(Continuous().(yields))
    func = par
    τ, result = __fit_NS(ns,func,yields,maturities,ns.τ_initial)
    return NelsonSiegelSvenssonCurve(result.param[1], result.param[2], result.param[3],result.param[4],  first(τ), last(τ))
end

function Forward(ns::NelsonSiegelSvensson,yields, maturities=eachindex(yields))
    yields = rate.(Continuous().(yields))
    func = forward
    τ, result = __fit_NS(ns,func,yields,maturities,ns.τ_initial)
    return NelsonSiegelSvenssonCurve(result.param[1], result.param[2], result.param[3],result.param[4],  first(τ), last(τ))
end

function Base.zero(nss::NelsonSiegelSvenssonCurve, t)
    if iszero(t)
        # zero rate is undefined for t = 0
        t += eps()
    end
    Continuous.(nss.β₀ .+ nss.β₁ .* (1.0 .- exp.(-t ./ nss.τ₁)) ./ (t ./ nss.τ₁) .+ nss.β₂ .* ((1.0 .- exp.(-t ./ nss.τ₁)) ./ (t ./ nss.τ₁) .- exp.(-t ./ nss.τ₁)) .+ nss.β₃ .* ((1.0 .- exp.(-t ./ nss.τ₂)) ./ (t ./ nss.τ₂) .- exp.(-t ./ nss.τ₂)))
end
discount(nss::NelsonSiegelSvenssonCurve, t) = discount.(zero.(nss,t),t)


text:
using Libdl

if Sys.islinux()
    version = "2.11.0"
    const has_driver = !isempty(Libdl.find_library(["libcassandra"]))
    const has_yum = try success(`yum --version`) catch e false end
    const has_apt = try success(`apt-get -v`) && success(`apt-cache -v`) catch e false end

    if has_driver
        println("Cassandra CPP driver already installed.")
    elseif has_yum
        cass_url = "http://downloads.datastax.com/cpp-driver/centos/7/cassandra/v" * version * "/"
        cass_file = "cassandra-cpp-driver-" * version * "-1.el7.centos.x86_64.rpm"
        cass_source = cass_url * cass_file
        cass_target = "/tmp/cassandra-cpp-driver.rpm"
        dl = try success(`wget -O $cass_target $cass_source`) catch e false end
        !dl && error("Unable to download CPP driver.")
        inst = try success(`sudo yum install -y $cass_target`) catch e false end
        !inst && error("Unable to install CPP driver.")
    elseif has_apt
        ubuntu_version = chomp(read(pipeline(`cat /etc/os-release`, `grep -Eo "VERSION_ID=\"[0-9\.]+\""`, `grep -Eo "[^\"]+"`, `grep -E "[0-9.]+"`), String))
        if (parse(Int64, ubuntu_version) > 10)
            cass_url = "http://downloads.datastax.com/cpp-driver/ubuntu/$(ubuntu_version)/cassandra/v" * version * "/"
            cass_file = "cassandra-cpp-driver_" * version * "-1_amd64.deb"
            cass_source = cass_url * cass_file
            cass_target = "/tmp/cassandra-cpp-driver.deb"
            libuv_inst = success(`sudo apt install libuv1`)
            !libuv_inst && error("Unable to install libuv driver.")
            dl = try success(`wget -O $cass_target $cass_source`) catch e false end
            !dl && error("Unable to download CPP driver.", cass_source)
            inst = try success(`sudo dpkg -i $cass_target`) catch e false end
            !inst && error("Unable to install CPP driver.")
        else
            println("Auto installation is not available for this OS version.")
        end
    else
        error("This package requires cassandra-cpp-driver to be installed, but the build system only understands apt and yum.")
    end
elseif Sys.iswindows()
    version = "2.11.0"

    const has_driver = !isempty(Libdl.find_library(["cassandra"]))

    if has_driver
        println("Cassandra CPP driver already installed.")
    else
        using ZipFile
        cass_url = "https://downloads.datastax.com/cpp-driver/windows/cassandra/v" * version * "/"
        if typeof(1) == Int64
            cass_file = "cassandra-cpp-driver-" * version * "-win64-msvc141.zip"
        else
            cass_file = "cassandra-cpp-driver-" * version * "-win32-msvc140.zip"
        end
        cass_source = cass_url * cass_file
        cass_target = tempdir() * "/cassandra-cpp-driver.zip"
        dl = download(cass_source, cass_target)
        rd = ZipFile.Reader(dl)
        for f in rd.files
            if f.name == "bin/cassandra.dll"
                write("C:\\Temp\\cassandra.dll", read(f, String))
            end
        end
    end
else
    error("This package does not support OSX")
end