text:
# Dictionary which returns (and sets) a default value for a requested item not
# already in to the dictionary

# DefaultDictBase is the main class used to in Default*Dicts.
#
# Each related (immutable) Default*Dict class contains a single
# DefautlDictBase object as a member, and delegates almost all
# functions to this object.
#
# The main rationale for doing this instead of using type aliases is
# that this way, we can have actual class names and constructors for
# each of the DefaultDictBase "subclasses", in some sense getting
# around the Julia limitation of not allowing concrete classes to be
# subclassed.
#

immutable DefaultDictBase{K,V,F,D<:Associative{K,V}} <: Associative{K,V}
    default::F
    d::D

    DefaultDictBase(x::F, kv::AbstractArray{(K,V)}) = new(x, D(kv))
    DefaultDictBase(x::F, d::DefaultDictBase) = DefaultDictBase(x, d.d)
    DefaultDictBase(x::F, d::D=D()) = new(x, d)
    DefaultDictBase(x, ks, vs) = new(x, D(ks,vs))

end

# Constructors

DefaultDictBase() = error("no default specified")
DefaultDictBase(k,v) = error("no default specified")

# TODO: these mimic similar Dict constructors, but may not be needed
DefaultDictBase{K,V,F}(default::F, ks::AbstractArray{K}, vs::AbstractArray{V}) =
    DefaultDictBase{K,V,F,Dict{K,V}}(default,ks,vs)
DefaultDictBase{F}(default::F,ks,vs) = DefaultDictBase{Any,Any,F,Dict}(default, ks, vs)

# syntax entry points
DefaultDictBase{F}(default::F) = DefaultDictBase{Any,Any,F,Dict}(default)
DefaultDictBase{K,V,F}(::Type{K}, ::Type{V}, default::F) = DefaultDictBase{K,V,F,Dict}(default)
DefaultDictBase{K,V,F}(default::F, kv::AbstractArray{(K,V)}) = DefaultDictBase{K,V,F,Dict}(default, kv)

DefaultDictBase{F,D<:Associative}(default::F, d::D) = ((K,V)=eltype(d); DefaultDictBase{K,V,F,D}(default, d))

# Functions

# most functions are simply delegated to the wrapped dictionary
@delegate DefaultDictBase.d [ sizehint, empty!, setindex!, get, haskey,
                              getkey, pop!, delete!, start, done, next,
                              isempty, length ]

similar{K,V,F}(d::DefaultDictBase{K,V,F}) = DefaultDictBase{K,V,F}()
in{T<:DefaultDictBase}(key, v::Base.KeyIterator{T}) = key in keys(v.dict.d)
next{T<:DefaultDictBase}(v::Base.KeyIterator{T}, i) = (v.dict.d.keys[i], Base.skip_deleted(v.dict.d,i+1))
next{T<:DefaultDictBase}(v::Base.ValueIterator{T}, i) = (v.dict.d.vals[i], Base.skip_deleted(v.dict.d,i+1))

getindex(d::DefaultDictBase, key) = get!(d.d, key, d.default)

# TODO: remove these if/when https://github.com/JuliaLang/julia/pull/5519 is committed
if !applicable(get!, (Dict,))
    global getindex
    function getindex{K,V,F<:Base.Callable}(d::DefaultDictBase{K,V,F,Dict}, key)
        if !haskey(d.d, key)
            return (d.d[key] = d.default())
        end
        return d.d[key]
    end

    function getindex{K,V,F}(d::DefaultDictBase{K,V,F,Dict}, key)
        if !haskey(d.d, key)
            return (d.d[key] = d.default)
        end
        return d.d[key]
    end
end

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

# Here begins the actual definition of the DefaultDict and
# DefaultOrderedDict classes.  As noted above, these are simply
# wrappers around a DefaultDictBase object, and delegate all functions
# to that object

for (DefaultDict,O) in [(:DefaultDict, :Unordered), (:DefaultOrderedDict, :Ordered)]
    @eval begin
        immutable $DefaultDict{K,V,F} <: Associative{K,V}
            d::DefaultDictBase{K,V,F,HashDict{K,V,$O}}

            $DefaultDict(x, kv::AbstractArray{(K,V)}) = new(DefaultDictBase{K,V,F,HashDict{K,V,$O}}(x, kv))
            $DefaultDict(x, d::$DefaultDict) = $DefaultDict(x, d.d)
            $DefaultDict(x, d::HashDict) = new(DefaultDictBase{K,V,F,HashDict{K,V,$O}}(x, d))
            $DefaultDict(x) = new(DefaultDictBase{K,V,F,HashDict{K,V,$O}}(x))
            $DefaultDict(x, ks, vs) = new(DefaultDictBase{K,V,F,HashDict{K,V,$O}}(x,ks,vs))
        end

        ## Constructors

        $DefaultDict() = error("$DefaultDict: no default specified")
        $DefaultDict(k,v) = error("$DefaultDict: no default specified")

        # TODO: these mimic similar Dict constructors, but may not be needed
        $DefaultDict{K,V,F}(default::F, ks::AbstractArray{K}, vs::AbstractArray{V}) = $DefaultDict{K,V,F}(default,ks,vs)
        $DefaultDict{F}(default::F,ks,vs) = $DefaultDict{Any,Any,F}(default, ks, vs)

        # syntax entry points
        $DefaultDict{F}(default::F) = $DefaultDict{Any,Any,F}(default)
        $DefaultDict{K,V,F}(::Type{K}, ::Type{V}, default::F) = $DefaultDict{K,V,F}(default)
        $DefaultDict{K,V,F}(default::F, kv::AbstractArray{(K,V)}) = $DefaultDict{K,V,F}(default, kv)

        $DefaultDict{F}(default::F, d::Associative) = ((K,V)=eltype(d); $DefaultDict{K,V,F}(default, HashDict(d)))

        ## Functions

        # Most functions are simply delegated to the wrapped DefaultDictBase object
        @delegate $DefaultDict.d [ sizehint, empty!, setindex!,
                                   getindex, get, get!, haskey,
                                   getkey, pop!, delete!, start, next,
                                   done, next, isempty, length]

        similar{K,V,F}(d::$DefaultDict{K,V,F}) = $DefaultDict{K,V,F}()
        in{T<:$DefaultDict}(key, v::Base.KeyIterator{T}) = key in keys(v.dict.d.d)
    end
end

## If/when a SortedDict becomes available, this should be uncommented to provide a DefaultSortedDict

# immutable DefaultSortedDict{K,V,F} <: Associative{K,V}
#     d::DefaultDictBase{K,V,F,SortedDict{K,V}}

#     DefaultSortedDict(x, kv::AbstractArray{(K,V)}) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x, kv))
#     DefaultSortedDict(x, d::DefaultSortedDict) = DefaultSortedDict(x, d.d)
#     DefaultSortedDict(x, d::SortedDict) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x, d))
#     DefaultSortedDict(x) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x))
#     DefaultSortedDict(x, ks, vs) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x,ks,vs))
# end

## Constructors

# DefaultSortedDict() = error("DefaultSortedDict: no default specified")
# DefaultSortedDict(k,v) = error("DefaultSortedDict: no default specified")

# # TODO: these mimic similar Dict constructors, but may not be needed
# DefaultSortedDict{K,V,F}(default::F, ks::AbstractArray{K}, vs::AbstractArray{V}) = DefaultSortedDict{K,V,F}(default,ks,vs)
# DefaultSortedDict{F}(default::F,ks,vs) = DefaultSortedDict{Any,Any,F}(default, ks, vs)

# # syntax entry points
# DefaultSortedDict{F}(default::F) = DefaultSortedDict{Any,Any,F}(default)
# DefaultSortedDict{K,V,F}(::Type{K}, ::Type{V}, default::F) = DefaultSortedDict{K,V,F}(default)
# DefaultSortedDict{K,V,F}(default::F, kv::AbstractArray{(K,V)}) = DefaultSortedDict{K,V,F}(default, kv)

# DefaultSortedDict{F}(default::F, d::Associative) = ((K,V)=eltype(d); DefaultSortedDict{K,V,F}(default, SortedDict(d)))

## Functions

## Most functions are simply delegated to the wrapped DefaultDictBase object

# @delegate DefaultSortedDict.d [ sizehint, empty!, setindex!,
#                            getindex, get, get!, haskey,
#                            getkey, pop!, delete!, start, next,
#                            done, next, isempty, length]

# similar{K,V,F}(d::DefaultSortedDict{K,V,F}) = DefaultSortedDict{K,V,F}()
# in{T<:DefaultSortedDict}(key, v::Base.KeyIterator{T}) = key in keys(v.dict.d.d)


text:
# Hill-Function
# https://www.wikiwand.com/en/Dose%E2%80%93response_relationship
# https://www.wikiwand.com/en/Hill_coefficient

# according to 10.1111/j.1600-0773.1997.tb00286.x:
hill_factory(n, ec50) = A -> 1 - A^n/(A^n + ec50^n)

# The Med and Expmed are used to specify medication usage in the simulation
struct Med
    time::Real        # time in minutes at wich the medication is added
    type::String      # type of medication (if type is equal to a key in the params dict that value will be changed to the value of dose)
    dose::Float64     # dose of medication (varies with type -> refere to handle_medication function)
end

mutable struct ExpMed # as above but with exponential decline
    time::Real
    type::String
    dose::Float64
    fade::Float64     # time * fade^(current_time - time) -> sets decline
end

# struct to encompasse a generic type of medication
# with potentially non linear active substance impact
@with_kw mutable struct Activa
    time :: Real
    duration :: Number = Inf
    dose :: Float64
    name :: String = ""
    fade :: Float64 = 1
    func :: Function = id
    param :: String
end

Meds = Union{Med,ExpMed,Activa}

gen_default_activa(time, dose; fade = 1.0, duration = Inf) = Activa(
    time = time,
    duration = duration,
    dose = dose,
    fade = fade,
    func = id,
    param = ""
)

function Dz(time, dose; fade = 1.0, duration = Inf)
    activa = gen_default_activa(time, dose, fade = fade, duration = duration)
    activa.name = "Diaxozide"
    activa.func = (x,_) -> params["gkatpbar"] + 25 * x
    activa.param = "gkatpbar"
    return activa
end

function Tolb(time, dose; fade = 1.0, duration = Inf, hill = true)
    activa = gen_default_activa(time, dose, fade = fade, duration = duration)
    activa.name = "Tolbutamid"
    if hill
        activa.func = (x,_) -> params["gkatpbar"] * (1 - x^1/(x^1 + 1.35^1))
    else
        activa.func = (x,_) -> params["gkatpbar"] - 25 * x
    end
    activa.param = "gkatpbar"
    return activa
end

function Tg(time, dose; fade = 1.0, duration = Inf)
    activa = gen_default_activa(time, dose, fade = fade, duration = duration)
    activa.name =  "Thapsigargin"
    activa.func = (x,_) -> params["kserca"] - 0.04x
    activa.param = "kserca"
    return activa
end

function KCl(time, dose; fade = 1.0, duration = Inf)
    activa = gen_default_activa(time, dose, fade = fade, duration = duration)
    activa.name = "KCl"
    activa.func = (x,_) -> vkf((x > 4.8) ? x : params["vk"], 130)
    activa.param = "vk"
    return activa
end

function Glucose(time, dose; fade = 1.0, duration = Inf)
    activa = gen_default_activa(time, dose, fade = fade, duration = duration)
    activa.name = "Glucose"
    activa.func = (x,_) -> x / 1000
    activa.param = "Jgk"
    return activa
end

# calculates the effects of medications
# by calling the func in the supplied activa
function handle_medication(med :: Activa, i :: Any)
    if isapprox(med.time + med.duration, i.t / 6000, rtol = 0.01) || med.dose ≤ 0
        i.p[med.param] = params[med.param]
    else
        i.p[med.param] = med.func(med.dose, i)
    end
end

# legacy medication handling
function handle_medication(med :: Union{Med, ExpMed}, i::Any)
    if med.type ≡ "tolbutamid"
        if med.dose > 0
                i.p["gkatpbar"] = params["gkatpbar"] - 0.02 * 1250 * med.dose
            else
                i.p["gkatpbar"] = params["gkatpbar"]
            end
    elseif med.type ≡ "Dz"
        if med.dose > 0
            i.p["gkatpbar"] = params["gkatpbar"] + 0.02 *  1250 * med.dose
        else
            i.p["gkatpbar"] = params["gkatpbar"]
        end
    elseif med.type ≡ "Tg"
        if med.dose > 0
                i.p["kserca"] = params["kserca"] - 0.04 * med.dose
            else
                i.p["kserca"] = params["kserca"]
            end
    elseif med.type ≡ "KCl"
            i.p["vk"] = vkf((med.dose > 4.8) ? med.dose : params["vk"], 130)
    elseif med.type ≡ "nifedipin"
        Bool(med.dose) ? i.p["gca"] = 500 : i.p["gca"] = 1000
    elseif med.type ≡ "glucose"
        i.p["Jgk"] = med.dose / 1000
    elseif med.type ⊆ keys(params)
        i.p[med.type] = med.dose
    end
end

# wrapper functions for handle_medication wich determain
function calc_dosage(m::Med, i)
    handle_medication(m, i)
    return true
end

# this version determains if the expoentialy falling
# dosage is still over the minimal threshold
function calc_dosage(m::ExpMed, i)
    m.dose *= m.fade^((i.t - m.time*6000)/6000)
    handle_medication(m, i)
    return m.dose ≤ 0.01 ? true : false
end

function calc_dosage(m::Activa, i)
    m.dose *= m.fade^((i.t - m.time*6000)/6000)
    handle_medication(m, i)
    return m.dose ≤ 0.01 ? true : false
end

function parse_medications(i)
    for m ∈ i.p["meds"]
        if m.time ≤ i.t / 6000
            if m isa Activa && m.time + m.duration ≥ i.t / 6000
                calc_dosage(m,i)
            elseif m isa Union{Med, ExpMed} && m.time ≤ i.t / 6000
                if calc_dosage(m, i);
                    filter!(n -> n ≢ m, i.p["meds"]) # remove "empty" Activa
                end
            end
        end
    end
end

ext: jl


text:
# Utils for modifying plot style variables

function PushStyleColor(idx, col::UInt32)
    ImPlot_PushStyleColor_U32(idx, col)
end

function PushStyleColor(idx, col::ImVec4)
    ImPlot_PushStyleColor_Vec4(idx, col)
end

function PushStyleVar(idx, val::AbstractFloat)
    ImPlot_PushStyleVarFloat(idx, Float32(val))
end

function PushStyleVar(idx, val::Integer)
    ImPlot_PushStyleVarInt(idx, Cint(val))
end

function PopStyleColor(count::Integer = 1)
    ImPlot_PopStyleColor(count)
end

function PopStyleVar(count::Integer = 1)
    ImPlot_PopStyleVar(count)
end


text:
module LassoTest

const  NUMBER_REPEAT = 1

using CoordinateDescent
using Test
using ProximalBase
using Random
using LinearAlgebra
using SparseArrays

Random.seed!(1)

##############################################
#
#  Lasso
#
##############################################

@testset "lasso" begin

  @testset "zero" begin
    n = 100
    p = 10

    X = randn(n, p)
    Y = X * ones(p) + 0.1 * randn(n)
    Xy = X' * Y / n

    lambda = maximum(abs.(Xy)) + 0.1
    out = lasso(X, Y, lambda)
    @test out.x == SparseIterate(p)
  end

  @testset "non-zero" begin
    for i=1:NUMBER_REPEAT
      n = 100
      p = 10
      s = 5

      X = randn(n, p)
      Y = X[:,1:s] * ones(s) + 0.1 * randn(n)

      λ = fill(0.3, p)
      beta = lasso(X, Y, 1., λ, CDOptions(;optTol=1e-12))

      f = CDQuadraticLoss(X'X/n, -X'Y/n)
      g = ProximalBase.ProxL1(1., λ)
      x1 = SparseIterate( p )
      coordinateDescent!(x1, f, g, CDOptions(;optTol=1e-12))
      @test beta.x ≈ x1 atol=1e-5

      @test (maximum(abs.(X'*(Y - X*beta.x) / n)) - 0.3) / 0.3 ≈ 0. atol=1e-5
    end
  end

  @testset "different interfaces" begin
    n = 500
    p = 500
    s = 50

    X = randn(n, p)
    β = randn(s)
    Y = X[:,1:s] * β + randn(n)

    λ = 0.1
    x1 = lasso(X, Y, λ)
    x2 = lasso(X, Y, λ, ones(p))

    @test Vector(x1.x) ≈ Vector(x2.x) atol=1e-5
  end

end

@testset "cd lasso" begin

  for i=1:NUMBER_REPEAT
    n = 200
    p = 50
    s = 10

    X = randn(n, p)
    β = randn(s)
    Y = X[:,1:s] * β + 0.1 * randn(n)

    g = ProximalBase.ProxL1(0.2)
    f1 = CDQuadraticLoss(X'X/n, -X'Y/n)
    f2 = CDLeastSquaresLoss(Y, X)

    x1 = SparseIterate(p)
    x2 = SparseIterate(p)
    coordinateDescent!(x1, f1, g, CDOptions(;optTol=1e-12))
    coordinateDescent!(x2, f2, g, CDOptions(;optTol=1e-12))

    @test maximum(abs.(x1 - x2)) ≈ 0. atol=1e-5
    @test (maximum(abs.(X'*(Y - X*x1) / n)) - 0.2) / 0.2 ≈ 0. atol=1e-5
    @test (maximum(abs.(X'*(Y - X*x2) / n)) - 0.2) / 0.2 ≈ 0. atol=1e-5
  end

end

@testset "cd sqrt-lasso" begin

  @testset "kkt" begin
    for i=1:NUMBER_REPEAT
      n = 100
      p = 50
      s = 5

      X = randn(n, p)
      β = randn(s)
      Y = X[:,1:s] * β + randn(n)

      λ = 2.8
      g = ProximalBase.ProxL1(λ)
      f = CDSqrtLassoLoss(Y, X)

      x1 = SparseIterate(p)
      coordinateDescent!(x1, f, g, CDOptions(;maxIter=5000, optTol=1e-8))

      @test max(0, maximum(abs.(X'*(Y - X*x1) / norm(Y - X*x1))) - λ) / λ  ≈ 0. atol=1e-3
    end
  end

  @testset "interfaces" begin
    for i=1:NUMBER_REPEAT
      n = 500
      p = 500
      s = 50

      X = randn(n, p)
      β = randn(s)
      Y = X[:,1:s] * β + randn(n)

      opt1 = CDOptions(;maxIter=5000, optTol=1e-10, warmStart=true, randomize=false)
      opt2 = CDOptions(;maxIter=5000, optTol=1e-10, warmStart=true, randomize=true)
      opt3 = CDOptions(;maxIter=5000, optTol=1e-10, warmStart=false, randomize=false)
      opt4 = CDOptions(;maxIter=5000, optTol=1e-10, warmStart=false, randomize=true)

      x1 = SparseIterate(sprand(p, 0.6))
      x2 = SparseIterate(sprand(p, 0.6))
      x3 = SparseIterate(sprand(p, 0.6))
      x4 = SparseIterate(sprand(p, 0.6))

      λ = 1.5
      g = ProximalBase.ProxL1(λ)
      f = CDSqrtLassoLoss(Y, X)

      coordinateDescent!(x1, f, g, opt1)
      coordinateDescent!(x2, f, g, opt2)
      coordinateDescent!(x3, f, g, opt3)
      coordinateDescent!(x4, f, g, opt4)

      @test Vector(x1) ≈ Vector(x2) atol=1e-4
      @test Vector(x3) ≈ Vector(x2) atol=1e-4
      @test Vector(x4) ≈ Vector(x2) atol=1e-4

      y1 = sqrtLasso(X, Y, λ, opt1, standardizeX=false)
      y2 = sqrtLasso(X, Y, λ, opt2, standardizeX=false)
      y3 = sqrtLasso(X, Y, λ, opt3, standardizeX=false)
      y4 = sqrtLasso(X, Y, λ, opt4, standardizeX=false)

      @test Vector(y1.x) ≈ Vector(x2) atol=1e-4
      @test Vector(y2.x) ≈ Vector(x2) atol=1e-4
      @test Vector(y3.x) ≈ Vector(x2) atol=1e-4
      @test Vector(y4.x) ≈ Vector(x2) atol=1e-4

      z1 = sqrtLasso(X, Y, λ, ones(p), opt1)
      z2 = sqrtLasso(X, Y, λ, ones(p), opt2)
      z3 = sqrtLasso(X, Y, λ, ones(p), opt3)
      z4 = sqrtLasso(X, Y, λ, ones(p), opt4)

      @test Vector(z1.x) ≈ Vector(x2) atol=1e-4
      @test Vector(z2.x) ≈ Vector(x2) atol=1e-4
      @test Vector(z3.x) ≈ Vector(x2) atol=1e-4
      @test Vector(z4.x) ≈ Vector(x2) atol=1e-4

    end
  end

end

@testset "scaled lasso" begin

  for i=1:NUMBER_REPEAT
    n = 1000
    p = 500
    s = 50

    X = randn(n, p)
    β = randn(s)
    Y = X[:,1:s] * β + randn(n)

    λ = rand() / 5.

    opt1 = IterLassoOptions(;maxIter=100, optTol=1e-8, optionsCD=CDOptions(;maxIter=5000, optTol=1e-8))
    opt2 = IterLassoOptions(;maxIter=100, optTol=1e-8, initProcedure=:InitStd, σinit=2., optionsCD=CDOptions(;maxIter=5000, optTol=1e-8))

    x1 = SparseIterate(p)
    x2 = SparseIterate(p)
    sol1 = scaledLasso!(x1, X, Y, λ, ones(p), opt1)
    sol2 = scaledLasso!(x2, X, Y, λ, ones(p), opt2)

    σ1 = sol1.σ
    σ2 = sol2.σ

    @test max.((maximum(abs.(X'*(Y - X*x1) / n)) - λ*σ1), 0.) / (σ1*λ) ≈ 0. atol=1e-4
    @test max.((maximum(abs.(X'*(Y - X*x2) / n)) - λ*σ2), 0.) / (σ2*λ) ≈ 0. atol=1e-4
    @test Vector(x1) ≈ Vector(x2) atol=1e-4
  end

end

@testset "lasso path" begin

  @testset "standardizeX = false" begin
    n = 1000
    p = 500
    s = 50

    X = randn(n, p)
    β = randn(s)
    Y = X[:,1:s] * β + randn(n)

    λ1 = 0.3
    λ2 = 0.1
    opt = CDOptions(;maxIter=5000, optTol=1e-8)

    x1 = lasso(X, Y, λ1)
    x2 = lasso(X, Y, λ2)

    λpath = [λ1, λ2]
    path = LassoPath(X, Y, λpath, opt; standardizeX=false)

    @test typeof(path) == LassoPath{Float64}
    @test Vector(path.βpath[1]) ≈ Vector(x1.x) atol=1e-5
    @test Vector(path.βpath[2]) ≈ Vector(x2.x) atol=1e-5

    S1 = findall(!iszero, x1.x)
    S2 = findall(!iszero, x2.x)
    rf = refitLassoPath(path, X, Y)

    @test rf[S1] ≈ X[:,S1] \ Y atol=1e-5
    @test rf[S2] ≈ X[:,S2] \ Y atol=1e-5
  end

  @testset "standardizeX = true" begin
    n = 1000
    p = 500
    s = 50

    X = randn(n, p)
    β = randn(s)
    Y = X[:,1:s] * β + randn(n)

    loadingX = Array{Float64}(undef, p)
    CoordinateDescent._stdX!(loadingX, X)

    λ1 = 0.3
    λ2 = 0.1
    opt = CDOptions(;maxIter=5000, optTol=1e-8)

    x1 = lasso(X, Y, λ1, loadingX, opt)
    x2 = lasso(X, Y, λ2, loadingX, opt)

    λpath = [λ1, λ2]
    path = LassoPath(X, Y, λpath, opt)

    @test typeof(path) == LassoPath{Float64}
    @test Vector(path.βpath[1]) ≈ Vector(x1.x) atol=1e-5
    @test Vector(path.βpath[2]) ≈ Vector(x2.x) atol=1e-5

    S1 = findall(!iszero, x1.x)
    S2 = findall(!iszero, x2.x)
    rf = refitLassoPath(path, X, Y)

    @test rf[S1] ≈ X[:,S1] \ Y atol=1e-5
    @test rf[S2] ≈ X[:,S2] \ Y atol=1e-5
  end

end

end


text:
module TestDef

using Test
using InitialValues: InitialValue, hasinitialvalue, isknown

module CleanNameSpace
    using InitialValues: @def, @disambiguate
    add(x, y) = x + y
    got(x) = string("Got: ", repr(x))
    @def add got(x)
    @disambiguate add Missing
end

@testset "CleanNameSpace" begin
    add = CleanNameSpace.add
    @test !isdefined(CleanNameSpace, :InitialValues)
    @test hasinitialvalue(add)
    @test isknown(InitialValue(add))
    @test add(InitialValue(add), :x) == "Got: :x"
    @test add(InitialValue(add), missing) === missing
end

module NonFunction
    using InitialValues: @def, @disambiguate
    struct Add end
    const add = Add()
    add(x, y) = x + y
    @def add
    @disambiguate add Missing
end

@testset "NonFunction" begin
    add = NonFunction.add
    @test !(add isa Function)
    @test hasinitialvalue(add)
    @test isknown(InitialValue(add))
    @test add(InitialValue(add), :x) == :x
    @test add(InitialValue(add), missing) === missing
end

end  # module


text:
function Base.show(io::IO, ::MIME"text/html", x::T) where {T <: Polynomials.AbstractPolynomial}
#     #write(io, "
")
     write(io, "
")
     show(io, "text/latex", x)
     write(io, "
")
 end


text:
function mp_contignecy(data_ip, Total_sample, total_yr, br_no, base_weight_ip, weights)
    Cont_list = Any[]
    base_list = Any[] #(base_nw, cont_nw, cont_br)
    year_base = Any[]
    base_weight = base_weight_ip
    st_nw = 0
    end_nw = 0
    for i = 1:Total_sample*total_yr
         if i == 1
            st_nw = 2
            end_nw = 2
         else
             st_nw = end_nw + 2
             end_nw = end_nw + 2
        end
        push!(base_list, st_nw - 1)
        base_weight[st_nw - 1] = weights[i]
        for n = st_nw:end_nw
            push!(Cont_list, (st_nw - 1, n, br_no))
            data_ip["nw"]["$n"]["branchdc"]["$br_no"]["rateA"] = 0
        end
    end

    return data_ip, Cont_list, base_list, base_weight
end

function mp_contignecy_nocl(data_ip, Total_sample, br_no)
    Cont_list = Any[]
    base_list = Any[] #(base_nw, cont_nw, cont_br)
    st_nw = 0
    end_nw = 0
    for i = 1:Total_sample
         if i == 1
            st_nw = 2
            end_nw = 2
         else
             st_nw = end_nw + 2
            end_nw = end_nw + 2
        end
        push!(base_list, st_nw - 1)
        for n = st_nw:end_nw
            push!(Cont_list, (st_nw - 1, n, br_no))
            data_ip["nw"]["$n"]["branchdc"]["$br_no"]["rateA"] = 0
         end
    end
    return data_ip, Cont_list, base_list
end

function year_base_networks(Total_sample, total_yr, base_list)
    base_list_no = 1
    year_base =  Dict([(yr, []) for yr = 1: total_yr])
    for k = 1:total_yr
            for t = 1:Total_sample
                     push!(year_base[k], base_list[base_list_no])
                    base_list_no += 1
            end
    end
    return year_base
end

function mp_datainputs(data_ip,Total_sample, year, year_num,file)
    kk= 1
    for yr = 1:length(year)
        fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output.mat")
        vars = matread(fname)
        fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output_MC.mat")
        vars1 = matread(fname1)
        fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
        vars2 = matread(fname2)

        file_cl = string("cluster_", year_num[yr], ".jld2")
        cl = load(file_cl)
    for i = 1:Total_sample
        # sample = Int(round(8760*(1-rand())))
          sample = cl["scenarios"][i]
         for tt = 1:2
            data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = vars["Heq"][sample]
            if data_ip["nw"]["$kk"]["reserves"]["2"]["H"] < 0.5
                data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = 0.5
            end
            # data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2] = maximum(vars1["MC1"])*100*scenarios[1,i] for kmeans
            data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2] = vars1["MC1"][sample]*100
            rating = deepcopy(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
            data_ip["nw"]["$kk"]["gen"]["1"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
            data_ip["nw"]["$kk"]["gen"]["2"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
            if occursin("6bus", file)
                data_ip["nw"]["$kk"]["gen"]["4"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
            end
            data_ip["nw"]["$kk"]["gen"]["3"]["pmax"] = vars["Ptotal"][sample]/100
            data_ip["nw"]["$kk"]["load"]["1"]["pd"] = vars["Ptotal"][sample]/100
            # data_ip["nw"]["$kk"]["gen"]["3"]["pmax"] = maximum(vars["Ptotal"])*scenarios[2,i]/100 # for kmeans
            # data_ip["nw"]["$kk"]["load"]["1"]["pd"]  = maximum(vars["Ptotal"])*scenarios[2,i]/100 # for k means
            # what about load?
            kk +=1
        end
    end
    end
    return data_ip
end

    function mp_datainputs_nocl(data_ip, Total_sample, year, year_num, scenario, file)
        kk= 1
        for yr = 1:length(year)
            fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output.mat")
            vars = matread(fname)
            fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output_MC.mat")
            vars1 = matread(fname1)
            fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
            vars2 = matread(fname2)

        for i = 1:Total_sample
            # sample = Int(round(8760*(1-rand())))
            sample = scenario[i]
            for tt = 1:2
                data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = vars["Heq"][sample]
                if data_ip["nw"]["$kk"]["reserves"]["2"]["H"] < 0.5
                    data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = 0.5
                end
                data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2] = vars1["MC1"][sample]*100
                rating = deepcopy(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
                data_ip["nw"]["$kk"]["gen"]["1"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                data_ip["nw"]["$kk"]["gen"]["2"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                if occursin("6bus", file)
                    data_ip["nw"]["$kk"]["gen"]["4"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                end
                data_ip["nw"]["$kk"]["gen"]["3"]["pmax"] = vars["Ptotal"][sample]/100
                data_ip["nw"]["$kk"]["load"]["1"]["pd"] = vars["Ptotal"][sample]/100
                # what about load?
                kk +=1
            end
        end
        end
        return data_ip
    end

    function mp_datainputs_nocl_PLdim(data_ip, Total_sample, year, year_num, scenario, file)
        kk= 1
        for yr = 1:length(year)
            fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output.mat")
            vars = matread(fname)
            fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output_MC.mat")
            vars1 = matread(fname1)
            fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
            vars2 = matread(fname2)

        for i = 1:Total_sample
            # sample = Int(round(8760*(1-rand())))
            sample = scenario[i]
            for tt = 1:2
                data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = vars["Heq"][sample]
                if data_ip["nw"]["$kk"]["reserves"]["2"]["H"] < 0.5
                    data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = 0.5
                end
                data_ip["nw"]["$kk"]["gen"]["2"]["cost"][2] = vars1["MC1"][sample]*100
                # rating = deepcopy(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
                # data_ip["nw"]["$kk"]["gen"]["1"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                if occursin("6bus", file)
                    data_ip["nw"]["$kk"]["gen"]["4"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                end
                data_ip["nw"]["$kk"]["gen"]["2"]["pmax"] = vars["Ptotal"][sample]/100
                data_ip["nw"]["$kk"]["load"]["1"]["pd"] = vars["Ptotal"][sample]/100
                # what about load?
                kk +=1
            end
        end
        end
        return data_ip
    end

    function mp_contignecy_ne(data_ip, Total_sample, base_weight_ip, weights)
        Cont_list = Any[]
        base_list = Any[] #(base_nw, cont_nw, cont_br)
        base_weight = base_weight_ip
        cont_no = 1
        for i = 1:Total_sample
             if i == 1
                global st_nw = 2
                global end_nw = length(data_ip["nw"]["1"]["branchdc_ne"]) + 1
             else
                st_nw = end_nw + 2
                end_nw = (st_nw - 1) + length(data_ip["nw"]["1"]["branchdc_ne"])
            end
            push!(base_list, st_nw - 1)
            base_weight[st_nw - 1] = weights[i]
            br_no = 1
            for n = st_nw:end_nw
                push!(Cont_list, (st_nw - 1, n, br_no))
                data_ip["nw"]["$n"]["branchdc_ne"]["$br_no"]["rateA"] = 0
                br_no += 1
                cont_no += 1
            end
        end
        return data_ip, Cont_list, base_list, base_weight
    end

    function mp_contignecy_ne_PL(data_ip, Total_sample, br_no, base_weight_ip, weights)
        Cont_list = Any[]
        base_list = Any[] #(base_nw, cont_nw, cont_br)
        base_weight = base_weight_ip
        cont_no = 1
        for i = 1:Total_sample
             if i == 1
                global st_nw = 2
                global end_nw = 2
             else
                st_nw = end_nw + 2
                end_nw = end_nw + 2
            end
            push!(base_list, st_nw - 1)
            base_weight[st_nw - 1] = weights[i]
            for n = st_nw:end_nw
                push!(Cont_list, (st_nw - 1, n, br_no))
                display(n)
                data_ip["nw"]["$n"]["branchdc_ne"]["$br_no"]["rateA"] = 0
                cont_no += 1
            end
        end
        return data_ip, Cont_list, base_list, base_weight
    end

    function mp_datainputs_ne(data_ip,Total_sample, year, year_num,file)
        kk= 1
        cost = Any[]
        for yr = 1:length(year)
            fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output.mat")
            vars = matread(fname)
            fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output_MC.mat")
            vars1 = matread(fname1)
            fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
            vars2 = matread(fname2)
            file_cl = string("cluster_", year_num[yr], ".jld2")
            cl = load(file_cl)
        for i = 1:Total_sample
            # sample = Int(round(8760*(1-rand())))
            sample = cl["scenarios"][i]
             for tt = 1:(length(data_ip["nw"]["1"]["branchdc_ne"])+1)
               data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = vars["Heq"][sample]
                if data_ip["nw"]["$kk"]["reserves"]["2"]["H"] < 0.5
                    data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = 0.5
                end
                data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2] = vars1["MC1"][sample]*100
                rating = deepcopy(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
                data_ip["nw"]["$kk"]["gen"]["1"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                data_ip["nw"]["$kk"]["gen"]["2"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                if occursin("6bus", file)
                    data_ip["nw"]["$kk"]["gen"]["4"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                end
                data_ip["nw"]["$kk"]["gen"]["3"]["pmax"] = vars["Ptotal"][sample]/100
                data_ip["nw"]["$kk"]["load"]["1"]["pd"] = vars["Ptotal"][sample]/100
                push!(cost,  data_ip["nw"]["$kk"]["gen"]["3"]["pmax"]*data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2])
                 kk +=1
            end
        end
        end
        return data_ip, cost
    end

    function mp_datainputs_ne_PL(data_ip,Total_sample, year, year_num,file)
        kk= 1
        cost = Any[]
        for yr = 1:length(year)
            fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output.mat")
            vars = matread(fname)
            fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output_MC.mat")
            vars1 = matread(fname1)
            fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
            vars2 = matread(fname2)
            file_cl = string("cluster_", year_num[yr], ".jld2")
            cl = load(file_cl)
        for i = 1:Total_sample
            # sample = Int(round(8760*(1-rand())))
            sample = cl["scenarios"][i]
            for tt = 1:2
                data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = vars["Heq"][sample]
                if data_ip["nw"]["$kk"]["reserves"]["2"]["H"] < 0.5
                    data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = 0.5
                end
                data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2] = vars1["MC1"][sample]*100
                rating = deepcopy(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
                data_ip["nw"]["$kk"]["gen"]["1"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                data_ip["nw"]["$kk"]["gen"]["2"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                if occursin("6bus", file)
                    data_ip["nw"]["$kk"]["gen"]["4"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                end
                data_ip["nw"]["$kk"]["gen"]["3"]["pmax"] = vars["Ptotal"][sample]/100
                data_ip["nw"]["$kk"]["load"]["1"]["pd"] = vars["Ptotal"][sample]/100
                # what about load?
                # display("Pgmax")
                display(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
                push!(cost,  data_ip["nw"]["$kk"]["gen"]["1"]["pmax"]* data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2])
                 kk +=1
            end
        end
        end
        return data_ip, cost
    end

    function mp_datainputs_ne_r1(data_ip,Total_sample)
        kk= 1
        cost = Any[]
        for yr = 1:length(year)
            fname = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output.mat")
            vars = matread(fname)
            fname1 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\",year[yr],"\\output_MC.mat")
            vars1 = matread(fname1)
            fname2 = string("C:\\Users\\djaykuma\\OneDrive - Energyville\\Freq_TNEP_paper\\MATLAB\\NAT_Results\\wind_sample.mat")
            vars2 = matread(fname2)
            file_cl = string("cluster_", year_num[yr], ".jld2")
            cl = load(file_cl)
        for i = 1:Total_sample
            # sample = Int(round(8760*(1-rand())))
            sample = cl["scenarios"][i]
            # display("sample:$sample")
            for tt = 1:(length(data_sp["branchdc_ne"])+1)
               data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = vars["Heq"][sample]
                if data_ip["nw"]["$kk"]["reserves"]["2"]["H"] < 0.5
                    data_ip["nw"]["$kk"]["reserves"]["2"]["H"] = 0.5
                end
                data_ip["nw"]["$kk"]["gen"]["1"]["cost"][2] = vars1["MC1"][sample]*100
                data_ip["nw"]["$kk"]["gen"]["2"]["cost"][2] = vars1["MC1"][sample]*100
                rating = deepcopy(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
                data_ip["nw"]["$kk"]["gen"]["1"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                data_ip["nw"]["$kk"]["gen"]["2"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                if occursin("6bus", file)
                    data_ip["nw"]["$kk"]["gen"]["4"]["pmax"] = vars2["winddata"][sample]/maximum(vars2["winddata"])*rating
                end
                data_ip["nw"]["$kk"]["gen"]["3"]["pmax"] = vars["Ptotal"][sample]/100
                data_ip["nw"]["$kk"]["load"]["1"]["pd"] = vars["Ptotal"][sample]/100
                # what about load?
                # display(data_ip["nw"]["$kk"]["gen"]["1"]["pmax"])
                push!(cost,  data_ip["nw"]["$kk"]["gen"]["3"]["pmax"]*data_ip["nw"]["$kk"]["gen"]["3"]["cost"][2])
                 kk +=1
            end
        end
        end
        return data_ip, cost
    end

    function mp_contignecy_nocl_sensitivity(data_ip, Total_sample, br_no)
        Cont_list = Any[]
        base_list = Any[] #(base_nw, cont_nw, cont_br)
        st_nw = 0
        end_nw = 0
        for i = 1:Total_sample
             if i == 1
                st_nw = 2
                end_nw = 2
             else
                 st_nw = end_nw + 2
                end_nw = end_nw + 2
            end
            push!(base_list, st_nw - 1)
            for n = st_nw:end_nw
                push!(Cont_list, (st_nw - 1, n, br_no))
                data_ip["nw"]["$n"]["branchdc"]["$br_no"]["rateA"] = 0
             end
        end
        return data_ip, Cont_list, base_list
    end


text:
using LinearAlgebra
using SparseArrays

export expmv

"""
    expmv!(f, t, A, b; <keyword arguments>)

Returns the matrix-vector product ``exp(t A) b`` where `A` is a ``n × n`` sparse
real or complex matrix or matrix-free object, `b` is a vector of ``n`` real or
complex elements, and `t` is a time step parameter (or a `StepRangeLen` object
representing a range of values).

# Arguments
* `f`: an `n`-vector for preallocating the output
* `t`: `Number` or `StepRangeLen` object
* `A`: a `n × n` real or complex sparse matrix
* `b`: an `n`-vector
* `M = []`: manually set the degree of the Taylor expansion
* `precision = :double`: can be `:double`, `:single` or `:half`.
* `shift = $DEFAULT_SHIFT`: set to `true` to apply a shift in order to reduce the norm of A
        (see Sec. 3.1 of the paper)
* `full_term = false`: set to `true` to evaluate the full Taylor expansion instead
        of truncating when reaching the required precision
* `check_positive = false`: set to `true` to check if `A` has strictly positive entries
        for a faster evaluation of the 1-norm
"""
function expmv!(
        f::AbstractVecOrMat, t::Number, A, b::AbstractVecOrMat,
        M = nothing, norm = LinearAlgebra.norm, opnorm = LinearAlgebra.opnorm,
        b1 = similar(b), b2 = similar(b);
        precision = :double, shift = DEFAULT_SHIFT, full_term = false, check_positive = false
    )

    if shift == true && !hasmethod(tr, Tuple{typeof(A)})
        shift = false
        @warn "Shift set to false as $(typeof(A)) doesn't support tr"
    end

    n = size(A, 2)
    Tr = real(eltype(b))

    mu = zero(Tr)
    if shift
        mu = tr(A)/n
        A = diagshift!(A, mu) #A = A - mu * I
    end

    if M == nothing
        tt = 1
        (M,alpha,unA) = select_taylor_degree(t*A, size(b,2); precision=precision, shift=shift, check_positive=check_positive)
    else
        tt = t
    end

    tol =
      if precision == :double
          Tr(2.0^(-53))
      elseif precision == :single
          Tr(2.0^(-24))
      elseif precision == :half
          Tr(2.0^(-10))
      end

    s = 1

    if t == 0
        m = 0
    else
        (m_max,p) = size(M)
        U = diagm(0 => 1:m_max)
        C = ((ceil.(abs.(tt)*M))'*U )

        C[C .== 0] .= Inf

        cost, idx = findmin(C')
        # idx is a CarthesianIndex if C' is a Matrix, or a scalar if C' is a row
        # vector. idx[1] extract the first index, i.e. row, of the CarthesianIndex
        m = idx[1]

        if cost == Inf
            cost = 0
        end

        s = max(cost/m,1)
    end

    eta = 1

    if shift
        eta = exp(t*mu/s)
    end

    f = copyto!(f, b)
    b1 = copyto!(b1, b)
    b2 = copyto!(b2, b)
    c1 = c2 = Tr(Inf)

    @inbounds for i = 1:s
        if !full_term
            c1 = norm(b1, Inf) # only need to update if !full_term
        end

        for k = 1:m
            b2 = mul!(b2, A, b1)
            b2 .*= (t / (s * k)) # b = (t/(s*k))*(A*b)
            b1, b2 = b2, b1
            f .+= b1 # f = axpy!(1, b1, f)

            if !full_term
                c2 = norm(b1, Inf) # only need to update if !full_term
                if c1 + c2 <= tol * norm(f, Inf)
                    break
                end
                c1 = c2
            end

        end

        if shift
            f .*= eta
        end

        b1 = copyto!(b1, f)
    end

    return f
end

function expmv(
        t::Number, A, b::AbstractVecOrMat,
        M = nothing, norm = LinearAlgebra.norm, opnorm = LinearAlgebra.opnorm,
        b1 = copy(b), b2 = similar(b);
        precision = :double, shift = DEFAULT_SHIFT, full_term = false, check_positive = false
    )

    return expmv!(
        similar(b), t, A, b, M, norm, opnorm, b1, b2;
        precision = precision, shift = shift, full_term = full_term, check_positive = check_positive
    )
end


text:
# using POMDPs

# import Pkg
# Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaPOMDP/DiscreteValueIteration.jl"))
# POMDPs.add("DiscreteValueIteration", false)


text:
using BrainFlow

# specify logging library to use
BrainFlow.enable_dev_logger(BrainFlow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = BrainFlow.BoardShim(BrainFlow.SYNTHETIC_BOARD, params)
sampling_rate = BrainFlow.get_sampling_rate(BrainFlow.SYNTHETIC_BOARD)
nfft = BrainFlow.get_nearest_power_of_two(sampling_rate)

BrainFlow.prepare_session(board_shim)
BrainFlow.start_stream(board_shim)
sleep(5)
BrainFlow.stop_stream(board_shim)
data = BrainFlow.get_board_data(board_shim)
BrainFlow.release_session(board_shim)

eeg_channels = BrainFlow.get_eeg_channels(BrainFlow.SYNTHETIC_BOARD)

bands = BrainFlow.get_avg_band_powers(data, eeg_channels, sampling_rate, true)
println(bands[1])
println(bands[2])


text:
using DataStructures
using Distributions

include("tree_utils.jl")
include("../dynamics_model/dubins_model.jl")

function model_check!(tree::KDTREE; max_iterations = 200, belres = 1e-6, γ = 1.0)
    for i = 1:max_iterations
        residual = 0.0

        # Stacks to go through tree
        lb_s = Stack{Vector{Float64}}()
        ub_s = Stack{Vector{Float64}}()
        s = Stack{Union{LEAFNODE, KDNODE}}()

        push!(lb_s, tree.lbs)
        push!(ub_s, tree.ubs)
        push!(s, tree.root_node)

        while !isempty(s)
            curr = pop!(s)
            curr_lbs = pop!(lb_s)
            curr_ubs = pop!(ub_s)

            if typeof(curr) == LEAFNODE
                min_cte = curr_lbs[1]
                max_cte = curr_ubs[1]

                if (min_cte < -10.0) || (max_cte > 10.0)
                    curr.prob = 1.0
                else
                    old_prob = curr.prob

                    next_lbs, next_ubs = reachable_cell(curr_lbs, curr_ubs, curr.min_control, curr.max_control)
                    next_nodes = get_overlapping_nodes(tree.root_node, next_lbs, next_ubs)
                    new_prob = γ * maximum([node.prob for node in next_nodes])
                    curr.prob = new_prob

                    change = abs(old_prob - new_prob)
                    change > residual && (residual = change)
                end
            else
                # Traverse tree and keep track of bounds
                dim = curr.dim
                split = curr.split
                # Go left, upper bounds will change
                left_ubs = copy(curr_ubs)
                left_ubs[dim] = split

                push!(lb_s, curr_lbs)
                push!(ub_s, left_ubs)
                push!(s, curr.left)

                # Go right, lower bounds will change
                right_lbs = copy(curr_lbs)
                right_lbs[dim] = split

                push!(lb_s, right_lbs)
                push!(ub_s, curr_ubs)
                push!(s, curr.right)
            end
        end
        println("[Iteration $i] residual: $residual")
        residual < belres && break
    end
end

function forward_reach(init_tree::KDTREE; max_iter = 50, verbose = false)
    # This function assumes that the first tree is labeled already with the states you
    # want to allow the system to start in

    trees = [init_tree]

    for i = 1:max_iter
        tree = trees[end]
        next_tree = copy(init_tree)
        zero_out!(next_tree)

        # Stacks to go through tree
        lb_s = Stack{Vector{Float64}}()
        ub_s = Stack{Vector{Float64}}()
        s = Stack{Union{LEAFNODE, KDNODE}}()

        push!(lb_s, tree.lbs)
        push!(ub_s, tree.ubs)
        push!(s, tree.root_node)

        while !isempty(s)
            curr = pop!(s)
            curr_lbs = pop!(lb_s)
            curr_ubs = pop!(ub_s)

            if typeof(curr) == LEAFNODE
                if curr.prob == 1
                    next_lbs, next_ubs = reachable_cell(curr_lbs, curr_ubs, curr.min_control, curr.max_control)
                    #curr_nodes = get_overlapping_nodes(tree.root_node, next_lbs, next_ubs)
                    next_nodes = get_overlapping_nodes(next_tree.root_node, next_lbs, next_ubs)
                    for node in next_nodes
                        node.prob = curr.prob
                    end
                end
             else
                # Traverse tree and keep track of bounds
                dim = curr.dim
                split = curr.split
                # Go left, upper bounds will change
                left_ubs = copy(curr_ubs)
                left_ubs[dim] = split

                push!(lb_s, curr_lbs)
                push!(ub_s, left_ubs)
                push!(s, curr.left)

                # Go right, lower bounds will change
                right_lbs = copy(curr_lbs)
                right_lbs[dim] = split

                push!(lb_s, right_lbs)
                push!(ub_s, curr_ubs)
                push!(s, curr.right)
            end
        end
        push!(trees, next_tree)
        # Check convergence
        converged = tree == next_tree
        verbose && println("[Iteration $i] converged: $converged")
        converged && break
    end
    return trees
end

function forward_reach!(trees::Vector{KDTREE}; verbose = false)
    # This function assumes that the first tree is labeled already with the states you
    # want to allow the system to start in

    for i = 1:length(trees) - 1
        tree = trees[i]
        next_tree = trees[i + 1]

        # Stacks to go through tree
        lb_s = Stack{Vector{Float64}}()
        ub_s = Stack{Vector{Float64}}()
        s = Stack{Union{LEAFNODE, KDNODE}}()

        push!(lb_s, tree.lbs)
        push!(ub_s, tree.ubs)
        push!(s, tree.root_node)

        while !isempty(s)
            curr = pop!(s)
            curr_lbs = pop!(lb_s)
            curr_ubs = pop!(ub_s)

            if typeof(curr) == LEAFNODE
                next_lbs, next_ubs = reachable_cell(curr_lbs, curr_ubs, curr.min_control, curr.max_control)
                curr_nodes = get_overlapping_nodes(tree.root_node, next_lbs, next_ubs)
                next_nodes = get_overlapping_nodes(next_tree.root_node, next_lbs, next_ubs)
                for node in next_nodes
                    node.prob = curr.prob
                end
             else
                # Traverse tree and keep track of bounds
                dim = curr.dim
                split = curr.split
                # Go left, upper bounds will change
                left_ubs = copy(curr_ubs)
                left_ubs[dim] = split

                push!(lb_s, curr_lbs)
                push!(ub_s, left_ubs)
                push!(s, curr.left)

                # Go right, lower bounds will change
                right_lbs = copy(curr_lbs)
                right_lbs[dim] = split

                push!(lb_s, right_lbs)
                push!(ub_s, curr_ubs)
                push!(s, curr.right)
            end
        end
        # Check convergence
        converged = tree == next_tree
        verbose && println("[Iteration $i] converged: $converged")
    end
end

function forward_reach(total_tree::KDTREE, init_xy_tree::KDTREE; max_iter = 50, verbose = false)
    # This function assumes that the first tree is labeled already with the states you
    # want to allow the system to start in

    total_tree_next = copy(total_tree)
    zero_out!(total_tree_next)

    label_xy!(init_xy_tree, total_tree)
    trees = [init_xy_tree]

    for i = 1:max_iter
        # Stacks to go through tree
        lb_s = Stack{Vector{Float64}}()
        ub_s = Stack{Vector{Float64}}()
        s = Stack{Union{LEAFNODE, KDNODE}}()

        push!(lb_s, total_tree.lbs)
        push!(ub_s, total_tree.ubs)
        push!(s, total_tree.root_node)

        while !isempty(s)
            curr = pop!(s)
            curr_lbs = pop!(lb_s)
            curr_ubs = pop!(ub_s)

            if typeof(curr) == LEAFNODE
                if curr.prob == 1.0 # Only need to update from cells that are currently reachable
                    next_lbs, next_ubs = reachable_cell_dtp(curr_lbs, curr_ubs, curr.min_control, curr.max_control)
                    next_nodes = get_overlapping_nodes(total_tree_next.root_node, next_lbs, next_ubs)
                    for node in next_nodes
                        node.prob = curr.prob
                    end
                end
             else
                # Traverse tree and keep track of bounds
                dim = curr.dim
                split = curr.split
                # Go left, upper bounds will change
                left_ubs = copy(curr_ubs)
                left_ubs[dim] = split

                push!(lb_s, curr_lbs)
                push!(ub_s, left_ubs)
                push!(s, curr.left)

                # Go right, lower bounds will change
                right_lbs = copy(curr_lbs)
                right_lbs[dim] = split

                push!(lb_s, right_lbs)
                push!(ub_s, curr_ubs)
                push!(s, curr.right)
            end
        end
        curr_xy_tree = trees[end]
        next_xy_tree = copy(curr_xy_tree)
        label_xy!(next_xy_tree, total_tree_next)
        push!(trees, next_xy_tree)

        copy_labels!(total_tree, total_tree_next)
        zero_out!(total_tree_next)

        verbose && (println("[Iteration:] ", i))
    end
    return trees
end

# NOTE: I actually do not think this works
function forward_reach_runway(total_tree::KDTREE; max_iter = 50, verbose = false)
    # This function assumes that the first tree is labeled already with the states you
    # want to allow the system to start in

    total_tree_next = copy(total_tree)
    zero_out!(total_tree_next)

    init_lbs, init_ubs = get_xy_bounds(total_tree)
    lbs = [init_lbs]
    ubs = [init_ubs]

    for i = 1:max_iter
        # Stacks to go through tree
        lb_s = Stack{Vector{Float64}}()
        ub_s = Stack{Vector{Float64}}()
        s = Stack{Union{LEAFNODE, KDNODE}}()

        push!(lb_s, total_tree.lbs)
        push!(ub_s, total_tree.ubs)
        push!(s, total_tree.root_node)

        while !isempty(s)
            curr = pop!(s)
            curr_lbs = pop!(lb_s)
            curr_ubs = pop!(ub_s)

            if typeof(curr) == LEAFNODE
                if curr.prob == 1.0 # Only need to update from cells that are currently reachable
                    next_lbs, next_ubs = reachable_cell_dtp(curr_lbs, curr_ubs, curr.min_control, curr.max_control)
                    next_nodes = get_overlapping_nodes(total_tree_next.root_node, next_lbs, next_ubs)
                    for node in next_nodes
                        node.prob = curr.prob
                    end
                end
             else
                # Traverse tree and keep track of bounds
                dim = curr.dim
                split = curr.split
                # Go left, upper bounds will change
                left_ubs = copy(curr_ubs)
                left_ubs[dim] = split

                push!(lb_s, curr_lbs)
                push!(ub_s, left_ubs)
                push!(s, curr.left)

                # Go right, lower bounds will change
                right_lbs = copy(curr_lbs)
                right_lbs[dim] = split

                push!(lb_s, right_lbs)
                push!(ub_s, curr_ubs)
                push!(s, curr.right)
            end
        end
        xy_lbs, xy_ubs = get_xy_bounds(total_tree_next)
        push!(lbs, xy_lbs)
        push!(ubs, xy_ubs)

        copy_labels!(total_tree, total_tree_next)
        zero_out!(total_tree_next)

        verbose && (println("[Iteration:] ", i))
    end
    return lbs, ubs
end

""" Sampling-based Sanity Checks """
function sim_traj(x0, tree::KDTREE; nsteps = 50)
    traj = [x0]

    x = copy(x0)

    for i = 1:nsteps
        curr_leaf, curr_lbs, curr_ubs = get_leaf_and_bounds(tree, x)
        next_lbs, next_ubs = reachable_cell(curr_lbs, curr_ubs, curr_leaf.min_control, curr_leaf.max_control)
        x = rand.(Uniform.(next_lbs, next_ubs))
        push!(traj, x)
    end
    return traj
end

function sim_traj_true_next(x0, tree::KDTREE; nsteps = 50)
    traj = [x0]

    x = copy(x0)

    for i = 1:nsteps
        curr_leaf = get_leaf(tree.root_node, x)
        a = rand.(Uniform(curr_leaf.min_control, curr_leaf.max_control))
        xnext, θnext = next_state(x[1], x[2], a)
        x = [xnext, θnext]
        push!(traj, x)
    end
    return traj
end

function find_failures(tree::KDTREE, lbs, ubs; nsteps = 50, npoints = 1000)
    failures = 0
    for i = 1:npoints
        x0 = rand.(Uniform.(lbs, ubs))
        traj = sim_traj(x0, tree, nsteps = nsteps)
        xvals = [traj[i][1] for i = 1:length(traj)]
        if any(xvals .< -10) || any(xvals .> 10)
            failures += 1
        end
    end
    return failures
end

""" Labeling functions """

function label_tree_failures!(tree)
    # Stacks to go through tree
    lb_s = Stack{Vector{Float64}}()
    ub_s = Stack{Vector{Float64}}()
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(lb_s, tree.lbs)
    push!(ub_s, tree.ubs)
    push!(s, tree.root_node)

    while !(isempty(s))
        curr = pop!(s)
        curr_lbs = pop!(lb_s)
        curr_ubs = pop!(ub_s)

        if typeof(curr) == LEAFNODE
            # Check if off runway
            min_cte = curr_lbs[1]
            max_cte = curr_ubs[1]

            curr.prob = (min_cte < -10.0) || (max_cte > 10.0) ? 1.0 : 0.0
        else
            # Traverse tree and keep track of bounds
            dim = curr.dim
            split = curr.split
            # Go left, upper bounds will change
            left_ubs = copy(curr_ubs)
            left_ubs[dim] = split

            push!(lb_s, curr_lbs)
            push!(ub_s, left_ubs)
            push!(s, curr.left)

            # Go right, lower bounds will change
            right_lbs = copy(curr_lbs)
            right_lbs[dim] = split

            push!(lb_s, right_lbs)
            push!(ub_s, curr_ubs)
            push!(s, curr.right)
        end
    end
end

function label_tree_observable!(tree; control_gains = [-0.74, -0.44])
    # Stacks to go through tree
    lb_s = Stack{Vector{Float64}}()
    ub_s = Stack{Vector{Float64}}()
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(lb_s, tree.lbs)
    push!(ub_s, tree.ubs)
    push!(s, tree.root_node)

    while !(isempty(s))
        curr = pop!(s)
        curr_lbs = pop!(lb_s)
        curr_ubs = pop!(ub_s)

        if typeof(curr) == LEAFNODE
            curr.min_control = control_gains' * curr_ubs
            curr.max_control = control_gains' * curr_lbs

            # if curr.min_control > curr.max_control
            #     println(curr_ubs)
            #     println(curr_lbs)
            # end

            # Check if off runway
            min_cte = curr_lbs[1]
            max_cte = curr_ubs[1]

            curr.prob = (min_cte < -10.0) || (max_cte > 10.0) ? 1.0 : 0.0
        else
            # Traverse tree and keep track of bounds
            dim = curr.dim
            split = curr.split
            # Go left, upper bounds will change
            left_ubs = copy(curr_ubs)
            left_ubs[dim] = split

            push!(lb_s, curr_lbs)
            push!(ub_s, left_ubs)
            push!(s, curr.left)

            # Go right, lower bounds will change
            right_lbs = copy(curr_lbs)
            right_lbs[dim] = split

            push!(lb_s, right_lbs)
            push!(ub_s, curr_ubs)
            push!(s, curr.right)
        end
    end
end

function label_tree_with_prob!(tree, prob; control_gains = [-0.74, -0.44])
    # Stacks to go through tree
    lb_s = Stack{Vector{Float64}}()
    ub_s = Stack{Vector{Float64}}()
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(lb_s, tree.lbs)
    push!(ub_s, tree.ubs)
    push!(s, tree.root_node)

    while !(isempty(s))
        curr = pop!(s)
        curr_lbs = pop!(lb_s)
        curr_ubs = pop!(ub_s)

        if typeof(curr) == LEAFNODE
            curr.min_control = control_gains' * curr_ubs
            curr.max_control = control_gains' * curr_lbs

            curr.prob = prob
        else
            # Traverse tree and keep track of bounds
            dim = curr.dim
            split = curr.split
            # Go left, upper bounds will change
            left_ubs = copy(curr_ubs)
            left_ubs[dim] = split

            push!(lb_s, curr_lbs)
            push!(ub_s, left_ubs)
            push!(s, curr.left)

            # Go right, lower bounds will change
            right_lbs = copy(curr_lbs)
            right_lbs[dim] = split

            push!(lb_s, right_lbs)
            push!(ub_s, curr_ubs)
            push!(s, curr.right)
        end
    end
end

function copy_labels!(copyee, copyer)
    # Stacks to go through tree
    s1 = Stack{Union{LEAFNODE, KDNODE}}()
    s2 = Stack{Union{LEAFNODE, KDNODE}}()

    push!(s1, copyee.root_node)
    push!(s2, copyer.root_node)

    while !(isempty(s1))
        curr_copyee = pop!(s1)
        curr_copyer = pop!(s2)

        if typeof(curr_copyee) == LEAFNODE
            curr_copyee.prob = curr_copyer.prob
        else
            # Traverse trees
            push!(s1, curr_copyee.left)
            push!(s1, curr_copyee.right)
            push!(s2, curr_copyer.left)
            push!(s2, curr_copyer.right)
        end
    end
end

function one_out!(tree)
    # Stacks to go through tree
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(s, tree.root_node)

    while !(isempty(s))
        curr = pop!(s)

        if typeof(curr) == LEAFNODE
            curr.prob = 1.0
        else
            # Traverse tree
            push!(s, curr.left)
            push!(s, curr.right)
        end
    end
end

function zero_out!(tree)
    # Stacks to go through tree
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(s, tree.root_node)

    while !(isempty(s))
        curr = pop!(s)

        if typeof(curr) == LEAFNODE
            curr.prob = 0.0
        else
            # Traverse tree
            push!(s, curr.left)
            push!(s, curr.right)
        end
    end
end

function label_xy!(xy_tree, total_tree)
    # Stacks to go through tree
    lb_s = Stack{Vector{Float64}}()
    ub_s = Stack{Vector{Float64}}()
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(lb_s, xy_tree.lbs)
    push!(ub_s, xy_tree.ubs)
    push!(s, xy_tree.root_node)

    while !(isempty(s))
        curr = pop!(s)
        curr_lbs = pop!(lb_s)
        curr_ubs = pop!(ub_s)

        if typeof(curr) == LEAFNODE
            total_lbs = [curr_lbs[1]; total_tree.lbs[2]; curr_lbs[2]]
            total_ubs = [curr_ubs[1]; total_tree.ubs[2]; curr_ubs[2]]
            nodes = get_overlapping_nodes(total_tree.root_node, total_lbs, total_ubs)

            max_prob = maximum([node.prob for node in nodes])
            if max_prob == 1.0
                curr.prob = 1.0
            end
        else
            # Traverse tree and keep track of bounds
            dim = curr.dim
            split = curr.split
            # Go left, upper bounds will change
            left_ubs = copy(curr_ubs)
            left_ubs[dim] = split

            push!(lb_s, curr_lbs)
            push!(ub_s, left_ubs)
            push!(s, curr.left)

            # Go right, lower bounds will change
            right_lbs = copy(curr_lbs)
            right_lbs[dim] = split

            push!(lb_s, right_lbs)
            push!(ub_s, curr_ubs)
            push!(s, curr.right)
        end
    end
end

# Modifies control outputs
function label_xyθ!(total_tree, xθ_tree)
    # Stacks to go through tree
    lb_s = Stack{Vector{Float64}}()
    ub_s = Stack{Vector{Float64}}()
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(lb_s, total_tree.lbs)
    push!(ub_s, total_tree.ubs)
    push!(s, total_tree.root_node)

    while !(isempty(s))
        curr = pop!(s)
        curr_lbs = pop!(lb_s)
        curr_ubs = pop!(ub_s)

        if typeof(curr) == LEAFNODE
            point = (curr_ubs[1:2] .+ curr_lbs[1:2]) ./ 2
            leaf = get_leaf(xθ_tree.root_node, point)
            curr.min_control = leaf.min_control
            curr.max_control = leaf.max_control
        else
            # Traverse tree and keep track of bounds
            dim = curr.dim
            split = curr.split
            # Go left, upper bounds will change
            left_ubs = copy(curr_ubs)
            left_ubs[dim] = split

            push!(lb_s, curr_lbs)
            push!(ub_s, left_ubs)
            push!(s, curr.left)

            # Go right, lower bounds will change
            right_lbs = copy(curr_lbs)
            right_lbs[dim] = split

            push!(lb_s, right_lbs)
            push!(ub_s, curr_ubs)
            push!(s, curr.right)
        end
    end
end

function label_start_states!(tree, lbs, ubs)
    zero_out!(tree)
    nodes = get_overlapping_nodes(tree.root_node, lbs, ubs)
    for node in nodes
        node.prob = 1.0
    end
end

function get_xy_bounds(total_tree)
     # Stacks to go through tree
    lb_s = Stack{Vector{Float64}}()
    ub_s = Stack{Vector{Float64}}()
    s = Stack{Union{LEAFNODE, KDNODE}}()

    push!(lb_s, total_tree.lbs)
    push!(ub_s, total_tree.ubs)
    push!(s, total_tree.root_node)

    lbs = total_tree.ubs[[1, 3]]
    ubs = total_tree.lbs[[1, 3]]

    while !(isempty(s))
        curr = pop!(s)
        curr_lbs = pop!(lb_s)
        curr_ubs = pop!(ub_s)

        if typeof(curr) == LEAFNODE
            if curr.prob == 1.0
                xy_lbs = curr_lbs[[1, 3]]
                xy_ubs = curr_ubs[[1, 3]]

                replace_lb_dims = findall(xy_lbs .< lbs)
                replace_ub_dims = findall(xy_ubs .> ubs)

                if length(replace_lb_dims) > 0
                    lbs[replace_lb_dims] = xy_lbs[replace_lb_dims]
                end

                if length(replace_ub_dims) > 0
                    ubs[replace_ub_dims] = xy_ubs[replace_ub_dims]
                end
            end
        else
            # Traverse tree and keep track of bounds
            dim = curr.dim
            split = curr.split
            # Go left, upper bounds will change
            left_ubs = copy(curr_ubs)
            left_ubs[dim] = split

            push!(lb_s, curr_lbs)
            push!(ub_s, left_ubs)
            push!(s, curr.left)

            # Go right, lower bounds will change
            right_lbs = copy(curr_lbs)
            right_lbs[dim] = split

            push!(lb_s, right_lbs)
            push!(ub_s, curr_ubs)
            push!(s, curr.right)
        end
    end

    return lbs, ubs
end

max_stars_repo_name: sisl/VerifyGAN

max_issues_repo_name: sisl/VerifyGAN

max_forks_repo_name: sisl/VerifyGAN


text:
# jni_md.h
const jint = Cint
#ifdef _LP64 /* 64-bit Solaris */
# typedef long jlong;
const jlong = Clonglong
const jbyte = Cchar

# jni.h

const jboolean = Cuchar
const jchar = Cushort
const jshort = Cshort
const jfloat = Cfloat
const jdouble = Cdouble
const jsize = jint
jprimitive = Union{jboolean, jchar, jshort, jfloat, jdouble, jint, jlong}

struct JavaMetaClass{T}
    ptr::Ptr{Nothing}
end

#The metaclass, sort of equivalent to a the
JavaMetaClass(T, ptr) = JavaMetaClass{T}(ptr)

mutable struct JavaObject{T}
    ptr::Ptr{Nothing}

    #This below is ugly. Once we stop supporting 0.5, this can be replaced by
    # function JavaObject{T}(ptr) where T
    function JavaObject{T}(ptr) where T
        j = new{T}(ptr)
        finalizer(deleteref, j)
        return j
    end

    #replace with: JavaObject{T}(argtypes::Tuple, args...) where T
    JavaObject{T}(argtypes::Tuple, args...) where {T} = jnew(T, argtypes, args...)
end

JavaObject(T, ptr) = JavaObject{T}(ptr)
JavaObject{T}() where {T} = JavaObject{T}((),)

function deleteref(x::JavaObject)
    if x.ptr == C_NULL; return; end
    if (penv==C_NULL); return; end
    #ccall(:jl_,Nothing,(Any,),x)
    ccall(jnifunc.DeleteLocalRef, Nothing, (Ptr{JNIEnv}, Ptr{Nothing}), penv, x.ptr)
    x.ptr=C_NULL #Safety in case this function is called direcly, rather than at finalize
    return
end

"""
```
isnull(obj::JavaObject)
```
Checks if the passed JavaObject is null or not

### Args
* obj: The object of type JavaObject

### Returns
true if the passed object is null else false
"""
isnull(obj::JavaObject) = obj.ptr == C_NULL

"""
```
isnull(obj::JavaMetaClass)
```
Checks if the passed JavaMetaClass is null or not

### Args
* obj: The object of type JavaMetaClass

### Returns
true if the passed object is null else false
"""
isnull(obj::JavaMetaClass) = obj.ptr == C_NULL

const JClass = JavaObject{Symbol("java.lang.Class")}
const JObject = JavaObject{Symbol("java.lang.Object")}
const JMethod = JavaObject{Symbol("java.lang.reflect.Method")}
const JThread = JavaObject{Symbol("java.lang.Thread")}
const JClassLoader = JavaObject{Symbol("java.lang.ClassLoader")}
const JString = JavaObject{Symbol("java.lang.String")}

function JString(str::AbstractString)
    jstring = ccall(jnifunc.NewStringUTF, Ptr{Nothing}, (Ptr{JNIEnv}, Ptr{UInt8}), penv, String(str))
    if jstring == C_NULL
        geterror()
    else
        return JString(jstring)
    end
end

# jvalue(v::Integer) = int64(v) << (64-8*sizeof(v))
jvalue(v::Integer) = Int64(v)
jvalue(v::Float32) = jvalue(reinterpret(Int32, v))
jvalue(v::Float64) = jvalue(reinterpret(Int64, v))
jvalue(v::Ptr) = jvalue(Int(v))

function _jimport(juliaclass)
    for str ∈ [" ", "(", ")"]
        juliaclass = replace(juliaclass, str=>"")
    end
    :(JavaObject{Symbol($juliaclass)})
end

macro jimport(class::Expr)
    juliaclass = sprint(Base.show_unquoted, class)
    _jimport(juliaclass)
end
macro jimport(class::Symbol)
    juliaclass = string(class)
    _jimport(juliaclass)
end
macro jimport(class::AbstractString)
    _jimport(class)
end

function jnew(T::Symbol, argtypes::Tuple, args...)
    assertroottask_or_goodenv()
    sig = method_signature(Nothing, argtypes...)
    jmethodId = ccall(jnifunc.GetMethodID, Ptr{Nothing},
                      (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(T),
                      String("<init>"), sig)
    if jmethodId == C_NULL
        throw(JavaCallError("No constructor for $T with signature $sig"))
    end
    return  _jcall(metaclass(T), jmethodId, jnifunc.NewObjectA, JavaObject{T}, argtypes, args...)
end

# Call static methods
function jcall(typ::Type{JavaObject{T}}, method::AbstractString, rettype::Type, argtypes::Tuple,
               args... ) where T
    assertroottask_or_goodenv()
    sig = method_signature(rettype, argtypes...)
    jmethodId = ccall(jnifunc.GetStaticMethodID, Ptr{Nothing},
                      (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(T),
                      String(method), sig)
    jmethodId==C_NULL && geterror(true)
    _jcall(metaclass(T), jmethodId, C_NULL, rettype, argtypes, args...)
end

# Call instance methods
function jcall(obj::JavaObject, method::AbstractString, rettype::Type, argtypes::Tuple, args... )
    assertroottask_or_goodenv()
    sig = method_signature(rettype, argtypes...)
    jmethodId = ccall(jnifunc.GetMethodID, Ptr{Nothing},
                      (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(obj),
                      String(method), sig)
    jmethodId==C_NULL && geterror(true)
    _jcall(obj, jmethodId, C_NULL, rettype,  argtypes, args...)
end

function jfield(typ::Type{JavaObject{T}}, field::AbstractString, fieldType::Type) where T
    assertroottask_or_goodenv()
    jfieldID  = ccall(jnifunc.GetStaticFieldID, Ptr{Nothing},
                      (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(T),
                      String(field), signature(fieldType))
    jfieldID==C_NULL && geterror(true)
    _jfield(metaclass(T), jfieldID, fieldType)
end

function jfield(obj::JavaObject, field::AbstractString, fieldType::Type)
    assertroottask_or_goodenv()
    jfieldID  = ccall(jnifunc.GetFieldID, Ptr{Nothing},
                      (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}), penv, metaclass(obj),
                      String(field), signature(fieldType))
    jfieldID==C_NULL && geterror(true)
    _jfield(obj, jfieldID, fieldType)
end

for (x, y, z) in [(:jboolean, :(jnifunc.GetBooleanField), :(jnifunc.GetStaticBooleanField)),
                  (:jchar, :(jnifunc.GetCharField), :(jnifunc.GetStaticCharField)),
                  (:jbyte, :(jnifunc.GetByteField), :(jnifunc.GetStaticBypeField)),
                  (:jshort, :(jnifunc.GetShortField), :(jnifunc.GetStaticShortField)),
                  (:jint, :(jnifunc.GetIntField), :(jnifunc.GetStaticIntField)),
                  (:jlong, :(jnifunc.GetLongField), :(jnifunc.GetStaticLongField)),
                  (:jfloat, :(jnifunc.GetFloatField), :(jnifunc.GetStaticFloatField)),
                  (:jdouble, :(jnifunc.GetDoubleField), :(jnifunc.GetStaticDoubleField)) ]

    m = quote
        function _jfield(obj, jfieldID::Ptr{Nothing}, fieldType::Type{$(x)})
            callmethod = ifelse( typeof(obj)<:JavaObject, $y , $z )
            result = ccall(callmethod, $x, (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{Nothing}), penv, obj.ptr,
                           jfieldID)
            result==C_NULL && geterror()
            return convert_result(fieldType, result)
        end
    end
    eval(m)
end

function _jfield(obj, jfieldID::Ptr{Nothing}, fieldType::Type)
    callmethod = ifelse( typeof(obj)<:JavaObject, jnifunc.GetObjectField , jnifunc.GetStaticObjectField )
    result = ccall(callmethod, Ptr{Nothing}, (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{Nothing}), penv, obj.ptr,
                   jfieldID)
    result==C_NULL && geterror()
    return convert_result(fieldType, result)
end

#Generate these methods to satisfy ccall's compile time constant requirement
#_jcall for primitive and Nothing return types
for (x, y, z) in [ (:jboolean, :(jnifunc.CallBooleanMethodA), :(jnifunc.CallStaticBooleanMethodA)),
                  (:jchar, :(jnifunc.CallCharMethodA), :(jnifunc.CallStaticCharMethodA)),
                  (:jbyte, :(jnifunc.CallByteMethodA), :(jnifunc.CallStaticByteMethodA)),
                  (:jshort, :(jnifunc.CallShortMethodA), :(jnifunc.CallStaticShortMethodA)),
                  (:jint, :(jnifunc.CallIntMethodA), :(jnifunc.CallStaticIntMethodA)),
                  (:jlong, :(jnifunc.CallLongMethodA), :(jnifunc.CallStaticLongMethodA)),
                  (:jfloat, :(jnifunc.CallFloatMethodA), :(jnifunc.CallStaticFloatMethodA)),
                  (:jdouble, :(jnifunc.CallDoubleMethodA), :(jnifunc.CallStaticDoubleMethodA)),
                  (:Nothing, :(jnifunc.CallVoidMethodA), :(jnifunc.CallStaticVoidMethodA)) ]
    m = quote
        function _jcall(obj, jmethodId::Ptr{Nothing}, callmethod::Ptr{Nothing}, rettype::Type{$(x)},
                        argtypes::Tuple, args...)
            if callmethod == C_NULL #!
                callmethod = ifelse( typeof(obj)<:JavaObject, $y , $z )
            end
            @assert callmethod != C_NULL
            @assert jmethodId != C_NULL
            isnull(obj) && throw(JavaCallError("Attempt to call method on Java NULL"))
            savedArgs, convertedArgs = convert_args(argtypes, args...)
            result = ccall(callmethod, $x , (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}), penv, obj.ptr, jmethodId, convertedArgs)
            result==C_NULL && geterror()
            result == nothing && (return)
            return convert_result(rettype, result)
        end
    end
    eval(m)
end

#_jcall for Object return types
#obj -- receiver - Class pointer or object prointer
#jmethodId -- Java method ID
#callmethod -- the C method pointer to call
function _jcall(obj, jmethodId::Ptr{Nothing}, callmethod::Ptr{Nothing}, rettype::Type, argtypes::Tuple,
                args...)
    if callmethod == C_NULL
        callmethod = ifelse(typeof(obj)<:JavaObject, jnifunc.CallObjectMethodA ,
                            jnifunc.CallStaticObjectMethodA)
    end
    @assert callmethod != C_NULL
    @assert jmethodId != C_NULL
    isnull(obj) && error("Attempt to call method on Java NULL")
    savedArgs, convertedArgs = convert_args(argtypes, args...)
    result = ccall(callmethod, Ptr{Nothing}, (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}),
                   penv, obj.ptr, jmethodId, convertedArgs)
    result==C_NULL && geterror()
    return convert_result(rettype, result)
end

global const _jmc_cache = Dict{Symbol, JavaMetaClass}()

function _metaclass(class::Symbol)
    jclass=javaclassname(class)
    jclassptr = ccall(jnifunc.FindClass, Ptr{Nothing}, (Ptr{JNIEnv}, Ptr{UInt8}), penv, jclass)
    jclassptr == C_NULL && throw(JavaCallError("Class Not Found $jclass"))
    return JavaMetaClass(class, jclassptr)
end

function metaclass(class::Symbol)
    if !haskey(_jmc_cache, class)
        _jmc_cache[class] = _metaclass(class)
    end
    return _jmc_cache[class]
end

metaclass(::Type{JavaObject{T}}) where {T} = metaclass(T)
metaclass(::JavaObject{T}) where {T} = metaclass(T)

javaclassname(class::Symbol) = replace(string(class), "."=>"/")

function geterror(allow=false)
    isexception = ccall(jnifunc.ExceptionCheck, jboolean, (Ptr{JNIEnv},), penv )

    if isexception == JNI_TRUE
        jthrow = ccall(jnifunc.ExceptionOccurred, Ptr{Nothing}, (Ptr{JNIEnv},), penv)
        jthrow==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        ccall(jnifunc.ExceptionDescribe, Nothing, (Ptr{JNIEnv},), penv ) #Print java stackstrace to stdout
        ccall(jnifunc.ExceptionClear, Nothing, (Ptr{JNIEnv},), penv )
        jclass = ccall(jnifunc.FindClass, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{UInt8}), penv,
                       "java/lang/Throwable")
        jclass==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        jmethodId=ccall(jnifunc.GetMethodID, Ptr{Nothing},
                        (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}), penv, jclass, "toString",
                        "()Ljava/lang/String;")
        jmethodId==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        res = ccall(jnifunc.CallObjectMethodA, Ptr{Nothing},
                    (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}), penv, jthrow, jmethodId,
                    C_NULL)
        res==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        msg = unsafe_string(JString(res))
        ccall(jnifunc.DeleteLocalRef, Nothing, (Ptr{JNIEnv}, Ptr{Nothing}), penv, jthrow)
        throw(JavaCallError(string("Error calling Java: ",msg)))
    else
        if allow==false
            return #No exception pending, legitimate NULL returned from Java
        else
            throw(JavaCallError("Null from Java. Not known how"))
        end
    end
end

#get the JNI signature string for a method, given its
#return type and argument types
function method_signature(rettype, argtypes...)
    s=IOBuffer()
    write(s, "(")
    for arg in argtypes
        write(s, signature(arg))
    end
    write(s, ")")
    write(s, signature(rettype))
    return String(take!(s))
end

#get the JNI signature string for a given type
function signature(arg::Type)
    if arg === jboolean
        return "Z"
    elseif arg === jbyte
        return "B"
    elseif arg === jchar
        return "C"
    elseif arg === jshort
        return "S"
    elseif arg === jint
        return "I"
    elseif arg === jlong
        return "J"
    elseif arg === jfloat
        return "F"
    elseif arg === jdouble
        return "D"
    elseif arg === Nothing
        return "V"
    elseif arg <: Array
        dims = "[" ^ ndims(arg)
        return string(dims, signature(eltype(arg)))
    end
end

signature(arg::Type{JavaObject{T}}) where {T} = string("L", javaclassname(T), ";")


text:
module OPFLearn

import Dates
using DelimitedFiles: readdlm, writedlm
import Distributed
import LinearAlgebra: diag, eigvals, I, cholesky, norm, Hermitian, dot, norm, logdet

import OnlineStatsBase: fit!, CovMatrix, Variance
import OnlineStats: Series
import StatsBase: sample

import JuMP
import MathOptInterface
const MOI = MathOptInterface
import PowerModels
import PowerModels: ref, ids, var
const PM = PowerModels

import Random

import Ipopt

include("create_samples.jl")
include("create_samples_distributed.jl")
include("tasks.jl")

include("formulations.jl")
include("run_ac_opf.jl")

include("network.jl")
include("results.jl")
include("io.jl")
include("sample.jl")
include("export.jl")

const TOL = 1e-6
const INFEASIBILITY_CERT_SHIFT = 0.0  # Tune this parameter based on willingness to cut out large loads
const FNFP_SHIFT = 1e-2  # Increase to find more AC-OPF feasible points as FNFP problem
const R_TOLERANCE = 1e-5  # 10x tolerance used for solving Nearest Feasible load problem

const NW_DEFAULT = PM.nw_id_default # The first network index for a PowerModels network.
				      # OPFLearn currently only supports single network models.

const SOLVED_STATUSES = [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_OPTIMAL, MOI.ALMOST_LOCALLY_SOLVED]

const EPS = 1e-5
const STAT_ORDER = ["pd", "qd", "pg", "vm_gen"]
const DEFAULT_INPUTS = ["pd", "qd"]
const DEFAULT_OUTPUTS = ["p_gen", "q_gen", "vm_gen", "vm_bus", "va_bus",
						 "p_to", "q_to", "p_fr", "q_fr"]
# "vm" and "va" can be added for voltage magnitudes and angles at each bus
const DEFAULT_DUALS = ["v_min", "v_max", "pg_min", "pg_max", "qg_min", "qg_max",
					   "pto_max", "qto_max", "pfr_max", "qfr_max"]

const DEFAULT_SAVE_ORDER = vcat(DEFAULT_INPUTS, DEFAULT_OUTPUTS, DEFAULT_DUALS)

const ELEMENT_LABELS = Dict(zip(append!(DEFAULT_SAVE_ORDER, ["v_bus", "total_cost"]),
				    ["load", "load",
					 "gen", "gen", "gen", "bus", "bus",
					 "line", "line", "line", "line",
					 "bus", "bus", "gen", "gen", "gen", "gen",
					 "line", "line", "line", "line",
					 "bus", "net"]))

const VAR_SIZES = Dict([
					 ("v_min", [:n_buses]),
					 ("v_max", [:n_buses]),
					 ("pg_min", [:n_gens, :n_ext]),
					 ("pg_max", [:n_gens, :n_ext]),
					 ("qg_min", [:n_gens, :n_ext]),
					 ("qg_max", [:n_gens, :n_ext]),
					 ("pfr_max", [:n_branches]),
					 ("pto_max", [:n_branches]),
					 ("qfr_max", [:n_branches]),
					 ("qto_max", [:n_branches]),
					 ("pd", [:n_loads]),
					 ("qd", [:n_loads]),
					 ("p_gen", [:n_gens, :n_ext]),
					 ("q_gen", [:n_gens, :n_ext]),
					 ("v_bus", [:n_buses]),
					 ("vm_gen", [:n_gens, :n_ext]),
					 ("vm_bus", [:n_buses]),
					 ("va_bus", [:n_buses]),
					 ("p_fr", [:n_branches]),
					 ("p_to", [:n_branches]),
					 ("q_fr", [:n_branches]),
					 ("q_to", [:n_branches]),
					 ("q_to", [:n_branches]),
					 ("total_cost", [:n_nets]),
					 ])

end


text:
#!/usr/bin/env julia

# active project directory
import Pkg
Pkg.activate(dirname(@__FILE__))

# load modules
using ArgParse
import Pkg.TOML
import LibGit2
import UUIDs
import Pkg.Types: VersionRange, VersionBound

include(joinpath(@__DIR__, "registry_parsing.jl"))

const NAME = "name"
const REPO = "repo"
const DEPS = "deps"
const PKGS = "packages"
const COMPAT = "compat"
const PROJECT_FILE  = "Project.toml"
const REGISTRY_FILE = "Registry.toml"
const USER_DEPO = first(Pkg.depots())

function isrepo(pkgdir::String)
    try
        repo = LibGit2.GitRepo(pkgdir)
        close(repo)
    catch
        return false
    end
    return true
end

function readregistry(regpath::String)
    regfile = joinpath(regpath, REGISTRY_FILE)
    !isfile(regfile) && error("Registry configuration is not found in $regpath")
    reg = TOML.parsefile(regfile)
    return TOML.parsefile(regfile)
end

function readproject(pkgdir::String)
    prjfile = joinpath(pkgdir, PROJECT_FILE)
    !isfile(prjfile) && error("Package must have `Project.toml` file")
    return TOML.parsefile(prjfile)
end

function genpackage(pkgdir::String, project::Dict{T,Any}) where T<:AbstractString
    package = Dict{T, T}()
    package[NAME] = if haskey(project, NAME)
        project[NAME]
    else
        basename(PKGSRC)
    end
    package["uuid"] = if haskey(project, "uuid")
        project["uuid"]
    else
        string(UUIDs.uuid1())
    end
    package["repo"] = LibGit2.with(LibGit2.GitRepo(pkgdir)) do repo
        LibGit2.with(LibGit2.get(LibGit2.GitRemote, repo, "origin")) do remote
            LibGit2.url(remote)
        end
    end
    return package
end

function genversions(pkgdir)
    versions = Dict{String, Any}()
    LibGit2.with(LibGit2.GitRepo(pkgdir)) do repo
        tags = LibGit2.tag_list(repo)
        for tag in tags
            ignore_version = false
            tagref = LibGit2.GitReference(repo, "refs/tags/$tag")
            tree = LibGit2.peel(LibGit2.GitTree, tagref)
            try
                cont = LibGit2.content(tree["Project.toml"])
            catch ex
                if !isa(ex, KeyError)
                    rethrow(ex)
                else
                    ignore_version = true
                end
            end
            if ignore_version
                println("Ignoring version $tag, no Project.toml")
                continue
            end
            treesha = LibGit2.GitHash(tree)
            v = string(VersionNumber(tag))
            versions[v] = Dict("git-tree-sha1" => string(treesha))
            println("Adding version $v")
        end
    end
    return versions
end

function group(A::Vector{T}, sentinel::T=last(A)) where T
    p = sortperm(A)
    q = A[p]
    res = Vector{Vector{Int}}()
    grp = Vector{Int}()
    first = true
    last = sentinel
    for (i,v) in enumerate(q)
        if !first && last != v
            push!(res, grp)
            grp = [p[i]]
        else
            push!(grp, p[i])
        end
        last = v
        first = false
    end
    push!(res, grp)
    return res
end

"""Generate dependencies from `REQUIRE` file"""
function gendependenciesreq(pkgdir::String, versions::Dict{T,Any}) where T<:AbstractString
    # read all dependencies from repo
    depvers = Dict{T,Vector{VersionNumber}}()
    verhash = [k => versions[k]["git-tree-sha1"] for k in sort(collect(keys(versions)))]
    LibGit2.with(LibGit2.GitRepo(pkgdir)) do repo
        for (ver, hash) in verhash
            LibGit2.with(LibGit2.GitTree(repo, "$hash^{tree}")) do tree
                try
                    cont = LibGit2.content(tree["REQUIRE"])
                    for req in split(cont, '\n', keepempty=false)
                        dep = first(split(req, ' '))
                        #dep == "julia" && continue
                        !haskey(depvers, dep) && setindex!(depvers, VersionNumber[], dep)
                        push!(depvers[dep], VersionNumber(ver))
                    end
                catch y
                    if !isa(y, KeyError)
                        rethrow(y)
                    end
                end
            end
        end
    end

    # aggregate dependencies by versions
    deps = Dict{T,Dict{T,T}}()

    # get registry context
    ctx = Pkg.Types.Context()
    Pkg.Types.Context!(ctx)

    # group dependencies by min-max versions
    verdeps = [extrema(vers) => dep for (dep, vers) in depvers if dep != "julia"]
    gidxs = group(map(first, verdeps))
    for g in gidxs
        # find uuids of packages
        pkgspecs = map(Pkg.Types.PackageSpec, map(last, verdeps[g]))
        Pkg.Types.project_deps_resolve!(ctx.env, pkgspecs)
        Pkg.Types.registry_resolve!(ctx.env, pkgspecs)
        # form version range
        vv = first(first(verdeps[g]))
        vrange = VersionRange(map(VersionBound, vv)...) |> string
        deps[vrange] = Dict{T,T}()
        for ps in pkgspecs
            deps[vrange][ps.name] = string(ps.uuid)
        end
    end

    return deps, depvers
end

function gendependencies(pkgdir::String, versions::Dict{T,Any}) where T<:AbstractString
    # read all dependencies from repo
    depvers = Dict{T,Any}()
    verhash = [k => versions[k]["git-tree-sha1"] for k in sort(collect(keys(versions)))]
    LibGit2.with(LibGit2.GitRepo(pkgdir)) do repo
        for (ver, hash) in verhash
            LibGit2.with(LibGit2.GitTree(repo, "$hash^{tree}")) do tree
                try
                    cont = LibGit2.content(tree["Project.toml"])
                    prj = TOML.parse(cont)
                    depvers[ver] = prj[DEPS]
                catch y
                    if !isa(y, KeyError)
                        rethrow(y)
                    end
                end
            end
        end
    end
    return depvers
end

function aggregateverdeps(verdeps)
    # aggregate dependencies by versions
    cdeps = Dict{String, Any}()

    # group dependencies by min-max versions
    gidxs = group(map(first, verdeps))
    for g in gidxs
        vv = first(first(verdeps[g])) # form version range
        vrange = VersionRange(map(VersionBound, vv)...) |> string
        cdeps[vrange] = Dict{String,String}()
        for i in g
            pkg = last(verdeps[i])
            cdeps[vrange][first(pkg)] = last(pkg)
        end
    end
    return cdeps
end

function compactdeps(deps::Dict{T,Any}) where T<:AbstractString
    alldeps = vcat(([ (d,uuid) for (d,uuid) in ds] for (v, ds) in deps)...) |> unique
    depvers = [ d => [ VersionNumber(k) for (k, depid) in deps if haskey(depid, first(d))] for d in alldeps]
    verdeps = [extrema(vers) => dep for (dep, vers) in depvers]
    return aggregateverdeps(verdeps)
end

function gencompatibility(pkgdir::String, versions::Dict{T,Any}) where T<:AbstractString
    # read all dependencies from repo
    compatvers = Dict{T,Any}()
    verhash = [k => versions[k]["git-tree-sha1"] for k in sort(collect(keys(versions)))]
    LibGit2.with(LibGit2.GitRepo(pkgdir)) do repo
        for (ver, hash) in verhash
            LibGit2.with(LibGit2.GitTree(repo, "$hash^{tree}")) do tree
                try
                    prjfile = LibGit2.content(tree["Project.toml"])
                    prj = TOML.parse(prjfile)
                    compatvers[ver] = prj[COMPAT]
                catch y
                    if !isa(y, KeyError)
                        rethrow(y)
                    end
                end
            end
        end
    end
    return compatvers
end

function compactcompats(compats::Dict{T,Any}) where T<:AbstractString
    cptverschange = vcat(([dver => v for dver in dvers] for (v, dvers) in compats)...)
    length(cptverschange) == 0 && return Dict{String,Any}() # no compat info
    gidxs = group(map(first, cptverschange))
    vercpts = [extrema(map(p->VersionNumber(last(p)), cptverschange[g])) => first(cptverschange[first(g)]) for g in gidxs]

    return aggregateverdeps(vercpts)
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "init"
            help = "initialize a registry"
            action = :command
        "update"
            help = "update a registry"
            action = :command
        "add"
            help = "add a package to a registry"
            action = :command
        "addurl"
            help = "add a package to a registry from url"
            action = :command
    end

    @add_arg_table s["init"] begin
        "path"
            arg_type = String
            help = "a registry location"
            required = true
    end

    @add_arg_table s["add"] begin
        "registry"
            arg_type = String
            help = "a registry location"
            required = true
        "pkg"
            arg_type = String
            help = "a package location"
            required = true
    end

    @add_arg_table s["addurl"] begin
        "registry"
            arg_type = String
            help = "a registry location"
            required = true
        "pkgurl"
            arg_type = String
            help = "a package location"
            required = true
    end

    @add_arg_table s["update"] begin
        "registry"
            arg_type = String
            help = "a registry location"
            required = true
    end

    return parse_args(s)
end

isregistrydir(path::String) = isdir(path) && isrepo(path) && isfile(joinpath(path, REGISTRY_FILE))

function getregistrypath(path::String)
    if isdirpath(path)
        return path, ""
    else
        return joinpath(USER_DEPO, "registries", path), path
    end
end

function writeto(f, fpath)
    io = open(fpath, "w")
    try
        f(io)
    finally
        close(io)
    end
end

function saveregistryfile(path, reg)
    name = reg[NAME]
    uuid = reg["uuid"]
    repo = reg[REPO]
    description = reg["description"]
    packages = reg[PKGS]
    data =
    regdata = RegistryData(name, UUID(uuid), repo, description, packages, Dict())
    writeto(io->TOML.print(io, regdata), path)
end

function init(path::String)
    regpath, regname = getregistrypath(path)

    if !isregistrydir(regpath)
        # input registry properties
        reg = Dict{String, Any}()
        reg["uuid"] = UUIDs.uuid4()
        print("Enter registry name$(isempty(regname) ? ":" : " ["*regname*"]:") ")
        tmpname = readline(stdin)
        reg[NAME] = isempty(tmpname) ? regname : tmpname
        print("Enter registry repository url: ")
        reg[REPO] = readline(stdin)
        reg[PKGS] = Dict{String,Any}()
        print("Enter registry description: ")
        reg["description"] = readline(stdin)
        print("Create registry [Y/n]: ")
        readline(stdin) == "n" && return

        # create directory
        mkpath(regpath)
        # create repo
        repo = LibGit2.init(regpath)
        try
            # generate Registry.toml
            regfilepath = joinpath(regpath, REGISTRY_FILE)
            saveregistryfile(regfilepath, reg)

            LibGit2.add!(repo, REGISTRY_FILE)
            LibGit2.commit(repo, "Registry created.")
            rmt = LibGit2.GitRemote(repo, "origin", reg[REPO])
            try
                LibGit2.add_push!(repo, rmt, "refs/heads/master")
                LibGit2.add_fetch!(repo, rmt, "refs/heads/master")
            finally
                close(rmt)
            end
            LibGit2.set_remote_url(repo, "origin", reg[REPO])

            println("Registry is created in $regpath")
        catch ex
            rm(regpath, recursive=true, force=true)
            rethrow(ex)
        finally
            close(repo)
        end
    else
        error("Registry is already created in $path.")
    end
end

function addjuliavers!(compats, vers)
    svers = sort!(collect(keys(vers)))
    # if no compatibility information exists all versions are possible
    lastver = svers[end]
    if !haskey(compats, lastver)
        compats[lastver] = Dict{String,Any}()
    end
    if !haskey(compats[lastver], "julia")
        compats[lastver]["julia"] = "*"
    end
end

function addurl(regdir::String, pkgurl::String)
    # clone to tempdir first
    tmp = mktempdir()
    Base.shred!(LibGit2.CachedCredentials()) do creds
        LibGit2.with(Pkg.GitTools.clone(pkgurl, tmp; header = "Pkg from $(repr(pkgurl))", credentials = creds)) do repo
        end
    end

    # verify it looks like a package
    if !isfile(joinpath(tmp, "Project.toml"))
        Pkg.Types.pkgerror("no `Project.toml` file in cloned package")
    end

    addpkg(regdir, tmp)
end

function addpkg(regdir::String, pkgdir::String)
    Pkg.Registry.update(regdir)
    regpath, regname = getregistrypath(regdir)

    # load registry & package config
    reg = readregistry(regpath)
    prj = readproject(pkgdir)
    prjname = prj[NAME]
    prjid = prj["uuid"]

    # get package versions
    vers = genversions(pkgdir)
    if length(vers) == 0
        @info "Package $prjname doesn't have any version git tags."
        return
    end

    # create registry record
    prjpath = joinpath(regpath, string(prjname[1])*"/"*prjname)
    # isdir(prjpath) && error("Package `$(prjname)` is already added to the registry `$(reg[NAME])`")
    mkpath(prjpath)

    # write package description
    writeto(joinpath(prjpath, "Package.toml")) do io
        prjdesc = Dict{String, Any}()
        prjdesc[NAME] = prjname
        prjdesc["uuid"] = prjid

        repo = LibGit2.GitRepo(pkgdir)
        try
            rmt = LibGit2.lookup_remote(repo, "origin")
            repourl = LibGit2.url(rmt)
            close(rmt)
            if startswith(repourl, "git") # convert to http
                repourl = replace(repourl, ":" =>"/")
                repourl = replace(repourl, "git@" =>"https://")
            end
            print("Confirm package remote repository location [$repourl]: ")
            tmpurl = readline(stdin)
            prjdesc[REPO] = isempty(tmpurl) ? repourl : tmpurl
        catch err
            rethrow(err)
        finally
            close(repo)
        end
        TOML.print(io, prjdesc)
    end

    # write package versions
    writeto(joinpath(prjpath, "Versions.toml")) do io
        TOML.print(io, vers, sorted=true)
    end

    # write compatibility reqs
    writeto(joinpath(prjpath, "Compat.toml")) do io
        compats = gencompatibility(pkgdir, vers)
        addjuliavers!(compats, vers)
        TOML.print(io, compactcompats(compats), sorted=true)
    end

    # write dependencies
    writeto(joinpath(prjpath, "Deps.toml")) do io
        deps = gendependencies(pkgdir, vers)
        TOML.print(io, compactdeps(deps))
    end

    # generate record in a Registry.toml
    if !haskey(reg, PKGS)
        reg[PKGS] = Dict{String, Any}()
    end
    reg[PKGS][prjid] = Dict{String, Any}()
    reg[PKGS][prjid]["name"] = prjname
    reg[PKGS][prjid]["path"] = string(prjname[1])*"/"*prjname
    saveregistryfile(joinpath(regpath, REGISTRY_FILE), reg)

    # commit changes
    repo = LibGit2.init(regpath)
    try
        LibGit2.add!(repo, string(prjname[1])*"/"*prjname)
        LibGit2.add!(repo, REGISTRY_FILE)
        LibGit2.commit(repo, "Added package $prjname to the registry.")
        LibGit2.push(repo, remote="origin", refspecs=["refs/heads/master"])
    finally
        close(repo)
    end

    println("$prjname added to registry in $regpath")
end

function main()
    args = parse_commandline()

    # process commands
    cmd = args["%COMMAND%"]
    if cmd == "init"
        init(args["init"]["path"])
    elseif cmd == "add"
        addargs = args["add"]
        addpkg(addargs["registry"], addargs["pkg"])
    elseif cmd == "addurl"
        addargs = args["addurl"]
        addurl(addargs["registry"], addargs["pkgurl"])
    elseif cmd == "update"
        updreg(args["update"]["registry"])
    else
        error("Unknow command `$cmd`")
    end
end

!isinteractive() && main()

max_stars_repo_name: sisl/registratura

max_issues_repo_name: sisl/registratura

max_forks_repo_name: sisl/registratura


text:
@assert 1 == 1
@assert 1 === 1.0 # ERROR: AssertionError: 1 === 1.0

@time sleep(1)

@code_llvm 1 == 1.0

@code_native 1 == 1.0


text:
## We compare here three different solvers : SCF, direct minimization and newton
## algorithm.

using DFTK, LinearAlgebra

a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]]

model = model_LDA(lattice, atoms)
basis = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3])
tol = 1e-12

## SCF
scfres_scf = self_consistent_field(basis, tol=tol,
                                   is_converged=DFTK.ScfConvergenceDensity(tol))

## Direct minimization
scfres_dm = direct_minimization(basis, tol=tol)

## Newton algorithm
# start not too far from the solution to ensure convergence : we use here the
# solution of a single iteration SCF
scfres_start = self_consistent_field(basis, maxiter=1)
# remove virtual orbitals
ψ, _ = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation)
scfres_newton = newton(basis, ψ, tol=tol)

println("|ρ_newton - ρ_scf| = ", norm(scfres_newton.ρ - scfres_scf.ρ))
println("|ρ_newton - ρ_dm| = ", norm(scfres_newton.ρ - scfres_dm.ρ))
println("|ρ_scf - ρ_dm| = ", norm(scfres_scf.ρ - scfres_dm.ρ))


text:
# linalg tests take the longest - start them off first
testnames = [
    "linalg", "core", "keywordargs", "numbers", "strings",
    "collections", "hashing", "remote", "iobuffer", "arrayops", "reduce", "reducedim",
    "simdloop", "blas", "fft", "dsp", "sparse", "bitarray", "random", "math",
    "functional", "bigint", "sorting", "statistics", "spawn",
    "backtrace", "priorityqueue", "arpack", "file", "suitesparse", "version",
    "resolve", "pollfd", "mpfr", "broadcast", "complex", "socket",
    "floatapprox", "readdlm", "regex", "float16", "combinatorics",
    "sysinfo", "rounding", "ranges", "mod2pi", "euler", "show",
    "lineedit", "replcompletions", "repl", "test", "examples", "goto"
]
@unix_only push!(testnames, "unicode")

# parallel tests depend on other workers - do them last
push!(testnames, "parallel")

tests = ARGS==["all"] ? testnames : ARGS

if "linalg" in tests
    # specifically selected case
    filter!(x -> x != "linalg", tests)
    prepend!(tests, ["linalg1", "linalg2", "linalg3", "linalg4"])
end

net_required_for = ["socket", "parallel"]
net_on = true
try
    getipaddr()
catch
    warn("Networking unavailable: Skipping tests [" * join(net_required_for, ", ") * "]")
    net_on = false
end

n = 1
if net_on
    n = min(8, CPU_CORES, length(tests))
    n > 1  && addprocs(n; exeflags=`--check-bounds=yes`)
    blas_set_num_threads(1)
else
    filter!(x -> !(x in net_required_for), tests)
end

@everywhere include("testdefs.jl")

reduce(propagate_errors, nothing, pmap(runtests, tests; err_retry=false, err_stop=true))

@unix_only n > 1 && rmprocs(workers(), waitfor=5.0)
println("    \033[32;1mSUCCESS\033[0m")

max_stars_repo_name: Aerlinger/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia