text:
using SyntheticGrids
using Test
using JSON

@testset "SyntheticGrids" begin
    SEED = 666
    LATLIM = (38, 40)
    LONGLIM = (-89, -88)
    GENPATH = joinpath(@__DIR__, "..", "data", "GenData.json")
    CENSUSPATH = joinpath(@__DIR__, "..", "data", "Census_data.dat")
    GENCOORDPATH = joinpath(@__DIR__, "..", "data", "Generator_coord.dat")
    GENDATAPATH = joinpath(@__DIR__, "..", "data", "Generator_data.dat")
    DUMMYPATH = joinpath(@__DIR__, "..", "data", "dummy.json")
    DUMMYPPC = joinpath(@__DIR__, "..", "data", "dummy.p")
    GRIDPATH = joinpath(@__DIR__, "..", "data", "testgrid.json")

    function baseline_grid(seed=666, latlim=LATLIM, longlim=LONGLIM)
        grid = Grid(seed)
        place_loads_from_zips!(grid; latlim=latlim, longlim=longlim)
        place_gens_from_data!(grid; latlim=latlim, longlim=longlim)
        connect!(grid)
        create_lines!(grid)
        return grid
    end

    @testset "Grids" begin
        @testset "Load bus length" begin
            grid = Grid(SEED)
            place_loads_from_zips!(grid; latlim=LATLIM, longlim=LONGLIM)
            @test length(buses(grid)) == 130
        end

        @testset "Gen bus length" begin
            grid = Grid(SEED)
            place_gens_from_data!(grid; latlim=LATLIM, longlim=LONGLIM)
            @test length(buses(grid)) == 17
        end

        @testset "Adjacency sum" begin
            grid = baseline_grid()
            @test sum(adjacency(grid)) == 358
        end

        @testset "Bus type count" begin
            grid = baseline_grid()
            count = SyntheticGrids.count_bus_type(grid)
            @test count == (130, 17)
        end

        @testset "Trans lines" begin
            grid = baseline_grid()
            @test grid.trans_lines[1].capacity == 4900
            @test length(trans_lines(grid)) == 179
            @test test_connectivity(grid.bus_conn, false)
            @test total_links(grid.bus_conn) == 179
            @test mean_node_deg(grid.bus_conn) ≈ 2.435374149659
            @test mean_shortest_path(adjacency(grid)) ≈ 7.693646165949373
            @test mean_shortest_path(adjacency(grid), distance(buses(grid))) ≈ 133.34235584030938
            @test robustness_line(grid.bus_conn, 10) > 1
            @test robustness_node(grid.bus_conn, 10) > 1
        end

        @testset "Add substation" begin
            grid = Grid(true)
            @test size(grid.bus_conn) == (137, 137)
            @test size(grid.sub_conn) == (43, 43)

            add_gen!(grid, (11., 11.), 100, [11], ["ki"], reconnect = false)
            add_load!(grid, (19.,12.), 100, 100, 100, reconnect = true)
            add_substation!(
                grid,
                (11., 18.),
                [11],
                1000,
                6000,
                500,
                Set(),
                [grid.buses[139], grid.buses[138]]
            )
            @test size(grid.bus_conn) == (139, 139)
            @test size(grid.sub_conn) == (44, 44)
        end

        @testset "Clustering grids" begin
            grid = baseline_grid()
            count = SyntheticGrids.count_bus_type(grid)
            cluster!(
                grid,
                round(Int, 0.3 * count[1]),
                round(Int, 0.05 * count[2]),
                round(Int, 0.5 * count[2])
            )
            @test length(substations(grid)) == 47
            @test grid.substations[1].population == 6193
            @test grid.substations[47].generation == 16.1
            @test test_connectivity(grid.sub_conn, false)
            @test total_links(grid.sub_conn) == 70
            @test mean_shortest_path(sub_connectivity(grid) .> 0) ≈ 4.1620642824807605
            @test mean_shortest_path(
                (sub_connectivity(grid) .> 0),
                distance(substations(grid))
                ) ≈ 111.6094877525516
            @test robustness_line(grid.sub_conn, 10) > 1
            @test robustness_node(grid.sub_conn, 10) > 1
        end

        @testset "Test seeding and object comparison" begin
            grid = baseline_grid()
            grid1 = baseline_grid()

            @test grid == grid1
        end
    end

    @testset "Input and Output" begin
        @testset "Pandapower interface" begin
            grid = baseline_grid()
            pgrid = to_pandapower(grid, DUMMYPPC)
            @test length(pgrid.trafo) == 31
            pgrid = SyntheticGrids.load_pp_grid(DUMMYPPC)
            @test length(pgrid.trafo) == 31
            rm(DUMMYPPC)
        end

        @testset "Generator data" begin
            SyntheticGrids.prepare_gen_data(GENDATAPATH, GENCOORDPATH, DUMMYPATH)
            genjson = JSON.parsefile(GENPATH)
            dummyjson = JSON.parsefile(DUMMYPATH)
            @test genjson == dummyjson
            rm(DUMMYPATH)
        end

        @testset "Save and Load Grid" begin
            grid = baseline_grid()
            save(grid, GRIDPATH)
            lgrid = load_grid(GRIDPATH)
            @test adjacency(grid) == adjacency(lgrid)
            @test sub_connectivity(grid) == sub_connectivity(lgrid)
            lgrid = 0 # trying to solve issues with Windows locking permissions to GRIDPATH
            rm(GRIDPATH)
        end
    end
end


text:
function _true_data(model_args)
  P0 = 50.0
  Q0 = 5.0
  r = 0.6
  K = 100.0
  s = 1.2
  a = 25.0
  u = 0.5
  v = 0.3
  parameter = PredatorPreyParameter(P0, Q0, r, K, s, a, u, v, model_args[3:4]..., nothing)
  stable!(parameter)
  generate(predator_prey_model, model_args)[1]
end

_model_args = (0.1, 110.0, 0.0, 50.0, 20, 1.0)
_trace = _true_data(_model_args)
data = choicemap(map(p -> (:measurements => p[1] => :z, p[2]), enumerate(_trace[:measurements]))...)

ext: jl

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
deval{N}(::Val{N}) = N
padding{IT<:BSpline}(::Type{IT}) = Val{0}()
@generated function padding{IT}(t::Type{IT})
    pad = [deval(padding(IT.parameters[d])) for d = 1:length(IT.parameters)]
    t = tuple(pad...)
    :(Val{$t}())
end

function padded_index{N,pad}(sz::NTuple{N,Int}, ::Val{pad})
    szpad = ntuple(i->sz[i]+2padextract(pad,i), N)::NTuple{N,Int}
    ind = Array(UnitRange{Int},N)
    for i in 1:N
        p = padextract(pad,i)
        ind[i] = 1+p:szpad[i]-p
    end
    ind,szpad
end

copy_with_padding{IT}(A, ::Type{IT}) = copy_with_padding(eltype(A), A, IT)
function copy_with_padding{TC,IT<:DimSpec{InterpolationType}}(::Type{TC}, A, ::Type{IT})
    Pad = padding(IT)
    ind,sz = padded_index(size(A), Pad)
    if sz == size(A)
        coefs = copy!(Array(TC, size(A)), A)
    else
        coefs = zeros(TC, sz...)
        coefs[ind...] = A
    end
    coefs, Pad
end

prefilter!{TWeights, IT<:BSpline, GT<:GridType}(::Type{TWeights}, A, ::Type{IT}, ::Type{GT}) = A
prefilter{TWeights, TC, IT<:BSpline, GT<:GridType}(::Type{TWeights}, ::Type{TC}, A, ::Type{IT}, ::Type{GT}) = prefilter!(TWeights, copy!(Array(TC, size(A)), A), IT, GT), Val{0}()

function prefilter{TWeights,TC,IT<:Union{Cubic,Quadratic},GT<:GridType}(
    ::Type{TWeights}, ::Type{TC}, A::AbstractArray, ::Type{BSpline{IT}}, ::Type{GT}
    )
    ret, Pad = copy_with_padding(TC, A, BSpline{IT})
    prefilter!(TWeights, ret, BSpline{IT}, GT), Pad
end

function prefilter{TWeights,TC,IT<:Tuple{Vararg{Union{BSpline,NoInterp}}},GT<:DimSpec{GridType}}(
    ::Type{TWeights}, ::Type{TC}, A::AbstractArray, ::Type{IT}, ::Type{GT}
    )
    ret, Pad = copy_with_padding(TC, A, IT)
    prefilter!(TWeights, ret, IT, GT), Pad
end

function prefilter!{TWeights,TCoefs<:AbstractArray,IT<:Union{Quadratic,Cubic},GT<:GridType}(
    ::Type{TWeights}, ret::TCoefs, ::Type{BSpline{IT}}, ::Type{GT}
    )
    local buf, shape, retrs
    sz = size(ret)
    first = true
    for dim in 1:ndims(ret)
        M, b = prefiltering_system(TWeights, eltype(TCoefs), sz[dim], IT, GT)
        A_ldiv_B_md!(ret, M, ret, dim, b)
    end
    ret
end

function prefilter!{TWeights,TCoefs<:AbstractArray,IT<:Tuple{Vararg{Union{BSpline,NoInterp}}},GT<:DimSpec{GridType}}(
    ::Type{TWeights}, ret::TCoefs, ::Type{IT}, ::Type{GT}
    )
    local buf, shape, retrs
    sz = size(ret)
    first = true
    for dim in 1:ndims(ret)
        it = iextract(IT, dim)
        if it != NoInterp
            M, b = prefiltering_system(TWeights, eltype(TCoefs), sz[dim], bsplinetype(it), iextract(GT, dim))
            if M != nothing
                A_ldiv_B_md!(ret, M, ret, dim, b)
            end
        end
    end
    ret
end

prefiltering_system(::Any, ::Any, ::Any, ::Any, ::Any) = nothing, nothing

"""
    M, b = prefiltering_system{T,TC,GT<:GridType,D<:Degree}m(::T, ::Type{TC}, n::Int, ::Type{D}, ::Type{GT})

Given element types (`T`, `TC`) and interpolation scheme
(`GT`, `D`) as well the number of rows in the data input (`n`), compute the
system used to prefilter spline coefficients. Boundary conditions determine the
values on the first and last rows.

Some of these boundary conditions require that these rows have off-tridiagonal
elements (e.g the `[1,3]` element of the matrix). To maintain the efficiency of
solving tridiagonal systems, the [Woodbury matrix identity](https://en.wikipedia.org/wiki/Woodbury_matrix_identity)
is used to add additional elements off the main 3 diagonals.

The filtered coefficients are given by solving the equation system

    M * c = v + b

where `c` are the sought coefficients, and `v` are the data points.
"""
prefiltering_system

"""
    dl, d, du = inner_system_diags{T,IT}(::Type{T}, n::Int, ::Type{IT})

Helper function to generate the prefiltering equation system: generates the diagonals
for a `n`-by-`n` tridiagonal matrix with eltype `T` corresponding to the interpolation
type `IT`.

`dl`, `d`, and `du` are intended to be used e.g. as in `M = Tridiagonal(dl, d, du)`
"""
inner_system_diags


text:
function _quicksort!(a::Vector, leftindex, rightindex, isless)
    # Note on chosen optimization:
    # Shaffer recommends to *just* return the small partitions and do a single final insertionsort. We follow Sedgewick
    # here, who recommends to do a insertionsort for the many small partitions. This is the same technique applied to
    # mergesort (see https://algs4.cs.princeton.edu/23quicksort/).
    if (rightindex + 1 - leftindex)  <= 9
        insertionsort!(view(a, leftindex:rightindex), isless)
        return
    end
    pivotindex = _findpivot(a, leftindex, rightindex)
    swap!(a, pivotindex, rightindex)
    newpivotindex = _partition!(a, leftindex-1, rightindex, a[rightindex], isless)
    swap!(a, newpivotindex, rightindex)
    _quicksort!(a, leftindex, newpivotindex-1, isless)
    _quicksort!(a, newpivotindex+1, rightindex, isless)
end

@inline function _findpivot(a, leftindex, rightindex)
    return (leftindex + rightindex) ÷ 2
end

@inline function _partition!(a, leftindex, rightindex, pivot, isless)
    let l = leftindex, r = rightindex
        while true
            l += 1
            while isless(a[l], pivot)
                l += 1
            end
            l < r || break
            r -= 1
            while isless(pivot, a[r])
                l < r || break
                r -= 1
            end
            swap!(a, l, r) # gets bound error
            l < r || break
        end
        return l
    end
end

function quicksort!(a::Vector, isless=isless)
    _quicksort!(a,  1, length(a), isless)
    a
end


text:
# Autogenerated wrapper script for NCCL_jll for powerpc64le-linux-gnu-cuda+10.2
export libnccl

using CUDA_loader_jll
JLLWrappers.@generate_wrapper_header("NCCL")
JLLWrappers.@declare_library_product(libnccl, "libnccl.so.2")
function __init__()
    JLLWrappers.@generate_init_header(CUDA_loader_jll)
    JLLWrappers.@init_library_product(
        libnccl,
        "lib/libnccl.so",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"Probability Space indexed with values of type I"
abstract type Ω{I} <: Random.AbstractRNG end

"This is base Omega - Sample Space Object"
abstract type ΩBase{I} <: Ω{I} end

idtype(::Type{OT}) where {I, OT <: Ω{I}} = I

# Must implement:
# * Base.==
# * Base.isempty
# * Base.merge! ( -- currently unused)
# ? others
"menrand(ωπ, args...).  Memoized `rand`.  If ωπ.i ∈ ωπ.ω, return ωπ.ω[ωπ.i] otherwise call `rand(args...)`"
function memrand end

"Linearize ω into flat vector"
function linearize end

"Inverse of `linearize`, structure vector into ω"
function unlinearize end

function nelem end

## Rand

"Random ω ∈ Ω"
Base.rand(x::Type{O}) where O <: ΩBase = O()
Base.rand(rng::AbstractRNG, ::Type{T}) where {T <: Ω} = tagrng(rand(T), rng)

Random.rng_native_52(ω::Ω) = Random.rng_native_52(Random.GLOBAL_RNG)


text:
module projection

using LinearAlgebra

export distance_matrix_to_absolute_coordinates

function distance_matrix_to_absolute_coordinates( D::Array{Float64,2} , bSquared::Bool=false , n_dimensions::Int64=2 )
    """
    Same function can be found in my impetuous-gfa repo here:
    https://github.com/richardtjornhammar/impetuous/blob/master/src/impetuous/clustering.py
    around line 121 (same name)
    """
    if !bSquared
        D = D.^2.
    end
    N,M = size(D)
    DIM = n_dimensions
    DIJ = D.*0.
    for i in 1:M
        for j in 1:M
            DIJ[i,j] = 0.5* (D[i,end]+D[j,end]-D[i,j])
	end
    end
    D = DIJ
    U,S,Vt = svd( D )
    S[DIM+1:end] *= 0.
    Z  = diagm(S.^0.5)[:,1:DIM]
    xr = Vt*Z
    return ( xr )
end

end

max_stars_repo_name: richardtjornhammar/ingenuous

max_issues_repo_name: richardtjornhammar/ingenuous

max_forks_repo_name: richardtjornhammar/ingenuous


text:
# include the package
using LatticePhysics

# obtain the unitcell of the honeycomb lattice
unitcell = getUnitcellHoneycomb()

# build a lattice for visual confirmation that this is the correct unitcell
lattice = getLatticeInBox(unitcell, [10.0, 8.0])
# plot the lattice for visual confirmation
plotLattice2D(lattice, openfile=true)

# get a path for band structure calculations, use the default path
path = getDefaultPathTriangular()
# calculate the band structure (using the unitcell)
plotBandstructure(unitcell, path, showPlot=true)


text:
@doc raw"""
    PrimalDualProblem {T, mT <: Manifold, nT <: Manifold} <: PrimalDualProblem} <: Problem{T}

Describes a Problem for the linearized or exact Chambolle-Pock algorithm.[^BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez2020][^ChambollePock2011]

# Fields

All fields with !! can either be mutating or nonmutating functions, which should be set
depenting on the parameter `T <: AbstractEvaluationType`.

* `M`, `N` – two manifolds ``\mathcal M``, ``\mathcal N``
* `cost` ``F + G(Λ(⋅))`` to evaluate interims cost function values
* `linearized_forward_operator!!` linearized operator for the forward operation in the algorthm ``DΛ``
* `linearized_adjoint_operator!!` The adjoint differential ``(DΛ)^* : \mathcal N → T\mathcal M``
* `prox_F!!` the proximal map belonging to ``f``
* `prox_G_dual!!` the proximal map belonging to ``g_n^*``
* `Λ!!` – (`fordward_operator`) the  forward operator (if given) ``Λ: \mathcal M → \mathcal N``

Either ``DΛ`` (for the linearized) or ``Λ`` are required usually.

# Constructor

    LinearizedPrimalDualProblem(M, N, cost, prox_F, prox_G_dual, adjoint_linearized_operator;
        linearized_forward_operator::Union{Function,Missing}=missing,
        Λ::Union{Function,Missing}=missing,
        evaluation::AbstractEvaluationType=AllocatingEvaluation()
    )

The last optional argument can be used to provide the 4 or 5 functions as allocating or
mutating (in place computation) ones.
Note that the first argument is always the manifold under consideration, the mutated one is
the second.

[^BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez2020]:
    > R. Bergmann, R. Herzog, M. Silva Louzeiro, D. Tenbrinck, J. Vidal-Núñez:
    > _Fenchel Duality Theory and a Primal-Dual Algorithm on Riemannian Manifolds_,
    > Foundations of Computational Mathematics, 2021.
    > doi: [10.1007/s10208-020-09486-5](http://dx.doi.org/10.1007/s10208-020-09486-5)
    > arXiv: [1908.02022](http://arxiv.org/abs/1908.02022)
[^ChambollePock2011]:
    > A. Chambolle, T. Pock:
    > _A first-order primal-dual algorithm for convex problems with applications to imaging_,
    > Journal of Mathematical Imaging and Vision 40(1), 120–145, 2011.
    > doi: [10.1007/s10851-010-0251-1](https://dx.doi.org/10.1007/s10851-010-0251-1)
"""
mutable struct PrimalDualProblem{T,mT<:Manifold,nT<:Manifold} <: Problem{T}
    M::mT
    N::nT
    cost::Function
    prox_F!!::Function
    prox_G_dual!!::Function
    linearized_forward_operator!!::Union{Function,Missing}
    adjoint_linearized_operator!!::Function
    Λ!!::Union{Function,Missing}
end
function PrimalDualProblem(
    M::mT,
    N::nT,
    cost,
    prox_F,
    prox_G_dual,
    adjoint_linearized_operator;
    linearized_forward_operator::Union{Function,Missing}=missing,
    Λ::Union{Function,Missing}=missing,
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
) where {mT<:Manifold,nT<:Manifold}
    return PrimalDualProblem{typeof(evaluation),mT,nT}(
        M,
        N,
        cost,
        prox_F,
        prox_G_dual,
        linearized_forward_operator,
        adjoint_linearized_operator,
        Λ,
    )
end

@doc raw"""
    y = get_primal_prox(p::PrimalDualProblem, σ, x)
    get_primal_prox!(p::PrimalDualProblem, y, σ, x)

Evaluate the proximal map of ``F`` stored within [`PrimalDualProblem`](@ref)

```math
\operatorname{prox}_{σF}(x)
```

which can also be computed in place of `y`.
"""
get_primal_prox(::PrimalDualProblem, ::Any...)

function get_primal_prox(p::PrimalDualProblem{AllocatingEvaluation}, σ, x)
    return p.prox_F!!(p.M, σ, x)
end
function get_primal_prox(p::PrimalDualProblem{MutatingEvaluation}, σ, x)
    y = allocate_result(p.M, get_primal_prox, x)
    return p.prox_F!!(p.M, y, σ, x)
end
function get_primal_prox!(p::PrimalDualProblem{AllocatingEvaluation}, y, σ, x)
    return copyto!(y, p.prox_F!!(p.M, σ, x))
end
function get_primal_prox!(p::PrimalDualProblem{MutatingEvaluation}, y, σ, x)
    return p.prox_F!!(p.M, y, σ, x)
end

@doc raw"""
    y = get_dual_prox(p::PrimalDualProblem, n, τ, ξ)
    get_dual_prox!(p::PrimalDualProblem, y, n, τ, ξ)

Evaluate the proximal map of ``G_n^*`` stored within [`PrimalDualProblem`](@ref)

```math
\operatorname{prox}_{τG_n^*}(ξ)
```

which can also be computed in place of `y`.
"""
get_dual_prox(::PrimalDualProblem, ::Any...)

function get_dual_prox(p::PrimalDualProblem{AllocatingEvaluation}, n, τ, ξ)
    return p.prox_G_dual!!(p.N, n, τ, ξ)
end
function get_dual_prox(p::PrimalDualProblem{MutatingEvaluation}, n, τ, ξ)
    η = allocate_result(p.N, get_dual_prox, ξ)
    return p.prox_G_dual!!(p.N, η, n, τ, ξ)
end
function get_dual_prox!(p::PrimalDualProblem{AllocatingEvaluation}, η, n, τ, ξ)
    return copyto!(η, p.prox_G_dual!!(p.N, n, τ, ξ))
end
function get_dual_prox!(p::PrimalDualProblem{MutatingEvaluation}, η, n, τ, ξ)
    return p.prox_G_dual!!(p.N, η, n, τ, ξ)
end
@doc raw"""
    Y = linearized_forward_operator(p::PrimalDualProblem, m X)
    linearized_forward_operator!(p::PrimalDualProblem, Y, m, X)

Evaluate the linearized operator (differential) ``DΛ(m)[X]`` stored within
the [`PrimalDualProblem`](@ref) (in place of `Y`).
"""
linearized_forward_operator(::PrimalDualProblem, ::Any...)

function linearized_forward_operator(p::PrimalDualProblem{AllocatingEvaluation}, m, X)
    return p.linearized_forward_operator!!(p.M, m, X)
end
function linearized_forward_operator(p::PrimalDualProblem{MutatingEvaluation}, m, X)
    y = random_point(p.N)
    forward_operator!(p, y, m)
    Y = zero_tangent_vector(p.N, y)
    return p.linearized_forward_operator!!(p.M, Y, m, X)
end
function linearized_forward_operator!(p::PrimalDualProblem{AllocatingEvaluation}, Y, m, X)
    return copyto!(Y, p.linearized_forward_operator!!(p.M, m, X))
end
function linearized_forward_operator!(p::PrimalDualProblem{MutatingEvaluation}, Y, m, X)
    return p.linearized_forward_operator!!(p.M, Y, m, X)
end

@doc raw"""
    y = forward_operator(p::PrimalDualProblem, x)
    forward_operator!(p::PrimalDualProblem, y, x)

Evaluate the forward operator of ``Λ(x)`` stored within the [`PrimalDualProblem`](@ref)
(in place of `y`).
"""
forward_operator(::PrimalDualProblem{AllocatingEvaluation}, ::Any...)

function forward_operator(p::PrimalDualProblem{AllocatingEvaluation}, x)
    return p.Λ!!(p.M, x)
end
function forward_operator(p::PrimalDualProblem{MutatingEvaluation}, x)
    y = random_point(p.N)
    return p.Λ!!(p.M, y, x)
end
function forward_operator!(p::PrimalDualProblem{AllocatingEvaluation}, y, x)
    return copyto!(y, p.Λ!!(p.M, x))
end
function forward_operator!(p::PrimalDualProblem{MutatingEvaluation}, y, x)
    return p.Λ!!(p.M, y, x)
end

@doc raw"""
    X = adjoint_linearized_operator(p::PrimalDualProblem, m, n, Y)
    adjoint_linearized_operator(p::PrimalDualProblem, X, m, n, Y)

Evaluate the adjoint of the linearized forward operator of ``(DΛ(m))^*[Y]`` stored within
the [`PrimalDualProblem`](@ref) (in place of `X`).
Since ``Y∈T_n\mathcal N``, both ``m`` and ``n=Λ(m)`` are necessary arguments, mainly because
the forward operator ``Λ`` might be `missing` in `p`.
"""
adjoint_linearized_operator(::PrimalDualProblem{AllocatingEvaluation}, ::Any...)

function adjoint_linearized_operator(p::PrimalDualProblem{AllocatingEvaluation}, m, n, Y)
    return p.adjoint_linearized_operator!!(p.N, m, n, Y)
end
function adjoint_linearized_operator(p::PrimalDualProblem{MutatingEvaluation}, m, n, Y)
    X = zero_tangent_vector(p.M, m)
    return p.adjoint_linearized_operator!!(p.N, X, m, n, Y)
end
function adjoint_linearized_operator!(
    p::PrimalDualProblem{AllocatingEvaluation}, X, m, n, Y
)
    return copyto!(X, p.adjoint_linearized_operator!!(p.N, m, n, Y))
end
function adjoint_linearized_operator!(p::PrimalDualProblem{MutatingEvaluation}, X, m, n, Y)
    return p.adjoint_linearized_operator!!(p.N, X, m, n, Y)
end

@doc raw"""
    PrimalDualOptions

A general type for all primal dual based options to be used within primal dual
based algorithms
"""
abstract type PrimalDualOptions <: Options end

@doc raw"""
    ChambollePockOptions <: PrimalDualOptions

stores all options and variables within a linearized or exact Chambolle Pock.
The following list provides the order for the constructor, where the previous iterates are
initialized automatically and values with a default may be left out.

* `m` - base point on ``\mathcal M``
* `n` - base point on ``\mathcal N``
* `x` - an initial point on ``x^{(0)} ∈\mathcal M`` (and its previous iterate)
* `ξ` - an initial tangent vector ``ξ^{(0)}∈T^*\mathcal N`` (and its previous iterate)
* `xbar` - the relaxed iterate used in the next dual update step (when using `:primal` relaxation)
* `ξbar` - the relaxed iterate used in the next primal update step (when using `:dual` relaxation)
* `Θ` – factor to damp the helping ``\tilde x``
* `primal_stepsize` – (`1/sqrt(8)`) proximal parameter of the primal prox
* `dual_stepsize` – (`1/sqrt(8)`) proximnal parameter of the dual prox
* `acceleration` – (`0.`) acceleration factor due to Chambolle & Pock
* `relaxation` – (`1.`) relaxation in the primal relaxation step (to compute `xbar`)
* `relax` – (`_primal`) which variable to relax (`:primal` or `:dual`)
* `stop` - a [`StoppingCriterion`](@ref)
* `type` – (`exact`) whether to perform an `:exact` or `:linearized` Chambolle-Pock
* `update_primal_base` (`(p,o,i) -> o.m`) function to update the primal base
* `update_dual_base` (`(p,o,i) -> o.n`) function to update the dual base
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use
* `inverse_retraction_method` - (`LogarithmicInverseRetraction()`) an inverse retraction to use.
* `vector_transport_method` - (`ParallelTransport()`) a vector transport to use

where for the last two the functions a [`Problem`](@ref)` p`,
[`Options`](@ref)` o` and the current iterate `i` are the arguments.
If you activate these to be different from the default identity, you have to provide
`p.Λ` for the algorithm to work (which might be `missing` in the linearized case).

# Constructor
    ChambollePockOptions(m::P, n::Q, x::P, ξ::T, primal_stepsize::Float64, dual_stepsize::Float64;
        acceleration::Float64 = 0.0,
        relaxation::Float64 = 1.0,
        relax::Symbol = :primal,
        stopping_criterion::StoppingCriterion = StopAfterIteration(300),
        variant::Symbol = :exact,
        update_primal_base::Union{Function,Missing} = missing,
        update_dual_base::Union{Function,Missing} = missing,
        retraction_method = ExponentialRetraction(),
        inverse_retraction_method = LogarithmicInverseRetraction(),
        vector_transport_method = ParallelTransport(),
    )
"""
mutable struct ChambollePockOptions{
    P,
    Q,
    T,
    RM<:AbstractRetractionMethod,
    IRM<:AbstractInverseRetractionMethod,
    VTM<:AbstractVectorTransportMethod,
} <: PrimalDualOptions
    m::P
    n::Q
    x::P
    xbar::P
    ξ::T
    ξbar::T
    primal_stepsize::Float64
    dual_stepsize::Float64
    acceleration::Float64
    relaxation::Float64
    relax::Symbol
    stop::StoppingCriterion
    variant::Symbol
    update_primal_base::Union{Function,Missing}
    update_dual_base::Union{Function,Missing}
    retraction_method::RM
    inverse_retraction_method::IRM
    vector_transport_method::VTM

    function ChambollePockOptions(
        m::P,
        n::Q,
        x::P,
        ξ::T,
        primal_stepsize::Float64=1 / sqrt(8),
        dual_stepsize::Float64=1 / sqrt(8);
        acceleration::Float64=0.0,
        relaxation::Float64=1.0,
        relax::Symbol=:primal,
        stopping_criterion::StoppingCriterion=StopAfterIteration(300),
        variant::Symbol=:exact,
        update_primal_base::Union{Function,Missing}=missing,
        update_dual_base::Union{Function,Missing}=missing,
        retraction_method::RM=ExponentialRetraction(),
        inverse_retraction_method::IRM=LogarithmicInverseRetraction(),
        vector_transport_method::VTM=ParallelTransport(),
    ) where {
        P,
        Q,
        T,
        RM<:AbstractRetractionMethod,
        IRM<:AbstractInverseRetractionMethod,
        VTM<:AbstractVectorTransportMethod,
    }
        return new{P,Q,T,RM,IRM,VTM}(
            m,
            n,
            x,
            deepcopy(x),
            ξ,
            deepcopy(ξ),
            primal_stepsize,
            dual_stepsize,
            acceleration,
            relaxation,
            relax,
            stopping_criterion,
            variant,
            update_primal_base,
            update_dual_base,
            retraction_method,
            inverse_retraction_method,
            vector_transport_method,
        )
    end
end
get_solver_result(o::ChambollePockOptions) = o.x
@doc raw"""
    primal_residual(p, o, x_old, ξ_old, n_old)

Compute the primal residual at current iterate ``k`` given the necessary values ``x_{k-1},
ξ_{k-1}``, and ``n_{k-1}`` from the previous iterate.

```math
\Bigl\lVert
\frac{1}{σ}\operatorname{retr}^{-1}_{x_{k}}x_{k-1} -
V_{x_k\gets m_k}\bigl(DΛ^*(m_k)\bigl[V_{n_k\gets n_{k-1}}ξ_{k-1} - ξ_k \bigr]
\Bigr\rVert
```
where ``V_{⋅\gets⋅}`` is the vector transport used in the [`ChambollePockOptions`](@ref)
"""
function primal_residual(p::PrimalDualProblem, o::ChambollePockOptions, x_old, ξ_old, n_old)
    return norm(
        p.M,
        o.x,
        1 / o.primal_stepsize *
        inverse_retract(p.M, o.x, x_old, o.inverse_retraction_method) -
        vector_transport_to(
            p.M,
            o.m,
            adjoint_linearized_operator(
                p,
                o.m,
                o.n,
                vector_transport_to(p.N, n_old, ξ_old, o.n, o.vector_transport_method) -
                o.ξ,
            ),
            o.x,
            o.vector_transport_method,
        ),
    )
end
@doc raw"""
    dual_residual(p, o, x_old, ξ_old, n_old)

Compute the dual residual at current iterate ``k`` given the necessary values ``x_{k-1},
ξ_{k-1}``, and ``n_{k-1}`` from the previous iterate. The formula is slightly different depending
on the `o.variant` used:

For the `:lineaized` it reads
```math
\Bigl\lVert
\frac{1}{τ}\bigl(
V_{n_{k}\gets n_{k-1}}(ξ_{k-1})
- ξ_k
\bigr)
-
DΛ(m_k)\bigl[
V_{m_k\gets x_k}\operatorname{retr}^{-1}_{x_{k}}x_{k-1}
\bigr]
\Bigr\rVert
```

and for the `:exact` variant

```math
\Bigl\lVert
\frac{1}{τ} V_{n_{k}\gets n_{k-1}}(ξ_{k-1})
-
\operatorname{retr}^{-1}_{n_{k}}\bigl(
Λ(\operatorname{retr}_{m_{k}}(V_{m_k\gets x_k}\operatorname{retr}^{-1}_{x_{k}}x_{k-1}))
\bigr)
\Bigr\rVert
```

where in both cases ``V_{⋅\gets⋅}`` is the vector transport used in the [`ChambollePockOptions`](@ref).
"""
function dual_residual(p::PrimalDualProblem, o::ChambollePockOptions, x_old, ξ_old, n_old)
    if o.variant === :linearized
        return norm(
            p.N,
            o.n,
            1 / o.dual_stepsize *
            (vector_transport_to(p.N, n_old, ξ_old, o.n, o.vector_transport_method) - o.ξ) -
            linearized_forward_operator(
                p,
                o.m,
                vector_transport_to(
                    p.M,
                    o.x,
                    inverse_retract(p.M, o.x, x_old, o.inverse_retraction_method),
                    o.m,
                    o.vector_transport_method,
                ),
            ),
        )
    elseif o.variant === :exact
        return norm(
            p.N,
            o.n,
            1 / o.dual_stepsize *
            (vector_transport_to(p.N, n_old, ξ_old, o.n, o.vector_transport_method) - o.n) -
            inverse_retract(
                p.N,
                o.n,
                forward_operator(
                    p,
                    retract(
                        p.M,
                        o.m,
                        vector_transport_to(
                            p.M,
                            o.x,
                            inverse_retract(p.M, o.x, x_old, o.inverse_retraction_method),
                            o.m,
                            o.vector_transport_method,
                        ),
                        o.retraction_method,
                    ),
                ),
                o.inverse_retraction_method,
            ),
        )
    else
        throw(
            DomainError(
                o.variant,
                "Unknown Chambolle–Pock variant, allowed are `:exact` or `:linearized`.",
            ),
        )
    end
end
#
# Special Debuggers
#
@doc raw"""
    DebugDualResidual <: DebugAction

A Debug action to print the dual residual.
The constructor accepts a printing function and some (shared) storage, which
should at least record `:x`, `:ξ` and `:n`.
"""
mutable struct DebugDualResidual <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugDualResidual(
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)), io::IO=stdout
    )
        return new(io, "Dual Residual: ", a)
    end
    function DebugDualResidual(
        values::Tuple{P,T,Q},
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)),
        io::IO=stdout,
    ) where {P,T,Q}
        update_storage!(a, Dict(k => v for (k, v) in zip((:x, :ξ, :n), values)))
        return new(io, "Dual Residual: ", a)
    end
end
function (d::DebugDualResidual)(p::PrimalDualProblem, o::ChambollePockOptions, i::Int)
    if all(has_storage.(Ref(d.storage), [:x, :ξ, :n])) && i > 0 # all values stored
        xOld, ξOld, nOld = get_storage.(Ref(d.storage), [:x, :ξ, :n]) #fetch
        print(d.io, d.prefix * string(dual_residual(p, o, xOld, ξOld, nOld)))
    end
    return d.storage(p, o, i)
end

@doc raw"""
    DebugPrimalResidual <: DebugAction

A Debug action to print the primal residual.
The constructor accepts a printing function and some (shared) storage, which
should at least record `:x`, `:ξ` and `:n`.
"""
mutable struct DebugPrimalResidual <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugPrimalResidual(
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)), io::IO=stdout
    )
        return new(io, "Primal Residual: ", a)
    end
    function DebugPrimalResidual(
        values::Tuple{P,T,Q},
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)),
        io::IO=stdout,
    ) where {P,T,Q}
        update_storage!(a, Dict(k => v for (k, v) in zip((:x, :ξ, :n), values)))
        return new(io, "Primal Residual: ", a)
    end
end
function (d::DebugPrimalResidual)(
    p::P, o::ChambollePockOptions, i::Int
) where {P<:PrimalDualProblem}
    if all(has_storage.(Ref(d.storage), [:x, :ξ, :n])) && i > 0 # all values stored
        xOld, ξOld, nOld = get_storage.(Ref(d.storage), [:x, :ξ, :n]) #fetch
        print(d.io, d.prefix * string(primal_residual(p, o, xOld, ξOld, nOld)))
    end
    return d.storage(p, o, i)
end
@doc raw"""
    DebugPrimalDualResidual <: DebugAction

A Debug action to print the primaldual residual.
The constructor accepts a printing function and some (shared) storage, which
should at least record `:x`, `:ξ` and `:n`.
"""
mutable struct DebugPrimalDualResidual <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugPrimalDualResidual(
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)), io::IO=stdout
    )
        return new(io, "PD Residual: ", a)
    end
    function DebugPrimalDualResidual(
        values::Tuple{P,T,Q},
        a::StoreOptionsAction=StoreOptionsAction((:x, :ξ, :n)),
        io::IO=stdout,
    ) where {P,Q,T}
        update_storage!(a, Dict(k => v for (k, v) in zip((:x, :ξ, :n), values)))
        return new(io, "PD Residual: ", a)
    end
end
function (d::DebugPrimalDualResidual)(
    p::P, o::ChambollePockOptions, i::Int
) where {P<:PrimalDualProblem}
    if all(has_storage.(Ref(d.storage), [:x, :ξ, :n])) && i > 0 # all values stored
        xOld, ξOld, nOld = get_storage.(Ref(d.storage), [:x, :ξ, :n]) #fetch
        print(
            d.io,
            d.prefix * string(
                (
                    primal_residual(p, o, xOld, ξOld, nOld) +
                    dual_residual(p, o, xOld, ξOld, nOld)
                ) / manifold_dimension(p.M),
            ),
        )
    end
    return d.storage(p, o, i)
end

#
# Debugs
#
"""
    DebugPrimalChange(opts...)

Print the change of the primal variable by using [`DebugChange`](@ref),
see their constructors for detail.
"""
function DebugPrimalChange(
    a::StoreOptionsAction=StoreOptionsAction((:x,)), prefix="Primal Change: ", io::IO=stdout
)
    return DebugChange(a, prefix, io)
end

"""
    DebugPrimalIterate(opts...)

Print the change of the primal variable by using [`DebugIterate`](@ref),
see their constructors for detail.
"""
DebugPrimalIterate(opts...) = DebugIterate(opts...)

"""
    DebugDualIterate(e)

Print the dual variable by using [`DebugEntry`](@ref),
see their constructors for detail.
This method is further set display `o.ξ`.
"""
DebugDualIterate(opts...) = DebugEntry(:ξ, "ξ:", opts...)

"""
    DebugDualChange(opts...)

Print the change of the dual variable, similar to [`DebugChange`](@ref),
see their constructors for detail, but with a different calculation of the change,
since the dual variable lives in (possibly different) tangent spaces.
"""
mutable struct DebugDualChange <: DebugAction
    io::IO
    prefix::String
    storage::StoreOptionsAction
    function DebugDualChange(
        a::StoreOptionsAction=StoreOptionsAction((:ξ, :n)), io::IO=stdout
    )
        return new(io, "Dual Change: ", a)
    end
    function DebugDualChange(
        values::Tuple{T,P},
        a::StoreOptionsAction=StoreOptionsAction((:ξ, :n)),
        io::IO=stdout,
    ) where {P,T}
        update_storage!(a, Dict{Symbol,Any}(k => v for (k, v) in zip((:ξ, :n), values)))
        return new(io, "Dual Change: ", a)
    end
end
function (d::DebugDualChange)(
    p::P, o::ChambollePockOptions, i::Int
) where {P<:PrimalDualProblem}
    if all(has_storage.(Ref(d.storage), [:ξ, :n])) && i > 0 # all values stored
        ξOld, nOld = get_storage.(Ref(d.storage), [:ξ, :n]) #fetch
        print(
            d.io,
            d.prefix * string(
                norm(
                    p.N,
                    o.n,
                    vector_transport_to(p.N, nOld, ξOld, o.n, o.vector_transport_method) -
                    o.ξ,
                ),
            ),
        )
    end
    return d.storage(p, o, i)
end

"""
    DebugDualBaseIterate(io::IO=stdout)

Print the dual base variable by using [`DebugEntry`](@ref),
see their constructors for detail.
This method is further set display `o.n`.
"""
DebugDualBaseIterate(io::IO=stdout) = DebugEntry(:n, "n:", io)
"""
    DebugDualChange(a=StoreOptionsAction((:ξ)),io::IO=stdout)

Print the change of the dual base variable by using [`DebugEntryChange`](@ref),
see their constructors for detail, on `o.n`.
"""
function DebugDualBaseChange(a::StoreOptionsAction=StoreOptionsAction((:n)), io::IO=stdout)
    return DebugEntryChange(
        :n, (p, o, x, y) -> distance(p.N, x, y), a, "Dual Base Change:", io
    )
end

"""
    DebugPrimalBaseIterate(io::IO=stdout)

Print the primal base variable by using [`DebugEntry`](@ref),
see their constructors for detail.
This method is further set display `o.m`.
"""
DebugPrimalBaseIterate(io::IO=stdout) = DebugEntry(:m, "m:", io)

"""
    DebugPrimalBaseChange(a::StoreOptionsAction=StoreOptionsAction((:m)),io::IO=stdout)

Print the change of the primal base variable by using [`DebugEntryChange`](@ref),
see their constructors for detail, on `o.n`.
"""
function DebugPrimalBaseChange(
    a::StoreOptionsAction=StoreOptionsAction((:m)), io::IO=stdout
)
    return DebugEntryChange(
        :m, (p, o, x, y) -> distance(p.M, x, y), a, "Primal Base Change:", io
    )
end

#
# Records
#

# Primals are just the entries
"""
    RecordPrimalChange(a)

Create an [`RecordAction`](@ref) that records the primal value change,
i.e. [`RecordChange`](@ref), since we just redord the change of `o.x`.
"""
RecordPrimalChange() = RecordChange()

"""
    RecordDualBaseIterate(x)

Create an [`RecordAction`](@ref) that records the dual base point,
i.e. [`RecordIterate`](@ref), i.e. `o.x`.
"""
RecordPrimalIterate(x) = RecordIterate(x)

"""
    RecordDualIterate(ξ)

Create an [`RecordAction`](@ref) that records the dual base point,
i.e. [`RecordEntry`](@ref) of `o.ξ`, so .
"""
RecordDualIterate(ξ) = RecordEntry(ξ, :ξ)

"""
    RecordDualChange()

Create the action either with a given (shared) Storage, which can be set to the
`values` Tuple, if that is provided).
"""
RecordDualChange() = RecordEntryChange(:ξ, (p, o, x, y) -> distance(p.N, x, y))

"""
    RecordDualBaseIterate(n)

Create an [`RecordAction`](@ref) that records the dual base point,
i.e. [`RecordEntry`](@ref) of `o.n`.
"""
RecordDualBaseIterate(n) = RecordEntry(n, :n)

"""
    RecordDualBaseChange(e)

Create an [`RecordAction`](@ref) that records the dual base point change,
i.e. [`RecordEntryChange`](@ref) of `o.n` with distance to the last value to store a value.
"""
RecordDualBaseChange() = RecordEntryChange(:n, (p, o, x, y) -> distance(p.N, x, y))

"""
    RecordPrimalBaseIterate(x)

Create an [`RecordAction`](@ref) that records the primal base point,
i.e. [`RecordEntry`](@ref) of `o.m`.
"""
RecordPrimalBaseIterate(m) = RecordEntry(m, :m)
"""
    RecordPrimalBaseChange()

Create an [`RecordAction`](@ref) that records the primal base point change,
i.e. [`RecordEntryChange`](@ref) of `o.m` with distance to the last value to store a value.
"""
RecordPrimalBaseChange() = RecordEntryChange(:m, (p, o, x, y) -> distance(p.M, x, y))


text:
using NonstationaryProcesses
using Plots

tom = thomasCyclicallySymmetricSim(
    X0 = [0.1, 0.0, 0.0],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (0.18,),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.01,
    tmax = 2000.0,
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-12))

𝔱𝔬𝔪 = corruptphase(tom, ramp, (0.0, 1.0, 0.0, 2000.0))()

a = animatespectrum(𝔱𝔬𝔪, downsample=100, trail=20000, nperseg=2000)
gif(a, fps=24)


text:
abstract type Shape end

struct Circle <: Shape
    r::Real
end

struct Rectangle <: Shape
    l::Real
    w::Real
end

struct Square <: Shape
    l::Real
end

area(a::Circle) = π*a.r^2
area(a::Rectangle) = a.l * a.w
area(a::Square) = a.l^2

function main()

    # numbers_example()
    shapes_example()

end

function numbers_example()

    c = 1 + 2im
    r = 5

    println(c + c)
    println(c + r)
    @which c + c

end

function shapes_example()

    c = Circle(3)
    r = Rectangle(4, 5)
    s = Square(5)

    println(area(c))
    println(area(r))
    println(area(s))

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
###
  Name     : c7_22_optim.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 3/15/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using Optim
function g!(s, x)
 s[1] = -2.0*(1.0-x[1])-400.0*(x[2]-x[1]^2)*x[1]
 s[2] = 200.0*(x[2]-x[1]^2)
end
lower = [1.25, -2.1]
upper = [Inf, Inf]
initial_x = [2.0, 2.0]
od = OnceDifferentiable(f, g!, initial_x)
results = optimize(od, initial_x,lower,upper,Fminbox{GradientDescent}())


text:
function convert_arguments(LS::Type{<:LineSegments},
                           rg::RectilinearGrid{X,Y,0,T}) where {X,Y,T}
    lines = Vector{LineSegment{Point{2,T}}}(undef, X+Y)
    ylower = ymin(rg)
    yupper = ymax(rg)
    xlower = xmin(rg)
    xupper = xmax(rg)
    for i = 1:X
        lines[i] = LineSegment(Point(rg.x[i], ylower), Point(rg.x[i], yupper))
    end
    for i = 1:Y
        lines[X+i] = LineSegment(Point(xlower, rg.y[i]), Point(xupper, rg.y[i]))
    end
    return convert_arguments(LS, lines)
end

# 3D grid plot quadrilaterals
function convert_arguments(LS::Type{<:LineSegments},
                           rg::RectilinearGrid{X,Y,Z,T}) where {X,Y,Z,T}
    lines = Vector{LineSegment{Point{3,T}}}(undef, Z*(X+Y) + Y*(X+Z))

    xlower = xmin(rg)
    xupper = xmax(rg)
    ylower = ymin(rg)
    yupper = ymax(rg)
    zlower = zmin(rg)
    zupper = zmax(rg)
    ictr = 1
    # XY planes
    for iz = 1:Z
        for i = 1:X
            lines[ictr] = LineSegment(Point(rg.x[i], ylower, rg.z[iz]),
                                      Point(rg.x[i], yupper, rg.z[iz]))
            ictr += 1
        end
        for i = 1:Y
            lines[ictr] = LineSegment(Point(xlower, rg.y[i], rg.z[iz]),
                                      Point(xupper, rg.y[i], rg.z[iz]))
            ictr += 1
        end
    end
    # XZ planes
    for iy = 1:Y
        for i = 1:X
            lines[ictr] = LineSegment(Point(rg.x[i], rg.y[iy], zlower),
                                      Point(rg.x[i], rg.y[iy], zupper))
            ictr += 1
        end
        for i = 1:Z
            lines[ictr] = LineSegment(Point(xlower, rg.y[iy], rg.z[i]),
                                      Point(xupper, rg.y[iy], rg.z[i]))
            ictr += 1
        end
    end

    return convert_arguments(LS, lines)
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
Constraint handler solving the GBD subproblem (fixed z).
"""
mutable struct GBDSubHdlr <: SCIP.AbstractConstraintHandler
    scip::SCIP.Optimizer # of master problem
    subsolver            # for subproblems
    finaltime::Float64 # seconds
    addcritpath::Bool
    addnogoods::Bool
    counter::Int
    primal_bound::Float64
    best_solution::Union{CandSol, Nothing}

    function GBDSubHdlr(scip, subsolver, finaltime, addcritpaths, addnogoods)
        return new(scip, subsolver, finaltime, addcritpaths, addnogoods,
                   0, Inf, nothing)
    end
end

"""
Constraint data, referencing problem data and master variables.
"""
mutable struct GBDSubCons <: SCIP.AbstractConstraint{GBDSubHdlr}
    inst::Instance
    topo::Topology
    master::Master
end

function solve_gbd_sub(ch::GBDSubHdlr, cons::GBDSubCons,
                       sol::Ptr{SCIP.SCIP_SOL}=C_NULL; enforce=true)
    ch.counter += 1

    # fetch solution values
    ysol = SCIP.sol_values(ch.scip, cons.master.y, sol)
    zsol = SCIP.sol_values(ch.scip, cons.master.z, sol)
    qsol = SCIP.sol_values(ch.scip, cons.master.q, sol)
    ϕsol = SCIP.sol_values(ch.scip, cons.master.ϕ, sol)

    cand = CandSol(zsol .>= 0.5, qsol, ϕsol)

    # solve subproblem (from scratch, no warmstart)
    submodel, π, Δl, Δu, ploss, plb, pub =
        make_sub(cons.inst, cons.topo, cand, ch.subsolver)
    # TODO: solver.writemodels && writeLP(submodel, "sub_$(counter).lp", genericnames=false)
    settimelimit!(submodel, ch.subsolver, ch.finaltime - time())
    JuMP.optimize!(submodel)
    substatus = JuMP.termination_status(submodel)
    @assert substatus == MOI.OPTIMAL "Slack model is always feasible"
    totalslack = JuMP.objective_value(submodel)

    if totalslack ≈ 0.0
        # maybe only the relaxation is feasible, we have to check also the
        # "exact" subproblem with equations constraints.
        submodel2, _ = make_sub(cons.inst, cons.topo, cand, ch.subsolver,
                                relaxed=false)
        # TODO: solver.writemodels && writeLP(submodel2, "sub_exact_iter$(iter).lp", genericnames=false)
        settimelimit!(submodel2, ch.subsolver, ch.finaltime - time())
        JuMP.optimize!(submodel2)
        substatus2 = JuMP.termination_status(submodel2)
        @assert substatus2 == MOI.OPTIMAL "Slack model is always feasible"
        totalslack2 = JuMP.objective_value(submodel2)

        if totalslack2 ≈ 0.0
            # check whether this candidate is actually an improvement
            L = pipelengths(cons.topo)
            A = cons.topo.arcs
            D = cons.inst.diameters
            c = [diam.cost for diam in D]
            obj = [c[i] * L[a] for a=1:length(A), i=1:length(D)]
            candprimal = sum(obj .* zsol)
            if candprimal < ch.primal_bound
                @debug "  found improving solution :-)"
                ch.best_solution = cand
                ch.primal_bound = candprimal
            else
                @debug "  cand's obj val $(candprimal) no improvement over $(ch.primal_bound)!"
            end

            return SCIP.SCIP_FEASIBLE
        else
            @debug "  subproblem/relaxation gap!"
            if !enforce
                return SCIP.SCIP_INFEASIBLE
            end

            # cut off candidate with no-good on z
            active = zsol .> 0.5
            coefs = 2.0 * active .- 1.0
            @cb_constraint(ch.scip, coefs ⋅ cons.master.z ≤ sum(active) - 1)

            return SCIP.SCIP_CONSADDED
        end
    else # totalslack > 0
        if !enforce
            return SCIP.SCIP_INFEASIBLE
        end

        # generate cuts and add to master
        if ch.addnogoods
            # cut off candidate with no-good on z
            active = zsol .> 0.5
            coefs = 2.0 * active .- 1.0
            @cb_constraint(ch.scip, coefs ⋅ cons.master.z ≤ sum(active) - 1)
        end

        if ch.addcritpath
            dualsol = SubDualSol(
                JuMP.dual.(ploss), JuMP.dual.(plb), JuMP.dual.(pub))
            critpathcuts(cons.inst, cons.topo, cons.master, cand,
                         subdualsol, ch.subsolver, cb=ch.scip)
        end

        @debug "  added some(?) cuts."
        return SCIP.SCIP_CONSADDED
    end
end

function SCIP.check(ch::GBDSubHdlr,
                    constraints::Array{Ptr{SCIP.SCIP_CONS}},
                    sol::Ptr{SCIP.SCIP_SOL},
                    checkintegrality::Bool,
                    checklprows::Bool,
                    printreason::Bool,
                    completely::Bool)
    @assert length(constraints) == 1
    cons::GBDSubCons = SCIP.user_constraint(constraints[1])
    return solve_gbd_sub(ch, cons, sol, enforce=false)
end

function SCIP.enforce_lp_sol(ch::GBDSubHdlr,
                             constraints::Array{Ptr{SCIP.SCIP_CONS}},
                             nusefulconss::Cint,
                             solinfeasible::Bool)
    @assert length(constraints) == 1
    cons::GBDSubCons = SCIP.user_constraint(constraints[1])
    return solve_gbd_sub(ch, cons, enforce=true)
end

function SCIP.enforce_pseudo_sol(ch::GBDSubHdlr,
                                 constraints::Array{Ptr{SCIP.SCIP_CONS}},
                                 nusefulconss::Cint,
                                 solinfeasible::Bool,
                                 objinfeasible::Bool)
    @assert length(constraints) == 1
    cons::GBDSubCons = SCIP.user_constraint(constraints[1])
    return solve_gbd_sub(ch, cons, enforce=true)
end

function SCIP.lock(ch::GBDSubHdlr,
                   constraint::Ptr{SCIP.SCIP_CONS},
                   locktype::SCIP.SCIP_LOCKTYPE,
                   nlockspos::Cint,
                   nlocksneg::Cint)
    cons::GBDSubCons = SCIP.user_constraint(constraint)
    master::Master = cons.master
    for v in [master.y; master.z[:]; master.q; master.ϕ]
        # TODO: understand why lock is called during SCIPfree, after the
        # constraint should have been deleted already. Does it mean we should
        # implement CONSTRANS?
        var_::Ptr{SCIP.SCIP_VAR} = SCIP.var(ch.scip, JuMP.index(v))
        var_ != C_NULL || continue  # avoid segfault!

        SCIP.@SC SCIP.SCIPaddVarLocksType(
            ch.scip, var_, locktype, nlockspos + nlocksneg, nlockspos + nlocksneg)
    end
end


text:
<<<<<<< HEAD
using GLFW
using GLTF
using ModernGL
using CSyntax

using Logging
using Dates

include("src/globVar.jl")
include("src/fileIO.jl")

using .errVar
using .glfwVar

include("src/exit.jl")
include("src/init.jl")

mainInit()
modelLoad("test.gltf")
#=
=======
include("src/exit.jl")
include("src/init.jl")

using GLFW
using ModernGL
using CSyntax

mainInit()

# hard-coded shaders
const vert_source = """
#version 410 core
in vec3 vp;
void main(void)
{
    gl_Position = vec4(vp, 1.0);
}"""
const frag_source = """
#version 410 core
out vec4 frag_colour;
void main(void)
{
    frag_colour = vec4(0.5, 0.0, 0.5, 1.0);
}"""

# compile shaders
vert_shader = glCreateShader(GL_VERTEX_SHADER)
glShaderSource(vert_shader, 1, Ptr{GLchar}[pointer(vert_source)], C_NULL)
glCompileShader(vert_shader)
frag_shader = glCreateShader(GL_FRAGMENT_SHADER)
glShaderSource(frag_shader, 1, Ptr{GLchar}[pointer(frag_source)], C_NULL)
glCompileShader(frag_shader)

# create and link shader program
program = glCreateProgram()
glAttachShader(program, vert_shader)
glAttachShader(program, frag_shader)
glLinkProgram(program)

# vertex data
points = GLfloat[ 0.0,  0.5, 0.0,
                  0.5, -0.5, 0.0,
                 -0.5, -0.5, 0.0]

# create buffers located in the memory of graphic card
vbo = GLuint(0)
@c glGenBuffers(1, &vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)

# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

>>>>>>> 5f9517fcfc507bc9c83bc0ca5755ab4a4d01f653
# render
while !GLFW.WindowShouldClose(win)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glUseProgram(program)
    glBindVertexArray(vao)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    # check and call events
    GLFW.PollEvents()
    # swap the buffers
    GLFW.SwapBuffers(win)
end

mainExit()
<<<<<<< HEAD
=#
=======
>>>>>>> 5f9517fcfc507bc9c83bc0ca5755ab4a4d01f653


text:
using LocalAnisotropies
using GeoStats
using Test
import LocalAnisotropies: rotmat

@testset "LocalAnisotropies.jl" begin
    # convert data to LocalAnisotropy
    dummy = georef((az=1:10, r1=1:10, r2=1:10), PointSet(rand(2,10)))
    pars  = localanisotropies(dummy, [:az], [:r1,:r2], :GSLIB)

    # convert between different rotation conventions
    angs1 = convertangles([30,30,30], :GSLIB, :Datamine)
    angs2 = convertangles(angs1, :Datamine, :GSLIB)
	rmat  = rotmat([1,1,1], [30,30,30], :GSLIB)
    @test all(rmat .≈ rotmat([1,1,1], angs1, :Datamine))
    @test all(rmat .≈ rotmat([1,1,1], angs2, :GSLIB))
    angs_ = convertangles.(pars.rotation, :GSLIB)

    # pars_ = convertpars(pars, convention=:Datamine)
    # exportpars("C:\\Users\\test.csv", pars, :GSLIB)

	# local anisotropies from pointset coordinates
	data  = rmat[1] * vcat(rand(2,100),zeros(1,100))
	pset  = PointSet(data)
	nhood = KNearestSearch(pset, 10)
    gpars = localanisotropies(Geometric, nhood, simplify=true)
    gpars = localanisotropies(Geometric, nhood, simplify=false)

    grid2d = (10,10)
    grid3d = (10,10,5)

    for dims in (grid2d, grid3d)
        # reference scenario for tests
		R = length(dims)
		D = if R==2
    		georef((P=[sin(i)+j for i in 1:dims[1], j in 1:dims[2]],))
    	else
    		georef((P=[sin(i)+j+k for i in 1:dims[1], j in 1:dims[2], k in 1:dims[3]],))
    	end

        n = round(Int, 0.2*prod(dims))
        S = sample(D, n, replace=false)
        S = georef(values(S), PointSet(centroid.(domain(S))))
        G = CartesianGrid(dims...)
        searcher = KNearestSearch(G, 10)

        # get local anisotropies
        lpars = localanisotropies(Gradients, D, :P, 3)

        # rescale magnitude and interpolate local anisotropies
        lpars = rescale_magnitude(lpars, r2=(0.2,1.0))
        lpars = reference_magnitude(lpars, :Y)
        lpars = smoothpars(lpars, searcher)

		spars = nnpars(lpars, D, S)

        # interpolate in a coarser grid
        G_ = if R==2
			CartesianGrid((5,5),(0.5,0.5),(2.0,2.0))
		else
			CartesianGrid((5,5,3),(0.5,0.5,0.5),(2.0,2.0,2.0))
		end
        lpars_ = idwpars(lpars, searcher, G_)

        # Estimation problem
        P = EstimationProblem(S, G, :P)
        γ = NuggetEffect(0.1) + 0.9*ExponentialVariogram(range=60.0)

        # LocalKriging (MW)
        MW = LocalKriging(:P => (variogram=γ, localaniso=lpars, method=:MovingWindows))
        s1 = solve(P, MW)

        # LocalKriging (KC)
        KC = LocalKriging(:P => (variogram=γ, localaniso=lpars, method=:KernelConvolution))
        s2 = solve(P, KC)

        # Spatial deformation: anisotropic distances
        Sd1, Dd1 = deformspace(S, G, lpars, AnisoDistance, anchors=1500)
        Pd1 = EstimationProblem(Sd1, Dd1, :P)
        s3 = solve(Pd1, Kriging(:P => (variogram=γ,)))
		x3 = to_3d(s3)

        # Spatial deformation: anisotropic variogram distances
        Sd2, Dd2 = deformspace(S, G, lpars, KernelVariogram, γ, anchors=1500)
        Pd2 = EstimationProblem(Sd2, Dd2, :P)
        s4 = solve(Pd2, Kriging(:P => (variogram=γ,)))
		x4 = to_3d(s4)

        # Spatial deformation: geodesic anisotropic distances
        LDa = graph(S, G, lpars, AnisoDistance, searcher)
        Sd3, Dd3 = deformspace(LDa, GraphDistance, anchors=1500)
        Pd3 = EstimationProblem(Sd3, Dd3, :P)
        s5 = solve(Pd3, Kriging(:P => (variogram=γ,)))
		x5 = to_3d(s5)

        # Spatial deformation: geodesic anisotropic variogram distances
        LDv = graph(S, G, lpars, KernelVariogram, γ, searcher)
        Sd4, Dd4 = deformspace(LDv, GraphDistance, anchors=1500)
        Pd4 = EstimationProblem(Sd4, Dd4, :P)
        s6 = solve(Pd4, Kriging(:P => (variogram=γ,)))
		x6 = to_3d(s6)
    end
end


text:
using CUDA, Lux, Random, Test

include("utils.jl")

@testset "Elementwise Operation Dispatches" begin
    rng = Random.default_rng()
    Random.seed!(rng, 0)
    custom_activation(x) = abs(x)

    for T in [Float64, Float32, ComplexF64, ComplexF32]
        x = randn(rng, T, 10, 5, 2)
        y = randn(rng, T, 10, 1, 2)

        # On CPU the fallback should always work
        @test Lux.elementwise_add(x, y) == x .+ y
        @test Lux.elementwise_mul(x, y) == x .* y
        @test Lux.applyactivation(tanh, x) == tanh.(x)
        @test Lux.applyactivation(custom_activation, x) == custom_activation.(x)

        if T <: Real
            # Gradient for complex outputs are not defined
            test_gradient_correctness_fdm(sum ∘ Lux.elementwise_add, x, y)
        end

        # On GPU try to use CUDNN
        if CUDA.functional()
            x_g = x |> gpu
            y_g = y |> gpu

            @test Lux.elementwise_add(x_g, y_g) == x_g .+ y_g
            @test Lux.elementwise_mul(x_g, y_g) == x_g .* y_g
            if T <: Real
                @test Lux.applyactivation(tanh, x_g) == tanh.(x_g)
            else
                ## See https://github.com/FluxML/NNlibCUDA.jl/issues/47
                @test_broken Lux.applyactivation(tanh, x_g) == tanh.(x_g)
            end
            # Custom Activation test
            @test Lux.applyactivation(custom_activation, x_g) == custom_activation.(x_g)
        end
    end
end


text:
using Test
using ColorForwardDiff
using ForwardDiff
using ForwardDiff: derivative, gradient
using ColorTypes

@testset "derivative" begin
    @test derivative(x->RGB(0.5x, x^2, sin(x)), 0.5) ≈ RGB{Float64}(0.5, 1.0, cos(0.5))
end

@testset "gradient" begin
    @test gradient(x->RGB(0.5 * x[1], 0.2, 0.3), [0.5])[1] ≈ RGB{Float64}(0.5, 0.0, 0.0)
    @test gradient(x->RGB(0.1, x[1]^2, 0.3), [0.5])[1] ≈ RGB{Float64}(0.0, 1.0, 0.0)
    @test gradient(x->RGB(0.1, 0.2, sin(x[1])), [0.5])[1] ≈ RGB{Float64}(0.0, 0.0, cos(0.5))

    rgb_xy = gradient(x->RGB(x[1], x[2]^2, 0.3), [0.5, 0.5])
    @test rgb_xy[1] ≈ RGB{Float64}(1.0, 0.0, 0.0)
    @test rgb_xy[2] ≈ RGB{Float64}(0.0, 1.0, 0.0)

    rgb_xyz = gradient(x->RGB(x[1] + x[3], x[2]^2 - x[1], x[2] * x[3]), [0.5, 0.5, 0.5])
    @test rgb_xyz[1] ≈ RGB{Float64}(1.0, -1.0, 0.0)
    @test rgb_xyz[2] ≈ RGB{Float64}(0.0, 1.0, 0.5)
    @test rgb_xyz[3] ≈ RGB{Float64}(1.0, 0.0, 0.5)
end

@testset "traits" begin
    @test gradient(x->red(RGB(0.5 * x[1], 0.2, 0.3)) + x[1], [0.5])[1] ≈ 1.5
end


text:
# ----------------------------------------------------------------
# Radiative forcing from methane (Etminan et al. 2016 equations).
# ----------------------------------------------------------------

@defcomp rf_ch4_etminan begin

    a₃         = Parameter()             # Methane forcing constant.
    b₃         = Parameter()             # Methane forcing constant.
    N₂O_0      = Parameter()             # Initial (pre-industrial) nitrous oxide concentration (ppb).
    scale_CH₄  = Parameter()             # Scaling factor for uncertainty in methane radiative forcing.
    CH₄_0      = Parameter()             # Initial (pre-industrial) methane concentration (ppb).
    N₂O        = Parameter(index=[time]) # Atmospheric nitrous oxide concentration (ppb).
    CH₄        = Parameter(index=[time]) # Atmospheric methane concentration (ppb).

    rf_CH₄     = Variable(index=[time])  # Direct forcing from atmospheric methane concentrations (Wm⁻²).

    function run_timestep(p, v, d, t)

        if is_first(t)
            # Set initial radiative forcing to 0.
            v.rf_CH₄[t] = 0.0

        else
            # Create temporary averaging variables.
            M̄ = 0.5 * (p.CH₄[t] + p.CH₄_0)
            N̄ = 0.5 * (p.N₂O[t] + p.N₂O_0)

            # Direct methane radiative forcing.
            v.rf_CH₄[t] = ((p.a₃ * M̄ + p.b₃ * N̄ + 0.043) * (sqrt(p.CH₄[t]) - sqrt(p.CH₄_0))) * p.scale_CH₄
        end
    end
end


text:
using RungeKutta.Tableaus: get_radau_1_nodes, get_radau_1_weights, get_radau_1_coefficients,
                           get_radau_2_nodes, get_radau_2_weights, get_radau_2_coefficients

@testset "$(rpad("Radau Tableaus",80))" begin

    @test_throws ErrorException get_radau_1_nodes(1)
    @test_throws ErrorException get_radau_1_weights(1)
    @test_throws ErrorException get_radau_1_coefficients(1)

    @test_throws ErrorException get_radau_2_nodes(1)
    @test_throws ErrorException get_radau_2_weights(1)
    @test_throws ErrorException get_radau_2_coefficients(1)

    function _TableauRadauIA2(T=Float64)
        a = [[1//4  -1//4  ]
             [1//4   5//12 ]]
        b = [1//4, 3//4]
        c = [0,    2//3]

        Tableau{T}(:RadauIA2, 3, a, b, c)
    end

    function _TableauRadauIA3(T=Float64)
        a = [
                [ 1//9     (- 1 -    √6)/18    (- 1 +    √6)/18  ]
                [ 1//9     ( 88 +  7*√6)/360   ( 88 - 43*√6)/360 ]
                [ 1//9     ( 88 + 43*√6)/360   ( 88 -  7*√6)/360 ]
            ]
        b = [1/9,         (16+√6)/36,  (16-√6)/36 ]
        c = [0,           ( 6-√6)/10,  ( 6+√6)/10 ]

        Tableau{T}(:RadauIA3, 5, a, b, c)
    end

    function _TableauRadauIB2(T=Float64)
        a = [[1//8  -1//8  ]
             [7//24  3//8  ]]
        b = [1//4, 3//4]
        c = [0,    2//3]

        Tableau{T}(:RadauIB2, 3, a, b, c)
    end

    function _TableauRadauIB3(T=Float64)
        a = [
                [            1//18    (- 1 -     √6)/36    (- 1 +     √6)/36   ]
                [ ( 52 +  3*√6)/450   ( 16 +     √6)/72    (472 - 217*√6)/1800 ]
                [ ( 52 -  3*√6)/450   (472 + 217*√6)/1800  ( 16 -     √6)/72   ]
            ]
        b = [1/9,         (16+√6)/36,  (16-√6)/36 ]
        c = [0,           ( 6-√6)/10,  ( 6+√6)/10 ]

        Tableau{T}(:RadauIB3, 5, a, b, c)
    end

    function _TableauRadauIIA2(T=Float64)
        a = [[5//12  -1//12]
             [3//4    1//4 ]]
        b = [3//4, 1//4]
        c = [1//3, 1//1]

        Tableau{T}(:RadauIIA2, 3, a, b, c)
    end

    function _TableauRadauIIA3(T=Float64)
        a = [
                [ 11/45 -  7*√6/360    37/225-169*√6/1800   -2/225+√6/75 ]
                [ 37/225+169*√6/1800   11/45 +  7*√6/360    -2/225-√6/75 ]
                [  4/9  -    √6/36      4/9  +    √6/36      1/9         ]
            ]
        b = [4/9-√6/36,   4/9+√6/36,  1/9 ]
        c = [2/5-√6/10,   2/5+√6/10,  1   ]

        Tableau{T}(:RadauIIA3, 5, a, b, c)
    end

    function _TableauRadauIIB2(T=Float64)
        a = [[3//8   -1//24]
             [7//8    1//8 ]]
        b = [3//4, 1//4]
        c = [1//3, 1//1]

        Tableau{T}(:RadauIIB2, 3, a, b, c)
    end

    function _TableauRadauIIB3(T=Float64)
        a = [
                [( 16 -     √6)/72    (328 - 167*√6)/1800   (-2 + 3*√6)/450 ]
                [(328 + 167*√6)/1800  ( 16 +     √6)/72     (-2 - 3*√6)/450 ]
                [( 85 -  10*√6)/180   ( 85 +  10*√6)/180              1/18  ]
            ]
        b = [4/9-√6/36,   4/9+√6/36,  1/9 ]
        c = [2/5-√6/10,   2/5+√6/10,  1   ]

        Tableau{T}(:RadauIIB3, 5, a, b, c)
    end

    @test_throws ErrorException TableauRadauIA(1)
    @test_throws ErrorException TableauRadauIIA(1)

    @test TableauRadauIA(2)  ≈ _TableauRadauIA2()
    @test TableauRadauIA(3)  ≈ _TableauRadauIA3()

    @test TableauRadauIB(2)  ≈ _TableauRadauIB2()
    @test TableauRadauIB(3)  ≈ _TableauRadauIB3()

    @test TableauRadauIIA(2) ≈ _TableauRadauIIA2()
    @test TableauRadauIIA(3) ≈ _TableauRadauIIA3()

    @test TableauRadauIIB(2) ≈ _TableauRadauIIB2()
    @test TableauRadauIIB(3) ≈ _TableauRadauIIB3()

end


text:
import RowEchelon
import LinearAlgebra

#=
Homework 7
Section 3.1
Problems 13
=#

A = [ 4 0 -7 3 -5 ;
     0 0 2 0 0;
     7 3 -6 4 -8;
     5 0 5 2 -3;
     0 0 9 -1 2]

detA = LinearAlgebra.det(A)

println("Original Matrix:")
println(A)
println("Determinate:")
println(detA)


text:
using DynamicGrids, Test

# Reducing functions may return the original grid,
# not a copy, so we need to be careful with them.

@testset "Reducing functions over NamedTuple" begin
    rule = Cell{:a}() do data, state, I
        state + 10
    end
    @testset "single named grid" begin
        init = (a=[1 3],)
        transformed_output = TransformedOutput(sum, init; tspan=1:3)
        @test length(transformed_output) == 3
        @test transformed_output[1] == [1 3]
        @test transformed_output[2] == [0 0]
        @test transformed_output[3] == [0 0]

        sim!(transformed_output, rule)
        @test transformed_output[1] == [1 3]
        @test transformed_output[2] == [11 13]
        @test transformed_output[3] == [21 23]
    end

    @testset "multiple named grids" begin
        init = (a=[1 3], b=[5 5],)
        transformed_output = TransformedOutput(sum, init; tspan=1:3)
        @test length(transformed_output) == 3
        @test transformed_output[1] == [6 8]
        @test transformed_output[2] == [0 0]
        @test transformed_output[3] == [0 0]

        sim!(transformed_output, rule)
        @test transformed_output[1] == [6 8]
        @test transformed_output[2] == [16 18]
        @test transformed_output[3] == [26 28]
    end
end

@testset "Reducing functions over Array" begin
    rule = Cell() do data, state, I
        state + 10.0
    end
    init = [1 3]
    transformed_output = TransformedOutput(sum, init; tspan=1:3)
    @test length(transformed_output) == 3
    @test transformed_output[1] == 4
    @test transformed_output[2] == 0
    @test transformed_output[3] == 0

    sim!(transformed_output, rule)
    @test transformed_output[1] == 4
    @test transformed_output[2] == 24
    @test transformed_output[3] == 44
end

@testset "Padded grids" begin
    init =  Bool[
             0 0 0 0 0 0 0
             0 0 0 0 1 1 1
             0 0 0 0 0 0 1
             0 0 0 0 0 1 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ]
    test2 = Bool[
             0 0 0 0 0 1 0
             0 0 0 0 0 1 1
             0 0 0 0 1 0 1
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ]
    test3 = Bool[
             0 0 0 0 0 1 1
             0 0 0 0 1 0 1
             0 0 0 0 0 0 1
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
             0 0 0 0 0 0 0
            ]

    # Sum the last column in the grid
    output = TransformedOutput(A -> sum(view(A, :, 7)), init; tspan=1:3)
    sim!(output, Life())
    output == [2, 2, 3]

    # Sum the first row in the :a grid of the NamedTuple
    output = TransformedOutput(gs -> sum(view(gs[:a], 1, :)), (a=init,); tspan=1:3)
    sim!(output, Life{:a}())
    output == [0, 1, 2]
end

@testset "Defining new function" begin
    ruleAR = Cell() do data, state, I
        state + 10.0
    end
    ruleNT = Cell{:a, :a}() do data, state, I
        state + 10.0
    end

    # Defining 2 functions giving square of all single values
    f_square(AR::AbstractArray) = map(x -> x^2, AR)
    function f_square(NT::NamedTuple)
        v = [map(x -> x^2, el) for el in values(NT)]
        return NamedTuple{keys(NT)}(v)
    end

    @testset "Array" begin
        init = [1 3]
        f_square(init) == [1 9]

        transformed_output = TransformedOutput(f_square, init; tspan=1:3)
        sim!(transformed_output, ruleAR)

        @test transformed_output[1] == [1 3] .^2
        @test transformed_output[2] == ([1 3] .+ 10) .^2
        @test transformed_output[3] == ([1 3] .+ 10 .+ 10) .^2
    end

    @testset "Single NamedTuple" begin
        init = (a = [1 3],)
        f_square(init) == (a = [ 1 9],)

        transformed_output = TransformedOutput(f_square, init; tspan=1:3)
        sim!(transformed_output, ruleNT)

        @test transformed_output[1] == (a=[1 3] .^2 ,)
        @test transformed_output[2] == (a=([1 3] .+ 10) .^2 ,)
        @test transformed_output[3] == (a=([1 3] .+ 10 .+ 10) .^2 ,)

    end

    @testset "Multiple NamedTuple" begin
        init = (a = [1 3], b = [2 5])
        @test f_square(init) == (a = [1 9], b = [4 25])

        transformed_output = TransformedOutput(f_square, init; tspan=1:3)
        sim!(transformed_output, ruleNT)

        @test transformed_output[1] == (a=[1 3] .^2 , b=[2 5] .^2)
        @test transformed_output[2] == (a=([1 3] .+ 10) .^2 , b=[2 5] .^2)
        @test transformed_output[3] == (a=([1 3] .+ 10 .+ 10) .^2 , b=[2 5] .^2)

    end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SLICPartition(k, m; tol=1e-4, maxiter=10, vars=nothing)

A method for partitioning spatial data into approximately `k`
clusters using Simple Linear Iterative Clustering (SLIC).
The method produces clusters of samples that are spatially
connected based on a distance `dₛ` and that, at the same
time, are similar in terms of `vars` with distance `dᵥ`.
The tradeoff is controlled with a hyperparameter parameter
`m` in an additive model `dₜ = √(dᵥ² + m²(dₛ/s)²)`.

## Parameters

* `k`       - Approximate number of clusters
* `m`       - Hyperparameter of SLIC model
* `tol`     - Tolerance of k-means algorithm (default to `1e-4`)
* `maxiter` - Maximum number of iterations (default to `10`)
* `vars`    - Variables (or features) to consider (default to all)

## References

* Achanta et al. 2011. [SLIC superpixels compared to state-of-the-art
  superpixel methods](https://ieeexplore.ieee.org/document/6205760)
"""
struct SLICPartition <: PartitionMethod
  k::Int
  m::Float64
  tol::Float64
  maxiter::Int
  vars::Union{Vector{Symbol},Nothing}
end

SLICPartition(k::Int, m::Real; tol=1e-4, maxiter=10, vars=nothing) =
  SLICPartition(k, m, tol, maxiter, vars)

function partition(sdata, method::SLICPartition)
  # variables used for clustering
  datavars = collect(name.(variables(sdata)))
  vars = isnothing(method.vars) ? datavars : method.vars

  @assert vars ⊆ datavars "SLIC features not found in spatial data"

  # SLIC hyperparameter
  m = method.m

  # initial spacing of clusters
  s = slic_spacing(sdata, method)

  # initialize cluster centers
  c = slic_initialization(sdata, s)

  # ball neighborhood search
  searcher = NeighborhoodSearch(sdata, BallNeighborhood(s))

  # pre-allocate memory for label and distance
  l = fill(0, nelms(sdata))
  d = fill(Inf, nelms(sdata))

  # performance parameters
  tol     = method.tol
  maxiter = method.maxiter

  # k-means algorithm
  err, iter = Inf, 0
  while err > tol && iter < maxiter
    o = copy(c)

    slic_assignment!(sdata, searcher, vars, m, s, c, l, d)
    slic_update!(sdata, c, l)

    err = norm(c - o) / norm(o)
    iter += 1
  end

  subsets = [findall(isequal(k), l) for k in 1:length(c)]

  SpatialPartition(sdata, subsets)
end

function slic_spacing(sdata, method)
  V = volume(boundbox(sdata))
  d = ncoords(sdata)
  k = method.k
  (V/k) ^ (1/d)
end

function slic_initialization(sdata, s)
  # efficient neighbor search
  searcher = KNearestSearch(sdata, 1)

  # bounding box properties
  bbox = boundbox(sdata)
  lo, up = extrema(bbox)

  # cluster centers
  clusters = Vector{Int}()
  neighbor = Vector{Int}(undef, 1)
  ranges = [(l+s/2):s:u for (l, u) in zip(lo, up)]
  for x in Iterators.product(ranges...)
    search!(neighbor, SVector(x), searcher)
    push!(clusters, neighbor[1])
  end

  unique(clusters)
end

function slic_assignment!(sdata, searcher, vars, m, s, c, l, d)
  for (k, cₖ) in enumerate(c)
    xₖ = coordinates(sdata, [cₖ])
    inds = search(vec(xₖ), searcher)

    # distance between coordinates
    X  = coordinates(sdata, inds)
    dₛ = pairwise(Euclidean(), X, xₖ, dims=2)

    # distance between variables
    𝒮ᵢ = view(sdata, inds, vars)
    𝒮ₖ = view(sdata, [cₖ], vars)
    V  = Tables.matrix(values(𝒮ᵢ))
    vₖ = Tables.matrix(values(𝒮ₖ))
    dᵥ = pairwise(Euclidean(), V, vₖ, dims=1)

    # total distance
    dₜ = @. √(dᵥ^2 + m^2 * (dₛ/s)^2)

    @inbounds for (i, ind) in enumerate(inds)
      if dₜ[i] < d[ind]
        d[ind] = dₜ[i]
        l[ind] = k
      end
    end
  end
end

function slic_update!(sdata, c, l)
  for k in 1:length(c)
    inds = findall(isequal(k), l)
    X  = coordinates(sdata, inds)
    μ  = mean(X, dims=2)
    dₛ = pairwise(Euclidean(), X, μ, dims=2)
    @inbounds c[k] = inds[argmin(vec(dₛ))]
  end
end


text:
using Plots; pyplot(xtickfontsize=13, ytickfontsize=13, xguidefontsize=16, yguidefontsize=16, legendfontsize=12)
default(palette = palette(:okabe_ito))
using DelimitedFiles
using Printf
#using StatsFuns, SpecialFunctions
using nevanlinna_ac

"""
function spectral_density(J::Real, ω::Real, β::Real; η::Float64=0.05,Γ::Real=1.)
    # 2Imχ(ω)
    if J/Γ == 0.
        return 2π*sinh(Γ*β)/cosh(Γ*β)*(delta(ω-2Γ,η) - delta(ω+2Γ,η))
    elseif J/Γ == 1.
        # Ref: PRX.4.031008(2014) A6
        g0 = 0.858714569; zc = J^(-1/4)
        T = 1/β
        up = zc * g0 * β^(3/4)
        down = 2^(1/4) * √π * gamma(1/8) * gamma(5/8)
        fac = up/down
        res = sinh(ω/(2*T)) * abs(gamma(1/8 - 1im*ω/(2π*T)))^2
        return 2*fac*res
    else
        @error "No exact results. J/Γ should be 1. or 0."
    end
end
"""

function chi_div_w(w::Vector, path::String, n::Int64)
	η = 0.001
	g = make_input(path)
	gs = MasubaraGF(n, g.GF[length(g.GF)-n+1:end])
	return [spectrum_density(i,η,gs) for i in w]
end

#setups
#o = [-i for i in range(-2π,2π,step = 4π/400)]
o = [-i for i in range(-1,1,step = 0.01)]
#brange = [10.0, 20.0, 30.0, 40.0]
#de = [spectral_density(1.0, i, 20) for i in o]

g = 1.0
β = 20

#load Zi-Long Li data
#pathl = @sprintf "./data/TFIsing-Li/Aw_%.1f_b20.txt" g
#dl = readdlm(pathl)

#n = 5 #number of Masubara frequencies used
nrange = [6, 7, 8]
chi = Array{Vector}(undef,3,2)
for  i = 1:3
    path1 = @sprintf "./data/gt/gt8_g_%.1f_b_%i.txt" g β
    path2 = @sprintf "./data/gt/gt28_g_%.1f_b_%i.txt" g β
    chi[i,1] = chi_div_w(o, path1, nrange[i])
    chi[i,2] = chi_div_w(o, path2, nrange[i])
end

#sc = 0.01 .* sum.(chi)
#for i = 1:4
#    s = @sprintf "%.5f" sc[i]
#    println(s)
#end

#v = 2 * abs(g-1)
plot()
#plot!(o, de ./ o, line=(:black,2),label="theory")
#plot!([-v,v], seriestype="vline",line=(:dash, :black), label=false)
#plot!(dl[:,1], dl[:,2] ./ dl[:,1], line=(:black,2), label="Li")
for j=1:3
    plot!(o, chi[j,1], line=(2), label=@sprintf "D=8,n=%i" nrange[j])
end
for j=1:3
    plot!(o, chi[j,2], line=(:dash, 1),
    marker=(:circle, 2, stroke(0)),
    label=@sprintf "D=2×8,n=%i" nrange[j])
end

plot!(xlim=(-1, 1))
plot!(title = @sprintf "g = %.1f, β=20, η = 0.001" g)
plot!(xlabel="ω", ylabel="χ''(ω)/ω")
spath = @sprintf "./notes/chi_g_%.1f_b_20_dchange.pdf" g
savefig(spath)


text:
import FLOWFarm; const ff = FLOWFarm
using DelimitedFiles
using CSV
using PyPlot

# based on IEA case study 3

# set initial turbine x and y locations
layout_file_name = "./inputfiles/shiloh_layout.csv"
layout_data = CSV.read(layout_file_name)

# println(layout_data)
# println(layout_data.elevation[10])

nturbines = length(layout_data.elevation)
lat = zeros(nturbines)
long = zeros(nturbines)
hub_height = zeros(nturbines)
model = layout_data.model
global m
m = 0

for i in 1:nturbines
    lat[i] = layout_data.latitude[i]
    long[i] = layout_data.longitude[i]
    # println(layout_data.model[i])
    if ismissing(layout_data.elevation[i])
        hub_height[i] = (layout_data.elevation[i-1] + layout_data.elevation[i+1])/2.0 + layout_data.hub_height[i]
        global m = m + 1
    else
        hub_height[i] = layout_data.elevation[i] + layout_data.hub_height[i]
    end
end
utm_zone = 10
turbine_x, turbine_y = ff.latlong_to_xy(lat, long, utm_zone)
# turbine_x = [0.0]
# turbine_y = [0.0]

# calculate the number of turbines
nturbines = length(turbine_x)

# set turbine base heights
turbine_z = zeros(nturbines)

# set turbine yaw values
turbine_yaw = zeros(nturbines)

# load ge turb data
turbge1p5file = string("./inputfiles/ge15-77_thrust_power.csv")

gedata = CSV.read(turbge1p5file)

ngepoints = length(gedata.id)
gespeed = zeros(ngepoints)
gecp= zeros(ngepoints)
gethrust = zeros(ngepoints)
gerd = gedata.rd[1]
gemodel = gedata.model[1]
for i in 1:ngepoints
    gespeed[i] = gedata.wind_speed[i]
    gecp[i] = gedata.power[i]
    gethrust[i] = gedata.thrust[i]
end

# https://www.thewindpower.net/turbine_en_56_ge-energy_1.5sl.php
gecutin = 3.0
gecutout = 25.0
geratedspeed =  15.0
geratedpower =  1.5E6

# load mm turb data
turbmm92file = string("./inputfiles/mm92_thrust_power.csv")

mmdata = CSV.read(turbmm92file)

nmmpoints = length(mmdata.id)
mmspeed = zeros(nmmpoints)
mmcp = zeros(nmmpoints)
mmthrust = zeros(nmmpoints)
mmrd = mmdata.rd[1]
mmmodel = mmdata.model[1]
for i in 1:nmmpoints
    mmspeed[i] = mmdata.wind_speed[i]
    mmcp[i] = mmdata.power[i]
    mmthrust[i] = mmdata.thrust[i]
end

# https://www.thewindpower.net/turbine_en_15_repower_mm92.php
mmcutin = 3.0
mmcutout = 25.0
mmratedspeed =  10.5
mmratedpower =  2.0E6

rotor_diameter = zeros(nturbines) # m
cut_in_speed = zeros(nturbines)   # m/s
cut_out_speed = zeros(nturbines)   # m/s
rated_speed = zeros(nturbines)  # m/s
rated_power = zeros(nturbines)  # W
generator_efficiency = zeros(nturbines) .+ 1.0

for i in 1:nturbines

    if layout_data.model[i] == gemodel
        rotor_diameter[i] = gerd    # m
        cut_in_speed[i] = gecutin   # m/s
        cut_out_speed[i] = gecutout   # m/s
        rated_speed[i] = geratedspeed  # m/s
        rated_power[i] = geratedpower  # W
    elseif layout_data.model[i] == mmmodel
        rotor_diameter[i] = mmrd    # m
        cut_in_speed[i] = mmcutin   # m/s
        cut_out_speed[i] = mmcutout   # m/s
        rated_speed[i] = mmratedspeed  # m/s
        rated_power[i] = mmratedpower  # W
    end
end

# rotor swept area sample points (normalized by rotor radius)
rotor_points_y = [0.0]
rotor_points_z = [0.0]

fig, ax = plt.subplots()
plot([minimum(turbine_x),minimum(turbine_x),maximum(turbine_x), maximum(turbine_x),
minimum(turbine_x)],[minimum(turbine_y),maximum(turbine_y),maximum(turbine_y),
minimum(turbine_y),minimum(turbine_y)])

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[i]/2.0, fill=false,color="C1", linestyle="--"))
end

# println(nturbines)
# println(m)
# println(minimum(layout_data.hub_height))
# # set up and show plot
# axis("square")
plt.show()

# set flow parameters
wind_rose_file_name = string("./inputfiles/wind_rose_jepson_prarie.txt")
winddata = readdlm(wind_rose_file_name, ',', skipstart=9)
# println(winddata)
nstates = length(winddata[:,1])*6

speeds = [3.0, 5.0, 7.0, 9.0, 11.0, 13.0]
winddirections = zeros(nstates)
windspeeds = zeros(nstates)
windprobabilities = zeros(nstates)
for i in 1:length(winddata[:,1])
    dirstr = split(winddata[i,1] , r"[-]")
    for j in 1:6
        winddirections[(i-1)*6 + j] = (0.5*(parse(Float64, dirstr[2])+parse(Float64,dirstr[1])))
        windspeeds[(i-1)*6 + j] = speeds[j]
        windprobabilities[(i-1)*6 + j] = winddata[i,j+2] / 100.0
    end
end

winddirections *= pi/180.0

air_density = 1.1716  # kg/m^3
shearexponent = 0.15
ambient_tis = zeros(nstates) .+ 0.1
measurementheight = zeros(nstates) .+ sum(hub_height)/nturbines

# initialize power model
gepower_model = ff.PowerModelCpPoints(gespeed, gecp)
mmpower_model = ff.PowerModelCpPoints(mmspeed, mmcp)
power_models = Vector{typeof(gepower_model)}(undef, nturbines)
for i = 1:nturbines
    if model[i] == gemodel
        power_models[i] = gepower_model
    elseif model[i] == mmmodel
        power_models[i] = mmpower_model
    else
        throw(UndefVarError(:model))
    end
end

# initialize thurst model
gect_model = ff.ThrustModelCtPoints(gespeed, gethrust)
mmct_model = ff.ThrustModelCtPoints(mmspeed, mmthrust)
ct_models = Vector{typeof(gect_model)}(undef, nturbines)
for i = 1:nturbines
    if model[i] == gemodel
        ct_models[i] = gect_model
    elseif model[i] == mmmodel
        ct_models[i] = mmct_model
    else
        throw(UndefVarError(:model))
    end
end

# initialize wind shear model
wind_shear_model = ff.PowerLawWindShear(shearexponent)

# get sorted indecies
sorted_turbine_index = sortperm(turbine_x)

# initialize the wind resource definition
windresource = ff.DiscretizedWindResource([winddirections[1]], [windspeeds[1]], [windprobabilities[1]], [measurementheight[1]], air_density, [ambient_tis[1]], wind_shear_model)

# set up wake and related models
wakedeficitmodel = ff.GaussYawVariableSpread()

wakedeflectionmodel = ff.GaussYawVariableSpreadDeflection()
wakecombinationmodel = ff.LinearLocalVelocitySuperposition()
localtimodel = ff.LocalTIModelMaxTI()

# initialize model set
model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

struct SubProblem{T <: AbstractFloat, F <: AbstractFeasibilityAlgorithm, I <: AbstractIntegerAlgorithm}
    id::Int
    probability::T
    model::JuMP.Model
    optimizer::MOI.AbstractOptimizer
    feasibility_algorithm::F
    integer_algorithm::I
    linking_constraints::Vector{MOI.ConstraintIndex}
    masterterms::Vector{Vector{Tuple{Int, Int, T}}}

    function SubProblem(model::JuMP.Model,
                        id::Integer,
                        π::AbstractFloat,
                        feasibility_strategy::AbstractFeasibilityStrategy,
                        integer_strategy::AbstractIntegerStrategy)
        T = typeof(π)
        # Get optimizer backend
        optimizer = backend(model)
        # Collect all constraints with known decision occurances
        constraints, terms =
            collect_linking_constraints(model,
                                        T)
        # Instantiate feasibility algorithm
        feasibility_algorithm = worker(feasibility_strategy, constraints, optimizer)
        F = typeof(feasibility_algorithm)
        # Instantiate integer algorithm
        integer_algorithm = worker(integer_strategy, T)
        I = typeof(integer_algorithm)

        subproblem =  new{T,F,I}(id,
                                 π,
                                 model,
                                 optimizer,
                                 feasibility_algorithm,
                                 integer_algorithm,
                                 constraints,
                                 terms)
        # Initialize integer algorithm
        initialize_integer_algorithm!(subproblem)
        return subproblem
    end
end

# Subproblem methods #
# ========================== #
function collect_linking_constraints(model::JuMP.Model,
                                     ::Type{T}) where T <: AbstractFloat
    linking_constraints = Vector{MOI.ConstraintIndex}()
    masterterms = Vector{Vector{Tuple{Int, Int, T}}}()
    master_indices = index.(all_known_decision_variables(model, 1))
    # Parse single rows
    F = AffineDecisionFunction{T}
    for S in [MOI.EqualTo{Float64}, MOI.LessThan{Float64}, MOI.GreaterThan{Float64}]
        for ci in MOI.get(backend(model), MOI.ListOfConstraintIndices{F,S}())
            coeffs = Vector{Tuple{Int, Int, T}}()
            f = MOI.get(backend(model), MOI.ConstraintFunction(), ci)::AffineDecisionFunction{T}
            aff = JuMP.jump_function(model, f)::DecisionAffExpr{T}
            for (coef, kvar) in linear_terms(aff.decisions)
                # Map known decisions to master decision,
                # assuming sorted order
                if state(kvar) == Known
                    col = master_indices[index(kvar).value].value
                    push!(coeffs, (1, col, T(coef)))
                end
            end
            if !isempty(coeffs)
                push!(masterterms, coeffs)
                push!(linking_constraints, ci)
            end
        end
    end
    # Parse vector rows
    F = VectorAffineDecisionFunction{T}
    for S in [MOI.Zeros, MOI.Nonpositives, MOI.Nonnegatives]
        for ci in MOI.get(backend(model), MOI.ListOfConstraintIndices{F,S}())
            coeffs = Vector{Tuple{Int, Int, T}}()
            f = MOI.get(backend(model), MOI.ConstraintFunction(), ci)::VectorAffineDecisionFunction{T}
            affs = JuMP.jump_function(model, f)::Vector{DecisionAffExpr{T}}
            for (row, aff) in enumerate(affs)
                for (coef, kvar) in linear_terms(aff.decisions)
                    # Map known decisions to master decision,
                    # assuming sorted order
                    if state(kvar) == Known
                        col = master_indices[index(kvar).value].value
                        push!(coeffs, (row, col, T(coef)))
                    end
                end
            end
            if !isempty(coeffs)
                push!(masterterms, coeffs)
                push!(linking_constraints, ci)
            end
        end
    end
    return linking_constraints, masterterms
end

function update_subproblem!(subproblem::SubProblem)
    update_known_decisions!(subproblem.optimizer)
    return nothing
end

function restore_subproblem!(subproblem::SubProblem)
    restore!(subproblem.optimizer, subproblem.integer_algorithm)
    restore!(subproblem.optimizer, subproblem.feasibility_algorithm)
    return nothing
end

function (subproblem::SubProblem)(x::AbstractVector, metadata)
    return solve_subproblem(subproblem,
                            metadata,
                            subproblem.feasibility_algorithm,
                            subproblem.integer_algorithm,
                            x)
end

function solve_subproblem(subproblem::SubProblem, x::AbstractVector)
    MOI.optimize!(subproblem.optimizer)
    status = MOI.get(subproblem.optimizer, MOI.TerminationStatus())
    if status ∈ AcceptableTermination
        return OptimalityCut(subproblem, x)
    elseif status == MOI.INFEASIBLE
        return Infeasible(subproblem, x)
    elseif status == MOI.DUAL_INFEASIBLE
        return Unbounded(subproblem)
    else
        error("Subproblem $(subproblem.id) was not solved properly, returned status code: $status")
    end
end

# Cuts #
# ========================== #
function OptimalityCut(subproblem::SubProblem{T}, x::AbstractVector) where T <: AbstractFloat
    π = subproblem.probability
    nterms = if isempty(subproblem.masterterms)
        nterms = 0
    else
        nterms = mapreduce(+, subproblem.masterterms) do terms
            length(terms)
        end
    end
    cols = zeros(Int, nterms)
    vals = zeros(T, nterms)
    j = 1
    for (i, ci) in enumerate(subproblem.linking_constraints)
        λ = MOI.get(subproblem.optimizer, MOI.ConstraintDual(), ci)
        for (row, col, coeff) in subproblem.masterterms[i]
            cols[j] = col
            vals[j] = π * λ[row] * coeff
            j += 1
        end
    end
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    # Create sense-corrected optimality cut
    δQ = sparsevec(cols, vals, length(x))
    q = correction * π * MOI.get(subproblem.optimizer, MOI.ObjectiveValue()) + δQ⋅x
    return OptimalityCut(δQ, q, subproblem.id)
end

function FeasibilityCut(subproblem::SubProblem{T}, x::AbstractVector) where T <: AbstractFloat
    nterms = if isempty(subproblem.masterterms)
        nterms = 0
    else
        nterms = mapreduce(+, subproblem.masterterms) do terms
            length(terms)
        end
    end
    cols = zeros(Int, nterms)
    vals = zeros(T, nterms)
    j = 1
    for (i, ci) in enumerate(subproblem.linking_constraints)
        λ = MOI.get(subproblem.optimizer, MOI.ConstraintDual(), ci)
        for (row, col, coeff) in subproblem.masterterms[i]
            cols[j] = col
            vals[j] = λ[row] * coeff
            j += 1
        end
    end
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    # Create sense-corrected optimality cut
    G = sparsevec(cols, vals, length(x))
    g = correction * MOI.get(subproblem.optimizer, MOI.ObjectiveValue()) + G⋅x
    return FeasibilityCut(G, g, subproblem.id)
end

function Infeasible(subproblem::SubProblem{T}, x) where T
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    return Infeasible(correction * Inf, subproblem.id)
end

function Unbounded(subproblem::SubProblem)
    # Get sense
    sense = MOI.get(subproblem.optimizer, MOI.ObjectiveSense())
    correction = (sense == MOI.MIN_SENSE || sense == MOI.FEASIBILITY_SENSE) ? 1.0 : -1.0
    return Unbounded(correction * -Inf, subproblem.id)
end


text:
using Reactive
using FactCheck

# Stop the runner task

facts("Queue runner") do
    @fact istaskdone(Reactive.runner_task[]) --> false
    Reactive.stop()
    @fact istaskdone(Reactive.runner_task[]) --> true
end

step() = Reactive.run(1)
queue_size() = Base.n_avail(Reactive._messages)
number() = round(Int, rand()*1000)

include("basics.jl")
include("push_to_non_input.jl")
# include("gc.jl")
include("node_order.jl")
include("call_count.jl")
include("flatten.jl")
include("time.jl")
include("async.jl")
include("queue_runner.jl")
FactCheck.exitstatus()


text:
using DataFrames
using MixedModels
using Test

# convert to DataFrame to modify it
slp = DataFrame(MixedModels.dataset(:sleepstudy))
allowmissing!(slp, :days)
slp[1,:days] = missing

# TODO: re-enable this test when better missing support has landed in StatsModels
# @testset "No impact from missing on schema" begin
#     f = @formula(Y ~ 1 + U + (1|G))
#     contrasts =  Dict{Symbol,Any}()
#     form = apply_schema(f, schema(f, dat[:sleepstudy], contrasts), LinearMixedModel)
#     form_missing = apply_schema(f, schema(f, slp, contrasts), LinearMixedModel)
#
#     @test form.lhs == form_missing.lhs
#     @test form.rhs == form_missing.rhs
# end

@testset "Missing Omit" begin
    @testset "Missing from unused variables" begin
        # missing from unused variables should have no impact
        m1 = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), MixedModels.dataset(:sleepstudy), progress=false)
        m1_missing = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), slp, progress=false)
        @test isapprox(m1.θ, m1_missing.θ, rtol=1.0e-12)
    end

    @testset "Missing from used variables" begin
        m1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), MixedModels.dataset(:sleepstudy), progress=false)
        m1_missing = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), slp, progress=false)
        @test nobs(m1) - nobs(m1_missing) == 1
    end
end