text:
export
    MetaAgentCBS

# CBS submodule
module MetaAgentCBS

using ..CRCBS
using Parameters, LightGraphs, DataStructures

@with_kw struct State{S}
    states::Vector{S} = Vector{S}()
end
State(states::T) where {T<:Tuple} = State([states...])
Base.string(s::State) = string("(",prod(map(s->"$(string(s)),",s.states)),")")
@with_kw struct Action{A}
    actions::Vector{A} = Vector{A}()
end
Action(actions::T) where {T<:Tuple} = Action([actions...])
Base.hash(s::State{S}) where {S} = hash(s.states)
Base.:(==)(s1::S,s2::S) where {S<:State} = hash(s1) == hash(s2)
Base.string(a::Action) = string("(",prod(map(a->"$(string(a)),",a.actions)),")")

get_start_index(path::Path{MS,A,C}) where {S<:AbstractGraphState,MS<:State{S},A,C} = get_t(get_initial_state(path).states[1])

abstract type AbstractMetaEnv{S,A,T,C} <: AbstractLowLevelEnv{State{S},Action{A},MetaCostModel{T,C}} end
@with_kw struct LowLevelEnv{S,A,T,C<:AbstractCostModel{T},E<:AbstractLowLevelEnv{S,A,C}} <: AbstractMetaEnv{S,A,T,C}
    envs::Vector{E}             = Vector{E}()
    agent_idxs::Vector{Int}     = Vector{Int}()
    cost_model::MetaCostModel{T,C} = MetaCostModel(
        FullCostModel(sum,TravelTime()),length(envs))
end
# CRCBS.state_type(env::AbstractMetaEnv{S,A,T,C}) where {S,A,T,C} = State{S}
function CRCBS.compute_heuristic_cost(env::MetaAgentCBS.AbstractMetaEnv,cost,sp) where {E<:AbstractLowLevelEnv}
# function compute_heuristic_cost(m::MetaCostModel,::AbstractCostModel,env::E,cost,sp) where {E<:AbstractLowLevelEnv}
    m = get_cost_model(env)
    # h_cost = get_heuristic_cost(env,sp)
    costs = map(i->CRCBS.compute_heuristic_cost(
        get_envs(env)[i],
        cost.independent_costs[i],
        sp.states[i],
        # h_cost[i],
        ),1:m.num_agents)
    CRCBS.MetaCost(costs, CRCBS.aggregate_costs_meta(m.model, costs))
end

"""
    TeamMetaEnv{S,A,T,C,E}

Represents an env where a team of agents must move in rigid formation. The
function `CRCBS.get_possible_actions` must be overridden for all environments to
reflect this rigidity.
"""
@with_kw struct TeamMetaEnv{S,A,T,C<:AbstractCostModel{T},E<:AbstractLowLevelEnv{S,A,C}} <: AbstractMetaEnv{S,A,T,C}
    envs::Vector{E}             = Vector{E}()
    agent_idxs::Vector{Int}     = Vector{Int}()
    cost_model::MetaCostModel{T,C} = MetaCostModel(
        FullCostModel(sum,TravelTime()),length(envs))
end
get_envs(env::AbstractMetaEnv)                 = env.envs
get_agent_idxs(env::AbstractMetaEnv)           = env.agent_idxs
CRCBS.get_cost_model(env::AbstractMetaEnv)     = env.cost_model

function construct_meta_env(envs::Vector{E},idxs::Vector{Int}) where {S,A,T,C<:AbstractCostModel{T},E <: AbstractLowLevelEnv{S,A,C}}
    LowLevelEnv{S,A,T,C,E}(envs=envs,agent_idxs=idxs)
end
function construct_meta_env(envs::Vector{E},idxs::Vector{Int},c_model::C) where {S,A,T,C<:AbstractCostModel{T},E <: AbstractLowLevelEnv{S,A,C}}
    model = MetaCostModel(c_model,length(envs))
    LowLevelEnv{S,A,T,C,E}(envs=envs,agent_idxs=idxs,cost_model=model)
end
function construct_team_env(envs::Vector{E},idxs::Vector{Int}) where {S,A,T,C<:AbstractCostModel{T},E <: AbstractLowLevelEnv{S,A,C}}
    TeamMetaEnv{S,A,T,C,E}(envs=envs,agent_idxs=idxs)
end
function construct_team_env(envs::Vector{E},idxs::Vector{Int},c_model::C) where {S,A,T,C<:AbstractCostModel{T},E <: AbstractLowLevelEnv{S,A,C}}
    model = MetaCostModel(c_model,length(envs))
    TeamMetaEnv{S,A,T,C,E}(envs=envs,agent_idxs=idxs,cost_model=model)
end

function split_meta_path_node(path_node::PathNode{State{S},Action{A}}) where {S,A,C}
    return map(args->PathNode(args...), zip(
        get_s(path_node).states,
        get_a(path_node).actions,
        get_sp(path_node).states
        ))
end
function construct_meta_path_node(nodes::Vector{N}) where {N<:PathNode}
    PathNode(
        State(map(get_s,nodes)),
        Action(map(get_a,nodes)),
        State(map(get_sp,nodes)),
    )
end
construct_meta_path_node(nodes::Tuple) = construct_meta_path_node([nodes...])

# # TODO implement MetaPath
# @with_kw struct MetaPath{S,A,C} <: AbstractPath
#     paths       ::Vector{Path{S,A,C}}   = Vector{Path{S,A,C}}()
#     s0          ::State                 = get_s(get(path_nodes, 1, PathNode{S,A}()))
#     cost        ::MetaCost{C}           = MetaCost(map(get_cost,paths),typemax(C))
# end
# state_type(p::MetaPath{S,A,C}) where {S,A,C}    = State{S}
# action_type(p::MetaPath{S,A,C}) where {S,A,C}   = Action{A}
# cost_type(p::MetaPath{S,A,C}) where {S,A,C}     = MetaCost{C}
# Base.get(p::MetaPath,i,default=node_type(p)()) = construct_meta_path_node(map(p->get(p,i),p.paths))
# Base.getindex(p::MetaPath,i)    = construct_meta_path_node(map(p->getindex(p,i),p.paths))
# function Base.cat(p::P,x::N,args...) where {P<:MetaPath,N<:PathNode}
#     new_paths = typeof(p.paths)()
#     for (path,node) in zip(p.paths,split_path_node(x))
#         setindex!(path,node,i)
#         push!(new_paths,cat(path,node,args...))
#     end
#     P(s0=p.s0,paths=new_paths,cost=p.cost)
# end
# function Base.setindex!(p::P,x,i) where {P<:MetaPath}
#     for (path,node) in zip(p.paths,split_path_node(x))
#         setindex!(path,node,i)
#     end
# end
# Base.length(p::MetaPath) = maximum(map(length,p.paths))
# function Base.push!(p::MetaPath{S,A,C},n::PathNode{State{S},Action{A}}) where {S,A,C}
#     for (path,node) in zip(p.paths,split_path_node(n))
#         push!(path,node)
#     end
#     p
# end
# get_cost(p::MetaPath) = p.cost
# function set_cost!(p::MetaPath,meta_cost::MetaCost)
#     for (path,cost) in zip(p.paths,meta_cost.independent_costs)
#         set_cost!(path)
#     end
#     p.cost = meta_cost.cost
# end
# Base.copy(p::MetaPath) = MetaPath(s0=p.s0,paths=map(copy,p.paths),cost=p.cost)
# get_initial_state(path::MetaPath) = State(map(get_initial_state,path.paths))
# get_final_state(path::MetaPath) = State(map(get_final_state,path.paths))

function construct_meta_path(paths::Vector{P},cost) where {S,A,C,P<:Path{S,A,C}}
    @assert all(map(length,paths) .- maximum(map(length,paths)) .== 0) "$(map(length,paths))"
    starts = map(get_initial_state,paths)
    costs = map(get_cost, paths)
    path_nodes = Vector{PathNode{State{S},Action{A}}}()
    for nodes in zip(map(p->p.path_nodes,paths)...)
        push!(path_nodes, construct_meta_path_node(nodes))
    end
    meta_path = Path(
        path_nodes = path_nodes,
        s0=State(starts),
        cost = MetaCost(costs,cost)
    )
end
function CRCBS.build_env(solver::MetaAgentCBS_Solver, mapf, node, i)
    group = node.solution.group_idxs[i]
    construct_meta_env(
        [build_env(mapf, node, j) for j in group],
        group,
        get_cost_model(mapf.env)
        )
end
CRCBS.get_start(mapf::MAPF,env::AbstractMetaEnv,i) = State([get_start(mapf,e,j) for (e,j) in zip(get_envs(env), get_agent_idxs(env))])
################################################################################
################################################################################
################################################################################

"""
    `split_path(path::Path{State{S},Action{A}}) where {S,A}`

    Helper for breaking a meta-agent path into a set of single-agent paths
"""
function split_path(path::Path{State{S},Action{A},MetaCost{C}}) where {S,A,C}
    N = length(get_s(get_path_node(path, 1)).states)
    paths = [Path{S,A,C}(s0=path.s0.states[i],cost=get_cost(path).independent_costs[i]) for i in 1:N]
    for t in 1:length(path)
        path_node = get_path_node(path, t)
        for (p,n) in zip(paths,split_meta_path_node(path_node))
            push!(p,n)
        end
    end
    return paths
end
function CRCBS.set_solution_path!(solution::MetaSolution, meta_path, group_idx)
    for (idx,path) in zip(solution.group_idxs[group_idx],split_path(meta_path))
        set_solution_path!(solution.solution, path, idx)
    end
end
function CRCBS.set_path_cost!(solution::MetaSolution, meta_cost, group_idx)
    for (idx,cost) in zip(solution.group_idxs[group_idx],meta_cost.independent_costs)
        set_path_cost!(solution.solution, cost, idx)
    end
end
function CRCBS.get_heuristic_cost(env::E,state::State) where {E<:AbstractMetaEnv}
    c = [get_heuristic_cost(e,s) for (e,s) in zip(get_envs(env), state.states)]
end
function CRCBS.states_match(s1::State,s2::State)
    for (state1, state2) in zip(s1.states,s2.states)
        if !states_match(state1,state2)
            return false
        end
    end
    return true
end
function CRCBS.is_goal(env::E,state::State) where {E<:AbstractMetaEnv}
    for (e,s) in zip(get_envs(env), state.states)
        if !is_goal(e,s)
            return false
        end
    end
    return true
end
CRCBS.wait(s::State) = Action(map(CRCBS.wait, s.states))
CRCBS.wait(env::AbstractMetaEnv,s::State) = Action(map(args->CRCBS.wait(args...),zip(get_envs(env),s.states)))
# get_possible_actions
struct ActionIter{A}
    action_lists::Vector{Vector{A}}
end
const ActionIterState = Vector{Int}
function Base.iterate(it::ActionIter)
    iter_state = ActionIterState(ones(Int, length(it.action_lists)))
    if length(iter_state) > 0
        iter_state[end] = 0
        return iterate(it, iter_state)
    else
        return nothing
    end
end
function Base.iterate(iter::ActionIter,iter_state::ActionIterState)
    iter_state = deepcopy(iter_state)
    i = length(iter_state)
    while i > 0
        if iter_state[i] < length(iter.action_lists[i])
            iter_state[i] += 1
            action = Action([iter.action_lists[j][idx] for (j,idx) in enumerate(iter_state)])
            return action, iter_state
        else
            iter_state[i] = 1
            i -= 1
        end
    end
    return nothing
end
function CRCBS.get_possible_actions(env::E,s::State) where {E<:LowLevelEnv}
    return ActionIter([collect(get_possible_actions(env_,s_)) for (env_,s_) in zip(get_envs(env), s.states)])
end
Base.length(iter::ActionIter) = product([length(a) for a in iter.action_lists])
CRCBS.get_next_state(s::State,a::Action) = State(
    [get_next_state(s_,a_) for (s_,a_) in zip(s.states,a.actions)]
    )
CRCBS.get_next_state(env::E,s::State,a::Action) where {E<:AbstractMetaEnv} = State(
    [get_next_state(e_,s_,a_) for (e_,s_,a_) in zip(get_envs(env),s.states,a.actions)]
    )
CRCBS.get_transition_cost(env::E,s::State,a::Action,sp::State) where {E<:AbstractMetaEnv} = [
    get_transition_cost(e_,s_,a_,sp_) for (e_,s_,a_,sp_) in zip(get_envs(env),s.states,a.actions,sp.states)
]
function CRCBS.violates_constraints(env::E, s::State, a::Action, sp::State) where {E<:AbstractMetaEnv}
    for (i, (e_,s_,a_,sp_)) in enumerate(zip(get_envs(env),s.states,a.actions,sp.states))
        if violates_constraints(e_,s_,a_,sp_)
            return true
        end
        # check for self conflict
        P1 = PathNode(s_,a_,sp_)
        for (j, (s2_,a2_,sp2_)) in enumerate(zip(s.states,a.actions,sp.states))
            if i != j
                P2 = PathNode(s2_,a2_,sp2_)
                if detect_state_conflict(P1,P2) || detect_action_conflict(P1,P2)
                    return true
                end
            end
        end
    end
    return false
end

### Reservation Table Interface
function CRCBS.reserve!(table::PIBTReservationTable,env::MetaAgentCBS.AbstractMetaEnv,s,a,sp,t=-1)
    valid = true
    for args in zip(get_envs(env),s.states,a.actions,sp.states)
        valid &= reserve!(table,args...)
    end
    return valid
end
function CRCBS.is_reserved(table::PIBTReservationTable,env::MetaAgentCBS.AbstractMetaEnv,s,a,sp,t=-1)
    for args in zip(get_envs(env),s.states,a.actions,sp.states)
        if is_reserved(table,args...)
            return true
        end
    end
    return false
end
function CRCBS.create_reservations(env::AbstractMetaEnv,s,a,sp,t=-1)
    vcat(map(args->create_reservations(args...,t), zip(get_envs(env),s.states,a.actions,sp.states))...)
end

end


text:
module DistributedStwdLDA

using Distributions, StatsBase, Statistics, SpecialFunctions
using LinearAlgebra, Distributed
using SharedArrays,DataFrames,CSV
using FileIO,MCMCChains,ProgressMeter
using Random: Random

export Word,getCorp,getSynthData
export getBitarrays
export initShared
export mat2docind, vec2procind, initterms!,initcorp,initcorpBeta,initcorpChain
export getZ,getJ,getW
export nj2token,nw2token
export chunkSharedSV!,sampleSharedSV!,runChainSharedSV
export splitApprox,epochIndex,filterTokens,getIndRand
export findGenePrefix,filterCells,filterGenes,getParamDict,filterPbulkparams,dsLDA_E_step,procsN,getTheta!

include("corpus.jl") # corpus data and methods to get bit arrays from
include("bitarrays.jl") # methods on bit arrays
include("partition.jl") # partitioning and allocation
include("sampling.jl") # sampler code
include("chains.jl") # sampling diagnostics
include("durian.jl") # dslda stuff for durian
end

ext: jl


text:
const thesevecnames = @compat Dict(
    :ncurve => :curve_names,
    :nmultimesh=>:multimesh_names,
    :nmultimeshadj=>:multimeshadj_names,
    :nmultivar=>:multivar_names,
    :nmultimat=>:multimat_names,
    :nmultimatspecies=>:multimatspecies_names,
    :ncsgmesh=>:csgmesh_names,
    :ncsgvar=>:csgvar_names,
    :ndefvars=>:defvars_names,
    :nqmesh=>:qmesh_names,
    :nqvar=>:qvar_names,
    :nucdmesh=>:ucdmesh_names,
    :nucdvar=>:ucdvar_names,
    :nptmesh=>:ptmesh_names,
    :nptvar=>:ptvar_names,
    :nmat=>:mat_names,
    :nmatspecies=>:matspecies_names,
    :nvar=>:var_names,
    :nobj=>:obj_names,
    :ndir=>:dir_names,
    :narray=>:array_names,
    :nmrgtree=>:mrgtree_names,
    :ngroupelmap=>:groupelmap_names,
    :nmrgvar=>:mrgvar_names)

type JToc
  curve_names::Vector{String}
  multimesh_names::Vector{String}
  multimeshadj_names::Vector{String}
  multivar_names::Vector{String}
  multimat_names::Vector{String}
  multimatspecies_names::Vector{String}
  csgmesh_names::Vector{String}
  csgvar_names::Vector{String}
  defvars_names::Vector{String}
  qmesh_names::Vector{String}
  qvar_names::Vector{String}
  ucdmesh_names::Vector{String}
  ucdvar_names::Vector{String}
  ptmesh_names::Vector{String}
  ptvar_names::Vector{String}
  mat_names::Vector{String}
  matspecies_names::Vector{String}
  var_names::Vector{String}
  obj_names::Vector{String}
  dir_names::Vector{String}
  array_names::Vector{String}
  mrgtree_names::Vector{String}
  groupelmap_names::Vector{String}
  mrgvar_names::Vector{String}

  namemapping::Dict{Symbol,Symbol}
  function JToc(dbfile::Ptr{DBfile})
    inttoc = new()

    rawnames = unsafe_load(DBGetToc(dbfile))
    thesenums = [keyname => Int32(rawnames.(keyname)) for keyname in keys(thesevecnames)]

    for keyname in keys(thesevecnames)
      if thesenums[keyname] == 0
        inttoc.(thesevecnames[keyname]) = Array(String, 0,)
      else
        inttoc.(thesevecnames[keyname]) = map(bytestring,pointer_to_array(rawnames.(thesevecnames[keyname]), thesenums[keyname]))
      end
  end
    inttoc.namemapping = thesevecnames
  return inttoc
  end #function
end #type

function Base.show(io::IO, tocin::JToc)
  for name in values(tocin.namemapping)
    print(io, "Has $(length(tocin.name)) $(name): ", join(tocin.(name), ", "))
  end
end

## Functionality related to file names introspection
function Base.names(dbfile::Ptr{DBfile})
  tocin = JToc(dbfile)
  return tocin.curve_names
end
function Base.names(dbfile::silofile)
  if isvalid(dbfile)
    namesarr = names(dbfile.dbfile)
  else
    error("Passed silofile object is not valid.")
  end
  namesarr
end
function Base.names(file_name::String)
  dbfile = siloopen(file_name, "r")
  namesarr = names(dbfile)
  siloclose(dbfile)
  return namesarr
end
#
# export silofile, isvalid, silo_is_silofile, siloopen, siloclose
#
## Thin Julia type wrapper around the internal DBfile pointer
# tracks whether pointer is internally allocated, and the associated
# filename.
type silofile
  dbfile::Ptr{DBfile}
  filemode::String
  file_name::String
  ptropen::Bool
end
# function silofile(dbfile::Ptr{DBfile}, filemode::String, file_name::String, ptropen::Bool, toclose::Bool)
#   f = silofile(dbfile, filemode, file_name, ptropen)
#   toclose && finalizer(f, siloclose)
#   return f
# end

isvalid(f::silofile) = f.ptropen
Base.convert(::Type{Ptr}, f::silofile) = f.dbfile
Base.show(io::IO, f::silofile) = isvalid(f) ? print(io, "Silo file: ", f.file_name, " open with mode ", f.filemode) : print(io, "Closed Silo file: ", f.file_name)
#
#
#
## High-level Julia open/close functionality
function siloopen(file_name::String,
                  mode::String;
                  filetype::Int = DB_PDB,
                  fileinfo::String = "",
                  filecreatemode::Int = DB_CLOBBER,
                  filecreatetarget::Int = DB_LOCAL,)
  if mode=="r"
    if !isfile(file_name)
      error("File $file_name does not exist.")
    end
    if !is_silofile(file_name)
      error("File $file_name is not a Silo file.")
    end
    dbfile = DBOpenReal(file_name, filetype, DB_READ)
  elseif mode=="w"
    if isfile(file_name)
      rm(file_name)
    end
    dbfile =  DBCreate(file_name, filecreatemode, filecreatetarget, fileinfo, Silo.DB_PDB)
  elseif mode=="a"
    if !isfile(file_name)
      dbfile = DBCreate(file_name, filecreatemode, filecreatetarget, fileinfo, Silo.DB_PDB)
    else
      dbfile = DBOpen(file_name, filetype, DB_APPEND)
    end
  else
    error("Mode $mode not defined for Silo files.")
  end
  if dbfile == C_NULL
    error("File $file_name could not be opened in mode $mode.")
  end
  return silofile(dbfile, mode, file_name, true)
end

function siloclose(dbfile::Ptr{DBfile})
  DBClose(dbfile)
end
function siloclose(dbfile::silofile)
  DBClose(dbfile.dbfile)
  dbfile.ptropen = false
end

include("write.jl")

include("read.jl")


text:
################################################################################
#
#    NfOrd/TorsionUnits.jl : Torsion units in generic number field orders
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# 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.
#
# 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.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

################################################################################
#
#  Torsion unit test
#
################################################################################

doc"""
***
    is_torsion_unit(x::NfOrdElem, checkisunit::Bool = false) -> Bool

> Returns whether $x$ is a torsion unit, that is, whether there exists $n$ such
> that $x^n = 1$.
>
> If `checkisunit` is `true`, it is first checked whether $x$ is a unit of the
> maximal order of the number field $x$ is lying in.
"""
function is_torsion_unit(x::NfOrdElem, checkisunit::Bool = false)
  return is_torsion_unit(x.elem_in_nf, checkisunit)
end

################################################################################
#
#  Order of a single torsion unit
#
################################################################################

doc"""
***
    torsion_unit_order(x::NfOrdElem, n::Int)

> Given a torsion unit $x$ together with a multiple $n$ of its order, compute
> the order of $x$, that is, the smallest $k \in \mathbb Z_{\geq 1}$ such
> that $x^k = 1$.
>
> It is not checked whether $x$ is a torsion unit.
"""
function torsion_unit_order(x::NfOrdElem, n::Int)
  return torsion_unit_order(x.elem_in_nf, n)
end

################################################################################
#
#  Torsion units via lattice enumeration
#
################################################################################

doc"""
***
    torsion_units(O::NfOrd) -> Array{NfOrdElem, 1}

> Given an order $O$, compute the torsion units of $O$.
"""
function torsion_units(O::NfOrd)
  ar, g = _torsion_units(O)
  return ar
end

doc"""
***
    torsion_units(O::NfOrd) -> NfOrdElem

> Given an order $O$, compute a generator of the torsion units of $O$.
"""
function torsion_units_gen(O::NfOrd)
  ar, g = _torsion_units(O)
  return g
end

function _torsion_units(O::NfOrd)
   if isdefined(O, :torsion_units)
    return O.torsion_units
  end

  z = _torsion_units_lifting(O)
  O.torsion_units = z
  return O.torsion_units
end

function _torsion_units_lattice_enum(O::NfOrd)
  n = degree(O)
  K = nf(O)
  r1, r2 = signature(K)
  B = basis(O)

  if r1 > 0
    return [ O(1), -O(1) ], -O(1)
  end

  function _t2_pairing(x, y, p)
    local i
    v = minkowski_map(x, p)
    w = minkowski_map(y, p)

    t = zero(parent(v[1]))

    for i in 1:r1
      t = t + v[i]*w[i]
    end

    for i in (r1 + 1):(r1 + 2*r2)
      t = t + v[i]*w[i]
    end

    return t
  end

  p = 64

  gram_found = false

  could_enumerate = false

  A = ArbField(p)
  M = ArbMatSpace(A, n, n)()

  while true
    A = ArbField(p)
    M = ArbMatSpace(A, n, n)()

    gram_found = true

    @vprint :UnitGroup 1 "Computing Gram matrix with precision $p ... \n"

    for i in 1:n, j in 1:n
      M[i,j] = _t2_pairing(B[i], B[j], p)
      if !isfinite(M[i, j])
        p = 2*p
        gram_found = false
        break
      end
    end

    if gram_found
      break
    end
  end

  @vprint :UnitGroup 1 "Gram matrix is $M\n"

  @vprint :UnitGroup 1 "Enumerating elements with T_2 bounded by $n ... \n"
  l = enumerate_using_gram(M, A(n))

  R = Array{NfOrdElem, 1}()

  for i in l
    if O(i) == zero(O)
      continue
    end
    if is_torsion_unit(O(i))
      push!(R, O(i))
    end
  end

  i = 0

  for i in 1:length(R)
    if torsion_unit_order(R[i], length(R)) == length(R)
      break
    end
  end

  return R, R[i]
end

################################################################################
#
#  Torsion units via root lifting
#
################################################################################

function _torsion_group_order_divisor(O::NfOrd, N::Int = 5)
  p = 1
  m = fmpz(0)
  m_old = fmpz(0)
  stable = 0

  while true
    p = next_prime(p)
    if isramified(O, p)
      continue
    end
    lp = prime_decomposition(O, p)
    m_new = m_old
    for (P, e) in lp
      m_new = gcd(m_new, norm(P) - 1)
    end

    if m_new == m_old
      stable += 1
    end

    if stable == 5
      return m_new
    end

    m_old = m_new
  end
end

function _torsion_units_lifting(O::NfOrd)
  m = _torsion_group_order_divisor(O)
  Oy, y = O["y"]
  f = y^Int(m) - 1
  R = _roots_hensel(f)

  i = 1
  for i in 1:length(R)
    if torsion_unit_order(R[i], length(R)) == length(R)
      break
    end
  end

  return R, R[i]
end


text:
# This script is modified from:
#   `https://github.com/evetion/GeoArrays.jl/blob/master/src/plot.jl`
# Copyright (c) 2018 Maarten Pronk, MIT license

using RecipesBase

@recipe function f(ga::AbstractGeoArray; band=1, fact=nothing)
    xflip --> false
    yflip --> false
    aspect_ratio --> 1
    seriestype := :heatmap
    seriescolor := :viridis

    # 默认的分辨率为全球0.25°
    if fact === nothing
        n = prod(size(ga.A)[1:2])
        cellsize = 0.5
        fact = ceil(Int, sqrt(n / (360 / cellsize * 180 / cellsize)))
    end

    if fact > 1
        println("[st_plot]: fact = $fact\n")
        ga = resample(ga, fact)
    end
    # coords = st_coordinates(ga)
    # @show coords
    x = st_coordinates(ga, :x)
    y = st_coordinates(ga, :y)
    # x = map(x->x[1], coords[:, 1])
    # y = map(x->x[2], coords[end, :])
    z = ga.A[:, :, band]'

    # Can't use x/yflip as x/y coords
    # have to be sorted for Plots
    if ga.f.linear[1] < 0
        z = reverse(z, dims=2)
        reverse!(x)
    end
    if ga.f.linear[4] < 0
        z = reverse(z, dims=1)
        reverse!(y)
    end

    xlims --> (extrema(x))
    ylims --> (extrema(y))

    x, y, z
end


text:
module WikstromTerelius

using Infiltrator

include("generators.jl")
include("verifier.jl")

export verify, prove

end # module


text:
module ProtoBufTestApis
using ProtoBuf
using Compat
using Base.Test
import ProtoBuf.meta

print_hdr(tname) = println("testing $tname...")

type TestType
    a::AbstractString
    b::Bool
    TestType() = (o=new(); fillunset(o); o)
end #type TestType
meta(t::Type{TestType}) = meta(t, Symbol[:a], Int[], Dict{Symbol,Any}())

function test_apis()
    t = TestType()

    @test !has_field(t, :a)
    @test !has_field(t, :b)

    @test false == try get_field(t, :a); true; catch; false; end

    set_field!(t, :b, true)
    @test has_field(t, :b)
    @test (get_field(t, :b) == true)

    @test !isinitialized(t)
    t.a = "hello"
    @test !isinitialized(t)
    set_field!(t, :a, "hello world")
    @test isinitialized(t)
    @test (get_field(t, :a) ==  "hello world")

    clear(t, :b)
    @test isinitialized(t)
    clear(t)
    @test !isinitialized(t)

    t = protobuild(TestType, @compat Dict(:a => "hello", :b => false))
    @test t.a == "hello"
    @test t.b == false
end
end # module ProtoBufTestApis

print_hdr("utility api methods")
ProtoBufTestApis.test_apis()


text:
using Ramnet.Encoders
using Ramnet

using Random

using ..Buffers: MultiStepDynamicBuffer, add!
using ..Environments

mutable struct MCDifferentialDiscriminatorAgent{O <: AbstractVector,A <: Real,T <: Real,E <: AbstractEncoder{T}} <: AbstractAgent{O,A}
    actions::UnitRange{A}
    n::Int
    obs_size::Int
    α::Float64
    γ::Float64
    ϵ::Float64
    epochs::Int
    encoder::E
    Q̂::Dict{A,DifferentialDiscriminator{T,E}}
    buffer::MultiStepDynamicBuffer{O,A}
    observation::Union{Nothing,O}
    done::Bool
    action::Union{Nothing,A}
    rng::MersenneTwister

  # TODO: Validate actions and encoder and ϵ
    function MCDifferentialDiscriminatorAgent{O,A,T,E}(actions, n, obs_size, α, γ, ϵ, epochs, encoder::E; seed::Union{Nothing,Int}=nothing) where {O <: AbstractVector,A <: Real,T <: Real,E <: AbstractEncoder{T}}
        !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
        rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

        n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))
        epochs < 1 && throw(DomainError(epochs, "Number of epochs must be at least 1"))

        Q̂ = Dict(a => DifferentialDiscriminator(obs_size, n, encoder; α) for a in actions)
        buffer = MultiStepDynamicBuffer{O,A}()

        new(actions, n, obs_size, α, γ, ϵ, epochs, encoder, Q̂, buffer, nothing, false, nothing, rng)
    end
end

function MCDifferentialDiscriminatorAgent(
  env::AbstractEnvironment{O,A}, encoder::E; n, α::Float64=1.0, γ::Float64=1.0, ϵ::Float64=0.01, epochs::Int=1, seed::Union{Nothing,Int}=nothing) where {O <: AbstractVector,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    MCDifferentialDiscriminatorAgent{O,A,T,E}(action_set(env), n, observation_length(env), α, γ, ϵ, epochs, encoder; seed)
end

function observe!(agent::MCDifferentialDiscriminatorAgent{O,A,T,E}, observation::O) where {O,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    agent.observation = observation
    agent.done = false

    nothing
end

# TODO: All observe! methods are the same for all agents. Generalize.
# TODO: This method does not need to take in the action
function observe!(agent::MCDifferentialDiscriminatorAgent{O,A,T,E}, action::A, reward::Float64, observation::O, done::Bool) where {O <: AbstractVector,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    add!(agent.buffer, agent.observation, agent.action, reward)

    agent.observation = observation
    agent.done = done

    nothing
end

function update!(agent::MCDifferentialDiscriminatorAgent{O,A,T,E}) where {O <: AbstractVector,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    if agent.done
        G = 0.0
        len = length(agent.buffer)
        Gs = Vector{Float64}(undef, len)
        obs = O[]
        idx_by_action = Dict{A,Vector{Int}}()

        for (i, transition) in Iterators.enumerate(Iterators.reverse(agent.buffer))
            G = transition.reward + agent.γ * G
            Gs[i] = G
            # push!(idx_by_action[transition.action], i)
            push!(get!(idx_by_action, transition.action, Int[]), i)
            push!(obs, transition.observation)
          # train!(
          #     agent.Q̂[transition.action],
          #     transition.observation,
          #     G
          # )
        end

        # @show Gs
        # @show idx_by_action
        for a in keys(idx_by_action)
            train!(
              agent.Q̂[a],
              hcat(obs[idx_by_action[a]]...),
              Gs[idx_by_action[a]];
              epochs=agent.epochs
            )
        end

        reset!(agent.buffer)
    end

    nothing
end

function Agents.reset!(agent::MCDifferentialDiscriminatorAgent{O,A,T,E}) where {O <: AbstractVector,A <: Real,T <: Real,E <: AbstractEncoder{T}}
    for action in agent.actions
        agent.Q̂[action] = DifferentialDiscriminator(agent.obs_size, agent.n, agent.encoder; α=agent.α)
    end

    reset!(agent.buffer)

    nothing
end


text:
@testset "HLADataset(::Val{:Test})" begin
    @test Escape.HLADataset("Test") isa Escape.HLADataset
end


text:
module Presentation

import Pandoc
import REPL
import REPL: Terminals
import Base: read

using Crayons
using Highlights
using Highlights.Tokens
using Highlights.Lexers
using TerminalExtensions
using Markdown

export render, next, previous, current_slide

abstract type PandocMarkdown end
abstract type JuliaMarkdown end

include("lexers.jl")

include("utils.jl")
include("slideshow.jl")
include("renderer.jl")

TERMINAL = nothing # Contains reference to the built in Terminal
SLIDES = nothing # Contains reference to a `Slides` object

render() = render(SLIDES)
next() = next(SLIDES)
previous() = previous(SLIDES)
current_slide() = current_slide(SLIDES)
filename() = filename(SLIDES)

function __init__()
    global TERMINAL
    TERMINAL = REPL.Terminals.TTYTerminal(get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)
end

end


text:
using RandomTweaks
using Test

@info "Testing Tweaks"
include("load_test.jl")
include("curve_test.jl")


text:
#!/usr/bin/env julia

using Primitiv
using Base.Test

this_file = basename(@__FILE__)

function test_dir(dir)
    jl_files = sort(filter(f -> ismatch(r"^.+\.jl$", f), readdir(dir)),
                    by = fn -> stat(joinpath(dir, fn)).mtime)
    map(reverse(jl_files)) do file
        file == this_file && return
        include(joinpath(dir, file))
    end
end

@testset "Primitiv Test" begin
    test_dir(dirname(@__FILE__))
end


text:
"""
    convertvideo(file; pathout, filenameout)

Convert the tar-file obtained from a simulation to an mp4-video.

# Available kwargs
* `pathout`:  The output path.
* `filename`: The new file name.
"""
function convertvideo(file; pathout::String=splitdir(file)[1], filenameout::String=splitext(splitdir(file)[2])[1])
    file = abspath(file)
    pathout = abspath(pathout)
    filenameout = filenameout*".mp4"
    println("Creating " * filenameout * ".")
    MeshCat.convert_frames_to_video(file, joinpath(pathout, filenameout), overwrite=true)
    return
end

# Fix for Windows video creation
# function convert_frames_to_video(tar_file_path::AbstractString, output_path::AbstractString="output.mp4"; framerate=60, overwrite=false)
#     output_path = abspath(output_path)
#     if !isfile(tar_file_path)
#         error("Could not find the input file $tar_file_path")
#     end
#     if isfile(output_path) && !overwrite
#         error("The output path $output_path already exists. To overwrite that file, you can pass `overwrite=true` to this function")
#     end

#     mktempdir() do tmpdir
#         run(unpack_cmd(tar_file_path, tmpdir, ".tar", nothing))
#         cmd = ["-r", string(framerate), "-i", "%07d.png", "-vcodec", "libx264", "-preset", "slow", "-pix_fmt", "yuv420p", "-crf", "18", "-filter_complex", "[0]split=2[bg][fg];[bg]drawbox=c=white@1:replace=1:t=fill[bg];[bg][fg]overlay=format=auto"]
#         if overwrite
#             push!(cmd, "-y")
#         end
#         push!(cmd, output_path)

#         cd(tmpdir) do
#             FFMPEG.exe(cmd...)
#         end
#     end
#     @info("Saved output as $output_path")
#     return output_path
# end


text:
# Functional implementation of https://egraphs-good.github.io/
# https://dl.acm.org/doi/10.1145/3434304

using DataStructures

"""
Abstract type representing an [`EGraph`](@ref) analysis,
attaching values from a join semi-lattice domain to
an EGraph
"""
abstract type AbstractAnalysis end

const ClassMem = Dict{Int64,Vector{Any}}
const HashCons = Dict{Any,Int64}
const Parent = Tuple{Any,Int64} # parent enodes and eclasses
const ParentMem = Dict{Int64,Vector{Parent}}
const AnalysisData = Dict{Int64,Any}
const Analyses = Vector{AbstractAnalysis}

"""
A concrete type representing an [`EGraph`].
See the [egg paper](https://dl.acm.org/doi/pdf/10.1145/3434304)
for implementation details
"""
mutable struct EGraph
    """stores the equality relations over e-class ids"""
    U::IntDisjointSets
    """map from eclass id to eclasses"""
    M::ClassMem             #
    H::HashCons             # hashcons
    parents::ParentMem
    """worklist for ammortized upwards merging"""
    dirty::Vector{Int64}
    root::Int64
    """A vector of analyses associated to the EGraph"""
    analyses::Analyses
end

EGraph() = EGraph(
    IntDisjointSets(0),
    ClassMem(),
    HashCons(),
    ParentMem(),
    Vector{Int64}(),
    0,
    Analyses(),
)

function EGraph(e)
    G = EGraph()
    rootclass = addexpr!(G, e)
    G.root = rootclass.id
    G
end

function addparent!(G::EGraph, a::Int64, parent::Parent)
    @assert isenode(parent[1])
    if !haskey(G.parents, a)
        G.parents[a] = [parent]
    else
        union!(G.parents[a], [parent])
    end
end

"""
Inserts an e-node in an [`EGraph`](@ref)
"""
function add!(G::EGraph, n)::EClass
    @debug("adding ", n)
    canonicalize!(G.U, n)
    if haskey(G.H, n)
        return find(G, G.H[n]) |> EClass
    end
    @debug(n, " not found in H")
    id = push!(G.U) # create new singleton eclass
    !haskey(G.parents, id) && (G.parents[id] = [])
    getfunargs(n) .|> x -> addparent!(G, x.id, (n, id))

    G.H[n] = id
    G.M[id] = [n]

    # make analyses for new enode
    for analysis ∈ G.analyses
        if !islazy(analysis)
            analysis[id] = make(analysis, n)
            modify!(analysis, id)
        end
    end

    return EClass(id)
end
add!(g::EGraph, c::EClass)::EClass = EClass(find(g, c))

"""
Recursively traverse an [`Expr`](@ref) and insert terms into an
[`EGraph`](@ref). If `e` is not an [`Expr`](@ref), then directly
insert the literal into the [`EGraph`](@ref).
"""
function addexpr!(G::EGraph, e)::EClass
    e = cleanast(e)
    df_walk((x -> add!(G, x)), e; skip_call = true)
end

function mergeparents!(G::EGraph, from::Int64, to::Int64)
    !haskey(G.parents, from) && (G.parents[from] = []; return)
    !haskey(G.parents, to) && (G.parents[to] = [])

    union!(G.parents[to], G.parents[from])
    delete!(G.parents, from)
end

"""
Given an [`EGraph`](@ref) and two e-class ids, set
the two e-classes as equal.
"""
function Base.merge!(G::EGraph, a::Int64, b::Int64)::Int64
    id_a = find(G, a)
    id_b = find(G, b)
    id_a == id_b && return id_a
    id_u = union!(G.U, id_a, id_b)

    @debug "merging" id_a id_b

    from, to = if (id_u == id_a)
        id_b, id_a
    elseif (id_u == id_b)
        id_a, id_b
    else
        error("egraph invariant maintenance error")
    end

    push!(G.dirty, id_u)

    clean(t) = begin
        delete!(G.H, t)
        canonicalize!(G.U, t)
        G.H[t] = to
        t
    end

    G.M[from] = map(clean, G.M[from])
    G.M[to] = map(clean, G.M[to])
    G.M[to] = G.M[from] ∪ G.M[to]

    if from == G.root
        G.root = to
    end

    delete!(G.M, from)
    mergeparents!(G, from, to)
    for analysis ∈ G.analyses
        if haskey(analysis, from) && haskey(analysis, to)
            analysis[to] = join(analysis, analysis[from], analysis[to])
            delete!(analysis, from)
        end
    end

    return id_u
end

"""
Returns the canonical e-class id for a given e-class.
"""
find(G::EGraph, a::Int64)::Int64 = find_root!(G.U, a)
find(G::EGraph, a::EClass)::Int64 = find_root!(G.U, a.id)

"""
This function restores invariants and executes
upwards merging in an [`EGraph`](@ref). See
the [egg paper](https://dl.acm.org/doi/pdf/10.1145/3434304)
for more details.
"""
function rebuild!(G::EGraph)
    while !isempty(G.dirty)
        todo = unique([find(G, id) for id ∈ G.dirty])
        empty!(G.dirty)
        foreach(todo) do x
            repair!(G, x)
        end
    end

    if G.root != 0
        G.root = find(G, G.root)
    end
end

function repair!(G::EGraph, id::Int64)
    id = find(G, id)
    @debug "repairing " id

    for (p_enode, p_eclass) ∈ G.parents[id]
        #old_id = G.H[p_enode]
        #delete!(G.M, old_id)
        delete!(G.H, p_enode)
        @debug "deleted from H " p_enode
        n = canonicalize(G.U, p_enode)
        n_id = find(G, p_eclass)
        G.H[n] = n_id
    end

    new_parents = OrderedDict{Any,Int64}()

    for (p_enode, p_eclass) ∈ G.parents[id]
        canonicalize!(G.U, p_enode)
        # deduplicate parents
        if haskey(new_parents, p_enode)
            @debug "merging classes" p_eclass (new_parents[p_enode])
            merge!(G, p_eclass, new_parents[p_enode])
        end
        new_parents[p_enode] = find(G, p_eclass)
    end
    G.parents[id] = collect(new_parents) .|> Tuple
    @debug "updated parents " id G.parents[id]

    # Analysis invariant maintenance
    for an ∈ G.analyses
        haskey(an, id) && modify!(an, id)
        # modify!(an, id)
        id = find(G, id)
        for (p_enode, p_eclass) ∈ G.parents[id]
            # p_eclass = find(G, p_eclass)
            if !islazy(an) && !haskey(an, p_eclass)
                an[p_eclass] = make(an, p_enode)
            end
            if haskey(an, p_eclass)
                new_data = join(an, an[p_eclass], make(an, p_enode))
                if new_data != an[p_eclass]
                    an[p_eclass] = new_data
                    push!(G.dirty, p_eclass)
                end
            end
        end
    end
end

"""
Recursive function that traverses an [`EGraph`](@ref) and
returns a vector of all reachable e-classes from a given e-class id.
"""
function reachable(g::EGraph, id::Int64; hist=Int64[])
    id = find(g, id)
    hist = hist ∪ [id]
    for n ∈ g.M[id]
        if n isa Expr
            for child_eclass ∈ getfunargs(n)
                c_id = child_eclass.id
                if c_id ∉ hist
                    hist = hist ∪ reachable(g, c_id; hist=hist)
                end
            end
        end
    end

    return hist
end


text:
using ImageView, Images   # TestImages  just a reminder poop stinks!
using Gtk.ShortNames

grid, frames, canvases = canvasgrid((1,1))  # 1 row, 2 columns
img = fill(RGBA(0.65625,0.796875,1,1), 500, 800)
imshow(canvases[1,1], img)
win = Window(grid)
showall(win)
# ==============================================================================
# This draws a pixel with Alpha chanel
# ==============================================================================
function pixel(image, X, Y, color )
    a = alpha(color)
    width, height = size(image)
    if X>0 && X < width && Y>0 && Y < height
        image[ X,Y ] = image[ X,Y ]*(1-a) + (color*a)
    end
end
# ==============================================================================
# This draws a rectrangle to a picture!
# ==============================================================================
function rect(image, left, top, width, height, color)
    for y in 1:height, x in 1:width
        pixel(image, top+y, left+x, color )
    end
end
# ==============================================================================
# This draws random rectrangles!
# ==============================================================================
function RandomPaint()
    left, top = rand( 1:800), rand( 1:500)
    rect(img, left, top,
              rand( 1:300), rand( 1:300),
              rand(RGBA{Float32})  )
    imshow(canvases[1,1], img)
end

u1 = Timer(u1 -> RandomPaint(), 0.01, 0.1)

#imshow(canvases[1,2], testimage("mandrill"))
#win = Window(grid)
# showall(win)

max_stars_repo_name: itsKriptunium/MaelEngine

max_issues_repo_name: itsKriptunium/MaelEngine

max_forks_repo_name: itsKriptunium/MaelEngine


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""TokenReviewStatus is the result of the token authentication request.

    IoK8sApiAuthenticationV1beta1TokenReviewStatus(;
        audiences=nothing,
        authenticated=nothing,
        error=nothing,
        user=nothing,
    )

    - audiences::Vector{String} : Audiences are audience identifiers chosen by the authenticator that are compatible with both the TokenReview and token. An identifier is any identifier in the intersection of the TokenReviewSpec audiences and the token's audiences. A client of the TokenReview API that sets the spec.audiences field should validate that a compatible audience identifier is returned in the status.audiences field to ensure that the TokenReview server is audience aware. If a TokenReview returns an empty status.audience field where status.authenticated is \"true\", the token is valid against the audience of the Kubernetes API server.
    - authenticated::Bool : Authenticated indicates that the token was associated with a known user.
    - error::String : Error indicates that the token couldn't be checked
    - user::IoK8sApiAuthenticationV1beta1UserInfo : User is the UserInfo associated with the provided token.
"""
mutable struct IoK8sApiAuthenticationV1beta1TokenReviewStatus <: SwaggerModel
    audiences::Any # spec type: Union{ Nothing, Vector{String} } # spec name: audiences
    authenticated::Any # spec type: Union{ Nothing, Bool } # spec name: authenticated
    error::Any # spec type: Union{ Nothing, String } # spec name: error
    user::Any # spec type: Union{ Nothing, IoK8sApiAuthenticationV1beta1UserInfo } # spec name: user

    function IoK8sApiAuthenticationV1beta1TokenReviewStatus(;audiences=nothing, authenticated=nothing, error=nothing, user=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1beta1TokenReviewStatus, Symbol("audiences"), audiences)
        setfield!(o, Symbol("audiences"), audiences)
        validate_property(IoK8sApiAuthenticationV1beta1TokenReviewStatus, Symbol("authenticated"), authenticated)
        setfield!(o, Symbol("authenticated"), authenticated)
        validate_property(IoK8sApiAuthenticationV1beta1TokenReviewStatus, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        validate_property(IoK8sApiAuthenticationV1beta1TokenReviewStatus, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type IoK8sApiAuthenticationV1beta1TokenReviewStatus

const _property_map_IoK8sApiAuthenticationV1beta1TokenReviewStatus = Dict{Symbol,Symbol}(Symbol("audiences")=>Symbol("audiences"), Symbol("authenticated")=>Symbol("authenticated"), Symbol("error")=>Symbol("error"), Symbol("user")=>Symbol("user"))
const _property_types_IoK8sApiAuthenticationV1beta1TokenReviewStatus = Dict{Symbol,String}(Symbol("audiences")=>"Vector{String}", Symbol("authenticated")=>"Bool", Symbol("error")=>"String", Symbol("user")=>"IoK8sApiAuthenticationV1beta1UserInfo")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1beta1TokenReviewStatus }) = collect(keys(_property_map_IoK8sApiAuthenticationV1beta1TokenReviewStatus))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1beta1TokenReviewStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthenticationV1beta1TokenReviewStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1beta1TokenReviewStatus }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1beta1TokenReviewStatus[property_name]

function check_required(o::IoK8sApiAuthenticationV1beta1TokenReviewStatus)
    true
end

function validate_property(::Type{ IoK8sApiAuthenticationV1beta1TokenReviewStatus }, name::Symbol, val)
end


text:
"""
    MCTerm(i, j, coeff, operator, source_orbital, integrals=[])

Represents one term in the multi-configurational expansion. `i` and
`j` are indices in the mixing-coefficient vector c (which is subject
to optimization, and thus has to be referred to), `coeff` is an
additional coefficient, and `integrals` is a list of indices into the
vector of common integrals, the values of which should be multiplied
to form the overall coefficient.
"""
struct MCTerm{T,QO,O}
    i::Int
    j::Int
    coeff::T
    operator::QO
    source_orbital::O
    integrals::Vector{Int}
end

"""
    OrbitalEquation(orbital, equation,
                    one_body, direct_terms, exchange_terms, source_terms)

Represents the integro-differential equation for `orbital`, expressed
as a linear combination of the different terms, with pointers to the
list of common integrals that is stored by the encompassing
[`MCEquationSystem`](@ref) object.
"""
struct OrbitalEquation{O,Equation}
    orbital::O
    equation::Equation
    terms::Vector{Pair{Int,Vector{MCTerm}}}
end

function Base.show(io::IO, oeq::OrbitalEquation)
    write(io, "OrbitalEquation($(oeq.orbital)): ")
    show(io, oeq.equation)
    write(io, "\n")
end

"""
    MCEquationSystem(equations, integrals)

Represents a coupled system of integro-differential `equations`,
resulting from the variation of a multi-configurational
[`EnergyExpression`](@ref), with respect to all constituent
orbitals. All `integrals` that are in common between the `equations`
need only be computed once per iteration, for efficiency.
"""
struct MCEquationSystem
    equations
    integrals
end

"""
    pushifmissing!(vector, element)

Push `element` to the end of `vector`, if not already present. Returns
the index of `element` in `vector`.
"""
function pushifmissing!(v::Vector, element)
    i = findfirst(isequal(element), v)
    isnothing(i) || return i
    push!(v, element)
    length(v)
end

"""
    orbital_equation(E::EnergyExpression, orbital, integrals::Vector)

Generate the [`OrbitalEquation`](@ref) governing `orbital` by varying
the [`EnergyExpression`](@ref) `E`, and storing common expressions in
`integrals`.
"""
function orbital_equation(E::EM, orbital::O, integrals::Vector) where {EM<:EnergyExpression,O}
    equation = diff(E, orbital)

    # Complementary orbital
    comp_orbital = orbital isa Conjugate ? orbital.orbital : Conjugate(orbital)

    terms = Dict{Int,Vector{MCTerm}}()

    # Invert matrix coordinate -> equation mapping, i.e. gather all
    # coordinates for which a specific NBodyEquation appears.

    for (i,j,eq) in zip(findnz(equation)...)
        for subeq = eq.equations
            operator = subeq.operator
            source_orbital = subeq.orbital
            coeff = MCTerm(i,j,subeq.factor.coeff,
                           operator, source_orbital,
                           map(factor -> pushifmissing!(integrals, factor),
                               subeq.factor.factors) |> Vector{Int})

            # If the operator isa ContractedOperator, an integral has
            # to be performed, which can potentially be reused among
            # common terms. This is only possible, however, if the
            # orbital under consideration is not inside the integrand;
            # if it is, we are dealing with an integral operator,
            # which has to be reevaluated each time it is applied to
            # an orbital.
            integral = (operator isa ContractedOperator && comp_orbital ∉ operator) ?
                pushifmissing!(integrals, operator) : 0
            terms[integral] = push!(get(terms, integral, MCTerm[]), coeff)
        end
    end

    OrbitalEquation(comp_orbital, equation, collect(pairs(terms)))
end

"""
    diff(energy_expression, orbitals)

Derive the integro-differential equations for all `orbitals`, from
`energy_expression`. Returns a [`MCEquationSystem`](@ref), that
gathers information on which integrals are common to all equations,
for efficient equation solving.
"""
function Base.diff(E::EM, orbitals::VO; verbosity=0) where {EM<:EnergyExpression, O,VO<:AbstractVector{O}}
    # Vector of common integrals
    integrals = []

    norb = length(orbitals)
    p = if verbosity > 0
        @info "Deriving equations for $(norb) orbitals"
        Progress(norb)
    end
    equations = map(orbitals) do orbital
        eq = orbital_equation(E, orbital, integrals)
        isnothing(p) || ProgressMeter.next!(p)
        eq
    end

    MCEquationSystem(equations, integrals)
end


text:
using Test

using CUDA

md = CuModuleFile(joinpath(@__DIR__, "vadd.ptx"))
vadd = CuFunction(md, "kernel_vadd")

dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)
c = similar(a)

d_a = CuArray(a)
d_b = CuArray(b)
d_c = CuArray(c)

len = prod(dims)
cudacall(vadd, Tuple{CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}}, d_a, d_b, d_c; threads=len)

@test a+b ≈ Array(d_c)


text:
using MCIntegrals
const P = MCIntegrals
using Test
using StaticArrays
using LinearAlgebra
using Random
using Setfield
using Cuba: vegas
using HCubature

function isconsistent(truth, est; nstd=6, kw_approx...)
    val = est.value
    Δ = nstd * est.std
    if ≈(val, truth; kw_approx...)
        true
    else
        truth - Δ <= val <= truth + Δ
    end
end

ALGS = [Vegas(), MCVanilla()]

@testset "RNG Reproducibility" begin
    for alg in ALGS
        @set! alg.rng = MersenneTwister(1)
        res1 = ∫(cos, (0,1), alg)
        @set! alg.rng = MersenneTwister(1)
        res2 = ∫(cos, (0,1), alg)
        @test res1 === res2
        res3 = ∫(cos, (0,1), alg)
        @test !(res3 === res2)
    end
end

@testset "Domain" begin
    dom = Domain((1f0, 2f0))
    @test P.pointtype(typeof(dom)) === SVector{1, Float32}
    @test P.volume(dom) === 1f0
end

@testset "exotic types $(typeof(alg))" for alg in [
        MCVanilla(10),
        Vegas(10),
        CubaAlg(vegas),
       ]
    est = ∫(identity, (0f0, 1f0), alg)
    @test typeof(est.value) === Float32
    @test typeof(est.std) === Float32

    est = ∫(_ -> true, (0f0, 1f0), alg)
    @test typeof(est.value) === Float64
    @test typeof(est.std) === Float64

    est = ∫(identity, (0, 1), alg)
    @test typeof(est.value) === Float64
    @test typeof(est.std) === Float64

    est = ∫(identity, (big"0", big"1"), alg)
    @test typeof(est.value) === BigFloat
    @test typeof(est.std) === BigFloat

    v = @SVector[1.,2.]
    est = ∫(_ -> v, (0,1), alg)
    @test typeof(est.value) === typeof(v)
    @test typeof(est.std)   === typeof(est.value)
end

@testset "std bernoulli" begin
    p = 0.1 + 0.8*rand()
    f = x -> 0 <= x <= p

    neval = 1000
    true_value = p
    true_std = sqrt(p*(1-p)/neval)

    est = ∫(f, (0,1), MCVanilla(neval))
    @test est.std ≈ true_std rtol=0.2
    @test isconsistent(true_value, est)

    est = ∫(f, (0,1), Vegas(neval))
    @test est.std < 0.2*true_std
    @test isconsistent(true_value, est)
end

@testset "constant $(typeof(alg))" for alg in [
    Vegas(10^1),
    MCVanilla(10^1),
    CubaAlg(vegas),
    HCubatureAlg(),
    ]
    for dim in 1:4
        for _ in 1:10
            lower = @SVector randn(dim)
            upper = lower + @SVector rand(dim)
            val = randn()
            f = _ -> val
            dom = Domain(lower, upper)
            est = ∫(f, dom, alg)
            vol = prod(upper .- lower)
            truth = val * vol
            @test isconsistent(truth, est)
            @test est.std < 1e-3
        end
    end
end

@testset "∫x^2" begin
    est_vegas = ∫(x->x^2, (0,1), Vegas(10000))
    @test isconsistent(1/3, est_vegas)
    @test est_vegas.std < 1e-3

    est_vanilla = ∫(x->x^2, (0,1), MCVanilla(1000))
    @test isconsistent(1/3, est_vanilla)
    @test est_vanilla.std < 1e-2
end

@testset "∫(f, (a,b), $(typeof(alg)))" for alg in [
        MCVanilla(10^4),
        Vegas(10^4),
       ]

    funs = [
        (f = x -> x^2, F = x -> x^3/3),
        (f = cos, F = sin),
        let
            a = randn()
            (f = _ -> a, F = x -> a*x)
        end
    ]

    for (f,F) in funs
        for _ in 1:10
            a = 10randn()
            b = a + 10rand()
            est = ∫(f, (a,b), alg)
            truth = F(b) - F(a)
            @test isconsistent(truth, est)
        end
    end
end

@testset "Ball $(typeof(alg))" for alg in [
    MCVanilla(10^4),
    Vegas(10^4, rng=MersenneTwister(1))]

    f(x) = norm(x) < 1 ? 1. : 0.
    est = ∫(f, (-2,2), alg)
    @test isconsistent(2,est)

    est = ∫(f, ((-1,1),(-1,1)), alg)
    @test isconsistent(π, est)

    est = ∫(f, ((-1,1),(-1,1),(-1,1)), alg)
    @test isconsistent(4π/3, est)
end

@testset "Fubini $(typeof(alg))" for alg in [
    MCVanilla(10^4),
    Vegas(10^4)]
    for _ in 1:10
        a1,a2 = randn(2)
        b1 = a1 + 10rand()
        b2 = a2 + 10rand()
        dom1 = Domain((a1, b1))
        dom2 = Domain((a2, b2))
        dom = Domain(((a1, b1), (a2, b2)))
        f1 = x -> cos(x[1])
        F1 = x -> sin(x[1])
        f2 = x -> x[1]^2
        F2 = x -> x[1]^3/3
        f  = v -> f1(v[1]) * f2(v[2])
        est1 = ∫(f1, dom1, alg)
        est2 = ∫(f2, dom2, alg)
        est  = ∫(f , dom , alg)
        @test isconsistent(est1.value*est2.value, est)
    end
end

@testset "Simple Vegas Example" begin
    a = -0.237267860990952
    b = -0.19099487884978464
    val = 0.665843778891144
    truth = (b - a) * val
    f = _ -> val
    Ω = Domain((a,b))

    alg = Vegas(5)
    est = ∫(f, Ω, alg)
    @test est.value ≈ truth
    @test 0 <= est.std <= sqrt(eps(Float64))

    iq = P.initvr(f, Ω, alg)
    est = ∫(f, iq, alg)
    @test est.value ≈ truth
    @test 0 <= est.std <= sqrt(eps(Float64))

    iq2 = P.tune(f, Ω, alg)
    est = ∫(f, iq2, alg)
    @test est.value ≈ truth
    @test 0 <= est.std <= sqrt(eps(Float64))
    walls = iq2.boundaries[1]
    @test walls ≈ range(a,stop=b, length=length(walls))
end


text:
using DistributedArrays
@everywhere importall DistributedArrays

@everywhere function get_data_plus_gc(domain, nguard, ngrid)
    if myid() in procs(domain)
        li = localindexes(domain)
        lp = localpart(domain)

        s = size(lp)
        data_plus_gc = zeros(s[1]+2*nguard, s[2]+2*nguard)
        for j in 1:s[2]
            for i in 1:s[1]
                data_plus_gc[i+nguard, j+nguard] = lp[i,j]
            end
        end

        xstart = li[1][1]
        xend   = li[1][end]
        ystart = li[2][1]
        yend   = li[2][end]

        for g in 1:nguard
            xsg = (xstart-1-g + ngrid) % ngrid + 1
            xeg = (xend-1+g) % ngrid + 1

            ysg = (ystart-1-g + ngrid) % ngrid + 1
            yeg = (yend-1+g) % ngrid + 1

            for i in 1+nguard:s[1]+nguard
                @inbounds data_plus_gc[i, nguard+1-g] = domain[i-nguard+xstart-1, ysg]
                @inbounds data_plus_gc[i, s[2]+nguard+g] = domain[i-nguard+xstart-1, yeg]
            end
            for j in 1+nguard:s[2]+nguard
                @inbounds data_plus_gc[nguard+1-g, j] = domain[xsg, j-nguard+ystart-1]
                @inbounds data_plus_gc[s[1]+nguard+g, j] = domain[xeg, j-nguard+ystart-1]
            end
        end
    end
    return data_plus_gc
end

@everywhere function advect_data(dens, nguard, ngrid, velx, vely, dx, dy, dt)
    locdens = get_data_plus_gc(dens, nguard, ngrid)

    s = size(locdens)
    nx = s[1] - 2*nguard
    ny = s[2] - 2*nguard

    gradx = zeros(locdens)
    grady = zeros(locdens)

    for j in 1+nguard:ny+nguard
        for i in 1+nguard:nx+nguard
            if velx > 0
                @inbounds gradx[i,j] = (3locdens[i,j] - 4locdens[i-1,j] + locdens[i-2,j])/(2dx)
            else
                @inbounds gradx[i-2,j] = (-locdens[i,j] + 4locdens[i-1,j] - 3locdens[i-2,j])/(2dx)
            end

            if vely > 0
                @inbounds grady[i,j] = (3locdens[i,j] - 4locdens[i,j-1] + locdens[i,j-1])/(2dy)
            else
                @inbounds grady[i,j-2] = (-locdens[i,j] + 4locdens[i,j-1] - 3locdens[i,j-2])/(2dy)
            end
        end
    end

    for j in 1+nguard:ny+nguard
        for i in 1+nguard:nx+nguard
            @inbounds localpart(dens)[i-nguard, j-nguard] -= dt*(velx*gradx[i,j] + vely*grady[i,j])
        end
    end
end

function initial_conditions(ngrid, dx, dy, posx, posy, sigma)
    dens = zeros(ngrid, ngrid)
    for j in 1:ngrid
        y = (j-1)*dy
        for i in 1:ngrid
            x = (i-1)*dx
            @inbounds dens[i,j] = exp(-((x-posx)^2 + (y-posy)^2)/(sigma^2))
        end
    end
    distribute(dens)
end

function output_csv(filename, dens)
    open(filename, "w") do f
    denslocal = convert(Array, dens)
        for j in 1:ngrid
            for i in 1:ngrid
                @inbounds d = dens[i,j]
                write(f,"$i, $j, $d\n")
            end
         end
    end
end

function timestep(dens, nguard, ngrid, velx, vely, dx, dy, dt)
    ps = procs(dens)
    refs = [(@spawnat p advect_data(dens, nguard, ngrid, velx, vely, dx, dy, dt)) for p in ps]
    pmap(fetch, refs)
end

ngrid = 1000
ntimesteps = 500
output = false

nargs = length(ARGS)
if nargs > 0
    ngrid = parse(Int32, ARGS[1])
end
if nargs > 1
    ntimesteps = parse(Int32, ARGS[2])
end
if nargs > 2
    output = parse(Bool, ARGS[2])
end

posx = 0.3
posy = 0.3
sigma = 0.15
dx = 1./ngrid
dy = 1./ngrid
dens = initial_conditions(ngrid, dx, dy, posx, posy, sigma)

if output
    output_csv("init_jl.csv", dens)
end

velx = 1.0
vely = 1.0
velmag = sqrt(velx*velx + vely*vely)
cfl = 0.125
dt = cfl*min(dx, dy)/velmag

for t in 1:ntimesteps
    timestep(dens, 2, ngrid, velx, vely, dx, dy, dt)
end

if output
    output_csv("final_jl.csv", dens)
end


text:
######################################
using DataFrames, CSV, CodecZlib, Mmap

#= Polytope Data from text files: 29×3 DataFrame
│ Row │ Database │ Polytopes │ cumulative_sum │
│     │ Int8     │ Int32     │ Int64          │
├─────┼──────────┼───────────┼────────────────┤
│ 1   │ 5        │ 1561      │ 1561           │
│ 2   │ 6        │ 24189     │ 25750          │
│ 3   │ 7        │ 177446    │ 203196         │
│ 4   │ 8        │ 834638    │ 1037834        │
│ 5   │ 9        │ 2867955   │ 3905789        │
│ 6   │ 10       │ 7725801   │ 11631590       │
│ 7   │ 11       │ 16608387  │ 28239977       │
│ 8   │ 12       │ 29270253  │ 57510230       │
│ 9   │ 13       │ 43458000  │ 100968230      │
│ 10  │ 14       │ 56060584  │ 157028814      │
│ 11  │ 15       │ 64085869  │ 221114683      │
│ 12  │ 16       │ 65615931  │ 286730614      │
│ 13  │ 17       │ 59972682  │ 346703296      │
│ 14  │ 18       │ 48703033  │ 395406329      │
│ 15  │ 19       │ 34847821  │ 430254150      │
│ 16  │ 20       │ 21913680  │ 452167830      │
│ 17  │ 21       │ 12070919  │ 464238749      │
│ 18  │ 22       │ 5826221   │ 470064970      │
│ 19  │ 23       │ 2450720   │ 472515690      │
│ 20  │ 24       │ 898929    │ 473414619      │
│ 21  │ 25       │ 284696    │ 473699315      │
│ 22  │ 26       │ 78468     │ 473777783      │
│ 23  │ 27       │ 18417     │ 473796200      │
│ 24  │ 28       │ 3781      │ 473799981      │
│ 25  │ 29       │ 647       │ 473800628      │
│ 26  │ 30       │ 114       │ 473800742      │
│ 27  │ 31       │ 23        │ 473800765      │
│ 28  │ 32       │ 8         │ 473800773      │
│ 29  │ 33       │ 2         │ 473800775      │   =#

const DIM = 4

# Construct polytope from the polytope coordinates in verts and append to DataFrame
# of all polytopes (Δ²). Note that the length of the polytope and the classification
# number are embedded in the DataFrame.
function MakePolytope!(Δ²::DataFrame, verts, ontology)
    Δ = hcat(verts...)
    size(Δ,1) < size(Δ,2) ? Δ = Δ' : nothing;
    len = convert(Int32,size(Δ, 1))
    delim = zeros(len)
    delim[1] += ontology.v
    delim[2] += ontology.p
    delim[3] += ontology.v_dual
    delim[4] += ontology.p_dual
    Δ = hcat(delim, Δ)
    header = zeros(1 + 4)
    header[1] = len
    header[2] = ontology.h11
    header[3] = ontology.h21
    header[4] = ontology.euler
    Δ = vcat(header', Δ)
    Δ = convert(DataFrame, Δ)
    println(Δ)
    run(`quit`)
    append!(Δ², Δ)
end

# Write Δ² to a .gzip file and compress. Each sub-database
# has a simple integer naming scheme (sub_db_id).
function ExportSubDB(Δ²::DataFrame, sub_db_id::Integer)
    write_path = "/media/share/Dev/CalabiYau/data/polytopes_db_4d/$(sub_db_id).gz"
    open(write_path, "w") do io
        stream = GzipCompressorStream(io)
        CSV.write(stream, Δ²)
        close(stream)
        stream = nothing
        GC.gc()
    end #close .gzf
end

# Grab feature data
function GetOntology(ln::String)
    info =  split(strip(ln))
    Δ⁰_pnts = parse(Int,split(info[3], ":")[2])
    Δ⁰_vrts = parse(Int,info[4])
    Δ_pnts = parse(Int,split(info[5], ":")[2])
    Δ_vrts = parse(Int,info[6])
    H = split(split(info[7], ":")[2], ",")
    h11 = parse(Int,H[1])
    h21 = parse(Int,H[2])
    euler = parse(Int,split(split(info[end], "[")[2],"]")[1])
    return (p_dual = Δ⁰_pnts, v_dual = Δ⁰_vrts,p = Δ_pnts, v =  Δ_vrts,h11= h11,h21= h21,euler= euler)
end

function DoTheThing()#(a::Integer, b::Integer, polytopeCount::Integer, info::DataFrame)
    Δ² = DataFrame()::DataFrame
    ontology = ()
    polytopeCount = 0::Integer
    verts = []
    sub_db_id = 1::Integer
    polyPerDB = 0::Integer
    info = DataFrame(Database = Int8[], Polytopes = Int32[])

    for id in 5:33
        # Initialize path to current database file
        db = "$id"
        if length(db) == 1  db = "0"*db; end
        data_path = "/media/share/Dev/CalabiYau/data/polytopes_db_4d/src/v$db"
        # Open current database file
        open(data_path) do file
            firstLine = true
            polyPerDB = 0
            for ln in eachline(file)
                # Are we on the first line?
                if firstLine
                    ontology = GetOntology(ln)
                    firstLine = false
                    continue
                end
                # Are we on the last line?
                if occursin("#", ln)
                    MakePolytope!(Δ², verts, ontology)
                    verts = nothing
                    polytopeCount += 1
                    polyPerDB += 1
                    verts = []
                    break
                end
                # Have we finished reading a polytope?
                if occursin("[", ln)
                    #= Construct polytope from the polytope coordinates in verts, append
                    to DataFrame of all polytopes (Δ²), and reset relevant variables       =#
                    MakePolytope!(Δ², verts, ontology)
                    verts = nothing
                    polytopeCount += 1
                    polyPerDB += 1
                    verts = []
                    #= If Δ² has reached 1,000,000 polytopes, then write Δ² to file. After
                    compression, iterate the sub-database identifer and parse the next
                    classification number (ontology) from the current line.             =#
                    if polytopeCount == 1000000
                        ExportSubDB(Δ², sub_db_id)
                        Δ² = nothing
                        sub_db_id += 1
                        polytopeCount = 0
                        Δ² = DataFrame()
                    end #end polytope count check
                    ontology = GetOntology(ln)
                else # if reading a coordinate line
                    row = [x for x in parse.(Int16, split(strip(ln)))]
                    push!(verts, row)
                end #end current line
            end #end eachline
            current_db_info = DataFrame(Database = id, Polytopes = polyPerDB)
            append!(info, current_db_info)
            run(`clear`)
            println(info)
        end #close file
    end
    ## Export final (partial) database
    ExportSubDB(Δ², sub_db_id)

    info.cumulative_sum = cumsum(info.Polytopes)
    run(`clear`)
    println(info)
end

DoTheThing()


text:
"""
Implements the reaction network detailed in https://www.pnas.org/content/pnas/suppl/2003/05/27/1132112100.DC1/2112SuppText.pdf

Accompanying article is here https://www.pnas.org/content/100/12/7051
"""

function CircadianOscillator(input=nothing)
println("default input is t -> 1.5. If CircadianOscillator() is called with no arguments, the default is used")
@parameters t
D = Differential(t)
@variables  M_P(t) M_C(t) M_B(t) P_C(t) C_C(t) P_CP(t) C_CP(t) PC_C(t) PC_N(t) PC_CP(t) PC_NP(t) B_C(t) B_CP(t) B_N(t) B_NP(t) I_N(t)

@parameters k_1 k_2 k_3 k_4 k_5 k_6 k_7 k_8 K_AP K_AC K_IB k_dmb k_dmc  k_dmp k_dn  k_dnc K_d K_dp K_p  K_mB K_mC K_mP k_sB k_sC k_sP m n V_1B V_1C V_1P V_1PC V_2B V_2C V_2P V_2PC V_3B V_3PC V_4B V_4PC V_phos v_dBC v_dBN v_dCC v_dIN v_dPC v_dPCC v_dPCN v_mB v_mC v_mP v_sB v_sC #v_sP

#also have v_phos = 0.4 in the paper

#v_sP is the transcription rate of PER mRNA. In the paper, the effect of light is modelled through changes in this parameter.

v_sP = input
if v_sP === nothing
    v_sP = 1.5
end

eqs =
[ D(M_P) ~ v_sP*(B_N^n/(K_AP^n + B_N^n)) - v_mP*(M_P/(K_mP + M_P)) - k_dmp*M_P,
  D(M_C) ~ v_sC*B_N^n/(K_AC^n + B_N^n) - v_mC*(M_C/(K_mC + M_C)) - k_dmc*M_C,
  D(M_B) ~ v_sB*K_IB^m/(K_IB^m + B_N^m ) - v_mB*(M_B/(K_mB + M_B)) - k_dmb*M_B,
  D(P_C) ~ k_sP*M_P - V_1P*P_C/(K_p + P_C) + V_2P*P_CP/(K_dp + P_CP) + k_4*PC_C - k_3*P_C*C_C - k_dn*P_C,
  D(C_C) ~ k_sC*M_C - V_1C*C_C/(K_p + C_C) + V_2C*C_CP/(K_dp + C_CP) + k_4*PC_C - k_3*P_C*C_C - k_dnc*C_C,
  D(P_CP) ~ V_1P*P_C/(K_p+P_C) - V_2P*P_CP/(K_dp + P_CP) - v_dPC*P_CP/(K_d + P_CP) - k_dn*P_CP,
  D(C_CP) ~ V_1C*C_C/(K_p +C_C) - V_2C*C_CP/(K_dp+C_CP) - v_dCC*C_CP/(K_d + C_CP) - k_dn*C_CP,
  D(PC_C) ~ -V_1PC*PC_C/(K_p + PC_C) + V_2PC*PC_CP/(K_dp + PC_CP) - k_4*PC_C + k_3*P_C*C_C + k_2*PC_N - k_1*PC_C - k_dn*PC_C,
  D(PC_N) ~ -V_3PC*PC_N/(K_p + PC_N) + V_4PC*PC_NP/(K_dp + PC_NP) - k_2*PC_N + k_1*PC_C - k_7*B_N*PC_N + k_8*I_N - k_dn*PC_N,
  D(PC_CP) ~ V_1PC*PC_C/(K_p + PC_C) - V_2PC*PC_CP/(K_dp + PC_CP) - v_dPCC*PC_CP/(K_d + PC_CP) - k_dn*PC_CP,
  D(PC_NP) ~ V_3PC*PC_N/(K_p + PC_N) - V_4PC*PC_NP/(K_dp + PC_NP) - v_dPCN*PC_NP/(K_d + PC_NP) - k_dn*PC_NP,
  D(B_C) ~ k_sB*M_B - V_1B*B_C/(K_p + B_C) + V_2B*B_CP/(K_dp + B_CP) - k_5*B_C + k_6*B_N - k_dn*B_C,
  D(B_CP) ~ V_1B*B_C/(K_p + B_C) - V_2B*B_CP/(K_dp + B_CP) - v_dBC*B_CP/(K_d + B_CP) - k_dn*B_CP,
  D(B_N) ~ -V_3B*B_N/(K_p + B_N) + V_4B*B_NP/(K_dp + B_NP) + k_5*B_C - k_6*B_N - k_7*B_N*PC_N + k_8*I_N - k_dn*B_N,
  D(B_NP) ~ V_3B*B_N/(K_p + B_N) - V_4B*B_NP/(K_dp + B_NP) - v_dBN*B_NP/(K_d + B_NP) - k_dn*B_NP,
  D(I_N) ~ -k_8*I_N + k_7*B_N*PC_N - v_dIN*I_N/(K_d + I_N) - k_dn*I_N
  ]

tspan = (0.,100.)

ic = [M_P=> 2.391615385357426 ,M_C=> 1.7547090952244027 ,M_B=> 8.561088205594936 ,P_C=> 0.9030978254420003 ,C_C=> 3.5790996574043907 ,P_CP=> 0.11556396765994004 ,C_CP=> 0.6301918743084862 ,PC_C=> 0.8525514643541497 ,PC_N=> 0.09761528221027084 ,PC_CP=> 0.1755706805316031 ,PC_NP=> 0.06289128872787488 ,B_C=> 2.2764826266300875 ,B_CP=> 0.9063725609755553 ,B_N=> 1.763837030511116 ,B_NP=> 0.32028080421972577 ,I_N=> 0.024197413185517866 ]

ps = [k_1=>0.4, k_2=>0.2, k_3=>0.4, k_4=>0.2, k_5=>0.4, k_6=>0.2, k_7=>0.5, k_8=>0.1, K_AP=>0.7, K_AC=>0.6, K_IB=>2.2, k_dmb=>0.01, k_dmc => 0.01, k_dmp=>0.01, k_dn => 0.01, k_dnc=>0.12, K_d=>0.3, K_dp=>0.1, K_p => 0.1, K_mB=>0.4, K_mC=>0.4, K_mP=>0.31, k_sB=>0.12, k_sC=>1.6, k_sP=>0.6, m=>2, n=>4, V_1B=>0.5, V_1C=>0.6, V_1P=>0.4, V_1PC=>0.4, V_2B=>0.1, V_2C=>0.1, V_2P=>0.3, V_2PC=>0.1, V_3B=>0.5, V_3PC=>0.4, V_4B=>0.2, V_4PC=>0.1, v_dBC=>0.5, v_dBN=>0.6, v_dCC=>0.7, v_dIN=>0.8, v_dPC=>0.7, v_dPCC=>0.7, v_dPCN=>0.7, v_mB=>0.8, v_mC=>1.0, v_mP=>1.1, v_sB=>1.0, v_sC=>1.1] #, v_sP=>1.5]

od = ODESystem(eqs, t, first.(ic), first.(ps))
# of = ODEFunction(od)
tspan = (0.,100.)
return od, ic, tspan, ps
end


text:
@testset "NumField/NfAbs/NfAbs" begin
  cyclo_expl = function(n, m)
    Fn, zn = CyclotomicField(n)
    Fnm, znm = CyclotomicField(n*m)
    x = zn
    x_up = Hecke.force_coerce_cyclo(Fnm, x)
    x_down = Hecke.force_coerce_cyclo(Fn, x_up)
    return (x, x_up, x_down)
  end

  res = cyclo_expl(3, 4)
  @test (res[1]^3, res[2]^3) == (1, 1)

  res = cyclo_expl(3, 2)
  z6 = gen(parent(res[2]))
  # Test that z3 is mapped to z6^2
  @test z6^2 == res[2]
end

@testset "Splitting Field" begin

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^3-2
  K = splitting_field(f)
  @test typeof(K) == AnticNumberField
  K1 = number_field([x^3-2, x^2+x+1])[1]
  K1abs = simple_extension(K1)[1]
  @test isisomorphic(K, K1abs)[1]
  K, R = splitting_field(f, do_roots = true)
  for r in R
    @test iszero(f(r))
  end

  K, a = number_field(x^2+1)
  Kt, t = PolynomialRing(K, "t")
  g = t^4-2
  L = splitting_field(g)
  @test typeof(L) == Hecke.NfRel{nf_elem}
  @test absolute_degree(L) == 8
end

@testset "simplify" begin

  Qx, x = PolynomialRing(FlintQQ, "x")
  K = number_field(x^5 - x^4 - x^3 - 220*x^2 - 360*x - 200)[1]
  @test typeof(K) == AnticNumberField
  L = simplify(K, canonical = true)[1]
  @test L.pol == x^5 - x^4 + x^3 - 19*x^2 + 35*x + 77
end

@testset "simplify-Gunter" begin

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^18 - x^16 - 6*x^14 - 4*x^12 - 4*x^10 + 2*x^8 + 6*x^6 - 4*x^4 + 3*x^2 - 1
  g = x^18 - 3*x^16 + 4*x^14 - 6*x^12 - 2*x^10 + 4*x^8 + 4*x^6 + 6*x^4 + x^2 - 1
  h = x^18 + x^16 - x^14 - 8*x^12 - 3*x^8 + 27*x^6 - 25*x^4 + 8*x^2 - 1
  @test simplify(number_field(f)[1], canonical = true)[1].pol == g
  @test simplify(number_field(g)[1], canonical = true)[1].pol == g
  @test simplify(number_field(h)[1], canonical = true)[1].pol == g
end


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ 8bf5ade3-3bd9-425a-a9e1-3c1cf749e1c1
using Pkg, DrWatson

# ╔═╡ 8ea8eabf-b063-40d1-b83c-4ed24cea4ef6
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ 0fbfaf77-6389-4ac1-8c9e-ae80d339336a
md" ### j-02-04.jl"

# ╔═╡ d0fbe0f6-d839-48f7-a1e0-14c785e4625d
begin
	a = Float64[0 0 16; 0 4 5; 8 -4 22]
	b = Float64[4, 2, 5]
end

# ╔═╡ d89b625f-aa4c-4e39-af34-b8a1dfcdaa3b
md" ##### Cholesky factorization using banded storage"

# ╔═╡ 3b0869a9-ae43-42d8-b8e6-cb21107aeef3
md" ##### Banded Coefficient Matrix:"

# ╔═╡ ca47bbfe-a36f-4214-98a3-7b72c1041980
a

# ╔═╡ 7b594403-1a12-46eb-ba83-d0c8dfdef682
md" ##### Right Hand side:"

# ╔═╡ 00fa6a6f-297f-4fd7-8b2c-8c0f7a1571d5
b

# ╔═╡ 44d0ef6b-9a6a-4e40-8f96-fbc4f33e468e
cholin!(a)

# ╔═╡ 59b640c3-e769-42d2-aa48-a6a00a0bd0cf
md" ##### L in Band Form:"

# ╔═╡ 5f3f4836-8270-4ba2-9eab-1b83816b56f4
a

# ╔═╡ f072e05c-9086-4f9c-bd2f-fe2bd9d262ec
chobac!(a, b)

# ╔═╡ f9b02427-0328-42e5-831b-23c77194b56f
md" ##### Solution Vector:"

# ╔═╡ 4c504c4e-7804-44af-ba2e-8f4cb84956c7
b

# ╔═╡ 15bac539-5e4f-4baa-b001-799feaf0f9a5
md" ### End of j-02-04.jl"

# ╔═╡ Cell order:
# ╟─0fbfaf77-6389-4ac1-8c9e-ae80d339336a
# ╠═8bf5ade3-3bd9-425a-a9e1-3c1cf749e1c1
# ╠═8ea8eabf-b063-40d1-b83c-4ed24cea4ef6
# ╠═d0fbe0f6-d839-48f7-a1e0-14c785e4625d
# ╠═d89b625f-aa4c-4e39-af34-b8a1dfcdaa3b
# ╠═3b0869a9-ae43-42d8-b8e6-cb21107aeef3
# ╠═ca47bbfe-a36f-4214-98a3-7b72c1041980
# ╠═7b594403-1a12-46eb-ba83-d0c8dfdef682
# ╠═00fa6a6f-297f-4fd7-8b2c-8c0f7a1571d5
# ╠═44d0ef6b-9a6a-4e40-8f96-fbc4f33e468e
# ╠═59b640c3-e769-42d2-aa48-a6a00a0bd0cf
# ╠═5f3f4836-8270-4ba2-9eab-1b83816b56f4
# ╠═f072e05c-9086-4f9c-bd2f-fe2bd9d262ec
# ╠═f9b02427-0328-42e5-831b-23c77194b56f
# ╠═4c504c4e-7804-44af-ba2e-8f4cb84956c7
# ╟─15bac539-5e4f-4baa-b001-799feaf0f9a5


text:
"Export a [`ClimateModel`](@ref) instance to a `.json` file.

This `.json` file can be shared with others, you can import it using "
function export_state(filename::AbstractString, m::ClimateModel)
    open(filename, "w") do io
        JSON2.write(io, m)
    end
end

function import_state(filename::String)::ClimateModel
    open(filename, "r") do io
        return JSON2.read(io, ClimateModel)
    end
end

function export_parameters(filename::AbstractString, params::ClimateModelParameters)
    open(filename, "w") do io
        JSON2.write(io, params)
    end
end

function import_parameters(filename::AbstractString)::ClimateModelParameters
    open(filename, "r") do io
        return JSON2.read(io, ClimateModelParameters)
    end
end

"""The [`ClimateModelParameters`](@ref) included with this package.

Currently `included_configurations["default"]` is the only included set."""
const included_configurations = let
    # find the config dir relative to this .jl file
    config_dir = joinpath(@__DIR__, "..", "..", "configurations")
    config_files = [file for file in readdir(config_dir) if occursin(".json", file)]
    loaded = [import_parameters(joinpath(config_dir, file)) for file in config_files]
    Dict(p.name => p for p in loaded)
end


text:
__precompile__()
module Weave
import Highlights
using Compat
using Requires

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" Base.include(Main, "plots.jl")
    @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" Base.include(Main, "gadfly.jl")
end

"""
`list_out_formats()`

List supported output formats
"""
function list_out_formats()
  for format = keys(formats)
      println(string(format,": ",  formats[format].description))
  end
end

"""
`tangle(source ; out_path=:doc, informat="noweb")`

Tangle source code from input document to .jl file.

* `informat`: `"noweb"` of `"markdown"`
* `out_path`: Path where the output is generated. Can be: `:doc`: Path of the source document, `:pwd`: Julia working directory,  `"somepath"`, directory name as a string e.g `"/home/mpastell/weaveout"`
or filename as string e.g. ~/outpath/outfile.jl.
"""
function tangle(source ; out_path=:doc, informat=:auto)
    doc = read_doc(source, informat)
    doc.cwd = get_cwd(doc, out_path)

    outname = get_outname(out_path, doc, ext = "jl")

    open(outname, "w") do io
    for chunk in doc.chunks
      if typeof(chunk) == CodeChunk
          options = merge(rcParams[:chunk_defaults], chunk.options)
          if options[:tangle]
            write(io, chunk.content*"\n")
          end
      end
    end
  end
  doc.cwd == pwd()  && (outname = basename(outname))
  @info("Writing to file $outname")
end

"""
    weave(source ; doctype = :auto,
        informat=:auto, out_path=:doc, args = Dict(),
        mod::Union{Module, Symbol} = Main,
        fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache=:off,
        template = nothing, highlight_theme = nothing, css = nothing,
        pandoc_options = "",
        latex_cmd = "xelatex")

Weave an input document to output file.

* `doctype`: :auto = set based on file extension or specify one of the supported formats.
  See `list_out_formats()`
* `informat`: :auto = set based on file extension or set to  `"noweb"`, `"markdown"` or  `script`
* `out_path`: Path where the output is generated. Can be: `:doc`: Path of the source document, `:pwd`:
   Julia working directory, `"somepath"`: output directory as a String e.g `"/home/mpastell/weaveout"` or filename as
   string e.g. ~/outpath/outfile.tex.
* `args`: dictionary of arguments to pass to document. Available as WEAVE_ARGS
* `mod`: Module where Weave `eval`s code. Defaults to `:sandbox`
   to create new sandbox module, you can also pass a module e.g. `Main`.
* `fig_path`: where figures will be generated, relative to out_path
* `fig_ext`: Extension for saved figures e.g. `".pdf"`, `".png"`. Default setting depends on `doctype`.
* `cache_path`: where of cached output will be saved.
* `cache`: controls caching of code: `:off` = no caching, `:all` = cache everything,
  `:user` = cache based on chunk options, `:refresh`, run all code chunks and save new cache.
* `throw_errors` if `false` errors are included in output document and the whole document is
    executed. if `true` errors are thrown when they occur.
* `template` : Template (file path) for md2html or md2tex formats.
* `highlight_theme` : Theme (Highlights.AbstractTheme) for used syntax highlighting
* `css` : CSS (file path) used for md2html format
* `pandoc_options` = String array of options to pass to pandoc for `pandoc2html` and
   `pandoc2pdf` formats e.g. ["--toc", "-N"]
* `latex_cmd` the command used to make pdf from .tex

**Note:** Run Weave from terminal and not using IJulia, Juno or ESS, they tend to mess with capturing output.
"""
function weave(source ; doctype = :auto,
        informat=:auto, out_path=:doc, args = Dict(),
        mod::Union{Module, Symbol} = :sandbox,
        fig_path = "figures", fig_ext = nothing,
        cache_path = "cache", cache=:off,
        throw_errors = false,
        template = nothing, highlight_theme = nothing, css = nothing,
        pandoc_options = String[]::Array{String},
        latex_cmd = "xelatex")

    doc = read_doc(source, informat)
    highlight_theme != nothing && (doc.highlight_theme = highlight_theme)
    #theme != nothing && (doc.theme = theme) #Reserved for themes
    css != nothing && (doc.css = css)
    template != nothing && (doc.template = template)

    try
      doc = run(doc, doctype = doctype,
              mod = mod,
              out_path=out_path, args = args,
              fig_path = fig_path, fig_ext = fig_ext, cache_path = cache_path, cache=cache,
              throw_errors = throw_errors)
      formatted = format(doc)

      outname = get_outname(out_path, doc)

      open(outname, "w") do io
          write(io, formatted)
      end

      #Special for that need external programs
      if doc.doctype == "pandoc2html"
          mdname = outname
          outname = get_outname(out_path, doc, ext = "html")
          pandoc2html(formatted, doc, outname, pandoc_options)
          rm(mdname)
      elseif doc.doctype == "pandoc2pdf"
          mdname = outname
          outname = get_outname(out_path, doc, ext = "pdf")
          pandoc2pdf(formatted, doc, outname, pandoc_options)
          rm(mdname)
      elseif doc.doctype == "md2pdf"
          success = run_latex(doc, outname, latex_cmd)
          success && rm(doc.fig_path, force = true, recursive = true)
          success || return
          outname = get_outname(out_path, doc, ext = "pdf")
      end

      doc.cwd == pwd() && (outname = basename(outname))
      @info("Report weaved to $outname")
    #    catch err
    #    @warn("Something went wrong during weaving")
    #    println(e)
    finally
        doctype == :auto && (doctype = detect_doctype(doc.source))
        if occursin("pandoc2pdf", doctype) && cache == :off
            rm(doc.fig_path, force = true, recursive = true)
        elseif occursin("2html", doctype)
            rm(doc.fig_path, force = true, recursive = true)
        end
    end
end

function weave(doc::AbstractString, doctype::AbstractString)
    weave(doc, doctype=doctype)
end

"""
  notebook(source::String, out_path=:pwd)

Convert Weave document `source` to Jupyter notebook and execute the code
using nbconvert. Requires IJulia. **Ignores** all chunk options

* `out_path`: Path where the output is generated. Can be: `:doc`: Path of the source document,
   `:pwd`: Julia working directory, `"somepath"`: Path as a
    String e.g `"/home/mpastell/weaveout"`
* nbconvert cell timeout in seconds. Defaults to -1 (no timeout)
"""
function notebook(source::String, out_path=:pwd, timeout=-1)
  doc = read_doc(source)
  converted = convert_doc(doc, NotebookOutput())
  doc.cwd = get_cwd(doc, out_path)
  outfile = get_outname(out_path, doc, ext="ipynb")

  open(outfile, "w") do f
    write(f, converted)
  end

  @info("Running nbconvert")
  eval(Meta.parse("using IJulia"))
  out = read(`$(IJulia.jupyter)-nbconvert --ExecutePreprocessor.timeout=$timeout --to notebook --execute $outfile --output $outfile`, String)
end

"""
    include_weave(doc, informat=:auto)

Include code from Weave document calling `include_string` on
all code from doc. Code is run in the path of the include document.
"""
function include_weave(source, informat=:auto)
  old_path = pwd()
  doc = read_doc(source, informat)
  cd(doc.path)
  try
    code = join([x.content for x in
      filter(x -> isa(x,Weave.CodeChunk), doc.chunks)], "\n")
    include_string(code)
  catch e
    cd(old_path)
    throw(e)
  end
end

#Hooks to run before and after chunks, this is form IJulia,
#but note that Weave hooks take the chunk as input
const preexecute_hooks = Function[]
push_preexecute_hook(f::Function) = push!(preexecute_hooks, f)
pop_preexecute_hook(f::Function) = splice!(preexecute_hooks, findfirst(preexecute_hooks, f))

const postexecute_hooks = Function[]
push_postexecute_hook(f::Function) = push!(postexecute_hooks, f)
pop_postexecute_hook(f::Function) = splice!(postexecute_hooks, findfirst(postexecute_hooks, f))

include("config.jl")
include("chunks.jl")
include("display_methods.jl")
include("readers.jl")
include("run.jl")
include("cache.jl")
include("formatters.jl")
include("Markdown2HTML.jl")
include("format.jl")
include("pandoc.jl")
include("writers.jl")

export weave, list_out_formats, tangle, convert_doc, notebook,
        set_chunk_defaults, get_chunk_defaults, restore_chunk_defaults,
        include_weave
end


text:
mvs = []
z = []

function go_to_corner_Sud_West!(r::Robot)
    while !isborder(r, Sud) || !isborder(r, West)
        while !isborder(r, Sud)
            move!(r, Sud)
        end
        while !isborder(r, West)
            move!(r, West)
        end
    end
end

function sledstroki!(r::Robot)
    while mvs != []
        putmarker!(r)
        pop!(mvs)
        move!(r, Ost)
    end
    while !isborder(r, West)
        push!(mvs, Ost)
        move!(r, West)
    end
end

function pervaiastroka!(r::Robot)
    while !isborder(r, Ost)
        putmarker!(r)
        push!(mvs, Ost)
        push!(z, Ost)
        move!(r, Ost)
    end
    putmarker!(r)
    push!(mvs, Ost)
    push!(z, Ost)
    go_to_corner_Sud_West!(r)
end

function main(r::Robot)
    go_to_corner_Sud_West!(r)
    pervaiastroka!(r)
    if !isborder(r, Nord)
        move!(r, Nord)
    end
    while mvs != 0
        pop!(mvs)
        sledstroki!(r)
        if !isborder(r, Nord)
            move!(r, Nord)
        else
            break
        end
    end
end

main(r)


text:
#Capacity Market Clearing Module

"""
This function does nothing if the product is not Capacity
"""
function update_capacity_supply_curve!(capacity_supply_curve::Vector{Vector{Union{String, Float64}}},
                                       product::T,
                                       project::P) where {T <: Product, P <: Project{<: BuildPhase}}

   return capacity_supply_curve
end

"""
This function pushes project bids in the capacity supply curve
"""
function update_capacity_supply_curve!(capacity_supply_curve::Vector{Vector{Union{String, Float64}}},
                                       product::Capacity,
                                       project::P) where P <: Project{<: BuildPhase}

   # Each Element of the Supply Curve:
   #   [1] - Project Name
   #   [2] - Project Size
   #   [3] - Project Capacity Bid
   #   [4] - Project Derating Factor

   push!(capacity_supply_curve, [get_name(project),
                                 get_maxcap(project),
                                 get_project_capacity_market_bid(project),
                                 get_project_derating(project)])
   return capacity_supply_curve
end

"""
This function creates the capacity market demand curve for actual capacity market clearing.
"""
function create_capacity_demand_curve(input_file::String,
                                      system_peak_load::Float64,
                                      capacity_mkt_bool::Bool)

   # Gather parameter data
   capacity_demand_params = read_data(input_file)[1, :]
   fpr = capacity_demand_params["Forecast Pool Req"] #Forecast Pool Req
   rel_req = system_peak_load * fpr * capacity_mkt_bool # Reliability Requirement
   irm = capacity_demand_params["IRM"] # Installed Reserve Margin
   irm_perc_points = parse.(Float64, split(capacity_demand_params["IRM perc points"], ";")) #Installed Reserve Margin percentage points

   net_CONE = capacity_demand_params["Net CONE per day"] * 365 * capacity_mkt_bool # Net CONE
   net_CONE_perc_points = parse.(Float64, split(capacity_demand_params["Net CONE perc points"], ";")) #Net CONE percentage points

   @assert length(net_CONE_perc_points) == length(irm_perc_points)

   gross_CONE_value = capacity_demand_params["Gross CONE per day"] * 365 * capacity_mkt_bool  # Gross CONE
   gross_CONE_points = zeros(length(net_CONE_perc_points))
   gross_CONE_points[1] = gross_CONE_value

   max_clear = capacity_demand_params["Max Clear"]

   # Construct demand curve
   demand_curve_break_points = zeros(length(net_CONE_perc_points) + 2)
   demand_curve_price_points = zeros(length(net_CONE_perc_points) + 2)

   demand_curve_break_points[1] = 0.0
   demand_curve_price_points[1] = max(net_CONE * net_CONE_perc_points[1], gross_CONE_points[1])

   demand_curve_break_points[end] = rel_req * max_clear
   demand_curve_price_points[end] = 0.0

   # Demand curve points based on PJM capacity market
   for i in 1:(length(net_CONE_perc_points))
      demand_curve_break_points[i + 1] = rel_req * ( 1 + irm + irm_perc_points[i]) / ( 1 + irm)
      demand_curve_price_points[i + 1] = max(net_CONE * net_CONE_perc_points[i], gross_CONE_points[i])
   end

   capacity_demand_curve = CapacityMarket(demand_curve_break_points, demand_curve_price_points)
   return capacity_demand_curve
end

"""
This function models actual capacity market clearing.
"""
function capacity_market_clearing(demand_curve::CapacityMarket,
                                  supply_curve::Vector{Vector{Union{String, Float64}}},
                                  solver::JuMP.MOI.OptimizerWithAttributes)

   # Each Element of the Supply Curve:
   #   [1] - Project Name
   #   [2] - Project Size
   #   [3] - Project Capacity Bid
   #   [4] - Project Derating Factor

   demand_segmentsize, demand_segmentgrad, demand_pricepoints = make_capacity_demand(demand_curve)

   #Number of segments of capacity supply and demand curves
   n_demand_seg = length(demand_segmentsize);
   n_supply_seg = length(supply_curve);

   #----------Capacity Market Clearing Problem---------------------------------------------------------------------------------
   cap_mkt = JuMP.Model(solver);

   #Define the variables
   JuMP.@variables(cap_mkt, begin
          Q_supply[s=1:n_supply_seg] >= 0 # Quantity of cleared capacity supply offers
          Q_demand[d=1:n_demand_seg] >= 0 # Quantity of cleared capacity demand from the demand curve
       end)

   #Functions----------------------------------------------------------------------------------------

   #Cost of procuring capacity supply for each segment
   supply_cost(Q_seg, s) = supply_curve[s][3] * Q_seg;

   #Welfare from meeting the capacity resource requirement for each segment
   demand_welfare(Q_seg, d) = demand_pricepoints[d] * Q_seg + 0.5 * demand_segmentgrad[d] * (Q_seg^2);

   #Expressions---------------------------------------------------------------------------------------

   #Totoal Cost of procuring capacity supply
   JuMP.@expression(cap_mkt, total_supply_cost, sum(supply_cost(Q_supply[s], s) for s = 1:n_supply_seg))

   #Total Welfare from meeting the capacity resource requirement
   JuMP.@expression(cap_mkt, total_dem_welfare, sum(demand_welfare(Q_demand[d], d) for d = 1:n_demand_seg))

   #Constraints--------------------------------------------------------------------------------------

   #Cleared capacity supply limit for each segment
   JuMP.@constraint(cap_mkt, [s=1:n_supply_seg], Q_supply[s] <= supply_curve[s][2])

   #Cleared demand limit for each segment
   JuMP.@constraint(cap_mkt, [d=1:n_demand_seg], Q_demand[d] <= demand_segmentsize[d])

   #Total cleared capacity supply should meet the total cleared demand
   JuMP.@constraint(cap_mkt, mkt_clear, sum(Q_supply[s] * supply_curve[s][4] for s in 1:n_supply_seg) == sum(Q_demand))

   #Define Objective Function - Social Welfare Maxmization
   JuMP.@objective(cap_mkt, Max, total_dem_welfare - total_supply_cost)

   println("Actual Capacity Market Clearing:")
   JuMP.optimize!(cap_mkt)
   println(JuMP.termination_status(cap_mkt))
   println(JuMP.objective_value(cap_mkt))

   #Capacity Market Clearing Price is the shadow variable of the capacity balance constraint
   cap_price = AxisArrays.AxisArray(reshape([JuMP.dual(mkt_clear)], 1,), [1])
   cap_accepted_bid = Dict(supply_curve[s][1] => value.(Q_supply[s]) / supply_curve[s][2]  for s in 1:n_supply_seg)
#------------------------------------------------------------------------------------------------
   return cap_price, cap_accepted_bid

end