text:
import CSV
using DataFrames

function read_spatial_df(data_path; x_col::Symbol=:x, y_col::Symbol=:y, gene_col::Union{Symbol, Nothing}=:gene, filter_cols::Bool=false)
    df_spatial = DataFrame!(CSV.File(data_path));

    for (cn, co) in zip((:x, :y, :gene), (x_col, y_col, gene_col))
        if co == nothing
            continue
        end

        if (cn in propertynames(df_spatial)) & (cn != co)
            cr = Symbol(String(cn) * "_reserved")
            if cr in propertynames(df_spatial)
                DataFrames.select!(df_spatial, DataFrames.Not(cr))
            end
            DataFrames.rename!(df_spatial, cn => cr);
        end

        DataFrames.rename!(df_spatial, co => cn);
    end

    if gene_col === nothing
        if filter_cols
            df_spatial = df_spatial[:, [:x, :y]]
        end
    else
        if filter_cols
            df_spatial = df_spatial[:, [:x, :y, :gene]]
        end
        df_spatial[!, :gene] = ["$g" for g in df_spatial.gene]
    end

    return df_spatial
end

function convert_segmentation_to_counts(genes::Vector{Int}, cell_assignment::Vector{Int}; drop_empty_labels::Bool=false, gene_names::Union{Vector{String}, Nothing}=nothing)
    if drop_empty_labels
        if minimum(cell_assignment) == 0
            cell_assignment = denserank(cell_assignment) .- 1
        else
            cell_assignment = denserank(cell_assignment)
        end
    end

    cm = zeros(Int, maximum(genes), maximum(cell_assignment))
    for i in 1:length(genes)
        if cell_assignment[i] == 0
            continue
        end
        cm[genes[i], cell_assignment[i]] += 1
    end

    if gene_names !== nothing
        cm = DataFrame(cm, [Symbol("$c") for c in 1:size(cm, 2)])
        cm[!, :gene] = gene_names
        cm = cm[:, vcat(end, 1:end-1)]
    end

    return cm
end

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
module TemperatureProfiles

import DocStringExtensions
const DSE = DocStringExtensions

export TemperatureProfile,
    IsothermalProfile, DecayingTemperatureProfile, DryAdiabaticProfile

import CLIMAParameters
const CP = CLIMAParameters
const APS = CP.AbstractParameterSet

import ..InternalClimaParams
const ICP = InternalClimaParams

"""
    TemperatureProfile

Specifies the temperature or virtual temperature profile for a reference state.

Instances of this type are required to be callable objects with the following signature

    T,p = (::TemperatureProfile)(param_set::APS, z::FT) where {FT}

where `T` is the temperature or virtual temperature (in K), and `p` is the pressure (in Pa).
"""
abstract type TemperatureProfile{FT} end

"""
    IsothermalProfile(param_set, T_virt)
    IsothermalProfile(param_set, ::Type{FT<:AbstractFloat})

A uniform virtual temperature profile, which is implemented
as a special case of [`DecayingTemperatureProfile`](@ref).
"""
IsothermalProfile(param_set::APS, T_virt::FT) where {FT} =
    DecayingTemperatureProfile{FT}(param_set, T_virt, T_virt)

function IsothermalProfile(param_set::APS, ::Type{FT}) where {FT}
    T_virt = FT(ICP.T_surf_ref(param_set))
    return DecayingTemperatureProfile{FT}(param_set, T_virt, T_virt)
end

"""
    DryAdiabaticProfile{FT} <: TemperatureProfile{FT}

A temperature profile that has uniform dry potential temperature `θ`

# Fields

$(DSE.FIELDS)
"""
struct DryAdiabaticProfile{FT} <: TemperatureProfile{FT}
    "Surface temperature (K)"
    T_surface::FT
    "Minimum temperature (K)"
    T_min_ref::FT
    function DryAdiabaticProfile{FT}(
        param_set::APS,
        T_surface::FT = FT(ICP.T_surf_ref(param_set)),
        _T_min_ref::FT = FT(ICP.T_min_ref(param_set)),
    ) where {FT}
        return new{FT}(T_surface, _T_min_ref)
    end
end

"""
    (profile::DryAdiabaticProfile)(
        param_set::APS,
        z::FT,
    ) where {FT}

Returns dry adiabatic temperature and pressure profiles
with zero relative humidity. The temperature is truncated
to be greater than or equal to `profile.T_min_ref`.
"""
function (profile::DryAdiabaticProfile)(param_set::APS, z::FT) where {FT}

    R_d::FT = ICP.R_d(param_set)
    cp_d::FT = ICP.cp_d(param_set)
    grav::FT = ICP.grav(param_set)
    MSLP::FT = ICP.MSLP(param_set)

    # Temperature
    Γ = grav / cp_d
    T = max(profile.T_surface - Γ * z, profile.T_min_ref)

    # Pressure
    p = MSLP * (T / profile.T_surface)^(grav / (R_d * Γ))
    if T == profile.T_min_ref
        z_top = (profile.T_surface - profile.T_min_ref) / Γ
        H_min = R_d * profile.T_min_ref / grav
        p *= exp(-(z - z_top) / H_min)
    end
    return (T, p)
end

"""
    DecayingTemperatureProfile{F} <: TemperatureProfile{FT}

A virtual temperature profile that decays smoothly with height `z`, from
`T_virt_surf` to `T_min_ref` over a height scale `H_t`. The default height
scale `H_t` is the density scale height evaluated with `T_virt_surf`.

```math
T_{\\text{v}}(z) = \\max(T_{\\text{v, sfc}} − (T_{\\text{v, sfc}} - T_{\\text{v, min}}) \\tanh(z/H_{\\text{t}})
```

# Fields

$(DSE.FIELDS)
"""
struct DecayingTemperatureProfile{FT} <: TemperatureProfile{FT}
    "Virtual temperature at surface (K)"
    T_virt_surf::FT
    "Minimum virtual temperature at the top of the atmosphere (K)"
    T_min_ref::FT
    "Height scale over which virtual temperature drops (m)"
    H_t::FT
    function DecayingTemperatureProfile{FT}(
        param_set::APS,
        _T_virt_surf::FT = FT(ICP.T_surf_ref(param_set)),
        _T_min_ref::FT = FT(ICP.T_min_ref(param_set)),
        H_t::FT = FT(ICP.R_d(param_set)) * _T_virt_surf /
                  FT(ICP.grav(param_set)),
    ) where {FT}
        return new{FT}(_T_virt_surf, _T_min_ref, H_t)
    end
end

function (profile::DecayingTemperatureProfile)(param_set::APS, z::FT) where {FT}
    R_d::FT = ICP.R_d(param_set)
    grav::FT = ICP.grav(param_set)
    MSLP::FT = ICP.MSLP(param_set)

    # Scale height for surface temperature
    H_sfc = R_d * profile.T_virt_surf / grav
    H_t = profile.H_t
    z′ = z / H_t
    tanh_z′ = tanh(z′)

    ΔTv = profile.T_virt_surf - profile.T_min_ref
    Tv = profile.T_virt_surf - ΔTv * tanh_z′

    ΔTv′ = ΔTv / profile.T_virt_surf
    p = -H_t * (z′ + ΔTv′ * (log(1 - ΔTv′ * tanh_z′) - log(1 + tanh_z′) + z′))
    p /= H_sfc * (1 - ΔTv′^2)
    p = MSLP * exp(p)
    return (Tv, p)
end

end

ext: jl


text:
abstract type AbstractConstruct{T} end
######################################
# Base functionalities for AbstractConstruct
######################################

"""
    Base.eltype(::AbstractConstruct{T}) where {T} = T

Returns the type that is generated by the iterator on an `AbstractConstruct{T}` structure
"""
Base.eltype(::AbstractConstruct{T}) where {T} = T

"""
    Base.getindex(Construct::AbstractConstruct,i::Int)

Returns the module on the i-th possition in the construct
Output of type `Mod`
"""
Base.getindex(c::AbstractConstruct,i::Int) = c.c[i]

"""
    Base.getindex(Construct::AbstractConstruct,i::UnitRange)

Returns a sub-part of the original construct.
Using the indexes of the `UnitRange` object. Extention to allow UnitRange indexing on a `AbstractConstruct`.
The type of output type is simular as the larger input construct
"""

Base.getindex(c::AbstractConstruct,i::UnitRange) = c.c[i] |> typeof(c)

"""
    Base.length(Construct::AbstractConstruct)

Returns the lenght of the input construct, the number of modules that are used
"""
Base.length(c::AbstractConstruct) = length(c.c)

"""
    Base.lastindex(Construct::AbstractConstruct)

Returns the lenght of the input construct, the number of modules that are used
"""
Base.lastindex(c::AbstractConstruct) =  last(eachindex(IndexLinear(), c.c))

"""
    Base.size(Construct::AbstractConstruct)

Returns the size a construct, a tuple of two dimensions. The first is the length of the construct the second equals to one (1D constructs)
Allows to use of `end` when indexing on a construct to obtain the last module.

"""
Base.size(c::AbstractConstruct) = (length(c.c),1)

"""
    Base.isequal(c1::AbstractConstruct,c2::AbstractConstruct)

Evaluates if two constructs are equal. If the the type of to constructs differces
If the type of two constructs aren't equal, it returns `false`, even if the same modules are present.
"""
Base.isequal(c1::AbstractConstruct,c2::AbstractConstruct) = typeof(c1) == typeof(c2) ? isequal(c1,c2) : false

"""
    Base.iterate(Construct::AbstractConstruct,state=1)

Iterator to loop over all modules in a give construct.
"""
function Base.iterate(Construct::AbstractConstruct,state=1)
    if  state <= length(Construct)
        mod = Construct[state]
        state += 1
        return (mod,state)
    else
        return
    end
end

##################
#New structs
##################

"""
    Ordered_Construct{T}  <: AbstractConstruct{T}

Structure that contains an ordered construct. In ordered construct, the order of the different modules is important.
For example: `Ordered_Construct([a,b,c])` is different compared to  `Ordered_Construct([c,b,a])`.
"""

struct Ordered_Construct{T<:Mod}  <: AbstractConstruct{T}
    c::Array{T}
end

"""
    Base.isequal(c1::Ordered_Construct,c2::Ordered_Construct)

Evaluates if two Ordered_Construct constructs are equal.
"""

Base.isequal(c1::Ordered_Construct,c2::Ordered_Construct) = isequal(c1.c,c2.c)

"""
    Base.*(m1::Mod,m2::Mod)

Returns a `Ordered_Construct` containing both modules `m1` and `m2`.Order of input argument is maintained in the returned construct.
To allow efficient on the fly calculation using Kronecker product the `Base.*`should be defined between two `Mod`
"""
*(m1::Mod,m2::Mod) = Ordered_Construct([m1,m2])

"""
    Base.*(c::Ordered_Construct ,m::Mod)

Returns a new `Ordered_Construct` where the module `m` is concatenated to the end of the input construct `c`.

"""
*(c::Ordered_Construct ,m::Mod) = push!(c.c,m) |> Ordered_Construct

"""
    *(m::Mod, c::Ordered_Construct)

Returns a new `Ordered_Construct` where the module `m` is concatenated in front of the input construct `c`.
Probably less efficient than `*(c::Ordered_Construct ,m::Mod)` because of a `firstpush!`
"""
*(m::Mod,  c::Ordered_Construct) = pushfirst!(c.c,m) |> Ordered_Construct

"""
    *(c1::Ordered_Construct,c2::Ordered_Construct)

Concatenated two ordered constructs `c1` and `c2`,
where the construct that is in the second argument position will be appended to the end of the construct in the first argument position

```jldoctest
m1 = Mod("a")
m2 = Mod("b")
m3 = Mod("c")
c1 = Ordered_Construct([m1,m2])
c2 = Ordered_Construct([m1,m3])
c1*c2

# output
Ordered_Construct([m1,m2,m1,m4])
```
"""

*(c1::Ordered_Construct ,c2::Ordered_Construct) = vcat(c1.c,c2.c) |> Ordered_Construct

"""
    Unordered_Construct{T} <: AbstractConstruct{T}

Structure that contains an unordered construct. In unordered construct, the order of the different modules has no meaning .
For example: `Unordered_Construct([a,b,c])` is considered equal to `Unordered_Construct([c,b,a])`.
"""
# currently the underlying structure is an array, a Set is can be a more logical alternative
# still for the first implementation and test, an array structure was preferred. An array can always be transformed Set where needed.
# The users shouldn't worry too much about this more to run everything smoothly internally

struct Unordered_Construct{T<:Mod} <: AbstractConstruct{T}
    c::Array{T}
end

"""
    +(m1::Mod,m2::Mod)

Returns a `Unordered_Construct` containing both modules `m1` and `m2`. Order of input argument is maintained in the returned construct.
The order of the modules doesn't have useful meaning, but it has some advantages internally to keep the order.
"""

+(m1::Mod,m2::Mod) = Unordered_Construct([m1,m2])
"""
    Base.*(c::Unordered_Construct ,m::Mod)

Returns a new `Unordered_Construct` where the module `m` is concatenated to the end of the input construct `c`.

"""
+(c::Unordered_Construct ,m::Mod) = push!(c.c,m) |> Unordered_Construct

"""
    Base.*(c::Ordered_Construct ,m::Mod)

Returns a new `Unordered_Construct` where the module `m` is concatenated in front of the input construct `c`

"""
+(m::Mod, c::Unordered_Construct) = pushfirst!(c.c,m) |> Unordered_Construct

"""
    *(c1::Unordered_Construct,c2::Unordered_Construct)

Concatenated two `Unordered_Construct` `c1` and `c2`,
where the construct that is in the second argument position will be appended to the end of the construct in the first argument position
"""

+(c1::Unordered_Construct,c2::Unordered_Construct) = vcat(c1.c,c2.c) |> Unordered_Construct

"""
    Base.isequal(c1::Unordered_Construct,c2::Unordered_Construct)

Evaluates if two `Unordered_Construct` are equal.
"""

Base.isequal(c1::Unordered_Construct,c2::Unordered_Construct) = isequal(Set(c1.c),Set(c2.c))


text:
function _simplify_angle{R}(angle::R)
    two_pi = R(2) * π
    ϕ = angle % two_pi
    ϕ < zero(R) ? two_pi + ϕ : ϕ
end

function _expand_size{R<:Real}(height, width, ϕ::R)
    two_pi = R(2) * π
    half_pi = π / R(2)

    θ = if ϕ > R(3) * half_pi
        two_pi - ϕ
    elseif ϕ > π
        ϕ - π
    elseif ϕ > half_pi
        π - ϕ
    else
        convert(typeof(two_pi), ϕ)
    end

    h1 = sin(θ) * width
    h2 = sin(half_pi - θ) * height
    height_new = floor(Int, h1 + h2)

    w1 = cos(θ) * width
    w2 = cos(half_pi - θ) * height
    width_new = floor(Int, w1 + w2)
    height_new, width_new
end

function _ta_rotate{T}(A::AbstractMatrix{T}, ϕ)
    itp = interpolate(A, BSpline(Linear()), OnGrid())
    etp = extrapolate(itp, zero(T))

    tfm = tformrotate(ϕ)
    TransformedArray(etp, tfm)
end

function rotate_expand{T}(A::AbstractMatrix{T}, angle)
    ϕ = _simplify_angle(angle)

    if ϕ == 0
        A
    elseif ϕ == π/2
        rotr90(A)
    elseif ϕ == 3π/2
        rotl90(A)
    elseif ϕ == 1π
        Images.reflect(A)
    else
        A_tfm = _ta_rotate(A, ϕ)
        h, w = _expand_size(size(A,1), size(A,2), ϕ)
        B = similar(A, h, w)
        transform!(B, A_tfm)
    end
end

function rotate_expand(img::AbstractImage, angle)
    assert2d(img)
    shareproperties(img, rotate_expand(data(img), isxfirst(img) ? angle : -angle))
end

function rotate_crop{T}(A::AbstractMatrix{T}, angle)
    ϕ = _simplify_angle(angle)

    if ϕ == 0
        A
    elseif ϕ == 1π
        Images.reflect(A)
    else
        A_tfm = _ta_rotate(A, ϕ)
        B = similar(A)
        transform!(B, A_tfm)
    end
end

function rotate_crop(img::AbstractImage, angle)
    assert2d(img)
    shareproperties(img, rotate_crop(data(img), isxfirst(img) ? angle : -angle))
end


text:
import RData
using RNGPool
using CoupledConditionalSMC

include("demo.jl")
include("lgModel.jl")

# data generated with A = 0.95 (not 0.9)
dataJLS = RData.load("demo/ar1data.RData")
ys = vec(dataJLS["observations"])

# they start at time 0 with observations from time 1,2,...
# so we have an initial distribution of N(0,1+0.9^2)
theta = LGTheta(0.9, 1.0, 1.0, 1.0, 0.0, 1.81)
lM = LinearGaussian.makelM(theta)

setRNGs(12345)

# one can couple in a few hundred iterations using 4096 particles with n=500
# but with n=10000 one would require millions or billions of particles
# and this would exhaust main memory
model = makeLGModel(theta, ys[1:8000])
visualizeCCSMC(model, 2^12, true, true)
CoupledConditionalSMC.couplingTimes(model, 2^12, 4000, 10, true, true, 2000)

model = makeLGModel(theta, ys[1:500])
visualizeCCSMC(model, lM, 128, :BS, true, true)
visualizeCCSMC(model, lM, 128, :AS, true, true)
CoupledConditionalSMC.couplingTimes(model, lM, 128, 100, 10, :BS, true, true, 2000)
tmp = CoupledConditionalSMC.couplingTimes(model, lM, 128, 100, 100, :AS, true, true, 2000)
maximum(tmp)

model = makeLGModel(theta, ys[1:100])
CoupledConditionalSMC.couplingTimes(model, lM, 128, 100, 10, :BS, true, true, 2000)

model = makeLGModel(theta, ys[1:500])
visualizeCCSMC(model, lM, 128, :BS, true, true)
visualizeCCSMC(model, lM, 128, :AS, true, true)

println("\nUsing a linear Gaussian Model with n = 1000, obsVariance = 100.0")
model, lM, _ = setupLGModel(1000, 100.0)
visualizeCCSMC(model, lM, 512, :BS, true, false)
visualizeCCSMC(model, 512, true, false, 100)
visualizeCCSMC(model, lM, 512, :AS, true, false, 100)

visualizeCCSMC(model, lM, 512, :BS, true, true)
visualizeCCSMC(model, 512, true, true, 100)
visualizeCCSMC(model, lM, 512, :AS, true, true, 100)


text:
using JuMP, GLPK, Ipopt
using DualDecomposition
using DRMSMIP
using Random

const DD = DualDecomposition
const parallel = DD.parallel

const rng = Random.MersenneTwister(1234)

"""
a: interest rate
π: unit stock price
ρ: unit dividend price

K: number of stages
L: number of stock types
2^L scenarios in each stage
2^L^(K-1)=16 scenarios in total
ρ = 0.05 * π
bank: interest rate 0.01
stock1: 1.03 or 0.97
stock2: 1.06 or 0.94
...

b_k: initial asset (if k=1) and income (else)
B_k: money in bank
x_{k,l}: number of stocks to buy/sell (integer)
y_{k,l}: total stocks

deterministic model:

    max     B_K+∑_{l=1}^{L}π_{K,l}y_{K,l}

    s.t.    B_1+∑_{l=1}^{L}π_{1,l}x_{1,l} = b_1

            b_k+(1+a)B_{k-1}+∑_{l=1}^{L}ρ_{k,l}y_{k-1,l} = B_k+∑_{l=1}^{L}π_{k,l}x_{k,l}, ∀ k=2,…,K

            y_{1,l} = x_{1,l}, ∀ l=1,…,L

            y_{k-1,l}+x_{k,l} = y_{k,l}, ∀ k=2,…,K, l=1,…,L

            x_{k,l} ∈ ℤ , ∀ k=1,…,K, l=1,…,L

            y_{k,l} ≥ 0, ∀ k=1,…,K, l=1,…,L

            B_k ≥ 0, ∀ k=1,…,K.
"""
const K = 3
const L = 2
const a = 0.01
const b_init = 100  # initial capital
const b_in = 30   # income

"""
In each node, we have Np=10 samples from a log-normal distribution
"""
const Np = 10 # number of samples

function generate_sample(π::Array{Float64})::Array{DRMSMIP.Sample}
    # generates random samples following a lognormal distribution
    ret = Array{DRMSMIP.Sample}(undef, Np)
    for ii in 1:Np
        π_samp = Array{Float64}(undef, L)
        for l in 1:L
            sig = sqrt( log( 0.5+sqrt( ( 0.03*l )^2+0.25 ) ) )
            rnd = sig * randn(rng) .+ log(π[l])
            π_samp[l] = exp(rnd)
        end
        ret[ii] = DRMSMIP.Sample(Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}(:π => π_samp), 1/Np)
    end
    return ret
end

# iteratively add nodes
# root node
function create_nodes()::DD.Tree
    ξ = Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}(:π => ones(L))
    ξ_samp = generate_sample(ξ[:π])
    set = DRMSMIP.WassersteinSet(ξ_samp, 1.0, DRMSMIP.norm_L1)
    tree = DD.Tree(ξ, set)

    #subproblem formulation
    function subproblem_builder(mdl::JuMP.Model, node::DD.SubTreeNode)
        @variable(mdl, x[l=1:L], DD.ControlInfo, subnode = node, ref_symbol = :x)

        @variable(mdl, y[l=1:L] >= 0, DD.OutStateInfo, subnode = node, ref_symbol = :y)

        @variable(mdl, B >= 0, DD.OutStateInfo, subnode = node, ref_symbol = :B)

        π = DD.get_scenario(node)[:π]
        @constraints(mdl,
            begin
                B + sum( π[l] * x[l] for l in 1:L) == b_init
                [l=1:L], y[l] - x[l] == 0
            end
        )
        DD.set_stage_objective(node, 0.0)
    end

    DD.set_stage_builder!(tree, 1, subproblem_builder)

    create_nodes!(tree, 1)
    return tree
end

# child nodes
function create_nodes!(tree::DD.Tree, pt::Int)
    for scenario = 1:2^L
        prob = 1/2^L
        π = get_realization(DD.get_scenario(tree, pt)[:π], scenario)
        ξ = Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}(:π => π)
        if DD.get_stage(tree, pt) != K-1
            ξ_samp = generate_sample(ξ[:π])
            set = DRMSMIP.WassersteinSet(ξ_samp, 1.0, DRMSMIP.norm_L1)
            id = DD.add_child!(tree, pt, ξ, set)
        else
            id = DD.add_child!(tree, pt, ξ, nothing)
        end

        #subproblem formulation
        function subproblem_builder(mdl::JuMP.Model, node::DD.SubTreeNode)
            @variable(mdl, x[l=1:L], DD.ControlInfo, subnode = node, ref_symbol = :x)

            @variable(mdl, y[l=1:L] >= 0, DD.OutStateInfo, subnode = node, ref_symbol = :y)

            @variable(mdl, B >= 0, DD.OutStateInfo, subnode = node, ref_symbol = :B)

            @variable(mdl, y_[l=1:L] >= 0, DD.InStateInfo, subnode = node, ref_symbol = :y)

            @variable(mdl, B_ >= 0, DD.InStateInfo, subnode = node, ref_symbol = :B)

            π = DD.get_scenario(node)[:π]
            pt = DD.get_parent(node)
            ρ = DD.get_scenario(tree, pt)[:π] * 0.05
            @constraint(mdl, B + sum( π[l] * x[l] - ρ[l] * y_[l] for l in 1:L) - (1+a) * B_ == b_in)
            @constraint(mdl, [l=1:L], y[l] - x[l] - y_[l] == 0)

            #dummy bound for input variables to avoid subproblem unboundedness
            @constraint(mdl, [l=1:L], y_[l] <= 500)
            @constraint(mdl, B_ <= 500)
            if DD.get_stage(node) < K
                DD.set_stage_objective(node, 0.0)
            else
                DD.set_stage_objective(node, -(B + sum( π[l] * y[l] for l in 1:L )))
            end
        end

        DD.set_stage_builder!(tree, id, subproblem_builder)
        if DD.get_stage(tree, id) < K
            create_nodes!(tree, id)
        end
    end
end

# construct realization event
function get_realization(ξ::Array{Float64,1}, scenario::Int)::Array{Float64,1}
    ret = ones(L)
    multipliers = digits(scenario - 1, base=2, pad=L)*2 - ones(L)
    for l = 1:L
        ret[l] = ξ[l] * (1 + multipliers[l] * l * 0.03)
    end
    return ret
end

tree = create_nodes()

det = DRMSMIP.create_Wasserstein_deterministic!(tree)
set_optimizer(det.model, GLPK.Optimizer)
JuMP.optimize!(det.model)
println(JuMP.objective_value(det.model))

#node_cluster = DD.decomposition_not(tree)
node_cluster = DD.decomposition_scenario(tree)
#node_cluster = DD.decomposition_temporal(tree) #There is a DUAL_INFEASIBLE issue

# Number of block components
NS = length(node_cluster)

# Initialize MPI
parallel.init()

# Create DualDecomposition instance.
algo = DRMSMIP.DR_LagrangeDual(tree)

# partition scenarios into processes
parallel.partition(NS)

coupling_variables = Vector{DD.CouplingVariableRef}()
models = Dict{Int,JuMP.Model}()

for block_id in parallel.getpartition()
    nodes = node_cluster[block_id]
    subtree = DD.create_subtree!(block_id, coupling_variables, nodes)
    set_optimizer(subtree.model, GLPK.Optimizer)
    DD.add_block_model!(algo, block_id, subtree.model)
    models[block_id] = subtree.model
end

# Set nonanticipativity variables as an array of symbols.
DD.set_coupling_variables!(algo, coupling_variables)

bundle_init = DRMSMIP.initialize_bundle(tree, algo)
#println(bundle_init)
# for (id, node) in tree.nodes
#    println(node.cost)
# end

# Lagrange master method
#LM = DD.BundleMaster(BM.ProximalMethod, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
LM = DD.BundleMaster(BM.TrustRegionMethod, GLPK.Optimizer)

# Solve the problem with the solver; this solver is for the underlying bundle method.
DD.run!(algo, LM, bundle_init)

# Write timing outputs to files
DD.write_all(algo)

# Finalize MPI
parallel.finalize()


text:
using Dates

function hfun_datetext(date=locvar("rss_pubdate"))
    datetext = "$(monthname(date)) $(day(date)), $(year(date))"
    return datetext
end

# Based on https://github.com/JuliaLang/www.julialang.org/blob/master/utils.jl
function hfun_blogposts()
    curyear = year(Dates.today())
    io = IOBuffer()
    for year in curyear:-1:2020
        ys = "$year"
        for month in 12:-1:1
            ms = "0"^(month < 10) * "$month"
            base = joinpath("blog", ys, ms)
            isdir(base) || continue
            posts = filter!(p -> endswith(p, ".md"), readdir(base))
            days = zeros(Int, length(posts))
            lines = Vector{String}(undef, length(posts))
            for (i, post) in enumerate(posts)
                ps = splitext(post)[1]
                url = "/blog/$ys/$ms/$ps/"
                surl = strip(url, '/')
                title = pagevar(surl, :title)
                rss = pagevar(surl, :rss)
                date = pagevar(surl, :rss_pubdate)
                datetext = hfun_datetext(date)
                days[i] = day(date)
                lines[i] =
                    """
                    ~~~

                            $title

$rss

$datetext

                    ~~~
                    """
            end
            # Sort by day
            foreach(line -> write(io, line), lines[sortperm(days, rev=true)])
        end
    end
    # markdown conversion adds `<p>` beginning and end but
    # we want to avoid this to avoid an empty separator
    r = Franklin.fd2html(String(take!(io)), internal=true)
    return r
end


text:
# ====================================================================
# Component structure
#
mutable struct Lorentzian_1D <: AbstractComponent
    norm::Parameter
    center::Parameter
    fwhm::Parameter

    function Lorentzian_1D(norm::Number, center::Number, fwhm::Number)
        @assert norm  > 0
        @assert fwhm > 0

        out = new(Parameter(norm), Parameter(center), Parameter(fwhm))
        out.norm.low = 0
        out.fwhm.low = 0
        return out
    end
end

mutable struct Lorentzian_2D <: AbstractComponent
    norm::Parameter
    centerX::Parameter
    centerY::Parameter
    fwhmX::Parameter
    fwhmY::Parameter

    function Lorentzian_2D(norm::Number, centerX::Number, centerY::Number, fwhmX::Number, fwhmY::Number)
        @assert norm   > 0
        @assert fwhmX > 0
        @assert fwhmY > 0

        out = new(Parameter(norm), Parameter(centerX), Parameter(centerY), Parameter(fwhmX), Parameter(fwhmY))
        out.norm.low = 0
        out.fwhmX.low = 0
        out.fwhmY.low = 0

        return out
    end
end

Lorentzian(norm, center, fwhm) = Lorentzian_1D(norm, center, fwhm)
Lorentzian(norm, centerX, centerY, fwhmX, fwhmY) = Lorentzian_2D(norm, centerX, centerY, fwhmX, fwhmY)

# ====================================================================
# Allocate the output array and the component `cdata`
struct Lorentzian_cdata <: AbstractComponentData; end

cdata(comp::Lorentzian_1D, domain::Domain_1D) = Lorentzian_cdata()
cdata(comp::Lorentzian_2D, domain::Domain_2D) = Lorentzian_cdata()

# ====================================================================
# Evaluate component
function evaluate!(cdata::Lorentzian_cdata, output::Vector{Float64}, domain::Domain_1D,
                   norm, center, fwhm)
    x = domain[1]

    @. (output = norm /
        (1. +
         ((x - center) / fwhm)^2.
        ))
    return output
end

function evaluate!(cdata::Lorentzian_cdata, output::Vector{Float64}, domain::Domain_2D,
                   norm, centerX, centerY, fwhmX, fwhmY)
    x = domain[1]
    y = domain[2]

    @. (output = norm /
        (1. +
         ((x - centerX) / fwhmX)^2. +
         ((y - centerY) / fwhmY)^2.
         ))

    return output
end


text:
"""
    wiggle(wav; taxis=1:size(wav,1), zaxis=1:size(wav,2), sc=1, EdgeColor=:black, FaceColor=:black, Orient=:across, Overlap=true, ZDir=:normal)

Plot a set of shaded wiggles
# Arguments
* 'wav::Array': matrix of waveform columns.
* 'taxis::Array=1:size(wav,1)': time axis vector
* 'zaxis::Array=1:size(wav,2)': space axis vector
* 'sc::Float=1': scale factor/magnification.
* 'EdgeColor::Symbol=:black': Sets edge of wiggles color.
* 'FaceColor::Symbol=:black': Sets shading color of wiggles.
* 'Overlap::bool=true': How signals are scaled.
        true  - Signals overlap (default);
        false - Signals are scaled so they do not overlap.
* 'Orient::Symbol=:across': Controls orientation of wiggles.
        :across - from left to right
        :down   - from top to down
* 'ZDir::Symbol=:normal': Direction of space axis.
        :normal  - First signal at bottom (default)
        :reverse - First signal at top.

Translated by Nicholas Hausch -- MATLAB file provided by Naoki Saito
The previous MATLAB version contributors are:
 Anthony K. Booer (SLB) and Bradley Marchand (NSWC-PC)
Revised by Naoki Saito, Feb. 05, 2018
"""
function wiggle(wav; taxis=1:size(wav,1), zaxis=1:size(wav,2), sc=1, EdgeColor=:black, FaceColor=:black, Overlap=true, Orient=:across, ZDir=:normal)

# Set axes
(n,m) = size(wav)

# Sanity check
if length(taxis) != n
    error("Inconsistent taxis dimension!")
end
if length(zaxis) != m
    error("Inconsistent zaxis dimension!")
end

# For calculation purposes
maxrow = zeros(m); minrow = zeros(m)
for k = 1:m
  maxrow[k] = maximum(wav[:,k]); minrow[k] = minimum(wav[:,k])
end

# Scale the data for plotting
wamp = deepcopy(wav)
dt = mean(diff(taxis))
dz = mean(diff(zaxis))
    if Overlap
  wamp *= 2 * dz * (sc/maximum(maxrow-minrow))
else
  wmax = maximum(maxrow); wmin = minimum(minrow);
  if wmax<=0
    wmax = 0
  end
  if wmin>=0
    wmin = 0
  end
    wamp = sc*wav/(wmax-wmin)
end

# Set initial plot
t0 = minimum(taxis)
t1 = maximum(taxis)
z0 = minimum(zaxis)
z1 = maximum(zaxis)
if Orient == :down
 plot(xlims=(z0-dz,z1+dz), ylims=(t0,t1), yflip=true, legend=:none)
else
 plot(xlims=(t0,t1), ylims=(z0-dz,z1+dz), legend=:none)
end
if ZDir == :reverse
    wamp = flipdim(wamp,2)
end

# Plot each wavelet
for k = 1:m
  sig = wamp[:,k]
  t = deepcopy(taxis)
  w_sign = sign.(sig)
  for j=1:n-1
    if (w_sign[j]!=w_sign[j+1] && w_sign[j]!=0 && w_sign[j+1]!=0)
      sig = [sig; 0]
      t = [t; t[j]-sig[j]*(t[j+1]-t[j])/(sig[j+1]-sig[j])]
    end
  end
  IX = sortperm(t)
  t = t[IX]
  sig = sig[IX]
  len = length(t)
  len1 = collect(len:-1:1)
  indperm = [1:len;len1]
  inputx = t[indperm]
  inputy = zaxis[k] .+ [sig;min.(sig[len1],0)]
    # In the plot! functions below, theoretically speaking, either
    # fillrange = zaxis[k] or fillrange=[zaxis[k], zaxis[k]+dz] should be used.
    # However, those do not generate the desired plots as of O3/19/2018.
    # Somehow, the relative value of 0, i.e., fillrange=0, works well,
    # which is used temporarily.
  if Orient == :down
    plot!(inputy, inputx, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor, orientation=:v)
  else
    plot!(inputx, inputy, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor)
  end
end
plot!() # flushing the display.
end

"""
    wiggle!(wav; taxis=1:size(wav,1), zaxis=1:size(wav,2), sc=1, EdgeColor=:black, FaceColor=:black, Orient=:across, Overlap=true, ZDir=:normal)

Plot a set of shaded wiggles on the current displayed graphics
# Arguments
* 'wav::Array': matrix of waveform columns.
* 'taxis::Array=1:size(wav,1)': time axis vector
* 'zaxis::Array=1:size(wav,2)': space axis vector
* 'sc::Float=1': scale factor/magnification.
* 'EdgeColor::Symbol=:black': Sets edge of wiggles color.
* 'FaceColor::Symbol=:black': Sets shading color of wiggles.
* 'Overlap::bool=true': How signals are scaled.
        true  - Signals overlap (default);
        false - Signals are scaled so they do not overlap.
* 'Orient::Symbol=:across': Controls orientation of wiggles.
        :across - from left to right
        :down   - from top to down
* 'ZDir::Symbol=:normal': Direction of space axis.
        :normal  - First signal at bottom (default)
        :reverse - First signal at top.

Translated by Nicholas Hausch -- MATLAB file provided by Naoki Saito
The previous MATLAB version contributors are:
 Anthony K. Booer (SLB) and Bradley Marchand (NSWC-PC)
Revised by Naoki Saito, Feb. 05, 2018
"""
function wiggle!(wav; taxis=1:size(wav,1), zaxis=1:size(wav,2), sc=1, EdgeColor=:black, FaceColor=:black, Overlap=true, Orient=:across, ZDir=:normal)

# Set axes
(n,m) = size(wav)

# Sanity check
if length(taxis) != n
    error("Inconsistent taxis dimension!")
end
if length(zaxis) != m
    error("Inconsistent zaxis dimension!")
end

# For calculation purposes
maxrow = zeros(m); minrow = zeros(m)
for k = 1:m
  maxrow[k] = maximum(wav[:,k]); minrow[k] = minimum(wav[:,k])
end

# Scale the data for plotting
wamp = deepcopy(wav)
dt = mean(diff(taxis))
dz = mean(diff(zaxis))
    if Overlap
  wamp *= 2 * dz * (sc/maximum(maxrow-minrow))
else
  wmax = maximum(maxrow); wmin = minimum(minrow);
  if wmax<=0
    wmax = 0
  end
  if wmin>=0
    wmin = 0
  end
    wamp = sc*wav/(wmax-wmin)
end

# Set initial plot
t0 = minimum(taxis)
t1 = maximum(taxis)
z0 = minimum(zaxis)
z1 = maximum(zaxis)
if Orient == :down
 plot!(xlims=(z0-dz,z1+dz), ylims=(t0,t1), yflip=true, legend=:none)
else
 plot!(xlims=(t0,t1), ylims=(z0-dz,z1+dz), legend=:none)
end
if ZDir == :reverse
    wamp = flipdim(wamp,2)
end

# Plot each wavelet
for k = 1:m
  sig = wamp[:,k]
  t = deepcopy(taxis)
  w_sign = sign.(sig)
  for j=1:n-1
    if (w_sign[j]!=w_sign[j+1] && w_sign[j]!=0 && w_sign[j+1]!=0)
      sig = [sig; 0]
      t = [t; t[j]-sig[j]*(t[j+1]-t[j])/(sig[j+1]-sig[j])]
    end
  end
  IX = sortperm(t)
  t = t[IX]
  sig = sig[IX]
  len = length(t)
  len1 = collect(len:-1:1)
  indperm = [1:len;len1]
  inputx = t[indperm]
  inputy = zaxis[k] .+ [sig;min.(sig[len1],0)]
    # In the plot! functions below, theoretically speaking, either
    # fillrange = zaxis[k] or fillrange=[zaxis[k], zaxis[k]+dz] should be used.
    # However, those do not generate the desired plots as of O3/19/2018.
    # Somehow, the relative value of 0, i.e., fillrange=0, works well,
    # which is used temporarily.
  if Orient == :down
    plot!(inputy, inputx, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor, orientation=:v)
  else
    plot!(inputx, inputy, fillrange=0, fillalpha=0.75, fillcolor=FaceColor, linecolor=EdgeColor)
  end
end
plot!() # flushing the display.
end


text:
## Sliceplorer definition stuff
using Sobol

function func1d(f, pt, d)
  function(x)
    xx = copy(pt)
    xx[d] = x
    f(xx)
  end
end

function func_samples(f, rng, pt, d) :: SliceplorerDimSample
  f1d = func1d(f, pt, d)
  samples = range(rng[1], stop=rng[2], length=51) # 50 samples for now
  SliceplorerDimSample([SliceplorerPoint((x, f1d(x))) for x = samples])
end

function sample_dim(f, rng, fps, d) :: SliceplorerDim
  SliceplorerDim([func_samples(f, rng, fp, d) for fp = fps])
end

@doc raw"""
    slice1d(f, spec)

Compute the slices for a given function `f` given the ProblemSpec
range `spec`.

See also: [`plot`](@ref)

# Examples

This will plot the function $f(x) = \sum x^2$.
```julia-repl
julia> f = x -> sum(x .* x)
julia> spec = ProblemSpec("x1" => (-1.,1.), "x2" => (-1.,1.), "x3" => (-1.,1.))
julia> d = slice1d(f, spec)
julia> plot(d)
```
"""
# TODO: add type signature
function slice1d(f, spec; n=50) :: SliceSet1D
  focuspoints = gen_fps(spec, n)
  samps = sample_fps(f, focuspoints, spec, n)
  SliceSet1D(spec, focuspoints, samps)
end

function gen_fps(spec, n)
  # Create focus points for slicing
  dmins = [x[1] for (_,x) in spec]
  dmaxs = [x[2] for (_,x) in spec]
  seq = SobolSeq(dmins, dmaxs)
  hcat([next!(seq) for i = 1:n])
end

function sample_fps(f, focuspoints, spec, n) :: SliceplorerSamples
  samps = OrderedDict()
  for (i,d) in enumerate(keys(spec))
    samps[d] = sample_dim(f, spec[d], focuspoints, i)
  end
  SliceplorerSamples(samps)
end

@doc raw"""
    sliceplorer(f, spec)

Slice and plot the function `f` given the ProblemSpec range `spec`

# Examples

This will plot the function $f(x) = \sum x^2$.
```julia-repl
julia> f = x -> sum(x .* x)
julia> spec = ProblemSpec("x1" => (-1.,1.), "x2" => (-1.,1.), "x3" => (-1.,1.))
julia> sliceplorer(f, spec)
"""
function sliceplorer(f, spec; n=50)
  slices = slice1d(f, spec; n=n);
  plot(slices)
end


text:
using Test
using SWIP
using Distributions

function testintervention()
  # Normally distributed random variable with id 1
  x = 1 ~ ω -> Normal(ω, 0, 1)

  # Normally distributed random variable with id 2 and x as mean
  y = 2 ~ (ω -> Normal(ω, x(ω), 1)) <| (x,)
  x_ = 0.1
  y_ = 0.3

  # An omega object -- like a trace
  ω = Ω(Dict((1,) => x_, (2,) => y_))

  # model -- tuple-valued random variable ω -> (x(ω), y(ω)), becase we want joint pdf
  m = rt(x, y)

  # Log density of model wrt ω
  l = logpdf(m, ω)

  # Check it is what it should be
  @test l == logpdf(Normal(0, 1), x_) + logpdf(Normal(x_, 1), y_)

  # Intervened model
  v_ = 100.0

  # y had x been v_
  yⁱ = y | had(x => v_)

  # new model with Intervened variables
  mⁱ = rt(x, yⁱ)

  # log pdf of interved model on same ω
  lⁱ = logpdf(mⁱ, ω)
  logpdf(x, ω)

  @test lⁱ == logpdf(Normal(0, 1), x_) + logpdf(Normal(v_, 1), y_)
  @test lⁱ < l
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule TDLib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("TDLib")
JLLWrappers.@generate_main_file("TDLib", UUID("6cfef11d-eb86-5696-bb50-1df673c3ffa1"))
end  # module TDLib_jll


text:
# Assert the conservation of species throughout the experiment.
function assertConservation(y)
    # These are all 0-indexed, so add 1
    consList = [
        [1, 4, 5, 7, 8, 11, 12, 14, 15],  # IL2Rb
        [0, 3, 5, 6, 8],  # IL2Ra
        [9, 10, 12, 13, 15],  # IL15Ra
        [16, 17, 18],  # IL7Ra
        [2, 6, 7, 8, 13, 14, 15, 18],
    ] #gc

    # Check for conservation of species sum
    for ii in range(1, stop = length(consList))
        diff = sum(y[consList[ii] .+ 1])
        @test isapprox(diff, 0.0, atol = 1.0e-12)
    end
end

@testset "Model tests." begin
    @testset "Reaction model mass conservation." begin
        dy = zeros(gcSolver.halfL)

        gcSolver.dYdT(dy, copy(dy), rxntfR, ones(3))

        # Check for conservation of each surface receptor
        assertConservation(dy)
    end

    @testset "Full model mass conservation." begin
        rr = copy(rxntfR)
        rr[24:end] .= 0.0
        dy = zeros(gcSolver.Nspecies)

        gcSolver.fullDeriv(dy, copy(dy), rr, 0.0)

        # Check for conservation of each surface receptor
        assertConservation(dy)
        # Check for conservation of each endosomal receptor
        assertConservation(dy[(gcSolver.halfL + 1):end])
    end

    @testset "Reproducibility." begin
        output = runCkine(tps, rxntfR)

        @test ndims(output) == 2
        @test output == runCkine(tps, rxntfR)
    end

    @testset "Steady-state at t=0." begin
        rr = copy(rxntfR)
        rr[1:3] .= 0.0

        out = gcSolver.runCkine([0.0, 1.0], rr)

        @test all(out .>= 0.0)
        @test norm(out[1, :] - out[2, :]) < 1.0e-4
    end

    @testset "Equilibrium." begin
        out = runCkine([100000.0], rxntfR)

        dy = ones(gcSolver.Nspecies)

        gcSolver.fullDeriv(dy, vec(out), rxntfR, 0.0)

        @test all(out .>= 0.0)

        @test norm(dy) < 1.0e-6
    end

    @testset "Detailed balance using no-trafficking model." begin
        prob = gcSolver.runCkineSetup([1000000.0], deepcopy(rxntfR))
        prob.p[20:30] .= 0.0 # Set everything outside of binding reactions to 0

        out = vec(runCkine([1000000.0], prob))

        J = ForwardDiff.jacobian((y, x) -> gcSolver.fullDeriv(y, x, prob.p, 0.0), ones(gcSolver.Nspecies), out)

        # Slice out just the surface species
        GK = J[1:(2gcSolver.halfL), 1:(2gcSolver.halfL)] * diagm(vec(out[1:(2gcSolver.halfL)]))

        @test norm(GK - transpose(GK)) < 1.0e-8
    end

    @testset "Test that there is at least 1 non-zero species at T=0." begin
        temp = runCkine(tps, rxntfR)
        @test any(temp[1, :] .> 0.0)
    end
end


text:
function get_index(quad::Quadtree, shape::Shape)
    x_max = maximum(shape.bounds[:,1])
    y_max = maximum(shape.bounds[:,2])
    x_min = minimum(shape.bounds[:,1])
    y_min = minimum(shape.bounds[:,2])

    left_quadrants = x_min > quad.x0 && x_max < (quad.x1 + quad.x0)/2
    right_quadrants = x_min > (quad.x1 + quad.x0)/2 && x_max < quad.x1

    top_quadrants = y_min > quad.y0 && y_max < (quad.y1 + quad.y0)/2
    bottom_quadrants = y_min > (quad.y1 + quad.y0)/2 && y_max < quad.y1

    if(right_quadrants && top_quadrants)
        return 1
    end

    if(left_quadrants && top_quadrants)
        return 2
    end

    if(left_quadrants && bottom_quadrants)
        return 3
    end

    if(right_quadrants && bottom_quadrants)
        return 4
    end

    return -1
end


text:
@testset "Statistics" begin
  @testset "Solution" begin
    grid2D = RegularGrid(3,3)
    realizations = Dict(:value => [i*ones(nelms(grid2D)) for i in 1:3])
    solution2D = SimulationSolution(grid2D, realizations)

    # mean
    mean2D = mean(solution2D)
    @test mean2D[:value] == 2.0*ones(nelms(mean2D))
    @test domain(mean2D) == solution2D.domain

    # variance
    variance2D = var(solution2D)
    @test variance2D[:value] == 1.0*ones(nelms(variance2D))
    @test domain(variance2D) == solution2D.domain

    # quantile (scalar)
    p = 0.5
    quantile2D = quantile(solution2D, p)
    @test quantile2D[:value] == 2.0*ones(nelms(quantile2D))
    @test domain(quantile2D) == solution2D.domain

    # quantile (vector)
    ps = [0.0, 0.5, 1.0]
    quantiles2D = quantile(solution2D, ps)
    @test quantiles2D[2][:value] == quantile2D[:value]
  end

  @testset "Data" begin
    # load data with bias towards large values (gold mine)
    sdata = readgeotable(joinpath(datadir,"clustered.csv"), coordnames=(:x,:y))

    # spatial mean
    μn = mean(sdata[:Au])
    μs = mean(sdata, :Au)
    @test abs(μn - 0.5) > abs(μs - 0.5)
    @test mean(sdata)[:Au] ≈ μs

    # spatial variance
    σn = var(sdata[:Au])
    σs = var(sdata, :Au)
    @test σn ≤ σs
    @test var(sdata)[:Au] ≈ σs

    # spatial quantile
    qn = quantile(sdata[:Au], 0.5)
    qs = quantile(sdata, :Au, 0.5)
    @test qn ≥ qs
    @test quantile(sdata, 0.5)[:Au] ≈ qs
  end
end


text:
# QCQP example
#    maximize x + y
#
#    s.t.  x, y >= 0
#          x^2 + y^2 <= 1
#
#    solution: (0.71, 0.71) objv = 1.414

using Gurobi

env = Gurobi.Env()

model = Gurobi.Model(env, "qcqp_01", :maximize)

add_cvars!(model, [1., 1.], 0., Inf)
update_model!(model)

 # add_qpterms!(model, linearindices, linearcoeffs, qrowinds, qcolinds, qcoeffs, sense, rhs)
add_qconstr!(model, [], [], [1, 2], [1, 2], [1, 1.], '<', 1.0)
update_model!(model)

println(model)

optimize(model)

println("sol = $(get_solution(model))")
println("obj = $(get_objval(model))")


text:
using JetPack, Jets, LinearAlgebra, Test

@testset "JopDifference, linearity, n=$n, T=$T" for n in ( (11,), (11,22), (11,22,33) ), T in (Float32,Float64)
    A = JopDifference(JetSpace(T,n))
    lhs, rhs = linearity_test(A)
    @test lhs ≈ rhs
end

@testset "JopDifference, dot-product, n=$n, T=$T" for n in ( (11,), (11,22), (11,22,33) ), T in (Float32,Float64)
    A = JopDifference(JetSpace(T,n))
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopDifference, correctness, n=$n" for n in ( (11,), (11,12), (11,12,13) )
    ndims = length(n)
    ops = [JopDifference(JetSpace(Float64,n), i) for i=1:ndims]
    α = [rand(Float64) for i=1:ndims]
    x = [mapreduce(i->α[i]*I.I[i], +, 1:length(n)) for I in CartesianIndices(n)]
    y = [A*x for A in ops]
    δ = [abs(sum(y[i] .- α[i]))/length(y[i]) for i=1:ndims]
    map(i->@test(δ[i]<0.1), 1:length(n))
end


text:
# MD5 HMAC tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/hmac-test.c
for (key,text,true_digest) in [
    (
        "",
        "",
        "74e6f7298a9c2d168935f58c001bad88"
    ),(
        "key",
        "The quick brown fox jumps over the lazy dog",
        "80070713463e7749b90c2dc24911e275"
    ),(
        hex2bytes(
        "0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b"),
        "Hi There",
        "9294727a3638bb1c13f48ef8158bfc9d"
    ),(
        "Jefe",
        "what do ya want for nothing?",
        "750c783e6ab0b503eaa86e310a5db738"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        hex2bytes(
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddd"),
        "56be34521d144c88dbb8c733f0e8b3f6"
    ),(
        hex2bytes(
        "0102030405060708090a0b0c0d0e0f10" *
        "111213141516171819"),
        hex2bytes(
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcd"),
        "697eaf0aca3a3aea3a75164746ffaa79"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key - Hash Key First",
        "6b1ab7fe4bd7bf8f0b62e6ce61b9d0cd"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key and Larger " *
        "Than One Block-Size Data",
        "6f630fad67cda0ee1fb1f562db3aa53e"
    )
]
    @test digest("md5", key, text) == hex2bytes(true_digest)
end

# SHA1 HMAC tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/hmac-test.c
for (key,text,true_digest) in [
    (
        "",
        "",
        "fbdb1d1b18aa6c08324b7d64b71fb76370690e1d"
    ),(
        "key",
        "The quick brown fox jumps over the lazy dog",
        "de7c9b85b8b78aa6bc8a7a36f70a90701c9db4d9"
    ),(
        hex2bytes(
        "0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b"),
        "Hi There",
        "b617318655057264e28bc0b6fb378c8ef146be00"
    ),(
        "Jefe",
        "what do ya want for nothing?",
        "effcdf6ae5eb2fa2d27416d5f184df9c259a7c79"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        hex2bytes(
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddd"),
        "125d7342b9ac11cd91a39af48aa17b4f63f175d3"
    ),(
        hex2bytes(
        "0102030405060708090a0b0c0d0e0f10" *
        "111213141516171819"),
        hex2bytes(
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcd"),
        "4c9007f4026250c6bc8414f9bf50c86c2d7235da"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key - Hash Key First",
        "aa4ae5e15272d00e95705637ce8a3b55ed402112"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"),
        "Test Using Larger Than Block-Size Key and Larger " *
        "Than One Block-Size Data",
        "e8e99d0f45237d786d6bbaa7965c7808bbff1a91"
    )
]
    @test hexdigest("sha1", key, text) == true_digest
end

# SHA256 HMAC tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/hmac-test.c
for (key,text,true_digest) in [
   (
        "",
        "",
        "b613679a0814d9ec772f95d778c35fc5ff1697c493715653c6c712144292c5ad"
    ),(
        "key",
        "The quick brown fox jumps over the lazy dog",
        "f7bc83f430538424b13298e6aa6fb143ef4d59a14946175997479dbc2d1a3cd8"
    ),(
        hex2bytes(
        "0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b0b" *
        "0b0b0b0b"),
        "Hi There",
        "b0344c61d8db38535ca8afceaf0bf12b" *
        "881dc200c9833da726e9376c2e32cff7"
    ),(
        "Jefe",
        "what do ya want for nothing?",
        "5bdcc146bf60754e6a042426089575c7" *
        "5a003f089d2739839dec58b964ec3843"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaa"),
        hex2bytes(
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddddddddddddddddddddddddddddddd" *
        "dddd"),
        "773ea91e36800e46854db8ebd09181a7" *
        "2959098b3ef8c122d9635514ced565fe"
    ),(
        hex2bytes(
        "0102030405060708090a0b0c0d0e0f10" *
        "111213141516171819"),
        hex2bytes(
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcdcdcdcdcdcdcdcdcdcdcdcdcdcdcd" *
        "cdcd"),
        "82558a389a443c0ea4cc819899f2083a" *
        "85f0faa3e578f8077a2e3ff46729665b"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaa"),
        "Test Using Larger Than Block-Size Key - Hash Key First",
        "60e431591ee0b67f0d8a26aacbf5b77f" *
        "8e0bc6213728c5140546040f0ee37f54"
    ),(
        hex2bytes(
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa" *
        "aaaaaa"),
        "This is a test using a larger than block-size ke" *
        "y and a larger than block-size data. The key nee" *
        "ds to be hashed before being used by the HMAC al" *
        "gorithm.",
        "9b09ffa71b942fcb27635fbcd5b0e944" *
        "bfdc63644f0713938a7f51535c3a35e2"
    )
]
    local h = HMACState("sha256", key)
    update!(h, text)
    @test hexdigest!(h) == true_digest
end

# Test invalid parameters
@test_throws ArgumentError HMACState("this is not a cipher", "")

# Test show methods
println("Testing HMAC show methods:")
println(HMACState("SHA256", ""))


text:
"""
Provides a domain specific language for representing HTML documents.

# Examples

```julia
using Documenter.Utilities.DOM

# `DOM` does not export any HTML tags. Define the ones we actually need.
@tags div p em strong ul li

div(
    p("This ", em("is"), " a ", strong("paragraph."),
    p("And this is ", strong("another"), " one"),
    ul(
        li("and"),
        li("an"),
        li("unordered"),
        li("list")
    )
)
```

*Notes*

All the arguments passed to a node are flattened into a single vector rather
than preserving any nested structure. This means that passing two vectors of
nodes to a `div` will result in a `div` node with a single vector of children
(the concatenation of the two vectors) rather than two vector children. The
only arguments that are not flattened are nested nodes.

String arguments are automatically converted into text nodes. Text nodes do not
have any children or attributes and when displayed the string is escaped using
[`escapehtml`](@ref).

# Attributes

As well as plain nodes shown in the previous example, nodes can have attributes
added to them using the following syntax.

```julia
div[".my-class"](
    img[:src => "foo.jpg"],
    input[\"#my-id\", :disabled]
)
```

In the above example we add a `class = "my-class"` attribute to the `div` node,
a `src = "foo.jpg"` to the `img`, and `id = "my-id" disabled` attributes to the
`input` node.

The following syntax is supported within `[...]`:

```julia
tag[\"#id\"]
tag[".class"]
tag[\".class#id\"]
tag[:disabled]
tag[:src => "foo.jpg"]
# ... or any combination of the above arguments.
```

# Internal Representation

The [`@tags`](@ref) macro defines named [`Tag`](@ref) objects as follows

```julia
@tags div p em strong
```

expands to

```julia
const div, p, em, strong = Tag(:div), Tag(:p), Tag(:em), Tag(:strong)
```

These [`Tag`](@ref) objects are lightweight representations of empty HTML
elements without any attributes and cannot be used to represent a complete
document. To create an actual tree of HTML elements that can be rendered we
need to add some attributes and/or child elements using `getindex` or `call`
syntax. Applying either to a [`Tag`](@ref) object will construct a new
[`Node`](@ref) object.

```julia
tag(...)      # No attributes.
tag[...]      # No children.
tag[...](...) # Has both attributes and children.
```

All three of the above syntaxes return a new [`Node`](@ref) object. Printing of
`Node` objects is defined using the standard Julia display functions, so only
needs a call to `print` to print out a valid HTML document with all nessesary
text escaped.
"""
module DOM

import ..Utilities

using Compat

tostr(p::Pair) = p

export @tags

#
# The following sets are based on:
#
# - https://developer.mozilla.org/en/docs/Web/HTML/Block-level_elements
# - https://developer.mozilla.org/en-US/docs/Web/HTML/Inline_elements
# - https://developer.mozilla.org/en-US/docs/Glossary/empty_element
#
const BLOCK_ELEMENTS = Set([
    :address, :article, :aside, :blockquote, :canvas, :dd, :div, :dl,
    :fieldset, :figcaption, :figure, :footer, :form, :h1, :h2, :h3, :h4, :h5,
    :h6, :header, :hgroup, :hr, :li, :main, :nav, :noscript, :ol, :output, :p,
    :pre, :section, :table, :tfoot, :ul, :video,
])
const INLINE_ELEMENTS = Set([
    :a, :abbr, :acronym, :b, :bdo, :big, :br, :button, :cite, :code, :dfn, :em,
    :i, :img, :input, :kbd, :label, :map, :object, :q, :samp, :script, :select,
    :small, :span, :strong, :sub, :sup, :textarea, :time, :tt, :var,
])
const VOID_ELEMENTS = Set([
    :area, :base, :br, :col, :command, :embed, :hr, :img, :input, :keygen,
    :link, :meta, :param, :source, :track, :wbr,
])
const ALL_ELEMENTS = union(BLOCK_ELEMENTS, INLINE_ELEMENTS, VOID_ELEMENTS)

#
# Empty string as a constant to make equality checks slightly cheaper.
#
const EMPTY_STRING = ""
const TEXT = Symbol(EMPTY_STRING)

"""
Represents a empty and attribute-less HTML element.

Use [`@tags`](@ref) to define instances of this type rather than manually
creating them via `Tag(:tagname)`.
"""
struct Tag
    name :: Symbol
end

Base.show(io::IO, t::Tag) = print(io, "<", t.name, ">")

"""
Define a collection of [`Tag`](@ref) objects and bind them to constants
with the same names.

# Examples

Defined globally within a module:

```julia
@tags div ul li
```

Defined within the scope of a function to avoid cluttering the global namespace:

```julia
function template(args...)
    @tags div ul li
    # ...
end
```
"""
macro tags(args...) esc(tags(args)) end
tags(s) = :(($(s...),) = $(map(Tag, s)))

const Attributes = Vector{Pair{Symbol, String}}

"""
Represents an element within an HTML document including any textual content,
children `Node`s, and attributes.

This type should not be constructed directly, but instead via `(...)` and
`[...]` applied to a [`Tag`](@ref) or another [`Node`](@ref) object.
"""
struct Node
    name :: Symbol
    text :: String
    attributes :: Attributes
    nodes :: Vector{Node}

    Node(name::Symbol, attr::Attributes, data::Vector{Node}) = new(name, EMPTY_STRING, attr, data)
    Node(text::AbstractString) = new(TEXT, text)
end

#
# Syntax for defining `Node` objects from `Tag`s and other `Node` objects.
#
(t::Tag)(args...) = Node(t.name, Attributes(), data(args))
(n::Node)(args...) = Node(n.name, n.attributes, data(args))
Base.getindex(t::Tag, args...) = Node(t.name, attr(args), Node[])
Base.getindex(n::Node, args...) = Node(n.name, attr(args), n.nodes)

#
# Helper methods for the above `Node` "pseudo-constructors".
#
data(args) = flatten!(nodes!, Node[], args)
attr(args) = flatten!(attributes!, Attributes(), args)

#
# Types that must not be flattened when constructing a `Node`'s child vector.
#
const Atom = Union{AbstractString, Node, Pair, Symbol}

"""
# Signatures

```julia
flatten!(f!, out, x::Atom)
flatten!(f!, out, xs)
```

Flatten the contents the third argument into the second after applying the
function `f!` to the element.
"""
flatten!(f!, out, x::Atom) = f!(out, x)
flatten!(f!, out, xs)      = (for x in xs; flatten!(f!, out, x); end; out)

#
# Helper methods for handling flattening children elements in `Node` construction.
#
nodes!(out, s::AbstractString) = push!(out, Node(s))
nodes!(out, n::Node)           = push!(out, n)

#
# Helper methods for handling flattening in construction of attribute vectors.
#
function attributes!(out, s::AbstractString)
    class, id = IOBuffer(), IOBuffer()
    for x in eachmatch(r"[#|\.]([\w\-]+)", s)
        print(startswith(x.match, '.') ? class : id, x.captures[1], ' ')
    end
    position(class) === 0 || push!(out, tostr(:class => rstrip(String(take!(class)))))
    position(id)    === 0 || push!(out, tostr(:id    => rstrip(String(take!(id)))))
    return out
end
attributes!(out, s::Symbol) = push!(out, tostr(s => ""))
attributes!(out, p::Pair)   = push!(out, tostr(p))

function Base.show(io::IO, n::Node)
    if n.name === Symbol("#RAW#")
        print(io, n.nodes[1].text)
    elseif n.name === TEXT
        print(io, escapehtml(n.text))
    else
        print(io, '<', n.name)
        for (name, value) in n.attributes
            print(io, ' ', name)
            isempty(value) || print(io, '=', repr(escapehtml(value)))
        end
        if n.name in VOID_ELEMENTS
            print(io, "/>")
        else
            print(io, '>')
            if n.name === :script || n.name === :style
                isempty(n.nodes) || print(io, n.nodes[1].text)
            else
                for each in n.nodes
                    show(io, each)
                end
            end
            print(io, "</", n.name, '>')
        end
    end
end

Base.show(io::IO, ::MIME"text/html", n::Node) = print(io, n)

"""
Escape characters in the provided string. This converts the following characters:

- `<` to `<`
- `>` to `>`
- `&` to `&`
- `'` to `'`
- `\"` to `"`

When no escaping is needed then the same object is returned, otherwise a new
string is constructed with the characters escaped. The returned object should
always be treated as an immutable copy and compared using `==` rather than `===`.
"""
function escapehtml(text::AbstractString)
    if occursin(r"[<>&'\"]", text)
        buffer = IOBuffer()
        for char in text
            char === '<'  ? write(buffer, "<")   :
            char === '>'  ? write(buffer, ">")   :
            char === '&'  ? write(buffer, "&")  :
            char === '\'' ? write(buffer, "'")  :
            char === '"'  ? write(buffer, """) : write(buffer, char)
        end
        String(take!(buffer))
    else
        text
    end
end

"""
A HTML node that wraps around the root node of the document and adds a DOCTYPE
to it.
"""
mutable struct HTMLDocument
    doctype :: String
    root    :: Node
end
HTMLDocument(root) = HTMLDocument("html", root)

function Base.show(io::IO, doc::HTMLDocument)
    println(io, "<!DOCTYPE $(doc.doctype)>")
    println(io, doc.root)
end

end


text:
module Generate

import Base.Pkg.PkgError
importall Base.LibGit2

copyright_year() = string(Dates.year(Dates.today()))
copyright_name(repo::GitRepo) = LibGit2.getconfig(repo, "user.name", "")

function package(pkg::AbstractString, description::AbstractString;
                    authors::Union{AbstractString, Array} = "",
                    authors_url::AbstractString          = "",
                    path::AbstractString                 = Pkg.Dir.path(),
                    github_name::AbstractString          = "",
                    github_ssh::Bool                     = true,
                    mkdocs::Bool                         = true
                )

    pkg_path = joinpath(path, pkg)
    ispath(pkg_path) && throw(PkgError("$pkg seems to exit: Path <$pkg_path>  To remove it use: `Pkg.rm`"))
    isempty(description) && throw(PkgError("`description` may not be empty."))

    year = copyright_year()
    try
        repo = Generate.init(pkg_path, github_name, github_ssh)
        LibGit2.transact(repo) do repo
            if isempty(authors)
                authors = copyright_name(repo)
            end

            if isempty(authors_url)
                authors_url = "https://github.com/$(copyright_name(repo))/"
            end

            if isempty(github_name)
                github_name = LibGit2.getconfig(repo, "user.name", "")
            end

            files = [Generate.license(pkg_path, year, authors, authors_url),
                     Generate.readme(pkg_path, description, github_name),
                     Generate.tests(pkg_path),
                     Generate.require(pkg_path),
                     Generate.gitignore(pkg_path),
                     Generate.entrypoint(pkg_path, description),
                     Generate.todo(pkg_path),
                     Generate.changelog(pkg_path),
                    ]
            mkdocs && push!(files, Generate.mkdocs(pkg_path, description, year, authors, authors_url, github_name))

            msg = """
            $pkg.jl "generated" files.

                authors:  $(join(vcat(authors),", "))
                authors info url: $authors_url
                year:    $year

            Julia Version $VERSION [$(Base.GIT_VERSION_INFO.commit_short)]
            """
            LibGit2.add!(repo, files..., flags = LibGit2.Consts.INDEX_ADD_FORCE)
            info("Committing $pkg generated files")
            LibGit2.commit(repo, msg)
        end
    catch
        rm(pkg_path, recursive=true)
        rethrow()
    end
    return
end

function init(pkg::AbstractString, github_name::AbstractString, github_ssh::Bool)
    ispath(pkg) && throw(PkgError("$pkg seems to exit: Path <$pkg_path"))

    pkg_name = basename(pkg)
    info("Initializing $pkg_name repo: $pkg")
    repo = LibGit2.init(pkg)
    try
        LibGit2.commit(repo, "initial empty commit")
    catch err
        throw(PkgError("Unable to initialize $pkg_name package: $err"))
    end

    if isempty(github_name)
        github_name = LibGit2.getconfig(repo, "user.name", "")
    end

    try
        if github_name != "NONE"
            if github_ssh
                repo_url = "git@github.com:$(github_name)/$(pkg_name).jl.git"
            else
                repo_url = "https://github.com/$(github_name)/$(pkg_name).jl/"
            end
            info("Origin: $repo_url")
            with(LibGit2.GitRemote, repo, "origin", repo_url) do rmt
                LibGit2.save(rmt)
            end
            LibGit2.set_remote_url(repo, repo_url)
        end
    end
    return repo
end

function genfile(f::Function, pkg::AbstractString, file::AbstractString)
    path = joinpath(pkg, file)
    if !ispath(path)
        info("Generating $file")
        mkpath(dirname(path))
        open(f, path, "w")
        return file
    end
    return ""
end

function copyright(year::AbstractString, authors::Array, authors_url::AbstractString)
    text = "> Copyright (c) $year - $year:"
    for author in authors
        text *= "\n>  * $author"
    end
    if authors_url != "NONE"
        text *= "\n>"
        text *= "\n> Authors info: $authors_url"
        text *= "\n>"
    end
    return text
end

function copyright(year::AbstractString, authors::AbstractString, authors_url::AbstractString)
    if authors_url != "NONE"
        return "> Copyright (c) $year - $year: **$authors**.   (<$authors_url>)"
    else
        return "> Copyright (c) $year - $year: **$authors**."
    end
end

function license(pkg::AbstractString,
                 year::AbstractString,
                 authors::Union{AbstractString,Array},
                 authors_url::AbstractString)
    pkg_name = basename(pkg)
    genfile(pkg,"LICENSE.md") do io
        println(io, "## Copyrights & Licenses")
        println(io)
        println(io, "The *$pkg_name.jl package* is licensed under the MIT \"Expat\" License:")
        println(io)
        println(io, copyright(year, authors, authors_url))
        lic=readall(normpath(dirname(@__FILE__), "..", "res", "licenses_MIT"))
        for l in split(lic,['\n','\r'])
            println(io, "> ", l)
        end
        println(io)
        println(io, "### Licenses for incorporated software")
        println(io)
        println(io, "The *$pkg_name.jl package* contains some code derived from the following sources, ",
                           "which have their own licenses:")
        println(io)
    end
end

function readme(pkg::AbstractString,
                description::AbstractString,
                github_name::AbstractString
                )
    pkg_name = basename(pkg)
    genfile(pkg,"README.md") do io
        println(io, "## $(pkg_name).jl")
        println(io)
        println(io, "$description")
        println(io)
        if github_name != "NONE"
            println(io, "### Web Presence")
            println(io)
            println(io, "* $(pkg_name).jl [web site and documentation]",
                        "(http://$(lowercase(github_name)).github.io/$(pkg_name).jl/)")
            println(io, "* $(pkg_name).jl [github repository](https://github.com/$(github_name)/$(pkg_name).jl/)")
            println(io)
        end
    end
end

function tests(pkg::AbstractString)
    pkg_name = basename(pkg)
    genfile(pkg,"test/runtests.jl") do io
        print(io, """
        using $pkg_name
        using Base.Test

        # write your own tests here
        @test 1 == 1
        """)
    end
end

function versionfloor(ver::VersionNumber)
    # return "major.minor" for the most recent release version relative to ver
    # for prereleases with ver.minor == ver.patch == 0, return "major-" since we
    # don't know what the most recent minor version is for the previous major
    if isempty(ver.prerelease) || ver.patch > 0
        return string(ver.major, '.', ver.minor)
    elseif ver.minor > 0
        return string(ver.major, '.', ver.minor - 1)
    else
        return string(ver.major, '-')
    end
end

function todo(pkg::AbstractString)
    genfile(pkg,"TODO.md") do io
        print(io, "## TODO")
        print(io)
    end
end

function require(pkg::AbstractString)
    genfile(pkg,"REQUIRE") do io
        print(io, """
        julia $(versionfloor(VERSION))
        """)
    end
end

function gitignore(pkg::AbstractString)
    genfile(pkg,".gitignore") do io
        print(io, """
        *.jl.cov
        *.jl.*.cov
        *.jl.mem

        # MkDocs build
        site/
        """)
    end
end

function entrypoint(pkg::AbstractString, description::AbstractString)
    pkg_name = basename(pkg)
    genfile(pkg,"src/$pkg_name.jl") do io
        print(io, """#__precompile__(true)

            "$description"
            module $pkg_name

            # package code goes here

            end # module
            """)
    end
end

function mkdocs_yml(pkg::AbstractString,
                    description::AbstractString,
                    year::AbstractString,
                    authors::Union{AbstractString, Array},
                    authors_url::AbstractString,
                    github_name::AbstractString
                    )
    pkg_name = basename(pkg)
    genfile(pkg,"mkdocs.yml") do io
        println(io, "site_name:        $(pkg_name).jl")
        github_name != "NONE" && println(io, "site_url:         https://$(github_name).github.io/$(pkg_name).jl/")
        println(io, "site_description: $(description)")
        println(io, "site_author:      $(authors)")
        println(io, "site_favicon:     images/$(lowercase(pkg_name))_favicon.ico")
        println(io)
        github_name != "NONE" && println(io, "repo_url:         https://github.com/$(github_name)/$(pkg_name).jl/")
        println(io)
        println(io, "pages:")
        println(io, "- Home: 'index.md'")
        println(io)
        println(io, "- About:")
        println(io, "    - Readme: about/readme.md")
        println(io, "    - License: about/license.md")
        println(io, "    - Changelog: about/changelog.md")
        println(io)
        println(io, "use_directory_urls: true")
        println(io, "copyright: Copyright (c) $(year) - $(year), $(authors) ",
                    "$(authors_url).")
        println(io)
        println(io, "theme: 'bootstrap'")
        println(io)
    end
end

function docs_index(pkg::AbstractString,
                    description::AbstractString,
                    year::AbstractString,
                    authors::Union{AbstractString, Array},
                    authors_url::AbstractString,
                    github_name::AbstractString
                    )
    pkg_name = basename(pkg)
    genfile(pkg,"docs/index.md") do io
        println(io, "# $(pkg_name).jl")
        println(io)
        println(io, "$(description)")
        println(io)
        if github_name != "NONE"
            println(io, "For *official* information visit the *$(pkg_name).jl package* [web site and documentation]",
                        "(http://$(lowercase(github_name)).github.io/$(pkg_name).jl/).")
            println(io)
        end
    end
end

function mkdocs(pkg::AbstractString,
                    description::AbstractString,
                    year::AbstractString,
                    authors::Union{AbstractString, Array},
                    authors_url::AbstractString,
                    github_name::AbstractString
                    )
    # dirs and links
    images_path = joinpath(pkg, "docs", "images")
    mkpath(images_path)
    touch(joinpath(images_path, ".gitkeep"))

    about_path = joinpath(pkg, "docs", "about")
    mkpath(joinpath(about_path, "readme", "images"))
    touch(joinpath(about_path, "readme", "images", ".gitkeep"))

    cd(about_path)
    symlink(joinpath(pkg, "README.md"), "readme.md")
    symlink(joinpath(pkg, "LICENSE.md"), "license.md")
    symlink(joinpath(pkg, "CHANGELOG.md"), "changelog.md")
    cd(normpath(dirname(@__FILE__)))

    #
    mkdocs_yml(pkg, description, year, authors, authors_url, github_name)
    docs_index(pkg, description, year, authors, authors_url, github_name)
end

function changelog(pkg::AbstractString)
    pkg_name = basename(pkg)
    genfile(pkg,"CHANGELOG.md") do io
        print(io, """
            ## v0.1.0 (xxxx-xx-xx)

            ### Summary

            -----------------------------------------------------------------------------------------------------------------------

            ## P-Versioning Based On [Semantic Versioning](http://semver.org/)

            **IMPORTANT DIFFERENCE** to the *Semantic Versioning 2.0.0*

            * A pre-release version MUST NOT be added.

            * Build metadata MUST comprise only ASCII alphanumerics [0-9A-Za-z] and MUST NOT contain any hyphen.

            ### Package Versioning

            1. **Software and related packages using this modified Semantic Versioning MUST declare a public API.** This API could
                be declared in the code itself or exist strictly in documentation. However it is done, it should be precise and
                comprehensive.

            2. A normal version number MUST take the form X.Y.Z where X, Y, and Z are non-negative integers, and MUST NOT contain
                leading zeroes. X is the major version, Y is the minor version, and Z is the patch version.
                Each element MUST increase numerically. For instance: 1.9.0 -> 1.10.0 -> 1.11.0.

            3. Once a versioned package has been released, the contents of that version MUST NOT be modified. Any modifications
                MUST be released as a new version.

            4. Major version zero (0.y.z) is for initial development. Anything may change at any time. The public API should not be
                considered stable.

            5. **Version 1.0.0 defines the public API. The way in which the version number is incremented after this release is
                dependent on this public API and how it changes.**

            6. Patch version Z (x.y.Z | x > 0) MUST be incremented if only backwards compatible bug fixes are introduced. A bug fix
                is defined as an internal change that fixes incorrect behavior.

            7. Minor version Y (x.Y.z | x > 0) MUST be incremented if new, backwards compatible functionality is introduced to the
                public API. It MUST be incremented if any public API functionality is marked as deprecated. It MAY be incremented
                if substantial new functionality or improvements are introduced within the private code. It MAY include patch level
                changes. Patch version MUST be reset to 0 when minor version is incremented.

            8. Major version X (X.y.z | X > 0) MUST be incremented if any backwards incompatible changes are introduced to the
                public API. It MAY include minor and patch level changes. Patch and minor version MUST be reset to 0 when major
                version is incremented.

            9. Build metadata MAY be denoted by appending a plus sign and a series of dot separated identifiers immediately
                following the patch version. Identifiers MUST comprise only ASCII alphanumerics [0-9A-Za-z] and MUST NOT contain
                hyphen. Identifiers MUST NOT be empty. Build metadata SHOULD be ignored when determining version precedence. Thus
                two versions that differ only in the build metadata, have the same precedence.

                Examples: 1.0.0+001, 1.0.0+20130313144700, 1.0.0+exp.sha.5114f85, 1.0.7+r128.g4560914.

            * What do I do if I accidentally release a backwards incompatible change as a minor version?

                As soon as you realize that you've broken the Semantic Versioning spec, fix the problem and release a new minor
                version that corrects the problem and restores backwards compatibility. Even under this circumstance, it is
                unacceptable to modify versioned releases. If it's appropriate, document the offending version and inform your
                users of the problem so that they are aware of the offending version.

            * How should I handle deprecating functionality?

                Deprecating existing functionality is a normal part of software development and is often required to make forward
                progress. When you deprecate part of your public API, you should do two things:

                1. update your documentation to let users know about the change,
                2. issue a new minor release with the deprecation in place. Before you completely remove the functionality in a new
                    major release there should be at least one minor release that contains the deprecation so that users can
                    smoothly transition to the new API.

            """)
    end
end

end # module


text:
struct BlocksInfoF90
    nblocks    :: Cint
    nelectrons :: Cint
    ncsfstotal :: Cint
    norbitals  :: Cint
    nvectotal  :: Cint
    nvecsize   :: Cint
end

struct BlockF90
    blockid       :: Cint
    ncsfs         :: Cint
    nevs          :: Cint
    iatjp         :: Cint
    iaspa         :: Cint
    eav           :: Cdouble
    eigenstates   :: Ptr{Cdouble}
    eigenenergies :: Ptr{Cdouble}
end

struct MixingBlock
    eav :: Float64
    energies :: Vector{Float64}
    states :: Matrix{Float64}
    _block_f90 :: BlockF90
end

struct MixingFile
    nelectrons :: Int
    blocks :: Vector{MixingBlock}
    _blocksinfo_f90 :: BlocksInfoF90
end

function read_rmix(filename)
    isfile(filename) || error("Unable to open $(filename)")

    blockinfo = Ref{BlocksInfoF90}()
    blocks = Ref{Ptr{BlockF90}}()
    status = ccall( (:mixread, libgrasp_so),
        Cint, (Cstring, Ref{BlocksInfoF90}, Ref{Ptr{BlockF90}}),
        filename, blockinfo, blocks
    )
    if status != zero(status)
        error("mixread returned with status $(status)")
    end
    blockinfo = blockinfo.x
    blocks = unsafe_wrap(Array{BlockF90}, blocks.x, blockinfo.nblocks, own=true)
    blocks = map(blocks) do block
        energies = unsafe_wrap(Array{Cdouble}, block.eigenenergies, (block.nevs,), own=true)
        states = unsafe_wrap(Array{Cdouble}, block.eigenstates, (block.ncsfs, block.nevs), own=true)
        MixingBlock(block.eav, energies, states, block)
    end
    MixingFile(blockinfo.nelectrons, blocks, blockinfo)
end


text:
function init_cond(type, N; dx=1/143, M = 5, α=500.0, β=1000.0, c01 = 1, BC=:Periodic, r₀ = 0.0, n_id = 0, n_vd = 0, γ = 50.0, f = 0.01)
    # Use odd sized pattern
    @assert M ÷ 2 + 1 == ceil(M/2)
    # Done
    if type==:NoisePatches
        u01 = c01*ones(Float64, N,N)
        n_rect = rand(20:50)
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        dy = dx
        y = collect(-N/2*dy:dy:(N/2-1)*dy) # Global Coordinate
        for i = 1:n_rect
            l = rand(1:N ÷ 5)
            w = rand(1:N ÷ 5)
            sx = rand(1:4*N ÷ 5)
            sy = rand(1:4*N ÷ 5)
            val = rand(Uniform(-c01*f, +c01*f))
            u01[sx:sx+l,sy:sy+w] .= c01 + val
            #val = rand()
            #u01[sx:sx+l,sy:sy+w] .= val
        end
    elseif type==:Square
        @assert M ≥ 3
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        y = collect(-N/2*dx:dx:(N/2-1)*dx) # Global Coordinate
        dy = dx
        len = x[end]/(M-1)*2 # The factor here is to make the whole MxM grid fit into the lattice
        basis = len*hcat([0,1],[1,0]) |> SMatrix{2,2}
        cell = HomogeneousCell([[0.0,0.0]])
        lattice = RegularLattice((M,M), basis, cell)
        u01 = zeros(Float64, N,N)
        for i in 1:M
            for j in 1:M
                ii, jj= (lattice[i,j,1] .- len*(M+1)/2)
                u01 .+= exp.(-α*((x .- ii).^2 .+ (y' .- jj).^2))
            end
        end
    elseif type==:Er2Mn2O7
        @assert M ≥ 5
        # Grid Stuff
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        dy = dx*sqrt(3)/2 # factor of 6/5??
        y = collect(-N/2*dy:dy:(N/2-1)*dy) # Global Coordinate
        #y = [y; y[end]*(1+dx)] # Hack Fix, idk why y sometimes has N-1 elements
        u01 = zeros(Float64, N,N)
        len = x[end]#/(M-3)*2
        len = 1.0
        # Declare Basis
        fbasis = len*hcat([1, 0],
                      [0.5, 0.5*sqrt(3)],
                     ) |> SMatrix{2,2}
        # Atom positions
        cell_vectors_raw1 = [[0.0, 0.0], [0.5, 0.0], [0.0, 0.5]]
        cell_vectors_raw2 = [[0.5, 0.0], [0.5, 0.5], [0.0, 0.5]]

        fcell = InhomogeneousCell([fbasis*vec for vec in cell_vectors_raw1], [fbasis*vec for vec in cell_vectors_raw2]; label = :ermno3)
        lattice = RegularLattice((M,M), fbasis, fcell; label = :hexagonal)

        #u01 = zeros(Float64, N,N)
        u01 = zeros(Float64, N, N)
        for i in 1:M
            for j in 1:M
                for ic in 1:length(cell_vectors_raw1)
                    ii, jj= lattice[i,j,ic,1]
                    #ii -= len*(M ÷ 2 + 1)
                    if BC != :Periodic
                      if isapprox(ii,x[end],atol=dx) || isapprox(ii,x[1],atol=dx) || isapprox(jj,y[1],atol=dx) || isapprox(jj,y[end],atol=dx)
                        #println("on boundary")
                        continue
                      else
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        u01 .+= exp.(-α*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                      end
                    else
                        u01 .+= exp.(-α*((x .- ii).^2 .+ (y' .- jj).^2))
                    end
                end
                for ic in 1:length(cell_vectors_raw2)
                    ii, jj= lattice[i,j,ic,2]
                    #ii -= len*(M ÷ 2 + 1)
                    if BC != :Periodic
                      if isapprox(ii,x[end],atol=dx) || isapprox(ii,x[1],atol=dx) || isapprox(jj,y[1],atol=dx) || isapprox(jj,y[end],atol=dx)
                        #println("on boundary")
                        continue
                      else
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        u01 .+= exp.(-β*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                      end
                    else
                        u01 .+= exp.(-β*((x .- ii).^2 .+ (y' .- jj).^2))
                    end
                end
            end
        end
    elseif type==:ErMnO3
        @assert M ≥ 5
        # Grid Stuff
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        dy = dx*sqrt(3)/2 # factor of 6/5??
        y = collect(-N/2*dy:dy:(N/2-1)*dy) # Global Coordinate
        #y = [y; y[end]*(1+dx)] # Hack Fix, idk why y sometimes has N-1 elements
        u01 = zeros(Float64, N,N)
        len = x[end]/(M-3)*2
        # Declare Basis
        fbasis = len*hcat([0.5, 0.5*sqrt(3)],
                      [0.5, -0.5*sqrt(3)],
                     ) |> SMatrix{2,2}
        # Atom positions
        cell_vectors_raw1 = [[0.0, 0.0], [1/3, 2/3], [2/3, 1/3]]
        cell_vectors_raw2 = [[0, 1/3], [0, 2/3], [1/3, 1/3], [1/3, 0], [2/3, 0], [2/3, 2/3]]

        fcell = InhomogeneousCell([fbasis*vec for vec in cell_vectors_raw1], [fbasis*vec for vec in cell_vectors_raw2]; label = :ermno3)
        lattice = RegularLattice((M,M), fbasis, fcell; label = :hexagonal)

        u01 = c01*ones(Float64, N,N)
        for i in 1:M
            for j in 1:M
                for ic in 1:length(cell_vectors_raw1)
                    ii, jj= lattice[i,j,ic,1]
                    ii -= len*(M ÷ 2 + 1)
                    if BC != :Periodic
                      if isapprox(ii,x[end],atol=dx) || isapprox(ii,x[1],atol=dx) || isapprox(jj,y[1],atol=dx) || isapprox(jj,y[end],atol=dx)
                        #println("on boundary")
                        continue
                      else
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        u01 .+= c01*exp.(-α*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                      end
                    else
                        u01 .+= c01*exp.(-α*((x .- ii).^2 .+ (y' .- jj).^2))
                    end
                end
                for ic in 1:length(cell_vectors_raw2)
                    ii, jj= lattice[i,j,ic,2]
                    ii -= len*(M ÷ 2 + 1)
                    if BC != :Periodic
                      if isapprox(ii,x[end],atol=1.1*dx) || isapprox(ii,x[1],atol=1.1*dx) || isapprox(jj,y[1],atol=1.1*dx) || isapprox(jj,y[end],atol=1.1*dx)
                        #println("on boundary")
                        continue
                      else
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        u01 .+= c01*exp.(-β*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                      end
                    else
                        u01 .+= c01*exp.(-β*((x .- ii).^2 .+ (y' .- jj).^2))
                    end
                end
            end
        end
    elseif type==:ErMnO3_ErOnly
        @assert M ≥ 9
        # Grid Stuff
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        dy = dx*sqrt(3)/2 # factor of 6/5??
        y = collect(-N/2*dy:dy:(N/2-1)*dy) # Global Coordinate
        #y = [y; y[end]*(1+dx)] # Hack Fix, idk why y sometimes has N-1 elements
        u01 = zeros(Float64, N,N)
        len = x[end]/(M-7)*2
        # Declare Basis
        fbasis = len*hcat([0.5, 0.5*sqrt(3)],
                      [0.5, -0.5*sqrt(3)],
                     ) |> SMatrix{2,2}
        # Atom positions
        cell_vectors_raw1 = [[0.0, 0.0], [1/3, 2/3], [2/3, 1/3]]
        cell_vectors_raw2 = [[0, 1/3], [0, 2/3], [1/3, 1/3], [1/3, 0], [2/3, 0], [2/3, 2/3]]

        fcell = InhomogeneousCell([fbasis*vec for vec in cell_vectors_raw1], [fbasis*vec for vec in cell_vectors_raw2]; label = :ermno3)
        lattice = RegularLattice((M,M), fbasis, fcell; label = :hexagonal)

        u01 = c01*ones(Float64, N,N)
        i_vd = rand(1:M, 1, n_vd)
        j_vd = rand(1:M, 1, n_vd)
        ic_vd = rand(1:length(cell_vectors_raw1), 1, n_vd)
        i_id = rand(1:M, 1, n_id)
        j_id = rand(1:M, 1, n_id)
        ic_id = rand(1:length(cell_vectors_raw2), 1, n_id)
        zv = zip(i_vd, j_vd, ic_vd)
        zi = zip(i_id, j_id, ic_id)
        for i in 1:M
            for j in 1:M
                for ic in 1:length(cell_vectors_raw1)
                    if (i, j, ic) ∈ zv
                        continue
                    end
                    ii, jj= lattice[i,j,ic,1]
                    ii -= len*(M ÷ 2 + 1) #+ 0.5*len
                    if BC != :Periodic
                      if isapprox(ii,x[end],atol=1.1*dx) || isapprox(ii,x[1],atol=1.1*dx) || isapprox(jj,y[1],atol=dx) || isapprox(jj,y[end],atol=dx)
                        #println("On Boundary")
                        continue
                      else
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        val = sign(rand(Uniform(-c01*f, +c01*f)))
                        u01 .+= f*c01*val*exp.(-α*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                      end
                      else
                        val = sign(rand(Uniform(-c01*f, +c01*f)))
                        u01 .+= f*c01*val*exp.(-α*((x .- ii).^2 .+ (y' .- jj).^2))
                    end
                end
                for ic in 1:length(cell_vectors_raw2)
                    ii, jj= lattice[i,j,ic,2]
                    ii -= len*(M ÷ 2 + 1)
                    if isapprox(ii,x[end],atol=dx) || isapprox(ii,x[1],atol=dx) || isapprox(jj,y[1],atol=dx) || isapprox(jj,y[end],atol=dx)
                        continue
                    end
                    if (i, j, ic) ∈ zi
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        val = sign(rand(Uniform(-c01*f, +c01*f)))
                        u01 .+= f*c01*val*exp.(-γ*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                    end
                end
            end
        end
    elseif type==:ErMnO3_MnOnly
        @assert M ≥ 9
        # Grid Stuff
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        dy = dx*sqrt(3)/2 # factor of 6/5??
        y = collect(-N/2*dy:dy:(N/2-1)*dy) # Global Coordinate
        #y = [y; y[end]*(1+dx)] # Hack Fix, idk why y sometimes has N-1 elements
        u01 = zeros(Float64, N,N)
        len = x[end]/(M-7)*2
        # Declare Basis
        fbasis = len*hcat([0.5, 0.5*sqrt(3)],
                      [0.5, -0.5*sqrt(3)],
                     ) |> SMatrix{2,2}
        # Atom positions
        cell_vectors_raw1 = [[0.0, 0.0], [1/3, 2/3], [2/3, 1/3]]
        cell_vectors_raw2 = [[0, 1/3], [0, 2/3], [1/3, 1/3], [1/3, 0], [2/3, 0], [2/3, 2/3]]

        fcell = InhomogeneousCell([fbasis*vec for vec in cell_vectors_raw1], [fbasis*vec for vec in cell_vectors_raw2]; label = :ermno3)
        lattice = RegularLattice((M,M), fbasis, fcell; label = :hexagonal)

        u01 = c01*ones(Float64, N,N)
        i_vd = rand(1:M, 1, n_vd)
        j_vd = rand(1:M, 1, n_vd)
        ic_vd = rand(1:length(cell_vectors_raw1), 1, n_vd)
        i_id = rand(1:M, 1, n_id)
        j_id = rand(1:M, 1, n_id)
        ic_id = rand(1:length(cell_vectors_raw2), 1, n_id)
        zv = zip(i_vd, j_vd, ic_vd)
        zi = zip(i_id, j_id, ic_id)
        for i in 1:M
            for j in 1:M
                for ic in 1:length(cell_vectors_raw2)
                    if (i, j, ic) ∈ zv
                        continue
                    end
                    ii, jj= lattice[i,j,ic,2]
                    ii -= len*(M ÷ 2 + 1)
                    if BC != :Periodic
                      if isapprox(ii,x[end],atol=dx) || isapprox(ii,x[1],atol=dx) || isapprox(jj,y[1],atol=dx) || isapprox(jj,y[end],atol=dx)
                        continue
                      else
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        val = sign(rand(Uniform(-c01*f, +c01*f)))
                        u01 .+= f*c01*val*exp.(-α*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                      end
                    else
                        val = sign(rand(Uniform(-c01*f, +c01*f)))
                        u01 .+= f*c01*val*exp.(-α*((x .- ii).^2 .+ (y' .- jj).^2))
                    end
                end
                for ic in 1:length(cell_vectors_raw2)
                    ii, jj= lattice[i,j,ic,2]
                    ii -= len*(M ÷ 2 + 1)
                    if isapprox(ii,x[end],atol=dx) || isapprox(ii,x[1],atol=dx) || isapprox(jj,y[1],atol=dx) || isapprox(jj,y[end],atol=dx)
                        continue
                    end
                    if (i, j, ic) ∈ zi
                        x₀ = r₀*x[end]*rand()
                        y₀ = r₀*x[end]*rand()
                        val = sign(rand(Uniform(-c01*f, +c01*f)))
                        u01 .+= f*c01*val*exp.(-γ*((x .- ii .- x₀).^2 .+ (y' .- jj .- y₀).^2))
                    end
                end
            end
        end
    elseif type==:HoneyComb
        @assert M ≥ 5
        # Grid Stuff
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        dy = dx*sqrt(3)/2*6/5 # factor of 6/5
        y = collect(-N/2*dy:dy:(N/2-1)*dy) # Global Coordinate
        #y = [y; y[end]*(1+dx)] # Hack Fix, idk why y sometimes has N-1 elements
        u01 = zeros(Float64, N,N)
        len = x[end]/(M-3)*2
        # Declare Basis
        fbasis = len*hcat([0.5, 0.5*sqrt(3)],
                      [0.5, -0.5*sqrt(3)],
                     ) |> SMatrix{2,2}
        # Atom positions
        cell_vectors_raw = [[2/3, 1/3], [1/3, 2/3]]
        pos = [fbasis*vec for vec in cell_vectors_raw]
        fcell = HomogeneousCell(pos)
        lattice = RegularLattice((M,M), fbasis, fcell; label = :hexagonal)

        u01 = zeros(Float64, N,N)
        for i in 1:M
            for j in 1:M
                for ic in 1:length(cell_vectors_raw)
                    ii, jj= lattice[i,j,ic]
                    ii -= len*(M ÷ 2 + 1)
                    u01 .+= exp.(-α*((x .- ii).^2 .+ (y' .- jj).^2))
                end
            end
        end
    elseif type==:Hexagonal
        @assert M ≥ 5
        # Grid Stuff
        x = collect(-N/2*dx:dx:(N/2-1)*dx) # Global coordinate
        dy = dx*sqrt(3)/2
        y = collect(-N/2*dy:dy:(N/2-1)*dy) # Global Coordinate
        #y = [y; y[end]*(1+dx)] # Hack Fix, idk why y sometimes has N-1 elements
        u01 = zeros(Float64, N,N)
        len = x[end]/(M-3)*2
        # Declare Basis
        fbasis = len*hcat([0.5, 0.5*sqrt(3)],
                      [0.5, -0.5*sqrt(3)],
                     ) |> SMatrix{2,2}
        # Atom positions
        cell_vectors_raw = [[0,0]]
        pos = [fbasis*vec for vec in cell_vectors_raw]
        fcell = HomogeneousCell(pos)
        lattice = RegularLattice((M,M), fbasis, fcell; label = :hexagonal)

        for i in 1:M
            for j in 1:M
                for ic in 1:length(cell_vectors_raw)
                    ii, jj= lattice[i,j, ic]
                    ii -= len*(M ÷ 2 + 1)
                    u01 .+= exp.(-α*((x .- ii).^2 .+ (y' .- jj).^2))
                end
            end
        end
    end
    return u01, x, y, dx, dy
end


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.

function MOIU.map_indices(index_map::Function, change::DecisionCoefficientChange)
    return DecisionCoefficientChange(index_map(change.decision), change.new_coefficient)
end
MOIU.map_indices(index_map::Function, change::DecisionStateChange) = change
MOIU.map_indices(index_map::Function, change::KnownValuesChange) = change

function MOIU.modify_function(f::AffineDecisionFunction, change::Union{MOI.ScalarConstantChange, MOI.ScalarCoefficientChange})
    return typeof(f)(MOIU.modify_function(f.variable_part, change),
                     copy(f.decision_part))
end

function MOIU.modify_function(f::QuadraticDecisionFunction{T}, change::Union{MOI.ScalarConstantChange, MOI.ScalarCoefficientChange}) where T
    return typeof(f)(
        MOIU.modify_function(f.variable_part, change),
        copy(f.decision_part),
        copy(f.cross_terms))
end

function MOIU.modify_function(f::VectorAffineDecisionFunction, change::Union{MOI.VectorConstantChange, MOI.MultirowChange})
    return typeof(f)(MOIU.modify_function(f.variable_part, change),
                     copy(f.decision_part))
end

function MOIU.modify_function(f::AffineDecisionFunction, change::DecisionCoefficientChange)
    return typeof(f)(copy(f.variable_part),
                     MOIU.modify_function(f.decision_part,
                                          MOI.ScalarCoefficientChange(change.decision, change.new_coefficient)))
end

function MOIU.modify_function(f::QuadraticDecisionFunction{T}, change::DecisionCoefficientChange) where T
    return typeof(f)(
        copy(f.variable_part),
        MOIU.modify_function(f.decision_part,
                             MOI.ScalarCoefficientChange(change.decision, change.new_coefficient)),
        copy(f.cross_terms))
end

function MOIU.modify_function(f::VectorAffineDecisionFunction, change::DecisionMultirowChange)
    return typeof(f)(copy(f.variable_part),
                     MOIU.modify_function(f.decision_part,
                                          MOI.MultirowChange(change.decision, change.new_coefficients)))
end

function MOIU.modify_function(f::Union{SingleDecision, AffineDecisionFunction, QuadraticDecisionFunction, VectorAffineDecisionFunction}, change::Union{DecisionStateChange, KnownValuesChange})
    # Nothing to do here, handled in bridges
    return f
end

# Can rely on scalarize bridge if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.ScalarizeBridge{T, AffineDecisionFunction{T}},
                    change::DecisionModification) where T
    for constraint in bridge.scalar_constraints
        MOI.modify(model, constraint, change)
    end
    return nothing
end

# Can rely on vectorize bridge if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.VectorizeBridge{T, VectorAffineDecisionFunction{T}},
                    change::DecisionModification) where T
    MOI.modify(model, bridge.vector_constraint, change)
    return nothing
end

# Can rely on flipsign bridges if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.FlipSignBridge,
                    change::DecisionModification)
    MOI.modify(model, bridge.constraint, change)
    return nothing
end

# Can rely on norm-bridges if modifications are passed along
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.NormInfinityBridge{T, VectorAffineDecisionFunction{T}},
                    change::DecisionModification) where T
    MOI.modify(model, bridge.constraint, change)
    return nothing
end
function MOI.modify(model::MOI.ModelLike,
                    bridge::MOIB.Constraint.NormOneBridge{T, VectorAffineDecisionFunction{T}},
                    change::DecisionModification) where T
    MOI.modify(model, bridge.nn_index, change)
    return nothing
end

function MOI.modify(uf::MOIU.UniversalFallback, ::CI{F,S}, ::Union{DecisionStateChange,KnownValuesChange}) where {F, S}
    return nothing
end

function MOI.get(model::MOIU.CachingOptimizer,
                 attr::DecisionIndex,
                 ci::CI)
    return MOI.get(model.optimizer, attr, model.model_to_optimizer_map[ci])
end

function MOI.get(b::MOIB.AbstractBridgeOptimizer,
                 attr::DecisionIndex, ci::MOI.ConstraintIndex)
    return MOIB.call_in_context(b, ci, bridge -> MOI.get(b, attr, bridge))
end