text:
# Rock Diagnosis
# Based on Araya-Lopez 2013 (Active Diagnosis Through Information-Lookahead Planning)

struct RockIPOMDP <: IPOMDP{Int, Symbol, Symbol}
    tdist::Array
    tarray::Array{Float64}
    odist::Array
    oarray::Array{Float64}
    inforeward::Vector{Vector{Float64}}
    discount::Float64
end

struct RockRIPOMDP <: RIPOMDP{Int, Symbol, Symbol}
    tarrayl::Array{Float64}
    tarrayu::Array{Float64}
    oarrayl::Array{Float64}
    oarrayu::Array{Float64}
    inforeward::Vector{Vector{Float64}}
    discount::Float64
    usize::Float64
end

# rock loation, starting location, rock status, and max sensor efficiency
const lrock = 2
const lstart = 1
const srock = :bad
const seff = 0.9

const discount_rock = 0.99
const states_rock = [10, 11, 20, 21] # location: 1 or 2, rock status: 0 (bad) or 1 (good)
const actions_rock = [:left, :right, :check]
const observations_rock = [:good, :bad, :none]
const inforeward_rock = [[1.0, -1/3, -1/3, -1/3],
                         [-1/3, 1.0, -1/3, -1/3],
                         [-1/3, -1/3, 1.0, -1/3],
                         [-1/3, -1/3, -1/3, 1.0]]
const initbelief_rock = [0.5, 0.5, 0.0, 0.0]
const initbelief_dist_rock = SparseCat(states_rock, initbelief_rock)
const initdist_rock = SparseCat(states_rock, [0.5, 0.5, 0.0, 0.0])
const pϵ_rock = 1e-6

# Nominal transition function distributions (deterministic)
const tdist_rock = [hcat(SparseCat(states_rock, [1.0, 0.0, 0.0, 0.0]),  #1B, L to 1B
                    SparseCat(states_rock, [0.0, 0.0, 1.0, 0.0]),  #1B, R to 2B
                    SparseCat(states_rock, [1.0, 0.0, 0.0, 0.0])); #1B, C to 1B
                    hcat(SparseCat(states_rock, [0.0, 1.0, 0.0, 0.0]),  #1G, L to 1G
                    SparseCat(states_rock, [0.0, 0.0, 0.0, 1.0]),  #1G, R to 2G
                    SparseCat(states_rock, [0.0, 1.0, 0.0, 0.0])); #1G, C to 1G
                    hcat(SparseCat(states_rock, [1.0, 0.0, 0.0, 0.0]),  #2B, L to 1B
                    SparseCat(states_rock, [0.0, 0.0, 1.0, 0.0]),  #2B, R to 2B
                    SparseCat(states_rock, [0.0, 0.0, 1.0, 0.0])); #2B, C to 2B
                    hcat(SparseCat(states_rock, [0.0, 1.0, 0.0, 0.0]),  #2G, L to 1G
                    SparseCat(states_rock, [0.0, 0.0, 0.0, 1.0]),  #2G, R to 2G
                    SparseCat(states_rock, [0.0, 0.0, 0.0, 1.0]))] #2G, C to 2G

# Nominal transition function array tarray_rock[s,a,s'] = Pr(s' | s ,a)
const tarray_rock = cat(3, [1.0 0.0 1.0; 0.0 0.0 0.0; 1.0 0.0 0.0; 0.0 0.0 0.0], # s' = 1B
                           [0.0 0.0 0.0; 1.0 0.0 1.0; 0.0 0.0 0.0; 1.0 0.0 0.0], # s' = 1G
                           [0.0 1.0 0.0; 0.0 0.0 0.0; 0.0 1.0 1.0; 0.0 0.0 0.0], # s' = 2B
                           [0.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0; 0.0 1.0 1.0]) # s' = 2G

# Nominal observation function distributions
function peff(xs::Int, ys::Int, xr::Int, yr::Int, maxeff::Float64)
   dist = sqrt((yr -ys)^2 + (xr - xs)^2)
   return maxeff * (1 + exp(-dist * 0.2)) / 2
end
peff(xs::Int, ys::Int, xr::Int, yr::Int) = peff(xs, ys, xr, yr, 1.0)
const ps1 = peff(1, 0, 2, 0, seff)
const ps2 = peff(2, 0, 2, 0, seff)

const odist_rock = [hcat(SparseCat(observations_rock, [0.0, 0.0, 1.0]),  # L, 1B
                   SparseCat(observations_rock, [0.0, 0.0, 1.0]),  # L, 1G
                   SparseCat(observations_rock, [0.0, 0.0, 1.0]),  # L, 2B
                   SparseCat(observations_rock, [0.0, 0.0, 1.0])); # L, 2G
                   hcat(SparseCat(observations_rock, [0.0, 0.0, 1.0]),  # R, 1B
                   SparseCat(observations_rock, [0.0, 0.0, 1.0]),  # R, 1G
                   SparseCat(observations_rock, [0.0, 0.0, 1.0]),  # R, 2B
                   SparseCat(observations_rock, [0.0, 0.0, 1.0])); # R, 2G
                   hcat(SparseCat(observations_rock, [1 - ps1, ps1, 0.0]),  # C, 1B
                   SparseCat(observations_rock, [ps1, 1 - ps1, 0.0]),  # C, 1G
                   SparseCat(observations_rock, [1 - ps2, ps2, 0.0]),  # C, 2B
                   SparseCat(observations_rock, [ps2, 1 - ps2, 0.0]))] # C, 2G

# Nominal observation function array oarray_rock[a,sp,z] = Pr(z|a,sp)
const oarray_rock = cat(3, [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 1-ps1 ps1 1-ps2 ps2], # z = good
                          [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; ps1 1-ps1 ps2 1-ps2], # z = bad
                          [1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0; 0.0 0.0 0.0 0.0])   # z = nothing

function buildor(array::Array{Float64}, ps1::Float64, ps2::Float64, uncsize::Float64)
  ps1_l = max(ps1 - 2 * pϵ_rock, 0.0 + pϵ_rock)
  ps1_u = min(ps1 + 2 * pϵ_rock, 1.0 - pϵ_rock)
  ps2_l = max(ps2 - uncsize, 0.0 + pϵ_rock)
  ps2_u = min(ps2 + uncsize, 1.0 - pϵ_rock)
  o_l = cat(3, [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 1-ps1_l ps1_l 1-ps2_l ps2_l],
               [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; ps1_l 1-ps1_l ps2_l 1-ps2_l],
               [1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0; 0.0 0.0 0.0 0.0])
  o_u = cat(3, [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 1-ps1_u ps1_u 1-ps2_u ps2_u],
               [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; ps1_u 1-ps1_u ps2_u 1-ps2_u],
               [1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0; 0.0 0.0 0.0 0.0])
  min.(o_l, o_u), max.(o_l, o_u)
end

function RockIPOMDP(alphas::Vector{Vector{Float64}}, disc::Float64)
    tdist = tdist_rock
    tarray = tarray_rock
    odist = odist_rock
    oarray = oarray_rock
    discount = disc
    RockIPOMDP(tdist, tarray, odist, oarray, alphas, discount)
end
RockIPOMDP(disc::Float64) = RockIPOMDP(inforeward_rock, disc)
RockIPOMDP(alphas::Vector{Vector{Float64}}) = RockIPOMDP(alphas, discount_rock)
RockIPOMDP() = RockIPOMDP(inforeward_rock)

function RockRIPOMDP(alphas::Vector{Vector{Float64}}, disc::Float64, usize::Float64)
    oarrayl, oarrayu = buildor(oarray_rock, ps1, ps2, usize)
    tarrayl = tarray_rock
    tarrayu = tarray_rock
    discount = disc
    RockRIPOMDP(tarrayl, tarrayu, oarrayl, oarrayu, alphas, discount, usize)
end
RockRIPOMDP(alphas::Vector{Vector{Float64}}) = RockRIPOMDP(alphas, discount_rock, 0.025)
RockRIPOMDP(err::Float64) = RockRIPOMDP(inforeward_rock, discount_rock, err)
RockRIPOMDP(alphas::Vector{Vector{Float64}}, disc::Float64) = RockRIPOMDP(alphas, disc, 0.025)
RockRIPOMDP(disc::Float64, err::Float64) = RockRIPOMDP(inforeward_rock, disc, err)
RockRIPOMDP() = RockRIPOMDP(0.025)

states(::Union{RockIPOMDP, RockRIPOMDP}) = states_rock
actions(::Union{RockIPOMDP, RockRIPOMDP}) = actions_rock
observations(::Union{RockIPOMDP, RockRIPOMDP}) = observations_rock

n_states(::Union{RockIPOMDP, RockRIPOMDP}) = 4
n_actions(::Union{RockIPOMDP, RockRIPOMDP}) = 3
n_observations(::Union{RockIPOMDP, RockRIPOMDP}) = 3

function state_index(::Union{RockIPOMDP, RockRIPOMDP}, s::Int)
    if s == 10
        return 1
    elseif s == 11
        return 2
    elseif s == 20
        return 3
    elseif s== 21
        return 4
    else
        return nothing
    end
end

function action_index(::Union{RockIPOMDP, RockRIPOMDP}, a::Symbol)
    if a == :left
        return 1
    elseif a == :right
        return 2
    elseif a == :check
        return 3
    else
        return nothing
    end
end

function observation_index(::Union{RockIPOMDP, RockRIPOMDP}, z::Symbol)
    if z == :good
        return 1
    elseif z == :bad
        return 2
    elseif z == :none
        return 3
    else
        return nothing
    end
end

obs_index(prob::Union{RockIPOMDP, RockRIPOMDP}, z::Symbol) = observation_index(prob, z)

initial_state_distribution(::Union{RockIPOMDP, RockRIPOMDP}) = initdist_rock
initial_belief(::Union{RockIPOMDP, RockRIPOMDP}) = initbelief_rock
initial_belief_distribution(::Union{RockIPOMDP, RockRIPOMDP}) = initbelief_dist_rock

# Nominal transitions
transition(prob::RockIPOMDP, s::Int, a::Symbol) = prob.tdist[state_index(prob, s), action_index(prob, a)]
transition(prob::RockIPOMDP, s::Int, a::Symbol, sp::Int) = prob.tarray[state_index(prob, s), action_index(prob, a), state_index(prob, sp)]

# Robust transitions
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PtLower = Plower(t|s,a) = t_lower[s,a,sp]
function transition(prob::RockRIPOMDP, s::Int, a::Symbol)
    si, ai = state_index(prob, s), action_index(prob, a)
    prob.tarrayl[si, ai, :], prob.tarrayu[si, ai, :]
end
function transition(prob::RockRIPOMDP, s::Int, a::Symbol, sp::Int)
    si, ai, spi = state_index(prob, s), action_index(prob, a), state_index(prob, sp)
    prob.tarrayl[si, ai, spi], prob.tarrayu[si, ai, spi]
end

# Nominal observations
observation(prob::RockIPOMDP, a::Symbol, sp::Int) = prob.odist[action_index(prob, a), state_index(prob, sp)]
observation(prob::RockIPOMDP, a::Symbol, sp::Int, z::Symbol) = prob.oarray[action_index(prob, a), state_index(prob, sp), observation_index(prob, z)]

# Robust observations
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PzLower = Plower(z|s,a,t) = o_lower[a,sp,z]
function observation(prob::RockRIPOMDP, a::Symbol, sp::Int)
    ai, spi = action_index(prob, a), state_index(prob, sp)
    prob.oarrayl[ai,spi,:], prob.oarrayu[ai,spi,:]
end
function observation(prob::RockRIPOMDP, a::Symbol, sp::Int, z::Symbol)
    ai, spi, zi = action_index(prob, a), state_index(prob, sp), observation_index(prob, z)
    prob.oarrayl[ai,spi,zi], prob.oarrayu[ai,spi,zi]
end

function dynamics(prob::RockIPOMDP)
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        pt = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pz = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        p[tind, zind, sind, aind] = pz * pt
    end
    p
end

function dynamics(prob::RockRIPOMDP)
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        ptl, ptu = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pzl, pzu = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        plower[tind, zind, sind, aind] = pzl * ptl
        pupper[tind, zind, sind, aind] = pzu * ptu
    end
    plower, pupper
end
function reward(prob::Union{RockIPOMDP, RockRIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    for α in prob.inforeward
        rmax = max(rmax, dot(α, b))
    end
    rmax
end
reward(prob::Union{RockIPOMDP,RockRIPOMDP}, b::Vector{Float64}, s::Int, a::Symbol, sp::Int) = reward(prob,b,a)

function rewardalpha(prob::Union{RockIPOMDP, RockRIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    ralpha = nothing
    for α in prob.inforeward
        rnext = dot(α, b)
        if rnext > rmax
            rmax = rnext
            ralpha = copy(α)
        end
    end
    ralpha
end

discount(p::Union{RockIPOMDP,RockRIPOMDP}) = p.discount

# Simulation Functions
function initial_state(prob::Union{RockIPOMDP,RockRIPOMDP}, rng::AbstractRNG)
    d = initial_state_distribution(prob)
    return rand(rng, d)
end

function generate_sor(prob::RockIPOMDP, b, s, a, rng::AbstractRNG)
    sp = rand(rng, transition(prob, s, a))
    o = rand(rng, observation(prob, a, sp))
    r = reward(prob, b, s, a, sp)
    sp, o, r
end

function generate_sor(prob::RockRIPOMDP, b, s, a, rng::AbstractRNG)
    tdist = SparseCat(states(prob), psample(transition(prob, s, a)...))
    sp = rand(rng, tdist)
    odist = SparseCat(observations(prob), psample(observation(prob, a, sp)...))
    o = rand(rng, odist)
    r = reward(prob, b, s, a, sp)
    sp, o, r
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using Documenter
using DocumenterMarkdown
using MXNet

"""
Return all files of a submodule

julia> listpages("ndarray")
15-element Array{String,1}:
 "ndarray.jl"
 "ndarray/activation.jl"
 "ndarray/arithmetic.jl"
 "ndarray/array.jl"
 ...
 "ndarray/statistic.jl"
 "ndarray/trig.jl"
 "ndarray/type.jl"
"""
listpages(x) =
  ["$x.jl"; joinpath.(x, readdir(joinpath(@__DIR__, "..", "src", x)))]

const api_pages = [
  "api/context.md",
  "api/ndarray.md",
  "api/symbolic-node.md",
  "api/model.md",
  "api/initializers.md",
  "api/optimizers.md",
  "api/callbacks.md",
  "api/metric.md",
  "api/io.md",
  "api/nn-factory.md",
  "api/executor.md",
  "api/kvstore.md",
  "api/visualize.md",
]

makedocs(
  sitename = "MXNet.jl",
  modules  = MXNet,
  doctest  = false,
  format   = Markdown(),
)

ext: jl


text:
export MIS2, MIS3C, MIS4, MIS4a, TVDMISA, TVDMISB

"""
    MultirateInfinitesimalStep

Multirate Infinitesimal Step (MIS) methods of [WKG2009](@cite) and [KW2014](@cite).

The available implementations are:

 - [`MIS2`](@ref)
 - [`MIS3C`](@ref)
 - [`MIS4`](@ref)
 - [`MIS4a`](@ref)
 - [`TVDMISA`](@ref)
 - [`TVDMISB`](@ref)
"""
abstract type MultirateInfinitesimalStep <: DistributedODEAlgorithm end

struct MultirateInfinitesimalStepTableau{Nstages, Nstages², RT}
  α::SArray{NTuple{2, Nstages}, RT, 2, Nstages²}
  β::SArray{NTuple{2, Nstages}, RT, 2, Nstages²}
  γ::SArray{NTuple{2, Nstages}, RT, 2, Nstages²}
  d::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
  c::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
  c̃::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
end
function MultirateInfinitesimalStepTableau(α,β,γ)
  d = SVector(sum(β, dims = 2)) # KW2014 (2)
  # c = (I - α - γ) \ d
  c = similar(d)
  for i in eachindex(c)
      c[i] = d[i]
      for j = 1:i-1
        c[i] += (α[i, j] + γ[i, j]) * c[j]
      end
  end
  c = SArray(c)
  c̃ = α * c
  MultirateInfinitesimalStepTableau(α,β,γ,d,c,c̃)
end

struct MultirateInfinitesimalStepCache{Nstages, A, T<:MultirateInfinitesimalStepTableau}
  "difference between stage and initial value ``U^{(i)} - u``"
  ΔU::NTuple{Nstages,A}
  "evaluated slow part of each stage ``f_slow(U^{(i)})``"
  F::NTuple{Nstages,A}
  tableau::T
end

nstages(::MultirateInfinitesimalStepCache{Nstages}) where {Nstages} = Nstages

function cache(
  prob::DiffEqBase.AbstractODEProblem{uType, tType, true},
  alg::MultirateInfinitesimalStep; kwargs...) where {uType,tType}

  tab = tableau(alg,eltype(prob.u0))
  N = length(tab.d)
  ΔU = ntuple(n -> similar(prob.u0), N)
  # at time i, contains
  #   ΔU[j] = U[j] - u    j < i
  #   ΔU[i] = U[i]
  #   ΔU[N] = offset vec
  F = ntuple(n -> similar(prob.u0), N)
  return MultirateInfinitesimalStepCache(ΔU, F, tab)
end

function init_inner(prob, outercache::MultirateInfinitesimalStepCache, dt)
  OffsetODEFunction(prob.f.f1, zero(dt), one(dt), one(dt), outercache.ΔU[end])
end

function update_inner!(innerinteg, outercache::MultirateInfinitesimalStepCache,
  f_slow, u, p, t, dt, i)

  f_offset = innerinteg.prob.f
  tab = outercache.tableau
  N = nstages(outercache)

  F = outercache.F
  ΔU = outercache.ΔU

  # F[i] = f_slow(U[i-1], p, t + c[i-1]*dt)
  u0 = i == 1 ? u : ΔU[i-1]
  t0 = i == 1 ? t : t + tab.c[i-1]*dt
  f_slow(F[i], u0, p, t0)

  # the (i+1)th stage of the paper
  innerinteg.u = i == N ? u : ΔU[i]

  groupsize = 256
  event = Event(array_device(u))
  event = mis_update!(array_device(u), groupsize)(
    u,
    ΔU,
    F,
    innerinteg.u,
    f_offset.x,
    tab,
    i,
    N,
    dt;
    ndrange = length(u),
    dependencies = (event,),
  )
  wait(array_device(u), event)

  # KW2014 (9)
  # evaluate f_fast(z(τ), p, t + c̃[i]*dt + (c[i]-c̃[i])/d[i] * τ)
  f_offset.α = t + tab.c̃[i] * dt
  f_offset.β = (tab.c[i] - tab.c̃[i]) / tab.d[i]

  innerinteg.t = zero(t)
  innerinteg.tstop = tab.d[i] * dt
end

@kernel function mis_update!(
  u,
  ΔU,
  F,
  innerinteg_u,
  f_offset_x,
  tab,
  i,
  N,
  dt,
)
  e = @index(Global, Linear)
  @inbounds begin
    if i > 1
      ΔU[i-1][e] -= u[e]
    end

    # KW2014 (1a)
    if i < N
      innerinteg_u[e] = u[e]
    end
    for j = 1:i-1
      innerinteg_u[e] += tab.α[i,j] * ΔU[j][e]
    end

    # KW2014 (1b) / (9)
    f_offset_x[e] = tab.β[i,i]/tab.d[i] .* F[i][e]
    for j = 1:i-1
      f_offset_x[e] += (tab.γ[i,j]/(tab.d[i]*dt)) * ΔU[j][e] + tab.β[i,j]/tab.d[i]  * F[j][e]
    end
  end
end

"""
    MIS2()

The MIS2 Multirate Infinitesimal Step (MIS) method of [KW2014](@cite).
"""
struct MIS2 <: MultirateInfinitesimalStep
end
function tableau(mis::MIS2, RT)
  α = @SArray RT[
    0 0 0
    0.536946566710 0 0
    0.480892968551 0.500561163566 0
  ]
  β = @SArray RT[
    0.126848494553 0 0
    -0.784838278826 1.37442675268 0
    -0.0456727081749 -0.00875082271190 0.524775788629
  ]
  γ = @SArray RT[
    0 0 0
    0.652465126004 0 0
    -0.0732769849457 0.144902430420 0
  ]
  MultirateInfinitesimalStepTableau(α,β,γ)
end

"""
    MIS3C()

The MIS3C Multirate Infinitesimal Step (MIS) method of [KW2014](@cite).
"""
struct MIS3C <: MultirateInfinitesimalStep
end
function tableau(mis::MIS3C, RT)
  α = @SArray RT[
    0 0 0
    0.589557277145 0 0
    0.544036601551 0.565511042564 0
  ]
  β = @SArray RT[
    0.397525189225 0 0
    -0.227036463644 0.624528794618 0
    -0.00295238076840 -0.270971764284 0.671323159437
  ]
  γ = @SArray RT[
    0 0 0
    0.142798786398 0 0
    -0.0428918957402 0.0202720980282 0
  ]
  MultirateInfinitesimalStepTableau(α,β,γ)
end

"""
    MIS4()

The MIS4 Multirate Infinitesimal Step (MIS) method of [KW2014](@cite).
"""
struct MIS4 <: MultirateInfinitesimalStep
end
function tableau(mis::MIS4, RT)
  α = @SArray RT[
    0 0 0 0
    0.914092810304 0 0 0
    1.14274417397 -0.295211246188 0 0
    0.112965282231 0.337369411296 0.503747183119 0
  ]
  β = @SArray RT[
    0.136296478423 0 0 0
    0.280462398979 -0.0160351333596 0 0
    0.904713355208 -1.04011183154 0.652337563489 0
    0.0671969845546 -0.365621862610 -0.154861470835 0.970362444469
  ]
  γ = @SArray RT[
    0 0 0 0
    0.678951983291 0 0 0
    -1.38974164070 0.503864576302 0 0
    -0.375328608282 0.320925021109 -0.158259688945 0
  ]
  MultirateInfinitesimalStepTableau(α,β,γ)
end

"""
    MIS4a()

The MIS4a Multirate Infinitesimal Step (MIS) method of [KW2014](@cite).
"""
struct MIS4a <: MultirateInfinitesimalStep
end
function tableau(mis::MIS4a, RT)
  α = @SArray RT[
    0 0 0 0
    0.52349249922385610 0 0 0
    1.1683374366893629 -0.75762080241712637 0 0
    -0.036477233846797109 0.56936148730740477 0.47746263002599681 0
  ]
  # β₅₁ in the paper is incorrect
  # the correct value (β[4,1]) is used below (from authors)
  β = @SArray RT[
    0.38758444641450318 0 0 0
    -0.025318448354142823 0.38668943087310403 0 0
    0.20899983523553325 -0.45856648476371231 0.43423187573425748 0
    -0.10048822195663100 -0.46186171956333327 0.83045062122462809 0.27014914900250392
  ]
  γ = @SArray RT[
    0 0 0 0
    0.13145089796226542 0 0 0
    -0.36855857648747881 0.33159232636600550 0 0
    -0.065767130537473045 0.040591093109036858 0.064902111640806712 0
  ]
  MultirateInfinitesimalStepTableau(α,β,γ)
end

"""
    TVDMISA()

The TVDMISA Total Variation Diminishing (TVD) Multirate Infinitesimal Step (MIS)
method of [KW2014](@cite).
"""
struct TVDMISA <: MultirateInfinitesimalStep
end
function tableau(mis::TVDMISA, RT)
  α = @SArray RT[
    0 0 0
    0.1946360605647457 0 0
    0.3971200136786614 0.2609434606211801 0
  ]
  β = @SArray RT[
    2//3 0 0
    -0.28247174703488398 4//9 0
    -0.31198081960042401 0.18082737579913699 9//16
  ]
  γ = @SArray RT[
    0 0 0
    0.5624048933209129 0 0
    0.4408467475713277 -0.2459300561692391 0
  ]
  MultirateInfinitesimalStepTableau(α,β,γ)
end

"""
    TVDMISB()

The TVDMISB Total Variation Diminishing (TVD) Multirate Infinitesimal Step (MIS)
method of [KW2014](@cite).
"""
struct TVDMISB <: MultirateInfinitesimalStep
end
function tableau(mis::TVDMISB, RT)
  α = @SArray RT[
    0 0 0
    0.42668232863311001 0 0
    0.26570779016173801 0.41489966891866698 0
  ]
  β = @SArray RT[
    2//3 0 0
    -0.25492859100078202 4//9 0
    -0.26452517179288798 0.11424084424766399 9//16
  ]
  γ = @SArray RT[
    0 0 0
    0.28904389120139701 0 0
    0.45113560071334202 -0.25006656847591002 0
  ]
  MultirateInfinitesimalStepTableau(α,β,γ)
end


text:
abstract type JLPSexp end

abstract type JLPType <: JLPSexp end
abstract type JLPAbstList <: JLPSexp end

abstract type JLPPrimType <: JLPType end

abstract type JLPNumber <: JLPPrimType end
abstract type JLPAbsChar <: JLPPrimType end

struct JLPSymbol <: JLPType
    val::String
end

struct JLPFloat <: JLPNumber
    val::Float32
end

struct JLPInteger <: JLPNumber
    val::Int32
end

struct JLPBool <: JLPNumber
    val::Bool
end

struct JLPChar <: JLPPrimType
    val::Char
end

struct JLPString <: JLPPrimType
    val::JLPSymbol
    str::String
end

struct JLPNil <: JLPAbstList
end

const jlpnil = JLPNil()

Base.iterate(iter::JLPNil) = nothing

struct JLPList <: JLPAbstList
    val::Union{JLPSymbol, JLPNil}
    car::JLPSexp
    cdr::JLPAbstList

    JLPList(a ,b ,c) = new(a, b ,c)
    JLPList(a ,b) = new(JLPSymbol(""), a, b)
end

Base.iterate(iter::JLPList) = (iter.car, iter.cdr)

function Base.iterate(iter::JLPList, state::Union{JLPList, JLPNil})
    if typeof(state) <: JLPNil
        nothing
    else
        (state.car, state.cdr)
    end
end

function add_list(l::JLPAbstList, target::JLPAbstList)
    if l == jlpnil
        return target
    elseif l.cdr == jlpnil
        return JLPList(l.val, l.car, target)
    else
        return JLPList(l.val, l.car, add_list(l.cdr, target))
    end
end


text:
export BrickletRealTimeClockDateTime
struct BrickletRealTimeClockDateTime
    year::Integer
    month::Integer
    day::Integer
    hour::Integer
    minute::Integer
    second::Integer
    centisecond::Integer
    weekday::Integer
end

export BrickletRealTimeClockAlarm
struct BrickletRealTimeClockAlarm
    month::Integer
    day::Integer
    hour::Integer
    minute::Integer
    second::Integer
    weekday::Integer
    interval::Integer
end

export BrickletRealTimeClockIdentity
struct BrickletRealTimeClockIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletRealTimeClock
"""
Battery-backed real-time clock
"""
mutable struct BrickletRealTimeClock <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletRealTimeClock(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_real_time_clock")
        deviceInternal = package.BrickletRealTimeClock(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_date_time
"""
    $(SIGNATURES)

Sets the current date (including weekday) and the current time.

If the backup battery is installed then the real-time clock keeps date and
time even if the Bricklet is not powered by a Brick.

The real-time clock handles leap year and inserts the 29th of February
accordingly. But leap seconds, time zones and daylight saving time are not
handled.
"""
function set_date_time(device::BrickletRealTimeClock, year, month, day, hour, minute, second, centisecond, weekday)
    device.deviceInternal.set_date_time(year, month, day, hour, minute, second, centisecond, weekday)
end

export get_date_time
"""
    $(SIGNATURES)

Returns the current date (including weekday) and the current time of the
real-time clock.
"""
function get_date_time(device::BrickletRealTimeClock)
    return device.deviceInternal.get_date_time()
end

export get_timestamp
"""
    $(SIGNATURES)

Returns the current date and the time of the real-time clock.
The timestamp has an effective resolution of hundredths of a
second and is an offset to 2000-01-01 00:00:00.000.
"""
function get_timestamp(device::BrickletRealTimeClock)
    return device.deviceInternal.get_timestamp()
end

export set_offset
"""
    $(SIGNATURES)

Sets the offset the real-time clock should compensate for in 2.17 ppm steps
between -277.76 ppm (-128) and +275.59 ppm (127).

The real-time clock time can deviate from the actual time due to the frequency
deviation of its 32.768 kHz crystal. Even without compensation (factory
default) the resulting time deviation should be at most ±20 ppm (±52.6
seconds per month).

This deviation can be calculated by comparing the same duration measured by the
real-time clock (``rtc_duration``) an accurate reference clock
(``ref_duration``).

For best results the configured offset should be set to 0 ppm first and then a
duration of at least 6 hours should be measured.

The new offset (``new_offset``) can be calculated from the currently configured
offset (``current_offset``) and the measured durations as follow::

  new_offset = current_offset - round(1000000 * (rtc_duration - ref_duration) / rtc_duration / 2.17)

If you want to calculate the offset, then we recommend using the calibration
dialog in Brick Viewer, instead of doing it manually.

The offset is saved in the EEPROM of the Bricklet and only needs to be
configured once.
"""
function set_offset(device::BrickletRealTimeClock, offset)
    device.deviceInternal.set_offset(offset)
end

export get_offset
"""
    $(SIGNATURES)

Returns the offset as set by :func:`Set Offset`.
"""
function get_offset(device::BrickletRealTimeClock)
    return device.deviceInternal.get_offset()
end

export set_date_time_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Date Time` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Date Time` Callback is only triggered if the date or time changed
since the last triggering.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_date_time_callback_period(device::BrickletRealTimeClock, period)
    device.deviceInternal.set_date_time_callback_period(period)
end

export get_date_time_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Date Time Callback Period`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_date_time_callback_period(device::BrickletRealTimeClock)
    return device.deviceInternal.get_date_time_callback_period()
end

export set_alarm
"""
    $(SIGNATURES)

Configures a repeatable alarm. The :cb:`Alarm` callback is triggered if the
current date and time matches the configured alarm.

Setting a parameter to -1 means that it should be disabled and doesn't take part
in the match. Setting all parameters to -1 disables the alarm completely.

For example, to make the alarm trigger every day at 7:30 AM it can be
configured as (-1, -1, 7, 30, -1, -1, -1). The hour is set to match 7 and the
minute is set to match 30. The alarm is triggered if all enabled parameters
match.

The interval has a special role. It allows to make the alarm reconfigure itself.
This is useful if you need a repeated alarm that cannot be expressed by matching
the current date and time. For example, to make the alarm trigger every 23
seconds it can be configured as (-1, -1, -1, -1, -1, -1, 23). Internally the
Bricklet will take the current date and time, add 23 seconds to it and set the
result as its alarm. The first alarm will be triggered 23 seconds after the
call. Because the interval is not -1, the Bricklet will do the same again
internally, take the current date and time, add 23 seconds to it and set that
as its alarm. This results in a repeated alarm that triggers every 23 seconds.

The interval can also be used in combination with the other parameters. For
example, configuring the alarm as (-1, -1, 7, 30, -1, -1, 300) results in an
alarm that triggers every day at 7:30 AM and is then repeated every 5 minutes.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_alarm(device::BrickletRealTimeClock, month, day, hour, minute, second, weekday, interval)
    device.deviceInternal.set_alarm(month, day, hour, minute, second, weekday, interval)
end

export get_alarm
"""
    $(SIGNATURES)

Returns the alarm configuration as set by :func:`Set Alarm`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_alarm(device::BrickletRealTimeClock)
    return device.deviceInternal.get_alarm()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletRealTimeClock)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletRealTimeClock, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using EDFPlus
using Test
# Run tests

@test include("readtest.jl")
@test include("writetest.jl")
@test include("cover.jl")


text:
module TermInterface

"""
    istree(x)

Returns `true` if `x` is a term. If true, `operation`, `arguments`
must also be defined for `x` appropriately.
"""
istree(x) = istree(typeof(x))
istree(x::Type{T}) where {T} = false
export istree

"""
    symtype(x)

Returns the symbolic type of `x`. By default this is just `typeof(x)`.
Define this for your symbolic types if you want `SymbolicUtils.simplify` to apply rules
specific to numbers (such as commutativity of multiplication). Or such
rules that may be implemented in the future.
"""
function symtype(x)
    typeof(x)
end
export symtype

"""
    issym(x)

Returns `true` if `x` is a symbol. If true, `nameof` must be defined
on `x` and must return a Symbol.
"""
issym(x) = issym(typeof(x))
issym(x::Type{T}) where {T} = false
export issym

"""
    exprhead(x)

If `x` is a term as defined by `istree(x)`, `exprhead(x)` must return a symbol,
corresponding to the head of the `Expr` most similar to the term `x`.
If `x` represents a function call, for example, the `exprhead` is `:call`.
If `x` represents an indexing operation, such as `arr[i]`, then `exprhead` is `:ref`.
Note that `exprhead` is different from `operation` and both functions should
be defined correctly in order to let other packages provide code generation
and pattern matching features.
"""
function exprhead end
export exprhead

"""
    operation(x)

If `x` is a term as defined by `istree(x)`, `operation(x)` returns the
head of the term if `x` represents a function call, for example, the head
is the function being called.
"""
function operation end
export operation

"""
    arguments(x)

Get the arguments of `x`, must be defined if `istree(x)` is `true`.
"""
function arguments end
export arguments

"""
    unsorted_arguments(x::T)

If x is a term satisfying `istree(x)` and your term type `T` orovides
and optimized implementation for storing the arguments, this function can
be used to retrieve the arguments when the order of arguments does not matter
but the speed of the operation does.
"""
unsorted_arguments(x) = arguments(x)
export unsorted_arguments

"""
    arity(x)

Returns the number of arguments of `x`. Implicitly defined
if `arguments(x)` is defined.
"""
arity(x) = length(arguments(x))
export arity

"""
    metadata(x)

Return the metadata attached to `x`.
"""
metadata(x) = nothing
export metadata

"""
    metadata(x, md)

Returns a new term which has the structure of `x` but also has
the metadata `md` attached to it.
"""
function metadata(x, data)
    error("Setting metadata on $x is not possible")
end

"""
    similarterm(x, head, args, symtype=nothing; metadata=nothing, exprhead=:call)

Returns a term that is in the same closure of types as `typeof(x)`,
with `head` as the head and `args` as the arguments, `type` as the symtype
and `metadata` as the metadata. By default this will execute `head(args...)`.
`x` parameter can also be a `Type`. The `exprhead` keyword argument is useful
when manipulating `Expr`s.
"""
function similarterm(x, head, args, symtype=nothing; metadata=nothing, exprhead=nothing)
    if exprhead === nothing
        similarterm(typeof(x), head, args, symtype; metadata=metadata)
    else
        similarterm(typeof(x), head, args, symtype; metadata=metadata, exprhead=exprhead)
    end
end

function similarterm(x::Type{T}, head, args, symtype=nothing; metadata=nothing, exprhead=:call) where T
    !istree(T) ? head : head(args...)
end
export similarterm

include("utils.jl")

include("expr.jl")

end # module


text:
using Parameters, JuMP

"""Define energy system model type as JuMP.Model."""
const EnergySystemModel = Model

"""Specifation for which constraints to include to the model. Constraints that
are not specified are included by default.
# Arguments
- `renewable_target::Bool`: Whether to include renewables target constraint.
- `storage::Bool`: Whether to include storage constraints.
- `ramping::Bool`: Whether to include ramping constraints.
- `voltage_angles::Bool`: Whether to include voltage angle constraints.
"""
@with_kw struct Specs
    transmission::Bool = true
    renewable_target::Bool = false
    carbon_cap::Bool = false
    nuclear_limit::Bool = false
    storage::Bool = false
    ramping::Bool = false
    voltage_angles::Bool = false
    hydro::Bool = false
    hydro_simple::Bool = false
end

"""Input indices and parameters for the model."""
@with_kw struct Params
    region_n::Array{AbstractString, 1}
    max_dem_n::Array{Float64, 1}
    technology_g::Array{AbstractString, 1}
    G::Array{Integer, 1}
    G_r::Array{Integer, 1}
    N::Array{Integer, 1}
    L::Array{Array{Integer, 1}, 1}
    L_ind::Array{Integer, 1}
    T::Array{Integer, 1}
    S::Array{Integer, 1}
    H::Array{Integer, 1}
    κ::AbstractFloat
    μ::AbstractFloat
    C::AbstractFloat
    C̄::AbstractFloat
    C_E::AbstractFloat
    R_E::AbstractFloat
    τ_t::Array{Integer, 1}
    Gmin_gn::Array{AbstractFloat, 2}
    Gmax_gn::Array{AbstractFloat, 2}
    A_gnt::Array{AbstractFloat, 3}
    D_nt::Array{AbstractFloat, 2}
    I_g::Array{AbstractFloat, 1}
    M_g::Array{AbstractFloat, 1}
    C_g::Array{AbstractFloat, 1}
    e_g::Array{AbstractFloat, 1}
    E_g::Array{AbstractFloat, 1}
    r⁻_g::Array{AbstractFloat, 1}
    r⁺_g::Array{AbstractFloat, 1}
    I_l::Array{AbstractFloat, 1}
    M_l::Array{AbstractFloat, 1}
    C_l::Array{AbstractFloat, 1}
    B_l::Array{AbstractFloat, 1}
    e_l::Array{AbstractFloat, 1}
    Tmin_l::Array{AbstractFloat, 1}
    Tmax_l::Array{AbstractFloat, 1}
    ξ_s::Array{AbstractFloat, 1}
    I_s::Array{AbstractFloat, 1}
    C_s::Array{AbstractFloat, 1}
    Smin_sn::Array{AbstractFloat, 2}
    Smax_sn::Array{AbstractFloat, 2}
    Wmax_hn::Array{AbstractFloat, 2}
    Wmin_hn::Array{AbstractFloat, 2}
    Hmax_hn::Array{AbstractFloat, 2}
    Hmin_hn::Array{AbstractFloat, 2}
    HRmax_n::Array{AbstractFloat, 1}
    Fmin_n::Array{AbstractFloat, 1}
    AH_nt::Array{AbstractFloat, 2}
    AR_nt::Array{AbstractFloat, 2}
    I_h::Vector{Float64}
    M_h::Vector{Float64}
    C_h::Vector{Float64}
    e_h::Vector{Float64}
    E_h::Vector{Float64}
    r⁻_h::Vector{Float64}
    r⁺_h::Vector{Float64}
end

"""Retrieving data from objects typed JuMP.Containers.DenseAxisArray"""
retrieve_data(a::Number) = a
retrieve_data(a::JuMP.Containers.DenseAxisArray) = (JuMP.value.(a)).data
retrieve_data(a::AffExpr) = JuMP.value.(a)
retrieve_data(a::Vector{AffExpr}) = JuMP.value.(a);

"""Extract objective values from a JuMP model.
# Arguments
- `model::EnergySystemModel`
- `JuMPObjDict::::Union{Dict{String, Float64}, Dict{String, Any}}`
"""
function JuMPObj(model::EnergySystemModel, JuMPObjDict::Dict{String, Any})
    dict = Dict(i => model[Symbol(i)] |> retrieve_data |> first for i in keys(JuMPObjDict))
    return dict
end

"""Extract objective values from a JuMP model.
# Arguments
- `model::EnergySystemModel`
- `JuMPObjDict::::Union{Dict{String, Float64}, Dict{String, Any}}`
"""
function JuMPVar(model::EnergySystemModel, JuMPVarDict::Dict{String, Any})
    dict = Dict(i => model[Symbol(i)] |> retrieve_data for i in keys(JuMPVarDict))
    return dict
end

"""Compute expression values from the results.
# Arguments
- `parameters::Params`
- `variables::Variables`
"""
function Expressions(parameters::Params, specs::Specs, variables::Dict{String, Array{Float64}})
    @unpack G, G_r, N, T, H, R_E, e_g, E_g, τ_t = parameters

    """Expression values:
    # Attributes
    - `κ′: Renewable generation share
    - `μ′: Hydro share
    - `C′_E: CO2 emission reduction
    """
    p_gnt = variables["p_gnt"]

    if (specs.hydro || specs.hydro_simple)
        h_hnt = variables["h_hnt"]
        hr_nt = variables["hr_nt"]

        κ′ = (sum(p_gnt[g,n,t]*τ_t[t] for g in G_r, n in N, t in T) + sum(h_hnt[h,n,t]*τ_t[t] for h in H, n in N, t in T) + sum(hr_nt[n,t]*τ_t[t] for n in N, t in T)) /
            (sum(p_gnt[g,n,t]*τ_t[t] for g in G, n in N, t in T) + sum(h_hnt[h,n,t]*τ_t[t] for h in H, n in N, t in T) + sum(hr_nt[n,t]*τ_t[t] for n in N, t in T))
        μ′ = sum(p_gnt[5,n,t]*τ_t[t] for n in N, t in T) /
            (sum(p_gnt[g,n,t]*τ_t[t] for g in G, n in N, t in T) + sum(h_hnt[h,n,t]*τ_t[t] for h in H, n in N, t in T) + sum(hr_nt[n,t]*τ_t[t] for n in N, t in T))
    else
        κ′ = (sum(p_gnt[g,n,t]*τ_t[t] for g in G_r, n in N, t in T)) /
            (sum(p_gnt[g,n,t]*τ_t[t] for g in G, n in N, t in T))
        μ′ = sum(p_gnt[5,n,t]*τ_t[t] for n in N, t in T) /
            (sum(p_gnt[g,n,t]*τ_t[t] for g in G, n in N, t in T))
    end

    C′_E = 1 - (sum(E_g[g] * (sum(p_gnt[g,n,t]*τ_t[t] for n in N, t in T)) / e_g[g] for g in G)) / R_E

    ExpressionsDict = Dict("κ′" => κ′, "μ′" => μ′, "C′_E" => C′_E)

    return ExpressionsDict
end

"""Creates the energy system model.
# Arguments
- `parameters::Params`
- `specs::Specs`
"""
function EnergySystemModel(parameters::Params, specs::Specs)
    @unpack max_dem_n, G, G_r, N, L, L_ind, T, S, H, κ, μ, C, C̄, C_E, R_E, τ_t, Gmin_gn, Gmax_gn, A_gnt, D_nt, I_g, M_g,
            C_g, e_g, E_g, r⁻_g, r⁺_g, I_l, M_l, C_l, B_l, e_l, Tmin_l, Tmax_l, ξ_s, I_s, C_s, Smin_sn, Smax_sn,
            Wmax_hn, Wmin_hn, Hmax_hn, Hmin_hn, HRmax_n, Fmin_n, AH_nt, AR_nt,
            I_h, M_h, r⁻_h, r⁺_h =
            parameters

    # TODO: include hydro environment constraints: i.e., use C_h, e_h, E_h

    ## Assumptions and caveats
    # i) Numerical scaling: divisions made through the model (usually per 1000) are done to improve numerical scale

    """Variable values."""
    VariablesDict = Dict{String, Any}()

    """Objective values:
    # Attributes
    - `f1: Generation investment and maintenance
    - `f2: Generation operational cost
    - `f3: Shedding cost
    - `f4: Transmission investment and maintenance
    - `f5: Transmission operational cost
    - `f6: Storage investment cost
    - `f7: Storage operational cost
    - `f8: Hydro investment
    """
    ObjectivesDict = Dict{String, Any}()

    # Indices of lines L: L_ind

    # Create an instance of JuMP model.
    model = EnergySystemModel()

    # -- Main variables --
    @variable(model, p_gnt[g in G, n in N, t in T] ≥ 0)
    VariablesDict["p_gnt"] = p_gnt
    @variable(model, Gmax_gn[g,n] ≥ p̄_gn[g in G, n in N] ≥ Gmin_gn[g,n])
    VariablesDict["p̄_gn"] = p̄_gn
    @variable(model, σ_nt[n in N, t in T] ≥ 0)
    VariablesDict["σ_nt"] = σ_nt

    # Transmission variables
    if specs.transmission
        @variable(model, f_lt[l in L_ind, t in T])
        VariablesDict["f_lt"] = f_lt
        @variable(model, f_abs_lt[l in L_ind, t in T] ≥ 0)
        VariablesDict["f_abs_lt"] = f_abs_lt
        @variable(model, Tmax_l[l] ≥ f̄_l[l in L_ind] ≥ Tmin_l[l])
        VariablesDict["f̄_l"] = f̄_l
    end

    if specs.storage
        # Storage variables
        @variable(model, b_snt[s in S, n in N, t in T] ≥ 0)
        VariablesDict["b_snt"] = b_snt
        @variable(model, Smax_sn[s,n] ≥ b̄_sn[s in S, n in N] ≥ Smin_sn[s,n])
        VariablesDict["b̄_sn"] = b̄_sn
        @variable(model, b⁺_snt[s in S, n in N, t in T] ≥ 0)
        VariablesDict["b⁺_snt"] = b⁺_snt
        @variable(model, b⁻_snt[s in S, n in N, t in T] ≥ 0)
        VariablesDict["b⁻_snt"] = b⁻_snt
    end

    if specs.voltage_angles
        # Voltage variables
        @variable(model, θ_nt[n in N, t in T] ≥ 0)
        VariablesDict["θ_nt"] = θ_nt
        @variable(model, θ′_nt[n in N, t in T] ≥ 0)
        VariablesDict["θ′_nt"] = θ′_nt
    end

    if specs.hydro
        # Hydro energy variables
        @variable(model, w_hnt[h in H, n in N, t in T]  ≥ Wmin_hn[h,n])
        VariablesDict["w_hnt"] = w_hnt
        @variable(model, h_hnt[h in H, n in N, t in T] ≥ 0)
        VariablesDict["h_hnt"] = h_hnt
        @variable(model, hr_nt[n in N, t in T] ≥ 0)
        VariablesDict["hr_nt"] = hr_nt
        @variable(model, h̄_hn[h in H, n in N] ≥ Hmin_hn[h,n])
        VariablesDict["h̄_hn"] = h̄_hn
    end

    if specs.hydro_simple
        if !(specs.hydro)
            # Hydro energy variables
            @variable(model, h_hnt[h in H, n in N, t in T] ≥ 0)
            VariablesDict["h_hnt"] = h_hnt
            @variable(model, h̄_hn[h in H, n in N] ≥ Hmin_hn[h,n])
            VariablesDict["h̄_hn"] = h̄_hn
            @variable(model, 0 ≤ hr_nt[n in N, t in T] ≤ 0)
            VariablesDict["hr_nt"] = hr_nt
        end
        # Compute maximal levels
        a_n = zeros(length(N))
        for n in N
            a_n[n] = sum(AH_nt[n,:])/(sum(Hmin_hn[:,n]) * length(T))
            isnan(a_n[n]) ? a_n[n] = 0 : a_n[n] = a_n[n]
        end
    end

    ## -- Objective --
    # Investment and maintenance of generation capacity
    @expression(model, f1,
        sum(I_g[g]*(p̄_gn[g,n]-Gmin_gn[g,n]) + M_g[g]*p̄_gn[g,n] for g in G, n in N))
    ObjectivesDict["f1"] = f1

    # Operational cost of generation dispatch
    @expression(model, f2,
        sum(C_g[g]*p_gnt[g,n,t]*τ_t[t] for g in G, n in N, t in T))
    ObjectivesDict["f2"] = f2

    # Shedding cost
    @expression(model, f3,
        sum(C*σ_nt[n,t]*τ_t[t] for n in N, t in T))
    ObjectivesDict["f3"] = f3

    if specs.transmission
        # Investment and maintenance cost of transmission cpacity
        @expression(model, f4,
            sum(I_l[l]*(f̄_l[l]-Tmin_l[l]) + M_l[l]*f̄_l[l] for l in L_ind))
        ObjectivesDict["f4"] = f4

        # Operational cost of transmission flow
        @expression(model, f5,
            sum(C_l[l]*f_abs_lt[l,t]*τ_t[t] for l in L_ind, t in T))
        ObjectivesDict["f5"] = f5
    end

    if specs.storage
        # Investment cost of storage capacity
        @expression(model, f6,
            sum(I_s[s]*(b̄_sn[s,n] - Smin_sn[s,n]) for s in S, n in N))
        ObjectivesDict["f6"] = f6

        # Operational cost of storage
        @expression(model, f7,
            sum(C_s[s]*(b⁺_snt[s,n,t] + b⁻_snt[s,n,t])*τ_t[t] for s in S, n in N, t in T))
        ObjectivesDict["f7"] = f7
    end

    # Investment cost of hydro capacity
    if (specs.hydro || specs.hydro_simple)
        @expression(model, f8,
            sum(I_h*(h̄_hn[h,n] - Hmin_hn[h,n]) + M_h*h̄_hn[h,n] for h in H, n in N))
        ObjectivesDict["f8"] = f8

        # TODO: add operational costs for hydro generation (i.e., C_h)
    end

    @objective(model, Min, sum(sum(flatten(ObjectivesDict[i])) for i in keys(ObjectivesDict))/10^6)

    ## -- Constraints --
    # Transmission lines to node n
    L⁻(n) = (l for (l,(i,j)) in zip(L_ind,L) if j==n)
    # Transmission lines from node n
    L⁺(n) = (l for (l,(i,j)) in zip(L_ind,L) if i==n)

    # TODO: add hydro generation efficiency dependencies (i.e., e_h)

    # Energy balance (dependent on the features selected)
    if specs.transmission && specs.storage && (specs.hydro || specs.hydro_simple)                   # Trans/Stor/Hydro
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t] +
        sum(e_l[l]*f_lt[l,t] for l in L⁻(n)) - sum(e_l[l]*f_lt[l,t] for l in L⁺(n)) +
        sum(ξ_s[s]*b⁻_snt[s,n,t] - b⁺_snt[s,n,t] for s in S) +
        sum(h_hnt[h,n,t] for h in H) + hr_nt[n,t])/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    elseif specs.transmission && specs.storage && !(specs.hydro || specs.hydro_simple)            # Trans/Stor
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t] +
        sum(e_l[l]*f_lt[l,t] for l in L⁻(n)) - sum(e_l[l]*f_lt[l,t] for l in L⁺(n)) +
        sum(ξ_s[s]*b⁻_snt[s,n,t] - b⁺_snt[s,n,t] for s in S))/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    elseif specs.transmission && !(specs.storage) && (specs.hydro || specs.hydro_simple)            # Trans/Hydro
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t] +
        sum(e_l[l]*f_lt[l,t] for l in L⁻(n)) - sum(e_l[l]*f_lt[l,t] for l in L⁺(n)) +
        sum(h_hnt[h,n,t] for h in H) + hr_nt[n,t] )/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    elseif specs.transmission && !(specs.storage) && !(specs.hydro || specs.hydro_simple)         # Trans
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t] +
        sum(e_l[l]*f_lt[l,t] for l in L⁻(n)) - sum(e_l[l]*f_lt[l,t] for l in L⁺(n)))/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    elseif !(specs.transmission) && !(specs.storage) && !(specs.hydro || specs.hydro_simple)      # -
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t])/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    elseif !(specs.transmission) && specs.storage && !(specs.hydro || specs.hydro_simple)         # Stor
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t] +
        sum(ξ_s[s]*b⁻_snt[s,n,t] - b⁺_snt[s,n,t] for s in S))/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    elseif !(specs.transmission) && !(specs.storage) && (specs.hydro || specs.hydro_simple)         # Hydro
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t] +
        sum(h_hnt[h,n,t] for h in H) + hr_nt[n,t] )/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    elseif !(specs.transmission) && specs.storage && (specs.hydro || specs.hydro_simple)            # Stor/Hydro
        @constraint(model,
        b1[n in N, t in T],
        (sum(p_gnt[g,n,t] for g in G) + σ_nt[n,t] +
        sum(ξ_s[s]*b⁻_snt[s,n,t] - b⁺_snt[s,n,t] for s in S) +
        sum(h_hnt[h,n,t] for h in H) + hr_nt[n,t])/1000
        == max_dem_n[n]*D_nt[n,t]/1000)
    end

    # Generation capacity
    @constraint(model,
        g1[g in G, n in N, t in T],
        p_gnt[g,n,t] ≤ A_gnt[g,n,t] * p̄_gn[g,n])

    # @constraint(model,
    #     g2[g in G, n in N],
    #     p̄_gn[g,n] ≤ Gmax_gn[g,n])

    # Minimum renewables share
    if specs.renewable_target
        if (specs.hydro || specs.hydro_simple)
            @constraint(model, g3,
                ((sum(p_gnt[g,n,t]*τ_t[t] for g in G_r, n in N, t in T) + sum(h_hnt[h,n,t]*τ_t[t] for h in H, n in N, t in T)) + sum(hr_nt[n,t]*τ_t[t] for n in N, t in T)) ≥
                κ * (sum(p_gnt[g,n,t]*τ_t[t] for g in G, n in N, t in T) + sum(h_hnt[h,n,t]*τ_t[t] for h in H, n in N, t in T) + sum(hr_nt[n,t]*τ_t[t] for n in N, t in T)))
        else
            @constraint(model, g3,
                sum(p_gnt[g,n,t]*τ_t[t] for g in G_r, n in N, t in T) ≥
                κ * sum(p_gnt[g,n,t]*τ_t[t] for g in G, n in N, t in T))
        end
    end

    # Maximum nuclear share
    if specs.nuclear_limit
        #TODO: include a parameter with the index(es) of nuclear sources as a parameter in io.jl
        @constraint(model, g4,
            sum(p_gnt[5,n,t] for n in N, t in T) / 1000 ≤ μ * (sum(p_gnt[g,n,t] for g in G, n in N, t in T) +
            sum(h_hnt[h,n,t] for h in H, n in N, t in T) + sum(hr_nt[n,t] for n in N, t in T)) / 1000)
    end

    #Carbon cap
    if specs.carbon_cap
        @constraint(model, g5,
            (sum(E_g[g] * sum(p_gnt[g,n,t]*τ_t[t] for n in N, t in T) / e_g[g] for g in G)) / 1000 ≤ (1-C_E) * R_E / 1000)
    end

    # Shedding upper bound
    # @constraint(model,
    #     g6[n in N, t in T],
    #     σ_nt[n,t] ≤ C̄ * max_dem_n[n]*D_nt[n,t])

    if specs.transmission
        # Transmission capacity
        @constraint(model,
            t1[l in L_ind, t in T],
            f_lt[l,t]/1000 ≤ f̄_l[l]/1000)

        @constraint(model,
            t2[l in L_ind, t in T],
            f_lt[l,t]/1000 ≥ -f̄_l[l]/1000)

        # Absolute value of f_lt
        @constraint(model,
            t3[l in L_ind, t in T],
            f_abs_lt[l,t]/1000 ≥ f_lt[l,t]/1000)

        @constraint(model,
            t4[l in L_ind, t in T],
            f_abs_lt[l,t]/1000 ≥ -f_lt[l,t]/1000)

        # @constraint(model,
        #     t6[l in L_ind, t in T],
        #     f̄_l[l] ≤ Tmax_l[l])
    end

    if specs.storage
        # Initial storage policy
        ## TODO: pass the battery storage factor via input argument
        @constraint(model,
            s0[s in S, n in N, t in [1]],
            b_snt[s,n,t]/1000 == 0.5*b̄_sn[s,n]/1000)
        # Storage capacity
        @constraint(model,
            s1[s in S, n in N, t in T],
            b_snt[s,n,t]/1000 ≤ b̄_sn[s,n]/1000)
        # Discharge limits (t = 1)
        @constraint(model,
            s2[s in S, n in N, t in [1]],
            τ_t[t]*ξ_s[s]*b⁻_snt[s,n,t]/1000 ≤ b_snt[s,n,t]/1000)
        # Discharge limits (t > 1)
        @constraint(model,
            s3[s in S, n in N, t in T[T.>1]],
            τ_t[t]*ξ_s[s]*b⁻_snt[s,n,t]/1000 ≤ τ_t[t-1]*b_snt[s,n,t-1]/1000)
        # Charge
        @constraint(model,
            s4[s in S, n in N, t in T],
            τ_t[t]*b⁺_snt[s,n,t] ≤ (b̄_sn[s,n] - b_snt[s,n,t]))
        # Storage balance
        @constraint(model,
            s5[s in S, n in N, t in T[T.>1]],
            b_snt[s,n,t] == (b_snt[s,n,t-1] + τ_t[t]*(b⁺_snt[s,n,t] - ξ_s[s]*b⁻_snt[s,n,t])))
        # Storage continuity
        @constraint(model,
            s6[s in S, n in N],
            b_snt[s,n,1] == b_snt[s,n,T[end]])
        # Storage capacity bounds
        # @constraint(model,
        #     s7[s in S, n in N],
        #     b̄_sn[s,n] ≤ Smax_sn[s,n])
    end

    if specs.ramping
        # Ramping limits
        @constraint(model,
            r1[g in G, n in N, t in T[T.>1]],
            (p_gnt[g,n,t]-p_gnt[g,n,t-1]) ≤ r⁺_g[g] * p̄_gn[g,n])
        @constraint(model,
            r2[g in G, n in N, t in T[T.>1]],
            (p_gnt[g,n,t]-p_gnt[g,n,t-1]) ≥ -r⁻_g[g] * p̄_gn[g,n])
    end

    if specs.voltage_angles
        # Voltage angles
        @constraint(model,
            v1[g in G, l in L_ind, n in N, n′ in N, t in T[T.>1]],
            (θ_nt[n,t] - θ′_nt[n′,t])*B_l[l] == p_gnt[g,n,t]-p_gnt[g,n′,t])
    end

    if specs.hydro
        # Hydro energy constraints

        # Full reservoir policy
        @constraint(model,
            h0[h in H, n in N],
            w_hnt[h,n,1]/1000 == Wmax_hn[h,n]/1000)
        # Maximum reservoir level
        @constraint(model,
            h1[h in H, n in N, t in T],
            w_hnt[h,n,t]/1000 ≤ Wmax_hn[h,n]/1000)
        # Reservoir balance
        @constraint(model,
            h2[n in N, t in T[T.>1]],
            sum(w_hnt[h,n,t] for h in H)/1000 ≤ (sum(w_hnt[h,n,t-1] for h in H) + (AH_nt[n,t-1] - sum(h_hnt[h,n,t-1] for h in H))*τ_t[t])/1000)
        # Reservoir temporal continuity
        @constraint(model,
            h3[n in N, h in H],
            w_hnt[h,n,1]/1000 == w_hnt[h,n,T[end]]/1000)
        # Define the minimum hydro flow possible to be used
        α_nt = zeros(length(N),length(T))
        for n in N, t in T
            AH_nt[n,t] + AR_nt[n,t] ≤ Fmin_n[n] ? α_nt[n,t] = 0 : α_nt[n,t] = AH_nt[n,t] + AR_nt[n,t] - Fmin_n[n]
        end
        # Minimum hydro flow constraint
        @constraint(model,
            h4[n in N, t in T],
            (sum(h_hnt[h,n,t] for h in H) + hr_nt[n,t])/1000 ≤ α_nt[n,t]/1000)
        # Maximum hydro flow (availability)
        @constraint(model,
            h5[n in N, t in T],
            sum(h_hnt[h,n,t] for h in H)/1000 ≤ AH_nt[n,t]/1000)
        # Maximum hydro flow (capacity)
        @constraint(model,
            h6[h in H, n in N, t in T],
            h_hnt[h,n,t]/1000 ≤ h̄_hn[h,n]/1000)
        # Maximum RoR hydro generation
        @constraint(model,
            h7[n in N, t in T],
            hr_nt[n,t]/1000 ≤ AR_nt[n,t]/1000)
        # Maximum hydro generation
        @constraint(model,
            h8[n in N, t in T],
            hr_nt[n,t]/1000 ≤ HRmax_n[n]/1000)
        # Maximum hydro installed capacity
        # @constraint(model,
        #     h9[h in H, n in N],
        #     h̄_hn[h,n] ≤ Hmax_hn[h,n])
        if specs.ramping
            # Hydro ramping up
            @constraint(model,
                hr1[h in H, n in N, t in T[T.>1]],
                (h_hnt[h,n,t]-h_hnt[h,n,t-1])/1000 ≤ (r⁺_h[h] * h̄_hn[h,n])/1000)
            # Hydro ramping down
            @constraint(model,
                hr2[h in H, n in N, t in T[T.>1]],
                (h_hnt[h,n,t]-h_hnt[h,n,t-1])/1000 ≥ (-r⁻_h[h] * h̄_hn[h,n])/1000)
        end
    end

    if specs.hydro_simple
        if !(specs.hydro)
        # Maximum hydro flow (capacity)
        @constraint(model,
            hs1[h in H, n in N, t in T],
            h_hnt[h,n,t]/1000 ≤ h̄_hn[h,n]/1000)
            if specs.ramping
                # Hydro ramping up
                @constraint(model,
                    hsr1[h in H, n in N, t in T[T.>1]],
                    (h_hnt[h,n,t]-h_hnt[h,n,t-1])/1000 ≤ (r⁺_h[h] * h̄_hn[h,n])/1000)
                # Hydro ramping down
                @constraint(model,
                    hsr2[h in H, n in N, t in T[T.>1]],
                    (h_hnt[h,n,t]-h_hnt[h,n,t-1])/1000 ≥ (-r⁻_h[h] * h̄_hn[h,n])/1000)
            end
        end
        # Maximum hydro flow (accumulated)
        @constraint(model,
        hs2[n in N],
        sum(τ_t[t]*h_hnt[h,n,t] for h in H, t in T)/1000 ≤ a_n[n]*sum(τ_t[t]*h̄_hn[h,n] for h in H, t in T)/1000)
    end

    return model, VariablesDict, ObjectivesDict
end


text:
module MyLib

using Dates

function isitfriday()
  Dates.dayofweek(Dates.now()) == Dates.Friday
end

end

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
#Check https://github.com/JuliaSmoothOptimizers/NLPModelsKnitro.jl/blob/master/src/NLPModelsKnitro.jl
#=
*General options*
algorithm: Indicates which algorithm to use to solve the problem
blasoption: Specifies the BLAS/LAPACK function library to use for basic vector
            and matrix computations
cg_maxit: Determines the maximum allowable number of inner conjugate gradient
          (CG) iterations
cg_pmem: Specifies number of nonzero elements per hessian column when computing
         preconditioner
cg_precond: Specifies whether or not to apply preconditioning during
            CG iterations in barrier algorithms
cg_stoptol: Relative stopping tolerance for CG subproblems
convex: Identify convex models and apply specializations often beneficial for convex models
delta: Specifies the initial trust region radius scaling factor
eval_fcga: Specifies that gradients are provided together with functions in one callback
honorbnds: Indicates whether or not to enforce satisfaction of simple variable bounds
initpenalty: Initial penalty value used in Knitro merit function
linesearch_maxtrials: Indicates the maximum allowable number of trial points during the linesearch
linesearch: Indicates which linesearch strategy to use for the Interior/Direct or SQP algorithm
linsolver_ooc: Indicates whether to use Intel MKL PARDISO out-of-core solve of linear systems
linsolver: Indicates which linear solver to use to solve linear systems arising in Knitro algorithms
linsolver_maxitref
linsolver_pivottol
objrange: Specifies the extreme limits of the objective function for purposes of determining unboundedness
          Default value: 1.0e20
presolve: Determine whether or not to use the Knitro presolver
presolve_initpt: Controls whether Knitro presolver can shift user-supplied initial point
restarts: Specifies whether to enable automatic restarts
restarts_maxit: Maximum number of iterations before restarting when restarts are enabled
scale: Specifies whether to perform problem scaling
soc: Specifies whether or not to try second order corrections (SOC)
strat_warm_start: Specifies whether or not to invoke a warm-start strategy

*Derivatives options*
bfgs_scaling: Specifies the initial scaling for the BFGS or L-BFGS Hessian approximation
derivcheck: Determine whether or not to perform a derivative check on the model
gradopt: Specifies how to compute the gradients of the objective and constraint functions
hessian_no_f: Determines whether or not to allow Knitro to request Hessian
              evaluations without the objective component included.
hessopt: Specifies how to compute the (approximate) Hessian of the Lagrangian
         1 (exact) User provides a routine for computing the exact Hessian. (Default)
         4 (product_findiff) Knitro computes Hessian-vector products using finite-differences.
         5 (product) User provides a routine to compute the Hessian-vector products.
         6 (lbfgs) Knitro computes a limited-memory quasi-Newton BFGS Hessian
                   (its size is determined by the option lmsize).
lmsize: Specifies the number of limited memory pairs stored when approximating the Hessian

*Termination options*
feastol: Specifies the final relative stopping tolerance for the feasibility error.
         Default value: 1.0e-6
feastol_abs: Specifies the final absolute stopping tolerance for the feasibility error.
             Default value: 1.0e-3
fstopval: Used to implement a custom stopping condition based on the objective function value
ftol: The optimization process will terminate if feasible and the relative change
      in the objective function is less than ftol
      Default value: 1.0e-15
ftol_iters: The optimization process will terminate if the relative change in
            the objective function is less than ftol for ftol_iters consecutive
            feasible iterations. Default value: 5

maxfevals: Specifies the maximum number of function evaluations before termination.
           Default value: -1 (unlimited)
maxit: Specifies the maximum number of iterations before termination
       0 is default value, let Knitro set it (10000 for LP/NLP)
maxtime_cpu: Specifies, in seconds, the maximum allowable CPU time before termination.
             Default value: 1.0e8
maxtime_real: Specifies, in seconds, the maximum allowable real time before termination.
              Default value: 1.0e8
opttol: Specifies the final relative stopping tolerance for the KKT (optimality) error
        Default value: 1.0e-6
opttol_abs: Specifies the final absolute stopping tolerance for the KKT (optimality) error
            Default value: 1.0e-3
xtol: The optimization process will terminate if the relative change of the
      solution point estimate is less than xtol. Default value: 1.0e-12
xtol_iters: Number of consecutive iterations where change of the solution point
            estimate is less than xtol before Knitro stops.
            Default is 1

*Output options*
out_hints: Print diagnostic hints (e.g. on user option settings) after solving
           0 no prints, and 1 prints (default in Knitro)
outlev: Controls the level of output produced by Knitro
        0 (none) Printing of all output is suppressed.
        1 (summary) Print only summary information.
        2 (iter_10) Print basic information every 10 iterations. (default in Knitro)
=#

@init begin
  @require KNITRO = "67920dd8-b58e-52a8-8622-53c4cffbe346" begin
    @require NLPModelsKnitro = "bec4dd0d-7755-52d5-9a02-22f0ffc7efcb" begin
      is_knitro_installed = true
      import NLPModelsKnitro: knitro

      """
          `knitro(stp::NLPStopping; subsolver_verbose::Int = 0, kwargs...)`

      Stopping-version of the `knitro` function from NLPModelsKnitro.jl.
      This function doesn't calls `stop!` after `knitro` call.
      Use the `KnitroSolver` structure to fill-in the gradient and objective in the state.

      `subsolver_verbose` corresponds to `print_level` argument in `knitro`.
      Other keyword arguments are passed to the `knitro` call.
      Selection of possible [options](https://www.artelys.com/docs/knitro/3_referenceManual/userOptions.html).

      It requires `using KNITRO, NLPModelsKnitro`.
      """
      function NLPModelsKnitro.knitro(
        stp::NLPStopping;
        convex::Int = -1, #let Knitro deal with it :)
        objrange::Real = stp.meta.unbounded_threshold,
        hessopt::Int = 1,
        feastol::Real = stp.meta.rtol,
        feastol_abs::Real = stp.meta.atol,
        opttol::Real = stp.meta.rtol,
        opttol_abs::Real = stp.meta.atol,
        maxfevals::Int = min(stp.meta.max_cntrs[:neval_sum], typemax(Int32)),
        maxit::Int = 0, #stp.meta.max_iter
        maxtime_real::Real = stp.meta.max_time,
        out_hints::Int = 0,
        subsolver_verbose::Int = 0, #1 to see everything
        algorithm::Int = 0, # *New* 2
        ftol::Real = 1.0e-15, # *New*
        ftol_iters::Int = 5, # *New*
        xtol::Real = stp.meta.atol, # 1.0e-12, # *New*
        xtol_iters::Int = 2, # 3 # *New*
        kwargs...,
      )
        @assert -1 ≤ convex ≤ 1
        @assert 1 ≤ hessopt ≤ 7
        @assert 0 ≤ out_hints ≤ 1
        @assert 0 ≤ subsolver_verbose ≤ 6
        @assert 0 ≤ maxit

        stp.meta.start_time = time()
        nlp = stp.pb
        #y0 = stp.current_state.lambda #si défini
        #z0 = stp.current_state.mu #si défini
        solver = NLPModelsKnitro.KnitroSolver(
          nlp,
          x0 = stp.current_state.x,
          objrange = objrange,
          feastol = feastol,
          feastol_abs = feastol_abs,
          opttol = opttol,
          opttol_abs = opttol_abs,
          maxfevals = maxfevals,
          maxit = maxit,
          maxtime_real = maxtime_real,
          out_hints = out_hints,
          algorithm = algorithm,
          ftol = ftol,
          ftol_iters = ftol_iters,
          xtol = xtol,
          xtol_iters = xtol_iters,
          outlev = subsolver_verbose;
          kwargs...,
        )
        stats = NLPModelsKnitro.knitro!(nlp, solver)
        #@show stats.status, stats.solver_specific[:internal_msg]
        #if stats.status ∉ (:unbounded, :exception, :unknown) #∈ (:first_order, :acceptable)
        stp.current_state.x = stats.solution
        stp.current_state.fx = stats.objective
        stp.current_state.gx = KNITRO.KN_get_objgrad_values(solver.kc)[2]
        #norm(stp.current_state.gx, Inf)#stats.dual_feas #TODO: this is for unconstrained problem!!
        stp.current_state.mu = stats.multipliers_L
        if stp.pb.meta.ncon > 0
          stp.current_state.cx = KNITRO.KN_get_con_values(solver.kc)
          stp.current_state.Jx = sparse(
            KNITRO.KN_get_jacobian_values(solver.kc)[1],
            KNITRO.KN_get_jacobian_values(solver.kc)[2],
            KNITRO.KN_get_jacobian_values(solver.kc)[3],
            stp.pb.meta.ncon,
            stp.pb.meta.nvar,
          )
          stp.current_state.lambda = stats.multipliers
        end
        stp.current_state.current_score = max(stats.dual_feas, stats.primal_feas)
        #end
        #Update the meta boolean with the output message
        stp = stats_status_to_meta!(stp, stats)
        #@show status(stp, list = true)
        if status(stp) == :Unknown
          @warn "Error in StoppingInterface statuses: return status is $(stats.status)"
          #print(stats)
        end

        stp.meta.nb_of_stop = stats.iter
        Stopping._update_time!(stp.current_state, time()) # stats.elapsed_time

        return stp #would be better to return the stats somewhere
      end
    end
  end
end


text:
module AuditoryBistabilityLE
using CorticalSpecrtralTemporalResponses
using DSP
using DataFrames
using Requires
using Statistics
using LinearAlgebra
using SparseArrays
using SpecialFunctions
using TOML
using MetaArrays
using Unitful

import CorticalSpecrtralTemporalResponses: Δt, Δf, times, freqs, scales, rates, asseconds, asHz

export adaptmi, drift, scale_weighting, ncomponents, nunits, CoherenceModel,
    fusion_ratio, object_SNR, mask, scene_object_ratio,
    object_SNR2, ab_match, mean_spect, mean_spect2, AdaptMI

using ProgressMeter
next!(x::Progress) = ProgressMeter.next!(x)
next!(x::Nothing) = nothing

# core model implementation
########################################
# implements inhibition, adaptation and noise
include("adaptmi.jl")
# implements inhibition interactions for each hierarchical level
include("cortmi.jl")
# utility functions for handling the windowing of a signal
include("windowing.jl")
# cohere is step 1 of the object-level analysis (temporal coherence)
include("cohere.jl")
# the parts of cohere specific to the NMF implementation
include("nmf.jl")
# general purpose functions for tracking (step 2 of the object-level analysis)
include("tracking.jl")
# function specific to using prior source models for object tracking
# (other approachs were examined)
include("prior_tracking.jl")
# function for considering multiple priors
# (i.e. the different β, κ parametesr from the paper)
include("multi_prior_tracking.jl")

# all the functions that use the model implementation to do stuff
########################################
# miscellaneous utility functions
include("util.jl")
# ABA stimulus generation
include("stim.jl")
# hueristic estimate of the number of sources present in an ABA_ stimulus
include("bandwidth_ratio.jl")
#  code to compute percept lengths using bandwidth ratio
include("lengths.jl")
# applies the bistable model to a specific set of conditions
include("biapply.jl")
# where everything comes together: all stages of the model are computed here
include("bimodel.jl")
# compress the model output so it takes up less space when stored
include("compress.jl")

# data visualization
include("plot_axes.jl")
function __init__()
  @require RCall="6f49c342-dc21-5d91-9882-a32aef131414" include("rplots.jl")
end

end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
using Test

include("sc_ball.jl")

@testset "Penetration Tests" begin
    @test all(getindex.(sc_ball.q,4) .- getindex.(sc_ball.q,3) .- 0.1 .>= -1e-3)
end

# @testset "Energy Tests" begin
#     energy(q, u) = 1/2 * mass * (u[1]^2 + u[2]^2) + mass * g * q[2]
#     initial_energy = energy(bouncing_ball.q[1], bouncing_ball.u[1])

#     @test all(energy.(bouncing_ball.q, bouncing_ball.u) .<= initial_energy)
# end

# include("three_bouncing_balls.jl")

# @testset "Multiple contact tests" begin
#     for i = 1:3
#         @test all(getindex.(three_bb.q,i) .>= -5e-3)
#     end
# end

# @testset "Contact force tests" begin
#     for i = 1:3
#         @test minimum(getindex.(three_bb.Λ,i)) >= -1e-10
#     end
# end


text:
module MimiRICE2010

using Mimi
using XLSX: readxlsx

include("helpers.jl")
include("marginaldamage.jl")
include("parameters.jl")

include("components/climatedynamics_component.jl")
include("components/co2cycle_component.jl")
include("components/damages_component.jl")
include("components/emissions_component.jl")
include("components/grosseconomy_component.jl")
include("components/neteconomy_component.jl")
include("components/radiativeforcing_component.jl")
include("components/slr_component.jl")
include("components/slrdamages_component.jl")
include("components/welfare_component.jl")

export constructrice, getrice

const model_years = 2005:10:2595

function constructrice(p)

    m = Model()
    set_dimension!(m, :time, model_years)
    set_dimension!(m, :regions, ["US", "EU", "Japan", "Russia", "Eurasia", "China", "India", "MidEast", "Africa", "LatAm", "OHI", "OthAsia"])

    #--------------------------------------------------------------------------
    # Add components in order
    #--------------------------------------------------------------------------

    add_comp!(m, grosseconomy, :grosseconomy)
    add_comp!(m, emissions, :emissions)
    add_comp!(m, co2cycle, :co2cycle)
    add_comp!(m, radiativeforcing, :radiativeforcing)
    add_comp!(m, climatedynamics, :climatedynamics)
    add_comp!(m, sealevelrise, :sealevelrise)
    add_comp!(m, sealeveldamages, :sealeveldamages)
    add_comp!(m, damages, :damages)
    add_comp!(m, neteconomy, :neteconomy)
    add_comp!(m, welfare, :welfare)

    #--------------------------------------------------------------------------
    # Make internal parameter connections
    #--------------------------------------------------------------------------

    # GROSS ECONOMY COMPONENT
    connect_param!(m, :grosseconomy, :I, :neteconomy, :I)

    # EMISSIONS COMPONENT
    connect_param!(m, :emissions, :YGROSS, :grosseconomy, :YGROSS)

    # CO2 CYCLE COMPONENT
    connect_param!(m, :co2cycle, :E, :emissions, :E)

    # RADIATIVE FORCING COMPONENT
    connect_param!(m, :radiativeforcing, :MAT, :co2cycle, :MAT)
    connect_param!(m, :radiativeforcing, :MATSUM, :co2cycle, :MATSUM)

    # CLIMATE DYNAMICS COMPONENT
    connect_param!(m, :climatedynamics, :FORC, :radiativeforcing, :FORC)

    # SEA LEVEL RISE COMPONENT(S)
    connect_param!(m, :sealevelrise, :TATM, :climatedynamics, :TATM)
    connect_param!(m, :sealeveldamages, :TOTALSLR, :sealevelrise, :TOTALSLR)
    connect_param!(m, :sealeveldamages, :YGROSS, :grosseconomy, :YGROSS)

    # DAMAGES COMPONENT
    connect_param!(m, :damages, :TATM, :climatedynamics, :TATM)
    connect_param!(m, :damages, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :damages, :SLRDAMAGES, :sealeveldamages, :SLRDAMAGES)

    # NET ECONOMY COMPONENT
    connect_param!(m, :neteconomy, :YGROSS, :grosseconomy, :YGROSS)
    connect_param!(m, :neteconomy, :DAMFRAC, :damages, :DAMFRAC)
    connect_param!(m, :neteconomy, :DAMAGES, :damages, :DAMAGES)
    connect_param!(m, :neteconomy, :ABATECOST, :emissions, :ABATECOST)

    # WELFARE COMPONENT
    connect_param!(m, :welfare, :CPC, :neteconomy, :CPC)

    #--------------------------------------------------------------------------
    # Set external parameter values
    #--------------------------------------------------------------------------
    for (name, value) in p
        set_param!(m, name, value)
    end

    return m
end #function

function get_model(;datafile=joinpath(@__DIR__, "..", "data", "RICE_2010_base_000.xlsm"))
    params = getrice2010parameters(datafile)

    m = constructrice(params)

    return m
end #function

getrice = get_model     # Maintain the old `getrice` function name in addition to the standard MimiRICE2010.get_model

end #module


text:
# TODO: code duplication could be reduced if both solvers supported the same interface
function estimate(
    solver::InverseKKTConstraintSolver;
    dataset,
    control_system,
    player_cost_models,
    solver_attributes,
    expected_observation = identity,
    estimator_name = (expected_observation === identity ? "Ours Full" : "Ours Partial"),
    solver_kwargs...,
)

    @showprogress pmap(enumerate(dataset)) do (observation_idx, d)
        observation_model = (; d.σ, expected_observation)

        converged, estimate, opt_model = solve_inverse_game(
            solver,
            expected_observation(d.x);
            control_system,
            observation_model,
            player_cost_models,
            solver_attributes,
            solver_kwargs...,
            # NOTE: This estimator does not use any information beyond the state observation!
        )
        converged || @warn "conKKT did not converge on observation $observation_idx."

        merge(estimate, (; converged, observation_idx, estimator_name))
    end
end

function estimate(
    solver::InverseKKTResidualSolver;
    dataset,
    control_system,
    player_cost_models,
    solver_attributes,
    expected_observation = identity,
    estimator_name = (expected_observation === identity ? "Baseline Full" : "Baseline Partial"),
)
    @showprogress pmap(enumerate(dataset)) do (observation_idx, d)
        observation_model = (; d.σ, expected_observation)
        smoothed_observation = let
            pre_solve_converged, pre_solve_solution = InversePreSolve.pre_solve(
                # pre-filter for baseline receives one extra state observation to avoid
                # unobservability of the velocity at the end-point.
                expected_observation([d.x d.x_extra]),
                nothing;
                control_system,
                observation_model,
                solver_attributes,
            )
            @assert pre_solve_converged
            # Filtered sequence is truncated to the original length to give all methods the same
            # number of data-points for inference.
            (;
                x = pre_solve_solution.x[:, 1:size(d.x, 2)],
                u = pre_solve_solution.u[:, 1:size(d.u, 2)],
            )
        end

        converged, estimate, opt_model = solve_inverse_game(
            solver,
            smoothed_observation.x,
            smoothed_observation.u;
            control_system,
            player_cost_models,
            solver_attributes,
        )
        converged || @warn "resKKT did not converge on observation $observation_idx."

        merge(estimate, (; converged, observation_idx, estimator_name, smoothed_observation))
    end
end


text:
abstract type SeaLevel end
abstract type Bintanja   end
abstract type Spratt450  end
abstract type Spratt800  end

struct SeaLevelOpts{T}
    path::String
end

sea_level_curve(sealevel::SeaLevelOpts{Spratt450}) = string(sealevel.path,"Spratt2016-450kyr.txt");
sea_level_curve(sealevel::SeaLevelOpts{Spratt800}) = string(sealevel.path,"Spratt2016-800kyr.txt");
sea_level_curve(sealevel::SeaLevelOpts{Bintanja}) = string(sealevel.path,"Bintanja-3Ma.txt");

sea_level_starting_time(sea_level_file::SeaLevelOpts{Bintanja}) = 3e6
sea_level_starting_time(sea_level_file::SeaLevelOpts{Spratt800}) = 800e3
sea_level_starting_time(sea_level_file::SeaLevelOpts{Spratt450}) = 450e3

function get_sea_level(sealevelopts::SeaLevelOpts)

    sl = sea_level_curve(sealevelopts)
    curves0 = readdlm(sl)::Matrix{Float64}
    id = @. isnan(curves0[:,2]) != 1 # get rid of NaNs
    sealevel = curves0[id,:]
    nn = size(sealevel,1)
    age = Vector{Float64}(undef,nn)
    sea_level = Vector{Float64}(undef,nn)
    correction = sealevel[1,2]
    @simd for ii = 1:nn
        @inbounds sea_level[ii] = sealevel[ii,2] - correction; # Sea Level,,,meters above present day,,climate reconstructions,,Scaled first principal component of seven sea level reconstructions (0-430 ka),N
        @inbounds age[ii] = sealevel[ii,1]; # now in [yrs], originally in [kyrs]
    end

    return sea_level, age

end

function sea_level_corrections(sea_lvl_curve, sea_age, starting_time)
    sea_age = sea_age[end:-1:1] # reverse age array -> 0 = oldest age
    id_time = sea_age .>= maximum(sea_age) .- starting_time
    sea_age = sea_age[id_time]
    sea_lvl_curve = sea_lvl_curve[id_time]
    fsea = interpolate((reverse(sea_age),), reverse(sea_lvl_curve), Gridded(Linear()))
    return sea_lvl_curve, sea_age, fsea
end

function read_uplift()
    fname = "/home/albert/Dropbox/Riverini/DevelopNoTouch/uplift_rates/accre1_VertVel_870-871.txt"
    num = readdlm(fname)
    # correct time serie
    t = @. num[:,1] - num[1,1]; # [year]
    # mm/year to m/year
    uplift = @. num[:,2] / 1e3; # [mm/year] -> [m.year]

    return t, uplift

end


text:
module test_julia_asyncevent_timedwait

using Test

@test timedwait(() -> true, 0.0) === :ok

state = false
@test timedwait(() -> state, 0.0) === :timed_out

end # module test_julia_asyncevent_timedwait


text:
## Test π0 estimator methods ##
module Test_pi0

using MultipleTesting
using Base.Test
using StatsBase

## test case: deterministic
p0 = collect(0.01:0.01:1)
p1 = p0 .^ 10
p = [p0; p1] ## unsorted
pi0 = length(p0) / length(p)

function unsort(x)
    y = copy(x)
    while issorted(y)
        sample!(x, y, replace = false)
    end
    return y
end

lambdas = collect(0.05:0.05:0.95)

## storey_pi0 ##
println(" ** ", "storey_pi0")
@test_throws MethodError storey_pi0()
@test_approx_eq storey_pi0(p, 0.2) 0.6
@test_approx_eq storey_pi0(p, 0.0) 1.0
#@test_throws DomainError storey_pi0(p, 1.) ## CHCK

p_unsort = unsort(p)
@test !issorted(p_unsort)
@test_approx_eq storey_pi0(p_unsort, 0.2) 0.6
@test !issorted(p_unsort)

@test issubtype(typeof(Storey()), Pi0Estimator)
@test issubtype(typeof(Storey(0.5)), Pi0Estimator)
@test_approx_eq estimate_pi0(p, Storey(0.2)) 0.6
@test_approx_eq estimate_pi0(p, Storey(0.0)) 1.0
@test_approx_eq estimate_pi0(p, Storey(1.0)) 1.0

@test_throws DomainError Storey(-0.1)
@test_throws DomainError Storey(1.1)

## bootstrap_pi0 ##
println(" ** ", "bootstrap_pi0")
@test_throws MethodError bootstrap_pi0()

@test_approx_eq bootstrap_pi0(p, lambdas) 0.6
@test_approx_eq bootstrap_pi0(p0, lambdas) 1.0
@test_approx_eq bootstrap_pi0(p1, lambdas) 0.15

## unsorted lambdas
lambdas_unsort = unsort(lambdas)
@test !issorted(lambdas_unsort)
@test_approx_eq bootstrap_pi0(p, lambdas_unsort) 0.6
@test !issorted(lambdas_unsort)

## with default 'lambdas'
@test_approx_eq bootstrap_pi0(p) 0.6
@test_approx_eq bootstrap_pi0(p0) 1.0
@test_approx_eq bootstrap_pi0(p1) 0.15

@test issubtype(typeof(StoreyBootstrap()), Pi0Estimator)
@test_approx_eq estimate_pi0(p, StoreyBootstrap()) 0.6
@test_approx_eq estimate_pi0(p0, StoreyBootstrap()) 1.0
@test_approx_eq estimate_pi0(p1, StoreyBootstrap()) 0.15

@test_throws DomainError StoreyBootstrap(lambdas, -0.1)
@test_throws DomainError StoreyBootstrap(lambdas, 1.1)
@test_throws MethodError StoreyBootstrap(0.5)
@test_throws MethodError StoreyBootstrap(lambdas)
@test_throws MethodError StoreyBootstrap(0.5, lambdas)

## lsl_pi0 ##
println(" ** ", "lsl_pi0")
@test_approx_eq lsl_pi0(p) MultipleTesting.lsl_pi0_vec(p)
@test_approx_eq lsl_pi0(p0) MultipleTesting.lsl_pi0_vec(p0)
@test_approx_eq lsl_pi0(p1) MultipleTesting.lsl_pi0_vec(p1)

## checked against structSSI::pi0.lsl
@test_approx_eq_eps lsl_pi0(p) 0.62 1e-2
@test_approx_eq_eps lsl_pi0(p0) 1.0 1e-2
@test_approx_eq_eps lsl_pi0(p1) 0.16 1e-2

@test issubtype(typeof(LeastSlope()), Pi0Estimator)
@test_approx_eq_eps estimate_pi0(p, LeastSlope()) 0.62 1e-2
@test_approx_eq_eps estimate_pi0(p0, LeastSlope()) 1.0 1e-2
@test_approx_eq_eps estimate_pi0(p1, LeastSlope()) 0.16 1e-2

@test_throws MethodError LeastSlope(0.1)

## unsorted p-values
p_unsort = unsort(p)
@test !issorted(p_unsort)
@test_approx_eq_eps lsl_pi0(p_unsort) 0.62 1e-2
@test !issorted(p_unsort)

p_unsort = unsort(p)
@test !issorted(p_unsort)
@test_approx_eq_eps MultipleTesting.lsl_pi0_vec(p_unsort) 0.62 1e-2
@test !issorted(p_unsort)

## Oracle ##
println(" ** ", "oracle")
@test estimate_pi0(p, Oracle(0.5)) == 0.5
@test estimate_pi0(p0, Oracle(0.6)) == 0.6
@test estimate_pi0(p1, Oracle()) == 1.0

## twostep_pi0 ##
println(" ** ", "twostep_pi0")

alpha = 0.05

## checked against mutoss::TSBKY_pi0_est
@test_approx_eq twostep_pi0(p, alpha) 0.665
@test_approx_eq twostep_pi0(p0, alpha) 1.0
@test_approx_eq twostep_pi0(p1, alpha) 0.29

@test issubtype(typeof(TwoStep(alpha)), Pi0Estimator)
@test_approx_eq estimate_pi0(p, TwoStep()) 0.665
@test_approx_eq estimate_pi0(p, TwoStep(alpha)) 0.665
@test_approx_eq estimate_pi0(p0, TwoStep(alpha)) 1.0
@test_approx_eq estimate_pi0(p1, TwoStep(alpha)) 0.29

@test_approx_eq estimate_pi0(p, TwoStep(0.1)) 0.63
@test_approx_eq estimate_pi0(p, TwoStep(0.0)) 1.0
@test_approx_eq estimate_pi0(p, TwoStep(1.0)) 0.415

## unsorted p-values
p_unsort = unsort(p)
@test !issorted(p_unsort)
@test_approx_eq twostep_pi0(p_unsort, alpha) 0.665
@test !issorted(p_unsort)

## rightboundary_pi0 ##
println(" ** ", "rightboundary_pi0")

# R code tested against (note this uses equidistant grid so we have to use
# λseq as in Storey for comparison

# library(pi0)
# p0 <- seq(0.01,1, by=0.01)
# histf1(p0, rightBoundary=TRUE)
# [1] 1
# p1 <- p0^10
# histf1(p1, rightBoundary=TRUE)
# [1] 0.1428571
# p <- c(p0,p1)# histf1(p, rightBoundary=TRUE)
# [1] 0.5714286

# only eps because pi0 package uses right closed histograms
@test_approx_eq_eps rightboundary_pi0(p, lambdas) 0.5714286 0.02
@test_approx_eq rightboundary_pi0(p0, lambdas) 1.0
@test_approx_eq_eps rightboundary_pi0(p1, lambdas) 0.1428571 10.0^(-7)

@test_approx_eq rightboundary_pi0(p, lambdas) rightboundary_pi0(p, unsort(lambdas))

@test issubtype(typeof(RightBoundary(lambdas)), Pi0Estimator)
@test_approx_eq_eps estimate_pi0(p, RightBoundary(lambdas)) 0.5714286 0.02
@test_approx_eq estimate_pi0(p0, RightBoundary(lambdas)) 1.0
@test_approx_eq_eps estimate_pi0(p1, RightBoundary(lambdas)) 0.1428571 10.0^(-7)

# not checked against R implementation but should hold (check for default lambda grid)
@test_approx_eq estimate_pi0(p0, RightBoundary()) 1.0

@test issubtype(typeof(RightBoundary()), Pi0Estimator)

p_unsort = unsort(p1)
@test !issorted(p_unsort)
@test_approx_eq_eps rightboundary_pi0(p_unsort, lambdas) 0.1428571 10.0^(-7)
@test !issorted(p_unsort)

## censoredBUM_pi0 ##
println(" ** ", "censoredBUM_pi0")

@test_approx_eq_eps estimate_pi0(p, CensoredBUM()) 0.55797 1e-5
@test_approx_eq estimate_pi0(p, CensoredBUM()) MultipleTesting.cbum_pi0_naive(p)[1]
@test_approx_eq_eps estimate_pi0(p0, CensoredBUM()) 1.0 1e-5
@test_approx_eq_eps estimate_pi0(p1, CensoredBUM()) 0.11608 2e-5
@test_approx_eq_eps estimate_pi0(ones(50), CensoredBUM()) 1.0 1e-5

## test case that does not converge
stderr_dump = redirect_stderr()
@test isnan(estimate_pi0(p, CensoredBUM(0.5, 0.05, 1e-6, 2)))

@test issubtype(typeof(CensoredBUM()), Pi0Estimator)
@test issubtype(typeof(CensoredBUM(0.2, 0.1)), Pi0Estimator)

@test_throws DomainError CensoredBUM(-0.5, 0.05)
@test_throws DomainError CensoredBUM(1.5, 0.05)
@test_throws DomainError CensoredBUM(0.5, -0.05)
@test_throws DomainError CensoredBUM(0.5, 1.05)

@test_throws DomainError CensoredBUM(0.5, 0.05, -1e-6, 100)
@test_throws DomainError CensoredBUM(0.5, 0.05, 1.1, 100)
@test_throws DomainError CensoredBUM(0.5, 0.05, 1e-6, -10)

f = fit(CensoredBUM(), p)
@test issubtype(typeof(f), CensoredBUMFit)
@test_approx_eq_eps f.π0 0.55797 1e-5

pi0_est, pars, is_converged = MultipleTesting.cbum_pi0(ones(50))
@test_approx_eq pi0_est 1.0
@test is_converged

pi0_est, pars, is_converged = MultipleTesting.cbum_pi0_naive(ones(50))
@test isnan(pi0_est)
@test !is_converged

## BUM_pi0 ##
## needs better test cases and reference values
println(" ** ", "BUM_pi0")

@test_approx_eq_eps estimate_pi0(p, BUM(0.5)) 0.55528 1e-5
@test_approx_eq_eps estimate_pi0(p, BUM()) 0.55528 1e-5
@test_approx_eq_eps estimate_pi0(p0, BUM()) 1.0 1e-5
@test_approx_eq_eps estimate_pi0(p1, BUM()) 0.10874 1e-5

## test case that does not converge
@test isnan(estimate_pi0(p, BUM(0.5, 1e-6, 2)))

@test issubtype(typeof(BUM()), Pi0Estimator)

@test_throws DomainError BUM(-0.5)
@test_throws DomainError BUM(1.5)

## Flat grenander ##
println(" ** ", "FlatGrenander_pi0")

FlatGrenander = MultipleTesting.FlatGrenander

@test issubtype(typeof(FlatGrenander()), Pi0Estimator)

pu = collect(0.1:0.05:0.9)
@test_approx_eq estimate_pi0(pu, FlatGrenander()) 1.0
@test_approx_eq estimate_pi0(pu.^0.5, FlatGrenander()) 1.0

@test_approx_eq estimate_pi0(p0, FlatGrenander()) 1.0
@test estimate_pi0(p1, FlatGrenander()) < 0.15
@test_approx_eq_eps estimate_pi0(p, FlatGrenander()) pi0 0.1

## longest constant interval: low level
lci = MultipleTesting.longest_constant_interval

p = [0:0.1:1;];
f = [1.5, 1.5, 1.5, 1.5, 1.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.1];
@test_approx_eq lci(p, f) 0.5

f = [1.5, 1.5, 1.5, 1.5, 1.5, 1.2, 1.2, 1.2, 0.2, 0.2, 0.1];
@test_approx_eq lci(p, f) 0.2

f = [0.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.1];
@test_approx_eq lci(p, f) 0.2

f = [0.5, 0.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.2, 0.2, 0.2, 0.1];
@test_approx_eq lci(p, f) 0.5

p = [0.1, 0.3, 0.5, 0.9];
f = [0.5, 0.5, 0.2, 0.2];
@test_approx_eq lci(p, f) 0.2

p = [0.1, 0.5, 0.7, 0.9];
f = [0.5, 0.5, 0.2, 0.2];
@test_approx_eq lci(p, f) 0.5

end


text:
export cse, jlcode, genfunc, sympi

using SymPy
import Base: sign #, zero, one
import LinearAlgebra: dot

sign(x::Sym) = sympy_meth(:sign, x)

# const _symzero = oftype(Sym, 0)
# zero(::Type{Sym}) = _symzero
#
# const _symone = oftype(Sym, 1)
# one(::Type{Sym}) = _symone

const sympi = SymPy.sympy[:pi]

# import SymPy.cse
#
# function cse(e::Sym; symbols=nothing, order="none")
#     cseout = sympy_meth(:cse, e.x, symbols=symbols, order=order)
#     (convert(Vector{Tuple{Sym,Sym}}, cseout[1]), cseout[2][1])
# end
#
# function cse(A::Array{Sym}; symbols=nothing, order="none")
#     A = convert(SymMatrix, A)
#     cseout = sympy_meth(:cse, A.x, symbols=symbols, order=order)
#     (convert(Vector{Tuple{Sym,Sym}}, cseout[1]), convert(Array{Sym}, cseout[2][1]))
# end
retype_cse(symcode) = (convert(Vector{Tuple{Sym,Sym}}, symcode[1]), [x for x in symcode[2]])
export retype_cse

function jlcode(expr::Sym)
    codestr = SymPy.sympy_meth(:ccode, expr)
    codestr = replace(codestr, "pow" => "^")
    Meta.parse(codestr)
end

jlcode(exprvec::Vector{Sym}) = Expr(:vcat, [jlcode(e) for e in exprvec]...)

jlcode(exprmat::Matrix{Sym}) =
    Expr(:vcat, [Expr(:row, [jlcode(exprmat[i,j]) for j in 1:size(exprmat,2)]...) for i in 1:size(exprmat,1)]...)

assigntoarray(array::Symbol, exprvec::Vector{Sym}) =
    [Expr(:(=), Expr(:ref, array, i), jlcode(exprvec[i])) for i in 1:length(exprvec)]

assigntoarray(array::Symbol, exprmat::Matrix{Sym}) =
    [Expr(:(=), Expr(:ref, array, i, j), jlcode(exprmat[i, j])) for j in 1:size(exprmat,2) for i in 1:size(exprmat,1)]

function genfunc(name, args::Vector, ret, ivs::Vector{Tuple{Sym,Sym}}; inbounds=false, overret=nothing)
    body_ivs = [Expr(:(=), jlcode(iv), jlcode(e)) for (iv,e) in ivs]
    if string(name)[end] == '!'
        out = Symbol(args[1])
        body_ret = assigntoarray(out, ret)
        push!(body_ret, Expr(:return, out))
    else
        if overret === nothing
            body_ret = [Expr(:return, jlcode(ret))]
        else
            body_ret = [Expr(:return, Expr(overret.head, overret.args..., jlcode(ret)))]
        end
    end
    body = Expr(:block, [body_ivs; body_ret]...)
    if inbounds
        body = Expr(:macrocall, Symbol("@inbounds"), body)
    end
    if name == :(->)
        Expr(:(->), Expr(:tuple, args...), body)
    else
        Expr(:function, Expr(:call, name, args...), body)
    end
end

genfunc(name, args::Vector, ret::Tuple{Vector{Tuple{Sym,Sym}}, T}; inbounds=false, overret=nothing) where T = genfunc(name, args::Vector, ret[2], ret[1]; inbounds=inbounds, overret=overret)


text:
using Test
using Flux: onehotbatch, σ

using Flux.Losses: mse, label_smoothing, crossentropy, logitcrossentropy, binarycrossentropy, logitbinarycrossentropy
using Flux.Losses: xlogx, xlogy

# group here all losses, used in tests
const ALL_LOSSES = [Flux.Losses.mse, Flux.Losses.mae, Flux.Losses.msle,
                    Flux.Losses.crossentropy, Flux.Losses.logitcrossentropy,
                    Flux.Losses.binarycrossentropy, Flux.Losses.logitbinarycrossentropy,
                    Flux.Losses.kldivergence,
                    Flux.Losses.huber_loss,
                    Flux.Losses.tversky_loss,
                    Flux.Losses.dice_coeff_loss,
                    Flux.Losses.poisson_loss,
                    Flux.Losses.hinge_loss, Flux.Losses.squared_hinge_loss]

@testset "xlogx & xlogy" begin
  @test iszero(xlogx(0))
  @test isnan(xlogx(NaN))
  @test xlogx(2) ≈ 2.0 * log(2.0)
  @inferred xlogx(2)
  @inferred xlogx(0)

  @test iszero(xlogy(0, 1))
  @test isnan(xlogy(NaN, 1))
  @test isnan(xlogy(1, NaN))
  @test isnan(xlogy(NaN, NaN))
  @test xlogy(2, 3) ≈ 2.0 * log(3.0)
  @inferred xlogy(2, 3)
  @inferred xlogy(0, 1)
end

# First, regression-style y's
y = [1, 1, 0, 0]
ŷ = [.9, .1, .1, .9]

@testset "mse" begin
  @test mse(ŷ, y) ≈ (.1^2 + .9^2)/2
end

@testset "mae" begin
  @test Flux.mae(ŷ, y) ≈ 1/2
end

@testset "huber_loss" begin
  @test Flux.huber_loss(ŷ, y) ≈ 0.20500000000000002
end

y = [123.0,456.0,789.0]
ŷ = [345.0,332.0,789.0]
@testset "msle" begin
  @test Flux.msle(ŷ, y) ≈ 0.38813985859136585
end

# Now onehot y's
y = onehotbatch([1, 1, 0, 0], 0:1)
y_smoothed = label_smoothing(y, 0.1)
ŷ = [.1 .9; .9 .1; .9 .1; .1 .9]'
v = log(.1 / .9)
logŷ = [v 0.0; 0.0 v; 0.0 v; v 0.0]'
lossvalue = 1.203972804325936
lossvalue_smoothed = 1.2039728043259348
yl = onehotbatch([1], 0:1)
sf = 0.1
yls = [sf (1-sf)]'  # Effective y after label smoothing
ylp = [0.9 0.1]'
logylp = [0.0 v]'

# Construct `sim`ilar and `dis`imilar versions of the dataset so we can test effect of smoothing
# smoothing should decrease loss on disimilar and increase the loss on similar, compared to
# the loss without smoothing
ya = onehotbatch([1, 1, 1, 0, 0], 0:1)
ya_smoothed = label_smoothing(ya, 2sf)
y_same = Float32.(ya)
y_sim = y_same .* (1-2*sf) .+ sf
y_dis = copy(y_sim)
y_dis[1,:], y_dis[2,:] = y_dis[2,:], y_dis[1,:]

@testset "crossentropy" begin
  @test crossentropy([0.1,0.0,0.9], [0.1,0.0,0.9]) ≈ crossentropy([0.1,0.9], [0.1,0.9])
  @test crossentropy(ŷ, y) ≈ lossvalue
  @test crossentropy(ŷ, y_smoothed) ≈ lossvalue_smoothed
  @test crossentropy(ylp, label_smoothing(yl, 2sf)) ≈ -sum(yls.*log.(ylp))
  @test crossentropy(ylp, yl) ≈ -sum(yl.*log.(ylp))
  @test iszero(crossentropy(y_same, ya, ϵ=0))
  @test iszero(crossentropy(ya, ya, ϵ=0))
  @test crossentropy(y_sim, ya) < crossentropy(y_sim, ya_smoothed)
  @test crossentropy(y_dis, ya) > crossentropy(y_dis, ya_smoothed)
end

@testset "logitcrossentropy" begin
  @test logitcrossentropy(logŷ, y) ≈ lossvalue
  @test logitcrossentropy(logylp, yl) ≈ -sum(yl.*logsoftmax(logylp))
  @test logitcrossentropy(logylp, label_smoothing(yl, 2sf)) ≈ -sum(yls.*logsoftmax(logylp))
end

logŷ, y = randn(3), rand(3)
yls = y.*(1-2sf).+sf

@testset "binarycrossentropy" begin
  @test binarycrossentropy.(σ.(logŷ), label_smoothing(y, 2sf; dims=0); ϵ=0) ≈ -yls.*log.(σ.(logŷ)) - (1 .- yls).*log.(1 .- σ.(logŷ))
  @test binarycrossentropy(σ.(logŷ), y; ϵ=0) ≈ mean(-y.*log.(σ.(logŷ)) - (1 .- y).*log.(1 .- σ.(logŷ)))
  @test binarycrossentropy(σ.(logŷ), y) ≈ mean(-y.*log.(σ.(logŷ) .+ eps.(σ.(logŷ))) - (1 .- y).*log.(1 .- σ.(logŷ) .+ eps.(σ.(logŷ))))
end

@testset "logitbinarycrossentropy" begin
  @test logitbinarycrossentropy.(logŷ, label_smoothing(y, 0.2)) ≈ binarycrossentropy.(σ.(logŷ), label_smoothing(y, 0.2); ϵ=0)
  @test logitbinarycrossentropy(logŷ, y) ≈ binarycrossentropy(σ.(logŷ), y; ϵ=0)
end

y = onehotbatch([1], 0:1)
yls = [0.1 0.9]'
@testset "label_smoothing" begin
  @test label_smoothing(y, 0.2) == yls
  @test label_smoothing(y, 0.2; dims=0) == label_smoothing.(y, 0.2; dims=0)
  @test_throws ArgumentError label_smoothing([0., 0., 1., 0.], 1.2)
  @test_throws ArgumentError label_smoothing([0., 0., 1., 0.], 0.)
end

y = [1 2 3]
ŷ = [4.0 5.0 6.0]

@testset "kldivergence" begin
  @test Flux.kldivergence([0.1,0.0,0.9], [0.1,0.0,0.9]) ≈ Flux.kldivergence([0.1,0.9], [0.1,0.9])
  @test Flux.kldivergence(ŷ, y) ≈ -1.7661057888493457
  @test Flux.kldivergence(y, y) ≈ 0
end

y = [1 2 3 4]
ŷ = [5.0 6.0 7.0 8.0]

@testset "hinge_loss" begin
  @test Flux.hinge_loss(ŷ, y) ≈ 0
  @test Flux.hinge_loss(y, 0.5 .* y) ≈ 0.125
end

@testset "squared_hinge_loss" begin
  @test Flux.squared_hinge_loss(ŷ, y) ≈ 0
  @test Flux.squared_hinge_loss(y, 0.5 .* y) ≈ 0.0625
end

y = [0.1 0.2 0.3]
ŷ = [0.4 0.5 0.6]

@testset "poisson_loss" begin
  @test Flux.poisson_loss(ŷ, y) ≈ 0.6278353988097339
  @test Flux.poisson_loss(y, y) ≈ 0.5044459776946685
end

y = [1.0 0.5 0.3 2.4]
ŷ = [0 1.4 0.5 1.2]

@testset "dice_coeff_loss" begin
  @test Flux.dice_coeff_loss(ŷ, y) ≈ 0.2799999999999999
  @test Flux.dice_coeff_loss(y, y) ≈ 0.0
end

@testset "tversky_loss" begin
  @test Flux.tversky_loss(ŷ, y) ≈ -0.06772009029345383
  @test Flux.tversky_loss(ŷ, y, β=0.8) ≈ -0.09490740740740744
  @test Flux.tversky_loss(y, y) ≈ -0.5576923076923075
end

@testset "no spurious promotions" begin
  for T in (Float32, Float64)
    y = rand(T, 2)
    ŷ = rand(T, 2)
    for f in ALL_LOSSES
      fwd, back = Flux.pullback(f, ŷ, y)
      @test fwd isa T
      @test eltype(back(one(T))[1]) == T
    end
  end
end


text:
@testset "1092.shortest-common-supersequence.jl" begin
    @test shortest_common_supersequence("abac", "cab") == "cabac"
    @test shortest_common_supersequence("a", "b") == "ba"
    @test shortest_common_supersequence("cgjmcmrapzvicukazaft", "bjyrpuldxazepiougskd") ==
        "bcgjymcmrpuldxazepzviocugskdazaft"
end


text:
#**************************************************************************************
# Input_Interface.jl
# =============== part of the GeoEfficiency.jl package.
#
# all the input either from the console or from the csv files to the package is handled by some function here.
#
#**************************************************************************************

const datafolder = "GeoEfficiency"
const datadir = joinpath(homedir(), datafolder); 			isdir(datadir) || mkdir(datadir)

const Detectors = "Detectors.csv";
const srcHeights = "srcHeights.csv";
const srcRhos = "srcRhos.csv";
const srcRadii = "srcRadii.csv";
const srcLengths = "srcLengths.csv";

"""
	input(prompt::AbstractString = "? ")
Prompt the user with the massage `prompt` defaults to `? `.
Return a string delimited by new line excluding the new line.
"""
function input(prompt::AbstractString = "? ", incolor::Symbol = :green)
    print_with_color(incolor, prompt)
    chomp(readline())
end # function

"""
	getfloat(prompt::AbstractString = "? ", from::Real = 0.0, to::Real = Inf)

Prompts the user with the massage `prompt` defaults to `? ` to input a numserical expression evaluate to a numerical value and asserts that the value is in the semi open interval [`from`, `to`[.

 > input from the `console` can be numerical expression not just a number.
 >Example:-
 > 5/2, 5//2, pi, e, 1E-2, 5.2/3, sin(1), pi/2/3
 > All are valid expressions.

`Note Please`
\n- a blank (just a return) input is interpreted as being `0.0`.
"""
function getfloat(prompt::AbstractString = "? ",
					from::Real = 0.0,
                    to::Real = Inf)
    value = input(prompt)
	"" == value	&&	return 0.0		# just pressing return is interapted as <0.0>
    try
        val = include_string(value)
		val = float(val)
        @assert from <= val < to
        return val

	catch err
		isa(err, AssertionError) ?
			info("provid a number in the semi open interval [$from, $to[.")	: begin
			info("provid a valid numerical value!")
		end #begin
        return getfloat(prompt, from, to)

    end #try
end	#function

"""

	read_from_csvFile()

read detectors data from predefined file and return its content as an array of detectors.
"""
function read_from_csvFile()
	Detector_info_array::Array{Float64,2} = Array(Float64,(0,0))
	println("INFO: opening '$(Detectors)'......")
	try
		Detector_info_array = readcsv(joinpath(datadir, Detectors),  header=true)[1];

	catch err
		warn("'$(Detectors)' can't be found in '$(datadir)'")
		return getDetectors()

	end #try
	return getDetectors(Detector_info_array)

end #function

"""

	read_from_csvFile(csv_data::AbstractString)

read data from a file and return its content as an array.
`csv_data`: filename of csv file containing data.
"""
function read_from_csvFile(csv_data::AbstractString)
	println("INFO: Opening '$(csv_data)'......")
	try
		return readcsv(joinpath(datadir, csv_data),  header=true)[1][:,1];

	catch err
		warn("'$(csv_data)' can't be found in '$(datadir)'")
		return Float64[0.0]

	end #try
end #function

"""

	read_batch_info()

read `detectors` and `sources` parameters from the predefined csv files.
Return a tuple

	(	Detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		ispoint
		)
"""
function read_batch_info()

	println("INFO: Starting the batch mode.....")
	ispoint = input("\n Is it a point source {Y|n} ? ") |> lowercase != "n"

	function batchfailure()
		input("\t----<( Press return: to treminated batch mode )>----", :blue)
		src = source(isPoint=ispoint)
		srcHeights_array, srcRhos_array = [src[1].Height], [src[1].Rho]
		srcRadii_array, srcLengths_array = [src[2]],[src[3]]
		nothing
	end

	Detectors_array::Array{RadiationDetector,1} = read_from_csvFile()
	srcHeights_array::Array{Float64,1} =  read_from_csvFile(srcHeights) |> sort
	srcRhos_array::Array{Float64,1}
	srcRadii_array::Array{Float64,1}
	srcLengths_array ::Array{Float64,1}

	if srcHeights_array == [0.0]
			batchfailure()

	elseif ispoint
		srcRadii_array  = [0.0]
		srcLengths_array  = [0.0]
		srcRhos_array =	read_from_csvFile(srcRhos) |> sort

	else
		srcRhos_array = [0.0]
		srcRadii_array	 = 	read_from_csvFile(srcRadii) |> sort
		if srcRadii_array == [0.0]
			batchfailure()

		else
			srcLengths_array = 	read_from_csvFile(srcLengths) |> sort

		end #if
	end #if
	println("\n Results log\n=============")
	return (
		Detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		ispoint
		)
end #fumction

"""

	getDetectors()

prompt the user to input detector parameters from the `console`.
Return `Detectors_array` an Array of the entered detectors.
"""
function getDetectors()
	Detectors_array = RadiationDetector[]
	input("----<( Press return: to provid detector specifiction from the console )>----", :blue);
	while(true)
		try
			push!(Detectors_array, detectorFactory())

		catch
			break

		end #try

		res = input("""\n
    	- To add a new detector press return\n
    	- To quit press 'q'|'Q' then return\n
			\n\tyour Choice: """, :blue) |> lowercase;
		res == "q" && break

	end #while
	return Detectors_array
end #function

"""

	getDetectors(Detector_info_array::Array{Float64,2})

convert detectors from the information in `Detector_info_array` and return `Detectors_array` an Array of successfully converted detectors. If the `Detector_info_array` is empty it will call `getDetectors()`.
"""
function getDetectors(Detector_info_array::Array{Float64,2})
	isempty(Detector_info_array) && return getDetectors()
	Detectors_array::Array{RadiationDetector,1} = RadiationDetector[]
	for i_th_line = 1:size(Detector_info_array)[1]
		try
			push!(Detectors_array,
					detectorFactory((Detector_info_array[i_th_line,:])...))

		catch err
			continue

		end #try
	end #for

	return Detectors_array
end #function