text:
"Generic PDDL state description."
mutable struct GenericState <: State
    types::Set{Compound} # Object type declarations
    facts::Set{Term} # Boolean-valued fluents
    values::Dict{Symbol,Any} # All other fluents
end

GenericState(types) = GenericState(types, Set{Term}(), Dict{Symbol,Any}())
GenericState(types, facts) = GenericState(types, facts, Dict{Symbol,Any}())

GenericState(state::GenericState) = copy(state)
function GenericState(state::State)
    types = Set{Compound}((Compound(ty, [o]) for (o, ty) in get_objtypes(state)))
    new = GenericState(types)
    for (term, val) in get_fluents(state)
        if val === false continue end
        set_fluent!(new, val, term)
    end
    return new
end

Base.copy(s::GenericState) =
    GenericState(copy(s.types), copy(s.facts), deepcopy(s.values))
Base.:(==)(s1::GenericState, s2::GenericState) =
    issetequal(s1.types, s2.types) && issetequal(s1.facts, s2.facts) &&
    isequal(s1.values, s2.values)
Base.hash(s::GenericState, h::UInt) =
    hash(s.values, hash(s.facts, hash(s.types, h)))
Base.issubset(s1::GenericState, s2::GenericState) =
    s1.types ⊆ s2.types && s1.facts ⊆ s2.facts

stateindex(domain::GenericDomain, state::GenericState) =
    hash(state)

get_objects(state::GenericState) =
    Const[ty.args[1] for ty in state.types]

get_objects(state::GenericState, type::Symbol) =
    Const[ty.args[1] for ty in state.types if ty.name == type]

get_objtypes(state::GenericState) =
    Dict(ty.args[1] => ty.name for ty in state.types)

get_facts(state::GenericState) =
    state.facts

function get_fluent(state::GenericState, term::Const)
    if term in state.facts
        return true
    else
        return get(state.values, term.name, false)
    end
end

function get_fluent(state::GenericState, term::Compound)
    if term in state.facts
        return true
    else
        d = get(state.values, term.name, nothing)
        return d === nothing ? false : d[Tuple(a.name for a in term.args)]
    end
end

get_fluent(state::GenericState, name::Symbol) =
    get_fluent(state, Const(name))

get_fluent(state::GenericState, name::Symbol, args...) =
    get_fluent(state, Compound(name, collect(args)))

function set_fluent!(state::GenericState, val::Bool, term::Compound)
    if val push!(state.facts, term) else delete!(state.facts, term) end
    return val
end

function set_fluent!(state::GenericState, val::Bool, term::Const)
    if val push!(state.facts, term) else delete!(state.facts, term) end
    return val
end

function set_fluent!(state::GenericState, val::Any, term::Const)
    state.values[term.name] = val
end

function set_fluent!(state::GenericState, val::Any, term::Compound)
    d = get!(state.values, term.name, Dict())
    d[Tuple(a.name for a in term.args)] = val
end

set_fluent!(state::GenericState, val, name::Symbol) =
    set_fluent!(state, val, Const(name))

set_fluent!(state::GenericState, val, name::Symbol, args...) =
    set_fluent!(state, val, Compound(name, collect(args)))

get_fluents(state::GenericState) =
    ((name => get_fluent(state, name)) for name in get_fluent_names(state))

function get_fluent_names(state::GenericState)
    f = x -> begin
        name, val = x
        if isa(val, Dict)
            return (Compound(name, Const.(collect(args))) for args in keys(val))
        else
            return (Const(name),)
        end
    end
    fluent_names = Iterators.flatten(Base.Generator(f, state.values))
    return Iterators.flatten((state.facts, fluent_names))
end

get_fluent_values(state::GenericState) =
    (get_fluent(state, name) for name in get_fluent_names(state))


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using DifferentialEquations: Tsit5, Rosenbrock23

t_start = 0.0
t_stop = 8.0

# case 1: ME-FMU with state events

myFMU = fmiLoad("SpringFrictionPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

solution = fmiSimulateME(fmuStruct, t_start, t_stop)
@test length(solution.states.u) > 0
@test length(solution.states.t) > 0

@test solution.states.t[1] == t_start
@test solution.states.t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test solution.states.u[1] == [0.5, 0.0]
@test sum(abs.(solution.states.u[end] - [1.06736, -1.03552e-10])) < 0.1
fmiUnload(myFMU)

# case 2: ME-FMU with state and time events

myFMU = fmiLoad("SpringTimeFrictionPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

### test without recording values

solution = fmiSimulateME(fmuStruct, t_start, t_stop; dtmax=0.001) # dtmax to force resolution
@test length(solution.states.u) > 0
@test length(solution.states.t) > 0

@test solution.states.t[1] == t_start
@test solution.states.t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test solution.states.u[1] == [0.5, 0.0]
@test sum(abs.(solution.states.u[end] - [1.05444, 1e-10])) < 0.01

### test with recording values (variable step record values)

solution= fmiSimulateME(fmuStruct, t_start, t_stop; recordValues="mass.f", dtmax=0.001) # dtmax to force resolution
dataLength = length(solution.states.u)
@test dataLength > 0
@test length(solution.states.t) == dataLength
@test length(solution.values.saveval) == dataLength
@test length(solution.values.t) == dataLength

@test solution.states.t[1] == t_start
@test solution.states.t[end] == t_stop
@test solution.values.t[1] == t_start
@test solution.values.t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test sum(abs.(solution.states.u[1] - [0.5, 0.0])) < 1e-4
@test sum(abs.(solution.states.u[end] - [1.05444, 1e-10])) < 0.01
@test abs(solution.values.saveval[1][1] - 0.75) < 1e-4
@test sum(abs.(solution.values.saveval[end][1] - -0.54435 )) < 0.015

### test with recording values (fixed step record values)

tData = t_start:0.1:t_stop
solution = fmiSimulateME(fmuStruct, t_start, t_stop; recordValues="mass.f", saveat=tData, dtmax=0.001) # dtmax to force resolution
@test length(solution.states.u) == length(tData)
@test length(solution.states.t) == length(tData)
@test length(solution.values.saveval) == length(tData)
@test length(solution.values.t) == length(tData)

@test solution.states.t[1] == t_start
@test solution.states.t[end] == t_stop
@test solution.values.t[1] == t_start
@test solution.values.t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test sum(abs.(solution.states.u[1] - [0.5, 0.0])) < 1e-4
@test sum(abs.(solution.states.u[end] - [1.05444, 1e-10])) < 0.01
@test abs(solution.values.saveval[1][1] - 0.75) < 1e-4
@test sum(abs.(solution.values.saveval[end][1] - -0.54435 )) < 0.015

fmiUnload(myFMU)

# case 3a: ME-FMU without events, but with input signal (explicit solver: Tsit5)

function extForce(t)
    [sin(t)]
end

myFMU = fmiLoad("SpringPendulumExtForce1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

solution = fmiSimulateME(fmuStruct, t_start, t_stop; inputValueReferences=["extForce"], inputFunction=extForce, solver=Tsit5(), dtmax=0.001) # dtmax to force resolution
@test length(solution.states.u) > 0
@test length(solution.states.t) > 0

@test solution.states.t[1] == t_start
@test solution.states.t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test solution.states.u[1] == [0.5, 0.0]
@test sum(abs.(solution.states.u[end] - [0.613371, 0.188633])) < 0.012
fmiUnload(myFMU)

# case 3b: ME-FMU without events, but with input signal (implicit solver: Rosenbrock23, autodiff)

myFMU = fmiLoad("SpringPendulumExtForce1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

# ToDo: autodiff=true not working currently!
# solution = fmiSimulateME(fmuStruct, t_start, t_stop; inputValueReferences=["extForce"], inputFunction=extForce, solver=Rosenbrock23(autodiff=true), dtmax=0.001) # dtmax to force resolution
# @test length(solution.states.u) > 0
# @test length(solution.states.t) > 0

# @test solution.states.t[1] == t_start
# @test solution.states.t[end] == t_stop

# # reference values from Simulation in Dymola2020x (Dassl)
# @test solution.states.u[1] == [0.5, 0.0]
# @test sum(abs.(solution.states.u[end] - [0.613371, 0.188633])) < 0.01
fmiUnload(myFMU)

# case 3c: ME-FMU without events, but with input signal (implicit solver: Rosenbrock23, no autodiff)

myFMU = fmiLoad("SpringPendulumExtForce1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

solution = fmiSimulateME(fmuStruct, t_start, t_stop; inputValueReferences=["extForce"], inputFunction=extForce, solver=Rosenbrock23(autodiff=false), dtmax=0.001) # dtmax to force resolution
@test length(solution.states.u) > 0
@test length(solution.states.t) > 0

@test solution.states.t[1] == t_start
@test solution.states.t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test solution.states.u[1] == [0.5, 0.0]
@test sum(abs.(solution.states.u[end] - [0.613371, 0.188633])) < 0.01
fmiUnload(myFMU)

# case 4: ME-FMU without events, but saving value interpolation

myFMU = fmiLoad("SpringPendulumExtForce1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

solution = fmiSimulateME(fmuStruct, t_start, t_stop; saveat=tData, recordValues=myFMU.modelDescription.stateValueReferences)
@test length(solution.states.u) == length(tData)
@test length(solution.states.t) == length(tData)
@test length(solution.values.saveval) == length(tData)
@test length(solution.values.t) == length(tData)

for i in 1:length(tData)
    @test sum(abs(solution.states.t[i] - solution.states.t[i])) < 1e-6
    @test sum(abs(solution.states.u[i][1] - solution.values.saveval[i][1])) < 1e-6
    @test sum(abs(solution.states.u[i][2] - solution.values.saveval[i][2])) < 1e-6
end

fmiUnload(myFMU)

# case 5: ME-FMU with different (random) start state

myFMU = fmiLoad("SpringFrictionPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct != nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

rand_x0 = rand(2)
solution = fmiSimulateME(fmuStruct, t_start, t_stop; x0=rand_x0)
@test length(solution.states.u) > 0
@test length(solution.states.t) > 0

@test solution.states.t[1] == t_start
@test solution.states.t[end] == t_stop

@test solution.states.u[1] == rand_x0
fmiUnload(myFMU)

ext: jl


text:
module Metis

using SparseArrays
using LinearAlgebra
import LightGraphs
import SimpleWeightedGraphs
using METIS_jll: libmetis

# Metis C API
include("metis_h.jl")
const options = fill(Cint(-1), METIS_NOPTIONS)
options[METIS_OPTION_NUMBERING] = 1

# Julia interface
"""
    Metis.Graph

1-based CSR representation of a graph as defined in
section 5.5 "Graph data structure" in the Metis manual.
"""
struct Graph
    nvtxs::idx_t
    xadj::Vector{idx_t}
    adjncy::Vector{idx_t}
    vwgt::Vector{idx_t}
    adjwgt::Vector{idx_t}
    Graph(nvtxs, xadj, adjncy) = new(nvtxs, xadj, adjncy)
    Graph(nvtxs, xadj, adjncy, vwgt) = new(nvtxs, xadj, adjncy, vwgt)
    Graph(nvtxs, xadj, adjncy, vwgt, adjwgt) = new(nvtxs, xadj, adjncy, vwgt, adjwgt)
end

"""
    Metis.graph(G::SparseMatrixCSC; check_hermitian=true)

Construct the 1-based CSR representation of the sparse matrix `G`.
If `check_hermitian` is `false` the matrix is not checked for being hermitian
before constructing the graph.
"""
function graph(G::SparseMatrixCSC; check_hermitian=true)
    if check_hermitian
        ishermitian(G) || throw(ArgumentError("matrix must be Hermitian"))
    end
    N = size(G, 1)
    xadj = Vector{idx_t}(undef, N+1)
    xadj[1] = 1
    adjncy = Vector{idx_t}(undef, nnz(G))
    adjncy_i = 0
    @inbounds for j in 1:N
        n_rows = 0
        for k in G.colptr[j] : (G.colptr[j+1] - 1)
            i = G.rowval[k]
            if i != j # don't include diagonal elements
                n_rows += 1
                adjncy_i += 1
                adjncy[adjncy_i] = i
            end
        end
        xadj[j+1] = xadj[j] + n_rows
    end
    resize!(adjncy, adjncy_i)
    return Graph(idx_t(N), xadj, adjncy)
end

"""
    graph(G::LightGraphs.AbstractSimpleGraph)

Construct the 1-based CSR representation of the `LightGraphs` graph `G`.
"""
function graph(G::LightGraphs.AbstractSimpleGraph)
    N = LightGraphs.nv(G)
    xadj = Vector{idx_t}(undef, N+1)
    xadj[1] = 1
    adjncy = Vector{idx_t}(undef, 2*LightGraphs.ne(G))
    adjncy_i = 0
    for j in 1:N
        ne = 0
        for i in LightGraphs.outneighbors(G, j)
            ne += 1
            adjncy_i += 1
            adjncy[adjncy_i] = i
        end
        xadj[j+1] = xadj[j] + ne
    end
    resize!(adjncy, adjncy_i)
    return Graph(idx_t(N), xadj, adjncy)
end

"""
    graph(G::SimpleWeightedGraphs.AbstractSimpleWeightedGraph)

Construct the 1-based CSR representation of the `SimpleWeightedGraphs` graph `G`.
The edge weights should be positive integers.
"""
function graph(G::SimpleWeightedGraphs.AbstractSimpleWeightedGraph)
    N = SimpleWeightedGraphs.nv(G)
    xadj = Vector{idx_t}(undef, N+1)
    xadj[1] = 1
    adjncy = Vector{idx_t}(undef, 2*SimpleWeightedGraphs.ne(G))
    vwgt = ones(idx_t, N)
    adjwgt = Vector{idx_t}(undef, 2*SimpleWeightedGraphs.ne(G))
    adjncy_i = 0
    for j in 1:N
        ne = 0
        for i in SimpleWeightedGraphs.outneighbors(G, j)
            ne += 1
            adjncy_i += 1
            adjncy[adjncy_i] = i
            edgeweight = SimpleWeightedGraphs.get_weight(G, j, i)
            if isinteger(edgeweight) && edgeweight > 0
                adjwgt[adjncy_i] = idx_t(edgeweight)
            else
                throw(ErrorException("edge weight should be a positive integer (edge ($j, $i) with weight $edgeweight)"))
            end
        end
        xadj[j+1] = xadj[j] + ne
    end
    resize!(adjncy, adjncy_i)
    return Graph(idx_t(N), xadj, adjncy, vwgt, adjwgt)
end

"""
    perm, iperm = Metis.permutation(G)

Compute the fill reducing permutation `perm`
and its inverse `iperm` of `G`.
"""
permutation(G) = permutation(graph(G))

function permutation(G::Graph)
    perm = Vector{idx_t}(undef, G.nvtxs)
    iperm = Vector{idx_t}(undef, G.nvtxs)
    vwgt = isdefined(G, :vwgt) ? G.vwgt : C_NULL
    METIS_NodeND(G.nvtxs, G.xadj, G.adjncy, vwgt, options, perm, iperm)
    return perm, iperm
end

"""
    Metis.partition(G, n; alg = :KWAY)

Partition the graph `G` in `n` parts.
The partition algorithm is defined by the `alg` keyword:
 - :KWAY: multilevel k-way partitioning
 - :RECURSIVE: multilevel recursive bisection
"""
partition(G, nparts; alg = :KWAY) = partition(graph(G), nparts, alg = alg)

function partition(G::Graph, nparts::Integer; alg = :KWAY)
    part = Vector{idx_t}(undef, G.nvtxs)
    vwgt = isdefined(G, :vwgt) ? G.vwgt : C_NULL
    adjwgt = isdefined(G, :adjwgt) ? G.adjwgt : C_NULL
    edgecut = fill(idx_t(0), 1)
    if alg === :RECURSIVE
        METIS_PartGraphRecursive(G.nvtxs, idx_t(1), G.xadj, G.adjncy, vwgt, C_NULL, adjwgt,
                                 idx_t(nparts), C_NULL, C_NULL, options, edgecut, part)
    elseif alg === :KWAY
        METIS_PartGraphKway(G.nvtxs, idx_t(1), G.xadj, G.adjncy, vwgt, C_NULL, adjwgt,
                            idx_t(nparts), C_NULL, C_NULL, options, edgecut, part)
    else
        throw(ArgumentError("unknown algorithm $(repr(alg))"))
    end
    return part
end

"""
    Metis.separator(G)

Compute a vertex separator of the graph `G`.
"""
separator(G) = separator(graph(G))

function separator(G::Graph)
    part = Vector{idx_t}(undef, G.nvtxs)
    sepsize = fill(idx_t(0), 1)
    vwgt = isdefined(G, :vwgt) ? G.vwgt : C_NULL
    # METIS_ComputeVertexSeparator segfaults with 1-based indexing
    xadj = G.xadj .- idx_t(1)
    adjncy = G.adjncy .- idx_t(1)
    METIS_ComputeVertexSeparator(G.nvtxs, xadj, adjncy, vwgt, options, sepsize, part)
    part .+= 1
    return part
end

end # module


text:
module Telegraph

using Reexport

@reexport using Revise
#@reexport using RandomWalker
#@reexport using DataFrames
#@reexport using Random
#@reexport using Distributions
#@reexport using Chain

export
	telegraph_points,
	telegraph_time_dist,
	update,
	get_all_indices,
	get_indices,
	get_times

include("telegraph_structs.jl")
include("update.jl")
include("times_indices.jl")

end


text:
include("Ndq.jl")
include("unique_nodes.jl")
include("combvec.jl")
include("quadtest.jl")


text:
module Scripting

using DrWatson

# Conventions
#   args        ::  NamedTuple  :   (x=1, y=2)
#   argdict     ::  Dict        :   Dict(:x => 1, :y => 2)
#   savename    ::  String      :   "x=1-y=2"
#   cmdstr      ::  String      :   "--x 1 --y 2"

function overwrite(args, override)
    argdict = ntuple2dict(args)
    for k in keys(override)
        if k in keys(argdict)
            v = override[k]
            if argdict[k] == v
                @warn "The values for key :$k are the same: $v; overwriting is skipped."
            else
                @info "Overriding :$k: $(args[k]) => $v"
            end
            argdict[k] = v
        else
            @warn "Key $k doesn't exist in `args`; overwriting is skipped."
        end
    end
    return dict2ntuple(argdict)
end

cmdstr(args::NamedTuple) = cmdstr(ntuple2dict(args))
function cmdstr(argdict::Dict)
    return replace(savename(Dict("--$k" => v for (k, v) in argdict); connector=" "), "=" => " ")
end

import Dates

const DATETIME_FMT = "ddmmyyyy-H-M-S"

function find_latestdir(parentdir)
    dates = []
    for dirstr in readdir(parentdir)
        try
            push!(dates, Dates.DateTime(dirstr, DATETIME_FMT))
        catch
            nothing
        end
    end
    return Dates.format(maximum(dates), DATETIME_FMT)
end

using Pkg.TOML

function parse_toml(tomlpath::String, tableinfo::Tuple{Vararg{Pair}})
    # Load TOML file into nested dictionary
    toml = TOML.parsefile(tomlpath)
    # Get the "common" section
    argdict = toml["common"]
    local table = toml
    for (k, v) in tableinfo
        table = table[v]
        argdict = merge(argdict, filter(p -> !(p.second isa Dict), table))
        argdict[string(k)] = v
    end
    # Check if there are sub-tables left to read
    @assert length(filter(p -> p.second isa Dict, table)) == 0 "There sub-table left in the current nesting table."
    # Convert keys to Symbol
    return Dict(Symbol(p.first) => p.second for p in argdict)
end

import ArgParse

function parse_argstr(argstr::AbstractString, settings; as_symbols::Bool=false)
    return ArgParse.parse_args(split(replace(argstr, r"\s+" => " "), " "), settings; as_symbols=as_symbols)
end

export overwrite, cmdstr, DATETIME_FMT, find_latestdir, parse_toml, parse_argstr

include("check.jl")
export @jupyter, @script, @tb, @wb, @checknumerics

end # module


text:
# operation on MATLAB engine sessions

###########################################################
#
#   Session open & close
#
###########################################################

# 64 K buffer should be sufficient to store the output text in most cases
const default_output_buffer_size = 64 * 1024

type MSession
    ptr::Ptr{Void}
    buffer::Vector{Uint8}
    bufptr::Ptr{Uint8}

    function MSession(bufsize::Integer)
        global libeng
        if libeng == C_NULL
            load_libeng()
        end
        @assert libeng != C_NULL

        ep = ccall(engfunc(:engOpen), Ptr{Void}, (Ptr{Uint8},), default_startcmd)
        if ep == C_NULL
            throw(MEngineError("Failed to open a MATLAB engine session."))
        end

        buf = Array(Uint8, bufsize)

        if bufsize > 0
            bufptr = pointer(buf)
            ccall(engfunc(:engOutputBuffer), Cint, (Ptr{Void}, Ptr{Uint8}, Cint),
                ep, bufptr, bufsize)
        else
            bufptr = convert(Ptr{Uint8}, C_NULL)
        end

        if OS_NAME == :Windows
            # Hide the MATLAB Command Window on Windows
            ccall(engfunc(:engSetVisible ), Cint, (Ptr{Void}, Cint), ep, 0)
        end

        println("A MATLAB session is open successfully")
        new(ep, buf, bufptr)
    end

    MSession() = MSession(default_output_buffer_size)
end

function close(session::MSession)
    # Close a MATLAB Engine session

    @assert libeng::Ptr{Void} != C_NULL

    r = ccall(engfunc(:engClose), Cint, (Ptr{Void},), session.ptr)
    if r != 0
        throw(MEngineError("Failed to close a MATLAB engine session (err = $r)"))
    end
end

# default session

default_msession = nothing

function restart_default_msession(bufsize::Integer)
    global default_msession
    if !(default_msession == nothing)
        close(default_msession)
    end
    default_msession = MSession(bufsize)
end

restart_default_msession() = restart_default_msession(default_output_buffer_size)

function get_default_msession()
    global default_msession
    if default_msession == nothing
        default_msession = MSession()
    end
    default_msession::MSession
end

function close_default_msession()
    global default_msession
    if !(default_msession == nothing)
        close(default_msession)
        default_msession = nothing
    end
end

###########################################################
#
#   communication with MATLAB session
#
###########################################################

function eval_string(session::MSession, stmt::ASCIIString)
    # Evaluate a MATLAB statement in a given MATLAB session

    @assert libeng::Ptr{Void} != C_NULL

    r::Cint = ccall(engfunc(:engEvalString), Cint,
        (Ptr{Void}, Ptr{Uint8}), session.ptr, stmt)

    if r != 0
        throw(MEngineError("Invalid engine session."))
    end

    bufptr::Ptr{Uint8} = session.bufptr
    if bufptr != C_NULL
        bs = bytestring(bufptr)
        if ~isempty(bs)
            print(bs)
        end
    end
end

eval_string(stmt::ASCIIString) = eval_string(get_default_msession(), stmt)

function put_variable(session::MSession, name::Symbol, v::MxArray)
    # Put a variable into a MATLAB engine session

    @assert libeng::Ptr{Void} != C_NULL

    r = ccall(engfunc(:engPutVariable), Cint,
        (Ptr{Void}, Ptr{Uint8}, Ptr{Void}), session.ptr, string(name), v.ptr)

    if r != 0
        throw(MEngineError("Failed to put the variable $(name) into a MATLAB session."))
    end
end

put_variable(session::MSession, name::Symbol, v) = put_variable(session, name, mxarray(v))

put_variable(name::Symbol, v) = put_variable(get_default_msession(), name, v)

function get_mvariable(session::MSession, name::Symbol)

    @assert libeng::Ptr{Void} != C_NULL

    pv = ccall(engfunc(:engGetVariable), Ptr{Void},
        (Ptr{Void}, Ptr{Uint8}), session.ptr, string(name))

    if pv == C_NULL
        throw(MEngineError("Failed to get the variable $(name) from a MATLAB session."))
    end
    MxArray(pv)
end

get_mvariable(name::Symbol) = get_mvariable(get_default_msession(), name)

get_variable(name::Symbol) = jvariable(get_mvariable(name))
get_variable(name::Symbol, kind) = jvariable(get_mvariable(name), kind)

###########################################################
#
#   macro to simplify syntax
#
###########################################################

function _mput_multi(vs::Symbol...)
    nv = length(vs)
    if nv == 1
        v = vs[1]
        :( MATLAB.put_variable($(Meta.quot(v)), $(v)) )
    else
        stmts = Array(Expr, nv)
        for i = 1 : nv
            v = vs[i]
            stmts[i] = :( MATLAB.put_variable($(Meta.quot(v)), $(v)) )
        end
        Expr(:block, stmts...)
    end
end

macro mput(vs...)
    esc( _mput_multi(vs...) )
end

function make_getvar_statement(v::Symbol)
    :( $(v) = MATLAB.get_variable($(Meta.quot(v))) )
end

function make_getvar_statement(ex::Expr)
    if !(ex.head == :(::))
        error("Invalid expression for @mget.")
    end
    v::Symbol = ex.args[1]
    k::Symbol = ex.args[2]

    :( $(v) = MATLAB.get_variable($(Meta.quot(v)), $(k)) )
end

function _mget_multi(vs::Union(Symbol, Expr)...)
    nv = length(vs)
    if nv == 1
        make_getvar_statement(vs[1])
    else
        stmts = Array(Expr, nv)
        for i = 1 : nv
            stmts[i] = make_getvar_statement(vs[i])
        end
        Expr(:block, stmts...)
    end
end

macro mget(vs...)
    esc( _mget_multi(vs...) )
end

macro matlab(ex)
    :( MATLAB.eval_string($(mstatement(ex))) )
end

###########################################################
#
#   mxcall
#
###########################################################

# MATLAB does not allow underscore as prefix of a variable name
_gen_marg_name(mfun::Symbol, prefix::ASCIIString, i::Int) = "jx_$(mfun)_arg_$(prefix)_$(i)"

function mxcall(session::MSession, mfun::Symbol, nout::Integer, in_args...)
    nin = length(in_args)

    # generate tempoary variable names

    in_arg_names = Array(ASCIIString, nin)
    out_arg_names = Array(ASCIIString, nout)

    for i = 1 : nin
        in_arg_names[i] = _gen_marg_name(mfun, "in", i)
    end

    for i = 1 : nout
        out_arg_names[i] = _gen_marg_name(mfun, "out", i)
    end

    # generate MATLAB statement

    buf = IOBuffer()
    if nout > 0
        if nout > 1
            print(buf, "[")
        end
        print(buf, join(out_arg_names, ", "))
        if nout > 1
            print(buf, "]")
        end
        print(buf, " = ")
    end

    print(buf, string(mfun))
    print(buf, "(")
    if nin > 0
        print(buf, join(in_arg_names, ", "))
    end
    print(buf, ");")

    stmt = bytestring(buf)

    # put variables to MATLAB

    for i = 1 : nin
        put_variable(session, symbol(in_arg_names[i]), in_args[i])
    end

    # execute MATLAB statement

    eval_string(session, stmt)

    # get results from MATLAB

    ret = if nout == 1
        jvariable(get_mvariable(session, symbol(out_arg_names[1])))
    elseif nout >= 2
        results = Array(Any, nout)
        for i = 1 : nout
            results[i] = jvariable(get_mvariable(session, symbol(out_arg_names[i])))
        end
        tuple(results...)
    else
        nothing
    end

    # clear temporaries from MATLAB workspace

    for i = 1 : nin
        eval_string(session, string("clear ", in_arg_names[i], ";"))
    end

    for i = 1 : nout
        eval_string(session, string("clear ", out_arg_names[i], ";"))
    end

    # return
    ret
end

mxcall(mfun::Symbol, nout::Integer, in_args...) = mxcall(get_default_msession(), mfun, nout, in_args...)


text:
using Distributed;

# add two local workers
addprocs(2)

# add TCP/IP workers (two per node)
addprocs([("grid61", 2)])
addprocs([("grid62", 2)])
addprocs([("grid63", 2)])

println(workers())

println("distributed workers initialised, starting tests")

@everywhere function test()
    id = myid()
    println(id)

    try
        open("/data/cache/$(id).txt", "w") do io
            println(io, "myid() = $id")
        end
    catch err
        println(err)
    end
end

ras = [@spawnat w test() for w in workers()]

println("ras: ", ras)

@time wait.(ras)


text:
module TestLift
    using Lift
    using Base.Test

    f(x::Int) = 5 * x
    g(x::Int, y::Int) = x + y
    x = Nullable(5)
    y = Nullable{Int}()
    X = NullableArray([1, 2, 3, 4, 5])

    @test isequal(@^(f(x), Int), Nullable(25))
    @test isequal(@^(f(y), Int), Nullable{Int}())
    @test_throws ArgumentError eval(macroexpand(:( @^f(y) )))
    @test_throws ArgumentError eval(macroexpand(:( @^ if x > 0; f(x) end )))

    @test isequal(@^(f(X[1]) + g(x, X[1]), Int), Nullable(11))
    @test isequal(@^(f(X[1]) + g(y, X[1]), Int), Nullable{Int}())
end


text:
export bpffunction, bpfasm

# Compilation support

Base.@kwdef struct BPFCompilerParams <: AbstractCompilerParams
    format::Symbol = :obj
    license::String = ""
    prog_section::String = "prog"
    btf::Bool = true
end

BPFCompilerJob = CompilerJob{BPFCompilerTarget,BPFCompilerParams}

GPUCompiler.runtime_module(::BPFCompilerJob) = BPFnative

## host-side API

"""
    bpffunction(f, tt=Tuple{}; kwargs...)

Compiles a function invocation into eBPF bytecode.

The following keyword arguments are supported:
- `name`: override the name that the kernel will have in the generated code
- `license::String=""`: license for the kernel source code and resulting object
- `prog_section::String=""`: ELF section that the kernel will be placed in
- `btf::Bool=true`: Whether to generate BTF debuginfo or not

The output of this function is automatically cached, i.e. you can simply call
`bpffunction` in a hot path without degrading performance. New code will be
generated automatically when the function changes or when different types or
keyword arguments are provided.
"""
function bpffunction(f::Core.Function, tt::Type=Tuple{}; name=nothing, kwargs...)
    source = FunctionSpec(f, tt, false, name)
    target = BPFCompilerTarget()
    params = BPFCompilerParams(; kwargs...)
    job = CompilerJob(target, source, params)
    GPUCompiler.cached_compilation(bpffunction_cache, job,
                                   bpffunction_compile,
                                   bpffunction_link)
end

const bpffunction_cache = Dict{UInt,Any}()

# actual compilation
function bpffunction_compile(@nospecialize(job::CompilerJob))
    params = job.params
    format = params.format
    btf = params.btf

    # compile to BPF
    method_instance, world = GPUCompiler.emit_julia(job)
    ir, kernel = GPUCompiler.emit_llvm(job, method_instance; libraries=false)

    if format == :llvm
        return collect.(codeunits.(string.((ir, kernel)))) # TODO: Make more efficient
    elseif format == :asm
        format = LLVM.API.LLVMAssemblyFile
    elseif format == :obj
        format = LLVM.API.LLVMObjectFile
    end
    code, _ = GPUCompiler.emit_asm(job, ir; format=format, validate=true, strip=!btf)
    return Vector{UInt8}(code)
end
bpffunction_link(@nospecialize(job::CompilerJob), exe) = exe

function GPUCompiler.optimize_module!(job::BPFCompilerJob, mod::LLVM.Module)
    # extra optimizations to kill uses of `gpu_signal_exception`
    # GPUCompiler does this slightly too late, causing an LLVM abort due to the
    # BPF target not supporting lowering of `llvm.trap`
    ModulePassManager() do pm
        global_optimizer!(pm)
        global_dce!(pm)
        strip_dead_prototypes!(pm)
        run!(pm, mod)
    end
    for func in LLVM.functions(mod)
        # validate no-throw
        if LLVM.name(func) == "gpu_signal_exception" && length(collect(LLVM.uses(func))) > 0
            throw(GPUCompiler.KernelError(job, "eBPF does not support exceptions"))
        end
    end
    mod
end
function GPUCompiler.finish_module!(job::BPFCompilerJob, mod::LLVM.Module, entry::LLVM.Function)
    params = job.params
    license = params.license
    prog_section = params.prog_section

    for func in LLVM.functions(mod)
        # Set entry section for loaders like libbpf
        LLVM.section!(func, prog_section)
    end

    # Set license
    if license != ""
        ctx = LLVM.context(mod)
        i8 = LLVM.Int8Type(ctx)
        glob = GlobalVariable(mod, LLVM.ArrayType(i8, length(license)+1), "_license")
        linkage!(glob, LLVM.API.LLVMExternalLinkage)
        constant!(glob, true)
        section!(glob, "license")
        str = ConstantArray(Vector{UInt8}(license*'\0'); ctx)
        @assert context(glob) == context(str) == ctx
        initializer!(glob, str)
    end

    # Set all map definitions as external linkage
    for gv in filter(x->(section(x)=="maps")||(section(x)==".maps"), collect(LLVM.globals(mod)))
        linkage!(gv, LLVM.API.LLVMExternalLinkage)
    end

    ModulePassManager() do pm
        if Base.JLOptions().debug_level > 1
            # Validate contexts, for my sanity
            add!(pm, ModulePass("BPFValidateContexts", validate_contexts!))
        end
        # Promote `@malloc` intrinsics
        add!(pm, FunctionPass("BPFHeapToStack", heap_to_stack!))
        run!(pm, mod)
    end

    entry
end

"Validates LLVM contexts of all the things."
function validate_contexts!(mod::LLVM.Module)
    ctx = LLVM.context(mod)
    for fn in LLVM.functions(mod)
        @assert context(fn) == ctx "Failed validation: $fn"
        for bb in LLVM.blocks(fn)
            for insn in LLVM.instructions(bb)
                @assert context(insn) == ctx "Failed validation: $insn"
                for op in LLVM.operands(insn)
                    @assert context(op) == ctx "Failed validation: $op"
                end
            end
        end
    end
    for gv in LLVM.globals(mod)
        @assert context(gv) == ctx "Failed validation: $gv"
    end
    false
end

"Promotes `@malloc` intrinsics to allocas."
function heap_to_stack!(fn::LLVM.Function)
    changed = false
    ctx = LLVM.context(fn)
    for bb in LLVM.blocks(fn)
        for insn in LLVM.instructions(bb)
            if insn isa LLVM.CallInst && LLVM.name(LLVM.called_value(insn)) == "malloc"
                sz = convert(Int64, LLVM.operands(insn)[1])
                T_i8 = LLVM.Int8Type(ctx)
                T_pi8 = LLVM.PointerType(T_i8)
                T_buf = LLVM.ArrayType(T_i8, sz)
                Builder(ctx) do builder
                    # Place alloca at beginning of entry
                    position!(builder, first(LLVM.instructions(first(LLVM.blocks(fn)))))
                    buf = alloca!(builder, T_buf)
                    # Replace malloc with bitcast'd alloca
                    position!(builder, insn)
                    new_insn = bitcast!(builder, buf, T_pi8)
                    replace_uses!(insn, new_insn)
                    unsafe_delete!(LLVM.parent(insn), insn)
                end
                changed = true
            end
        end
    end
    changed
end


text:
# Custom definitions of LinearIndices and CartesianIndices to take into account
# index permutations.
#
# In particular, when array dimensions are permuted, the default
# CartesianIndices do not iterate in memory order, making them suboptimal.
# We try to workaround that by adding a custom definition of CartesianIndices.
#
# (TODO Better / cleaner way to do this??)

# We make LinearIndices(::PencilArray) return a PermutedLinearIndices, which
# takes index permutation into account.
struct PermutedLinearIndices{
        N, L <: LinearIndices, Perm,
        Offsets <: Union{Nothing, Dims{N}},
    } <: AbstractArray{Int,N}
    data :: L  # indices in permuted order
    perm :: Perm
    offsets :: Offsets
    function PermutedLinearIndices(ind::LinearIndices{N},
                                   perm::Perm, offsets=nothing) where {N, Perm}
        L = typeof(ind)
        Off = typeof(offsets)
        new{N, L, Perm, Off}(ind, perm, offsets)
    end
end

Base.length(L::PermutedLinearIndices) = length(L.data)
Base.size(L::PermutedLinearIndices) = L.perm \ size(L.data)
Base.iterate(L::PermutedLinearIndices, args...) = iterate(L.data, args...)
Base.lastindex(L::PermutedLinearIndices) = lastindex(L.data)

@inline _apply_offset(I::CartesianIndex, ::Nothing) = I
@inline _apply_offset(I::CartesianIndex{N}, off::Dims{N}) where {N} =
    CartesianIndex(Tuple(I) .- off)

@inline function Base.getindex(L::PermutedLinearIndices, i::Integer)
    @boundscheck checkbounds(L.data, i)
    @inbounds L.data[i]
end

@inline function Base.getindex(
        L::PermutedLinearIndices{N}, I::CartesianIndex{N}) where {N}
    Ioff = _apply_offset(I, L.offsets)
    J = L.perm * Ioff
    @boundscheck checkbounds(L.data, J)
    @inbounds L.data[J]
end

Base.LinearIndices(A::PencilArray) =
    PermutedLinearIndices(LinearIndices(parent(A)), permutation(A))

function Base.LinearIndices(g::GlobalPencilArray)
    A = parent(g)
    off = g.offsets
    PermutedLinearIndices(LinearIndices(parent(A)), permutation(A), off)
end

# We make CartesianIndices(::PencilArray) return a PermutedCartesianIndices,
# which loops faster (in memory order) when there are index permutations.
struct PermutedCartesianIndices{
        N, C <: CartesianIndices{N}, Iperm,
        Offsets} <: AbstractArray{CartesianIndex{N}, N}
    data  :: C  # indices in permuted order
    iperm :: Iperm  # inverse permutation
    offsets :: Offsets
    function PermutedCartesianIndices(ind::CartesianIndices{N},
                                      perm::Perm, offsets=nothing) where {N, Perm}
        iperm = inv(perm)
        C = typeof(ind)
        Iperm = typeof(iperm)
        Off = typeof(offsets)
        new{N, C, Iperm, Off}(ind, iperm, offsets)
    end
end

Base.size(C::PermutedCartesianIndices) = C.iperm * size(C.data)

@inline function Base.iterate(C::PermutedCartesianIndices, args...)
    next = iterate(C.data, args...)
    next === nothing && return nothing
    I, state = next  # `I` has permuted indices
    J = C.iperm * I  # unpermute indices
    Joff = _apply_offset(J, C.offsets)
    Joff, state
end

# Get i-th Cartesian index in memory (permuted) order.
# Returns the Cartesian index in logical (unpermuted) order.
@inline function Base.getindex(
        C::PermutedCartesianIndices, i::Integer)
    @boundscheck checkbounds(C.data, i)
    @inbounds I = C.data[i]  # convert linear to Cartesian index (relatively slow...)
    J = C.iperm * I          # unpermute indices
    Joff = _apply_offset(J, C.offsets)
    Joff
end

Base.CartesianIndices(A::PencilArray) =
    PermutedCartesianIndices(CartesianIndices(parent(A)), permutation(A))

function Base.CartesianIndices(g::GlobalPencilArray)
    A = parent(g)
    off = g.offsets .* -1  # negative offsets
    PermutedCartesianIndices(CartesianIndices(parent(A)), permutation(A), off)
end


text:
# Timescale for diffusion across one cell
min_Δxyz(grid) = min(grid.Δx, grid.Δy, grid.Δz)
min_Δxy(grid) = min(grid.Δx, grid.Δy)
min_Δz(grid) = grid.Δz

cell_diffusion_timescale(model) = cell_diffusion_timescale(model.closure, model.diffusivities, model.grid)

function cell_diffusion_timescale(closure::ConstantIsotropicDiffusivity{V, <:NamedTuple{()}},
                                  diffusivities, grid) where V
    return min_Δxyz(grid)^2 / closure.ν
end

function cell_diffusion_timescale(closure::ConstantIsotropicDiffusivity, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_κ = maximum(closure.κ)
    return min(Δ^2 / closure.ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::ConstantAnisotropicDiffusivity{V, <:NamedTuple{()}, <:NamedTuple{()}},
                                  diffusivities, grid) where V
    Δh = min_Δxy(grid)
    Δz = min_Δz(grid)
    return min(Δh^2 / closure.νh, Δz^2 / closure.νv)
end

function cell_diffusion_timescale(closure::ConstantAnisotropicDiffusivity, diffusivities, grid)
    Δh = min_Δxy(grid)
    Δz = min_Δz(grid)
    max_κh = maximum(closure.κh)
    max_κv = maximum(closure.κv)
    return min(Δh^2 / closure.νh, Δz^2 / closure.νv,
               Δh^2 / max_κh, Δz^2 / max_κv)
end

function cell_diffusion_timescale(closure::AnisotropicBiharmonicDiffusivity{V, <:NamedTuple{()}, <:NamedTuple{()}},
                                  diffusivities, grid) where V
    Δh = min_Δxy(grid)
    Δz = min_Δz(grid)
    return min(Δh^4 / closure.νh, Δz^4 / closure.νv)
end

function cell_diffusion_timescale(closure::AnisotropicBiharmonicDiffusivity, diffusivities, grid)
    Δh = min_Δxy(grid)
    Δz = min_Δz(grid)
    max_κh = maximum(closure.κh)
    max_κv = maximum(closure.κv)
    return min(Δh^4 / closure.νh, Δz^4 / closure.νv,
               Δh^4 / max_κh, Δz^4 / max_κv)
end

function cell_diffusion_timescale(closure::SmagorinskyLilly{FT, P, <:NamedTuple{()}},
                                  diffusivities, grid) where {FT, P}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::BlasiusSmagorinsky{ML, FT, P, <:NamedTuple{()}},
                                  diffusivities, grid) where {ML, FT, P}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::AbstractSmagorinsky, diffusivities, grid)
    Δ = min_Δxyz(grid)
    min_Pr = minimum(closure.Pr)
    max_κ = maximum(closure.κ)
    max_νκ = maximum(diffusivities.νₑ.data.parent) * max(1, 1/min_Pr)
    return min(Δ^2 / max_νκ, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::RozemaAnisotropicMinimumDissipation{FT, <:NamedTuple{()}},
                                  diffusivities, grid) where FT
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::VerstappenAnisotropicMinimumDissipation{FT, PK, PN, <:NamedTuple{()}},
                                  diffusivities, grid) where {FT, PK, PN}
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

function cell_diffusion_timescale(closure::AbstractAnisotropicMinimumDissipation, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    max_κ = max(Tuple(maximum(κₑ.data.parent) for κₑ in diffusivities.κₑ)...)
    return min(Δ^2 / max_ν, Δ^2 / max_κ)
end

function cell_diffusion_timescale(closure::AbstractLeith, diffusivities, grid)
    Δ = min_Δxyz(grid)
    max_ν = maximum(diffusivities.νₑ.data.parent)
    return Δ^2 / max_ν
end

cell_diffusion_timescale(closure::Tuple, diffusivities, grid) =
    min(Tuple(cell_diffusion_timescale(c, diffusivities, grid) for c in closure)...)


text:
"""
LTESOM - Long Term Energy System Optimization Model
Copyright (C) 2021, Peking University, College of Engineering (COE),
Department of Industrial Engineering and Management
"""

function configure_glpk(solver_settings_path::String)

    solver_settings = YAML.load(open(solver_settings_path))

    # solver parameters ############################################
    TolBnd = 1e-8
    if (haskey(solver_settings, "Feasib_Tol"))
        TolBnd = solver_settings["Feasib_Tol"]
    end

    TolDj = 1e-8
    if (haskey(solver_settings, "Optimal_Tol"))
        TolDj = solver_settings["Optimal_Tol"]
    end

    Presolve = 0
    if (haskey(solver_settings, "PreSolve"))
        Presolve = solver_settings["PreSolve"]
    end

    TmLim = -1
    if (haskey(solver_settings, "TimeLimit"))
        TmLim = solver_settings["TimeLimit"]
    end

    Dual = 0
    if (haskey(solver_settings, "Dual"))
        Dual = solver_settings["Dual"]
    end
    ################################################################

    # OPTIMIZER = optimizer_with_attributes(
    #     GLPK.Optimizer,
    #     "TolBnd" => TolBnd,
    #     "TolDj" => TolDj,
    #     "Presolve" => Presolve,
    #     "TmLim" => TmLim,
    #     "Dual" => Dual
    # )
    OPTIMIZER = GLPK.Optimizer
    return OPTIMIZER
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module MCTS  #ExprSearch.MCTS

export MCTSESParams, MCTSESResult, mcts_search, exprsearch, SearchParams, SearchResult
export get_derivtree

include("DerivTreeMDPs.jl")

using Reexport
using ExprSearch
using RLESUtils, GitUtils, CPUTimeUtils, Observers, LogSystems
import RLESTypes.SymbolTable
@reexport using LinearDerivTrees  #pretty_string
using .DerivTreeMDPs
using GrammaticalEvolution
using GBMCTS

import .DerivTreeMDPs.get_fitness
import ..ExprSearch: SearchParams, SearchResult, exprsearch, ExprProblem, get_grammar, get_fitness,
    get_expr, get_derivtree

include("logdefs.jl")

type MCTSESParams <: SearchParams
  #tree params
  maxsteps::Int64

  #mdp params
  max_neg_reward::Float64
  step_reward::Float64

  #mcts
  n_iters::Int64
  searchdepth::Int64
  exploration_const::Float64
  maxmod::Bool
  q0::Float64
  seed::Int64

  logsys::LogSystem
  userargs::SymbolTable
end
MCTSESParams(maxsteps::Int64, max_neg_reward::Float64, step_reward::Float64, n_iters::Int64,
    searchdepth::Int64, exploration_const::Float64, maxmod::Bool, q0::Float64, seed::Int64,
    logsys::LogSystem=logsystem(); userargs::SymbolTable=SymbolTable()) =
        MCTSESParams(maxsteps, max_neg_reward, step_reward, n_iters, searchdepth,
        exploration_const, maxmod, q0, seed, logsys, userargs)

type MCTSESResult <: SearchResult
  tree::LinearDerivTree
  actions::Vector{Int64}
  fitness::Float64
  expr::Union{Symbol,Expr}
  best_at_eval::Int64
  totalevals::Int64
end

exprsearch(p::MCTSESParams, problem::ExprProblem) = mcts_search(p, problem::ExprProblem)

get_derivtree(result::MCTSESResult) = get_derivtree(result.tree)
get_expr(result::MCTSESResult) = result.expr
get_fitness(result::MCTSESResult) = result.fitness

function mcts_search(p::MCTSESParams, problem::ExprProblem)
  @notify_observer(p.logsys.observer, "verbose1", ["Starting MCTS search"])
  @notify_observer(p.logsys.observer, "computeinfo", ["starttime", string(now())])

  grammar = get_grammar(problem)
  tree_params = LDTParams(grammar, p.maxsteps)
  mdp_params = DerivTreeMDPParams(grammar, p.max_neg_reward, p.step_reward)

  tree = LinearDerivTree(tree_params)
  mdp = DerivTreeMDP(mdp_params, tree, problem, p.userargs)

  solver = MCTSSolver(n_iterations=p.n_iters, depth=p.searchdepth,
    exploration_constant=p.exploration_const, maxmod=p.maxmod, rng=MersenneTwister(p.seed))
  policy = MCTSPolicy(solver, mdp, q0=p.q0)

  initialize!(tree)
  s = create_state(mdp)

  tstart = CPUtime_start()
  i = 1
  while !GBMCTS.isexplored(policy.mcts.tree, s) && i <= p.n_iters
    @notify_observer(p.logsys.observer, "iteration", [i])

    simulate(policy, s, p.searchdepth) #FIXME: remove searchdepth??

    @notify_observer(p.logsys.observer, "elapsed_cpu_s", [i, CPUtime_elapsed_s(tstart)])
    best_state = policy.best_state
    @notify_observer(p.logsys.observer, "current_best", [i, -policy.best_reward,
        string(best_state.past_actions), string(get_expr(best_state))]) #report fitness instead of reward
    @notify_observer(p.logsys.observer, "mcts_tree", [i, policy.mcts.tree, s])

    i += 1
  end
  best_fitness = -policy.best_reward
  expr = get_expr(policy.best_state)
  best_actions = policy.best_state.past_actions
  @notify_observer(p.logsys.observer, "result", [best_fitness, string(expr), policy.best_at_eval, policy.totalevals])

  #meta info
  @notify_observer(p.logsys.observer, "computeinfo", ["endtime",  string(now())])
  @notify_observer(p.logsys.observer, "computeinfo", ["hostname", gethostname()])
  @notify_observer(p.logsys.observer, "computeinfo", ["gitSHA",  get_SHA(dirname(@__FILE__))])
  @notify_observer(p.logsys.observer, "computeinfo", ["cpu_time", CPUtime_elapsed_s(tstart)])
  @notify_observer(p.logsys.observer, "parameters", ["maxsteps", p.maxsteps])
  @notify_observer(p.logsys.observer, "parameters", ["max_neg_reward", p.max_neg_reward])
  @notify_observer(p.logsys.observer, "parameters", ["step_reward", p.step_reward])
  @notify_observer(p.logsys.observer, "parameters", ["discount", mdp_params.discount])
  @notify_observer(p.logsys.observer, "parameters", ["n_iters", p.n_iters])
  @notify_observer(p.logsys.observer, "parameters", ["searchdepth", p.searchdepth])
  @notify_observer(p.logsys.observer, "parameters", ["exploration_const", p.exploration_const])
  @notify_observer(p.logsys.observer, "parameters", ["maxmod", p.maxmod])
  @notify_observer(p.logsys.observer, "parameters", ["q0", p.q0])

  return MCTSESResult(tree, best_actions, best_fitness, expr, policy.best_at_eval, policy.totalevals)
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
using DataFrames, TensorToolbox, MatrixImpute

n=10^5
m=10^3
k=7
U=rand(Float64,n,k)
V=rand(Float64,k,m)
Afull=U*V
A=copy(Afull)
A = allowmissing(A)
A[rand(Float64,n,m).>0.05].=missing
t1 = time_ns()
Aopt = MatrixImpute.Impute(A,k,γ = 10000)
println("Completed Test for Matrix Completion without Side Information")

n=10^3
m=10^3
p=100
k=3
U=rand(Float64,n,k)
R=rand(Float64,k,p)
B=rand(Float64,p,m)
Afull=U*R*B
A=copy(Afull)
A = allowmissing(A)
A[rand(Float64,n,m).>0.05].=missing
t1 = time_ns()
Aopt = MatrixImpute.Impute(A,k,γ = 10000, B = B')
println("Completed Test for Matrix Completion with Side Information")

n=100
m=200
p=500
k=6
M=0.8
γ=1000000
stepsize=64
A=rand(n,m,p)
A=allowmissing(A)
A[rand(n,m,p).>1-M].=missing
t1 = time_ns()
Aopt = MatrixImpute.Impute(A,k,γ = γ, method = :fastImputeT)
println("Completed Test for Tensor Completion")


text:
export
    SWR_CH_MAX,
    SWR_FLAG_RESAMPLE,
    SwrDitherType,
    SWR_DITHER_NONE,
    SWR_DITHER_RECTANGULAR,
    SWR_DITHER_TRIANGULAR,
    SWR_DITHER_TRIANGULAR_HIGHPASS,
    SWR_DITHER_NS,
    SWR_DITHER_NS_LIPSHITZ,
    SWR_DITHER_NS_F_WEIGHTED,
    SWR_DITHER_NS_MODIFIED_E_WEIGHTED,
    SWR_DITHER_NS_IMPROVED_E_WEIGHTED,
    SWR_DITHER_NS_SHIBATA,
    SWR_DITHER_NS_LOW_SHIBATA,
    SWR_DITHER_NS_HIGH_SHIBATA,
    SWR_DITHER_NB,
    SwrEngine,
    SWR_ENGINE_SWR,
    SWR_ENGINE_SOXR,
    SWR_ENGINE_NB,
    SwrFilterType,
    SWR_FILTER_TYPE_CUBIC,
    SWR_FILTER_TYPE_BLACKMAN_NUTTALL,
    SWR_FILTER_TYPE_KAISER,
    SwrContext

const SWR_CH_MAX  =  32
const SWR_FLAG_RESAMPLE  =  1

# begin enum SwrDitherType
const SwrDitherType = UInt32
const SWR_DITHER_NONE  =  UInt32(0)
const SWR_DITHER_RECTANGULAR  =  UInt32(1)
const SWR_DITHER_TRIANGULAR  =  UInt32(2)
const SWR_DITHER_TRIANGULAR_HIGHPASS  =  UInt32(3)
const SWR_DITHER_NS  =  UInt32(64)
const SWR_DITHER_NS_LIPSHITZ  =  UInt32(65)
const SWR_DITHER_NS_F_WEIGHTED  =  UInt32(66)
const SWR_DITHER_NS_MODIFIED_E_WEIGHTED  =  UInt32(67)
const SWR_DITHER_NS_IMPROVED_E_WEIGHTED  =  UInt32(68)
const SWR_DITHER_NS_SHIBATA  =  UInt32(69)
const SWR_DITHER_NS_LOW_SHIBATA  =  UInt32(70)
const SWR_DITHER_NS_HIGH_SHIBATA  =  UInt32(71)
const SWR_DITHER_NB  =  UInt32(72)
# end enum SwrDitherType

# begin enum SwrEngine
const SwrEngine = UInt32
const SWR_ENGINE_SWR  =  UInt32(0)
const SWR_ENGINE_SOXR  =  UInt32(1)
const SWR_ENGINE_NB  =  UInt32(2)
# end enum SwrEngine

# begin enum SwrFilterType
const SwrFilterType = UInt32
const SWR_FILTER_TYPE_CUBIC  =  UInt32(0)
const SWR_FILTER_TYPE_BLACKMAN_NUTTALL  =  UInt32(1)
const SWR_FILTER_TYPE_KAISER  =  UInt32(2)
# end enum SwrFilterType

const SwrContext = Void

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
mape(y, yhat) = mean(100.0 * abs((y .- yhat) ./ y))
medape(y, yhat) = median(100.0 * abs((y .- yhat) ./ y))

function eval_sparse_reg(A, b, xtrue, xhat)
    bhat = A * xhat
    strue = Set(find(xtrue))
    shat = Set(find(xhat))
    missing = setdiff(strue, shat)
    common = setdiff(strue, missing)
    additional = setdiff(shat, strue)
    mapeval = mape(b, bhat)
    medapeval = medape(b, bhat)
    println("MAPE: $(round(mapeval,1))%, MEDAPE: $(round(medapeval,1))%\n",
        "  Support: common = $(length(common)), missing = $(length(missing)), additional = $(length(additional))")
    return [mapeval, medapeval, length(common), length(missing), length(additional)]
end


text:
using LinearAlgebra
using SparseArrays
using Random
Random.seed!(2)
include("../src/SocpSolver.jl")
using .SocpSolver

A = sparse(randn(2, 5))
c = rand(5)
b = A*randn(5)
G = -sparse(Diagonal(ones(5)))
h = zeros(5)
prob = LinearProgram(A, G, c, b, h)
sol, stats = solve(prob, Settings());


text:
# Copyright (c) 2019 Guido Kraemer
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

# struct Header
#     version           ::UInt32
#     previous_hash     ::UInt256 #Vector{UInt8} # length 32
#     merkle_root       ::UInt256 #Vector{UInt8} # length 32
#     timestamp         ::UInt32
#     difficulty_target ::UInt32
#     nonce            ::UInt32
# end
#
# The following implementation hopefully is
# more efficient, because no conversions at reading time have to be made:

"""
    Header

Data Structure representing the Header of a Block in the Bitcoin blockchain.

Data are store as an `NTuple{80, UInt8}` without parsinf per se.
The elements of the `Header` can be accessed by `header.element`.

```julia
header.version
header.prevhash
header.merkleroot
header.time
header.bits
header.nonce
```
"""
struct Header
    data :: Vector{UInt8}
end

"""
    Header(x::IO) -> Header

Parse `Header` from an `IO`
"""
Header(io::IO) = Header(read(io, 80))

@inline Base.getindex(x::Header, r) = x.data[r]

@inline function Base.getproperty(x::Header, d::Symbol)
    if     d == :version    ltoh(reinterpret(UInt32, x.data[1:4])[1])
    elseif d == :prevhash   x.data[5:36]
    elseif d == :merkleroot x.data[37:68]
    elseif d == :time       ltoh(reinterpret(UInt32, x.data[69:72])[1])
    elseif d == :bits       ltoh(reinterpret(UInt32, x.data[73:76])[1])
    elseif d == :nonce      ltoh(reinterpret(UInt32, x.data[77:80])[1])
    else getfield(x, d)
    end
end

function Base.propertynames(::Type{Header}, private = false)
    (:version, :prevhash, :merkleroot, :time, :bits, :nonce,
     fieldnames(Header)...)
end

function Base.show(io::IO, header::Header)
    println(io, "  Version:    " * string(header.version,    base = 16))
    println(io, "  Prev Hash:  " * bytes2hex(header.prevhash))
    println(io, "  Root:       " * bytes2hex(header.merkleroot))
    println(io, "  Time:       " * string(header.time,       base = 10))
    println(io, "  Difficulty: " * string(header.bits,       base = 16))
    println(io, "  Nonce:      " * string(header.nonce,      base = 10))
end

"""
    serialize(header::Header) -> Vector{UInt8}

Returns the 80 byte Vector{UInt8} for the block header
"""
function serialize(header::Header)
    result = Vector(reinterpret(UInt8, [htol(header.version)]))
    append!(result, header.prevhash)
    append!(result, header.merkleroot)
    append!(result, Vector(reinterpret(UInt8, [htol(header.time)])))
    append!(result, Vector(reinterpret(UInt8, [htol(header.bits)])))
    append!(result, Vector(reinterpret(UInt8, [htol(header.nonce)])))
end

"""
    bip9(header::Header) -> Bool

Returns whether this header is signaling readiness for BIP9

    BIP9 is signalled if the top 3 bits are 001
    remember version is 32 bytes so right shift 29 (>> 29) and see if
    that is 001
"""
bip9(header::Header) = header.version >> 29 == 0b001

"""
    bip91(header::Header) -> Bool

Returns whether this header is signaling readiness for BIP91

    BIP91 is signalled if the 5th bit from the right is 1
    shift 4 bits to the right and see if the last bit is 1
"""
bip91(header::Header) = header.version >> 4 & 1 == 1

"""
    bip141(header::Header) - > Bool

Returns whether this header is signaling readiness for BIP141

    BIP91 is signalled if the 2nd bit from the right is 1
    shift 1 bit to the right and see if the last bit is 1
"""
bip141(header::Header) = header.version >> 1 & 1 == 1

"""
    target(header::Header) -> BigInt

Returns the proof-of-work target based on the bits

    last byte is exponent
    the first three bytes are the coefficient in little endian
    the formula is: coefficient * 256**(exponent-3)
"""
function target(header::Header)
    exponent = header.bits >> 24
    coefficient = header.bits & 0x00ffffff
    return coefficient * big(256)^(exponent - 3)
end

"""
    difficulty(header::Header) -> BigInt

Returns the header difficulty based on the bits

    difficulty is (target of lowest difficulty) / (header's target)
    lowest difficulty has bits that equal 0xffff001d
"""
function difficulty(header::Header)
    lowest = 0xffff * big(256)^(0x1d - 3)
    return div(lowest, target(header))
end

"""
    check_pow(header::Header) -> Bool

Returns whether this header satisfies proof of work

    get the hash256 of the serialization of this header
    interpret this hash as a little-endian number
    return whether this integer is less than the target
"""
function check_pow(header::Header)
    block_hash = hash256(header)
    proof = to_int(block_hash, little_endian=true)
    return proof < target(header)
end


text:
@testset "Test quadrature rules" begin

using HDiscontinuousGalerkin

# Test GrundmannMoeller Quadrature
quad_rule = QuadratureRule{2,RefTetrahedron}(GrundmannMoeller(),0)
@test getweights(quad_rule) ≈ [0.5]
@test getpoints(quad_rule)[1] ≈ [1//3,1//3]
quad_rule = QuadratureRule{2,RefTetrahedron}(GrundmannMoeller(),1)
@test getweights(quad_rule) ≈ 0.5*[0.520833333333333 ,0.520833333333333,0.520833333333333 ,-0.562500000000000]
@test getpoints(quad_rule)[1] ≈ [1//5,1//5]
@test getpoints(quad_rule)[2] ≈ [3//5,1//5]
@test getpoints(quad_rule)[3] ≈ [1//5,3//5]
@test getpoints(quad_rule)[4] ≈ [1//3,1//3]

for i in 0:5
    quad_rule = QuadratureRule{2,RefTetrahedron}(GrundmannMoeller(),i)
    @test sum(getweights(quad_rule)) ≈ 0.5
end

# Test Strang Quadrature
for i in 1:5
    quad_rule = QuadratureRule{2,RefTetrahedron}(Strang(),i)
    @test sum(getweights(quad_rule)) ≈ 0.5
end

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
using Test
using JuMP
using Cbc
using JSON
using REoptLite

if "Xpress" in ARGS
    @testset "test_with_xpress" begin
        include("test_with_xpress.jl")
    end

elseif "CPLEX" in ARGS
    @testset "test_with_cplex" begin
        include("test_with_cplex.jl")
    end

else  # run Cbc tests
    #=
    Cbc can be very slow with integers: if ElectricUtility.allow_simultaneous_export_import == false for example (a
        test with Xpress and CPLEX) the problem does not solve in under five hours.
    So, we test some simple cases with Cbc to stay aware of solve times with a free solver.
    =#
    @testset "January Export Rates" begin
        model = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
        data = JSON.parsefile("./scenarios/monthly_rate.json")

        # create wholesale_rate with compensation in January > retail rate
        jan_rate = data["ElectricTariff"]["monthly_energy_rates"][1]
        data["ElectricTariff"]["wholesale_rate"] =
            append!(repeat([jan_rate + 0.1], 31 * 24), repeat([0.0], 8760 - 31*24))
        data["ElectricTariff"]["monthly_demand_rates"] = repeat([0], 12)

        s = Scenario(data)
        inputs = REoptInputs(s)
        results = run_reopt(model, inputs)

        @test results["PV"]["size_kw"] ≈ 70.3084 atol=0.01
        @test results["Financial"]["lcc_us_dollars"] ≈ 430747.0 rtol=1e-5 # with levelization_factor hack the LCC is within 5e-5 of REopt Lite API LCC
        @test all(x == 0.0 for x in results["PV"]["year_one_to_load_series_kw"][1:744])
    end

    @testset "Blended tariff" begin
        model = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
        results = run_reopt(model, "./scenarios/no_techs.json")
        @test results["ElectricTariff"]["year_one_energy_cost_us_dollars"] ≈ 1000.0
        @test results["ElectricTariff"]["year_one_demand_cost_us_dollars"] ≈ 136.99
    end

    @testset "Solar and Storage" begin
        model = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
        results = run_reopt(model, "./scenarios/pv_storage.json")

        @test results["PV"]["size_kw"] ≈ 216.6667 atol=0.01
        @test results["Financial"]["lcc_us_dollars"] ≈ 1.23887e7 rtol=1e-5
        @test results["Storage"]["size_kw"] ≈ 55.9 atol=0.1
        @test results["Storage"]["size_kwh"] ≈ 78.9 atol=0.1
    end

    @testset "Outage with Generator" begin
        model = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
        results = run_reopt(model, "./scenarios/generator.json")
        @test results["Generator"]["size_kw"] ≈ 8.12 atol=0.01
        @test (sum(results["Generator"]["year_one_to_load_series_kw"][i] for i in 1:9) +
            sum(results["Generator"]["year_one_to_load_series_kw"][i] for i in 13:8760)) == 0
        p = REoptInputs("./scenarios/generator.json")
        simresults = simulate_outages(results, p)
        @test simresults["resilience_hours_max"] == 11
    end

    # TODO test MPC with outages
    @testset "MPC" begin
        model = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
        r = run_mpc(model, "./scenarios/mpc.json")
        @test maximum(r["ElectricUtility"]["to_load_series_kw"][1:15]) <= 98.0
        @test maximum(r["ElectricUtility"]["to_load_series_kw"][16:24]) <= 97.0
        @test sum(r["PV"]["to_grid_series_kw"]) ≈ 0
    end

    @testset "Complex Incentives" begin
        """
        This test was compared against the API test:
            reo.tests.test_reopt_url.EntryResourceTest.test_complex_incentives
        when using the hardcoded levelization_factor in this package's REoptInputs function.
        The two LCC's matched within 0.00005%. (The Julia pkg LCC is 1.0971991e7)
        """
        model = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
        results = run_reopt(model, "./scenarios/incentives.json")
        @test results["Financial"]["lcc_us_dollars"] ≈ 1.1152536e7 atol=5e4
        # The Cbc LCC is 1.7% higher than the Xpress LCC ? Probably due to integer issues in Cbc
    end

    @testset "Wind" begin
        model = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
        results = run_reopt(model, "./scenarios/wind.json")
        @test results["Wind"]["size_kw"] ≈ 3752 atol=0.1
        @test results["Financial"]["lcc_us_dollars"] ≈ 8.591017e6 rtol=1e-5
        #=
        0.5% higher LCC in this package as compared to API ? 8,591,017 vs 8,551,172
        - both have zero curtailment
        - same energy to grid: 5,839,317 vs 5,839,322
        - same energy to load: 4,160,683 vs 4,160,677
        - same city: Boulder
        - same total wind prod factor

        REoptLite.jl has:
        - bigger turbine: 3752 vs 3735
        - net_capital_costs_plus_om_us_dollars: 8,576,590 vs. 8,537,480

        TODO: will these discrepancies be addressed once NMIL binaries are added?
        =#
    end

    ## much too slow with Cbc (killed after 8 hours)
    # @testset "Minimize Unserved Load" begin
    #     m = Model(optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0))
    #     results = run_reopt(m, "./scenarios/outage.json")

    #     @test results["expected_outage_cost"] ≈ 0
    #     @test results["total_unserved_load"] ≈ 0
    #     @test value(m[:binMGTechUsed]["Generator"]) == 1
    #     @test value(m[:binMGTechUsed]["PV"]) == 0
    #     @test value(m[:binMGStorageUsed]) == 1
    #     @test results["lcc"] ≈ 1.5291695e7
    # end
end

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite


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

JLLWrappers.@generate_main_file_header("CGAL")
JLLWrappers.@generate_main_file("CGAL", UUID("8fcd9439-76b0-55f4-a525-bad0597c05d8"))
end  # module CGAL_jll


text:
push!(LOAD_PATH,"../src/")

using Documenter, Main.LifeVal

makedocs(sitename="LifeVal Documentation")


text:
# Regression

## Auxiliary

function lreg_chkdims(X::AbstractMatrix, Y::AbstractVecOrMat, trans::Bool)
    mX, nX = size(X)
    dX = ifelse(trans, mX, nX)
    dY = ifelse(trans, nX, mX)
    size(Y, 1) == dY || throw(DimensionMismatch("Dimensions of X and Y mismatch."))
    return dX
end

lrsoltype(::AbstractVector{T}) where T = Vector{T}
lrsoltype(::AbstractMatrix{T}) where T = Matrix{T}

_vaug(X::AbstractMatrix{T}) where T = vcat(X, ones(T, 1, size(X,2)))::Matrix{T}
_haug(X::AbstractMatrix{T}) where T = hcat(X, ones(T, size(X,1), 1))::Matrix{T}

## Linear Least Square Regression

"""
    llsq(X, y; ...)

Solve the linear least square problem.

Here, `y` can be either a vector, or a matrix where each column is a response vector.

This function accepts two keyword arguments:

- `dims`: whether input observations are stored as rows (`1`) or columns (`2`). (default is `1`)
- `bias`: whether to include the bias term `b`. (default is `true`)

The function results the solution `a`. In particular, when `y` is a vector (matrix), `a` is also a vector (matrix). If `bias` is true, then the returned array is augmented as `[a; b]`.
"""
function llsq(X::AbstractMatrix{T}, Y::AbstractVecOrMat{T};
              trans::Bool=false, bias::Bool=true,
              dims::Union{Integer,Nothing}=nothing) where {T<:Real}
    if dims === nothing && trans
        Base.depwarn("`trans` argument is deprecated, use llsq(X, Y, dims=d) instead.", :trans)
        dims = 1
    end
    if dims == 2
        mX, nX = size(X)
        size(Y, 1) == nX || throw(DimensionMismatch("Dimensions of X and Y mismatch."))
        mX <= nX || error("mX <= nX is required when trans is false.")
    else
        mX, nX = size(X)
        size(Y, 1) == mX || throw(DimensionMismatch("Dimensions of X and Y mismatch."))
        mX >= nX || error("mX >= nX is required when trans is false.")
    end
    _ridge(X, Y, zero(T), dims == 2, bias)
end
llsq(x::AbstractVector{T}, y::AbstractVector{T}; kwargs...) where {T<:Real} =
    llsq(x[:,:], y; dims=1, kwargs...)

## Ridge Regression (Tikhonov regularization)

"""

    ridge(X, y, r; ...)

Solve the ridge regression problem.

Here, ``y`` can be either a vector, or a matrix where each column is a response vector.

The argument `r` gives the quadratic regularization matrix ``Q``, which can be in either of the following forms:

- `r` is a real scalar, then ``Q`` is considered to be `r * eye(n)`, where `n` is the dimension of `a`.
- `r` is a real vector, then ``Q`` is considered to be `diagm(r)`.
- `r` is a real symmetric matrix, then ``Q`` is simply considered to be `r`.

This function accepts two keyword arguments:

- `dims`: whether input observations are stored as rows (`1`) or columns (`2`). (default is `1`)
- `bias`: whether to include the bias term `b`. (default is `true`)

The function results the solution `a`. In particular, when `y` is a vector (matrix), `a` is also a vector (matrix). If `bias` is true, then the returned array is augmented as `[a; b]`.
"""
function ridge(X::AbstractMatrix{T}, Y::AbstractVecOrMat{T}, r::Union{Real, AbstractVecOrMat};
               trans::Bool=false, bias::Bool=true,
               dims::Union{Integer,Nothing}=nothing) where {T<:Real}
    if dims === nothing && trans
        Base.depwarn("`trans` argument is deprecated, use ridge(X, Y, r, dims=d) instead.", :trans)
        dims = 1
    end
    d = lreg_chkdims(X, Y, dims == 2)
    if isa(r, Real)
        r >= zero(r) || error("r must be non-negative.")
        r = convert(T <: AbstractFloat ? T : Float64, r)
    elseif isa(r, AbstractVector)
        length(r) == d || throw(DimensionMismatch("Incorrect length of r."))
    elseif isa(r, AbstractMatrix)
        size(r) == (d, d) || throw(DimensionMismatch("Incorrect size of r."))
    end
    _ridge(X, Y, r, dims == 2, bias)
end
ridge(x::AbstractVector{T}, y::AbstractVector{T}, r::Union{Real, AbstractVecOrMat};
      kwargs...) where {T<:Real} = ridge(x[:,:], y, r; dims=1, kwargs...)

### implementation

function _ridge(_X::AbstractMatrix{T}, _Y::AbstractVecOrMat{T},
                r::Union{Real, AbstractVecOrMat}, trans::Bool, bias::Bool) where {T<:Real}
    # convert integer data to Float64
    X = T <: AbstractFloat ? _X : convert(Array{Float64}, _X)
    Y = T <: AbstractFloat ? _Y : convert(Array{Float64}, _Y)
    if bias
        if trans
            X_ = _vaug(X)
            A = cholesky!(Hermitian(_ridge_reg!(X_ * transpose(X_), r, bias))) \ (X_ * Y)
        else
            X_ = _haug(X)
            A = cholesky!(Hermitian(_ridge_reg!(X_'X_, r, bias))) \ (X_'Y)
        end
    else
        if trans
            A = cholesky!(Hermitian(_ridge_reg!(X * X', r, bias))) \ (X * Y)
        else
            A = cholesky!(Hermitian(_ridge_reg!(X'X, r, bias))) \ (X'Y)
        end
    end
    return A::lrsoltype(Y)
end

function _ridge_reg!(Q::Matrix, r::Real, bias::Bool)
    if r > zero(r)
        n = size(Q, 1) - Int(bias)
        for i = 1:n
            @inbounds Q[i,i] += r
        end
    end
    return Q
end

function _ridge_reg!(Q::AbstractMatrix, r::AbstractVector, bias::Bool)
    n = size(Q, 1) - Int(bias)
    @assert length(r) == n
    for i = 1:n
        @inbounds Q[i,i] += r[i]
    end
    return Q
end

function _ridge_reg!(Q::AbstractMatrix, r::AbstractMatrix, bias::Bool)
    n = size(Q, 1) - Int(bias)
    @assert size(r) == (n, n)
    for j = 1:n, i = 1:n
        @inbounds Q[i,j] += r[i,j]
    end
    return Q
end

## Isotonic Regression

"""
    isotonic(x, y[, w])

Solve the isotonic regression problem using the pool adjacent violators algorithm[^1].

Here `x` is the regressor vector, `y` is response vector, and `w` is an optional
weights vector.

The function returns a prediction vector of the same size as the regressor vector `x`.
"""
function isotonic(x::AbstractVector{T}, y::AbstractVector{T},
                  w::AbstractVector{T} = ones(T, length(y))) where {T<:Real}
    n = length(x)
    n == length(y) || throw(DimensionMismatch("Dimensions of x and y mismatch."))

    idx = sortperm(x)

    # PVA algorithm
    J = map(i->(Float64(y[i]*w[i]), w[i], [i]), idx)
    i = 1
    B₀ = J[i]
    while i < length(J)
        B₊ = J[i+1]
        if B₀[1] <= B₊[1] # step 1
            B₀ = B₊
            i += 1
        else # step 2
            ww = B₀[2] + B₊[2]
            J[i] = ((B₀[1]*B₀[2]+B₊[1]*B₊[2])/ww, ww, append!(B₀[3], B₊[3]))
            deleteat!(J, i+1)
            B₀ = J[i]
            while i > 1 # step 2.1
                B₋ = J[i-1]
                if B₀[1] <= B₋[1]
                    ww = B₀[2] + B₋[2]
                    J[i] = ((B₀[1]*B₀[2]+B₋[1]*B₋[2])/ww, ww, append!(B₀[3], B₋[3]))
                    deleteat!(J, i-1)
                    i-=1
                else
                    break
                end
            end
        end
    end
    [y for (y,w,ii) in J for i in sort(ii)]
end


text:
using CompScienceMeshes, BEAST
Base.getindex(m::CompScienceMeshes.AbstractMesh, i::Int) = cells(m)[i]
cellvertices(m::CompScienceMeshes.AbstractMesh, i::Int) = vertices(m)[cells(m)[i]]

m = meshsphere(1.0, 0.35)
# m = meshcuboid(1.0, 1.0, 0.2, 0.2)

X = raviartthomas(m)
Y = buffachristiansen(m, ibscaled=true)

Id = BEAST.Identity()
Nx = BEAST.NCross()

E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=1.0)
G = -assemble(Id, Y, Y)

b = assemble((n×E)×n, Y)
x = G \ b

# x = zeros(Float64,numfunctions(Y))
# x[20] = 1
fcr, geo = facecurrents(x, Y)

import PlotlyJS
using LinearAlgebra
PlotlyJS.plot(patch(geo.mesh, norm.(fcr)))
# PlotlyJS.plot(patch(geo.mesh, rand(numcells(geo))))

edges = skeleton(m,1)
verts = skeleton(m,0)

vertex_list = [v[1] for v in cells(verts)]
Q = lagrangec0d1(m, vertex_list, Val{3})
Z = curl(Q)

Gzy = assemble(Id, Z, Y)

# vtoc, nc = vertextocellmap(edges)
D = connectivity(verts, edges)
@assert size(D) == (numcells(edges), numcells(verts))

los = geometry(Y)
# Choose a vertex:
for n_vert in 1:numcells(verts)
    @assert cellvertices(verts, n_vert)[1] ≈ Q.pos[n_vert]
# find the connected edges
    # n_edges = vtoc[n_vert,1:nc[n_vert]]
    using SparseArrays
    rows = rowvals(D)
    vals = nonzeros(D)
    for k in nzrange(D,n_vert)
        n_edge = rows[k]
        cht = chart(edges, edges[n_edge])
        ctr = cartesian(center(cht))
        vt1 = cht.vertices[1]
        lvt = los.mesh.vertices[numvertices(m)+n_edge]
        pos = Y.pos[n_edge]
        @assert norm(lvt - pos) ≤ 1e-8
        @assert norm(cross(pos-vt1, ctr-vt1)) ≤ 1e-8
    end
end

n_vert = 13
x = [float(i == n_vert) for i in 1:numfunctions(Q)]
fcr1, geo1 = facecurrents(x, Q)
colors = getindex.(fcr1,1)
PlotlyJS.plot(patch(m, colors))

divY = divergence(Y)
fcr2, geo2 = facecurrents(rand(numfunctions(divY)), divY)
colors = getindex.(fcr2,1)
PlotlyJS.plot(patch(geo2.mesh, colors))

for n_fn in 1:numfunctions(divY)
    n_fn = 21
    charges = zeros(numcells(los))
    for shape in divY.fns[n_fn]
        f = shape.cellid
        ch = chart(los, los[f])
        charges[f] += shape.coeff #* volume(ch)
    end
    nzs = (charges[(charges .≈ 0) .== false])
    xtr = extrema(nzs)
    @assert all( Vector{Bool}(nzs .≈ xtr[1]) .| Vector{Bool}(nzs .≈ xtr[2]))
end


text:
#=
    Now, we define the new data types.
    Model => ModelProto
    Graph => GraphProto
    Node  => NodeProto
    Attribute => AttributeProto
    The new types will consist of Julian attributes.

    The purpose of dealing with these newer type is to make
    the process simpler and easier to debug.
=#

module Types

mutable struct ValueInfo
    name::AbstractString
    doc_string::AbstractString
end

mutable struct Node
    input::Vector{AbstractString}
    output::Vector{AbstractString}
    name::AbstractString
    op_type::AbstractString             # Done!
    domain::AbstractString
    attribute::Dict{Any, Any}         # AttributeProto to Dict
    doc_string::AbstractString
end

mutable struct Graph
    node::Array{Any, 1}
    name::AbstractString
    initializer::Dict{Any ,Any}             #Storing the array data instead of the tensorproto vector.
    doc_string::AbstractString              #in Dict format.
    input::Array{ValueInfo ,1}              # ValueInfoProto -> ValueInfo
    output::Array{ValueInfo, 1}                    #
    value_info::Array{ValueInfo, 1}                # Done!
end

mutable struct Model
    ir_version::Int64
    opset_import::Array{Any, 1}              #OperatorSetIdProto to Dict
    producer_name::AbstractString
    producer_version::AbstractString            # Done!
    domain::AbstractString
    model_version::Int64
    doc_string::AbstractString
    graph::Graph
    metadata_props::Array{Any, 1}            #StringStringEntryProto to Dict
end

export Model, Graph, Node, ValueInfo
end


text:
# using LCMCore
# using CaesarLCMTypes
using Rotations, CoordinateTransformations
using MeshCat
using TransformUtils
using AprilTags

include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","racecar","cameraUtils.jl"))

vis = Visualizer()
open(vis)

setobject!(vis["home"], Triad(1.0))

detector = AprilTagDetector()
# setobject!(vis["tag01"], Triad(0.2))
# settransform!(vis, Translation(-0.5, -0.5, 0))

fx = fy = 80.0
cx,cy = (320,240)
k1,k2 = (0.0,0.0)
camK = [fx 0 cx; 0 fy cy; 0 0 1]

tagsize = 0.172

function tagPicture(vis::Visualizer,msg::apriltag_t, drawdict)
    # m=homography_to_pose(Array{Float64,2}(msg.homography),300,300,320,240)
    m = msg.pose
    tkey = "tag$(msg.id)"
    if !haskey(drawdict, tkey)
        setobject!(vis[tkey], Triad(0.2))
        drawdict[tkey] = 1
    end
    # wTc=LinearMap(RotX(pi/2))
    # cTt=Translation(m[1:3,4]...) ∘ LinearMap(Quat(m[1:3,1:3]))
    # settransform!(vis[tkey], wTc ∘ cTt)
    P = Vector{Vector{Float64}}()
    for i in 1:4
      push!(P, msg.P[i,:])
    end
    tag = AprilTag(msg.familyName, msg.id, msg.hammingDistance, Float32(0.0), P)
    T,Q, = getAprilTagTransform(tag,
                                camK,
                                k1,
                                k2,
                                tagsize )
    cTt = T ∘ LinearMap(Q)
    settransform!(vis[tkey], cTt)
    nothing
end

function typed_callback(channel::String, msg::apriltag_t, drawdict)
    tagPicture(vis,msg, drawdict)
end

drawdict = Dict{String, Int}()

lcm=LCM()
subscribe(lcm, "TEST_CHANNEL", (m,c)->typed_callback(m,c,drawdict), apriltag_t)

while true
  handle(lcm)
end


text:
module POMDPGym
    using POMDPs
    using POMDPModels
    import POMDPModelTools:render
    using POMDPModelTools
    using POMDPPolicies
    using PyCall
    import PyCall: hasproperty
    using Parameters
    using Random
    using ImageCore
    using ImageShow
    using Images
    using Distributions
    import Cairo
    using Compose
    using StaticArrays
    using LinearAlgebra
    using GridInterpolations

    include("pycalls.jl")

    export GymPOMDP, reset!, step!, finished, render, close, init_mujoco_render
    include("gym_pomdp.jl")

    export torgb, preproc_atari_frame, AtariPOMDP
    include("atari_helpers.jl")

    export GridWorldMDP, GWPos, LavaWorldMDP, PendulumMDP, PendulumPOMDP, InvertedPendulumPOMDP, InvertedPendulumMDP, ContinuousBanditMDP, random_lava, ContinuumWorldMDP, Circle, Vec2, Vec4, LagrangeConstrainedPOMDP, RewardModPOMDP, isfailure,
    CollisionAvoidanceMDP, OptimalCollisionAvoidancePolicy, EpisodicSafetyGym,
    CartPoleMDP
    include("extra/gridworld.jl")
    include("extra/lavaworld.jl")
    include("extra/pendulum.jl")
    include("extra/continuous_bandit.jl")
    include("extra/continuumworld.jl")
    include("extra/lagrange_constraints.jl")
    include("extra/reward_mod.jl")
    include("extra/collision_avoidance.jl")
    include("extra/episodic_safety_gym.jl")
    include("extra/cartpole.jl")

    export AdversarialMDP, AdversarialPOMDP, AdditiveAdversarialMDP, AdditiveAdversarialPOMDP, disturbances, disturbanceindex
    include("extra/adversarialmdp.jl")
end # module

max_stars_repo_name: ancorso/POMDPGym

max_issues_repo_name: ancorso/POMDPGym

max_forks_repo_name: ancorso/POMDPGym


text:
export Signal
export signals

"""
    $(TYPEDEF)

A signal type for indicating a particular signal should be processed.
"""
struct Signal{S} end
@inline Signal(S::Symbol) = Signal{S}()

"""
Issue a signal to all problem owners.

# Example

To signal that the correction step has finished (irrespective of success)
```
signal!(Signal(:post_correct), problem)
```
"""
function signal! end

"""
    $(SIGNATURES)

Return a list of documented signals. This may not be an exhaustive list of all possible
signals.
"""
function signals()
    signal = Symbol[]
    for m in methods(signals)
        if (m.nargs == 2) && (m.sig.parameters[2] <: Signal)
            par = m.sig.parameters[2]
            if par isa Union
                while par.b isa Union
                    push!(signal, par.a.parameters[1])
                    par = par.b
                end
                push!(signal, par.a.parameters[1])
                push!(signal, par.b.parameters[1])
            else
                push!(signal, par.parameters[1])
            end
        end
    end
    return signal
end

"""
    $(SIGNATURES)

Return the docstring for a particular signal.
"""
signals(signal::Symbol) = signals(Signal(signal))

"""
    Signal(:initial_state)

Indicates that the function [`get_initial_state`](@ref) has been called. Problems should use
this as an opportunity to update initial conditions and initial data.

# Signature

```
owner(signal, [indices], problem, u, data)
```

# Arguments

- `signal::Signal{:initial_state}`
- `indices::NTuple{N,Int64} where N`: (optional) indices of any [`Var`](@ref),
    [`Data`](@ref), [`Func`](@ref), or [`Problem`](@ref) requested using
    [`pass_indices`](@ref)
- `problem::ClosedProblem`: the underlying problem structure
- `u::Vector{Vector{T}} where {T<:Number}`: vector of initial values
- `t::Vector{Vector{T}} where {T<:Number}`: vector of initial tangents
- `data::Vector{Any}`: vector of initial data

# Notes

- While `u`, `t`, and `data` can be mutated at will, the values contained within them should
    not be since they are user supplied. Instead, they should be copied before mutating.

# Example

```
# Change the first element of the first initial value to be 1.0.

# BAD - mutating a user supplied value directly
u[1][1] = 1.0

# GOOD - copying a user supplied value before mutating
u[1] = copy(u[1])
u[1][1] = 1.0
```
"""
signals(signal::Signal{:initial_state}) = (@doc signals(::typeof(signal)))  # brackets are required due to special casing of @doc

"""
    Signal(:pre_init)
    Signal(:post_init)

These signals occur when [`init!`](@ref) is called with a [`ClosedProblem`](@ref).

# Signature

```
owner(signal, [indices], problem)
```

# Arguments

- `signal::Signal{:initial_state}`
- `indices::NTuple{N,Int64} where N`: (optional) indices of any [`Var`](@ref),
    [`Data`](@ref), [`Func`](@ref), or [`Problem`](@ref) requested using
    [`pass_indices`](@ref)
- `problem::ClosedProblem`: the underlying problem structure
"""
function signals(signal::Union{Signal{:pre_init},Signal{:post_init}})
    return (@doc signals(::typeof(signal)))
end  # brackets are required due to special casing of @doc


text:
using Documenter
using BeaData

makedocs(
    modules = [BeaData],
    format = Documenter.Formats.HTML,
    sitename = "BeaData.jl",
    pages = Any[
        "Home" => "index.md",
        "Index" => "lib/public.md"
    ],
    doctest = false
)

deploydocs(
    repo = "github.com/stephenbnicar/BeaData.jl.git",
    target = "build",
    julia  = "0.6",
    deps = nothing,
    make = nothing,
)


text:
using Mimi
using Distributions
using DataFrames
using CSV
using Random

# standard output folder:
global dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/FixSCC"
global scenario = "NDCs"

function set_default_bools()
    # use annual version
    global use_annual = true
    # use interannual variability
    global use_interannvar = true

    # set default monte carlo parameters to the full set of probability distributions
    global use_only_varMC = false
    global use_no_varMC = false
    # set default interpolation parameters
    global use_linear = false
    global use_logburke = true
    global use_logpopulation = false

    # set only when doing specific variability sensitivity analysis
    global globallyset_use_interannvar = false
end

# set default; interannual global mean standard deviation
global set_gvarsd = 0.11294 # default from https://github.com/jkikstra/climvar
# in case you so desire, set one default; interannual regional mean standard deviation [normally read from rvarsd_regionalvariabilitystandarddeviation.csv]
# global v_multiplier = 1 # introduced for sensitivity analysis for variability; default = 1 -> variability according to empirical results
set_default_bools()

# # Run.
# # with variability
# numMCruns = 20000
# include("getpagefunction.jl")
# include("utils/mctools.jl")
# include("mcs_scc_dist.jl")
#
# get_pre_scc_mcs(numMCruns, 2020, dir_output, scenario)

# Run.
# without variability
global use_interannvar = false
dir_output = "C:/Users/jarmo/Documents/GitHub/PAGEoutput/FixSCC/Novar"

numMCruns = 30000
include("getpagefunction.jl")
include("utils/mctools.jl")
include("mcs_scc_dist.jl")

get_pre_scc_mcs(numMCruns, 2020, dir_output, scenario)

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
module FunctionalResponse
  using Test
  using BioEnergeticFoodWebs

  foodchain = [0 1 0; 0 0 1; 0 0 0]
  omnivory = [0 1 1; 0 0 1; 0 0 0]
  init_biomass = [0.5, 0.3, 0.8]
  for a in [foodchain, omnivory]
      for i in [1.0, 2.0]
          for j in [0.0, 1.0]
              p = model_parameters(a, h = i, c = j)
              ω = a ./ sum(a, dims = 2)
              ω[isnan.(ω)] .= 0.0
              @test p[:w] == ω
              num = ω .* (init_biomass' .^ i)
              bm = zeros(eltype(init_biomass), (length(init_biomass), length(init_biomass)))
              BioEnergeticFoodWebs.fill_bm_matrix!(bm, init_biomass, p[:w], p[:A], p[:h])
              @test num == bm
              hsd = [0.5, 0.5, 0.0] .^ i
              @test p[:Γh] == hsd
              ip = j .* hsd .* init_biomass
              fa = sum(num, dims = 2)
              den = hsd .+ ip .+ fa
              f_calc = num./den
              f_calc[isnan.(f_calc)] .= 0.0
              f = zeros(eltype(init_biomass), (length(init_biomass), length(init_biomass)))
              BioEnergeticFoodWebs.fill_F_matrix!(f, num, init_biomass, hsd, j)
              @test f == f_calc

              p = model_parameters(a, h = i, c = j, vertebrates = [true, false, false])
              x = [0.88, 0.3141, 0.138] #metabolic rates
              y = [4.0, 8.0, 0.0] #max. consumption efficiency
              eff = [0.0 0.85 0.45; 0.0 0.0 0.45; 0.0 0.0 0.0] #max. efficiency

              inflows = x .* y .* init_biomass .* f_calc
              outflows = inflows ./ eff
              outflows[isnan.(outflows)] .= 0.0

              infl, outfl = BioEnergeticFoodWebs.consumption(p, init_biomass)
              @test infl == vec(sum(inflows, dims = 2))
              @test outfl == vec(sum(outflows, dims = 1))
          end
      end
  end
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
using MultivariateFunctions
using Random
using DataFrames
using Distributions
using DataStructures

Random.seed!(1992)
nObs = 1000
dd = DataFrame()
dd[:x] = rand( Normal(),nObs) + 0.1 .* rand( Normal(),nObs)
dd[:z] = rand( Normal(),nObs) + 0.1 .* rand( Normal(),nObs)
dd[:w] = (0.5 .* rand( Normal(),nObs)) .+ 0.7.*(dd[:z] .- dd[:x]) + 0.1 .* rand( Normal(),nObs)
dd[:y] = (dd[:x] .*dd[:w] ) .* (dd[:z] .- dd[:w]) .+ dd[:x] + rand( Normal(),nObs)
dd[7,:y] = 1.0
y = :y
x_variables = Set{Symbol}([:w, :x, :z])

dd2 = deepcopy(dd)
dd2[dd2[:x] .< 0.0, :y] = 2.7
dd2[(dd2[:x] .> 0.0) .& (dd2[:w] .> 0.0), :y] = 1.3
dd2[(dd2[:x] .> 0.0) .& (dd2[:w] .< 0.0), :y] = 1.1

# Recursive Partitioning
rp_1, rp_reg_1 = create_recursive_partitioning(dd, y, x_variables, 3; rel_tol = 1e-3)
rp_2, rp_reg_2 = create_recursive_partitioning(dd, y, x_variables, 3; rel_tol = 1e-10)
SSR_1 = sum((rp_reg_1.rr.mu .- rp_reg_1.rr.y) .^ 2)
rp_3, rp_reg_3 = create_recursive_partitioning(dd, y, x_variables, 4; rel_tol = 1e-3)
SSR_3 = sum((rp_reg_3.rr.mu .- rp_reg_3.rr.y) .^ 2)
SSR_3 <= SSR_1
rp_4, rp_reg_4 = create_recursive_partitioning(dd, y, x_variables, 7; rel_tol = 1e-3)
SSR_4 = sum((rp_reg_4.rr.mu .- rp_reg_4.rr.y) .^ 2)
SSR_4 <= SSR_3
# Where we have piecewise constant data.
rp_21, rp_reg_21 = create_recursive_partitioning(dd2, y, x_variables, 3; rel_tol = 1e-3)
sum(abs.(evaluate(rp_21, dd2) .- dd2[:y])) < 1e-10

# MARS Spline
rp_1, rp_reg_1 = create_mars_spline(dd, y, x_variables, 3; rel_tol = 1e-3)
dd3 = deepcopy(dd)
dd3[:y] = 8.0 .* max.(0.0, dd3[:x] .+ 3.2) .+ 2.0 .* max.(0.0, dd3[:z] .+ 2.4) .+ 2.7 .* max.(0.0, dd3[:w] .- 1.2)
rp_32, rp_reg_32 = create_mars_spline(dd3, y, x_variables, 7; rel_tol = 1e-3)
sum(abs.(evaluate(rp_32, dd3) .- dd3[:y])) < 1e-08
dd4 = deepcopy(dd)
dd4[:y] = 1.0 .* max.(0.0, dd4[:x] .+ 3.2) .+ 2.0 .* max.(0.0, dd4[:z] .+ 2.4) .* max.(0.0, dd4[:w] .- 1.2)
rp_41, rp_reg_41 = create_mars_spline(dd4, y, x_variables, 7; rel_tol = 1e-3)
rp_41_RSS = sum(abs.(evaluate(rp_41, dd4) .- dd4[:y]))
rp_41_RSS < 0.05 # Here we dont get a great result. Because first it splits but w and then z so both are out of whack.
rp_41_5, rp_reg_41_5 = create_mars_spline(dd4, y, x_variables, 3; rel_tol = 1e-3)
reg_41_5_RSS = sum(abs.(evaluate(rp_41_5, dd4) .- dd4[:y]))

# Trim number of functions
trimmed_rp_41_fin, trimmed_rp_reg_41_fin = trim_mars_spline(dd4, y, rp_41; final_number_of_functions = 5)
reg_41_RSS_trimmed_to_five = sum(abs.(evaluate(trimmed_rp_41_fin, dd4) .- dd4[:y]))
reg_41_RSS_trimmed_to_five < reg_41_5_RSS # Because building a bigger model and then trimming should probably lead to lower RSS
reg_41_RSS_trimmed_to_five > rp_41_RSS # Because trimmed model should have higher RSS
# Maximum increase in RSS
trimmed_rp_41_maxinc, trimmed_rp_reg_41_maxinc = trim_mars_spline(dd4, y, rp_41; maximum_increase_in_RSS = 5.0)
trimmed_rp_41_maxrss, trimmed_rp_reg_41_maxrss = trim_mars_spline(dd4, y, rp_41; maximum_RSS = 100.0)
trimmed_rp_41_maxrss_trimmed = sum(abs.(evaluate(trimmed_rp_41_maxrss, dd4) .- dd4[:y]))
trimmed_rp_41_maxrss_trimmed < 100.0