text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using JuMP, LinearAlgebra, Geodesy, Cbc, Clp, ProgressBars, Printf, DataStructures

function build_model(instance::Instance, graph::Graph, optimizer)::JuMP.Model
    model = Model(optimizer)
    model[:instance] = instance
    model[:graph] = graph
    create_vars!(model)
    create_objective_function!(model)
    create_shipping_node_constraints!(model)
    create_process_node_constraints!(model)
    return model
end

function create_vars!(model::JuMP.Model)
    graph, T = model[:graph], model[:instance].time
    model[:flow] =
        Dict((a, t) => @variable(model, lower_bound = 0) for a in graph.arcs, t = 1:T)
    model[:dispose] = Dict(
        (n, t) => @variable(
            model,
            lower_bound = 0,
            upper_bound = n.location.disposal_limit[n.product][t]
        ) for n in values(graph.plant_shipping_nodes), t = 1:T
    )
    model[:store] = Dict(
        (n, t) =>
            @variable(model, lower_bound = 0, upper_bound = n.location.storage_limit)
        for n in values(graph.process_nodes), t = 1:T
    )
    model[:process] = Dict(
        (n, t) => @variable(model, lower_bound = 0) for
        n in values(graph.process_nodes), t = 1:T
    )
    model[:open_plant] = Dict(
        (n, t) => @variable(model, binary = true) for n in values(graph.process_nodes),
        t = 1:T
    )
    model[:is_open] = Dict(
        (n, t) => @variable(model, binary = true) for n in values(graph.process_nodes),
        t = 1:T
    )
    model[:capacity] = Dict(
        (n, t) => @variable(
            model,
            lower_bound = 0,
            upper_bound = n.location.sizes[2].capacity
        ) for n in values(graph.process_nodes), t = 1:T
    )
    model[:expansion] = Dict(
        (n, t) => @variable(
            model,
            lower_bound = 0,
            upper_bound = n.location.sizes[2].capacity - n.location.sizes[1].capacity
        ) for n in values(graph.process_nodes), t = 1:T
    )
end

function slope_open(plant, t)
    if plant.sizes[2].capacity <= plant.sizes[1].capacity
        0.0
    else
        (plant.sizes[2].opening_cost[t] - plant.sizes[1].opening_cost[t]) /
        (plant.sizes[2].capacity - plant.sizes[1].capacity)
    end
end

function slope_fix_oper_cost(plant, t)
    if plant.sizes[2].capacity <= plant.sizes[1].capacity
        0.0
    else
        (plant.sizes[2].fixed_operating_cost[t] - plant.sizes[1].fixed_operating_cost[t]) /
        (plant.sizes[2].capacity - plant.sizes[1].capacity)
    end
end

function create_objective_function!(model::JuMP.Model)
    graph, T = model[:graph], model[:instance].time
    obj = AffExpr(0.0)

    # Process node costs
    for n in values(graph.process_nodes), t = 1:T

        # Transportation and variable operating costs
        for a in n.incoming_arcs
            c = n.location.input.transportation_cost[t] * a.values["distance"]
            add_to_expression!(obj, c, model[:flow][a, t])
        end

        # Opening costs
        add_to_expression!(
            obj,
            n.location.sizes[1].opening_cost[t],
            model[:open_plant][n, t],
        )

        # Fixed operating costs (base)
        add_to_expression!(
            obj,
            n.location.sizes[1].fixed_operating_cost[t],
            model[:is_open][n, t],
        )

        # Fixed operating costs (expansion)
        add_to_expression!(obj, slope_fix_oper_cost(n.location, t), model[:expansion][n, t])

        # Processing costs
        add_to_expression!(
            obj,
            n.location.sizes[1].variable_operating_cost[t],
            model[:process][n, t],
        )

        # Storage costs
        add_to_expression!(obj, n.location.storage_cost[t], model[:store][n, t])

        # Expansion costs
        if t < T
            add_to_expression!(
                obj,
                slope_open(n.location, t) - slope_open(n.location, t + 1),
                model[:expansion][n, t],
            )
        else
            add_to_expression!(obj, slope_open(n.location, t), model[:expansion][n, t])
        end
    end

    # Shipping node costs
    for n in values(graph.plant_shipping_nodes), t = 1:T

        # Disposal costs
        add_to_expression!(
            obj,
            n.location.disposal_cost[n.product][t],
            model[:dispose][n, t],
        )
    end

    @objective(model, Min, obj)
end

function create_shipping_node_constraints!(model::JuMP.Model)
    graph, T = model[:graph], model[:instance].time
    model[:eq_balance] = OrderedDict()
    for t = 1:T
        # Collection centers
        for n in graph.collection_shipping_nodes
            model[:eq_balance][n, t] = @constraint(
                model,
                sum(model[:flow][a, t] for a in n.outgoing_arcs) == n.location.amount[t]
            )
        end

        # Plants
        for n in graph.plant_shipping_nodes
            @constraint(
                model,
                sum(model[:flow][a, t] for a in n.incoming_arcs) ==
                sum(model[:flow][a, t] for a in n.outgoing_arcs) + model[:dispose][n, t]
            )
        end
    end

end

function create_process_node_constraints!(model::JuMP.Model)
    graph, T = model[:graph], model[:instance].time

    for t = 1:T, n in graph.process_nodes
        input_sum = AffExpr(0.0)
        for a in n.incoming_arcs
            add_to_expression!(input_sum, 1.0, model[:flow][a, t])
        end

        # Output amount is implied by amount processed
        for a in n.outgoing_arcs
            @constraint(
                model,
                model[:flow][a, t] == a.values["weight"] * model[:process][n, t]
            )
        end

        # If plant is closed, capacity is zero
        @constraint(
            model,
            model[:capacity][n, t] <= n.location.sizes[2].capacity * model[:is_open][n, t]
        )

        # If plant is open, capacity is greater than base
        @constraint(
            model,
            model[:capacity][n, t] >= n.location.sizes[1].capacity * model[:is_open][n, t]
        )

        # Capacity is linked to expansion
        @constraint(
            model,
            model[:capacity][n, t] <=
            n.location.sizes[1].capacity + model[:expansion][n, t]
        )

        # Can only process up to capacity
        @constraint(model, model[:process][n, t] <= model[:capacity][n, t])

        if t > 1
            # Plant capacity can only increase over time
            @constraint(model, model[:capacity][n, t] >= model[:capacity][n, t-1])
            @constraint(model, model[:expansion][n, t] >= model[:expansion][n, t-1])
        end

        # Amount received equals amount processed plus stored
        store_in = 0
        if t > 1
            store_in = model[:store][n, t-1]
        end
        if t == T
            @constraint(model, model[:store][n, t] == 0)
        end
        @constraint(
            model,
            input_sum + store_in == model[:store][n, t] + model[:process][n, t]
        )

        # Plant is currently open if it was already open in the previous time period or
        # if it was built just now
        if t > 1
            @constraint(
                model,
                model[:is_open][n, t] == model[:is_open][n, t-1] + model[:open_plant][n, t]
            )
        else
            @constraint(model, model[:is_open][n, t] == model[:open_plant][n, t])
        end

        # Plant can only be opened during building period
        if t ∉ model[:instance].building_period
            @constraint(model, model[:open_plant][n, t] == 0)
        end
    end
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
"""
    ShallowWater{T<:Real,Dim}

The shallow water equations with gravitational constant `g`
in `Dim` dimensions using `T` as scalar type.
"""
struct ShallowWater{T,Dim} <: AbstractBalanceLaw{1}
  g::T
end

function ShallowWater(g=1., Dim=1)
  ShallowWater{typeof(g), Dim}(g)
end

function Base.show(io::IO, model::ShallowWater{T,Dim}) where {T,Dim}
  print(io, "Shallow water equations with g=", model.g, " {T=", T, ", Dim=", Dim, "}")
end

"""
Conserved variables (h, hv) of the shallow water system in one space dimension.
"""
struct ShallowWaterVar1D{T} <: FieldVector{2,T}
  h ::T
  hv::T
end

function (::Type{ShallowWaterVar1D{T}})(val::Real) where {T}
  ShallowWaterVar1D{T}(val, val)
end

function similar_type(::ShallowWaterVar1D{T}) where {T}
  ShallowWaterVar1D{T}
end

@inline variables(model::ShallowWater{T,1}) where {T} = ShallowWaterVar1D{T}

@inline function primitive_variables(u::ShallowWaterVar1D)
  @unpack h, hv = u
  if h ≈ 0
    v = zero(h)
  else
    v = hv/h
  end
  h, v
end

@inline function primitive_variables(u::ShallowWaterVar1D, model::ShallowWater)
  primitive_variables(u)
end

@inline function entropy_variables(u::ShallowWaterVar1D, model::ShallowWater)
    entropy_variables(primitive_variables(u, model)..., model)
end

@inline function entropy_variables(h, v, model::ShallowWater)
    @unpack g = model
    SVector(g*h - v^2/2, v)
end

@inline function flux_potential(u::ShallowWaterVar1D, model::ShallowWater)
    flux_potential(primitive_variables(u, model)..., model)
end

@inline function flux_potential(h, v, model::ShallowWater)
    @unpack g = model
    g*h^2*v/2
end

@inline function flux(u::ShallowWaterVar1D{T}, model::ShallowWater{T,1}) where {T}
  h, v = primitive_variables(u)
  @unpack g = model

  SVector{2,T}(h*v, h*v^2 + g*h^2/2)
end

@inline function max_abs_speed(u::ShallowWaterVar1D, model::ShallowWater{T,1}) where {T}
  h, v = primitive_variables(u)
  @unpack g = model

  abs(v) + sqrt(g*h)
end

"""
    max_abs_speed(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)

Estimate the maximal speed in the solution of the Riemann problem with left and
right states `uₗ`, `uᵣ` for `model` via the two rarefaction approximation, see
Guermond and Popov (2016) Fast Estimation from above of the maximum wave speed in
the Riemann roblem for the Euler equations, Journal of Computational Physics 321,
pp. 908-926, and
Chen and Shu (2017) Entropy stable high order discontinuous Galerkin methods with
suitable quadrature rules for hyperbolic conservation laws, Journal of Computational
Physics 345, pp. 427-461.
"""
@inline function max_abs_speed(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)
    λ₋, λ₊ = min_max_speed(uₗ, uᵣ, model)

    max(abs(λ₋), abs(λ₊))
end

"""
    min_max_speed(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)

Estimate the maximal speed in the solution of the Riemann problem with left and
right states `uₗ`, `uᵣ` for `model` via the two rarefaction approximation, see
Guermond and Popov (2016) Fast Estimation from above of the maximum wave speed in
the Riemann roblem for the Euler equations, Journal of Computational Physics 321,
pp. 908-926, and
Chen and Shu (2017) Entropy stable high order discontinuous Galerkin methods with
suitable quadrature rules for hyperbolic conservation laws, Journal of Computational
Physics 345, pp. 427-461.
"""
@inline function min_max_speed(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)
    @unpack g = model
    hₗ, vₗ = primitive_variables(uₗ, model)
    aₗ = sqrt(g * hₗ)
    hᵣ, vᵣ = primitive_variables(uᵣ, model)
    aᵣ = sqrt(g * hᵣ)

    hₘ = (2*(aᵣ+aₗ) - (vᵣ-vₗ))^2 / 16g
    if hₘ <= hₗ
        # approximated by rarefaction wave
        λ₋ = vₗ - aₗ
    else
        # approximated by shock wave
        λ₋ = vₗ - hₘ * sqrt( g*(hₘ+hₗ) / (2*hₘ*hₗ) )
    end
    if hₘ <= hᵣ
        # approximated by rarefaction wave
        λ₊ = vᵣ + aᵣ
    else
        # approximated by shock wave
        λ₊ = vᵣ + hₘ * sqrt( g*(hₘ+hᵣ) / (2*hₘ*hᵣ) )
    end

    λ₋, λ₊
end

"""
    (::EnergyConservativeFlux1Param)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)

The one-parameter family of energy conservative numerical fluxes for the shallow
water equations, see
Ranocha (2017) Shallow water equations: Split-form, entropy stable, well-balanced,
and positivity-preserving numerical methods, GEM - International Journal on
Geomathematics 8(1), pp. 85-133.
"""
function (fvol::EnergyConservativeFlux1Param)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D,
                                              model::ShallowWater)
    h₋, v₋ = primitive_variables(uₗ, model)
    h₊, v₊ = primitive_variables(uᵣ, model)
    @unpack g = model
    @unpack a₁ = fvol

    fnum_h  = ( (3-a₁)*(h₋*v₋ + h₊*v₊) + (1+a₁)*(h₊*v₋ + h₋*v₊) ) / 8
    fnum_hv = g * ( (1+a₁)*(h₊^2 + h₋^2) + (2-2a₁)*(h₊*h₋) ) / 8 +
                ( (3-a₁)*(h₊*v₊^2 + h₋*v₋^2) + (1+a₁)*(h₊*v₋^2 + h₋*v₊^2) ) / 16 +
                (h₊*v₊*v₋ + h₋*v₊*v₋) / 4

    SVector(fnum_h, fnum_hv)
end

"""
    (::EnergyConservativeFlux2Param)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)

The two-parameter family of energy conservative numerical fluxes for the shallow
water equations, see
Ranocha (2017) Shallow water equations: Split-form, entropy stable, well-balanced,
and positivity-preserving numerical methods, GEM - International Journal on
Geomathematics 8(1), pp. 85-133.
"""
function (fvol::EnergyConservativeFlux2Param)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D,
                                              model::ShallowWater)
    h₋, v₋ = primitive_variables(uₗ, model)
    h₊, v₊ = primitive_variables(uᵣ, model)
    @unpack g = model
    @unpack a₁, a₂ = fvol

    fnum_h  = ( (3-a₁)*(h₋*v₋ + h₊*v₊) + (1+a₁)*(h₊*v₋ + h₋*v₊) ) / 8 +
                (a₁+3a₂-2) * (v₊^3 - v₊^2*v₋ - v₊*v₋^2 + v₋^3) / 16g
    fnum_hv = g * ( (1+a₁)*(h₊^2 + h₋^2) + (2-2a₁)*(h₊*h₋) ) / 8 -
                (2a₁+3a₂-5) * (h₊*v₊^2 + h₋*v₋^2) / 16 +
                (2a₁+3a₂-1) * (h₊*v₋^2 + h₋*v₊^2) / 16 +
                (h₊*v₊*v₋ + h₋*v₊*v₋) / 4 +
                (a₁+3a₂-2) * ((v₊^2)^2 - 2*v₊^2*v₋^2 + (v₋^2)^2)

    SVector(fnum_h, fnum_hv)
end

"""
    (diss::ScalarDissipation)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)

The scalar dissipation operator \$- \\frac{\\lambda}{2} R \\cdot R^T \\cdot  (w_r - w_l)\$
for the shallow water equations, where \$\\partial_w u = R \\cdot R^T\$ is evaluated
at the arithmetic mean values \$h = (h_- + h_+) / 2\$, \$v = (v_- + v_+) / 2\$.
"""
function (diss::ScalarDissipation)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)
    h₋, v₋ = primitive_variables(uₗ, model)
    h₊, v₊ = primitive_variables(uᵣ, model)
    @unpack g = model

    w₋ = entropy_variables(h₋, v₋, model)
    w₊ = entropy_variables(h₊, v₊, model)

    h = (h₋ + h₊) / 2
    v = (v₋ + v₊) / 2
    du_dw = SArray{Tuple{2,2}}(1, v, v, g*h+v^2) / g

    λ = diss.max_abs_speed(uₗ, uᵣ, model)

    -λ/2 * du_dw * (w₊ - w₋)
end

"""
    (diss::MatrixDissipation)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)

The scalar dissipation operator \$- \\frac{1}{2} R \\cdot |\\Lambda| \\cdot R^T \\cdot  (w_r - w_l)\$
for the shallow water equations, where \$f'(u) = R \\cdot \\Lambda \\cdot R^{-1}\$ is
evaluated at the arithmetic mean values \$h = (h_- + h_+) / 2\$, \$v = (v_- + v_+) / 2\$.
"""
function (diss::MatrixDissipation)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater)
    h₋, v₋ = primitive_variables(uₗ, model)
    h₊, v₊ = primitive_variables(uᵣ, model)
    @unpack g = model

    w₋ = entropy_variables(h₋, v₋, model)
    w₊ = entropy_variables(h₊, v₊, model)

    h = (h₋ + h₊) / 2
    v = (v₋ + v₊) / 2
    sqrt_gh = sqrt(g*h)
    R = SArray{Tuple{2,2}}(1, v-sqrt_gh, 1, v+sqrt_gh) / sqrt(2g)
    Λ = SArray{Tuple{2,2}}(abs(v-sqrt_gh), 0, 0, abs(v+sqrt_gh))

    -R * Λ * R' * (w₊ - w₋) / 2
end

"""
    (::SuliciuFlux)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater{T,1})

The Suliciu relaxation solver / numerical flux, see
Bouchut (2004) Nonlinear Stability of Finite Volume Methods for Hyperbolic
Conservation Laws and Well-Balanced Schemes for Sources.
"""
function (::SuliciuFlux)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater{T,1}) where T
    hl, vl = primitive_variables(uₗ, model)
    hr, vr = primitive_variables(uᵣ, model)
    @unpack g = model

    sqrt_ghl = sqrt(g*hl)
    sqrt_ghr = sqrt(g*hr)

    # compute speeds
    cl_hl = zero(T)
    cr_hr = zero(T)
    if hl <= hr && 0 < hr
      cl_hl = sqrt_ghl + 1.5*max(0, 0.5*g*(hr*hr-hl*hl)/(hr*sqrt_ghr) + vl - vr )
      cr_hr = sqrt_ghr + 1.5*max(0, 0.5*g*(hl*hl-hr*hr)/(cl_hl*hl) + vl - vr )
    elseif hr <= hl && 0 < hl
      cr_hr = sqrt_ghr + 1.5*max(0, 0.5*g*(hl*hl-hr*hr)/(hl*sqrt_ghl) + vl - vr )
      cl_hl = sqrt_ghl + 1.5*max(0, 0.5*g*(hr*hr-hl*hl)/(cr_hr*hr) + vl - vr )
    else
      cl_hl = cr_hr = zero(T)
    end

    # compute intermediate values
    cl = cl_hl*hl
    cr = cr_hr*hr
    vls = vrs = ( cl*vl+cr*vr+0.5*g*(hl*hl-hr*hr) ) / ( cl+cr )
    vls = vrs = ifelse(isnan(vls), zero(T), vls)
    πls = πrs = ( 0.5*g*(cr*hl*hl+cl*hr*hr)-cl*cr*(vr-vl) ) / (cl+cr)
    πls = πrs = ifelse(isnan(πls), zero(T), πls)
    hls = 1 / ( 1/hl + (cr*(vr-vl)+0.5*g*(hl*hl-hr*hr))/(cl*(cl+cr)) )
    hls = ifelse(isnan(hls), zero(T), hls)
    hrs = 1 / ( 1/hr + (cl*(vr-vl)+0.5*g*(hr*hr-hl*hl))/(cr*(cl+cr)) )
    hrs = ifelse(isnan(hrs), zero(T), hrs)

    # compute fluxes
    fnum_h  = zero(T)
    fnum_hv = zero(T)
    if 0 <= vl-cl_hl
      fnum_h  = hl*vl
      fnum_hv = hl*vl*vl + 0.5*g*hl*hl
    elseif 0 <= vls
      fnum_h  = hls*vls
      fnum_hv = hls*vls*vls + πls
    elseif 0 <= vr+cr_hr
      fnum_h  = hrs*vrs
      fnum_hv = hrs*vrs*vrs + πrs
    else
      fnum_h  = hr*vr
      fnum_hv = hr*vr*vr + 0.5*g*hr*hr
    end

    SVector(fnum_h, fnum_hv)
end

"""
    (::KineticFlux)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater{T,1})

The kinetic relaxation solver / numerical flux, see
Perthame and Simeoni (2001) A kinetic scheme for the Saint-Venant system with a
source term, Calcolo 38(4), pp. 201-231.
"""
function (::KineticFlux)(uₗ::ShallowWaterVar1D, uᵣ::ShallowWaterVar1D, model::ShallowWater{T,1}) where T
    hₗ, vₗ = primitive_variables(uₗ, model)
    hᵣ, vᵣ = primitive_variables(uᵣ, model)
    @unpack g = model

    # integrals for \xi \geq 0
    fnum_h, fnum_hv = kinetic_integrals_ξ_geq_0(hₗ, vₗ, g)
    # use symmetry for \xi \leq 0
    val_h, val_hv   = kinetic_integrals_ξ_geq_0(hᵣ, -vᵣ, g)
    fnum_h  -= val_h
    fnum_hv += val_hv

    SVector(fnum_h, fnum_hv)
end

# \int_{\xi \geq 0} (\xi, xi^2)  M_l(\xi) \dif \xi  for (fnum_h, fnum_hv)
function kinetic_integrals_ξ_geq_0(h, v, g)
    val_h  = zero(h)
    val_hv = zero(h)
    sqrt_2gh = sqrt(2g*h)

    if h > 0 && v >= sqrt_2gh
        val_h += h*v

        val_hv += h*v*v + g*h*h/2
    elseif h > 0 && sqrt_2gh+v > 0 && sqrt_2gh-v > 0
        sqrt_2ghmv2 = sqrt(2g*h-v^2)
        atan_v_sqrt_2ghmv2 = atan(v/sqrt_2ghmv2)

        val_h += h * v / 2
        val_h += h * sqrt_2ghmv2 * 2 / (3π)
        val_h += v^2 * sqrt_2ghmv2 / (6g*π)
        val_h += h*v * atan_v_sqrt_2ghmv2 / π

        val_hv += ( h*v*v + g*h*h/2 ) / 2
        val_hv += h*v*sqrt_2ghmv2 * 13 / (12π)
        val_hv += v^3*sqrt_2ghmv2 / (12g*π)
        val_hv += ( h*v^2 + g*h^2/2 ) * atan_v_sqrt_2ghmv2 / π
    end

    val_h, val_hv
end

################################################################################

struct ShallowWaterRiemannSolution{T,T1} <: AbstractRiemannSolution
  prob::RiemannProblem{ShallowWater{T,1},ShallowWaterVar1D{T},T1}
  uₘ::ShallowWaterVar1D{T} # middle state
  σ₁⁻::T # slow speed of first family
  σ₁⁺::T # fast speed of first family
  σ₂⁻::T # slow speed of second family
  σ₂⁺::T # fast speed of second family

  function ShallowWaterRiemannSolution{T,T1}(prob::RiemannProblem{ShallowWater{T,1},ShallowWaterVar1D{T},T1}) where {T,T1}
    uₘ, σ₁⁻, σ₁⁺, σ₂⁻, σ₂⁺ = compute_state_and_speeds(prob.uₗ, prob.uᵣ, prob.model)
    new(prob, uₘ, σ₁⁻, σ₁⁺, σ₂⁻, σ₂⁺)
  end
end

function solve(prob::RiemannProblem{ShallowWater{T,1},ShallowWaterVar1D{T},T1}) where {T,T1}
  ShallowWaterRiemannSolution{T,T1}(prob)
end

function minmax_speeds(sol::ShallowWaterRiemannSolution)
  sol.σ₁⁻, sol.σ₂⁺
end

function (sol::ShallowWaterRiemannSolution)(ξ::Real)
  @unpack uₘ, σ₁⁻, σ₁⁺, σ₂⁻, σ₂⁺ = sol
  @unpack uₗ, uᵣ = sol.prob
  @unpack g = sol.prob.model

  if ξ <= σ₁⁻
    uₗ
  elseif ξ <= σ₁⁺
    rarefaction_wave_1(ξ, uₗ, g)
  elseif ξ <= σ₂⁻
    uₘ
  elseif ξ <= σ₂⁺
    rarefaction_wave_2(ξ, uᵣ, g)
  else # σ₂⁺ < ξ
    uᵣ
  end
end

function rarefaction_wave_1(ξ, uₗ::ShallowWaterVar1D{T}, g) where {T}
  hₗ, vₗ = primitive_variables(uₗ)
  h  = (vₗ + 2*sqrt(g*hₗ) -  ξ)^2 / 9g
  hv = (vₗ + 2*sqrt(g*hₗ) + 2ξ)*h / 3
  ShallowWaterVar1D{T}(h, hv)
end

function rarefaction_wave_2(ξ, uₗ::ShallowWaterVar1D{T}, g) where {T}
  hₗ, vₗ = primitive_variables(uₗ)
  h  = (vₗ - 2*sqrt(g*hₗ) -  ξ)^2 / 9g
  hv = (vₗ - 2*sqrt(g*hₗ) + 2ξ)*h / 3
  ShallowWaterVar1D{T}(h, hv)
end

function (sol::ShallowWaterRiemannSolution)(t::Real, x::Real)
  @unpack x₀, t₀ = sol.prob

  sol((x-x₀)/(t-t₀))
end

function compute_state_and_speeds(uₗ::ShallowWaterVar1D{T}, uᵣ::ShallowWaterVar1D{T}, model::ShallowWater) where {T}
  @unpack g = model
  hₗ, vₗ = primitive_variables(uₗ)
  hᵣ, vᵣ = primitive_variables(uᵣ)

  if v_S1_SWE(hᵣ,hₗ,vₗ,g) <= vᵣ && vᵣ <= v_R2_SWE(hᵣ,hₗ,vₗ,g)
    # Case I: wave1 = shock, wave2 = rarefaction
    if hₗ ≈ 0
      hₘ = hₗ
      vₘ = v_R2_SWE(hₘ,hᵣ,vᵣ,g)
      uₘ = ShallowWaterVar1D{T}(hₘ, hₘ*vₘ)
      σ₁⁻ = σ₁⁺ = vₘ - sqrt(g*hₗ)
      σ₂⁻ = vₘ + sqrt(g*hₘ)
      σ₂⁺ = vᵣ + sqrt(g*hᵣ)
    else
      hₘ = Roots.fzero(h -> vᵣ-v_R2_SWE(hᵣ,h,v_S1_SWE(h,hₗ,vₗ,g),g), hₗ, hᵣ)
      # this is better than vₘ = v_S1_SWE(hₘ,hₗ,vₗ) for hₗ≈0
      vₘ = v_R2_SWE(hₘ,hᵣ,vᵣ,g)
      uₘ = ShallowWaterVar1D{T}(hₘ, hₘ*vₘ)
      # this is better than σ₁⁻ = σ₁⁺ = vₗ - sqrt(hₘ+hₘ^2/hₗ) / sqrt(2) for hₗ≈0
      σ₁⁻ = σ₁⁺ = vₘ - sqrt(hₗ+hₗ^2/hₘ) * sqrt(g/2)
      σ₂⁻ = vₘ + sqrt(g*hₘ)
      σ₂⁺ = vᵣ + sqrt(g*hᵣ)
    end
  elseif v_R1_SWE(hᵣ,hₗ,vₗ,g) <= vᵣ && v_R2_SWE(hᵣ,hₗ,vₗ,g) <= vᵣ
    # Case II: wave1 = rarefaction, wave2 = rarefaction
    if 2*(sqrt(g*hᵣ)+sqrt(g*hₗ)) >= vᵣ-vₗ
      # no vacuum
      hₘ = (2*(sqrt(g*hᵣ)+sqrt(g*hₗ)) -(vᵣ-vₗ))^2 / 16g
      vₘ = v_R1_SWE(hₘ,hₗ,vₗ,g)
      uₘ = ShallowWaterVar1D{T}(hₘ,hₘ*vₘ)
      σ₁⁻ = vₗ - sqrt(g*hₗ)
      σ₁⁺ = vₘ - sqrt(g*hₘ)
      σ₂⁻ = vₘ + sqrt(g*hₘ)
      σ₂⁺ = vᵣ + sqrt(g*hᵣ)
    else
      # vacuum
      uₘ = ShallowWaterVar1D{T}(0,0)
      σ₁⁻ = vₗ - sqrt(g*hₗ)
      σ₁⁺ = vₗ + 2*sqrt(g*hₗ)
      σ₂⁻ = vᵣ - 2*sqrt(g*hᵣ)
      σ₂⁺ = vᵣ + sqrt(g*hᵣ)
    end
  elseif v_S2_SWE(hᵣ,hₗ,vₗ,g) <= vᵣ && vᵣ <= v_R1_SWE(hᵣ,hₗ,vₗ,g)
    # Case III: wave1 = rarefaction, wave2 = shock
    if hᵣ ≈ 0
      hₘ = hᵣ
      vₘ = v_R1_SWE(hₘ,hₗ,vₗ,g)
      uₘ = ShallowWaterVar1D{T}(hₘ,hₘ*vₘ)
      σ₁⁻ = vₗ - sqrt(g*hₗ)
      σ₁⁺ = vₘ - sqrt(g*hₘ)
      σ₂⁻ = σ₂⁺ = vₘ + sqrt(g*hᵣ)
    else
      hₘ = Roots.fzero(h->vᵣ-v_S2_SWE(hᵣ,h,v_R1_SWE(h,hₗ,vₗ,g),g), hᵣ, hₗ)
      vₘ = v_R1_SWE(hₘ,hₗ,vₗ,g)
      uₘ = ShallowWaterVar1D{T}(hₘ,hₘ*vₘ)
      σ₁⁻ = vₗ - sqrt(g*hₗ)
      σ₁⁺ = vₘ - sqrt(g*hₘ)
      σ₂⁻ = σ₂⁺ = vₘ + sqrt(hᵣ+hᵣ^2/hₘ) * sqrt(g/2)
    end
  else #vᵣ <= v_S1_SWE(hᵣ,hₗ,vₗ) && vᵣ <= v_S2_SWE(hᵣ,hₗ,vₗ)
    # Case IV: wave1 = shock, wave2 = shock
    # TODO: Chose some interval?
    hₘ = Roots.fzero(h->vᵣ-v_S2_SWE(hᵣ,h,v_S1_SWE(h,hₗ,vₗ,g),g), max(hₗ,hᵣ))
    vₘ = v_S1_SWE(hₘ,hₗ,vₗ,g)
    uₘ = ShallowWaterVar1D{T}(hₘ,hₘ*vₘ)
    σ₁⁻ = σ₁⁺ = vₗ - hₘ * sqrt(1/hₘ+1/hₗ) * sqrt(g/2)
    σ₂⁻ = σ₂⁺ = vₘ + hᵣ * sqrt(1/hᵣ+1/hₘ) * sqrt(g/2)
  end

  uₘ, σ₁⁻, σ₁⁺, σ₂⁻, σ₂⁺
end

# wave curves parameterised by the height h
v_R1_SWE(h, hₗ, vₗ, g) = vₗ - 2*(sqrt(g*h) - sqrt(g*hₗ))
v_R2_SWE(h, hₗ, vₗ, g) = vₗ + 2*(sqrt(g*h) - sqrt(g*hₗ))
v_S1_SWE(h, hₗ, vₗ, g) = vₗ - (h-hₗ) * sqrt( (1/h+1/hₗ)*g/2 )
v_S2_SWE(h, hₗ, vₗ, g) = vₗ + (h-hₗ) * sqrt( (1/h+1/hₗ)*g/2 )

@recipe function f(ξu::Tuple{Ξ,Vector{ShallowWaterVar1D{T}}}) where {Ξ,T}
  ξ, u = ξu

  h  = mappedarray(u->u.h, u)
  hv = mappedarray(u->u.hv, u)

  size --> (1000, 400)
  layout --> (1,2)
  legend --> false

  @series begin
    subplot := 1
    ylabel --> L"h"
    label  --> L"h"
    ξ, h
  end

  @series begin
    subplot := 2
    ylabel --> L"hv"
    label  --> L"hv"
    ξ, hv
  end
end

ext: jl


text:
using JuliaFormatter
using JuliaFormatter: DefaultStyle, YASStyle, Options, CONFIG_FILE_NAME
using CSTParser
using Test

fmt1(s; i = 4, m = 80, kwargs...) =
    JuliaFormatter.format_text(s; kwargs..., indent = i, margin = m)
fmt1(s, i, m; kwargs...) = fmt1(s; kwargs..., i = i, m = m)

# Verifies formatting the formatted text
# results in the same output
function fmt(s; i = 4, m = 80, kwargs...)
    s1 = fmt1(s; kwargs..., i = i, m = m)
    return fmt1(s1; kwargs..., i = i, m = m)
end
fmt(s, i, m; kwargs...) = fmt(s; kwargs..., i = i, m = m)

function run_pretty(text::String; style = DefaultStyle(), opts = Options())
    d = JuliaFormatter.Document(text)
    s = JuliaFormatter.State(d, opts)
    x = CSTParser.parse(text, true)
    t = JuliaFormatter.pretty(style, x, s)
    t
end
run_pretty(text::String, max_margin::Int) =
    run_pretty(text, opts = Options(max_margin = max_margin))

function run_nest(text::String; opts = Options(), style = DefaultStyle())
    d = JuliaFormatter.Document(text)
    s = JuliaFormatter.State(d, opts)
    x = CSTParser.parse(text, true)
    t = JuliaFormatter.pretty(style, x, s)
    JuliaFormatter.nest!(style, t, s)
    t, s
end
run_nest(text::String, max_margin::Int) =
    run_nest(text, opts = Options(max_margin = max_margin))

yasfmt1(s, i, m; kwargs...) = fmt1(s; kwargs..., i = i, m = m, style = YASStyle())
yasfmt(s, i, m; kwargs...) = fmt(s; kwargs..., i = i, m = m, style = YASStyle())

include("default_style.jl")
include("yas_style.jl")
include("issues.jl")
include("options.jl")
include("config.jl")


text:
### SimpleMuseProblem

struct SimpleMuseProblem{X,S,Gθ,Pθ,GZ} <: AbstractMuseProblem
    x :: X
    sample_x_z :: S
    ∇θ_logLike :: Gθ
    logLike_and_∇z_logLike :: GZ
    logPriorθ :: Pθ
end

@doc doc"""
    SimpleMuseProblem(x, sample_x_z, logLike, logPriorθ=(θ->0); ad=AD.ForwardDiffBackend())

Specify a MUSE problem by providing the simulation and posterior
evaluation code by-hand. The argument `x` should be the observed data.
The function `sample_x_z` should have signature:

```julia
function sample_x_z(rng::AbstractRNG, θ)
    # ...
    return (;x, z)
end
```
and return a joint sample of data `x` and latent space `z`. The
function `logLike` should have signature:

```julia
function logLike(x, z, θ)
    # return log likelihood
end
```

and return the likelihood $\log\mathcal{P}(x,z\,|\,\theta)$ for
your problem. The optional function `logPriorθ` should have signature:

```julia
function logPriorθ(θ)
    # return log prior
end
```

and should return the prior $\log\mathcal{P}(\theta)$ for your
problem. The `autodiff` parameter should be either
`MuseInference.ForwardDiffBackend()` or
`MuseInference.ZygoteBackend()`, specifying which library to use for
automatic differenation through `logLike`.

All variables `(x, z, θ)` can be any types which support basic
arithmetic.

# Example

```julia
# 512-dimensional noisy funnel
prob = SimpleMuseProblem(
    rand(512),
    function sample_x_z(rng, θ)
        z = rand(rng, MvNormal(zeros(512), exp(θ)*I))
        x = rand(rng, MvNormal(z, I))
        (;x, z)
    end,
    function logLike(x, z, θ)
        -(1//2) * (sum((x .- z).^2) + sum(z.^2) / exp(θ) + 512*θ)
    end,
    function logPrior(θ)
        -θ^2/(2*3^2)
    end;
    autodiff = MuseInference.ZygoteBackend()
)

# get solution
muse(prob, (θ=1,))
```
"""
function SimpleMuseProblem(x, sample_x_z, logLike, logPriorθ=(θ->0); autodiff::AD.AbstractBackend=AD.ForwardDiffBackend())
    SimpleMuseProblem(
        x,
        sample_x_z,
        (x,z,θ) -> first(AD.gradient(autodiff, θ -> logLike(x,z,θ), θ)),
        (x,z,θ) -> first.(AD.value_and_gradient(autodiff, z -> logLike(x,z,θ), z)),
        logPriorθ
    )
end

∇θ_logLike(prob::SimpleMuseProblem, x, z, θ) = prob.∇θ_logLike(x, z, θ)
logPriorθ(prob::SimpleMuseProblem, θ) = prob.logPriorθ(θ)
logLike_and_∇z_logLike(prob::SimpleMuseProblem, x, z, θ) = prob.logLike_and_∇z_logLike(x, z, θ)
sample_x_z(prob::SimpleMuseProblem, rng::AbstractRNG, θ) = prob.sample_x_z(rng, θ)


text:
#--------------------------------------------------------------------------------------------------
using Test
using Mahler

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

@test Mahler.mexpansion((x)->2,    5) == [2,0,0,0,0]
@test Mahler.mexpansion((x)->x,    5) == [0,1,0,0,0]
@test Mahler.mexpansion((x)->x+1,  5) == [1,1,0,0,0]
@test Mahler.mexpansion((x)->2*x+1,5) == [1,2,0,0,0]

@test Mahler.mval(  9,Mahler.mexpansion((x) ->    +1, 5)) == 1
@test Mahler.mval(  0,Mahler.mexpansion((x) -> 2*x+1, 5)) == 1
@test Mahler.mval(0.5,Mahler.mexpansion((x) -> 2*x+1, 5)) == 2
@test Mahler.mval(  1,Mahler.mexpansion((x) -> 2*x+1, 5)) == 3
@test Mahler.mval(5.5,Mahler.mexpansion((x) -> 2*x+1, 5)) == 12

@test Mahler.mval(1,Mahler.mexpansion(log1p, 10)) ≈ log1p(1)
@test Mahler.mval(2,Mahler.mexpansion(log1p, 10)) ≈ log1p(2)
@test Mahler.mval(9,Mahler.mexpansion(log1p, 10)) ≈ log1p(9)

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

max_stars_repo_name: ayarodionov/Mahler

max_issues_repo_name: ayarodionov/Mahler

max_forks_repo_name: ayarodionov/Mahler


text:
# A JSON-LD processing implementation

# JSON-LD seems pretty weird at first glance

# Needs a good understanding of IRI's too

# The main idea for JSON-LD is that IRI's used in keys can instead be mapped to keyword values using contexts to allow easier processing
# JSON-LD also allows type metadata to be specified with the @type property that can be applied to context keyword definitions and to document nodes

# Relative IRIs get interesting as they are dependent on the document from which they were fetched

# application/ld+json is the mime-type

# Needs an implementation for transforming from IRIs to URIs and so-forth
struct IRI
  identifier::String
end

# TODO perform transformation from IRI to URI
as_uri(id::IRI) = id.identifier

# If the IRI references a valid JSON-LD context spec then types can be generated from the spec easily enough
macro LdContextType(iri::String)
  # Tries to get the spec remotely from the provided IRI
end

macro LdType(iri::String, properties::Expr...)
  # Define an LdType using ordinary type syntax
  if properties.head != :type
    throw("LdTypes must be declared with the same syntax as Julia types")
  end
  # LdType properties are wrapped to be Nullable while processed to a regular type
end

# Need to implement three primary algorithms

function dereference_uri(uri)
  # TODO Having this be a blocking operation is likely to be an issue
  r = HTTP.request("GET", uri, [("Accept", "application/ld+json")])
  return JSON.parse(String(r.body))
end

# === CONTEXT PROCESSING ALGORITHM ===

# Based off of the default implementation
# TODO Something that does processing based on native Julia structs would be nicer
# Probably would be much more efficient and readable too

# Some sub-algorithms too, kind of sucks for type stability
# remote_contexts is there to avoid cyclicality problems
function add_local_context(active_context::Dict{String, Any}, local_context, remote_contexts::Set{String} = Set{String}())
  # Wrap the local context as an array
  return add_local_context(active_context, Any[local_context], remote_contexts)
end

function add_local_context(active_context::Dict{String, Any}, local_context::Array{Any, 1}, remote_contexts::Set{String} = Set{String}())
  # Making it a copy is probably good to keep logic clean
  # May be strictly unnecessary unless multi-threading is brought in
  result = copy(active_context)
  for added_context in local_context
    result = add_single_local_context(result, added_context, remote_contexts)
  end
end

# Local contexts must be strings, Dicts, or nothing
function add_single_local_context(active_context::Dict{String, Any}, local_context, remote_contexts::Set{String} = Set{String}())
  throw("Invalid local context: An invalid local context was detected")
end

function add_single_local_context(active_context::Dict{String, Any}, local_context::String, remote_contexts::Set{String} = Set{String}())
  uri = iri_to_uri(local_context)
  if uri in remote_contexts
    throw("Recursive context inclusion error")
  end
  # Required to be a dictionary
  context = dereference_uri(iri)
  if !(isa(context, Dict) && ("@context" in keys(context)))
    throw("Invalid remote context")
  end

  # TODO Hmm, I'm pretty sure I recurse using the value from @context, even though their algorithm is a little shady on that point
  return add_local_context(active_context, context["@context"], Set([remote_contexts..., uri]))
end

function add_single_local_context(active_context::Dict{String, Any}, local_context::Dict{String, Any}, remote_contexts::Set{String} = Set{String}())
  result = copy(active_context)
  # Check for an @base key
  if "@base" in keys(local_context)
    value = local_context["@base"]
  end
  if "@version" in keys(local_context) && local_context["@version"] != 1.1
    throw("Invalid @version value: The @version key was used in a context with an out of range value")
    # TODO handle processing mode behavior
  end
  if "@vocab" in keys(local_context)
  end
  # Process languages
  if "@language" in keys(local_context)
    value = local_context["@language"]
    if value == nothing
      delete!(result, "@language")
    elseif isa(value, String)
      local_context["@language"] = lowercase(value)
    else
      throw("Invalid default language")
    end
  end
  return result
end

# TODO is the copy operation really necessary? Seems like it should be a nop.
add_single_local_context(active_context::Dict{String, Any}, local_context::Void, remote_contexts::Set{String} = Set{String}()) = copy(active_context)

function iri_to_uri(iri)
  # TODO actually resolve the iri to a uri per RFC3987, ugh
  return iri
end


text:
const POINTS         = 0x0000
const LINES          = 0x0001
const LINE_LOOP      = 0x0002
const LINE_STRIP     = 0x0003
const TRIANGLES      = 0x0004
const TRIANGLE_STRIP = 0x0005
const TRIANGLE_FAN   = 0x0006
const QUADS          = 0x0007

function DrawElements(mode::Integer, count::Integer, type_::Integer, ibo::Integer)
	ccall( (:glDrawElements, lib), Void,
		(GLenum, GLsizei, GLenum, Ptr{GLvoid}),
		mode, count, type_, ibo)
end

function DrawElements(mode::Integer, count::Integer, type_::Integer, indices::Ptr)
	ccall( (:glDrawElements, lib), Void,
		(GLenum, GLsizei, GLenum, Ptr{GLvoid}),
		mode, count, type_, indices)
end

DrawElements{T}(mode::Integer, indices::Array{T}) =
	DrawElements(mode, length(indices), GLtype[T], convert(Ptr{T}, indices))


text:
function rober!(du, u, p, t)
    y1, y2, y3 = u
    k1, k2, k3 = p

    du[1] = -k1*y1 + k3*y2*y3
    du[2] = k1*y1 - k2*y2^2. - k3*y2*y3
    du[3] = k2*y2^2.

end

rober_model = DEsteady(
    func = rober!,
    u0 = [1.0,0.0,0.0],
    p = [1.0,0.0,0.0]
)

function bistable_ode!(du, u, p ,t)
    s1, s2 = u
    K1, K2, k1, k2, k3, k4, n1 , n2  = p
    du[1] = k1 / (1 + (s2/K2)^n1) - k3*s1
    du[2] = k2/  (1 + (s1/K1)^n2) - k4*s2
end

bistable_model = DEsteady(
    func = bistable_ode!,
    p = [1., 1., 20., 20., 5., 5.,  4., 4.],
    u0 =  [3., 1.]
)


text:
#!/usr/bin/env julia

using FromFile

@from "Helper.jl" import Helper

function main()
    fname = "example.txt"
    # fname = "input.txt"

    lines = Helper.read_lines(fname)
    for line in lines
        println(line)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
using Tricks: compat_hasmethod
# TODO use traits for compat_hasmethod checks

# Scalar parsers

# content input
"""
    noedparse(type::Type, node)
    nodeparse(type::Type, content)
    noedparse(::Type{Vector{T}}, nodes)
    nodeparse(::Type{Vector{T}}, contents)
Parses a node content with the specidied type
"""
nodeparse(type::Type{<:Number}, content::String) = parse(type, content)
nodeparse(type::Type{<:AbstractString}, content::String) = content
function nodeparse(type::Type, content::String)
    # TODO: better specialized method detection
    # https://julialang.slack.com/archives/C6A044SQH/p1578442480438100
    if hasmethod(type, Tuple{String}) &&  Core.Compiler.return_type(type, Tuple{Node})=== Union{}
        return type(content)
    elseif compat_hasmethod(convert, Tuple{String, type})
        return convert(type, content)
    elseif compat_hasmethod(parse, Tuple{type, String})
        return parse(type, content)
    else
        error("Could not parse a String as type $type")
    end
end

# elm input
nodeparse(type::Type{<:Union{Number, AbstractString}}, elm::Node) = nodeparse(type, elm.content)
function nodeparse(type::Type, elm::Node)
    content = elm.content
    if hasmethod(type, Tuple{String}) &&  Core.Compiler.return_type(type, Tuple{Node})=== Union{}
        return type(content)
    elseif compat_hasmethod(convert, Tuple{String, type})
        return convert(type, content)
    elseif compat_hasmethod(parse, Tuple{type, String})
        return parse(type, content)
    else
        # should be the last to avoid invoking generic methods
        return type(elm)
    end
end

# Vector parsers
# both elm and content input

# TODO slow for other than String|Number but compact
numorstr(::Nothing) = [Number, AbstractString]
@transform function nodeparse(type::Type{<:numorstr()}, elmsOrContents::Vector{Node})
    return nodeparse.(type, elmsOrContents)
end
@transform function nodeparse(type::Type{<:numorstr()}, elmsOrContents::Vector{String})
    return nodeparse.(type, elmsOrContents)
end
# TODO Faster code, but has redundancy:

function nodeparse(::Type{T}, contents::Vector{String}) where {T}
    elms_typed = Vector{T}(undef, length(contents))
    i = 1
    # TODO https://github.com/oxinabox/Tricks.jl/issues/2#issuecomment-630450764
    if hasmethod(type, Tuple{String}) &&  Core.Compiler.return_type(type, Tuple{Node})=== Union{}
        for content in contents
            elms_typed[i] = type(content)
            i+=1
        end
    elseif compat_hasmethod(convert, Tuple{String, type})
        for content in contents
            elms_typed[i] = convert(type, content)
            i+=1
        end
    elseif compat_hasmethod(parse, Tuple{type, String})
        for content in contents
            elms_typed[i] = parse(type, content)
            i+=1
        end
    else
        error("Could not parse a String as type $type")
    end
    return elms_typed
end

function nodeparse(type::Type{T}, elms::Vector{Node}) where {T}
    elms_typed = Vector{T}(undef, length(elms))
    i = 1
    # TODO https://github.com/oxinabox/Tricks.jl/issues/2#issuecomment-630450764
    if hasmethod(type, Tuple{String}) &&  Core.Compiler.return_type(type, Tuple{Node})=== Union{}
        for elm in elms
            elms_typed[i] = type(elm.content)
            i+=1
        end
    elseif compat_hasmethod(convert, Tuple{String, type})
        for elm in elms
            elms_typed[i] = convert(type, elm.content)
            i+=1
        end
    elseif compat_hasmethod(parse, Tuple{type, String})
        for elm in elms
            elms_typed[i] = parse(type, elm.content)
            i+=1
        end
    else
        # should be the last to avoid invoking generic methods
        for elm in elms
            elms_typed[i] = type(elm)
            i+=1
        end
    end
    return elms_typed
end


text:
"""
    boundingbox(polygon)

Compute axis-aligned bounding box for polygon.
"""
function boundingbox(polygon)
    sum(map(boundingbox, polygon))
end

function boundingbox(triangle::Triangle)
    matrix = hcat(triangle.a, triangle.b, triangle.c)
    AABB(
        minimum(matrix[1,:]),
        maximum(matrix[1,:]),
        minimum(matrix[2,:]),
        maximum(matrix[2,:]),
        minimum(matrix[3,:]),
        maximum(matrix[3,:]),
    )
end

function boundingbox(vertex::Vertex)
    AABB(
        vertex[1],  # xmin
        vertex[1],  # xmax
        vertex[2],  # ymin
        vertex[2],  # ymax
        vertex[3],  # zmin
        vertex[3],  # zmax
    )
end

function Base.:+(a::AABB, b::AABB)
    AABB(
        min(a.xmin, b.xmin),
        max(a.xmax, b.xmax),
        min(a.ymin, b.ymin),
        max(a.ymax, b.ymax),
        min(a.zmin, b.zmin),
        max(a.zmax, b.zmax),
    )
end

function Base.issubset(a::AABB, b::AABB)
    a.xmin >= b.xmin &&
    a.xmax <= b.xmax &&
    a.ymin >= b.ymin &&
    a.ymax <= b.ymax &&
    a.zmin >= b.zmin &&
    a.zmax <= b.zmax
end


text:
# SPDX-License-Identifier: X11
# 2020-11-22
# Sum and Product (300pt)

function main()
  S,P = parse.(Int, (split ∘ readline)())

  det::BigInt = big(S)^2 - 4P
  if det < 0
    println("No")
    return
  end

  rtdet = √det
  if !isinteger(rtdet)
    println("No")
    return
  end

  x,y = S + Int(rtdet), S - Int(rtdet)
  println(x & 1 == 0 && x > 0 &&
          y & 1 == 0 && y > 0
          ? "Yes" : "No")
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
__precompile__()

module Elements

using LinearAlgebra, Printf
using Distributed, SparseArrays
using ProgressMeter, Dates, StatsBase

using ..adiff, ..Materials
import ..Materials.getϕ

p = Int64(nworkers())
function setp(x)
  global p = Int64(x)
end
# continous elements
struct C2D{M}
  nodes::Vector{<:Integer}
  Nx::NTuple{M,Array{Float64,1}}
  Ny::NTuple{M,Array{Float64,1}}
  wgt::NTuple{M,Float64}
  V::Real
  mat::Materials.Material
end
struct C3D{M}
  nodes::Vector{<:Integer}
  Nx::NTuple{M,Array{Float64,1}}
  Ny::NTuple{M,Array{Float64,1}}
  Nz::NTuple{M,Array{Float64,1}}
  wgt::NTuple{M,Float64}
  V::Real
  mat::Materials.Material
end
struct CAS{M}
  nodes::Vector{<:Integer}
  N0::NTuple{M,Array{Float64,1}}
  Nx::NTuple{M,Array{Float64,1}}
  Ny::NTuple{M,Array{Float64,1}}
  X0::NTuple{M,Float64}
  wgt::NTuple{M,Float64}
  V::Real
  mat::Materials.Material
end
#  structural elements
struct Rod
  nodes::Vector{<:Integer}         # node id
  r0::Array{Float64}          #
  l0::Float64
  A::Float64                  # area
  mat::Materials.Material     # material properties
end
struct Beam
  nodes::Vector{<:Integer}
  r0::Array{Float64,1}
  L::Real
  t::Real
  w::Real
  lgwx::Array{Tuple{Float64,Float64},1}
  lgwy::Array{Tuple{Float64,Float64},1}
  mat::Materials.Material
end
Elems   = Union{Rod, Beam, C2D, C3D, CAS}
CElems  = Union{C2D, C3D, CAS}
CartEls = Union{C2D, C3D}
# structure for constraint eqs
struct ConstEq
  func
  iDoFs::Array{Int64}
  D::Type
end
ConstEq(func, iDoFs) = ConstEq(func, iDoFs, adiff.D2)
# constructors
function Rod(nodes, p0, A; mat=Materials.Hooke())

  r0  = p0[2]-p0[1]
  l0  = norm(r0)
  Rod(nodes.|>Int64, r0.|>Float64, l0|>Float64, A|>Float64, mat)
end
function Beam(nodes, p0, t, w; mat=Materials.Hooke(1, 0.3), Nx = 5, Ny = 3)

  lgwx = lgwt(Nx)
  lgwy = lgwt(Ny, a=-0.5, b=0.5)

  d0  = p0[2]-p0[1]
  L   = norm(d0)
  r0  = d0/L

  Beam(nodes, r0, L, Float64(t), Float64(w), lgwx, lgwy, mat)
end
function Tria(nodes::Vector{<:Integer},
              p0::Vector{Vector{T}} where T<:Real ;
              mat=Materials.Hooke())
  (Nx,Ny,wgt,A) = begin
    (x1,x2,x3) = (p0[1][1],p0[2][1],p0[3][1])
    (y1,y2,y3) = (p0[1][2],p0[2][2],p0[3][2])
    Delta      = x1*y2-x2*y1-x1*y3+x3*y1+x2*y3-x3*y2
    Nx         = [y2-y3,y3-y1,y1-y2]./Delta
    Ny         = [x3-x2,x1-x3,x2-x1]./Delta
    A          = abs(Delta)/2
    ((Nx,),(Ny,),(A,), A)
  end

  C2D(nodes,Nx,Ny,wgt,A,mat)
end
function Quad(nodes::Vector{<:Integer},
              p0::Vector{Vector{T}} where T<:Real;
              mat=Materials.Hooke())

  (A,Nx,Ny,wgt) = begin
    r        = [-1, 1]*0.577350269189626 # √3/3
    N(ξ,η)   = [(1-ξ)*(1-η),(1+ξ)*(1-η),(1+ξ)*(1+η),(1-ξ)*(1+η)]/4

    Nx  = Array{Array{Float64,1},2}(undef,2,2)
    Ny  = Array{Array{Float64,1},2}(undef,2,2)
    wgt = Array{Float64,2}(undef,2,2)
    A   = 0
    for (ii, ξ) in enumerate(r), (jj, η) in enumerate(r)
      N0  = N(adiff.D2([ξ,η])...)
      p   = sum([N0[ii]p0[ii] for ii in 1:4])
      J   = [p[ii].g[jj] for jj in 1:2, ii in 1:2]
      Nxy = J\hcat(adiff.grad.(N0)...)

      Nx[ii,jj]  = Nxy[1,:]
      Ny[ii,jj]  = Nxy[2,:]
      wgt[ii,jj] = det(J)

      A +=wgt[ii,jj]
    end
    (A,tuple(Nx...),tuple(Ny...),tuple(wgt...))
  end

  C2D(nodes,Nx,Ny,wgt,A,mat)
end
function Tet04(nodes::Vector{<:Integer},
               p0::Vector{Vector{T}} where T<:Real;
               mat=Materials.Hooke())
  (V, Nx, Ny, Nz) = begin
    A        = ones(4,4)
    A[1,2:4] = p0[1]
    A[2,2:4] = p0[2]
    A[3,2:4] = p0[3]
    A[4,2:4] = p0[4]
    C        = inv(A)
    V        = det(A)/6
    (V,(C[2,:],),(C[3,:],),(C[4,:],))
  end
  C3D(nodes,Nx,Ny,Nz,(1.0,),V,mat)
end
function Tet10(nodes::Vector{<:Integer},
               p0::Vector{Vector{T}} where T<:Real;
               mat=Materials.Hooke())

  (V, Nx, Ny, Nz) = begin

    La = 0.585410196624968
    Lb = 0.138196601125010
    ξ  = [La, Lb, Lb, Lb]
    A  = hcat([[1,p[1],p[2],p[3],p[1]^2,p[2]^2,p[3]^2,p[2]p[3],p[1]p[3],p[1]p[2]] for p in p0]...)
    C  = inv(A)
    V  = det(A[1:4,1:4])/6

    Nx,Ny,Nz = zeros(10,4),zeros(10,4),zeros(10,4)

    for ii in 1:4
      ξ        = circshift([La, Lb, Lb, Lb], ii)
      p        = p0[1]ξ[1]+p0[2]ξ[2]+p0[3]ξ[3]+p0[4]ξ[4]
      Nx[:,ii] = C[:,2]+2C[:,5]p[1]+C[:,9]p[3]+C[:,10]p[2]
      Ny[:,ii] = C[:,3]+2C[:,6]p[2]+C[:,8]p[3]+C[:,10]p[1]
      Nz[:,ii] = C[:,4]+2C[:,7]p[3]+C[:,8]p[2]+C[:, 8]p[1]
    end

    (V,
     (tuple([Nx[:,ii] for ii in 1:4]...),),
     (tuple([Ny[:,ii] for ii in 1:4]...),),
     (tuple([Nz[:,ii] for ii in 1:4]...),))
  end

  C3D(nodes,Nx,Ny,Nz,V,mat)
end
function Hex08(nodes::Vector{<:Integer},
               p0::Vector{Vector{T}} where T<:Real;
               mat=Materials.Hooke())
  (V, Nx, Ny, Nz, wgt) = begin
    r        = [-1, 1]*0.577350269189626 # √3/3
    # N(ξ,η,ζ) = [(1-ξ)*(1-η)*(1+ζ),(1-ξ)*(1-η)*(1-ζ),(1-ξ)*(1+η)*(1-ζ),(1-ξ)*(1+η)*(1+ζ),
    #            (1+ξ)*(1-η)*(1+ζ),(1+ξ)*(1-η)*(1-ζ),(1+ξ)*(1+η)*(1-ζ),(1+ξ)*(1+η)*(1+ζ)]/8
    N(ξ,η,ζ) = [(1-ξ)*(1-η)*(1-ζ),(1+ξ)*(1-η)*(1-ζ), (1+ξ)*(1+η)*(1-ζ),(1-ξ)*(1+η)*(1-ζ),
                (1-ξ)*(1-η)*(1+ζ),(1+ξ)*(1-η)*(1+ζ), (1+ξ)*(1+η)*(1+ζ),(1-ξ)*(1+η)*(1+ζ)]/8

    Nx  = Array{Array{Float64,1},3}(undef,2,2,2)
    Ny  = Array{Array{Float64,1},3}(undef,2,2,2)
    Nz  = Array{Array{Float64,1},3}(undef,2,2,2)
    wgt = Array{Float64,3}(undef,2,2,2)
    V   = 0
    for (ii, ξ) in enumerate(r), (jj, η) in enumerate(r), (kk, ζ) in enumerate(r)
      N0   = N(adiff.D2([ξ,η,ζ])...)
      p    = sum([N0[ii]p0[ii] for ii in 1:8])
      J    = [p[ii].g[jj] for jj in 1:3, ii in 1:3]
      Nxyz = J\hcat(adiff.grad.(N0)...)

      Nx[ii,jj,kk]  = Nxyz[1,:]
      Ny[ii,jj,kk]  = Nxyz[2,:]
      Nz[ii,jj,kk]  = Nxyz[3,:]
      wgt[ii,jj,kk] = det(J)

      V +=wgt[ii,jj,kk]
    end
    (V,tuple(Nx...),tuple(Ny...),tuple(Nz...),tuple(wgt...))
  end
  C3D(nodes,Nx,Ny,Nz,wgt,V,mat)
end
function ASTria(nodes::Vector{<:Integer},
                p0::Vector{Vector{T}} where T<:Real;
                mat=Materials.Hooke())
  (N,Nx,Ny,X0,wgt,A) = begin
    (x1, x2, x3) = (p0[1][1], p0[2][1], p0[3][1])
    (y1, y2, y3) = (p0[1][2], p0[2][2], p0[3][2])

    Delta = x1*y2-x2*y1-x1*y3+x3*y1+x2*y3-x3*y2
    N     = [1, 1, 1]/3
    Nx    = [y2-y3, y3-y1, y1-y2]/Delta
    Ny    = [x3-x2, x1-x3, x2-x1]/Delta
    A     = abs(Delta)/2
    X0    = (x1+x2+x3)/3
    wgt   = A*2π*X0
    ((N,),(Nx,),(Ny,),(X0,),(wgt,),A)
  end
  CAS(nodes,N,Nx,Ny,X0,wgt,A,mat)
end
function ASQuad(nodes::Vector{<:Integer},
                p0::Vector{Vector{T}} where T<:Real;
                mat=Materials.Hooke())
  (V,N0,Nx,Ny,X0,wgt) = begin
    r        = [-1, 1]*0.577350269189626 # √3/3
    N(ξ,η)   = [(1-ξ)*(1-η),(1+ξ)*(1-η),(1+ξ)*(1+η),(1-ξ)*(1+η)]/4

    N0  = Array{Array{Float64,1},2}(undef,2,2)
    Nx  = Array{Array{Float64,1},2}(undef,2,2)
    Ny  = Array{Array{Float64,1},2}(undef,2,2)
    X0  = Array{Float64,2}(undef,2,2)
    wgt = Array{Float64,2}(undef,2,2)
    V   = 0
    for (ii, ξ) in enumerate(r), (jj, η) in enumerate(r)
      Nij = N(adiff.D2([ξ,η])...)
      p   = sum([Nij[ii]p0[ii] for ii in 1:4])
      J   = [p[ii].g[jj] for jj in 1:2, ii in 1:2]
      Nxy = J\hcat(adiff.grad.(Nij)...)

      N0[ii,jj]  = adiff.val.(Nij)
      Nx[ii,jj]  = Nxy[1,:]
      Ny[ii,jj]  = Nxy[2,:]
      X0[ii,jj]  = p[1].v
      wgt[ii,jj] = det(J)*2π*p[1].v

      V +=wgt[ii,jj]
    end
    (V,tuple(N0...),tuple(Nx...),tuple(Ny...),tuple(X0...),tuple(wgt...))
  end
  CAS(nodes,N0,Nx,Ny,X0,wgt,V,mat)
end
# elastic energy evaluation functions for elements
function getϕ(elem::Rod,  u::Array{<:Real,2})

  l   = norm(elem.r0+u[:,2]-u[:,1])
  F11 = l/elem.l0
  elem.A*elem.l0*getϕ(F11, elem.mat)

end
function getϕ(elem::Beam, u::Array{<:Real,2})

  L, r0, t, w = elem.L, elem.r0, elem.t, elem.w
  T    = [r0[1] r0[2]; -r0[2] r0[1]]
  u0   = vcat(T*u[1:2,1], u[3,1], T*u[1:2,2], u[3,2])
  u0_x = (u0[4]-u0[1])/L

  ϕ  = 0
  for (r,wr) in elem.lgwx
    x, dx     = r*L, wr*L

    v0_x  = (-6x/L^2 + 6x^2/L^3)u0[2] +
    (1 - 4x/L + 3x^2/L^2)u0[3] +
    (6x/L^2 - 6x^2/L^3)u0[5] +
    (-2x/L + 3x^2/L^2)u0[6]

    v0_xx = (-6/L^2 + 12x/L^3)u0[2] +
    (-4/L + 6x/L^2)u0[3] +
    (6/L^2 - 12x/L^3)u0[5] +
    (-2/L + 6x/L^2)u0[6]

    for (s,ws) in elem.lgwy
      y, dy  = s*elem.t, ws*elem.t

      dV   = dx*dy*elem.w
      C11 = (1+u0_x-v0_xx*y)^2 + v0_x^2
      ϕ   += getϕ(C11, elem.mat)*dV
    end
  end
  return ϕ
end
function getϕ(elem::T where T<:CElems, u::Array{<:Real,2})
  M = length(elem.wgt)
  if isa(u[1], adiff.D2)
    ϕ = sum([begin
               F = getF(elem,u,ii)
               ϕ = getϕ(adiff.D2(getfield.(F,:v)),elem.mat)
               elem.wgt[ii]cross(ϕ,F)
             end  for ii in 1:M])
  else
    ϕ = sum([elem.wgt[ii]getϕ(getF(elem,u,ii), elem.mat) for ii in 1:M])
  end
end
function cross(ϕ, F)
  N = length(F)
  g = sum([ϕ.g[ii]F[ii].g for ii in 1:N])
  h = sum([0.5ϕ.h[ii,jj]*(F[ii].g*F[jj].g+F[jj].g*F[ii].g) for jj=1:N for ii=1:N])
  adiff.D2(ϕ.v, g, h)
end
# methods for evaluating def. gradient
function getF(elem::C3D, u::Array{N} where N, ii::Int64)

  Nx, Ny, Nz = elem.Nx[ii], elem.Ny[ii], elem.Nz[ii]

  [Nx⋅u[1:3:end] Ny⋅u[1:3:end] Nz⋅u[1:3:end];
   Nx⋅u[2:3:end] Ny⋅u[2:3:end] Nz⋅u[2:3:end];
   Nx⋅u[3:3:end] Ny⋅u[3:3:end] Nz⋅u[3:3:end] ] + I
end
function getF(elem::C2D, u::Array{N} where N, ii::Int64)

  Nx, Ny = elem.Nx[ii], elem.Ny[ii]
  my0    = zero(u[1])

  # [Nx⋅u[1:2:end] Ny⋅u[1:2:end] my0;
  #  Nx⋅u[2:2:end] Ny⋅u[2:2:end] my0;
  #  my0           my0           my0] + I
  [Nx⋅u[1:2:end] Ny⋅u[1:2:end];
   Nx⋅u[2:2:end] Ny⋅u[2:2:end]] + I
end
function getF(elem::CAS, u::Array{N} where N, ii::Int64)
  Nx,  Ny   = elem.Nx[ii],   elem.Ny[ii]
  N0,  X0   = elem.N0[ii],   elem.X0[ii]
  u0x, u0y  = Nx⋅u[1:2:end], Ny⋅u[1:2:end]
  v0x, v0y  = Nx⋅u[2:2:end], Ny⋅u[2:2:end]
  w0z       = N0⋅u[1:2:end]/X0
  my0       = zero(u[1])

  [u0x  u0y   my0;
   v0x  v0y   my0;
   my0  my0   w0z] + I
end
function getJ(F)
  if (length(F)==9)
    F[1]F[5]F[9]-F[1]F[6]F[8]-F[2]F[4]F[9]+F[2]F[6]F[7]+F[3]F[4]F[8]-F[3]F[5]F[7]
  else
    F[1]F[4]-F[2]F[3]
  end
end
getV(elem,u)     = sum([elem.wgt[ii]getJ(elem,u,ii) for ii in 1:length(elem.wgt)])
getJ(elem,u,ii)  = getJ(getF(elem, u, ii))
getJ(elem,u)     = getV(elem,u)/elem.V
getI3(elem,u,ii) = getJ(getF(elem, u, ii))^2
getI3(elem,u)    = sum([elem.wgt[ii]getI3(elem,u,ii) for ii in 1:length(elem.wgt)])/elem.V
# elastic energy evaluation functions for models (lists of elements)
function getϕ(elems::Array, u)

  nElems = length(elems)
  if p==1 || nElems<=p
    (Φ, r, C) = getϕ(elems, u, 1:nElems)
  else
    nDoFs  = length(u)

    Φ = zeros(nElems)
    r = zeros(nDoFs)
    C = spzeros(nDoFs, nDoFs)

    chunks = split(nElems, Elements.p)
    procs  = [@spawn getϕ(elems, u, chunk)  for chunk in chunks]

    for (ii,chunk) in enumerate(chunks)
      retval =   fetch(procs[ii])
      Φ      .+= retval[1]
      r      .+= retval[2]
      C      .+= retval[3]
    end
  end
  (Φ, r, C)
end
function getϕ(elems::Array, u, chunk)

  nDoFs  = length(u)
  nElems = length(elems)

  Φ = zeros(size(elems))
  r = zeros(nDoFs)
  C = spzeros(nDoFs, nDoFs)
  for ii ∈ chunk
    elem           =  elems[ii]
    nodes          =  elem.nodes
    iDoFs          =  LinearIndices(u)[:,nodes][:]
    ϕ              =  getϕ(elem, adiff.D2(u[:,nodes]))
    Φ[ii]          =  adiff.val(ϕ)
    r[iDoFs]       += adiff.grad(ϕ)
    C[iDoFs,iDoFs] += adiff.hess(ϕ)
  end
  (Φ, r, C)
end
function getϕ(eqns::Array{ConstEq}, u::Array{Float64}, λ::Array{Float64})

  nEqs   = length(eqns)
  if p==1 || nEqs<=p
    (veqs, reqs, Keqs) = getϕ(eqns, u, λ, 1:nEqs)
  else
    nDoFs  = length(u)
    veqs   = zeros(nEqs)
    reqs   = spzeros(nDoFs, nEqs)
    Keqs   = spzeros(nDoFs, nDoFs)
    chunks = split(nEqs, Elements.p)
    procs  = [@spawn getϕ(eqns, u, λ, chunk)  for chunk in chunks]

    for ii in 1:p
      retval  = fetch(procs[ii])
      veqs .+= retval[1]
      reqs .+= retval[2]
      Keqs .+= retval[3]
    end
  end
  (veqs, reqs, Keqs)
end
function getϕ(eqns::Array{ConstEq}, u::Array{Float64}, λ::Array{Float64}, chunk)

  nEqs  = length(eqns)
  nDoFs = length(u)

  veqs  = zeros(nEqs)
  reqs  = spzeros(nDoFs, nEqs)
  Keqs  = spzeros(nDoFs, nDoFs)

  for ii ∈ chunk
    eqn               =  eqns[ii]
    iDoFs             =  eqn.iDoFs
    ϕ                 =  eqn.func(eqn.D(u[iDoFs]))
    veqs[ii]          =  adiff.val(ϕ)
    reqs[iDoFs,ii]    =  adiff.grad(ϕ)
    Keqs[iDoFs,iDoFs] += λ[ii]adiff.hess(ϕ)
  end
  (veqs, reqs, Keqs)
end
function getinfo(elem::Elems, u::Array{<:Real,2}; info=:detF)
  M = length(elem.Nx)
  F = sum([getF(elem, u, ii) for ii in 1:M])/M
  Materials.getinfo(F, elem.mat, info=info)
end
getinfo(elems::Array, u; info=:detF) =  [getinfo(elem, u[:,elem.nodes], info=info) for elem in elems]
# solver
function solve(elems, u;
               N          = 11,
               LF         = range(1e-4, stop=1, length=N),
               eqns       = [],
               λ          = zeros(length(eqns)),
               ifree      = isnan.(u),
               fe         = zeros(size(u)),
               bprogress  = false,
               becho      = false,
               dTol       = 1e-5,
               dTolu      = dTol,
               dTole      = 1e2dTol,
               dNoise     = 1e-12,
               maxiter    = 11,
               bechoi     = false,
               bprogressi = false,
               ballus     = true,
               bpredict   = true,
               maxupdt    = NaN)

  N     = length(LF)
  t0    = Base.time_ns()
  beqns = length(eqns)>0
  if bprogress; p    = ProgressMeter.Progress(length(LF)); end
  if ballus;    allu = [];  end

  fnew  = copy(fe)
  icnst = .!ifree
  unew  = copy(u)
  unew[ifree] .= 0
  uold  = zeros(size(unew))
  λnew  = copy(λ)

  for (ii,LF) in enumerate(LF)
    unew[icnst] .= u[icnst]*LF
    fnew[ifree] .= fe[ifree]*LF
    lastu       = copy(unew)
    lastλ       = copy(λnew)
    T           = @elapsed (bfailed, normr, iter) =
    try
      solvestep!(elems, uold, unew, ifree,
                 eqns      = eqns,
                 λ         = λnew,
                 fe        = fnew,
                 dTole     = dTole,
                 dTolu     = dTolu,
                 dNoise    = dNoise,
                 maxiter   = maxiter,
                 bprogress = bprogressi,
                 becho     = bechoi,
                 bpredict  = bpredict,
                 maxupdt   = maxupdt)
    catch
      (true, Inf, 0)
    end
    if bfailed
      @printf("\n!! failed at LF: %.3f, with normr/dTol: %.3e\n", LF, normr/dTol)
      unew = lastu
      λnew = lastλ
      break
    else
      uold[:] .= unew[:]
      if ballus
        if beqns
          push!(allu, (copy(unew), copy(fnew), copy(λnew)))
        else
          push!(allu, (copy(unew), copy(fnew)))
        end
      end
      bprogress && ProgressMeter.next!(p)
      becho     && @printf("step %3i/%i, LF=%.3f, done in %2i iter, after %.2f sec.\n",
                           ii,N,LF,iter,T)
    end
    becho && flush(stdout)
  end
  becho && @printf("completed in %s\n",(Base.time_ns()-t0)÷1e9|>
                   Dates.Second|>Dates.CompoundPeriod|>Dates.canonicalize)
  becho && flush(stdout)

  ballus ? allu : unew
end
function solvestep!(elems, uold, unew, bfreeu;
                    fe        = zeros(length(unew)),
                    eqns      = [],
                    λ         = zeros(length(eqns)),
                    dTol      = 1e-5,
                    dTolu     = dTol,
                    dTole     = 1e2dTol,
                    dNoise    = 1e-12,
                    maxiter   = 11,
                    becho     = false,
                    bprogress = false,
                    bpredict  = true,
                    maxupdt   = NaN)

  if bprogress
    p = ProgressMeter.ProgressThresh(dTolu)
  end

  ifreeu    = findall(bfreeu[:])
  icnstu    = findall(.!bfreeu[:])

  nEqs      = length(eqns)
  nfreeu    = length(ifreeu)
  ncnstu    = length(icnstu)
  nDoFs     = nfreeu + nEqs
  iius      = 1:nfreeu
  iieqs     = nfreeu .+ (1:nEqs)

  bdone     = false
  bfailed   = false
  iter      = 0
  normupdt  = 0
  normre    = NaN
  oldupdt   = zeros(nDoFs)
  updt      = zeros(nDoFs)
  if nEqs != 0
    H       = spzeros(nDoFs,nDoFs)
  end

  # predictor step
  if bpredict
    deltat = @elapsed begin
      Δucnst    = unew[icnstu]-uold[icnstu]
      (Φ,fi,Kt) = getϕ(elems, uold)
      if nEqs == 0
        res              = fi[ifreeu]-fe[ifreeu]
        res[:]         .-= Kt[ifreeu,icnstu]*Δucnst
        updt[:]          = Kt[ifreeu,ifreeu]\res
        unew[ifreeu]    .= uold[ifreeu] .+ updt
        normupdt         = maximum(abs.(updt))
      else
        (vEqs,rEqs,KEqs) = getϕ(eqns, uold, λ)
        resu             = fi[ifreeu]-fe[ifreeu]-rEqs[ifreeu,:]*λ
        resu[:]        .-= (Kt[ifreeu,icnstu]-KEqs[ifreeu,icnstu])*Δucnst
        rese             = -vEqs
        res              = vcat(resu, rese)

        H[iius,iius]     = Kt[ifreeu,ifreeu]-KEqs[ifreeu,ifreeu]
        H[iius,iieqs]    = -rEqs[ifreeu,:]
        H[iieqs,iius]    = transpose(H[iius,iieqs])
        H[iieqs,iieqs]   = spdiagm(0=>dNoise*randn(nEqs))
        updt[:]          = qr(H)\res
        unew[ifreeu]    .= uold[ifreeu] .+ updt[iius]
        λ              .-= updt[iieqs]
        normupdt         = maximum(abs.(updt[iius]))
      end
    end
    becho && @printf("\npredictor step done in %.2f sec., ", deltat)
    becho && @printf("with normupdt: %.2e, starting corrector loop\n", normupdt); flush(stdout)
  else
    unew[ifreeu] .= uold[ifreeu]
  end
  # corrector loop
  while !bdone & !bfailed
    global normru
    oldupdt = copy(updt)
    tic     = Base.time_ns()
    (Φ,fi,Kt) = getϕ(elems, unew)

    if nEqs == 0
      res    = fe[ifreeu]-fi[ifreeu]
      norm0  = ncnstu > 0     ? norm(fi[icnstu])/ncnstu : 0
      normru = norm0  > dTolu ? norm(res)/nfreeu/norm0  : norm(res)/nfreeu
      bdone  = (normru<dTolu)
    else
      (vEqs,rEqs,KEqs) = getϕ(eqns, unew, λ)
      resu   = fi[ifreeu]-fe[ifreeu]-rEqs[ifreeu,:]*λ
      rese   = -vEqs
      res    = -vcat(resu, rese)
      norm0  = ncnstu > 0     ? norm(fi[icnstu])/ncnstu : 0
      normru = norm0  > dTolu ? norm(res)/nfreeu/norm0  : norm(res)/nfreeu
      normre = maximum(abs.(rese))
      bdone  = (normru<dTolu) && (normre<dTole)
    end

    if bdone
      fe[:]   = nEqs==0 ? fi[:] : fi[:]-rEqs*λ
    elseif iter < maxiter
      if nEqs == 0
        updt[:]         = Kt[ifreeu,ifreeu]\res
        normupdt        = maximum(abs.(updt))
        if !isnan(maxupdt)
          if normupdt > maxupdt
            updt      .*= (maxupdt/normupdt)
            normupdt    = maximum(abs.(updt))
          end
        end
        unew[ifreeu]  .+= updt
      else
        H[iius,iius]    = Kt[ifreeu,ifreeu]-KEqs[ifreeu,ifreeu]
        H[iius,iieqs]   = -rEqs[ifreeu,:]
        H[iieqs,iius]   = transpose(H[iius,iieqs])
        H[iieqs,iieqs]  = spdiagm(0=>dNoise*randn(nEqs))
        updt[:]         = qr(H)\res
        normupdt        = maximum(abs.(updt[iius]))
        if !isnan(maxupdt)
          if normupdt > maxupdt
            updt      .*= (maxupdt/normupdt)
            normupdt    = maximum(abs.(updt[iius]))
          end
        end
        unew[ifreeu]  .+= updt[iius]
        λ             .+= updt[iieqs]
      end
    else
      bfailed = true
    end
    bprogress && ProgressMeter.update!(p, normru)
    if becho
      if (bdone | bfailed)
        @printf("iter: %2i, norm0: %.2e, normru: %.2e, normre: %.2e, eltime: %.2f sec.\n",
                iter, norm0, normru, normre, Int64(Base.time_ns()-tic)/1e9)
      else
        @printf("iter: %2i, norm0: %.2e, normru: %.2e, normre: %.2e, normupdt: %.2e, α: %6.3f, eltime: %.2f sec.\n",
                iter, norm0, normru, normre, normupdt,
                oldupdt⋅updt/norm(updt)/norm(oldupdt), Int64(Base.time_ns()-tic)/1e9)
      end
      flush(stdout)
    end
    iter  += 1
  end

  (bfailed, normru, iter)
end
# helper functions
# find the Gauss-Legendre quadrature points and weights
function lgwt(N::Integer; a=0, b=1)

  N, N1, N2 = N-1, N, N+1
  xu   = range(-1, stop=1,length=N1)
  y    = cos.((2collect(0:N) .+ 1)*pi/(2N+2)) .+ (0.27/N1)*sin.(π*xu*N/N2)
  L    = zeros(N1,N2)
  dTol = 1e-16
  y0   = 2

  while maximum(abs.(y.-y0)) > dTol
    L[:,1] .= 1
    L[:,2] .= y
    for k = 2:N1
      L[:,k+1]=((2k-1)*y.*L[:,k] .- (k-1)*L[:,k-1])/k
    end
    global Lp = N2*(L[:,N1] .- y .* L[:,N2])./(1 .- y.^2)
    y0 = y
    y  = y0 .- L[:,N2]./Lp
  end

  x = (a.*(1 .- y) .+ b.* (1 .+ y))./2
  w = (b-a)./((1 .- y.^2).*Lp.^2).*(N2/N1)^2

  return [(x[ii], w[ii]) for ii ∈ 1:N1]
  # (x,w)
end
function renumber!(nodes, elems)

  unodes   = sort(unique(vcat([elem.nodes for elem in elems]...)))
  el_nodes = hcat([elem.nodes for elem in elems]...)
  shift    = maximum(unodes) + 1

  el_nodes .+= shift
  for (ii, nodeid) in enumerate(unodes)
    nodeid += shift
    el_nodes[el_nodes .== nodeid] .= ii
  end
  nodes = nodes[unodes]

  for (ii,elem) in enumerate(elems)
    elem.nodes[:] .= el_nodes[:,ii]
  end
  (nodes, elems)
end
function split(N::Int64, p::Int64)
  n    = Int64(floor(N/p))
  nEls = ones(Int64, p)*n
  nEls[1:N-p*n] .+= 1

  slice = [ range(sum(nEls[1:ii-1])+1, length=nEls[ii])
           for ii in 1:p]
end

end


text:
module VariationalOnlineStats

using Statistics
using Base.Iterators

export MedianOnline

struct MedianOnline{T}
    med::T   # median
    mad::T   # mean absolute deviation
    n::Int
end
(::Type{MedianOnline{T}})(x) where T = MedianOnline{T}(x, zero(T), 1)

MedianOnline(x::T) where {T<:AbstractFloat} = MedianOnline{T}(x)

function MedianOnline{T}(v::AbstractVector{<:Number}) where {T<:AbstractFloat}
    m = MedianOnline{T}(first(v))
    for x in drop(v, 1)
        m = median(m, x)
    end
    m
end
MedianOnline(v::AbstractVector{T}) where {T<:AbstractFloat} = MedianOnline{T}(v)

@inline function Statistics.median(m::MedianOnline{T}, x) where {T}
    f = 1/(m.n + 1)  # fraction of contribution to mad
    Δx = x - m.med
    aΔx = abs(Δx)
    mad = (1-f) * m.mad + f * aΔx
    Δv = mad/m.n   # separation of the closest point among n from an
                   # exponential distribution with the given mad
    med = m.med + sign(Δx) * min(aΔx, Δv)
    MedianOnline{T}(med, mad, m.n+1)
end

(::Type{T})(m::MedianOnline{T}) where {T<:Number} = m.med
(::Type{T})(m::MedianOnline{T}) where T = m.med
(::Type{T})(m::MedianOnline{<:Number}) where {T<:Number} = convert(T, m.med)

end # module


text:
module NCDataFrame

import NCDatasets: Dataset, defDim, defVar, varbyattrib, name
import DataFrames: DataFrame

export writenc, readnc

"""
writenc(df, "path/to/data.nc")

Create netcdf file from DataFrame
"""
function writenc(df::DataFrame, path)
    ds = Dataset(path, "c")
    cols = names(df);
    for col in cols
        data = df[:, col];
        #ty = eltype(data);
        len = length(data);
        defDim(ds, col, len);
        defVar(ds, col, data, (col,));
    end
    close(ds)
end

"""
df = readnc("path/to/data.nc")

Read netcdf file to DataFrame
"""
function readnc(path)
    ds = Dataset(path, "r")
    cols = map(name, varbyattrib(ds))
    df = DataFrame()
    for col in cols
        df[!,col] = ds[col][:]
    end
    return df
end

end # module


text:
"""
    imlink(imgs...; gridsize=(1,length(imgs)), dims=(1,2))

Show multiple images in a single window, linking higher-dimensional
axes to shared GUI control(s).
"""
function imlink(imgs...; gridsize=imlink_grid(imgs), dims=(1,2))
    zr, slicedata = roi(first(imgs), dims)
    gd = imshow_gui((200, 200), slicedata, gridsize)
    guidata = Vector{Any}(length(imgs))
    for (img, g, i) in zip(imgs, CartesianRange(gridsize), 1:length(imgs))
        if isa(img, AbstractArray)
            guidata[i] = imshow(gd["frame"][g], gd["canvas"][g], img, nothing, zr, slicedata)
        else
            guidata[i] = imshow(gd["frame"][g], gd["canvas"][g], img, zr, slicedata)
        end
    end
    gd["guidata"] = guidata
    showall(gd["window"])
    gd
end

imlink_grid(imgs) = (1,length(imgs))


text:
using Revise
using ADSeismic
using ADCME
using PyPlot
using MAT
using DelimitedFiles
using JLD2
using PyCall
matplotlib.use("Agg")
close("all")

ADCME.options.customop.verbose = false

output_dir = "data/elastic"
if !ispath(output_dir)
  mkpath(output_dir)
end

################### Generate synthetic data #####################
reset_default_graph()
mpi_init()

# model_name = "models/marmousi2-model-true-large.mat"
# model_name = "models/marmousi2-model-true.mat"
# model_name = "models/BP-model-true.mat"
model_name = "models/marmousi2-model-smooth.mat"

param = load_params(model_name, "MPIElastic")
compute_PML_Params!(param)

vp = matread(model_name)["vp"]
vs = matread(model_name)["vs"]
rho = matread(model_name)["rho"]
##
vs = matread("models/marmousi2-model-true.mat")["vs"]
rho = matread("models/marmousi2-model-true.mat")["rho"]

@info "loading data..."
@load "data/data$(mpi_rank())-$(mpi_size()).jld2" Rs

vp_ = Variable(vp)
var_size = size(vp_)

vp = mpi_bcast(vp_)
vs = constant(vs)
rho = constant(rho)

vp = extract_local_patch(param, vp, tag = 0)
vs = extract_local_patch(param, vs, deps = sum(vp), tag = 810000)
rho = extract_local_patch(param, rho, deps = sum(vs), tag = 820000)

λ, μ, ρ = compute_lame_parameters(vp, vs, rho)

# src = load_elastic_source(model_name, use_mpi=true, param=param)[4:4]
# receiver = load_elastic_receiver(model_name, use_mpi=true, param=param)[4:4]
# propagator = x->MPIElasticPropagatorSolver(param, x, ρ, λ, μ)
# [MPISimulatedObservation!(propagator(src[i]), receiver[i]) for i = 1:length(src)]

src = load_elastic_source(model_name, use_mpi=true, param=param)
receiver = load_elastic_receiver(model_name, use_mpi=true, param=param)

Rs_ = []
Vs_ = []

# for i = 1:length(src)
local_loss = constant(0.0)
for i = 1:4

  rcvv = Rs[i]
  src_ = src[i]
  receiver_ = receiver[i]

  propagator = MPIElasticPropagatorSolver(param, src_, ρ, λ, μ; tag_offset = i*100000, dep = length(Vs_)==0 ? nothing : sum(Vs_[end]))
  MPISimulatedObservation!(propagator, receiver_)

  push!(Rs_, receiver_.rcvv)
  push!(Vs_, propagator.vy)

  if !isnothing(rcvv)
      global local_loss += sum((receiver_.rcvv - rcvv )^2)
  end

end

local_loss += sum(Vs_[end]) * 1e-20

loss = mpi_sum(local_loss) * 1e7
g = gradients(loss, vp_)

py"""
from tensorflow.core.protobuf import rewriter_config_pb2
import tensorflow as tf
config = tf.ConfigProto(inter_op_parallelism_threads=1)
off = rewriter_config_pb2.RewriterConfig.OFF
config.graph_options.rewrite_options.memory_optimization  = off
"""
session_conf = py"config"
sess = Session(config = session_conf); init(sess)

@info "running loss"
@info run(sess, loss)
@info "running gradients"
@info run(sess, g)

# using ADOPT
# @info "starting..."
# @info "initial local loss = $(run(sess, local_loss)) at rank = $(mpi_rank())"

# function calculate_loss(x)
#     x = reshape(x, var_size)
#     L = run(sess, loss, vp_=>x)
#     L
# end

# iter = 0
# iter_result = Array{Float64, 2}[]
# function calculate_gradients(G, x)
#     x = reshape(x, var_size)
#     G[:] = run(sess, g, vp_=>x)
#     if iter % 1 == 0
#       clf()
#       pcolormesh(x', cmap="jet")
#       gca().invert_yaxis()
#       colorbar()
#       savefig("test_figures/x-$(iter).png")
#     end
#     global iter += 1
# end

# losses = Float64[]
# function step_callback(x)
#     @info "loss = $x"
#     push!(losses, x)
# end

# initial_x = run(sess, vp_)
# options = Options()
# result = ADOPT.mpi_optimize(calculate_loss, calculate_gradients, initial_x[:], LBFGS(), options; step_callback = step_callback)

# if mpi_rank()==0 && mpi_size()==1
#   ii = 1 # source number
#   u = run(sess, propagator(src[ii]).vx)
#   p = visualize_wavefield(u, param)
#   saveanim(p, "marmousi-wavefield.gif")
# end

if mpi_size()>1
    mpi_finalize()
end


text:
using Evolutionary: initial_population
# Implementation of an elementary cellular automata
# according to https://mathworld.wolfram.com/ElementaryCellularAutomaton.html
# uses random initialization and uses wolfram codes to specify the rule
# synchronous update of the 1D lattice

using Agents, Random
using CairoMakie
using InteractiveDynamics
using CSV
using Distributions
using Evolutionary

RADIUS = 3
RULE_LENGTH = 2^(2*RADIUS+1)
N_CELLS = 149

"""
The automaton living in a 1D space
"""
mutable struct Cell <: AbstractAgent
    id::Int
    pos::Dims{1}
    status::Bool # either 0 or 1, where 1 is 'alive'
end

"""
Takes the status of a neighborhood and returns the corresponding
index in the model rule (0,0,0) -> 1 (0,1,0) -> 3
"""
function configuration_index(cell_statuses)
    # takes the tuples and forms a string (0,1,1) -> "011"
    binary_code = string(UInt8.(cell_statuses)...)
    # turns the string binary code into the integer
    index = parse(Int, binary_code, base=2) + 1
    return index
end

"""
Given a cell checks its neighbors and
decides the next status of the cell based on the model rule
"""
function next_status(cell, model)
    neighbors = collect(nearby_agents(cell, model, model.radius))
    cell_statuses = (c.status for c in vcat(neighbors[1:model.radius], cell, neighbors[model.radius+1:end]))
    index = configuration_index(cell_statuses)
    return model.rule[index]
end

"""
Initializes the ABM

# Arguments
- `n_cells::Int`: total of cell automaton
- `wolfram_code::Int`: wolfram code for the rule
- `seed::Int`: random seed
"""
function build_model(rule; n_cells = 100, cell_p=0.5, radius=3, seed = 30)
    space = GridSpace((n_cells,); metric=:chebyshev)

    properties = Dict(:rule => rule, :radius => radius)
    model = ABM(
        Cell,
        space;
        properties,
        rng = MersenneTwister(seed))

    dist = Bernoulli(cell_p)
    for x in 1:n_cells
        cell = Cell(nextid(model), (x,), rand(dist))
        add_agent_pos!(cell, model)
    end

    return model
end

"""
Dummy update of the cells
"""
function cell_step!(_, _)

end

"""
Performs a synchronous update of all cells
"""
function ca_step!(model)
    new_statuses = fill(false, nagents(model))
    for agent in allagents(model)
        new_statuses[agent.id] = next_status(agent, model)
    end

    for id in allids(model)
        model[id].status = new_statuses[id]
    end
end

allsame(x) = all(y -> y == first(x), x)

function predominant_status(model)
    statuses = (c.status for c in allagents(model))
    return mean(statuses) >= 0.5
end

function score_scenario(model, target)
    statuses = (c.status for c in allagents(model))
    return all(y -> y == target, statuses)
end

function evaluate_rule(rule; n_cells=149, n_scenarios=100, max_steps=2*n_cells)
    cell_ps = rand(Uniform(0,1), n_scenarios)
    results = falses(n_scenarios)

    for i in 1:n_scenarios
        model = build_model(rule; n_cells=n_cells, cell_p=cell_ps[i], radius=RADIUS)
        target = predominant_status(model)
        run!(model, cell_step!, ca_step!, max_steps);
        results[i] = score_scenario(model, target)
    end

    return 1 - mean(results) # the lower the better
end

function random_rule(p)
    dist = Bernoulli(p)
    return rand(dist, RULE_LENGTH) .== 1
end

#println("running test eval")
#rintln(evaluate_rule(random_rule(0.5)))

println("Running optimization")
POP_SIZE = 100
MAX_ITERATIONS = 100

initial_rules = [random_rule(p) for p in rand(Uniform(0,1), POP_SIZE)]

ga = GA(
    populationSize = POP_SIZE,
    mutationRate = 0.15,
    epsilon = 0.1,
    selection = susinv,
    crossover = twopoint,
    mutation = flip
)

population = Evolutionary.initial_population(ga, initial_rules)

#result = Evolutionary.optimize(
#    evaluate_rule,
#    population,
#    ga,
#    Evolutionary.Options(
#        iterations=MAX_ITERATIONS, store_trace=true, show_trace=true, parallelization=:thread))

function plot_rule(rule; n_cells=149, cell_p=0.5, max_steps=2*n_cells)
    model = build_model(rule; n_cells=n_cells, cell_p=cell_p, radius=RADIUS)
    data, _ = run!(model, cell_step!, ca_step!, max_steps; adata=[:status])

    fig, ax = heatmap(data.id, data.step, data.status, colormap=:Blues_3)
    ax.yreversed = true
    return fig
end

rule_decimals(rule_hex) = [parse(Int, c, base= 16) for c in  rule_hex]
rule_from_decimals(rule_decimals) = BitArray(vcat([reverse(digits(d, base=2, pad=4)) for d in rule_decimals]...))
rule_from_hex(rule_hex) = (rule_from_decimals ∘ rule_decimals)(rule_hex)

maj_rule = BitArray(mean(digits(i, base=2, pad=7)) .> 0.5 for i in 0:127)

par_rule_hex = "0504058705000f77037755837bffb77f"
par_rule = rule_from_hex(par_rule_hex)

#plot_rule(Evolutionary.minimizer(result))

# Initialize model
#model = build_model(random_rule(0.5); n_cells=149, cell_p=0.5, radius=3)

# Runs the model and collects data
#data, _ = run!(model, cell_step!, ca_step!, 100; adata=[:status]);

# The data contains the step, id/position and status of the cell (1/0)
#data

#CSV.write("ca_data.csv", data);

# Lets plot the time evolution in the y axis
#heatmap(data.id, data.step, data.status, colormap=:Blues_3)


text:
function icp(c::Point3D, ċ::FreeVector3D, ω::Number)
    @framecheck c.frame ċ.frame
    Point3D(c.frame, horizontal_projection((c + ċ / ω).v)..., 0)
end

function icp(state::MechanismState, ω::Number)
    c = center_of_mass(state)
    h = momentum(state)
    ċ = FreeVector3D(h.frame, linear(momentum(state)) / mass(state.mechanism))
    icp(c, ċ, ω)
end

function icp(state::MechanismState)
    c = center_of_mass(state)
    z = c.v[3]
    gz = norm(state.mechanism.gravitational_acceleration)
    ω = sqrt(abs(gz / z))
    icp(state, ω)
end


text:
using Graphene
using Test
using GeometricalPredicates: Point2D
using Statistics: mean

test_gatom = GAtom(1, 2., 3., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
test_gbond = GBond(11, 22., 33., Set([44, 55, 66, 77]), "GBOND_Signature", 88, "dataset_99")
test_gpolygon = GPolygon(111, 222., 333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 888, "dataset_999", 6)
test_gdefect = GDefect(1111, 2222., 3333., Set([4444, 5555, 6666, 7777]), "GDEFECT_Signature", 8888, "dataset_9999", 14, Set([1, 11, 111]), "Divacancy")

test_simple_gatom = GAtom(1, 2., 3., Set([]), "", 0, "dataset")
test_simple_gbond = GBond(11, 22., 33., Set([]), "", 0, "dataset")
test_simple_gpolygon = GPolygon(111, 222., 333., Set([]), "", 0, "dataset", 6)

@testset "GAtom" begin
    @test isa(test_gatom, GAtom)
    @test get_id(test_gatom) == 1
    @test get_x(test_gatom) == 2.
    @test get_y(test_gatom) == 3.
    @test get_relatives(test_gatom) == Set([4, 5, 6, 7])
    @test get_signature(test_gatom) == "GATOM_Signature"
    @test get_frame(test_gatom) == 8
    @test get_dataset(test_gatom) == "dataset_9"
    @test get_noa(test_gatom) == 1
    @test get_type(test_gatom) == "Atom"
    @test isatom(test_gatom)
    @test get_members(test_gatom) == Set([1])

    @test GAtom(1, 2., 3., 0, "dataset") == test_simple_gatom
    @test GAtom(1, 2., 3., 0) == test_simple_gatom
    @test GAtom(1, 2., 3., Set([])) == test_simple_gatom
    @test GAtom(1, 2., 3.) == test_simple_gatom
    @test GAtom(1, Point2D(2., 3.)) == test_simple_gatom
end

@testset "GBond" begin
    @test isa(test_gbond, GBond)
    @test get_id(test_gbond) == 11
    @test get_x(test_gbond) == 22.
    @test get_y(test_gbond) == 33.
    @test get_relatives(test_gbond) == Set([44, 55, 66, 77])
    @test get_signature(test_gbond) == "GBOND_Signature"
    @test get_frame(test_gbond) == 88
    @test get_dataset(test_gbond) == "dataset_99"
    @test get_noa(test_gbond) == 2
    @test get_type(test_gbond) == "Bond"
    @test isbond(test_gbond)
    @test get_members(test_gbond) == Set([11])

    @test GBond(11, 22., 33., 0, "dataset") == test_simple_gbond
    @test GBond(11, 22., 33., 0) == test_simple_gbond
    @test GBond(11, 22., 33., Set([])) == test_simple_gbond
    @test GBond(11, 22., 33.) == test_simple_gbond
end

@testset "GPolygon" begin
    @test isa(test_gpolygon, GPolygon)
    @test get_id(test_gpolygon) == 111
    @test get_x(test_gpolygon) == 222.
    @test get_y(test_gpolygon) == 333.
    @test get_relatives(test_gpolygon) == Set([444, 555, 666, 777])
    @test get_signature(test_gpolygon) == "GPOLYGON_Signature"
    @test get_frame(test_gpolygon) == 888
    @test get_dataset(test_gpolygon) == "dataset_999"
    @test get_noa(test_gpolygon) == 6
    @test get_type(test_gpolygon) == "Polygon"
    @test ispolygon(test_gpolygon)
    @test get_members(test_gpolygon) == Set([111])

    @test GPolygon(111, 222., 333., 0, "dataset", 6) == test_simple_gpolygon
    @test GPolygon(111, 222., 333., 0, 6) == test_simple_gpolygon
    @test GPolygon(111, 222., 333., 6) == test_simple_gpolygon
end

@testset "GDefect" begin
    @test isa(test_gdefect, GDefect)
    @test get_id(test_gdefect) == 1111
    @test get_x(test_gdefect) == 2222.
    @test get_y(test_gdefect) == 3333.
    @test get_relatives(test_gdefect) == Set([4444, 5555, 6666, 7777])
    @test get_signature(test_gdefect) == "GDEFECT_Signature"
    @test get_frame(test_gdefect) == 8888
    @test get_dataset(test_gdefect) == "dataset_9999"
    @test get_noa(test_gdefect) == 14
    @test get_members(test_gdefect) == Set([1, 11, 111])
    @test get_type(test_gdefect) == "Divacancy"
end

test_gatom1 = GAtom(1, 12., 13., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
test_gatom2 = GAtom(2, 22., 23., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
test_gbond1 = GBond(11, 122., 133., Set([44, 55, 66, 77]), "GBOND_Signature", 8, "dataset_9")
test_gbond2 = GBond(22, 222., 233., Set([44, 55, 66, 77]), "GBOND_Signature", 8, "dataset_9")
test_gpolygon1 = GPolygon(111, 1222., 1333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 8, "dataset_9", 5)
test_gpolygon2 = GPolygon(222, 2222., 2333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 8, "dataset_9", 7)
test_g_for_sort = [test_gatom1, test_gatom2, test_gpolygon2, test_gbond1, test_gpolygon1, test_gbond2]

@testset "isless and sort" begin
    @test isless(test_gatom1, test_gpolygon1)
    @test isless(test_gbond2, test_gpolygon2)
    @test isless(test_gpolygon1, test_gpolygon2)
end

link_relatives!(test_gatom1, test_gatom2)
@test get_relatives(test_gatom1) == Set([4, 5, 6, 7, 2])
@test get_relatives(test_gatom2) == Set([4, 5, 6, 7, 1])
link_relatives!(test_gbond1, test_gbond2)
@test get_relatives(test_gbond1) == Set([44, 55, 66, 77, 22])
@test get_relatives(test_gbond2) == Set([44, 55, 66, 77, 11])
link_relatives!(test_gpolygon1, test_gpolygon2)
@test get_relatives(test_gpolygon1) == Set([444, 555, 666, 777, 222])
@test get_relatives(test_gpolygon2) == Set([444, 555, 666, 777, 111])

link_relatives!(test_gatom1, test_gbond1)
@test get_relatives(test_gatom1) == Set([4, 5, 6, 7, 2, 11])
@test get_relatives(test_gbond1) == Set([44, 55, 66, 77, 22, 1])
link_relatives!(test_gbond2, test_gatom2)
@test get_relatives(test_gatom2) == Set([4, 5, 6, 7, 1, 22])
@test get_relatives(test_gbond2) == Set([44, 55, 66, 77, 11, 2])
link_relatives!(test_gatom1, test_gpolygon1)
@test get_relatives(test_gatom1) == Set([4, 5, 6, 7, 2, 11, 111])
@test get_relatives(test_gpolygon1) == Set([444, 555, 666, 777, 222, 1])
link_relatives!(test_gpolygon2, test_gatom2)
@test get_relatives(test_gatom2) == Set([4, 5, 6, 7, 1, 22, 222])
@test get_relatives(test_gpolygon2) == Set([444, 555, 666, 777, 111, 2])
link_relatives!(test_gbond1, test_gpolygon1)
@test get_relatives(test_gbond1) == Set([44, 55, 66, 77, 22, 1, 111])
@test get_relatives(test_gpolygon1) == Set([444, 555, 666, 777, 222, 1, 11])
link_relatives!(test_gpolygon2, test_gbond2)
@test get_relatives(test_gbond2) == Set([44, 55, 66, 77, 11, 2, 222])
@test get_relatives(test_gpolygon2) == Set([444, 555, 666, 777, 111, 2, 22])

test_gatom1 = GAtom(1, 12., 13., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
test_gbond1 = GBond(11, 122., 133., Set([44, 55, 66, 77]), "GBOND_Signature", 8, "dataset_9")
test_gpolygon1 = GPolygon(111, 1222., 1333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 8, "dataset_9", 6)
test_gpolygon2 = GPolygon(222, 2222., 2333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 8, "dataset_9", 6)

link_relatives!(test_gpolygon2, [test_gatom1, test_gbond1, test_gpolygon1])
@test get_relatives(test_gatom1) == Set([4, 5, 6, 7, 222])
@test get_relatives(test_gbond1) == Set([44, 55, 66, 77, 222])
@test get_relatives(test_gpolygon1) == Set([444, 555, 666, 777, 222])
@test get_relatives(test_gpolygon2) == Set([444, 555, 666, 777, 1, 11, 111])

test_gatom1 = GAtom(1, 12., 13., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
test_gbond1 = GBond(11, 122., 133., Set([44, 55, 66, 77]), "GBOND_Signature", 8, "dataset_9")
test_gbond2 = GBond(22, 222., 233., Set([44, 55, 66, 77]), "GBOND_Signature", 8, "dataset_9")
test_gpolygon1 = GPolygon(111, 1222., 1333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 8, "dataset_9", 6)

link_relatives!(test_gbond2, [test_gatom1, test_gbond1, test_gpolygon1])
@test get_relatives(test_gatom1) == Set([4, 5, 6, 7, 22])
@test get_relatives(test_gbond1) == Set([44, 55, 66, 77, 22])
@test get_relatives(test_gpolygon1) == Set([444, 555, 666, 777, 22])

test_gatom1 = GAtom(1, 12., 13., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
test_gpolygon1 = GPolygon(111, 1222., 1333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 8, "dataset_9", 6)
test_gpolygon2 = GPolygon(222, 2222., 2333., Set([444, 555, 666, 777]), "GPOLYGON_Signature", 8, "dataset_9", 6)

link_relatives!(test_gatom1, [test_gpolygon1, test_gpolygon2])
@test get_relatives(test_gatom1) == Set([4, 5, 6, 7, 111, 222])
@test get_relatives(test_gpolygon1) == Set([444, 555, 666, 777, 1])
@test get_relatives(test_gpolygon2) == Set([444, 555, 666, 777, 1])

@test make_gatom(collect(test_indexed_atoms)[5]) == GAtom(5, 15., 26.1, Set([]), "", 0, "dataset")
gatoms_vector = make_gatom.(collect(test_indexed_atoms))
@test length(gatoms_vector) == length(test_indexed_atoms)
@test gatoms_vector[5] == GAtom(5, 15., 26.1, Set([]), "", 0, "dataset")

test_gatom1 = GAtom(1, 12., 13., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
test_gatom2 = GAtom(2, 22., 23., Set([4, 5, 6, 7]), "GATOM_Signature", 8, "dataset_9")
@test make_gbond([test_gatom1, test_gatom2], 12) == GBond(12, 17., 18., Set([1, 2]), "", 0, "dataset")
@test make_gbond!([test_gatom1, test_gatom2], 12) == GBond(12, 17., 18., Set([1, 2]), "", 0, "dataset")
@test test_gatom1 == GAtom(1, 12., 13., Set([4, 5, 6, 7, 12]), "GATOM_Signature", 8, "dataset_9")
@test test_gatom2 == GAtom(2, 22., 23., Set([4, 5, 6, 7, 12]), "GATOM_Signature", 8, "dataset_9")

test_bondmatrix = make_bondmatrix(test_bonds)
@test count(!iszero, test_bondmatrix) == 30
@test sort(unique([x for x in test_bondmatrix if !iszero(x)])) == collect(1:15)

@test make_bondmatrix([(1,2)]) == [0 1; 1 0]
@test make_bondmatrix([(1,2), (1,4), (3,4)]) == [0 1 0 2;
                                                 1 0 0 0;
                                                 0 0 0 3;
                                                 2 0 3 0]

test_gpolygon1 = make_gpolygon([test_gatom1, test_gatom2], [test_gbond1, test_gbond2], 111)
test_gpolygon2 = make_gpolygon!([test_gatom1, test_gatom2], [test_gbond1, test_gbond2], 111)
@test get_x(test_gpolygon1) == get_x(test_gpolygon2) == mean(get_x.([test_gatom1, test_gatom2]))
@test get_y(test_gpolygon1) == get_y(test_gpolygon2) == mean(get_y.([test_gatom1, test_gatom2]))
@test get_id(test_gpolygon1) == get_id(test_gpolygon2)
@test get_relatives(test_gpolygon1) == Set([])
@test get_relatives(test_gpolygon2) == Set([1, 11, 2, 22])

test_gatom1 = GAtom(1, 12., 13., Set([66, 11, 111, 222]), "", 8, "dataset")
test_gatom2 = GAtom(2, 22., 23., Set([11, 22, 111, 222]), "", 8, "dataset")
test_gatom3 = GAtom(3, 32., 33., Set([22, 33, 111, 222]), "", 8, "dataset")
test_gatom4 = GAtom(4, 42., 43., Set([33, 44, 111, 222]), "", 8, "dataset")
test_gatom5 = GAtom(5, 52., 53., Set([44, 55, 111, 222]), "", 8, "dataset")
test_gatom6 = GAtom(6, 62., 63., Set([55, 66, 111, 222]), "", 8, "dataset")
test_gbond1 = GBond(11, 122., 133., Set([1, 2, 111, 222]), "", 8, "dataset")
test_gbond2 = GBond(22, 222., 233., Set([2, 3, 111, 222]), "", 8, "dataset")
test_gbond3 = GBond(33, 322., 333., Set([3, 4, 111, 222]), "", 8, "dataset")
test_gbond4 = GBond(44, 422., 433., Set([4, 5, 111, 222]), "", 8, "dataset")
test_gbond5 = GBond(55, 522., 533., Set([5, 6, 111, 222]), "", 8, "dataset")
test_gbond6 = GBond(66, 622., 633., Set([6, 1, 111, 222]), "", 8, "dataset")
test_gpolygon1 = GPolygon(111, 1222., 1333., Set([1, 2, 3, 4, 5, 6, 11, 22, 33, 44, 55, 66, 111]), "", 8, "dataset", 6)
test_gpolygon2 = GPolygon(222, 2222., 2333., Set([1, 2, 3, 4, 5, 6, 11, 22, 33, 44, 55, 66, 222]), "", 8, "dataset", 6)

test_graphene = vcat([test_gatom1, test_gatom2, test_gatom3, test_gatom4, test_gatom5, test_gatom6],
                      [test_gbond1, test_gbond2, test_gbond3, test_gbond4, test_gbond5, test_gbond6],
                      [test_gpolygon1, test_gpolygon2])
test_noa_dict = Dict([Pair(x._id, get_noa(x)) for x in test_graphene])
map(x -> make_signature!(x, test_noa_dict), test_graphene)

@test test_gatom1 == GAtom(1, 12., 13., Set([66, 11, 111, 222]), "6-2|", 8, "dataset")
@test test_gatom2 == GAtom(2, 22., 23., Set([11, 22, 111, 222]), "6-2|", 8, "dataset")
@test test_gatom3 == GAtom(3, 32., 33., Set([22, 33, 111, 222]), "6-2|", 8, "dataset")
@test test_gatom4 == GAtom(4, 42., 43., Set([33, 44, 111, 222]), "6-2|", 8, "dataset")
@test test_gatom5 == GAtom(5, 52., 53., Set([44, 55, 111, 222]), "6-2|", 8, "dataset")
@test test_gatom6 == GAtom(6, 62., 63., Set([55, 66, 111, 222]), "6-2|", 8, "dataset")
@test test_gbond1 == GBond(11, 122., 133., Set([1, 2, 111, 222]), "6-2|", 8, "dataset")
@test test_gbond2 == GBond(22, 222., 233., Set([2, 3, 111, 222]), "6-2|", 8, "dataset")
@test test_gbond3 == GBond(33, 322., 333., Set([3, 4, 111, 222]), "6-2|", 8, "dataset")
@test test_gbond4 == GBond(44, 422., 433., Set([4, 5, 111, 222]), "6-2|", 8, "dataset")
@test test_gbond5 == GBond(55, 522., 533., Set([5, 6, 111, 222]), "6-2|", 8, "dataset")
@test test_gbond6 == GBond(66, 622., 633., Set([6, 1, 111, 222]), "6-2|", 8, "dataset")
@test test_gpolygon1 == GPolygon(111, 1222., 1333., Set([1, 2, 3, 4, 5, 6, 11, 22, 33, 44, 55, 66, 111]), "6-1|", 8, "dataset", 6)
@test test_gpolygon2 == GPolygon(222, 2222., 2333., Set([1, 2, 3, 4, 5, 6, 11, 22, 33, 44, 55, 66, 222]), "6-1|", 8, "dataset", 6)

@testset "Make Graphene" begin
    test_graphene1 = make_graphene(test_xy, max_bondlength=8, frame=12, dataset="test_dataset")
    test_graphene2 = make_graphene(test_xy, max_bondlength=12, frame=12, dataset="test_dataset")
    test_graphene3 = make_graphene(test_xy[:, 2], max_bondlength=8, frame=12, dataset="test_dataset")
    test_graphene4 = make_graphene(test_xy[:, 2:4], max_bondlength=12, frame=12, dataset="test_dataset")
    test_graphene5 = make_graphene(test_xy[:, [2,4]], max_bondlength=8, frame=12, dataset="test_dataset")
    test_graphene6 = make_graphene(test_xy[:, [2,4]], max_bondlength=12, frame=12, dataset="test_dataset")
    @test length(test_graphene1) == 15
    @test length(test_graphene2) == 33
    @test length(test_graphene3) == 1
    @test length(test_graphene4) == 5
    @test length(test_graphene5) == 2
    @test length(test_graphene6) == 4

    @test get_id.(test_graphene1) == collect(1:length(test_graphene1))
    @test get_id.(test_graphene2) == collect(1:length(test_graphene2))
end

test_defect_image = import_stack("test_files/test_defect_6_bw.h5")
test_graphene = make_graphene(centroid2xy(make_centroids(test_defect_image)), max_bondlength=10)

test_flower = first(find_defect(test_graphene, "Flower"))
test_divacancy = first(find_defect(test_graphene, "V2(585)"))
test_divacancy14 = first(find_defect(test_graphene, "V2(14)"))
test_butterfly = first(find_defect(test_graphene, "Butterfly"))
test_stonewales = first(find_defect(test_graphene, "Stone-Wales"))
test_5775 = first(find_defect(test_graphene, "5775"))
test_defects = find_defect(test_graphene, "Flower", "Divacancy", "Butterfly", "Stone-Wales", "5775")
test_defects[1]._id = length(test_graphene) + 1
test_defects[2]._id = length(test_graphene) + 1
test_defects[3]._id = length(test_graphene) + 1
test_defects[4]._id = length(test_graphene) + 1
test_defects[5]._id = length(test_graphene) + 1
test_defects[6]._id = length(test_graphene) + 1
@test test_defects[1] == test_flower
@test test_defects[2] == test_divacancy
@test test_defects[3] == test_divacancy14
@test test_defects[4] == test_butterfly
@test test_defects[5] == test_stonewales
@test test_defects[6] == test_5775

@testset "Stone-Wales" begin
    @test get_id(test_stonewales) == length(test_graphene) + 1
    @test get_frame(test_stonewales) == 1
    @test get_dataset(test_stonewales) == "dataset"
    @test get_noa(test_stonewales) == 16
    @test length(get_relatives(test_stonewales)) == 10
    @test length(get_members(test_stonewales)) == 39
    @test get_type(test_stonewales) == "Stone-Wales"
end

@testset "5775" begin
    @test get_id(test_5775) == length(test_graphene) + 1
    @test get_frame(test_5775) == 1
    @test get_dataset(test_5775) == "dataset"
    @test get_noa(test_5775) == 18
    @test length(get_relatives(test_5775)) == 12
    @test length(get_members(test_5775)) == 43
    @test get_type(test_5775) == "5775"
end

@testset "Divacancy" begin
    @test get_id(test_divacancy) == length(test_graphene) + 1
    @test get_frame(test_divacancy) == 1
    @test get_dataset(test_divacancy) == "dataset"
    @test get_noa(test_divacancy) == 14
    @test length(get_relatives(test_divacancy)) == 10
    @test length(get_members(test_divacancy)) == 33
    @test get_type(test_divacancy) == "Divacancy"
end

@testset "Divacancy-14" begin
    @test get_id(test_divacancy14) == length(test_graphene) + 1
    @test get_frame(test_divacancy14) == 1
    @test get_dataset(test_divacancy14) == "dataset"
    @test get_noa(test_divacancy14) == 14
    @test length(get_relatives(test_divacancy14)) == 10
    @test length(get_members(test_divacancy14)) == 29
    @test get_type(test_divacancy14) == "Divacancy-14"
end

@testset "Flower" begin
    @test get_id(test_flower) == length(test_graphene) + 1
    @test get_frame(test_flower) == 1
    @test get_dataset(test_flower) == "dataset"
    @test get_noa(test_flower) == 22
    @test length(get_relatives(test_flower)) == 12
    @test length(get_members(test_flower)) == 55
    @test get_type(test_flower) == "Flower"
end

@testset "Butterfly" begin
    @test get_id(test_butterfly) == length(test_graphene) + 1
    @test get_frame(test_butterfly) == 1
    @test get_dataset(test_butterfly) == "dataset"
    @test get_noa(test_butterfly) == 30
    @test length(get_relatives(test_butterfly)) == 14
    @test length(get_members(test_butterfly)) == 77
    @test get_type(test_butterfly) == "Butterfly"
end


text:
lines = readlines(ARGS[1])

function check_boards(check_val, boards)
    for (i, o) in enumerate(boards)
        if board_wins(check_val, boards[i][1], boards[i][2])
            return boards[i][1], boards[i][2]
        else
            continue
        end
    end

    return nothing, nothing
end

function board_wins(check_val, check_board, result_board)
    for (row_i, row_v) in enumerate(check_board)
        for (col_i, col_v) in enumerate(row_v)
            if col_v == check_val
                result_board[row_i][col_i] = check_val
                break
            end
        end
    end
    if check_rows(result_board)
        return true
    elseif check_columns(result_board)
        return true
    else
        return false
    end
end

function check_columns(result_board)
    for (i, v) in enumerate(result_board[1])
        if !isnothing(v)
            if size([row[i] for row in result_board if row[i] != ""]) == (5,)
                return true
            end
        end
    end
    return false
end

function check_rows(result_board)
    for row in result_board
        if size([x for x in row if x != ""]) == (5,)
            return true
        end
    end
    return false
end

function solution1(lines)
    moves = nothing
    boards = []
    board = []
    for (iter_val, line) in enumerate(lines)
        if iter_val == 1
            moves = split(line, ",")
        elseif iter_val > 2
            if line == ""
                append!(boards, [(board, [["" for i in 1:5] for i in 1:5])])
                board = []
            else
                append!(board, [split(line)])
            end
        else
            continue
        end
    end

    (win_board, win_positions) = (nothing, nothing)
    skip = Int[]
    for move in moves
        for (i, o) in enumerate(boards)
            if i in skip
                continue
            end
            for (row_i, row_v) in enumerate(boards[i][1])
                for (col_i, col_v) in enumerate(row_v)
                    if col_v == move
                        boards[i][2][row_i][col_i] = move
                        break
                        break
                    end
                end
            end
            if check_rows(boards[i][2])
                (win_board, win_positions) = boards[i]
                append!(skip, i)
            elseif check_columns(boards[i][2])
                (win_board, win_positions) = boards[i]
                append!(skip, i)
            end
        end

        if isnothing(win_board)
            continue
        else
            sum_val = 0
            for (row_i, row_v) in enumerate(win_positions)
                for (col_i, col_v) in enumerate(row_v)
                    if col_v == ""
                        sum_val += parse(Int, win_board[row_i][col_i])
                    end
                end
            end
            println(sum_val*parse(Int, move))

        end
        (win_board, win_positions) = (nothing, nothing)

    end
    # moves = parse_moves(lines[1])
    # lines = Array(lines)
    # println(lines[3:size(lines)])

    # boards = parse_games(lines[3:-1])
    return nothing
end

# function solution2(lines)

# end

solution1(lines)


text:
using UUIDs

"""
`AbstractScheduler` is an abstract type for schedulers.

Schedulers are structs which are responsible of running
`Action` at given instants (according a `Trigger`).

Several kind of schedulers can implement `AbstractScheduler`.

The most simple scheduler is `BlockingScheduler` which is monothread.
"""
abstract type AbstractScheduler end

include("background.jl")
include("blocking.jl")

"""
    get_scheduler_id()

Return a random identifier for a scheduler.
"""
function get_scheduler_id()
    string(UUIDs.uuid4())[1:13]
end


text:
const HELICS_HAVE_MPI = 0
const HELICS_HAVE_ZEROMQ = 1
const BOOST_VERSION_LEVEL = 71
const HELICS_VERSION_MAJOR = 2
const HELICS_VERSION_MINOR = 3
const HELICS_VERSION_PATCH = 0
const HELICS_VERSION = "2.3.0"
const HELICS_VERSION_BUILD = ""
const HELICS_VERSION_STRING = "2.3.0 (2019-11-13)"
const HELICS_DATE = "2019-11-13"