text:
#=
   DONALD+GERALD=ROBERT in ConstraintSolver.jl

   Model created by Hakan Kjellerstrand, hakank@gmail.com
   See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
const CS = ConstraintSolver
include("constraints_utils.jl")

function donald_gerald_robert()
    model = Model(optimizer_with_attributes(CS.Optimizer, "all_solutions"=>true,"logging"=>[]))
    @variable(model, 0 <= x[1:10] <= 9, Int)
    d, o, n, a, l, g, e, r, b, t = x

    @constraint(model, x in CS.AllDifferent())
    @constraint(model, d != 0)
    @constraint(model, g != 0)
    @constraint(model, r != 0)
    @constraint(model,
                              100000*d + 10000*o + 1000*n + 100*a + 10*l + d +
                              100000*g + 10000*e + 1000*r + 100*a + 10*l + d ==
                              100000*r + 10000*o + 1000*b + 100*e + 10*r + t

                              # Different approach:
                              # sum([100000*d,10000*o,1000*n,100*a,10*l + d]) +
                              # sum([100000*g,10000*e,1000*r, 100*a, 10*l, d]) ==
                              # sum([100000*r,10000*o, 1000*b, 100*e, 10*r, t])

        )

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL
        x = convert.(Integer,JuMP.value.(x))
        println(x)
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
    end
end

donald_gerald_robert()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using PySide
reload(Pkg.dir("PySide", "examples", "DThree.jl"))

## This is an attempt to integrate D3 graphics with julia via QtWebView.
## The performance is *surpringly* poor due to really slow passing of JavaScript
## from julia to QtWebView via its evaluateJavaScript method. Oh well. Using the
## webstack being developed by the hacker school is the right place for this type
## of work anyways, as it is more general and does not rely on Qt, Python, PySide
## installations

using DThree
using JSON

## Two examples

## Follow example on bar charts: http://mbostock.github.io/d3/tutorial/bar-1.html
w = Widget()

style = "
.chart div {
  font: 10px sans-serif;
  background-color: steelblue;
  text-align: right;
  padding: 3px;
  margin: 1px;
  color: white;
}
"

d3 = D3(w.w, style)

push!(w, d3.web)
raise(w)

data = [4, 8, 15, 16, 23, 42]

d3.eval_js("var data = $(json(data));")
d3.select("body").append("div").attr("class", "chart").render("chart")
## initial part of demo
d3.receiver("chart").selectAll("div").data(AsIs("data")).enter().append("div").style("width", AsIs("function(d) { return d * 10 + \"px\"; }")).text(AsIs("function(d) { return d; }")).render()
## using a scale instead
d3.scale_linear().domain([0, maximum(data)]).range(["0px", "420px"]).render("x")
d3.receiver("chart").selectAll("div").data(AsIs("data")).enter().append("div").style("width",AsIs("x")).text(AsIs("String")).render()

## A scatterplot example following http://ofps.oreilly.com/titles/9781449339739/
## not quite right, but not worth fixing up

style = "
.axis path,
.axis line {
    fill: none;
    stroke: black;
    shape-rendering: crispEdges;
}

.axis text {
    font-family: sans-serif;
    font-size: 11px;
}
"

    ## Make scatter plot
    x = linspace(0, 2*pi, 250)
    y = cos(x)
    dataset = "[" * join(["[$x,$y]" for (x,y) in zip(x,y)], ", ") * "]"

        width, height = 600, 400
        padding = 30

    d3.select("body").append("svg").attr("width", width).attr("height", height).render("svg")
    ## scales
    d3.scale_linear().domain([minimum(x),maximum(x)]*1.1).range([0, width]).render("xScale")
    d3.scale_linear().domain([minimum(y),maximum(y)]*1.1).range([height, 0]).render("yScale")
    ## axes
    d3.svg_axis().scale(I"xScale").orient("bottom").render("xAxis")
    d3.receiver("svg").append("g").attr("class", "axis").attr("transform", "translate(0,$(height - padding))").call(I"xAxis").render()
    d3.svg_axis().scale(I"yScale").orient("left").ticks(5).render("yAxis")
    d3.receiver("svg").append("g").attr("class", "axis").attr("transform", "translate($padding, 0)").call(I"yAxis").render()

        ## add points. This takes seconds for even 250 points....

    d3.receiver("svg").selectAll("circle").data(AsIs(dataset)).enter().append("circle").
       attr("cx", AsIs("function(d) {return xScale(d[0])}")).
       attr("cy", AsIs("function(d) {return yScale(d[1])}")).
       attr("r", 5).render()

ext: jl


text:
@testset "simulation" begin
    @testset "simulate" begin
        Random.seed!(60)
        # use simulate function (Munthe-Kaas integrator)
        acrobot = parse_urdf(joinpath(@__DIR__, "urdf", "Acrobot.urdf"), remove_fixed_tree_joints=false)
        x = MechanismState(acrobot)
        rand!(x)
        total_energy_before = gravitational_potential_energy(x) + kinetic_energy(x)
        times, qs, vs = simulate(x, 0.1; Δt = 1e-2)
        set_configuration!(x, qs[end])
        set_velocity!(x, vs[end])
        total_energy_after = gravitational_potential_energy(x) + kinetic_energy(x)
        @test isapprox(total_energy_after, total_energy_before, atol = 1e-3)

        total_energy_before = gravitational_potential_energy(x) + kinetic_energy(x)
        result = DynamicsResult(acrobot)

        # use RingBufferStorage
        using RigidBodyDynamics.OdeIntegrators
        passive_dynamics! = (vd::AbstractArray, sd::AbstractArray, t, state) -> begin
            dynamics!(result, state)
            copyto!(vd, result.v̇)
            copyto!(sd, result.ṡ)
            nothing
        end
        storage = RingBufferStorage{Float64}(x, 3)
        integrator = MuntheKaasIntegrator(x, passive_dynamics!, runge_kutta_4(Float64), storage)
        integrate(integrator, 0.1, 1e-2, max_realtime_rate = 2.0)
        set_configuration!(x, storage.qs[storage.last_index])
        set_velocity!(x, storage.vs[storage.last_index])
        @test isapprox(gravitational_potential_energy(x) + kinetic_energy(x), total_energy_before, atol = 1e-3)
    end

    @testset "elastic ball drop" begin
        # Drop a single rigid body with a contact point at the center of mass
        # onto the floor with a conservative contact model. Check energy
        # balances and bounce heights.
        Random.seed!(61)
        world = RigidBody{Float64}("world")
        mechanism = Mechanism(world)
        bodyframe = CartesianFrame3D()
        body = RigidBody("body", rand(SpatialInertia{Float64}, bodyframe))
        floatingjoint = Joint("floating", QuaternionFloating{Float64}())
        attach!(mechanism, world, body, floatingjoint)

        com = center_of_mass(spatial_inertia(body))
        model = SoftContactModel(hunt_crossley_hertz(; α = 0.), ViscoelasticCoulombModel(0.5, 1e3, 1e3))
        contactpoint = ContactPoint(com, model)
        add_contact_point!(body, contactpoint)

        point = Point3D(root_frame(mechanism), zero(SVector{3}))
        normal = FreeVector3D(root_frame(mechanism), 0., 0., 1.)
        halfspace = HalfSpace3D(point, normal)
        add_environment_primitive!(mechanism, halfspace)

        state = MechanismState(mechanism)
        z0 = 0.05
        zero!(state)
        tf = transform_to_root(state, body)
        tf = Transform3D(frame_after(floatingjoint), frame_before(floatingjoint), rotation(tf), SVector(1., 2., z0 - com.v[3]))
        set_configuration!(state, floatingjoint, tf)

        energy0 = gravitational_potential_energy(state)
        com_world = transform_to_root(state, body) * com

        ts, qs, vs = simulate(state, 0.5; Δt = 1e-3)

        currentsign = 0.
        switches = 0
        for (t, q, v) in zip(ts, qs, vs)
            set_configuration!(state, q)
            set_velocity!(state, v)
            twist = twist_wrt_world(state, body)
            com_world = transform_to_root(state, body) * com
            velocity = point_velocity(twist, com_world)

            z = com_world.v[3]
            penetration = max(-z, zero(z))
            n = model.normal.n
            elastic_potential_energy = model.normal.k * penetration^(n + 1) / (n + 1)
            energy = elastic_potential_energy + kinetic_energy(state) + gravitational_potential_energy(state)
            @test isapprox(energy0, energy; atol = 1e-2)

            newsign = sign(velocity.v[3])
            (newsign != currentsign) && (switches += 1)
            currentsign = newsign
        end
        @test switches > 3
    end

    @testset "inclined plane" begin
        θ = 0.5 # plane angle
        μcrit = tan(θ) # μ > μcrit should show sticking behavior, μ < μcrit should show sliding behavior

        # set up point mass + inclined plane
        world = RigidBody{Float64}("world")
        mechanism = Mechanism(world)
        floatingjoint = Joint("floating", QuaternionFloating{Float64}())
        body = RigidBody("body", SpatialInertia(CartesianFrame3D("inertia"), SMatrix{3, 3}(1.0I), zero(SVector{3}), 2.))
        attach!(mechanism, world, body, floatingjoint)
        worldframe = root_frame(mechanism)
        inclinedplane = HalfSpace3D(Point3D(worldframe, zero(SVector{3})), FreeVector3D(worldframe, sin(θ), 0., cos(θ)))
        add_environment_primitive!(mechanism, inclinedplane)
        irrelevantplane = HalfSpace3D(Point3D(worldframe, 0., 0., -100.), FreeVector3D(worldframe, 0., 0., 1.)) # #211
        add_environment_primitive!(mechanism, irrelevantplane)

        # simulate inclined plane friction experiments
        normalmodel = hunt_crossley_hertz(k = 50e3; α = 1.)
        contactlocation = Point3D(default_frame(body), 0., 0., 0.)
        for μ in (μcrit + 1e-2, μcrit - 1e-2)
            frictionmodel = ViscoelasticCoulombModel(μ, 50e3, 1e4)
            m, b = deepcopy((mechanism, body))
            add_contact_point!(b, ContactPoint(contactlocation, SoftContactModel(normalmodel, frictionmodel)))
            state = MechanismState(m)
            simulate(state, 1., Δt = 1e-3) # settle into steady state
            x1 = transform(state, contactlocation, worldframe)
            simulate(state, 0.5, Δt = 1e-3)
            x2 = transform(state, contactlocation, worldframe)
            if μ > μcrit # should stick
                @test isapprox(x1, x2, atol = 1e-4)
            else # should slip
                @test !isapprox(x1, x2, atol = 5e-2)
            end
        end
    end

    @testset "four-bar linkage" begin
        # gravitational acceleration
        g = -9.81

        # link lengths
        l_0 = 1.10
        l_1 = 0.5
        l_2 = 1.20
        l_3 = 0.75

        # link masses
        m_1 = 0.5
        m_2 = 1.0
        m_3 = 0.75

        # link center of mass offsets from the preceding joint axes
        c_1 = 0.25
        c_2 = 0.60
        c_3 = 0.375

        # moments of inertia about the center of mass of each link
        I_1 = 0.333
        I_2 = 0.537
        I_3 = 0.4

        # scalar type
        T = Float64

        # Rotation axis: negative y-axis
        axis = SVector(zero(T), -one(T), zero(T))

        world = RigidBody{T}("world")
        mechanism = Mechanism(world; gravity = SVector(0., 0., g))
        rootframe = root_frame(mechanism)

        # link1 and joint1
        joint1 = Joint("joint1", Revolute(axis))
        inertia1 = SpatialInertia(CartesianFrame3D("inertia1_centroidal"), moment=I_1*axis*axis', com=zero(SVector{3, T}), mass=m_1)
        link1 = RigidBody(inertia1)
        before_joint1_to_world = one(Transform3D, frame_before(joint1), default_frame(world))
        c1_to_joint = Transform3D(inertia1.frame, frame_after(joint1), SVector(c_1, 0, 0))
        attach!(mechanism, world, link1, joint1, joint_pose = before_joint1_to_world, successor_pose = c1_to_joint)

        # link2 and joint2
        joint2 = Joint("joint2", Revolute(axis))
        inertia2 = SpatialInertia(CartesianFrame3D("inertia2_centroidal"), moment=I_2*axis*axis', com=zero(SVector{3, T}), mass=m_2)
        link2 = RigidBody(inertia2)
        before_joint2_to_after_joint1 = Transform3D(frame_before(joint2), frame_after(joint1), SVector(l_1, 0., 0.))
        c2_to_joint = Transform3D(inertia2.frame, frame_after(joint2), SVector(c_2, 0, 0))
        attach!(mechanism, link1, link2, joint2, joint_pose = before_joint2_to_after_joint1, successor_pose = c2_to_joint)

        # link3 and joint3
        joint3 = Joint("joint3", Revolute(axis))
        inertia3 = SpatialInertia(CartesianFrame3D("inertia3_centroidal"), moment=I_3*axis*axis', com=zero(SVector{3, T}), mass=m_3)
        link3 = RigidBody(inertia3)
        before_joint3_to_world = Transform3D(frame_before(joint3), default_frame(world), SVector(l_0, 0., 0.))
        c3_to_joint = Transform3D(inertia3.frame, frame_after(joint3), SVector(c_3, 0, 0))
        attach!(mechanism, world, link3, joint3, joint_pose = before_joint3_to_world, successor_pose = c3_to_joint)

        # loop joint between link2 and link3
        joint4 = Joint("joint4", Revolute(axis))
        before_joint4_to_joint2 = Transform3D(frame_before(joint4), frame_after(joint2), SVector(l_2, 0., 0.))
        joint3_to_after_joint4 = Transform3D(frame_after(joint3), frame_after(joint4), SVector(-l_3, 0., 0.))
        attach!(mechanism, link2, link3, joint4, joint_pose = before_joint4_to_joint2, successor_pose = joint3_to_after_joint4)

        # initial state
        set_initial_state! = function (state::MechanismState)
            # found through nonlinear optimization. Slightly inaccurate.
            set_configuration!(state, joint1, 1.6707963267948966) # θ
            set_configuration!(state, joint2, -1.4591054166649482) # γ
            set_configuration!(state, joint3, 1.5397303602625536) # ϕ
            set_velocity!(state, joint1, 0.5)
            set_velocity!(state, joint2, -0.47295)
            set_velocity!(state, joint3, 0.341)
        end

        # no stabilization
        state = MechanismState(mechanism)
        set_initial_state!(state)
        zero_velocity!(state)
        energy0 = gravitational_potential_energy(state)
        ts, qs, vs = simulate(state, 1., Δt = 1e-3, stabilization_gains=nothing)
        @test kinetic_energy(state) ≉ 0 atol=1e-2
        energy1 = gravitational_potential_energy(state) + kinetic_energy(state)
        @test energy0 ≈ energy1 atol=1e-8
        @test transform(state, Point3D(Float64, frame_before(joint4)), rootframe) ≈
              transform(state, Point3D(Float64, frame_after(joint4)), rootframe) atol=1e-10 # no significant separation after a short simulation

        # with default stabilization: start with some separation
        set_initial_state!(state)
        set_configuration!(state, joint1, 1.7)
        @test transform(state, Point3D(Float64, frame_before(joint4)), rootframe) ≉
              transform(state, Point3D(Float64, frame_after(joint4)), rootframe) atol=1e-2 # significant separation initially
        simulate(state, 15., Δt = 1e-3)
        @test transform(state, Point3D(Float64, frame_before(joint4)), rootframe) ≈
              transform(state, Point3D(Float64, frame_after(joint4)), rootframe) atol=1e-5 # reduced separation after 15 seconds
        energy15 = gravitational_potential_energy(state) + kinetic_energy(state)
        simulate(state, 10, Δt = 1e-3)
        energy20 = gravitational_potential_energy(state) + kinetic_energy(state)
        @test energy20 ≈ energy15 atol=1e-5 # stabilization doesn't significantly affect energy after converging
    end
end


text:
# -------------------------------------------------------------
# Repurposed from https://github.com/tbreloff/OnlineAI.jl
# -------------------------------------------------------------

"""
Calculate the amount `Δwᵢⱼ` to adjust the parameter.  Signatures:

```
    param_change!(state::ParameterUpdaterState, updater::ParameterUpdater, ploss::ParameterLoss, gradient::Real, param::Real)
    param_change!(state::ParameterUpdaterState, updater::ParameterUpdater, gradient::Real)
```

The first allows for an additional `ParameterLoss` to be added to the final gradient.
"""
function param_change!(state::ParameterUpdaterState, updater::ParameterUpdater, ploss::ParameterLoss, gradient::Real, param::Real)
    ∇ = gradient + deriv(ploss, param)
    param_change!(state, updater, ∇)
end

# -------------------------------------------------------------

"Stochastic Gradient Descent with Momentum"
type SGDUpdater <: ParameterUpdater
    η::Float64 # learning rate
    μ::Float64 # momentum
end
SGDUpdater(; η=0.1, μ=0.5) = SGDUpdater(η, μ)

type SGDState <: ParameterUpdaterState
    lastChange::Float64
end
SGDState() = SGDState(0.0)

function param_change!(state::SGDState, updater::SGDUpdater, gradient::Real)
    state.lastChange = -updater.η * gradient + updater.μ * state.lastChange
end

# -------------------------------------------------------------

"Adaptive Gradient"
type AdagradUpdater <: ParameterUpdater
    ε::Float64  # try 0.01?
    η::Float64 # base learning rate (numerator)
end
AdagradUpdater(; ε=1e-8, η=1.0) = AdagradUpdater(ε, η)

type AdagradState <: ParameterUpdaterState
    G::Float64
end
AdagradState() = AdagradState(0.0)

function param_change!(state::AdagradState, updater::AdagradUpdater, gradient::Real)
    state.G += gradient^2
    η = updater.η / sqrt(updater.ε + state.G)
    -η * gradient
end

# -------------------------------------------------------------

"""
See: ADADELTA: An Adaptive Learning Rate Method (Zeiler 2012)

Relatively parameter-free... can probably avoid changing ε and ρ
"""
type AdadeltaUpdater <: ParameterUpdater
    ε::Float64  # try 0.01?
    η::Float64
    ρ::Float64  # try 0.97?
end
AdadeltaUpdater(; ε=1e-8, η=0.1, ρ=0.95) = AdadeltaUpdater(ε, η, ρ)

type AdadeltaState <: ParameterUpdaterState
    dMean::Float64
    GMean::Float64
end
AdadeltaState() = AdadeltaState(0.0, 0.0)

function param_change!(state::AdadeltaState, updater::AdadeltaUpdater, gradient::Real)
    ε, ρ = updater.ε, updater.ρ

    # average g²
    state.GMean = ρ * state.GMean + (1.0 - ρ) * gradient^2

    # compute learning rate from previous average dw² and current average g²
    η = updater.η * sqrt(state.dMean + ε) / sqrt(state.GMean + ε)

    # compute change and update average dw²
    dij = -η * gradient
    state.dMean = ρ * state.dMean + (1.0 - ρ) * dij^2
    dij
end

# -------------------------------------------------------------

"""
see: ADAM: A method for Stochastic Optimization (Kingma and Ba 2015)

Tracks an exponential moving average of the first and second moments of the gradient,
adjusting for zero-bias.  The defaults are those suggested in the paper.

TODO: AdaMax is similar, using the p-norm as p -> ∞
"""
type AdamUpdater <: ParameterUpdater
    ε::Float64  # small number so we don't divide by 0
    η::Float64  # learning rate... (this is α in the paper) maybe use around 1e-3?
    ρ1::Float64 # decay for first moment (β₁ in the paper)
    ρ2::Float64 # decay for second moment (β₂ in the paper)
end
AdamUpdater(; ε=1e-8, η=1e-3, ρ1=0.9, ρ2=0.999) = AdamUpdater(ε, η, ρ1, ρ2)

immutable AdamState{T<:Number, M<:AbstractArray{T}} <: ParameterUpdaterState
    m::M      # average first moment
    v::M      # average second moment
    ρ1t::M    # β₁ᵗ from the paper... t-th power of β₁
    ρ2t::M    # β₂ᵗ from the paper... t-th power of β₂
end
AdamState() = AdamState(0.0, 0.0, 1.0, 1.0)

function param_change!(state::AdamState{T,M}, updater::AdamUpdater, gradient::M)
    ρ1, ρ2 = updater.ρ1, updater.ρ2
    state.m = ρ1 * state.m + (1.0 - ρ1) * gradient
    state.v = ρ2 * state.v + (1.0 - ρ2) * gradient^2
    state.ρ1t *= updater.ρ1
    state.ρ2t *= updater.ρ2
    ηt = updater.η * (sqrt(1.0 - state.ρ2t) / (1.0 - state.ρ1t))
    -ηt * state.m / (sqrt(state.v + updater.ε))
end

function param_change!(state::AdamState{T,M}, updater::AdamUpdater, gradient::T, i::Int...)

# -------------------------------------------------------------

"""
see: ADAM: A method for Stochastic Optimization (Kingma and Ba 2015)

AdaMax is similar to Adam, using the p-norm as p -> ∞
"""
type AdaMaxUpdater <: ParameterUpdater
    # ε::Float64  # small number so we don't divide by 0
    η::Float64  # learning rate... (this is α in the paper) maybe use around 1e-3?
    ρ1::Float64 # decay for first moment (β₁ in the paper)
    ρ2::Float64 # decay for second moment (β₂ in the paper)
end
AdaMaxUpdater(; η=1e-3, ρ1=0.9, ρ2=0.99) = AdaMaxUpdater(η, ρ1, ρ2)

type AdaMaxState <: ParameterUpdaterState
    m::Float64      # average first moment
    u::Float64      # average second moment
    ρ1t::Float64    # β₁ᵗ from the paper... t-th power of β₁
    # ρ2t::Float64  # β₂ᵗ from the paper... t-th power of β₂
end
AdaMaxState() = AdaMaxState(0.0, 0.0, 1.0)

function param_change!(state::AdaMaxState, updater::AdaMaxUpdater, gradient::Real)
    ρ1 = updater.ρ1
    mij = ρ1 * state.m + (1.0 - ρ1) * gradient
    state.m = mij
    uij = max(updater.ρ2 * state.u, abs(gradient))
    state.u = uij
    state.ρ1t *= ρ1
    -updater.η * mij / ((uij + 1e-10) * (1.0 - state.ρ1t[1]))
end

# -------------------------------------------------------------
# Constructors for updater states
# -------------------------------------------------------------

get_state_type(::Type{SGDUpdater})      = SGDState
get_state_type(::Type{AdagradUpdater})  = AdagradState
get_state_type(::Type{AdadeltaUpdater}) = AdadeltaState
get_state_type(::Type{AdamUpdater})     = AdamState
get_state_type(::Type{AdaMaxUpdater})   = AdaMaxState

# single state
ParameterUpdaterState{U<:ParameterUpdater}(::Type{U}) = get_state_type(U)()
ParameterUpdaterState(updater::ParameterUpdater) = ParameterUpdaterState(typeof(updater))

# arrays of states
function ParameterUpdaterState{U<:ParameterUpdater}(::Type{U}, dims::Integer...)
    T = get_state_type(U)
    states = Array(T, dims...)
    for i in eachindex(states)
        states[i] = T()
    end
    states
end
ParameterUpdaterState(updater::ParameterUpdater, dims::Integer...) = ParameterUpdaterState(typeof(updater), dims...)

# -------------------------------------------------------------

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
##############
# RTN | LVLH #
##############

export rRTNtoECI
"""
Compute the radial, along-track, cross-track (RTN) to Earth-centered inertial
rotation matrix. If applied to a position vector in the RTN frame, it will
transform that vector to into the equivalent position vector in the ECI frame.

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite

Returns:
- `R_rtn2eci::SMatrix{3,3}`: Rotation matrix transforming _from_ the RTN frame _to_ the ECI frame
"""
function rRTNtoECI(x::AbstractVector{<:Real})
    r = x[idx1t3]
    v = x[idx4t6]
    n = cross(r, v)

    R = normalize(r)
    N = normalize(n)
    T = cross(N, R)

    R_rtn2eci = hcat(R, T, N)

    return R_rtn2eci
end

export rECItoRTN
"""
Compute the Earth-centered inertial to radial, along-track, cross-track (RTN)
rotation matrix. If applied to a position vector in the ECI frame, it will
transform that vector into the equivalent position vector in the RTN frame.

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite

Returns:
- `R_eci2rtn::SMatrix{3,3}`: Rotation matrix transforming _from_ the ECI frame _to_ the RTN frame
"""
function rECItoRTN(x::AbstractVector{<:Real})
    r = x[idx1t3]
    v = x[idx4t6]
    n = cross(r, v)

    R = normalize(r)
    N = normalize(n)
    T = cross(N, R)

    R_eci2rtn = hcat(R, T, N)'

    return R_eci2rtn
end

export sECItoRTN
"""
Compute the radial, along-track, cross-track (RTN) coordinates of a target satellite in the primary satellites RTN frame.

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xt::AbstractVector{<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xrtn::AbstractVector{<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sECItoRTN(x::AbstractVector{<:Real}, xt::AbstractVector{<:Real})
    if length(xt) >= 6
        return sECItoRTN(x, xt[idx1t6])
    else
        return sECItoRTN(x, xt[idx1t3])
    end
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xt::SVector{3,<:Real}`: Inertial state (position) of the target satellite

Returns:
- `xrtn::SVector{3,<:Real}`: Position of the target relative of the observing satellite in the RTN.
"""
function sECItoRTN(x::AbstractVector{<:Real}, xt::SVector{3,<:Real})
    # Create RTN rotation matrix
    R_eci2rtn = rECItoRTN(x)

    # Transform Position
    r     = x[idx1t3]
    rho   = xt[idx1t3] - r
    r_rtn = R_eci2rtn*rho

    return r_rtn
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xt::SVector{6,<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xrtn::SVector{6,<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sECItoRTN(x::AbstractVector{<:Real}, xt::SVector{6,<:Real})
    # Create RTN rotation matrix
    R_eci2rtn = rECItoRTN(x)

    # Transform Position
    r     = x[idx1t3]
    rho   = xt[idx1t3] - r
    r_rtn = R_eci2rtn*rho

    v       = x[idx4t6]
    f_dot   = norm(cross(r, v)) / norm(r)^2
    omega   = SVector(0, 0, f_dot)
    rho_dot = xt[idx4t6] - v
    v_rtn   = R_eci2rtn*rho_dot - cross(omega, r_rtn)

    return vcat(r_rtn, v_rtn)
end

export sRTNtoECI
"""
Compute the Earth-center

The RTN frame is also commonly refered to as the local-vertical, local-horizontal (LVLH) frame.

Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xrtn::AbstractVector{<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xt::AbstractVector{<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sRTNtoECI(x::AbstractVector{<:Real}, xrtn::AbstractVector{<:Real})
    if length(xrtn) >= 6
        return sRTNtoECI(x, xrtn[idx1t6])
    else
        return sRTNtoECI(x, xrtn[idx1t3])
    end
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xrtn::SVector{3,<:Real}`: Inertial state (position) of the target satellite

Returns:
- `xt::SVector{3,<:Real}`: Position of the target relative of the observing satellite in the RTN.
"""
function sRTNtoECI(x::AbstractVector{<:Real}, xrtn::SVector{3,<:Real})
    # Create RTN rotation matrix
    R_rtn2eci = rRTNtoECI(x)

    # Transform position
    r   = x[idx1t3]
    r_t = R_rtn2eci*xrtn + r

    return r_t
end

"""
Arguments:
- `x::AbstractVector{<:Real}`: Inertial state (position and velocity) of primary (observing) satellite
- `xrtn::SVector{6,<:Real}`: Inertial state (position and velocity) of the target satellite

Returns:
- `xt::SVector{6,<:Real}`: Position and velocity of the target relative of the observing satellite in the RTN.
"""
function sRTNtoECI(x::AbstractVector{<:Real}, xrtn::SVector{6,<:Real})
    # Create RTN rotation matrix
    R_rtn2eci = rRTNtoECI(x)

    # Transform position
    r     = x[idx1t3]
    r_rtn = xrtn[idx1t3]
    r_t   = R_rtn2eci*r_rtn + r

    # Transform velocity
    v     = x[idx4t6]
    v_rtn = xrtn[idx4t6]
    f_dot = norm(cross(r, v)) / norm(r)^2
    omega = SVector(0, 0, f_dot)
    v_t   = R_rtn2eci*(v_rtn + cross(omega, r_rtn)) + v

    return vcat(r_t, v_t)
end

#######################################
# IAU 2010 | Inertial <-> Earth-Fixed #
#######################################

export bias_precession_nutation
"""
Computes the Bias-Precession-Nutation matrix transforming the GCRS to the
CIRS intermediate reference frame. This transformation corrects for the
bias, precession, and nutation of Celestial Intermediate Origin (CIO) with
respect to inertial space.

Arguments:
- `epc::Epoch`: Epoch of transformation

Returns:
- `rc2i::Matrix{<:Real}`: 3x3 Rotation matrix transforming GCRS -> CIRS
"""
function bias_precession_nutation(epc::Epoch)
    # Constants of IAU 2006A transofrmation
    DMAS2R =  4.848136811095359935899141e-6 / 1.0e3
    dx06   =  0.0001750*DMAS2R
    dy06   = -0.0002259*DMAS2R

    # Compute X, Y, s terms using low-precision series terms
    x, y, s = iauXys00b(MJD_ZERO, mjd(epc, tsys="TT"))

    # Apply IAU2006 Offsets
    x += dx06
    y += dy06

    # Compute transformation and return
    rc2i = iauC2ixys(x, y, s)

    return rc2i
end

export earth_rotation
"""
Computes the Earth rotation matrix transforming the CIRS to the TIRS
intermediate reference frame. This transformation corrects for the Earth
rotation.

Arguments:
- `epc::Epoch`: Epoch of transformation

Returns:
- `r::Matrix{<:Real}`: 3x3 Rotation matrix transforming CIRS -> TIRS
"""
function earth_rotation(epc::Epoch)
    # Compute Earth rotation angle
    era = iauEra00(MJD_ZERO, mjd(epc, tsys="UT1"))

    # Rotate Matrix and return
    r = iauRz(era, Matrix{Float64}(I, 3, 3))

    return r
end

export polar_motion
"""
Computes the Earth rotation matrix transforming the TIRS to the ITRF reference
frame.

# Arguments
- `epc::Epoch`: Epoch of transformation

# Returns
- `rpm::Matrix{<:Real}`: 3x3 Rotation matrix transforming TIRS -> ITRF
"""
function polar_motion(epc::Epoch)
    xp, yp = POLE_LOCATOR(mjd(epc, tsys="UTC"))

    # Compute transformation and return
    rpm = iauPom00(xp, yp, iauSp00(MJD_ZERO, mjd(epc, tsys="TT")))

    return rpm
end

export rECItoECEF
"""
Computes the combined rotation matrix from the inertial to the Earth-fixed
reference frame. Applies corrections for bias, precession, nutation,
Earth-rotation, and polar motion.

The transformation is accomplished using the IAU 2006/2000A, CIO-based
theory using classical angles. The method as described in section 5.5 of
the SOFA C transformation cookbook.

# Arguments
- `epc::Epoch`: Epoch of transformation

# Returns
- `r::Matrix{<:Real}`: 3x3 Rotation matrix transforming GCRF -> ITRF
"""
function rECItoECEF(epc::Epoch)
    # Compute intermediate transformations
    rc2i = bias_precession_nutation(epc)
    r    = earth_rotation(epc)
    rpm  = polar_motion(epc)

    return rpm * r * rc2i
end

export rECEFtoECI
"""
Computes the combined rotation matrix from the Earth-fixed to the inertial
reference frame. Applies corrections for bias, precession, nutation,
Earth-rotation, and polar motion.

The transformation is accomplished using the IAU 2006/2000A, CIO-based
theory using classical angles. The method as described in section 5.5 of
the SOFA C transformation cookbook.

# Arguments
- `epc::Epoch`: Epoch of transformation

# Returns
- `r::Matrix{<:Real}`: 3x3 Rotation matrix transforming ITRF -> GCRF
"""
function rECEFtoECI(epc::Epoch)
    # Compute intermediate transformations
    rc2i = bias_precession_nutation(epc)
    r    = earth_rotation(epc)
    rpm  = polar_motion(epc)

    return rc2i' * r' * rpm'
end

export sECItoECEF
"""
Transforms an Earth inertial state into an Earth fixed state

The transformation is accomplished using the IAU 2006/2000A, CIO-based
theory using classical angles. The method as described in section 5.5 of
the SOFA C transformation cookbook.

# Arguments
- `epc::Epoch`: Epoch of transformation
- `x::AbstractVector{<:Real}`: Inertial state (position, velocity) [m; m/s]

# Returns
- `x_ecef::Vector{<:Real}`: Earth-fixed state (position, velocity)
"""
function sECItoECEF(epc::Epoch, x::AbstractVector{<:Real})
    dim_x  = length(x)
    x_ecef = zeros(Float64, dim_x)

    # Extract State Components
    r_eci = x[1:3]

    if dim_x >= 6
        v_eci = x[4:6]
    end

    if dim_x == 9
        a_eci = x[7:9]
    end

    # Compute Sequential Transformation Matrices
    rc2i = bias_precession_nutation(epc)
    r    = earth_rotation(epc)
    pm   = polar_motion(epc)

    # Create Earth's Angular Rotation Vector
    omega_vec = [0, 0, OMEGA_EARTH] # Neglect LOD effect

    # Calculate ECEF State
    x_ecef[1:3] = pm *  r * rc2i * r_eci

    if dim_x == 6
        x_ecef[4:6] = pm * (r * rc2i * v_eci - cross(omega_vec, r * rc2i * r_eci))
    end

    if dim_x == 9
        x_ecef[7:9] = pm * (r * rc2i * a_eci - cross(omega_vec, cross(omega_vec, r * rc2i * r_eci))
                                         - 2 * cross(omega_vec, r * rc2i * v_eci))
    end

    return x_ecef
end

export sECEFtoECI
"""
Transforms an Earth fixed state into an Inertial state

The transformation is accomplished using the IAU 2006/2000A, CIO-based
theory using classical angles. The method as described in section 5.5 of
the SOFA C transformation cookbook.

# Arguments
- `epc::Epoch`: Epoch of transformation
- `x::AbstractVector{<:Real}`: Earth-fixed state (position, velocity) [m; m/s]

# Returns
- `x_eci::Vector{<:Real}`: Inertial state (position, velocity)
"""
function sECEFtoECI(epc::Epoch, x::AbstractVector{<:Real})
    # Set state variable size
    dim_x = length(x)
    x_eci = zeros(Float64, dim_x)

    # Extract State Components
    r_ecef = x[1:3]

    if dim_x >= 6
        v_ecef = x[4:6]
    end

    if dim_x == 9
        a_ecef = x[7:9]
    end

    # Compute Sequential Transformation Matrices
    bpn = bias_precession_nutation(epc)
    rot = earth_rotation(epc)
    pm  = polar_motion(epc)

    # Create Earth's Angular Rotation Vector
    omega_vec = [0, 0, OMEGA_EARTH] # Neglect LOD effect

    # Calculate ECEF State
    x_eci[1:3] = (pm * rot * bpn)' * r_ecef

    if dim_x >= 6
        x_eci[4:6] = (rot * bpn)' * ((pm)' * v_ecef + cross(omega_vec, (pm)' * r_ecef))
    end

    if dim_x >= 9
        x_eci[7:9] = (rot * bpn)' * ((pm)' * a_ecef + cross(omega_vec, cross(omega_vec, (pm)' * x_eci[1:3]))
                                 + 2 * cross(omega_vec, (pm)' * x_eci[4:6]))
    end

    return x_eci
end


text:
function solvesplits!(iters::Vector{Int64},Q::Vector{Float64},A::Vector{Float64},
    children::Vector{Int64},Qold::Vector{Float64},Aold::Vector{Float64},W::Vector{Float64},
    beta::Vector{Float64},A0::Vector{Float64},c0::Vector{Float64},rho::Float64,
    f::Function,J::Function,maxiter::Int16,maxval::Float64,epsJ::Float64,epsN::Float64)
    LegHemodynamics.newton!(iters,Q,A,children,Qold,Aold,W,beta,A0,c0,rho,f,J,
        maxiter,maxval,epsJ,epsN);
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# Autogenerated wrapper script for Zenith_jll for aarch64-apple-darwin
export zenith

JLLWrappers.@generate_wrapper_header("Zenith")
JLLWrappers.@declare_executable_product(zenith)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        zenith,
        "bin/zenith",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    PathGraph{T<:Integer}

A path graph, with each node linked to the previous and next one.
"""
struct PathGraph{T<:Integer} <: LG.AbstractGraph{T}
    nv::Int
    function PathGraph{T}(nv::Integer) where {T<:Integer}
        _nv = nv >= 0 ? Int(nv) : 0
        new{T}(_nv)
    end
end

PathGraph(nv::Integer) = PathGraph{Int}(nv)

LG.edgetype(::PathGraph) = LG.Edge{Int}
LG.is_directed(::Type{<:PathGraph}) = false
LG.nv(g::PathGraph) = g.nv
LG.ne(g::PathGraph) = LG.nv(g) - 1
LG.vertices(g::PathGraph) = 1:LG.nv(g)

LightGraphs.edges(g::PathGraph) = [LG.Edge(i, i+1) for i in 1:LG.nv(g)-1]

LightGraphs.has_vertex(g::PathGraph, v) = 1 <= v <= LG.nv(g)

function LightGraphs.outneighbors(g::PathGraph, v)
    LG.has_vertex(g, v) || return Int[]
    nv(g) > 1 || return Int[]
    if v == 1
        return [2]
    end
    if v == nv(g)
        return [nv(g)-1]
    end
    return [v-1, v+1]
end

LightGraphs.inneighbors(g::PathGraph, v) = outneighbors(g, v)

function LightGraphs.has_edge(g::PathGraph, v1, v2)
    if !has_vertex(g, v1) || !has_vertex(g, v2)
        return false
    end
    return abs(v1-v2) == 1
end


text:
function initialize!(integrator,cache::KuttaPRK2p5ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::KuttaPRK2p5ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α21,α31,α32,α41,α42,α43,α5_6 = cache
  @unpack β1,β3,β5,β6,c2,c3,c4,c5_6 = cache

  k1 = f(uprev, p, t)
  k2 = f(uprev + dt*α21*k1, p, t + c2*dt)
  k3 = f(uprev + dt*(α31*k1 + α32*k2), p, t + c3*dt)
  k4 = f(uprev + dt*(α41*k1 + α42*k2 + α43*k3), p, t + c4*dt)

  k5_6 = Array{typeof(k1)}(undef, 2)

  if integrator.alg.threading == false
    k5_6[1] = f(uprev + dt*(α5_6[1,1]*k1 + α5_6[1,2]*k2 + α5_6[1,3]*k3 + α5_6[1,4]*k4), p, t + c5_6[1]*dt)
    k5_6[2] = f(uprev + dt*(α5_6[2,1]*k1 + α5_6[2,2]*k2 + α5_6[2,3]*k3 + α5_6[2,4]*k4), p, t + c5_6[2]*dt)
  else
    let
      Threads.@threads for i in [1,2]
        k5_6[i] = f(uprev + dt*(α5_6[i,1]*k1 + α5_6[i,2]*k2 + α5_6[i,3]*k3 + α5_6[i,4]*k4), p, t + c5_6[i]*dt)
      end
    end
  end

  u = uprev + dt*(β1*k1 + β3*k3 + β5*k5_6[1] + β6*k5_6[2])
  k = f(u, p, t+dt)

  integrator.fsallast = k # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.k[1] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator,cache::KuttaPRK2p5Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::KuttaPRK2p5Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,k1,k2,k3,k4,k5_6,fsalfirst,tmp = cache
  @unpack α21,α31,α32,α41,α42,α43,α5_6 = cache.tab
  @unpack β1,β3,β5,β6,c2,c3,c4,c5_6 = cache.tab

  f( k1, uprev, p, t)

  @.. u = uprev + dt*α21*k1
  f( k2, u, p, t + c2*dt)

  @.. u = uprev + dt*(α31*k1 + α32*k2)
  f( k3, u, p, t + c3*dt)

  @.. u = uprev + dt*(α41*k1 + α42*k2 + α43*k3)
  f( k4, u, p, t + c4*dt)

  if integrator.alg.threading == false
    @.. u = uprev + dt*(α5_6[1,1]*k1 + α5_6[1,2]*k2 + α5_6[1,3]*k3 + α5_6[1,4]*k4)
    f( k5_6[1], u, p, t + c5_6[1]*dt)

    @.. u = uprev + dt*(α5_6[2,1]*k1 + α5_6[2,2]*k2 + α5_6[2,3]*k3 + α5_6[2,4]*k4)
    f( k5_6[2], u, p, t + c5_6[2]*dt)
  else
    tmps = (u, tmp)
    let
      Threads.@threads for i in [1,2]
        @.. tmps[i] = uprev + dt*(α5_6[i,1]*k1 + α5_6[i,2]*k2 + α5_6[i,3]*k3 + α5_6[i,4]*k4)
        f( k5_6[i], tmps[i], p, t + c5_6[i]*dt)
      end
    end
  end

  @.. u = uprev + dt*(β1*k1 + β3*k3 + β5*k5_6[1] + β6*k5_6[2])
  f( k,  u, p, t+dt)
end


text:
function puzzle17(input::Int=348)

    input=348
    arr = [0]
    for iter in 1:2017
        arr = circshift(arr, -input)
        insert!(arr, 2, iter)
        arr = circshift(arr, -1)
    end
    partone = arr[2]

    cur = 0
    val = 0
    for iter in 1:50_000_000
        cur = mod1(cur + input, iter)
        if cur == 1
            val = iter
        end
        cur += 1
    end
    parttwo = val

    return [partone, parttwo]
end


text:
#' ---
#' title: Benders Decomposition for Norm-1 Regression with ParameterJuMP.jl
#' author: Joaquim Dias Garcia
#' date: March 9th 2019
#' options:
#'     github:
#'         out_path : markdown
#'     md2html:
#'         out_path : html
#'     md2pdf:
#'         out_path : pdf
#' ---

#' # Introduction

#' This notebook is part of the talk on ParameterJuMP.jl in the
#' third annual JuMP-dev workshop, held in Santiago, Chile, 2019

#' The main purpose of this notebook is to show an application of
#' [ParameterJuMP.jl](https://github.com/JuliaStochOpt/ParameterJuMP.jl).
#' ParameterJuMP is well suited for Benders like decompositions
#' therefore we shall try to demonstrate the usage of the library
#' in one of the simplest problems that fits in the Benders
#' decomposition framework, Norm-1 regression, which is a particular
#' case of quantile regression.
#'
#' Note that this is NOT the standard technique to solve Norm-1
#' regressions. Taylor made methods are available
#' [here](https://cran.r-project.org/web/packages/quantreg/quantreg.pdf)
#' for instance.
#'

#' This notebook will require the following libraries:

#' ParameterJuMP itself
using ParameterJuMP

#' JuMP: the julia mathematical programming modeling tool
using JuMP

#' GLPK: A linear programming solver (other solvers could be used -
#' such as Clp, Xpress, Gurobi, CPLEX and so on)

#+ echo = false
# using Xpress
# const OPTIMIZER = Xpress.Optimizer;
#+

using GLPK
const OPTIMIZER = GLPK.Optimizer;

#' TimerOutputs: a time measuring library to demonstrate the advantage of using ParameterJuMP
using TimerOutputs

#' The following two julia default libraries
using LinearAlgebra # just use the dot function
using Random # to use random number generators

#' Plots library
using Plots
gr(); # plotting backend

#' # Norm-1 regression

#' We will apply Norm-1 regression to the
#' [Linear Regression](https://en.wikipedia.org/wiki/Linear_regression) problem.
#' Linear regression is a statistical tool to obtain the relation
#' between one **dependent variable** and other **explanatory variables**.
#' In other words, given a set of $n$ explanatory variables $X = \{ X_1, \dots, X_n \}$
#' we would like to obtain the best possible estimate for $Y$.
#' In order to accomplish such a task we make the hypothesis that $Y$
#' is approximately linear function of $X$:
#'
#' $Y = \sum_{j =1}^n \beta_j X_j + \varepsilon$
#'
#' where $\varepsilon$ is some random error.
#'
#' The estimation of the $\beta$ values relies on observations of the variables:
#' $\{y^i, x_1^i, \dots, x_n^i\}_i$
#'
#' In this notebook we will solve a problem where the explanatory variables
#' are sinusoids of differents frequencies

#' First, we define the number of explanatory variables and observations

#+ echo = false

# GLPK small
# const N_Candidates = 100
# const N_Observations = 600
# const N_Nodes = 100

# GLPK large
# const Candidates = 600
# const Observations = 5000
# const Nodes = 1000

# Xpress
# const N_Candidates = 700
# const N_Observations = 5000
# const N_Nodes = 1000

#+

const N_Candidates = 100
const N_Observations = 600
const N_Nodes = 100

const Observations = 1:N_Observations
const Candidates = 1:N_Candidates
const Nodes = 1:N_Nodes
;

#' Initialize a random number generator to keep results deterministic

rng = Random.MersenneTwister(123);

#' Building regressors (explanatory) sinusoids

const X = zeros(N_Candidates, N_Observations)
const time = [obs / N_Observations * 1 for obs in Observations]
for obs in Observations, cand in Candidates
    t = time[obs]
    f = cand
    X[cand, obs] = sin(2 * pi * f * t)
end

#' Define coefficients

β = zeros(N_Candidates)
for i ∈ Candidates
    if rand(rng) <= (1-i/N_Candidates)^2 && i<=100
        β[i] = 4*rand(rng)/i
    end
end
println("First coefs: $(β[1:min(10, N_Candidates)])")

#' Create noisy observations
const y = X' * β .+ 0.1*randn(rng, N_Observations)

plt = plot(time, y,
    xlabel = "Time (s)", ylabel = "Amplitude")
plot!(plt, time, X'[:,1])
plot!(plt, time, X'[:,3])
plot!(plt, time, X'[:,9])

#' The classic tool to estimate linear regression models is the
#' [Least Squares](https://en.wikipedia.org/wiki/Least_squares) method.
#'
#' The least squares method relies on solving the optimization problem:
#'
#' $\max \Bigg\{ \sum_{i \in Observations} \Big( y_i - \sum_{j \in Candidates} \beta_j x_{i,j} \Big) ^2 \Bigg\}$
#'
#' In Norm-1 regression, the quadratic functions are replaced by absolute values:
#'
#' $\max\Bigg\{ \sum_{i \in Observations} \Big| y_i - \sum_{j \in Candidates} \beta_j x_{i,j} \Big| \Bigg\}$
#'
#' This optimization problem can be recast as a [Linear Programming](https://en.wikipedia.org/wiki/Linear_programming) Problem:
#'
#' $$
#'\begin{align}
#'    & \min_{\varepsilon^{up}, \varepsilon^{dw}, \beta}  &&  \sum_{i \in Observations} {\varepsilon^{up}}_i + {\varepsilon^{dw}}_i  && \notag \\
#'    & \text{subject to}     &&  {\varepsilon^{up}}_i \geq + y_i - \sum_{j \in Candidates} \beta_j x_{i,j} && \forall i \in Observations \notag \\
#'    &                       &&  {\varepsilon^{dw}}_i \geq - y_i + \sum_{j \in Candidates} \beta_j x_{i,j} && \forall i \in Observations \notag \\
#'    &                       &&  {\varepsilon^{up}}_i, {\varepsilon^{dw}}_i \geq 0 && \forall i \in Observations \notag \\
#'\end{align}
#' $$
#'
#' Where $Observations$ is the set of all observations.
#'

#' This linear programming problem can be described in julia with JuMP
function full_model_regression()
    time_build = @elapsed begin # measure time to create a model

        # initialize a optimization model
        full_model = Model(with_optimizer(OPTIMIZER))

        # create optimization variables of the problem
        @variables(full_model, begin
            ɛ_up[Observations] >= 0
            ɛ_dw[Observations] >= 0
            β[1:N_Candidates]
            # 0 <= β[Candidates] <= 8
        end)

        # define constraints of the model
        @constraints(full_model, begin
            ɛ_up_ctr[i in Observations],
                ɛ_up[i] >= + sum(X[j,i] * β[j] for j ∈ Candidates) - y[i]
            ɛ_dw_ctr[i in Observations],
                ɛ_dw[i] >= - sum(X[j,i] * β[j] for j ∈ Candidates) + y[i]
        end)

        # construct the objective function to be minimized
        @objective(full_model, Min, sum(ɛ_up[i] + ɛ_dw[i] for i ∈ Observations))
    end

    # solve the problem
    time_solve = @elapsed optimize!(full_model)

    println("First coefficients in solution: $(value.(β)[1:min(10, N_Candidates)])")
    println("Objective value: $(objective_value(full_model))")
    println("Time in solve: $time_solve")
    println("Time in build: $time_build")

    return nothing
end

#' Now we execute the function that builds the model and solves it

N_Observations*N_Candidates < 10_000_000 && full_model_regression()

#' # Benders decomposition

#' Benders decomposition is used to solve large optimization problems
#' with some special characteristics.
#' LP's can be solved with classical linear optimization methods
#' such as the Simplex method or Interior point methods provided by
#' solvers like GLPK.
#' However, these methods do not scale linearly with the problem size.
#' In the Benders decomposition framework we break the problem in two pieces:
#' A master and a slave problem.
#'
#' Of course some variables will belong to both problems, this is where the
#' cleverness of Benders kicks in:
#' The master problem is solved and passes the shared variables to the slave.
#' The slave problem is solved with the shared variables FIXED to the values
#' given by the master problem. The solution of the slave problem can be used
#' to generate a constraint to the master problem to describe the linear
#' approximation of the cost function of the shared variables.
#' In many cases, like stochastic programming, the slave problems have a
#' interesting structure and might be broken in smaller problem to be solved
#' in parallel.

#' We will descibe the decomposition similarly to what is done in:
#' Introduction to Linear Optimization, Bertsimas & Tsitsiklis (Chapter 6.5):
#' Where the problem in question has the form
#'
#' $$
#'\begin{align}
#'    & \min_{x, y_k}     &&  c^T x && + f_1^T y_1 && + \dots && + f_n^T y_n  &&  \notag \\
#'    & \text{subject to} &&  Ax    &&             &&         &&              && = b \notag \\
#'    &                   &&  B_1 x && + D_1 y_1   &&         &&              && = d_1 \notag \\
#'    &                   &&  \dots &&             &&  \dots  &&              &&       \notag \\
#'    &                   &&  B_n x &&             &&         && + D_n y_n    && = d_n \notag \\
#'    &                   &&   x,   &&     y_1,    &&         &&       y_n    && \geq 0 \notag \\
#'\end{align}
#' $$
#'

#' ### Slave

#' Given a solution for the $x$ variables we can define the slave problem as
#'
#'
#' $$
#'\begin{align}
#' z_k(x) \ = \ & \min_{y_k}        &&  f_k^T y_k &&                  \notag \\
#'              & \text{subject to} &&  D_k y_k   &&  = d_k - B_k x \notag \\
#'              &                   &&  y_k       && \geq 0          \notag \\
#'\end{align}
#' $$
#'
#'
#' The $z_k(x)$ function represents the cost of the subproblem given a
#' solution for $x$. This function is a convex function because $x$
#' affects only the right hand side of the problem (this is a standard
#' results in LP theory).
#'
#' For the special case of the Norm-1 reggression the problem is written as:
#'
#' $$
#'\begin{align}
#' z_k(\beta) \ = \ & \min_{\varepsilon^{up}, \varepsilon^{dw}}         &&  \sum_{i \in ObsSet(k)} {\varepsilon^{up}}_i + {\varepsilon^{dw}}_i               && \notag \\
#'                  & \text{subject to}     &&  {\varepsilon^{up}}_i \geq + y_i - \sum_{j \in Candidates} \beta_j x_{i,j} && \forall i \in ObsSet(k) \notag \\
#'                  &                       &&  {\varepsilon^{dw}}_i \geq - y_i + \sum_{j \in Candidates} \beta_j x_{i,j} && \forall i \in ObsSet(k) \notag \\
#'                  &                       &&  {\varepsilon^{up}}_i, {\varepsilon^{dw}}_i \geq 0                             && \forall i \in ObsSet(k) \notag \\
#'\end{align}
#' $$
#'
#' The collection $ObsSet(k)$ is a sub-set of the N_Observations.
#' Any partition of the N_Observations collection is valid.
#' In this notebook we will partition with the function:

function ObsSet(K)
    obs_per_block = div(N_Observations, N_Nodes)
    return (1 + (K - 1) * obs_per_block):(K * obs_per_block)
end

#'
#' Which can be written in JuMP as follows.
#'
#' At this point we make a small detour to highlight the ParameterJuMP
#' application. Every time you a find a IF block with the flag `PARAM`
#' it means that we have two different implementations of the method:
#' one relying on ParameterJuMP and the other using pure JuMP.
#'
#'

function slave_model(PARAM, K)

    # initialize the JuMP model
    slave = if PARAM
        # special constructor exported by ParameterJuMP
        # to add the functionality to the model
        ModelWithParams(with_optimizer(OPTIMIZER))
    else
        # regular JuMP constructor
        Model(with_optimizer(OPTIMIZER))
    end

    # Define local optimization variables for norm-1 error
    @variables(slave, begin
        ɛ_up[ObsSet(K)] >= 0
        ɛ_dw[ObsSet(K)] >= 0
    end)

    # create the regression coefficient representation
    if PARAM
        # here is the main constructor of the Parameter JuMP packages
        # It will create model *parameters* instead of variables
        # Variables are added to the optimization model, while parameters
        # are not. Parameters are merged with LP problem constants and do not
        # increase the model dimensions.
        β = add_parameters(slave, zeros(N_Candidates))
    else
        # Create fixed variables
        @variables(slave, begin
            β[Candidates]
            β_fixed[1:N_Candidates] == 0
        end)
        @constraint(slave, β_fix[i in Candidates], β[i] == β_fixed[i])
    end

    # create local constraints
    # Note that *parameter* algebra is implemented just like variables
    # algebra. We can multiply parameters by constants, add parameters,
    # sum parameters and varaibles and so on.
    @constraints(slave, begin
        ɛ_up_ctr[i in ObsSet(K)],
            ɛ_up[i] >= + sum(X[j,i] * β[j] for j ∈ Candidates) - y[i]
        ɛ_dw_ctr[i in ObsSet(K)],
            ɛ_dw[i] >= - sum(X[j,i] * β[j] for j ∈ Candidates) + y[i]
    end)
    # ATTENTION β[j] * X[j,i] Is much slower

    # create local objective function
    @objective(slave, Min, sum(ɛ_up[i] + ɛ_dw[i] for i ∈ ObsSet(K)))

    # return the correct group of parameters
    if PARAM
        return (slave, β)#β_upper, β_lower)
    else
        return (slave, β, β_fixed, β_fix)#β_upper, β_lower)
    end
end

#' ### Master

#' Now that all pieces of the original problem can be representad by
#' the convex $z_k(x)$ functions we can recast the problem in the the equivalent form:
#'
#' $$
#'\begin{align}
#'  & \min_{x}          &&  c^T x + z_1(x) + \dots + z_n(x) && \notag \\
#'  & \text{subject to} &&  Ax = b                          && \notag \\
#'  &                   &&  x \geq 0                        && \notag \\
#'\end{align}
#' $$
#'
#' However we cannot pass a problem in this form to a linear programming
#' solver (it could be passed to other kinds of solvers).
#'
#' Another standart result of optimization theory is that a convex function
#' can be represented by its supporting hyper-planes:
#'
#' $$
#'\begin{align}
#'  z_k(x) \ = \ & \min_{z, x}       &&  z && \notag \\
#'               & \text{subject to} &&  z \geq \pi_k(\hat{x}) (x - \hat{x}) + z_k(\hat{x}), \ \forall \hat{x} \in dom(z_k) && \notag \\
#'\end{align}
#' $$
#'
#' Then we can re-write (again) the master problem as
#'
#' $$
#'\begin{align}
#'  & \min_{x, z_k}     &&  c^T x + z_1 + \dots + z_n \notag \\
#'  & \text{subject to} &&  z_i \geq \pi_i(\hat{x}) (x - \hat{x}) + z_i(\hat{x}), \ \forall \hat{x} \in dom(z_i), i \in \{1, \dots, n\} \notag \\
#'  &                   &&  Ax = b \notag \\
#'  &                   &&  x \geq 0 \notag \\
#'\end{align}
#' $$
#'
#' Which is a linear program!
#'
#' However, it has infinitely many constraints !!!
#'
#' We can relax the infinite constraints and write:
#'
#' $$
#'\begin{align}
#'  & \min_{x, z_k}     &&  c^T x + z_1 + \dots + z_n \notag \\
#'  & \text{subject to} &&  Ax = b \notag \\
#'  &                   &&  x \geq 0 \notag \\
#'\end{align}
#' $$
#'
#' But now its only an underestimated problem.
#' In the case of our problem it can be written as:
#'
#' $$
#'\begin{align}
#'  & \min_{\varepsilon, \beta} &&  \sum_{i \in Nodes} \varepsilon_i \notag \\
#'  & \text{subject to} &&  \varepsilon_i \geq 0 \notag \\
#'\end{align}
#' $$
#'
#' This model can be written in JuMP
#'

function master_model(PARAM)
    master = Model(with_optimizer(OPTIMIZER))
    @variables(master, begin
        ɛ[Nodes] >= 0
        β[1:N_Candidates]
        # 0 <= β[Candidates] <= 8
    end)
    @objective(master, Min, sum(ɛ[i] for i ∈ Nodes))
    sol = zeros(N_Candidates)
    return (master, ɛ, β, sol)
end

#' The method to solve the master problem and query its solution
#' is given here:

function master_solve(PARAM, master_model)
    model = master_model[1]
    β = master_model[3]
    optimize!(model)
    return (value.(β), objective_value(model))
end

#' ### Supporting Hyperplanes

#' With these building blocks in hand, we can start building the algorithm.
#'
#' So far we know how to:
#' - Solve the relaxed master problem
#' - Obtain the solution for the $\hat{x}$ (or $\beta$ in our case)
#'
#'
#' Now we can:
#' - Fix the values of $\hat{x}$ in the slave problems
#' - Solve the slave problems
#' - query the solution of the slave problems to obtain the supporting hyperplane
#'
#' the value of $z_k(\hat{x})$, which is the objective value of the slave problem
#'
#' and the derivative $\pi_k(\hat{x}) = \frac{d z_k(x)}{d x} \Big|_{x = \hat{x}}$
#' The derivative is the dual variable associated to the variable $\hat{x}$,
#' which results by applying the chain rule on the constraints duals.
#'
#' These new steps are executed by the function:

function slave_solve(PARAM, model, master_solution)
    β0 = master_solution[1]
    slave = model[1]

    # The first step is to fix the values given by the master problem
    @timeit "fix" if PARAM
        # *parameters* can be set to new values and the optimization
        # model will be automatically updated
        β = model[2]
        fix.(β, β0)
    else
        # JuMP also has the hability to fix variables to new values
        β_fixed = model[3]
        β_fix = model[4]
        fix.(β_fixed, β0)
    end

    # here the slave problem is solved
    @timeit "opt" optimize!(slave)

    # query dual variables, which are sensitivities
    # They represent the subgradient (almost a derivative)
    # of the objective function for infinitesimal variations
    # of the constants in the linear constraints
    @timeit "dual" if PARAM
        # we can query dual values of *parameters*
        π = dual.(β)
    else
        # or, in pure JuMP, we query the duals form
        # constraints that fix the values of our regression
        # coefficients
        π = dual.(β_fix)
    end

    # π2 = shadow_price.(β_fix)
    # @show sum(π .- π2)
    obj = objective_value(slave)
    rhs = obj - dot(π, β0)
    return (rhs, π, obj)
end

#'
#' Now that we have cutting plane in hand we can add them to the master problem:
#'

function master_add_cut(PARAM, master_model, cut_info, node)
    master = master_model[1]
    ɛ = master_model[2]
    β = master_model[3]

    rhs = cut_info[1]
    π = cut_info[2]

    @constraint(master,
        ɛ[node] >= sum(π[j] * β[j] for j ∈ Candidates) + rhs)
end

#' ### Algorithm wrap up

#'
#' The complete algorithm is
#'
#'
#' - Solve the relaxed master problem
#' - Obtain the solution for the $\hat{x}$ (or $\beta$ in our case)
#' - Fix the values of $\hat{x}$ in the slave problems
#' - Solve the slave problem
#' - query the solution of the slave problem to obtain the supporting hyperplane
#' - add hyperplane to master problem
#' - repeat
#'

#' Now we grab all the pieces that we built and we write the benders
#' algorithm by calling the above function in a proper order.
#'
#' The macros `@timeit` are use to time each step of the algorithm.

function decomposed_model(PARAM)
    reset_timer!() # reset timer fo comparision
    time_init = @elapsed @timeit "Init" begin
        println("Initialize decomposed model")

        # Create the master problem with no cuts
        println("Build master problem")
        @timeit "Master" master = master_model(PARAM)

        # initialize solution for the regression coefficients in zero
        println("Build initial solution")
        @timeit "Sol" solution = (zeros(N_Candidates), Inf)
        best_sol = deepcopy(solution)

        # Create the slave problems
        println("Build slave problems")
        @timeit "Slaves" slaves = [slave_model(PARAM, i) for i ∈ Candidates]

        # Save initial version of the slave problems and create
        # the first set of cuts
        println("Build initial cuts")
        @timeit "Cuts" cuts = [slave_solve(PARAM, slaves[i], solution) for i ∈ Candidates]
    end

    UB = +Inf
    LB = -Inf

    println("Initialize Iterative step")
    time_loop = @elapsed @timeit "Loop"  for k in 1:80

        # Add cuts generated from each slave problem to the master problem
        @timeit "add cuts" for i ∈ Candidates
            master_add_cut(PARAM, master, cuts[i], i)
        end

        # Solve the master problem with the new set of cuts
        # Obtain new solution candidate for the regression coefficients
        @timeit "solve master" solution = master_solve(PARAM, master)

        # Pass the new candidate solution to each of the slave problems
        # Solve the slave problems and obtain cutting planes
        # @show solution[2]
        @timeit "solve nodes" for i ∈ Candidates
            cuts[i] = slave_solve(PARAM, slaves[i], solution)
        end

        LB = solution[2]
        new_UB = sum(cuts[i][3] for i ∈ Candidates)
        if new_UB <= UB
            best_sol = deepcopy(solution)
        end
        UB = min(UB, new_UB)
        println("Iter = $k, LB = $LB, UB = $UB")

        if abs(UB - LB)/(abs(UB)+abs(LB)) < 0.05
            println("Converged!")
            break
        end
    end
    println("First coefficients in solution: $(solution[1][1:min(10, N_Candidates)])")
    println("Objective value: $(solution[2])")
    println("Time in loop: $time_loop")
    println("Time in init: $time_init")

    print_timer()

    return best_sol[1]
end

#' Run benders decomposition with pure JuMP

GC.gc()
β1 = decomposed_model(false);

#' Run benders decomposition with ParameterJuMP

GC.gc()
β2 = decomposed_model(true);

#' Plot resulting time series from the benders base estimations

const y1 = X' * β1
const y2 = X' * β2

plt = plot(time, y,
    xlabel = "Time (s)", ylabel = "Amplitude")
plot!(plt, time, y1)
plot!(plt, time, y2)

#' # Acknowledgments

#' ParameterJuMP was developed by Joaquim Dias Garcia (@joaquimg) and Benoît Legat (@blegat)

#+ echo=false
# using Weave; cd("C:\\Users\\joaquimgarcia\\.julia\\dev\\ParameterJuMP\\examples"); weave("benders_quatile_regression.jl"); weave("benders_quatile_regression.jl", doctype = "github");# weave("benders_quatile_regression.jl", doctype = "md2pdf")
#+


text:
export between, P, O, X, Y, xs, ys, magnitude, normalize, signed_angle_to

X(x::Real) = P(x, 0)
Y(y::Real) = P(0, y)
const O = P(0, 0)

xs(ps::Array{Point}) = [p.x for p in ps]
ys(ps::Array{Point}) = [p.y for p in ps]

Point(ang::Angle) = Point(cos(ang.rad), sin(ang.rad))
Point(ang::Angle, length::Real) = Point(cos(ang.rad) * length, sin(ang.rad) * length)

Base.convert(::Type{Point}, t::Tuple{S, T}) where {S<:Real,T<:Real} = Point(t[1], t[2])

Base.show(io::IO, p::Point) = print(io, "Point($(p.xy[1]), $(p.xy[2]))")

function needed_attributes(::Type{Point})
    (Visible, Fill, Stroke, Linestyle, Linewidth, Markersize, Marker)
end

magnitude(p::Point) = sqrt(sum(p.xy .^ 2))
normalize(p::Point) = p / magnitude(p)

function Base.getproperty(p::Point, sym::Symbol)
    if sym == :x
        return p.xy[1]
    elseif sym == :y
        return p.xy[2]
    else
        getfield(p, sym)
    end
end

Base.convert(::Type{SVector{2, Float64}}, p::Point) = p.xy

Base.:+(p1::Point, p2::Point) = Point(p1.xy + p2.xy)
Base.:-(p1::Point, p2::Point) = Point(p1.xy - p2.xy)
Base.:-(p::Point) = Point(-p.xy)
Base.:*(p::Point, factor::Real) = Point(p.xy .* factor)
Base.:*(factor::Real, p::Point) = p * factor
Base.:/(p::Point, r::Real) = Point(p.xy ./ r)
Base.:(==)(p1::Point, p2::Point) = p1.xy == p2.xy
Base.:(!=)(p1::Point, p2::Point) = p1.xy != p2.xy
from_to(p1::Point, p2::Point) = p2 - p1
const → = from_to
between(p1::Point, p2::Point, fraction::Real) = p1 + (p2 - p1) * fraction
cross(p1::Point, p2::Point) = p1.x * p2.y - p1.y * p2.x
dot(p1::Point, p2::Point) = p1.x * p2.x + p1.y * p2.y

Base.isapprox(p1::Point, p2::Point) = Base.isapprox(magnitude(p1 → p2), 0; atol=6)

function Base.angle(p::Point)
    Angle(atan(p.y, p.x))
end

function signed_angle_to(p1::Point, p2::Point)
    Angle(atan(cross(p1, p2), dot(p1, p2)))
end

function _rotation_matrix(angle::Angle)
    c = cos(angle)
    s = sin(angle)
    SMatrix{2, 2, Float64}(c, s, -s, c)
end

function rotate(p::Point, angle::Angle; around::Point=Point(0, 0))
    vector = from_to(around, p)
    rotated_vector = Point(_rotation_matrix(angle) * vector.xy)
    rotated_vector + around
end


text:
using EvolutionaryModelingTools
using EvolutionaryModelingTools: sample, collectargs, _parsefunc
using Test

@testset "EvolutionaryModelingTools" begin
    @testset "Examples" begin
        # Examples are test with hand-code simulations
        @testset "Logistic growth with noise" begin
            include("example/logistic.jl")
        end
        @testset "Competitive Dynamics with noise and mutation" begin
            include("example/competition.jl")
        end
        @testset "SIR model with noise" begin
            include("example/sir.jl")
        end
    end

    @testset "Sample" begin
        @test sample(1, 1) == CartesianIndex()
        @test sample(fill(1), 1) == CartesianIndex()
        @test sample([1, 2], 2) == CartesianIndex(2)
        @test sample([1 2], 2) == CartesianIndex(1, 2)
        @test sample(reshape(1:4, 2, 2), 1) == CartesianIndex(1, 1)
        @test sample(reshape(1:4, 2, 2), 2) == CartesianIndex(2, 1)
        @test sample(reshape(1:4, 2, 2), 4) == CartesianIndex(1, 2)
        @test sample(reshape(1:4, 2, 2), 7) == CartesianIndex(2, 2)
    end

    @testset "Tools" begin
        @testset "parse function" begin
            fn, args =
                _parsefunc(:(@inline f(rng, ind, v) = v[ind] = rand(rng)), (:rng, :ind))
            @test fn == :f
            @test args == [:rng, :ind, :(args.v)]
        end

        @testset "collectargs" begin
            @test collectargs(quote
                global exp
                iv = f(arg) # internal variable
                mpd = map(exp, iv) # global variable
                A[1] = 1 # setindex
                A[begin+1:end] # setindex
                A[(i = 1; begin+i:end-i)]
                i # i is a local variable defined in above expression
                A[begin+j:end-j] # j is a arguments needed
                arg.x
                g(x) = x + 1 # local function
                (x, y) -> x + y # anonymous function
                @inbounds iv[ind] # macro call
                rand(rng)::Float64 # type annotation
            end) == Set([:A, :arg, :j, :ind, :rng])
            @test collectargs(quote
                iv = f(arg) # internal variable
                mpd = map(exp, iv) # global variable
                A[1] = 1 # setindex
                A[begin+1:end] # setindex
                arg.x
                g(x) = x + 1 # local function
                (x, y) -> x + y # anonymous function
                @inbounds iv[ind] # macro call
                rand(rng)::Float64 # type annotation
            end) == Set([:A, :arg, :exp, :ind, :rng]) # if exp is not mark as global, it will be collected
            @test collectargs(:((x::Int, y) -> x * y)) == Set([:x, :y])
            @test collectargs(:((x, y) -> @ein C[i, j] := x[i] * y[j])) == Set([:x, :y])
        end
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct HydroPumpedStorage <: HydroGen
        name::String
        available::Bool
        bus::Bus
        active_power::Float64
        reactive_power::Float64
        rating::Float64
        base_power::Float64
        prime_mover::PrimeMovers
        active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
        ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        rating_pump::Float64
        active_power_limits_pump::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits_pump::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
        ramp_limits_pump::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        time_limits_pump::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        storage_capacity::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
        inflow::Float64
        outflow::Float64
        initial_storage::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
        storage_target::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
        operation_cost::OperationalCost
        pump_efficiency::Float64
        conversion_factor::Float64
        time_at_status::Float64
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `active_power::Float64`
- `reactive_power::Float64`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: `(0, nothing)`, action if invalid: `error`
- `base_power::Float64`: Base power of the unit in MVA, validation range: `(0, nothing)`, action if invalid: `warn`
- `prime_mover::PrimeMovers`: Prime mover technology according to EIA 923
- `active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`, validation range: `(0, nothing)`, action if invalid: `warn`
- `reactive_power_limits::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}`, action if invalid: `warn`
- `ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: ramp up and ramp down limits in MW (in component base per unit) per minute, validation range: `(0, nothing)`, action if invalid: `error`
- `time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: Minimum up and Minimum down time limits in hours, validation range: `(0, nothing)`, action if invalid: `error`
- `rating_pump::Float64`: Thermal limited MVA Power Withdrawl of the pump. <= Capacity, validation range: `(0, nothing)`, action if invalid: `error`
- `active_power_limits_pump::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `reactive_power_limits_pump::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}`, action if invalid: `warn`
- `ramp_limits_pump::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: ramp up and ramp down limits in MW (in component base per unit) per minute of pump, validation range: `(0, nothing)`, action if invalid: `error`
- `time_limits_pump::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: Minimum up and Minimum down time limits of pump in hours, validation range: `(0, nothing)`, action if invalid: `error`
- `storage_capacity::NamedTuple{(:up, :down), Tuple{Float64, Float64}}`: Maximum storage capacity in the upper and lower reservoirs (units can be p.u-hr or m^3)., validation range: `(0, nothing)`, action if invalid: `error`
- `inflow::Float64`: Baseline inflow into the upper reservoir (units can be p.u. or m^3/hr), validation range: `(0, nothing)`, action if invalid: `error`
- `outflow::Float64`: Baseline outflow from the lower reservoir (units can be p.u. or m^3/hr), validation range: `(0, nothing)`, action if invalid: `error`
- `initial_storage::NamedTuple{(:up, :down), Tuple{Float64, Float64}}`: Initial storage capacity in the upper and lower reservoir (units can be p.u-hr or m^3)., validation range: `(0, nothing)`, action if invalid: `error`
- `storage_target::NamedTuple{(:up, :down), Tuple{Float64, Float64}}`: Storage target of upper reservoir at the end of simulation as ratio of storage capacity.
- `operation_cost::OperationalCost`: Operation Cost of Generation [`OperationalCost`](@ref)
- `pump_efficiency::Float64`: Efficiency of pump, validation range: `(0, 1)`, action if invalid: `warn`
- `conversion_factor::Float64`: Conversion factor from flow/volume to energy: m^3 -> p.u-hr.
- `time_at_status::Float64`
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct HydroPumpedStorage <: HydroGen
    name::String
    available::Bool
    bus::Bus
    active_power::Float64
    reactive_power::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "Base power of the unit in MVA"
    base_power::Float64
    "Prime mover technology according to EIA 923"
    prime_mover::PrimeMovers
    active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
    "ramp up and ramp down limits in MW (in component base per unit) per minute"
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Minimum up and Minimum down time limits in hours"
    time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Thermal limited MVA Power Withdrawl of the pump. <= Capacity"
    rating_pump::Float64
    active_power_limits_pump::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits_pump::Union{Nothing, NamedTuple{(:min, :max), Tuple{Float64, Float64}}}
    "ramp up and ramp down limits in MW (in component base per unit) per minute of pump"
    ramp_limits_pump::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Minimum up and Minimum down time limits of pump in hours"
    time_limits_pump::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Maximum storage capacity in the upper and lower reservoirs (units can be p.u-hr or m^3)."
    storage_capacity::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
    "Baseline inflow into the upper reservoir (units can be p.u. or m^3/hr)"
    inflow::Float64
    "Baseline outflow from the lower reservoir (units can be p.u. or m^3/hr)"
    outflow::Float64
    "Initial storage capacity in the upper and lower reservoir (units can be p.u-hr or m^3)."
    initial_storage::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
    "Storage target of upper reservoir at the end of simulation as ratio of storage capacity."
    storage_target::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
    "Operation Cost of Generation [`OperationalCost`](@ref)"
    operation_cost::OperationalCost
    "Efficiency of pump"
    pump_efficiency::Float64
    "Conversion factor from flow/volume to energy: m^3 -> p.u-hr."
    conversion_factor::Float64
    time_at_status::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function HydroPumpedStorage(name, available, bus, active_power, reactive_power, rating, base_power, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, rating_pump, active_power_limits_pump, reactive_power_limits_pump, ramp_limits_pump, time_limits_pump, storage_capacity, inflow, outflow, initial_storage, storage_target=(up=1.0, down=1.0), operation_cost=TwoPartCost(0.0, 0.0), pump_efficiency=1.0, conversion_factor=1.0, time_at_status=INFINITE_TIME, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    HydroPumpedStorage(name, available, bus, active_power, reactive_power, rating, base_power, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, rating_pump, active_power_limits_pump, reactive_power_limits_pump, ramp_limits_pump, time_limits_pump, storage_capacity, inflow, outflow, initial_storage, storage_target, operation_cost, pump_efficiency, conversion_factor, time_at_status, services, dynamic_injector, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function HydroPumpedStorage(; name, available, bus, active_power, reactive_power, rating, base_power, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, rating_pump, active_power_limits_pump, reactive_power_limits_pump, ramp_limits_pump, time_limits_pump, storage_capacity, inflow, outflow, initial_storage, storage_target=(up=1.0, down=1.0), operation_cost=TwoPartCost(0.0, 0.0), pump_efficiency=1.0, conversion_factor=1.0, time_at_status=INFINITE_TIME, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    HydroPumpedStorage(name, available, bus, active_power, reactive_power, rating, base_power, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, rating_pump, active_power_limits_pump, reactive_power_limits_pump, ramp_limits_pump, time_limits_pump, storage_capacity, inflow, outflow, initial_storage, storage_target, operation_cost, pump_efficiency, conversion_factor, time_at_status, services, dynamic_injector, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function HydroPumpedStorage(::Nothing)
    HydroPumpedStorage(;
        name="init",
        available=false,
        bus=Bus(nothing),
        active_power=0.0,
        reactive_power=0.0,
        rating=0.0,
        base_power=0.0,
        prime_mover=PrimeMovers.HY,
        active_power_limits=(min=0.0, max=0.0),
        reactive_power_limits=nothing,
        ramp_limits=nothing,
        time_limits=nothing,
        rating_pump=0.0,
        active_power_limits_pump=(min=0.0, max=0.0),
        reactive_power_limits_pump=nothing,
        ramp_limits_pump=nothing,
        time_limits_pump=nothing,
        storage_capacity=(up=0.0, down=0.0),
        inflow=0.0,
        outflow=0.0,
        initial_storage=(up=0.0, down=0.0),
        storage_target=(up=0.0, down=0.0),
        operation_cost=TwoPartCost(nothing),
        pump_efficiency=0.0,
        conversion_factor=0.0,
        time_at_status=INFINITE_TIME,
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

InfrastructureSystems.get_name(value::HydroPumpedStorage) = value.name
"""Get [`HydroPumpedStorage`](@ref) `available`."""
get_available(value::HydroPumpedStorage) = value.available
"""Get [`HydroPumpedStorage`](@ref) `bus`."""
get_bus(value::HydroPumpedStorage) = value.bus
"""Get [`HydroPumpedStorage`](@ref) `active_power`."""
get_active_power(value::HydroPumpedStorage) = get_value(value, value.active_power)
"""Get [`HydroPumpedStorage`](@ref) `reactive_power`."""
get_reactive_power(value::HydroPumpedStorage) = get_value(value, value.reactive_power)
"""Get [`HydroPumpedStorage`](@ref) `rating`."""
get_rating(value::HydroPumpedStorage) = get_value(value, value.rating)
"""Get [`HydroPumpedStorage`](@ref) `base_power`."""
get_base_power(value::HydroPumpedStorage) = value.base_power
"""Get [`HydroPumpedStorage`](@ref) `prime_mover`."""
get_prime_mover(value::HydroPumpedStorage) = value.prime_mover
"""Get [`HydroPumpedStorage`](@ref) `active_power_limits`."""
get_active_power_limits(value::HydroPumpedStorage) = get_value(value, value.active_power_limits)
"""Get [`HydroPumpedStorage`](@ref) `reactive_power_limits`."""
get_reactive_power_limits(value::HydroPumpedStorage) = get_value(value, value.reactive_power_limits)
"""Get [`HydroPumpedStorage`](@ref) `ramp_limits`."""
get_ramp_limits(value::HydroPumpedStorage) = get_value(value, value.ramp_limits)
"""Get [`HydroPumpedStorage`](@ref) `time_limits`."""
get_time_limits(value::HydroPumpedStorage) = value.time_limits
"""Get [`HydroPumpedStorage`](@ref) `rating_pump`."""
get_rating_pump(value::HydroPumpedStorage) = get_value(value, value.rating_pump)
"""Get [`HydroPumpedStorage`](@ref) `active_power_limits_pump`."""
get_active_power_limits_pump(value::HydroPumpedStorage) = get_value(value, value.active_power_limits_pump)
"""Get [`HydroPumpedStorage`](@ref) `reactive_power_limits_pump`."""
get_reactive_power_limits_pump(value::HydroPumpedStorage) = get_value(value, value.reactive_power_limits_pump)
"""Get [`HydroPumpedStorage`](@ref) `ramp_limits_pump`."""
get_ramp_limits_pump(value::HydroPumpedStorage) = get_value(value, value.ramp_limits_pump)
"""Get [`HydroPumpedStorage`](@ref) `time_limits_pump`."""
get_time_limits_pump(value::HydroPumpedStorage) = value.time_limits_pump
"""Get [`HydroPumpedStorage`](@ref) `storage_capacity`."""
get_storage_capacity(value::HydroPumpedStorage) = get_value(value, value.storage_capacity)
"""Get [`HydroPumpedStorage`](@ref) `inflow`."""
get_inflow(value::HydroPumpedStorage) = get_value(value, value.inflow)
"""Get [`HydroPumpedStorage`](@ref) `outflow`."""
get_outflow(value::HydroPumpedStorage) = value.outflow
"""Get [`HydroPumpedStorage`](@ref) `initial_storage`."""
get_initial_storage(value::HydroPumpedStorage) = get_value(value, value.initial_storage)
"""Get [`HydroPumpedStorage`](@ref) `storage_target`."""
get_storage_target(value::HydroPumpedStorage) = value.storage_target
"""Get [`HydroPumpedStorage`](@ref) `operation_cost`."""
get_operation_cost(value::HydroPumpedStorage) = value.operation_cost
"""Get [`HydroPumpedStorage`](@ref) `pump_efficiency`."""
get_pump_efficiency(value::HydroPumpedStorage) = value.pump_efficiency
"""Get [`HydroPumpedStorage`](@ref) `conversion_factor`."""
get_conversion_factor(value::HydroPumpedStorage) = value.conversion_factor
"""Get [`HydroPumpedStorage`](@ref) `time_at_status`."""
get_time_at_status(value::HydroPumpedStorage) = value.time_at_status
"""Get [`HydroPumpedStorage`](@ref) `services`."""
get_services(value::HydroPumpedStorage) = value.services
"""Get [`HydroPumpedStorage`](@ref) `dynamic_injector`."""
get_dynamic_injector(value::HydroPumpedStorage) = value.dynamic_injector
"""Get [`HydroPumpedStorage`](@ref) `ext`."""
get_ext(value::HydroPumpedStorage) = value.ext

InfrastructureSystems.get_time_series_container(value::HydroPumpedStorage) = value.time_series_container
"""Get [`HydroPumpedStorage`](@ref) `internal`."""
get_internal(value::HydroPumpedStorage) = value.internal

InfrastructureSystems.set_name!(value::HydroPumpedStorage, val) = value.name = val
"""Set [`HydroPumpedStorage`](@ref) `available`."""
set_available!(value::HydroPumpedStorage, val) = value.available = val
"""Set [`HydroPumpedStorage`](@ref) `bus`."""
set_bus!(value::HydroPumpedStorage, val) = value.bus = val
"""Set [`HydroPumpedStorage`](@ref) `active_power`."""
set_active_power!(value::HydroPumpedStorage, val) = value.active_power = val
"""Set [`HydroPumpedStorage`](@ref) `reactive_power`."""
set_reactive_power!(value::HydroPumpedStorage, val) = value.reactive_power = val
"""Set [`HydroPumpedStorage`](@ref) `rating`."""
set_rating!(value::HydroPumpedStorage, val) = value.rating = val
"""Set [`HydroPumpedStorage`](@ref) `base_power`."""
set_base_power!(value::HydroPumpedStorage, val) = value.base_power = val
"""Set [`HydroPumpedStorage`](@ref) `prime_mover`."""
set_prime_mover!(value::HydroPumpedStorage, val) = value.prime_mover = val
"""Set [`HydroPumpedStorage`](@ref) `active_power_limits`."""
set_active_power_limits!(value::HydroPumpedStorage, val) = value.active_power_limits = val
"""Set [`HydroPumpedStorage`](@ref) `reactive_power_limits`."""
set_reactive_power_limits!(value::HydroPumpedStorage, val) = value.reactive_power_limits = val
"""Set [`HydroPumpedStorage`](@ref) `ramp_limits`."""
set_ramp_limits!(value::HydroPumpedStorage, val) = value.ramp_limits = val
"""Set [`HydroPumpedStorage`](@ref) `time_limits`."""
set_time_limits!(value::HydroPumpedStorage, val) = value.time_limits = val
"""Set [`HydroPumpedStorage`](@ref) `rating_pump`."""
set_rating_pump!(value::HydroPumpedStorage, val) = value.rating_pump = val
"""Set [`HydroPumpedStorage`](@ref) `active_power_limits_pump`."""
set_active_power_limits_pump!(value::HydroPumpedStorage, val) = value.active_power_limits_pump = val
"""Set [`HydroPumpedStorage`](@ref) `reactive_power_limits_pump`."""
set_reactive_power_limits_pump!(value::HydroPumpedStorage, val) = value.reactive_power_limits_pump = val
"""Set [`HydroPumpedStorage`](@ref) `ramp_limits_pump`."""
set_ramp_limits_pump!(value::HydroPumpedStorage, val) = value.ramp_limits_pump = val
"""Set [`HydroPumpedStorage`](@ref) `time_limits_pump`."""
set_time_limits_pump!(value::HydroPumpedStorage, val) = value.time_limits_pump = val
"""Set [`HydroPumpedStorage`](@ref) `storage_capacity`."""
set_storage_capacity!(value::HydroPumpedStorage, val) = value.storage_capacity = val
"""Set [`HydroPumpedStorage`](@ref) `inflow`."""
set_inflow!(value::HydroPumpedStorage, val) = value.inflow = val
"""Set [`HydroPumpedStorage`](@ref) `outflow`."""
set_outflow!(value::HydroPumpedStorage, val) = value.outflow = val
"""Set [`HydroPumpedStorage`](@ref) `initial_storage`."""
set_initial_storage!(value::HydroPumpedStorage, val) = value.initial_storage = val
"""Set [`HydroPumpedStorage`](@ref) `storage_target`."""
set_storage_target!(value::HydroPumpedStorage, val) = value.storage_target = val
"""Set [`HydroPumpedStorage`](@ref) `operation_cost`."""
set_operation_cost!(value::HydroPumpedStorage, val) = value.operation_cost = val
"""Set [`HydroPumpedStorage`](@ref) `pump_efficiency`."""
set_pump_efficiency!(value::HydroPumpedStorage, val) = value.pump_efficiency = val
"""Set [`HydroPumpedStorage`](@ref) `conversion_factor`."""
set_conversion_factor!(value::HydroPumpedStorage, val) = value.conversion_factor = val
"""Set [`HydroPumpedStorage`](@ref) `time_at_status`."""
set_time_at_status!(value::HydroPumpedStorage, val) = value.time_at_status = val
"""Set [`HydroPumpedStorage`](@ref) `services`."""
set_services!(value::HydroPumpedStorage, val) = value.services = val
"""Set [`HydroPumpedStorage`](@ref) `ext`."""
set_ext!(value::HydroPumpedStorage, val) = value.ext = val

InfrastructureSystems.set_time_series_container!(value::HydroPumpedStorage, val) = value.time_series_container = val


text:
export eventpaths, eventstream
import SeisIO: merge,merge!

"""
     eventpaths()

Use S3 to query SCEDC-pds event paths. Returns array of event S3 paths.

# Keywords
- `starttime::TimeType`: The start time/day of the event query.
- `endtime::TimeType`: The end time/day of the event query.
- `minlatitude::Real`: Minimum latitude in event query.
- `maxlatitude::Real`: Maximum latitude in event query.
- `minlongitude::Real`: Minimum longitude in event query.
- `maxlongitude::Real`: Maximum longitude in event query.
- `minmagnitude::Real`: Minimum magnitude in event query.
- `maxmagnitude::Real`: Maximum magnitude in event query.
- `mindepth::Real`: Minimum event depth in km.
- `maxdepth::Real`: Maximum event depth in km.
- `eventtype::String`: Event type: {"eq"=>"earthquake", "qb"=> "quarry blast",
"sn"=>"sonic boom', "nt"=>"nuclear blast", "uk"=>"not reported","*"=>"all events"}
"""
function eventpaths(
    aws::AWSConfig;
    starttime::TimeType=today(),
    endtime::TimeType=today(),
    minlatitude::Real=-90,
    maxlatitude::Real=90,
    minlongitude::Real=-180,
    maxlongitude::Real=180,
    minmagnitude::Real=-10,
    maxmagnitude::Real=10,
    mindepth::Real=-10,
    maxdepth::Real=100,
    eventtype::String="*",
)
    @assert starttime >= Date(1977,1,1) "Starttime must be greater or equal to $(Date(1977,1,1))"
    # quert possible events from catalog
    eventdf = catalogquery(
        aws,
        starttime=starttime,
        endtime=endtime,
        minlatitude=minlatitude,
        maxlatitude=maxlatitude,
        minlongitude=minlongitude,
        maxlongitude=maxlongitude,
        minmagnitude=minmagnitude,
        mindepth=mindepth,
        maxdepth=maxdepth,
        eventtype=eventtype
    )

    N = size(eventdf,1)
    epaths = ["event_waveforms/$(year(row.EVENTTIME))/$(year(row.EVENTTIME))_$(lpad(dayofyear(row.EVENTTIME),3,'0'))/$(row.EVID).ms" for row in eachrow(eventdf)]
    ppaths = ["event_phases/$(year(row.EVENTTIME))/$(year(row.EVENTTIME))_$(lpad(dayofyear(row.EVENTTIME),3,'0'))/$(row.EVID).phase" for row in eachrow(eventdf)]
    paths = [epaths[ii] for ii = 1:length(epaths) if (s3_exists(aws,"scedc-pds",epaths[ii])) & (s3_exists(aws,"scedc-pds",ppaths[ii]))]

    # check that events actually exist
    if isempty(paths)
        throw(ErrorException("No events for query."))
    end
    return paths
end
eventpaths(a...;b...) = eventpaths(global_aws_config(region="us-west-2"),a...;b...)
eventpaths(d::Dict,a...;b...) = eventpaths(global_aws_config(region=d[:region]),a...;b...)

function eventdownload()
end

function eventstream(
    aws::AWSConfig;
    starttime::TimeType=today(),
    endtime::TimeType=today(),
    minlatitude::Real=-90,
    maxlatitude::Real=90,
    minlongitude::Real=-180,
    maxlongitude::Real=180,
    minmagnitude::Real=-10,
    maxmagnitude::Real=10,
    mindepth::Real=-10,
    maxdepth::Real=100,
    eventtype::String="*",
)
    s3paths = eventpaths(
        aws,
        starttime=starttime,
        endtime=endtime,
        minlatitude=minlatitude,
        maxlatitude=maxlatitude,
        minlongitude=minlongitude,
        maxlongitude=maxlongitude,
        minmagnitude=minmagnitude,
        mindepth=mindepth,
        maxdepth=maxdepth,
        eventtype=eventtype
    )

    # grab event id and date from s3path
    evid = [parse(Int,replace(basename(s),".ms"=>"")) for s in s3paths]
    s3dates = date_yyyyddd.([basename(dirname(s)) for s in s3paths])
    s3phases = [phasequery(evid[ii],s3dates[ii]) for ii = 1:length(evid)]
    s3events = ec2stream(aws,"scedc-pds",s3paths,rtype=Array)

    for ii = 1:length(s3events)
        s3phases[ii] = merge(s3phases[ii],s3events[ii])
    end
    return s3phases
end
eventstream(a...;b...) = eventstream(global_aws_config(region="us-west-2"), a...; b...)
eventstream(d::Dict,a...;b...) = eventstream(global_aws_config(region=d[:region]), a...; b...)

function date_yyyyddd(yearday::String)
    @assert occursin(r"[1-2][0-9][0-9][0-9]_[0-3][0-6][0-9]",yearday)
    y,d = split(yearday,"_")
    yint = parse(Int,y)
    dint = parse(Int,d)
    @assert dint <= 366 "Input day must be less than or equal to 366"
    return Date(yint) + Day(dint-1)
end

function merge!(EC::EventChannel,SC::SeisChannel)
    if EC.id != SC.id
        throw(ErrorException("EventChannel ID $(EC.id) and SeisChannel ID $(SC.id) do not match."))
    end
    for f in SeisIO.datafields
        fval = getfield(SC,f)
        if !isempty(fval)
            setfield!(EC,f,fval)
        end
    end
    return nothing
end
merge!(SC::SeisChannel,EC::EventChannel,) = merge!(EC,SC)
merge(EC::EventChannel,SC::SeisChannel) = (U = deepcopy(EC);merge!(U,SC);return U)
merge(SC::SeisChannel,EC::EventChannel) = merge(EC,SC)

function merge(SE::SeisEvent,SD::SeisData)
    newEvent = SeisEvent()
    newEvent.hdr = SE.hdr
    newEvent.source = SE.source
    newEvent.data = merge(SE.data,SD)
    return newEvent
end
merge(SD::SeisData,SE::SeisEvent) = merge(SE,SD)

function merge(ETD::EventTraceData,SD::SeisData)
    N = length(SD.id)
    ETDdict = Dict(ETD.id[ii]=>ii for ii = 1:length(ETD.id))
    ETDkeys = keys(ETDdict)
    newETD = EventTraceData(N)
    for ii = 1:N
        if in(SD.id[ii],ETDkeys)
            ind = ETDdict[SD.id[ii]]
            newETD[ii] = merge(SD[ii],ETD[ind])
        else
            newETD[ii] = EventTraceData(SD[ii])[1] # need to make EventChannel(SeisChannel())
        end
    end
    return newETD
end
merge(SD::SeisData,ETD::EventTraceData) = merge(ETD,SD)