text:
export GammaMixture, GammaMixtureNode, GammaMixtureNodeMetadata

# Gamma Mixture Functional Form
struct GammaMixture{N} end

# Special node
# Generic FactorNode implementation does not work with dynamic number of inputs
# We need to reimplement the following set of functions
# functionalform(factornode::FactorNode)
# sdtype(factornode::FactorNode)
# interfaces(factornode::FactorNode)
# factorisation(factornode::FactorNode)
# localmarginals(factornode::FactorNode)
# localmarginalnames(factornode::FactorNode)
# metadata(factornode::FactorNode)
# get_pipeline_stages(factornode::FactorNode)
#
# setmarginal!(factornode::FactorNode, cname::Symbol, marginal)
# getmarginal!(factornode::FactorNode, localmarginal::FactorNodeLocalMarginal)
#
# functional_dependencies(factornode::FactorNode, iindex::Int)
# get_messages_observable(factornode, message_dependencies)
# get_marginals_observable(factornode, marginal_dependencies)
#
# score(::Type{T}, ::FactorBoundFreeEnergy, ::Stochastic, node::AbstractFactorNode, scheduler) where T
#
# Base.show
#
const GammaMixtureNodeFactorisationSupport = Union{MeanField, }

struct GammaMixtureNode{N, F <: GammaMixtureNodeFactorisationSupport, M, P} <: AbstractFactorNode
    factorisation :: F

    # Interfaces
    out    :: NodeInterface
    switch :: NodeInterface
    as  :: NTuple{N, IndexedNodeInterface}
    bs  :: NTuple{N, IndexedNodeInterface}

    meta     :: M
    pipeline :: P
end

functionalform(factornode::GammaMixtureNode{N}) where N = GammaMixture{N}
sdtype(factornode::GammaMixtureNode)                    = Stochastic()
interfaces(factornode::GammaMixtureNode)                = (factornode.out, factornode.switch, factornode.as..., factornode.bs...)
factorisation(factornode::GammaMixtureNode)             = factornode.factorisation
localmarginals(factornode::GammaMixtureNode)            = error("localmarginals() function is not implemented for GammaMixtureNode")
localmarginalnames(factornode::GammaMixtureNode)        = error("localmarginalnames() function is not implemented for GammaMixtureNode")
metadata(factornode::GammaMixtureNode)                  = factornode.meta
getpipeline(factornode::GammaMixtureNode)               = factornode.pipeline

setmarginal!(factornode::GammaMixtureNode, cname::Symbol, marginal)                = error("setmarginal() function is not implemented for GammaMixtureNode")
getmarginal!(factornode::GammaMixtureNode, localmarginal::FactorNodeLocalMarginal) = error("getmarginal() function is not implemented for GammaMixtureNode")

## activate!

struct GammaMixtureNodeFunctionalDependencies <: AbstractNodeFunctionalDependenciesPipeline end

default_functional_dependencies_pipeline(::Type{ <: GammaMixture }) = GammaMixtureNodeFunctionalDependencies()

function functional_dependencies(::GammaMixtureNodeFunctionalDependencies, factornode::GammaMixtureNode{N, F}, iindex::Int) where { N, F <: MeanField }
    message_dependencies = ()

    marginal_dependencies = if iindex === 1
        (factornode.switch, factornode.as, factornode.bs)
    elseif iindex === 2
        (factornode.out, factornode.as, factornode.bs)
    elseif 2 < iindex <= N + 2
        (factornode.out, factornode.switch, factornode.bs[iindex - 2])
    elseif N + 2 < iindex <= 2N + 2
        (factornode.out, factornode.switch, factornode.as[iindex - N - 2])
    else
        error("Invalid index in functional_dependencies for GammaMixtureNode")
    end

    return message_dependencies, marginal_dependencies
end

function get_messages_observable(factornode::GammaMixtureNode{N, F}, message_dependencies::Tuple{}) where { N, F <: MeanField }
    return nothing, of(nothing)
end

function get_marginals_observable(
    factornode::GammaMixtureNode{N, F},
    marginal_dependencies::Tuple{ NodeInterface, NTuple{N, IndexedNodeInterface}, NTuple{N, IndexedNodeInterface} }) where { N, F <: MeanField }

    varinterface = marginal_dependencies[1]
    asinterfaces = marginal_dependencies[2]
    bsinterfaces = marginal_dependencies[3]

    marginal_names = Val{ (name(varinterface), name(asinterfaces[1]), name(bsinterfaces[1])) }
    marginals_observable = combineLatest((
        getmarginal(connectedvar(varinterface), IncludeAll()),
        combineLatest(map((rate) -> getmarginal(connectedvar(rate), IncludeAll()), reverse(bsinterfaces)), PushNew()),
        combineLatest(map((shape) -> getmarginal(connectedvar(shape), IncludeAll()), reverse(asinterfaces)), PushNew()),
    ), PushNew()) |> map_to((
        getmarginal(connectedvar(varinterface), IncludeAll()),
        map((shape) -> getmarginal(connectedvar(shape), IncludeAll()), asinterfaces),
        map((rate) -> getmarginal(connectedvar(rate), IncludeAll()), bsinterfaces)
    ))

    return marginal_names, marginals_observable
end

function get_marginals_observable(
    factornode::GammaMixtureNode{N, F},
    marginal_dependencies::Tuple{ NodeInterface, NodeInterface, IndexedNodeInterface }) where { N, F <: MeanField }

    outinterface    = marginal_dependencies[1]
    switchinterface = marginal_dependencies[2]
    varinterface    = marginal_dependencies[3]

    marginal_names       = Val{ (name(outinterface), name(switchinterface), name(varinterface)) }
    marginals_observable = combineLatestUpdates((
        getmarginal(connectedvar(outinterface), IncludeAll()),
        getmarginal(connectedvar(switchinterface), IncludeAll()),
        getmarginal(connectedvar(varinterface), IncludeAll())
    ), PushNew())

    return marginal_names, marginals_observable
end

# FreeEnergy related functions

@average_energy GammaMixture (q_out::Any, q_switch::Any, q_a::NTuple{N, Any}, q_b::NTuple{N, GammaShapeRate}) where N = begin
    z_bar = probvec(q_switch)
    return mapreduce(+, 1:N, init = 0.0) do i
        return z_bar[i] * score(AverageEnergy(), GammaShapeRate, Val{ (:out, :α , :β) }, map((q) -> Marginal(q, false, false), (q_out, q_a[i], q_b[i])), nothing)
    end
end

function score(::Type{T}, objective::BetheFreeEnergy, ::FactorBoundFreeEnergy, ::Stochastic, node::GammaMixtureNode{N, MeanField}, scheduler) where { T <: InfCountingReal, N }

    skip_strategy = marginal_skip_strategy(objective)

    stream = combineLatest((
        getmarginal(connectedvar(node.out), skip_strategy) |> schedule_on(scheduler),
        getmarginal(connectedvar(node.switch), skip_strategy) |> schedule_on(scheduler),
        combineLatest(map((as) -> getmarginal(connectedvar(as), skip_strategy) |> schedule_on(scheduler), node.as), PushNew()),
        combineLatest(map((bs) -> getmarginal(connectedvar(bs), skip_strategy) |> schedule_on(scheduler), node.bs), PushNew())
    ), PushNew())

    mapping = let fform = functionalform(node), meta = metadata(node)
        (marginals) -> begin
            average_energy   = score(AverageEnergy(), fform, Val{ (:out, :switch, :a, :b) }, marginals, meta)

            out_entropy     = score(DifferentialEntropy(), marginals[1])
            switch_entropy  = score(DifferentialEntropy(), marginals[2])
            a_entropies = mapreduce((m) -> score(DifferentialEntropy(), m), +, marginals[3])
            b_entropies = mapreduce((m) -> score(DifferentialEntropy(), m), +, marginals[4])

            return convert(T, average_energy - (out_entropy + switch_entropy + a_entropies + b_entropies))
        end
    end

    return stream |> map(T, mapping)
end

as_node_functional_form(::Type{ <: GammaMixture }) = ValidNodeFunctionalForm()

# Node creation related functions

sdtype(::Type{ <: GammaMixture }) = Stochastic()

collect_factorisation(::Type{ <: GammaMixture }, factorisation) = factorisation

function ReactiveMP.make_node(::Type{ <: GammaMixture{N} }; factorisation::F = MeanField(), meta::M = nothing, pipeline::P = nothing) where { N, F, M, P }
    @assert N >= 2 "GammaMixtureNode requires at least two mixtures on input"
    @assert typeof(factorisation) <: GammaMixtureNodeFactorisationSupport "GammaMixtureNode supports only following factorisations: [ $(GammaMixtureNodeFactorisationSupport) ]"
    out    = NodeInterface(:out)
    switch = NodeInterface(:switch)
    as   = ntuple((index) -> IndexedNodeInterface(index, NodeInterface(:a)), N)
    bs   = ntuple((index) -> IndexedNodeInterface(index, NodeInterface(:b)), N)
    meta = collect_meta(GammaMixture, meta)
    return GammaMixtureNode{N, F, typeof(meta), P}(factorisation, out, switch, as, bs, meta, pipeline)
end

function ReactiveMP.make_node(::Type{ <: GammaMixture }, out::AbstractVariable, switch::AbstractVariable, as::NTuple{N, AbstractVariable}, bs::NTuple{N, AbstractVariable}; factorisation = MeanField(), meta = nothing, pipeline = nothing) where { N}
    node = make_node(GammaMixture{N}, factorisation = collect_factorisation(GammaMixture, factorisation), meta = collect_meta(GammaMixture, meta), pipeline = collect_pipeline(GammaMixture, pipeline))

    # out
    out_index = getlastindex(out)
    connectvariable!(node.out, out, out_index)
    setmessagein!(out, out_index, messageout(node.out))

    # switch
    switch_index = getlastindex(switch)
    connectvariable!(node.switch, switch, switch_index)
    setmessagein!(switch, switch_index, messageout(node.switch))

    # as
    foreach(zip(node.as, as)) do (ainterface, avar)
        shape_index = getlastindex(avar)
        connectvariable!(ainterface, avar, shape_index)
        setmessagein!(avar, shape_index, messageout(ainterface))
    end

    # bs
    foreach(zip(node.bs, bs)) do (binterface, bvar)
        rate_index = getlastindex(bvar)
        connectvariable!(binterface, bvar, rate_index)
        setmessagein!(bvar, rate_index, messageout(binterface))
    end

    return node
end

function ReactiveMP.make_node(fform::Type{ <: GammaMixture }, autovar::AutoVar, args::Vararg{ <: AbstractVariable }; kwargs...)
    var  = randomvar(getname(autovar))
    node = make_node(fform, var, args...; kwargs...)
    return node, var
end

ext: jl


text:
#coast.msh geometry, created by JIGSAW
#adapted by @miguelraz 28/03/2020 from
#https://github.com/dengwirda/jigsaw-matlab
#edge2=303
const lakes_edges = [(0,1)
(1,2)
(2,3)
(3,4)
(4,5)
(5,6)
(6,7)
(7,8)
(8,9)
(9,10)
(10,11)
(11,12)
(12,13)
(13,14)
(14,15)
(15,16)
(16,17)
(17,18)
(18,19)
(19,20)
(20,21)
(21,22)
(22,23)
(23,24)
(24,25)
(25,26)
(26,27)
(27,28)
(28,29)
(29,30)
(30,31)
(31,32)
(32,33)
(33,34)
(34,35)
(35,36)
(36,37)
(37,38)
(38,39)
(39,40)
(40,41)
(41,42)
(42,43)
(43,44)
(44,45)
(45,46)
(46,47)
(47,48)
(48,49)
(49,50)
(50,51)
(51,52)
(52,53)
(53,54)
(54,55)
(55,56)
(56,57)
(57,58)
(58,59)
(59,60)
(60,61)
(61,62)
(62,63)
(63,64)
(64,65)
(65,66)
(66,67)
(67,68)
(68,69)
(69,70)
(70,71)
(71,72)
(72,73)
(73,74)
(74,75)
(75,76)
(76,77)
(77,78)
(78,79)
(79,80)
(80,81)
(81,82)
(82,83)
(83,84)
(84,85)
(85,86)
(86,87)
(87,88)
(88,89)
(89,90)
(90,91)
(91,92)
(92,93)
(93,94)
(94,95)
(95,96)
(96,97)
(97,98)
(98,99)
(99,100)
(100,101)
(101,102)
(102,103)
(103,104)
(104,105)
(105,106)
(106,107)
(107,108)
(108,109)
(109,110)
(110,111)
(111,112)
(112,113)
(113,114)
(114,115)
(115,116)
(116,117)
(117,118)
(118,119)
(119,120)
(120,121)
(121,122)
(122,123)
(123,124)
(124,125)
(125,126)
(126,127)
(127,128)
(128,129)
(129,130)
(130,131)
(131,132)
(132,133)
(133,134)
(134,135)
(135,136)
(136,137)
(137,138)
(138,139)
(139,140)
(140,141)
(141,142)
(142,143)
(143,144)
(144,145)
(145,146)
(146,147)
(147,148)
(148,149)
(149,150)
(150,151)
(151,152)
(152,153)
(153,154)
(154,155)
(155,156)
(156,157)
(157,158)
(158,159)
(159,160)
(160,161)
(161,162)
(162,163)
(163,164)
(164,165)
(165,166)
(166,167)
(167,168)
(168,169)
(169,170)
(170,171)
(171,172)
(172,173)
(173,174)
(174,175)
(175,176)
(176,177)
(177,178)
(178,179)
(179,180)
(180,181)
(181,182)
(182,183)
(183,184)
(184,185)
(185,186)
(186,187)
(187,188)
(188,189)
(189,190)
(190,191)
(191,192)
(192,193)
(193,194)
(194,195)
(195,196)
(196,197)
(197,198)
(198,199)
(199,200)
(200,201)
(201,202)
(202,203)
(203,204)
(204,205)
(205,206)
(206,207)
(207,208)
(208,209)
(209,210)
(210,211)
(211,212)
(212,213)
(213,214)
(214,215)
(215,216)
(216,217)
(217,218)
(218,219)
(219,220)
(220,221)
(221,222)
(222,223)
(223,224)
(224,0)
(225,226)
(226,227)
(227,228)
(228,229)
(229,230)
(230,231)
(231,232)
(232,233)
(233,234)
(234,235)
(235,236)
(236,225)
(237,238)
(238,239)
(239,240)
(240,241)
(241,242)
(242,243)
(243,244)
(244,237)
(245,246)
(246,247)
(247,248)
(248,249)
(249,250)
(250,251)
(251,252)
(252,245)
(253,254)
(254,255)
(255,256)
(256,257)
(257,258)
(258,259)
(259,260)
(260,261)
(261,262)
(262,263)
(263,264)
(264,265)
(265,266)
(266,267)
(267,268)
(268,269)
(269,270)
(270,271)
(271,272)
(272,273)
(273,274)
(274,275)
(275,253)
(276,277)
(277,278)
(278,279)
(279,280)
(280,281)
(281,282)
(282,283)
(283,284)
(284,285)
(285,286)
(286,287)
(287,288)
(288,289)
(289,290)
(290,291)
(291,276)
(292,293)
(293,294)
(294,295)
(295,296)
(296,297)
(297,298)
(298,299)
(299,300)
(300,301)
(301,302)
(302,292)]


text:
export CSRGraph
export IndexType, LOCAL_INDICES, GLOBAL_INDICES
export StorageStatus, STORAGE_2D, STORAGE_1D_UNPACKED, STORAGE_1D_PACKED
export ProfileType, STATIC_PROFILE, DYNAMIC_PROFILE

"""
Can be used to differentiate global and local indices
"""
@enum IndexType LOCAL_INDICES GLOBAL_INDICES

"""
Status of the graph's or matrix's storage, when not in
a fill-complete state.
"""
@enum StorageStatus STORAGE_2D STORAGE_1D_UNPACKED STORAGE_1D_PACKED

"""
Allocation profile for matrix/graph entries
"""
@enum ProfileType STATIC_PROFILE DYNAMIC_PROFILE

#TODO document the type and constructors

mutable struct CSRGraph{GID <: Integer, PID <: Integer, LID <: Integer} <: RowGraph{GID, PID, LID}
    rowMap::BlockMap{GID, PID, LID}
    colMap::Union{BlockMap{GID, PID, LID}, Nothing}
    rangeMap::Union{BlockMap{GID, PID, LID}, Nothing}
    domainMap::Union{BlockMap{GID, PID, LID}, Nothing}

    #may be null if domainMap and colMap are the same
    importer::Union{Import{GID, PID, LID}, Nothing}
    #may be null if rangeMap and rowMap are the same
    exporter::Union{Export{GID, PID, LID}, Nothing}

    localGraph::LocalCSRGraph{LID, LID}

    #Local number of (populated) entries; must always be consistent
    nodeNumEntries::LID

    #Local number of (populated) diagonal entries.
    nodeNumDiags::LID

    #Local maximum of the number of entries in each row.
    nodeMaxNumRowEntries::LID

    #Global number of entries in the graph.
    globalNumEntries::GID

    #Global number of (populated) diagonal entries.
    globalNumDiags::GID

    #Global maximum of the number of entries in each row.
    globalMaxNumRowEntries::GID

    #Whether the graph was allocated with static or dynamic profile.
    pftype::ProfileType

    ## 1-D storage (Static profile) data structures ##
    localIndices1D::Array{LID, 1}
    globalIndices1D::Array{GID, 1}
    rowOffsets::Array{LID, 1}  #Tpetra: k_rowPts_

    ## 2-D storage (Dynamic profile) data structures ##
    localIndices2D::Array{Array{LID, 1}, 1}
    globalIndices2D::Array{Array{GID, 1}, 1}
    #may exist in 1-D storage if not packed
    numRowEntries::Array{LID, 1}

    storageStatus::StorageStatus

    indicesAllowed::Bool
    indicesType::IndexType
    fillComplete::Bool

    lowerTriangle::Bool
    upperTriangle::Bool
    indicesAreSorted::Bool
    noRedundancies::Bool
    haveLocalConstants::Bool
    haveGlobalConstants::Bool
    sortGhostsAssociatedWithEachProcessor::Bool

    plist::Dict

    nonLocals::Dict{GID, Array{GID, 1}}

    #Large ammounts of duplication between the constructors, so group it in an inner constructor
    function CSRGraph(
        rowMap::BlockMap{GID, PID, LID},
        colMap::Union{BlockMap{GID, PID, LID}, Nothing},
        rangeMap::Union{BlockMap{GID, PID, LID}, Nothing},
        domainMap::Union{BlockMap{GID, PID, LID}, Nothing},

        localGraph::LocalCSRGraph,

        nodeNumEntries::LID,

        pftype::ProfileType,
        storageStatus::StorageStatus,

        indicesType::IndexType,
        plist::Dict
    ) where {GID <: Integer, PID <: Integer, LID <: Integer}

        graph = new{GID, PID, LID}(
            rowMap,
            colMap,
            rangeMap,
            domainMap,

            nothing,
            nothing,

            localGraph,

            #Local number of (populated) entries; must always be consistent
            nodeNumEntries,

            #using 0 to indicate uninitiallized, since -1 isn't gareenteed to work
            0, #nodeNumDiags
            0, #nodeMaxNumRowEntries
            0, #globalNumEntries
            0, #globalNumDiags
            0, #globalMaxNumRowEntries

            #Whether the graph was allocated with static or dynamic profile.
            pftype,

            ## 1-D storage (Static profile) data structures ##
            LID[],
            GID[],
            LID[],

            ## 2-D storage (Dynamic profile) data structures ##
            Array{Array{LID, 1}, 1}(undef, 0),
            Array{Array{GID, 1}, 1}(undef, 0),
            LID[],

            storageStatus,

            false,
            indicesType,
            false,

            false,
            false,
            true,
            true,
            false,
            false,
            true,

            plist,

            Dict{GID, Array{GID, 1}}()
        )

        ## staticAssertions()
        #skipping sizeof checks
        #skipping max value checks related to size_t

        @assert(indicesType != LOCAL_INDICES || colMap != nothing,
            "Cannot have local indices with a null column Map")

        #ensure LID is a subset of GID (for positive numbers)
        if !(LID <: GID) && (GID != BigInt) && (GID != Integer)
            # all ints are assumed to be able to handle 1, up to their max
            if LID == BigInt || LID == Integer || typemax(LID) > typemax(GID)
                throw(InvalidArgumentError("The positive values of GID must "
                        * "be a superset of the positive values of LID"))
            end
        end

        graph
    end
end

#### Constructors #####

function CSRGraph(rowMap::BlockMap{GID, PID, LID}, maxNumEntriesPerRow::Integer,
        pftype::ProfileType; plist...) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, LID(maxNumEntriesPerRow),  pftype, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, maxNumEntriesPerRow::Integer,
        pftype::ProfileType, plist::Dict) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, nothing, LID(maxNumEntriesPerRow), pftype, plist)
end

#= Should be handled by the following constructors
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        maxNumEntriesPerRow::Integer, pftype::ProfileType; plist...) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, LID(maxNumEntriesPerRow), pftype, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        maxNumEntriesPerRow::Integer, pftype::ProfileType, plist::Dict) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, LID(maxNumEntriesPerRow), pftype, plist)
end
=#

function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::Union{BlockMap{GID, PID, LID}, Nothing},
        maxNumEntriesPerRow::Integer, pftype::ProfileType; plist...) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, maxNumEntriesPerRow, pftype, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::Union{BlockMap{GID, PID, LID}, Nothing},
        maxNumEntriesPerRow::Integer, pftype::ProfileType, plist::Dict) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    graph = CSRGraph(
        rowMap,
        colMap,
        nothing,
        nothing,

        LocalCSRGraph{LID, LID}(), #localGraph

        LID(0), #nodeNumEntries

        pftype,

        (pftype == STATIC_PROFILE ?
              STORAGE_1D_UNPACKED
            : STORAGE_2D),

        (colMap == nothing) ? GLOBAL_INDICES : LOCAL_INDICES,
        plist
    )
    allocateIndices(graph, graph.indicesType, LID(maxNumEntriesPerRow))

    resumeFill(graph, plist)
    checkInternalState(graph)

    graph
end

function CSRGraph(rowMap::BlockMap{GID, PID, LID}, numEntPerRow::AbstractArray{<:Integer, 1},
        pftype::ProfileType; plist...)  where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, numEntPerRow, pftype, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, numEntPerRow::AbstractArray{<:Integer, 1},
        pftype::ProfileType, plist::Dict)  where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, nothing, numEntPerRow, pftype, plist)
end

function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        numEntPerRow::AbstractArray{<:Integer, 1}, pftype::ProfileType;
        plist...)  where {GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, numEntPerRow, pftype, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        numEntPerRow::AbstractArray{<:Integer, 1}, pftype::ProfileType,
        plist::Dict)  where {GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, numEntPerRow, pftype, plist)
end

function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::Union{BlockMap{GID, PID, LID}, Nothing},
        numEntPerRow::AbstractArray{<:Integer, 1}, pftype::ProfileType;
        plist...)  where {GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, numEntPerRow, pftype, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::Union{BlockMap{GID, PID, LID}, Nothing},
        numEntPerRow::AbstractArray{<:Integer, 1}, pftype::ProfileType,
        plist::Dict)  where {GID <: Integer, PID <: Integer, LID <: Integer}
    graph = CSRGraph(
        rowMap,
        colMap,
        nothing,
        nothing,

        LocalCSRGraph{LID, LID}(), #localGraph

        LID(0), #nodeNumEntries

        #Whether the graph was allocated with static or dynamic profile.
        pftype,

        (pftype == STATIC_PROFILE ?
              STORAGE_1D_UNPACKED
            : STORAGE_2D),

        (colMap == nothing) ? GLOBAL_INDICES : LOCAL_INDICES,
        plist
    )

    localNumRows = numMyElements(rowMap)
    if length(numEntPerRow) != localNumRows
        throw(InvalidArgumentError("numEntPerRows has length $(length(numEntPerRow)) " *
                "!= the local number of rows $lclNumRows as spcified by the input row Map"))
    end

    allocateIndices(graph, graph.indicesType, Array{LID, 1}(numEntPerRow))

    resumeFill(graph, plist)
    checkInternalState(graph)

    graph
end

function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        rowPointers::AbstractArray{LID, 1}, columnIndices::Array{LID, 1};
        plist...) where {GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, rowPointers, columnIndices, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        rowPointers::AbstractArray{LID, 1}, columnIndices::Array{LID, 1},
        plist::Dict) where {GID <: Integer, PID <: Integer, LID <: Integer}
    graph = CSRGraph(
        rowMap,
        colMap,
        nothing,
        nothing,

        LocalCSRGraph{LID, LID}(), #localGraph

        LID(0), #nodeNumEntries

        STATIC_PROFILE,

        STORAGE_1D_PACKED,

        LOCAL_INDICES,
        plist
    )
    #seems to be already taken care of
    #allocateIndices(graph, LOCAL_INDICES, numEntPerRow)

    setAllIndicies(graph, rowPointers, columnIndicies)
    checkInternalState(graph)

    graph
end

function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        localGraph::LocalCSRGraph{LID, LID}; plist...) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    CSRGraph(rowMap, colMap, localGraph, Dict(plist))
end
function CSRGraph(rowMap::BlockMap{GID, PID, LID}, colMap::BlockMap{GID, PID, LID},
        localGraph::LocalCSRGraph{LID, LID}, plist::Dict) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    mapRowCount = numMyElements(rowMap)
    graph = CSRGraph(
        rowMap,
        colMap,
        rowMap,
        colMap,

        localGraph,

        localGraph.rowMap[mapRowCount+1], #nodeNumEntries

        STATIC_PROFILE,

        STORAGE_1D_PACKED,

        LOCAL_INDICES,
        plist
    )

    if numRows(localGraph) != numMyElements(rowMap)
        throw(InvalidArgumentError("input row map and input local "
                * "graph need to have the same number of rows.  The "
                * "row map claims $(numMyElements(rowMap)) row(s), "
                * "but the local graph claims $(numRows(localGraph)) "
                * "row(s)."))
    end

    #seems to be already taken care of
    #allocateIndices(graph, LOCAL_INDICES, numEntPerRow)

    makeImportExport(graph)

    d_inds = localGraph.entries
    graph.localIndices1D = d_inds

    d_ptrs = localGraph.rowMap
    graph.rowOffsets = d_ptrs

    #reset local properties
    graph.upperTriangle = true
    graph.lowerTriangle = true
    graph.nodeMaxNumRowEntries = 0
    graph.nodeNumDiags

    for localRow = 1:mapRowCount
        globalRow = gid(rowMap, localRow)
        rowLID = lid(colMap, globalRow)

        #possible that the local matrix has no entries in the column
        #corrisponding to the current row, in that case, the column map
        #might not contain that GID.  Hence, the index validity check
        if rowLID != 0
            if rowLID +1 > length(d_ptrs)
                throw(InvalidArgumentError("The given row Map and/or column Map "
                        * "is/are not compatible with the provided local graphs."))
            end
            if d_ptrs[rowLID] != d_ptr[rowLID+1]
                smallestCol = d_inds[d_ptrs[rowLID]]
                largestCol  = d_inds[d_ptrs[rowLID+1]-1]

                if smallestCol < localRow
                    graph.upperTriangle = false
                end
                if localRow < largestCol
                    graph.lowerTriangle = false
                end

                if rowLID in d_inds[d_ptrs[rowLID]:d_ptrs[rowLID]-1]
                    graph.nodeNumDiags += 1
                end
            end

            graph.nodeMaxNumRowEntries = max((d_ptrs[rowLID + 1] - d_ptrs[rowLID]),
                                            graph.nodeMaxNumRowEntries)
        end
    end

    graph.hasLocalConstants = true
    computeGlobalConstants(graph)

    graph.fillComplete = true
    checkInternalState(graph)

    graph
end


text:
mutable struct MarketModelPathwiseInverseFloater <: MarketModelPathwiseMultiProduct
  rateTimes::Vector{Float64}
  fixedAccruals::Vector{Float64}
  floatingAccruals::Vector{Float64}
  fixedStrikes::Vector{Float64}
  fixedMultipliers::Vector{Float64}
  floatingSpreads::Vector{Float64}
  paymentTimes::Vector{Float64}
  payer::Bool
  multiplier::Float64
  lastIndex::Int
  evolution::EvolutionDescription
  currentIndex::Int
end

function MarketModelPathwiseInverseFloater(rateTimes::Vector{Float64},
                                          fixedAccruals::Vector{Float64},
                                          floatingAccruals::Vector{Float64},
                                          fixedStrikes::Vector{Float64},
                                          fixedMultipliers::Vector{Float64},
                                          floatingSpreads::Vector{Float64},
                                          paymentTimes::Vector{Float64},
                                          payer::Bool)

  multiplier = payer ? -1.0 : 1.0
  lastIndex = length(rateTimes)

  check_increasing_times(paymentTimes)

  length(fixedAccruals) == lastIndex - 1 || error("Incorrect number of fixedAccruals given")
  length(floatingAccruals) == lastIndex - 1 || error("Incorrect number of floatingAccruals given")
  length(fixedStrikes) == lastIndex - 1 || error("Incorrect number of fixedStrikes given")
  length(fixedMultipliers) == lastIndex - 1 || error("Incorrect number of fixedMultipliers given")
  length(floatingSpreads) == lastIndex - 1 || error("Incorrect number of floatingSpreads given")
  length(paymentTimes) == lastIndex - 1 || error("Incorrect number of paymentTimes given")

  evolTimes = copy(rateTimes)
  pop!(evolTimes)

  evolution = EvolutionDescription(rateTimes, evolTimes)

  return MarketModelPathwiseInverseFloater(rateTimes, fixedAccruals, floatingAccruals, fixedStrikes, fixedMultipliers, floatingSpreads,
                                        paymentTimes, payer, multiplier, lastIndex, evolution, 1)
end

number_of_products(::MarketModelPathwiseInverseFloater) = 1
get_evolution(mm::MarketModelPathwiseInverseFloater) = mm.evolution
possible_cash_flow_times(mm::MarketModelPathwiseInverseFloater) = mm.paymentTimes
max_number_of_cashflows_per_product_per_step(mm::MarketModelPathwiseInverseFloater) = 1
already_deflated(::MarketModelPathwiseInverseFloater) = false

reset!(mm::MarketModelPathwiseInverseFloater) = mm.currentIndex = 1

function next_time_step!(mm::MarketModelPathwiseInverseFloater, currentState::CurveState, numberCashFlowsThisStep::Vector{Int},
                        cashFlowsGenerated::Vector{Vector{MarketModelPathWiseCashFlow}})
  numberCashFlowsThisStep[1] = 1
  @simd for i = 2:mm.lastIndex
    @inbounds cashFlowsGenerated[1][1].amount[i] = 0.0
  end

  liborRate = forward_rate(currentState, mm.currentIndex)
  inverseFloatingCoupon = max((mm.fixedStrikes[mm.currentIndex] - mm.fixedMultipliers[mm.currentIndex] * liborRate), 0.0) * mm.fixedAccruals[mm.currentIndex]
  floatingCoupon = (liborRate + mm.floatingSpreads[mm.currentIndex]) * mm.fixedAccruals[mm.currentIndex]

  cashFlowsGenerated[1][1].timeIndex = mm.currentIndex
  cashFlowsGenerated[1][1].amount[1] = mm.multiplier * (inverseFloatingCoupon - floatingCoupon)

  if inverseFloatingCoupon > 0.0
    cashFlowsGenerated[1][1].amount[mm.currentIndex + 1] = mm.multiplier * ( -mm.fixedMultipliers[mm.currentIndex] * mm.fixedAccruals[mm.currentIndex] -
                                                          mm.floatingAccruals[mm.currentIndex])
  else
    cashFlowsGenerated[1][1].amount[mm.currentIndex + 1] = -mm.multiplier * mm.floatingAccruals[mm.currentIndex]
  end

  mm.currentIndex += 1

  return mm.currentIndex == mm.lastIndex
end

## Clone ##
clone(mm::MarketModelPathwiseInverseFloater) = MarketModelPathwiseInverseFloater(copy(mm.rateTimes), copy(mm.fixedAccruals), copy(mm.floatingAccruals), copy(mm.fixedStrikes),
                                                copy(mm.fixedMultipliers), copy(mm.floatingSpreads), copy(mm.paymentTimes), mm.payer, mm.multiplier,
                                                mm.lastIndex, clone(mm.evolution), mm.currentIndex)


text:
module TestCLIBucket

using Base.Test
using SimpleTasks.Services.CLIBucket

import AWS
import SimpleTasks.Services.Bucket
import SimpleTasks.Services.CLIBucket
import SimpleTasks.Services.AWSCLIProvider
import SimpleTasks.Services.GCSCLIProvider

const TEST_FILE_NAME = "testfile$(rand())"
const BUCKET_NAME = "seunglab-test"
const MAX_INDEX = 100000

macro silent(expr::Expr)
    return :(try run(pipeline($expr, stdout=DevNull, stderr=DevNull)) catch end)
end

function upload_remote_test_file(provider::CLIBucket.Provider)
    create_local_test_file(TEST_FILE_NAME)
    try
        run(`$(provider.command) cp $TEST_FILE_NAME
            $(provider.prefix)/$BUCKET_NAME/$TEST_FILE_NAME`)
    finally
        delete_local_test_file()
    end
end

function delete_remote_test_file(provider::CLIBucket.Provider)
    run(`$(provider.command) rm $(provider.prefix)/$BUCKET_NAME/$TEST_FILE_NAME`)
end

function write_numbers(io::IO)
    for i in 1:MAX_INDEX
        write(io, "$i\n")
    end
end

function create_local_test_file(filename::AbstractString)
    file = open(filename, "w")
    write_numbers(file)
    close(file)
end

function delete_local_test_file()
    rm(TEST_FILE_NAME)
end

function test_creatable(provider::CLIBucket.Provider)
    bucket = CLIBucketService(provider, BUCKET_NAME)
    @test bucket != nothing
end

function test_bad_bucket(provider::CLIBucket.Provider)
    @test_throws ArgumentError CLIBucketService(provider,
        "$BUCKET_NAME-bad")
end

function test_download_IO(provider::CLIBucket.Provider)
    upload_remote_test_file(provider)

    bucket = CLIBucketService(provider, BUCKET_NAME)

    buffer = IOBuffer()
    Bucket.download(bucket, TEST_FILE_NAME, buffer)

    file_indices = split(takebuf_string(buffer), "\n")
    for index in 1:length(file_indices) - 1
        @test parse(Int, file_indices[index]) == index
    end
    # last split from \n is generates an empty string
    @test length(file_indices) - 1 == MAX_INDEX

    delete_remote_test_file(provider)
end

function test_download_file(provider::CLIBucket.Provider)
    upload_remote_test_file(provider)

    bucket = CLIBucketService(provider, BUCKET_NAME)

    download_filename = "new$TEST_FILE_NAME"
    Bucket.download(bucket, TEST_FILE_NAME, download_filename)

    file = open(download_filename, "r")
    index = 0
    try
        while !eof(file)
            index = index + 1
            @test parse(Int, readline(file)) == index
        end
    finally
        close(file)
        rm(download_filename)
    end
    @test index == MAX_INDEX

    delete_remote_test_file(provider)
end

function test_download_stream(provider::CLIBucket.Provider)
    upload_remote_test_file(provider)

    bucket = CLIBucketService(provider, BUCKET_NAME)

    stream = Bucket.download(bucket, TEST_FILE_NAME)

    index = 0
    try
        while !eof(stream)
            index = index + 1
            @test parse(Int, readline(stream)) == index
        end
    finally
        close(stream)
    end
    @test index == MAX_INDEX

    delete_remote_test_file(provider)
end

function test_upload_io(provider::CLIBucket.Provider)
    bucket = CLIBucketService(provider, BUCKET_NAME)

    io = IOBuffer()
    write_numbers(io)
    seekstart(io)
    # also tests parentheses commas and underscores
    upload_filename = "$(TEST_FILE_NAME)UpIO(,_)"

    # try to remove the file we are tring to upload from bucket
    @silent `$(provider.command) rm $(provider.prefix)/$BUCKET_NAME/$upload_filename`

    Bucket.upload(bucket, io, upload_filename)

    # Verify the uploaded file by manually download the uploaded file
    try
        run(`$(provider.command) cp $(provider.prefix)/$BUCKET_NAME/$upload_filename $upload_filename`)
    catch
        error("Unable to find downloaded file $upload_filename")
    end

    downloaded_file = open(upload_filename, "r")
    index = 0
    try
        while !eof(downloaded_file)
            index = index + 1
            @test parse(Int, readline(downloaded_file)) == index
        end
    finally
        close(downloaded_file)
        rm(upload_filename)
    end
    @test index == MAX_INDEX
    @silent `$(provider.command) rm $(provider.prefix)/$BUCKET_NAME/$upload_filename`
end

function test_upload_file(provider::CLIBucket.Provider)
    bucket = CLIBucketService(provider, BUCKET_NAME)

    upload_filename = "$(TEST_FILE_NAME)UpFile"
    # try to remove the file we are tring to upload from bucket
    @silent `$(provider.command) rm $(provider.prefix)/$BUCKET_NAME/$upload_filename`

    create_local_test_file(upload_filename)

    Bucket.upload(bucket, upload_filename, upload_filename)

    # Verify the uploaded file by manually download the uploaded file
    try
        run(`$(provider.command) cp
            $(provider.prefix)/$BUCKET_NAME/$upload_filename $upload_filename`)
    catch
        error("Unable to find downloaded file $upload_filename")
    end

    downloaded_file = open(upload_filename, "r")
    index = 0
    try
        while !eof(downloaded_file)
            index = index + 1
            @test parse(Int, readline(downloaded_file)) == index
        end
    finally
        close(downloaded_file)
        rm(upload_filename)
    end
    @test index == MAX_INDEX
    @silent `$(provider.command) rm $(provider.prefix)/$BUCKET_NAME/$upload_filename`
end

function test_delete_file_exists(provider::CLIBucket.Provider)
    upload_remote_test_file(provider)

    bucket = CLIBucketService(provider, BUCKET_NAME)
    Bucket.delete(bucket, TEST_FILE_NAME)

    # Verify that the remote file was deleted
    try
        @test_throws ErrorException run(`$(provider.command) ls
                $(provider.prefix)/$BUCKET_NAME/$TEST_FILE_NAME`)
    finally
        # Remove the test file if we had been unable to remove the file
        @silent `$(provider.command) rm $(provider.prefix)/$BUCKET_NAME/$TEST_FILE_NAME`
    end
end

function test_delete_file_not_exists(provider::CLIBucket.Provider)
    bucket = CLIBucketService(provider, BUCKET_NAME)
    # Verify that we throw an exception if the file doesn't exist
    @test_throws ArgumentError Bucket.delete(bucket, TEST_FILE_NAME)
end

function test_bucket(provider::CLIBucket.Provider)
    test_creatable(provider)
    test_bad_bucket(provider)

    test_download_IO(provider)
    test_download_file(provider)
    test_download_stream(provider)

    test_upload_io(provider)
    test_upload_file(provider)

    test_delete_file_exists(provider)
    test_delete_file_not_exists(provider)
end

function __init__()
    env = AWS.AWSEnv()
    aws_provider = AWSCLIProvider.Details(env)
    test_bucket(aws_provider)

    #=gcs_provider = GCSCLIProvider.Details()=#
    #=test_bucket(gcs_provider)=#
end

end # module TestCLIBucket

max_forks_repo_name: seung-lab/SimpleTasks


text:
# gen/common.jl
# https://developer.riotgames.com/docs/lol#data-dragon

data_version = "11.17.1"
locales = ("en_US", "ko_KR")
resources = (champion=:Champions, item=:Items, summoner=:SummonerSpells)


text:
# Autogenerated wrapper script for gh_cli_jll for armv7l-linux-gnueabihf
export gh

JLLWrappers.@generate_wrapper_header("gh_cli")
JLLWrappers.@declare_executable_product(gh)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        gh,
        "bin/gh",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
onfail(body, x) = error("I might have overlooked something: $x")
onfail(body, _::Test.Pass) = nothing
onfail(body, _::Tuple{Test.Fail,T}) where {T} = body()

function is_ci()
    get(ENV, "TRAVIS", "") == "true" ||
    get(ENV, "APPVEYOR", "") in ("true", "True") ||
    get(ENV, "CI", "") in ("true", "True")
end

absdiff(x, y) = x > y ? x - y : y - x

function imdiff(a, b)
    a_view = channelview(a)
    b_view = channelview(b)
    diffscale = max(maxabsfinite(a_view), maxabsfinite(b_view))
    d = sum(absdiff.(a_view, b_view))
    return d / (length(a) * diffscale)
end

_convert(C, T, xs::AbstractArray) =
    collect(colorview(C, map(i -> collect(reinterpret(T, collect(xs)[:, :, i])), 1:size(xs, 3))...))
_convert(C, T, xs::AbstractMatrix) = collect(colorview(C, collect(reinterpret(T, collect(xs)))))

_standardize_grayness(x) = x
_standardize_grayness(x::AbstractArray{<:Gray{Bool}}) = convert(Array{Gray{N0f8}}, x)
_standardize_grayness(x::AbstractArray{<:RGB}) = all(red.(x) .≈ green.(x) .≈ blue.(x)) ? Gray.(red.(x)) : x
_standardize_grayness(x::AbstractArray{<:RGBA}) = all(red.(x) .≈ green.(x) .≈ blue.(x)) ?
    convert(Array{GrayA}, colorview(GrayA, red.(x), alpha.(x))) :
    x

function plotdiffs(p, i)
    """
    Returns a Matrix of Colorants with following structure:
    --------------|---------------|----------------------|---------------------
    Orig Image p  | Orig Image i  | AbsDiff all channels | RelDiff all channels
    --------------|---------------|----------------------|---------------------
    1st channel p | 1st channel i | AbsDiff 1st channel  | RelDiff 1st channel
    --------------|---------------|----------------------|---------------------
    ...
    """
    _p = convert.(RGBA{Float64}, p)
    __p = collect(channelview(_p))
    __p[4, :, :] .= 0.0
    _i = convert.(RGBA{Float64}, i)
    __i = collect(channelview(_i))
    __i[4, :, :] .= 1.0 # make the difference in alpha be 1 so it is visible
    d = collect(colorview(RGBA{Float64}, absdiff.(__p, __i)))
    o = _plotdiffs(p, i)
    hborder = hcat(
        fill(RGBA{Float64}(0, 0, 0, 1), size(o, 1) + size(p, 1), 1),
        fill(RGBA{Float64}(1, 1, 1, 1), size(o, 1) + size(p, 1), 1),
        fill(RGBA{Float64}(0, 0, 0, 1), size(o, 1) + size(p, 1), 1),
    )

    collect(hcat(
        vcat(_p, _plotchannels(p)),
        hborder,
        vcat(_i, _plotchannels(i)),
        hborder,
        vcat(d, o),
        hborder,
        vcat(_rescale(d), _rescale(o)),
    ))
end

vborder(p) = vcat(
    fill(RGBA{Float64}(0, 0, 0, 1), 1, size(p, 2)),
    fill(RGBA{Float64}(1, 1, 1, 1), 1, size(p, 2)),
    fill(RGBA{Float64}(0, 0, 0, 1), 1, size(p, 2)),
)

function _plotdiffs(p, i)
    vcat(
        vborder(p),
        RGBA{Float64}.(absdiff.(red.(i), red.(p)), 0, 0, 1),
        vborder(p),
        RGBA{Float64}.(0, absdiff.(green.(i), green.(p)), 0, 1),
        vborder(p),
        RGBA{Float64}.(0, 0, absdiff.(blue.(i), blue.(p)), 1),
    )
end

function _plotdiffs(p::AbstractArray{<:RGBA}, i::AbstractArray{<:RGBA})
    a = absdiff.(alpha.(i), alpha.(p))
    vcat(
        vborder(p),
        RGBA{Float64}.(absdiff.(red.(i), red.(p)), 0, 0, 1),
        vborder(p),
        RGBA{Float64}.(0, absdiff.(green.(i), green.(p)), 0, 1),
        vborder(p),
        RGBA{Float64}.(0, 0, absdiff.(blue.(i), blue.(p)), 1),
        vborder(p),
        RGBA{Float64}.(a, a, a, 1),
    )
end

function _plotchannels(p)
    vcat(
        vborder(p),
        RGBA{Float64}.(red.(p), 0, 0, 1),
        vborder(p),
        RGBA{Float64}.(0, green.(p), 0, 1),
        vborder(p),
        RGBA{Float64}.(0, 0, blue.(p), 1),
    )
end

function _plotchannels(p::AbstractArray{<:RGBA})
    vcat(
        vborder(p),
        RGBA{Float64}.(red.(p), 0, 0, 1),
        vborder(p),
        RGBA{Float64}.(0, green.(p), 0, 1),
        vborder(p),
        RGBA{Float64}.(0, 0, blue.(p), 1),
        vborder(p),
        RGBA{Float64}.(alpha.(p), alpha.(p), alpha.(p), 1),
    )
end

function _rescale(x::AbstractArray)
    c = channelview(x)
    m = maximum(c)
    m = iszero(m) ? one(m) : m
    colorview(RGBA{Float64}, 100 * c ./ m)
end

## COPIED FROM Images.jl/src/algorithms.jl

"""
`m = maxabsfinite(A)` calculates the maximum absolute value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function maxabsfinite(A::AbstractArray{T}) where T
    ret = sentinel_min(typeof(abs(A[1])))
    for a in A
        ret = maxfinite_scalar(abs(a), ret)
    end
    ret
end

minfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b < a ? b : a) : b
maxfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b > a ? b : a) : b
minfinite_scalar(a::T, b::T) where {T<:Union{Integer,FixedPoint}} = b < a ? b : a
maxfinite_scalar(a::T, b::T) where {T<:Union{Integer,FixedPoint}} = b > a ? b : a
minfinite_scalar(a, b) = minfinite_scalar(promote(a, b)...)
maxfinite_scalar(a, b) = maxfinite_scalar(promote(a, b)...)

function minfinite_scalar(c1::C, c2::C) where C<:AbstractRGB
    C(minfinite_scalar(c1.r, c2.r),
      minfinite_scalar(c1.g, c2.g),
      minfinite_scalar(c1.b, c2.b))
end
function maxfinite_scalar(c1::C, c2::C) where C<:AbstractRGB
    C(maxfinite_scalar(c1.r, c2.r),
      maxfinite_scalar(c1.g, c2.g),
      maxfinite_scalar(c1.b, c2.b))
end

sentinel_min(::Type{T}) where {T<:Union{Integer,FixedPoint}} = typemax(T)
sentinel_max(::Type{T}) where {T<:Union{Integer,FixedPoint}} = typemin(T)
sentinel_min(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
sentinel_max(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
sentinel_min(::Type{C}) where {C<:AbstractRGB} = _sentinel_min(C, eltype(C))
_sentinel_min(::Type{C},::Type{T}) where {C<:AbstractRGB,T} = (s = sentinel_min(T); C(s,s,s))
sentinel_max(::Type{C}) where {C<:AbstractRGB} = _sentinel_max(C, eltype(C))
_sentinel_max(::Type{C},::Type{T}) where {C<:AbstractRGB,T} = (s = sentinel_max(T); C(s,s,s))
sentinel_min(::Type{C}) where {C<:AbstractGray} = _sentinel_min(C, eltype(C))
_sentinel_min(::Type{C},::Type{T}) where {C<:AbstractGray,T} = C(sentinel_min(T))
sentinel_max(::Type{C}) where {C<:AbstractGray} = _sentinel_max(C, eltype(C))
_sentinel_max(::Type{C},::Type{T}) where {C<:AbstractGray,T} = C(sentinel_max(T))


text:
"Molecular weight of air in kg/mol"
const Ma = 0.0289635

"Molecular weight of water in kg/mol"
const Mv = 0.01801528

"Universal gas constant kg m^2 /s^2 /K /mol"
const R = 8.314459848

"Melting point of water"
const T0 = 273.15

abstract type AbstractMoistAirModel <: ThermoState.Types.ThermoModel end
"""
Type used to model dry air

This uses the formulations presente in the ASHRAE handbook Psychrometrics: Theory and Practice, reference [5].
"""
struct ASHRAEDryAir <: AbstractMoistAirModel end

const DryAirModel = ASHRAEDryAir
molecular_weight(model::ASHRAEDryAir) = Ma
"""
Type used to model saturated water vapor

This uses the formulations present in the ASHRAE handbook Psychrometrics: Theory and Practice, reference [5].

"""
struct ASHRAEWaterVapor <:  AbstractMoistAirModel end

const WaterVaporModel = ASHRAEWaterVapor
molecular_weight(model::ASHRAEWaterVapor) = Mv

"""
Type used to model moist air

This uses the formulations presente in the ASHRAE handbook Psychrometrics: Theory and Practice, reference [5].

The temperature range is `173.15 < T < 474.15 K` with pressures `P < 5 MPa`.

It should be noted that the specific enthalpy, entropy and volume are calculated with respect to the mass of dry air. Therefore, 1 J/kg is actually 1 J per kg of dry air.
"""
struct ASHRAEMoistAir <:  AbstractMoistAirModel end
const MoistAirModel = ASHRAEMoistAir
molecular_weight(model::ASHRAEMoistAir) = SVector(1000*Ma,1000*Mv)

# package code goes here
include("utilities.jl")
include("docs.jl")
include("hyland83.jl")
include("hyland83a.jl")
include("wetbulb.jl")
include("moistair_impl.jl")


text:
using BitIntegers, Test
using Statistics: mean
using Serialization: serialize, deserialize

module TestBitIntegers

using BitIntegers, Test

BitIntegers.@define_integers 24
BitIntegers.@define_integers 200
BitIntegers.@define_integers 8  MyInt8 MyUInt8

# the following should throw, but is hard to test:
# BitIntegers.@define_integers 8 MyInt8

@testset "definitions" begin
    @test @isdefined Int24
    @test @isdefined UInt24
    @test sizeof(Int24) == sizeof(UInt24) == 3
    @test Int24  <: Signed
    @test UInt24 <: Unsigned

    @test @isdefined Int200
    @test @isdefined UInt200
    @test sizeof(Int200) == sizeof(UInt200) == 25
    @test Int200  <: Signed
    @test UInt200 <: Unsigned

    @test @isdefined MyInt8
    @test @isdefined MyUInt8
    @test sizeof(MyInt8) == sizeof(MyUInt8) == 1
    @test MyInt8  <: Signed
    @test MyUInt8 <: Unsigned
end

end # module TestBitIntegers

const BInts = Base.BitInteger_types
const XInts = (BitIntegers.BitInteger_types..., TestBitIntegers.UInt24, TestBitIntegers.Int24)
const Ints = (BInts..., XInts...)

# we don't include most Base-only type combinations:
const TypeCombos =
    [((X, Y) for X in Ints for Y in (X ∈ BInts ? XInts : Ints))...,
     (Int, Int), (Int, UInt), (UInt, Int), (UInt, UInt)]

nothing


text:
module Torch

using DocStringExtensions
using StructuredMechModels: AbstractRigidBody
using PyCall
using Zygote: Params

const torch = PyNULL()
const nn = PyNULL()
const optim = PyNULL()
const F = PyNULL()
const py_utils = PyNULL()
const structmechmod = PyNULL()

abstract type PyTorchRigidBody <: AbstractRigidBody end

"""
$(TYPEDEF)
"""
struct Naive{I<:Integer} <: PyTorchRigidBody
    nn::PyObject
    qdim::I
    udim::I
    thetamask
end

"""
Models the rigidbody as a Black-Box Neural Network

$(SIGNATURES)

**Arguments**:

- `qdim :: Int` Dimensions of position configuration
- `udim :: Int` Dimensions of active control
- `thetamask :: Vector{Int}` Binary Vector with ones corresponding to angles in the configuration space
"""
function Naive(qdim::Int, udim::Int, thetamask; activation="Tanh", hidden_sizes=(128, 128, 128))
    nn = structmechmod.rigidbody.NaiveRigidBody(qdim, udim, torch.from_numpy(thetamask), hidden_sizes)
    return Naive(nn, qdim, udim, thetamask)
end

struct DeLan{I<:Integer} <: PyTorchRigidBody
    nn::PyObject
    qdim::I
    udim::I
    thetamask
end

"""
Models the rigidbody as a Structured Mechanical Model with no constraints on the external forces.

$(SIGNATURES)

**Arguments**:

- `qdim`: Dimensions of position configuration
- `udim`: Dimensions of active control
- `thetamask`: Binary Vector with ones corresponding to angles in the configuration space

Useful when no prior knowledge about the active forces on a mechanical system.
"""
function DeLan(qdim, udim, thetamask; activation="Tanh", hidden_sizes=(32, 32, 32))
    @assert all(hidden_sizes .== hidden_sizes[1])
    nn = structmechmod.rigidbody.DeLan(qdim, hidden_sizes[1], length(hidden_sizes), torch.from_numpy(thetamask), udim=udim, activation=activation)
    return DeLan(nn, qdim, udim, thetamask)
end

struct ControlAffine{I<:Integer} <: PyTorchRigidBody
    nn::PyObject
    qdim::I
    udim::I
    thetamask
end

"""
Models the rigidbody as a Structured Mechanical Model with the constraints that active forces on the system are control affine.

$(SIGNATURES)

**Arguments**:

- `qdim`: Dimensions of position configuration
- `udim`: Dimensions of active control
- `thetamask`: Binary Vector with ones corresponding to angles in the configuration space

Almost all mechanical systems are control affine in the input forces for control.
"""
function ControlAffine(qdim, udim, thetamask; activation="Tanh", hidden_sizes=(32, 32, 32))
    @assert all(hidden_sizes .== hidden_sizes[1])
    forces = structmechmod.models.GeneralizedForces([structmechmod.models.ControlAffineForceNet(qdim, udim, hidden_sizes),
                                               structmechmod.models.QVForceNet(qdim, hidden_sizes)])
    nn = structmechmod.rigidbody.DeLan(qdim, hidden_sizes[1], length(hidden_sizes),
                                 torch.from_numpy(thetamask), udim=udim, activation=activation, forces=forces)
    return ControlAffine(nn, qdim, udim, thetamask)
end

"""

$(SIGNATURES)

**Arguments**:

- `x`: Input state (includes state and actions)
- `u`: Input control
"""
function (f::PyTorchRigidBody)(x, u)
    q = torch.from_numpy(permutedims(x[1:f.qdim, :], (2, 1)))
    v = torch.from_numpy(permutedims(x[f.qdim+1:2*f.qdim, :], (2, 1)))
    u = torch.from_numpy(permutedims(u, (2, 1)))

    qdot, vdot = f.nn(0, (q, v), u)
    qdot, vdot = map(x->permutedims(x, (2, 1)), (qdot.detach().numpy(), vdot.detach().numpy()))
    return vcat(qdot, vdot)
end

"""
Returns the parameters of `PyTorchRigidBody` as required by `Flux`.

$(SIGNATURES)
"""
function params(f::PyTorchRigidBody)
    ps = f.nn.get_params()
    return Params(ps)
end

"""

$(SIGNATURES)
"""
function train!(model::PyTorchRigidBody, train_data, valid_data, hparams)
    train_data = map(x->(permutedims(x, (2, 1))), (train_data.xs, train_data.us, train_data.xsp))
    valid_data = map(x->(permutedims(x, (2, 1))), (valid_data.xs, valid_data.us, valid_data.xsp))
    structmechmod.trainer.train(model.nn, train_data, valid_data, hparams)
end

"""
$(SIGNATURES)
"""
function load_state!(model::PyTorchRigidBody, path)
    model.nn.load_state_dict(torch.load(path))
    return model
end

function __init__()
    copy!(torch, pyimport("torch"))
    copy!(nn, pyimport("torch.nn"))
    copy!(optim, pyimport("torch.optim"))
    copy!(optim, pyimport("torch.nn.functional"))
    copy!(structmechmod, pyimport("structmechmod"))
end
end


text:
__precompile__()
module ResponseModels
using FilmDataStructures
export eVtoOmega, eps_InAs_struct, eps_InAsntype, eps_InAsptype, hStep, cstRsp, sicRsp, siDsc, dptDsc, siRsp, prmMSi
# Collection of material permittivity models for use in heatLayers code.
# Conversion factors.
# Microns to electron volts.
const muEv = 1.239
# Radial frequency to electron volts.
const hBEv = 6.582e-16
# Temperature to electron volts.
const blzK = 8.6173e-5
# Electron volts to Joules.
const evJ = 1.602176565e-19
const eps_0 = 8.8542*^(10,-12)
const hbar = 1.05457*^(10,-34) #m2kg/s

@inline function cstRsp(cst::ComplexF64, enr::Float64)::ComplexF64

	return cst
end
precompile(cstRsp, (ComplexF64, Float64))
## Heaviside function for testing
@inline function hStep(arg::Float64)::Float64

	if arg >= 0.0

		return 1.0
	else
		return 0.0
	end
end
precompile(hStep, (Float64, ))
### Silicon carbide.
"""

	sicRsp(enr::Float64)::ComplexF64

Simplified permittivity response model for silicon carbide.
Approximately valid for energies between 0.1 and 10 electron volts.
"""
# Reproduction of the model reported in J. Appl. Phys. 106, 044306 (2009).
@inline function sicRsp(enr::Float64)::ComplexF64

	# Energy values (in electron volts) for the longitudinal and transverse optical phonons,
	# and decay rate.
	enrL = 0.12014342321784009
	enrT = 0.09831711568536007
	enrD = 0.0005901528146448003

	return 6.7 * (1.0 + /(enrL^2.0 - enrT^2.0, enrT^2.0 - enr^2.0 - im * enrD * enr))
end

precompile(sicRsp, (Float64, ))
### Silicon
# Sample usage
# Boron p-type doping
# brnDpt = dptDsc(0.044, 1.0e18)
# Empty doping
# nllDpt = dptDsc(0.0, 0.0)
# Calculate silicon model parameters
# siModE = prmMSi(tmpLst[1], nllDpt, nllDpt)
# siRspE(enr) = siRsp(enr, siModE)
"""
Data structure holding computed model parameters for the permittivity response of silicon for
given dopant properties at a fixed temperature. Call prmMSi(tmp::Float64, dnr::dptDsc,
acp::dptDsc)::siDsc for construction.

# Arguments
.tmp: temperature at which parameters where calculated.
.dnrFac: strength of donor impurity Lorentzian.
.acpFac: strength of acceptor impurity Lorentzian.
.dnrDcy: decay rate (gamma) for donor impurity Lorentzian.
.acpDcy: decay rate (gamma) for acceptor impurity Lorentzian.
"""
struct siDsc

	tmp::Float64
	dnrFac::Float64
	acpFac::Float64
	dnrDcy::Float64
	acpDcy::Float64
end
# Include model code.
include("siMod.jl")
"""

	siRsp(enr::Float64, tmp::Float64, siVal::siDsc)::ComplexF64

Permittivity response model for doped silicon. See source file for additional details.
"""
@inline function siRsp(enr::Float64, siVal::siDsc)::ComplexF64

	return siRspBck(enr, siVal.tmp) - /(siVal.dnrFac, enr^2 + im * enr * siVal.dnrDcy) - /(siVal.acpFac, enr^2 + im * enr * siVal.acpDcy)
end
precompile(siRsp, (Float64, siDsc))

using eps_InAs_ntype
using eps_InAs_ptype
struct InAsDsc
    N0::Float64
    T::Float64
	mstar::Float64
	gamma_ntype::Float64
	gamma_ptype::Float64
	E0_T_InAs_value::Float64
end

@inline function eVtoOmega(energy) #energy in eV
    return(energy*evJ/hbar)
end
precompile(eVtoOmega, (Float64,))

@inline function eps_InAs_struct(N0,T)
	E0_T_InAs_value = E0_T_InAs(N0,T)
	gamma_ntype = Gamma_ntype(N0,T)
	gamma_ptype = Gamma_ptype(N0,T)
    mstar= m_star(N0)
    return InAsDsc(N0,T,mstar,gamma_ntype,gamma_ptype,E0_T_InAs_value)
end
precompile(eps_InAs_struct,(Float64,Float64))

@inline function eps_InAsntype(E_photon::Float64,InAsDsc_struct::InAsDsc)::ComplexF64 #N_base in m^(-3)
     omega = eVtoOmega(E_photon)
     #o_angular = 5.38*10^14   #hbar*5.38*10^14/eV = 0.35eV
    if(E_photon>0.35)
        return epsIB(omega,InAsDsc_struct.N0,InAsDsc_struct.T,InAsDsc_struct.E0_T_InAs_value) + epsFCL(omega,InAsDsc_struct.mstar,InAsDsc_struct.gamma_ntype,InAsDsc_struct.N0)
    else
        return epsFCL(omega,InAsDsc_struct.mstar,InAsDsc_struct.gamma_ntype,InAsDsc_struct.N0)
    end
end
precompile(eps_InAs_ntype,(Float64,InAsDsc))

@inline function eps_InAsptype(E_photon::Float64,InAsDsc_struct::InAsDsc)::ComplexF64 #N_base in m^(-3)
	omega = eVtoOmega(E_photon)
	#o_angular = 5.38*10^14   #hbar*5.38*10^14/eV = 0.35eV
   if(E_photon>0.35)
	   return epsIB(omega,InAsDsc_struct.N0,InAsDsc_struct.T,InAsDsc_struct.E0_T_InAs_value) + epsFCL(omega,InAsDsc_struct.mstar,InAsDsc_struct.gamma_ptype,InAsDsc_struct.N0)
   else
	   return epsFCL_ptype(omega,InAsDsc_struct.mstar,InAsDsc_struct.gamma_ptype,InAsDsc_struct.N0)
   end
end
precompile(eps_InAs_ptype,(Float64,InAsDsc))
end


text:
###################################################
"""
    getR_tile = getR_tile(tile::VisClaw.AMRGrid; length_unit="")

Get x and y ranges of a tile in AbstractString for -R option in GMT
"""
function getR_tile(tile::VisClaw.AMRGrid; length_unit="")
    xs = tile.xlow
    ys = tile.ylow
    xe = round(tile.xlow + tile.mx*tile.dx, digits=4)
    ye = round(tile.ylow + tile.my*tile.dy, digits=4)
    xyrange="$xs/$xe/$ys/$ye"
    isempty(length_unit) || (xyrange = xyrange*"+u"*length_unit)
    # return value
    return xyrange
end
###################################################
"""
    xyrange = getR(tiles::Vector{VisClaw.AMRGrid}; length_unit="")
    xyrange = getR(topo::VisClaw.AbstractTopo; length_unit="")
    xyrange = getR(region::VisClaw.Region; length_unit="")
    xyrange = getR(G::GMT.GMTgrid; length_unit="")

Get x and y ranges in AbstractString for -R option in GMT
"""
function getR(tiles::Vector{VisClaw.AMRGrid}; length_unit="")
    xs, xe, ys, ye = VisClaw.getlims(tiles)
    xyrange="$xs/$xe/$ys/$ye"
    isempty(length_unit) || (xyrange = xyrange*"+u"*length_unit)
    # return value
    return xyrange
end
###################################################
function getR(topo::VisClaw.AbstractTopo; length_unit="")
    xs=topo.x[1]
    xe=topo.x[end]
    ys=topo.y[1]
    ye=topo.y[end]
    xyrange="$xs/$xe/$ys/$ye"
    isempty(length_unit) || (xyrange = xyrange*"+u"*length_unit)
    return xyrange
end
###################################################
function getR(G::GMT.GMTgrid; length_unit="")
    x = extrema(G.x)
    y = extrema(G.y)
    xyrange="$(x[1])/$(x[2])/$(y[1])/$(y[2])"
    isempty(length_unit) || (xyrange = xyrange*"+u"*length_unit)
    return xyrange
end
###################################################
function getR(region::VisClaw.Region; length_unit="")
    xyrange="$(region.xlims[1])/$(region.xlims[2])/$(region.ylims[1])/$(region.ylims[2])"
    isempty(length_unit) || (xyrange = xyrange*"+u"*length_unit)
    return xyrange
end
###################################################

"""
    hwratio = axesratio(tiles::Vector{VisClaw.AMRGrid})
    hwratio = axesratio(topo::VisClaw.AbstractTopo)
    hwratio = axesratio(region::VisClaw.Region)
    hwratio = axesratio(G::GMT.GMTgrid)

Get height/width ratio
"""
function axesratio(tiles::Vector{VisClaw.AMRGrid})
    xs, xe, ys, ye = VisClaw.getlims(tiles)
    hwratio = (ye-ys)/(xe-xs)
    # return value
    return hwratio
end
###################################################
function axesratio(topo::VisClaw.AbstractTopo)
    xs=topo.x[1]
    xe=topo.x[end]
    ys=topo.y[1]
    ye=topo.y[end]
    hwratio = (ye-ys)/(xe-xs)
    # return value
    return hwratio
end
###################################################
function axesratio(G::GMT.GMTgrid)
    x = extrema(G.x)
    y = extrema(G.y)
    hwratio = (y[2]-y[1])/(x[2]-x[1])
    # return value
    return hwratio
end
###################################################
axesratio(region::VisClaw.Region) = (region.ylims[2]-region.ylims[1])/(region.xlims[2]-region.xlims[1])
###################################################

###################################################
"""
    G = geogrd(geo::VisClaw.Topo; kwargs...)
    G = geogrd(geo::VisClaw.DTopo, itime::Integer=0; kwargs...)

Generate grd (GMT) data
"""
function geogrd(geo::VisClaw.Topo; kwargs...)

    Δ = geo.dx
    R = VisClaw.getR(geo)
    xvec = repeat(geo.x, inner=(geo.nrows,1))
    yvec = repeat(geo.y, outer=(geo.ncols,1))

    G = GMT.surface([xvec[:] yvec[:] geo.elevation[:]]; R=R, I=Δ, kwargs...)

    return G
end
###################################################
function geogrd(geo::VisClaw.DTopo, itime::Integer=0; kwargs...)

    Δ = geo.dx
    R = VisClaw.getR(geo)
    xvec = repeat(geo.x, inner=(geo.my,1))
    yvec = repeat(geo.y, outer=(geo.mx,1))

    ( (itime < 0) || (geo.mt < itime) ) && error("Invalid time")
    if geo.mt == 1
        G = GMT.surface([xvec[:] yvec[:] geo.deform[:]]; R=R, I=Δ, kwargs...)
    elseif itime == 0
        G = GMT.surface([xvec[:] yvec[:] vec(geo.deform[:,:,end])]; R=R, I=Δ, kwargs...)
    else
        G = GMT.surface([xvec[:] yvec[:] vec(geo.deform[:,:,itime])]; R=R, I=Δ, kwargs...)
    end
    return G
end
###################################################

###################################################
"""
    proj = getJ(proj_base::AbstractString, hwratio::Real)

Correct J option
"""
function getJ(proj_base::AbstractString, hwratio)
    # find projection specifier
    J1 = match(r"^([a-zA-Z]+)", proj_base)
    J2 = match(r"([a-zA-Z]+).+?([a-zA-Z]+)", proj_base)

    J1 === nothing && error("Invald argument proj_base: $proj_base")

    # assign figure width
    # check whether variable proj_base contains any number
    regex = r"([+-]?(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?)"
    chkwidth = match(regex, proj_base)

    fwidth = chkwidth === nothing ? fwidth=10 : parse(Float64, chkwidth.captures[1])

    # assign figure height
    # check whether variable proj_base contains the height
    regex = r"([+-]?(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?).+?([+-]?(?:\d+\.?\d*|\.\d+)(?:[eE][+-]?\d+)?)"
    chkheight = match(regex, proj_base)

    fheight = chkheight === nothing ? hwratio*fwidth : parse(Float64, chkheight.captures[2])

    # generate J option
    if occursin("/",proj_base) && chkheight !== nothing
        proj = proj_base
    else
        proj = J2 === nothing ? J1.captures[1]*"$fwidth"*"/$fheight" : J1.captures[1]*"$fwidth"*J2.captures[2]*"/$fheight"*J2.captures[2]
    end
    # return value
    return proj
end
###################################################
getJ(proj_base::AbstractString, topo::VisClaw.Topo) = getJ(proj_base, VisClaw.axesratio(topo))
###################################################
getJ(proj_base::AbstractString, amr::Vector{VisClaw.AMRGrid}) = getJ(proj_base, VisClaw.axesratio(amr))
###################################################


text:
"""
    abstract type AbstractSingularityHandler{R}

Used for handling localized integrand singularities in `R`.
"""
abstract type AbstractSingularityHandler{T} end

"""
    struct IMT{A,P} <: AbstractSingularityHandler{ReferenceLine}

One-dimensional change of variables mapping `[0,1] -> [0,1]` with the property that
all derivatives vanish at the point `x=0`.

See [Davis and Rabinowitz](https://www.elsevier.com/books/methods-of-numerical-integration/davis/978-0-12-206360-2).
"""
struct IMT{A,P} <: AbstractSingularityHandler{ReferenceLine}
end
IMT(;a=1,p=1) = IMT{a,p}()

domain(::IMT) = ReferenceLine()
range(::IMT)  = ReferenceLine()

function (f::IMT{A,P})(x) where {A,P}
    exp(A*(1-1/x[1]^P))
end

derivative(f::IMT{A,P},x) where {A,P} = f(x) * A*P*1 / x[1]^(P+1)

jacobian(f::IMT,x) = derivative(f,x) |> SMatrix{1,1}

"""
    struct Kress{P} <: AbstractSingularityHandler{ReferenceLine}

Change of variables mapping `[0,1]` to `[0,1]` with the property that the first
`P` derivatives of the transformation vanish at `x=0`.
"""
struct Kress{P} <: AbstractSingularityHandler{ReferenceLine}
end
Kress(;order=5) = Kress{order}()

domain(k::Kress) = ReferenceLine()
range(k::Kress)  = ReferenceLine()

# NOTE: fastmath is needed here to allow for various algebraic simplifications
# which are not exact in floating arithmetic. Maybe reorder the operations *by
# hand* to avoid having to use fastmath? In any case, benchmark first.
function (f::Kress{P})(y) where {P}
    x = y[1]
    v = (x) -> (1/P - 1/2)*((1-x))^3 + 1/P*((x-1)) + 1/2
    return 2v(x)^P / (v(x)^P + v(2-x)^P)
end

function derivative(f::Kress{P},y) where {P}
    x = y[1]
    v = (x) -> (1/P - 1/2)*((1-x))^3 + 1/P*((x-1)) + 1/2
    vp = (x) -> -3*(1/P - 1/2)*((1-x))^2 + 1/P
    return 2 * (P*v(x)^(P-1)*vp(x) * (v(x)^P + v(2-x)^P) - (P*v(x)^(P-1)*vp(x) - P*v(2-x)^(P-1)*vp(2-x) ) * v(x)^P ) /
        (v(x)^P + v(2-x)^P)^2
end

jacobian(f::Kress,x) = derivative(f,x) |> SMatrix{1,1}

"""
    struct KressP{P} <: AbstractSingularityHandler{ReferenceLine}

Like [`Kress`](@ref), this change of variables maps the interval `[0,1]` onto
itself, but the first `P` derivatives of the transformation vanish at **both**
endpoints.

This change of variables can be used to *periodize* integrals in the following
sense. Suppose we wish to compute the integral of `f(x)` from `0` to `1` where
`f is not a `1-periodic function. If `ϕ` is an object of type `KressP`, then
using it as a change of variables in the integration yields a similar integral
from `0` to `1` (the interval `0≤0≤1` is mappend onto itself), but with
integrand given by `g(x) = f(ϕ(x))ϕ'(x)`. Since `ϕ'` vanishes (together with `P`
of its derivatives), the function `g(x)` is now periodic (up to derivatives of
order up to `P`) at the endpoints. Thus quadrature rules designed for periodic
functions like the [`TrapezoidalP`](@ref) can be used to obtain high order
convergence of `g`, which in turn yields a modified quadrature rule when viewed
as a quadrature rule for `f`.
"""
struct KressP{P} <: AbstractSingularityHandler{ReferenceLine}
end
KressP(;order=5) = KressP{order}()

domain(k::KressP) = ReferenceLine()
range(k::KressP)  = ReferenceLine()

@fastmath function (f::KressP{P})(y) where {P}
    x = y[1]
    v = (x) -> (1/P - 1/2)*((1-2x))^3 + 1/P*((2x-1)) + 1/2
    return v(x)^P / (v(x)^P + v(1-x)^P)
end

@fastmath function derivative(f::KressP{P},y) where {P}
    x = y[1]
    v =  (x) -> (1/P - 1/2)*((1-2x))^3 + 1/P*((2x-1)) + 1/2
    vp = (x) -> -6*(1/P - 1/2)*((1-2x))^2 + 2/P
    return (P*v(x)^(P-1)*vp(x) * (v(x)^P + v(1-x)^P) - (P*v(x)^(P-1)*vp(x) - P*v(1-x)^(P-1)*vp(1-x) ) * v(x)^P ) /
        (v(x)^P + v(1-x)^P)^2
end

jacobian(f::KressP,x) = derivative(f,x) |> SMatrix{1,1}

function (shand::KressP)(ent::ParametricCurve)
    f = parametrization(ent)
    new_f = v -> f(shand(v))
    return ParametricCurve(new_f)
end

function (shand::KressP)(ent::AbstractParametricBody)
    bnd = boundary(ent)
    map!(f->shand(f),bnd,bnd)
    return ent
end

"""
    struct Duffy <: AbstractSingularityHandler{RefereceTriangle}

Change of variables mapping the `ReferenceSquare` to the `RefereceTriangle` with
the property that the jacobian vanishes at the `(1,0)` vertex of the triangle.

Useful for integrating functions with a singularity on the `(1,0)` edge of the
reference triangle.
"""
struct Duffy <: AbstractSingularityHandler{ReferenceTriangle} end

domain(::Duffy) = ReferenceSquare()
range(::Duffy)  = ReferenceTriangle()

function (::Duffy)(u)
    SVector(u[1],(1-u[1])*u[2])
end

function jacobian(::Duffy,u)
    SMatrix{2,2,Float64}(1,0,-u[2][1],(1-u[1][1]))
end

# TODO: generalize to `N` dimensions
"""
    struct TensorProductSingularityHandler{S} <: AbstractSingularityHandler{ReferenceSquare}

A tensor product of two one-dimensional `AbstractSingularityHandler`s for performing
integration over the `ReferenceSquare`.
"""
struct TensorProductSingularityHandler{S} <: AbstractSingularityHandler{ReferenceSquare}
    shandler::S
end

domain(::TensorProductSingularityHandler) = ReferenceSquare()
range(::TensorProductSingularityHandler) = ReferenceSquare()

function TensorProductSingularityHandler(q...)
    TensorProductSingularityHandler(q)
end

TensorProductSingularityHandler{Tuple{P,Q}}() where {P,Q} = TensorProductSingularityHandler(P(),Q())

function (f::TensorProductSingularityHandler)(x)
    shandler = f.shandler
    @assert length(shandler) == length(x)
    N = length(shandler)
    svector(i->shandler[i](x[i]),N)
end

function jacobian(f::TensorProductSingularityHandler,x)
    shandler = f.shandler
    @assert length(shandler) == length(x)
    N = length(shandler)
    if N == 2
        jx = derivative(shandler[1],x[1])
        jy = derivative(shandler[2],x[2])
        return SMatrix{2,2}(jx,0,0,jy)
    else
        notimplemented()
    end
end

# # overload quadgk so that a singularity handler can be passed as first argument
# function quadgk(phi::SingularityHandler{ReferenceLine},f,args...;xs,kwargs...)
#     a,b = args[1], args[end]
#     quadgk(args...;kwargs...) do x
#         xi,dxi = phi(x,a,b,xs), derivative(phi,x,a,b,xs)
#         f(xi)*dxi
#     end
# end

# """
#     Window{A,B,S}

# Change of variables mapping `[0,1]` to `[0,1]` with the following properties:
# - smooth (infinitely differentiable) on `[0,1]`
# - all derivatives vanish on both `x=0` and `x=1`
# - derivative of the transoformation is exactly 1 between `A` and `B`
# """
# struct Window{A,B,S} <: AbstractSingularityHandler{ReferenceLine}
# end
# Window() = Window{1,1,5}()

# domain(::Window) = ReferenceLine()
# range(::Window)  = ReferenceLine()

# normalization(w::Window) = quadgk(t->_derivative(w,t),0,1,rtol=1e-16)[1]

# function (f::Window{A,B,S})(x) where {A,B,S}
#     x == 0 && return 0.0
#     I,_ = quadgk(t->derivative(f,t),0,x)
#     return I
#     # FIXME: the current way of computing the values of the `Window` change of
#     # variables is very innefficient. Could use a e.g. `cumsum`.
# end

# function derivative(f::Window,x)
#     x==0 && (return 0.0)
#     x==1 && (return 0.0)
#     _derivative(f,x) / normalization(f)
# end

# function _derivative(f::Window{A,B,S},x) where {A,B,S}
#     if  A ≤ x ≤ B
#         return 1
#     elseif 0 ≤ x < A
#         u = (A - x) / (A)
#         return exp(S*exp(-1/u)/(u-1))
#     elseif B < x ≤ 1
#         u = (x - B) / (1 - B)
#         return exp(S*exp(-1/u)/(u-1))
#     else
#         return 0
#     end
# end

# jacobian(f::Window,x) = derivative(f,x) |> SMatrix{1,1}

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
"""
CumulativeResourceDeadlineType

Whether resources in `CumulativeResource` constraint have deadlines:

* either there are no deadlines: `NO_DEADLINE_CUMULATIVE_RESOURCE`
* or deadlines are given as variables: `VARIABLE_DEADLINE_CUMULATIVE_RESOURCE`
"""
@enum CumulativeResourceDeadlineType begin
    NO_DEADLINE_CUMULATIVE_RESOURCE
    VARIABLE_DEADLINE_CUMULATIVE_RESOURCE
end

"""
    CumulativeResource{CRDT}(n_tasks::Int)

This set models most variants of task scheduling with cumulative resource
usage. Presence of deadlines can be indicated with the
`CumulativeResourceDeadlineType` enumeration.

## Without deadline

Each task is given by a minimum start time (the first `n_tasks` variables),
a duration (the next `n_tasks` variables), and the resource consumption
(the following `n_tasks` variables). The final variable is the maximum amount of
the resource available.

Also called [`cumulative`](https://sofdem.github.io/gccat/gccat/Ccumulative.html).
This version does not consider end deadlines for tasks.

## With variable deadline

Each task is given by a minimum start time (the first `n_tasks` variables),
a duration (the next `n_tasks` variables), a deadline (the following `n_tasks`
variables), and the resource consumption (the next `n_tasks` variables).
The final variable is the maximum amount of the resource available.

Also called [`cumulative`](https://sofdem.github.io/gccat/gccat/Ccumulative.html)
"""
struct CumulativeResource{CRDT} <: MOI.AbstractVectorSet
    n_tasks::Int
end

function CumulativeResource(n_tasks::Int)
    return CumulativeResource{NO_DEADLINE_CUMULATIVE_RESOURCE}(n_tasks)
end

MOI.dimension(set::CumulativeResource{NO_DEADLINE_CUMULATIVE_RESOURCE}) = 3 * set.n_tasks + 1
MOI.dimension(set::CumulativeResource{VARIABLE_DEADLINE_CUMULATIVE_RESOURCE}) = 4 * set.n_tasks + 1

"""
    NonOverlappingOrthotopesConditionalityType

Whether orthotopes in `NonOverlappingOrthotopes` constraint are considered:

* either all orthotopes must be considered: `UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES`
* or orthotopes can be disabled by variables: `CONDITIONAL_NONVERLAPPING_ORTHOTOPES`
"""
@enum NonOverlappingOrthotopesConditionalityType begin
    UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES
    CONDITIONAL_NONVERLAPPING_ORTHOTOPES
end

"""
    NonOverlappingOrthotopes{NOOCT}(n_orthotopes::Int, n_dimensions::Int)

This set corresponds to a guarantee that orthotopes do not overlap. Some
orthotopes can optionally be disabled for the constraint (guided by variables),
based on the value of `NonOverlappingOrthotopesConditionalityType`.

## Unconditional constraint

Guarantees that the `n_orthotopes` orthotopes do not overlap. The orthotopes
live in various dimensions: segments if `n_dimensions = 1`, rectangles if
`n_dimensions = 2`, rectangular parallelepiped if `n_dimensions = 3`,
hyperrectangles otherwise.

The variables are packed by orthotope:

* the first `n_dimensions` are the origin of the orthotope
* the next `n_dimensions` are the size of the orthotope in each dimension
* the last `n_dimensions` are the destination of the orthotope. These variables
  are automatically constrained to be `origin + size` (unlike other modelling
  layers, such as Gecode)

The set can be defined as:

``(o_1, s_1, d_1, o_2, s_2, d_2 \\dots o_\\mathtt{o}, s_\\mathtt{o}, d_\\mathtt{o}) \\in \\mathbb{R}^{3 \\times \\mathtt{o} \\times \\mathtt{d} }``

Also called [`diffn`](https://sofdem.github.io/gccat/gccat/Cdiffn.html),
`geost`, `nooverlap`, `diff2`, or `disjoint`.

### Example: two 2-D rectangles
    [x1, y1, w1, h1, x1e, y1e, x2, y2, w2, h2, x2e, y2e] in NonOverlappingOrthotopes(2, 2)
    # Enforces the following five constraints:
    #   OR(
    #     x1 + w1 <= x2,
    #     x2 + w2 <= x1,
    #     y1 + h1 <= y2,
    #     y2 + h2 <= y1
    #   )
    #   x1e = x1 + w1
    #   y1e = y1 + h1
    #   x2e = x2 + w2
    #   y2e = y2 + h2

## Conditional constraint

Guarantees that the `n_orthotopes` orthotopes do not overlap, with a binary
variable indicating whether a given orthotope must not overlap with other
orthotopes (if 1) or if it can be ignored (if 0). The orthotopes live in
various dimensions: segments if `n_dimensions = 1`, rectangles if
`n_dimensions = 2`, rectangular parallelepiped if `n_dimensions = 3`,
hyperrectangles otherwise.

The variables are packed by orthotope:

* the first `n_dimensions` are the origin of the orthotope
* the next `n_dimensions` are the size of the orthotope in each dimension
* the next `n_dimensions` are the destination of the orthotope. These variables
  are automatically constrained to be `origin + size` (unlike other modelling
  layers, such as Gecode)
* the last variable indicates whether the orthotope is mandatory (`true`) or
  optional (`false`)

The set can be defined as:

``(o_1, s_1, d_1, m_1, o_2, s_2, d_2, m_2 \\dots o_\\mathtt{o}, s_\\mathtt{o}, d_\\mathtt{o}, m_\\mathtt{o}) \\in \\prod_{i=1}^{\\mathtt{o}} (\\mathbb{R}^{3 \\times \\mathtt{d} \\times \\{0, 1\\}) }``

Also called [`diffn`](https://sofdem.github.io/gccat/gccat/Cdiffn.html),
`nooverlap`, or `disjointconditional`.
"""
struct NonOverlappingOrthotopes{NOOCT} <: MOI.AbstractVectorSet
    n_orthotopes::Int
    n_dimensions::Int
end

function NonOverlappingOrthotopes(n_orthotopes::Int, n_dimensions::Int)
    return NonOverlappingOrthotopes{UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES}(n_orthotopes, n_dimensions)
end

MOI.dimension(set::NonOverlappingOrthotopes{UNCONDITIONAL_NONVERLAPPING_ORTHOTOPES}) = 3 * set.n_orthotopes * set.n_dimensions
MOI.dimension(set::NonOverlappingOrthotopes{CONDITIONAL_NONVERLAPPING_ORTHOTOPES}) = 3 * set.n_orthotopes * set.n_dimensions + set.n_orthotopes

# isbits types, nothing to copy
function copy(
    set::Union{
        CumulativeResource,
        NonOverlappingOrthotopes,
    },
)
    return set
end


text:
using Plots, LaTeXStrings, Pkg, Revise
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

#pkg"activate ."
using FourierGPE

# ==== Units: ========================
# this example works in oscillator units

# ==== set simulation parameters ====
L = (20.0,20.0)
N = (128,128)
μ = 25.0

# ========= Initialize simulation ======
sim = Sim(L,N)
@pack! sim = μ
@unpack_Sim sim

# declare the potential function
import FourierGPE.V
V(x,y,t) = 0.5*(x^2 + y^2)

# useful TF state
using VortexDistributions
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

#make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
makevortex!(ψi,[0. 0. 1],x,y)
showpsi(x,y,ψi)
ϕi = kspace(ψi,sim)
@pack! sim = ϕi
sim
@unpack T = sim

T = Transforms(T,psi = ArrayPartition(ϕi,copy(ϕi),copy(ϕi)))

sim = Sim(sim,T=T)
ϕj = copy(ϕi)

function Lz(ϕ,sim)
    @unpack T,X,K = sim
    @unpack psi = T
    xky,ykx,L = psi.x[1],psi.x[2],psi.x[3]
    kx,ky = K
    x,y = X
    xky .= ϕ.*ky'
    ykx .= kx.*ϕ
    T.Tkx!*xky
    T.Tkx!*ykx
    @. L = x*xky - ykx*y'
    T.Txk!*L
    return L
end

@time Ltest = Lz(ϕj,sim)
Ltest == ϕj
ψjL = xspace(Ltest,sim)
sum(conj(ψi).*ψjL)/sum(abs2.(ψi))

# Ok, we have allocation free evaluation!!!

# ====== Evolve in k space ==========
sol = runsim(sim)
# ===================================

# pull out the ground state:
ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)

# Add a vortex off-axis
using VortexDistributions

# initial state: imprint a vortex inside Thomas-Fermi radius
healing(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
Rtf = sqrt(2*μ)
rv = 0.5*Rtf
xv,yv,cv = rv, 0.0, 1
initialvortex = [xv yv cv]
ξv = healing(xv,yv,μ,g)
ψv = copy(ψg)
makeallvortices!(ψv,initialvortex,x,y,ξv)
showpsi(x,y,ψv)

# In TF regim precession frequency is given analytically by:
# (see Fetter JLTP 2010)
ξ = 1/sqrt(μ)
Ωm = 3*log(Rtf/ξ/sqrt(2))/2/Rtf^2
Ωv = Ωm/(1-rv^2/Rtf^2)

#or a period of
Tv = 2*π/Ωv

# ==== set simulation parameters ====
γ = 0.0
tf = Tv
t = LinRange(ti,tf,Nt)
ϕi = kspace(ψv,sim)
reltol = 1e-7
alg = DP5()
@pack! sim = tf,t,γ,ϕi,reltol,alg
# initsim!(sim)

@unpack_Sim sim

# ====== Evolve in k space ==========
solv = runsim(sim)
# ===================================

ϕf = solv[100]
ψf = xspace(ϕf,sim)
showpsi(x,y,ψf)

#trim last few frames to show one orbit
# analytical result is within 10%
anim = @animate for i=1:Nt-6
    ψ = xspace(solv[i],sim)
    showpsi(x,y,ψ)
end

gif(anim,"./examples/vortex.gif",fps=30)


text:
### LiS module ###
### Linnette Teo 2022.01.07 ###
### Thermodynamic model to solve for equilibrium solutions for lithium sulfur system ###
### Three defined regions based on differing solids present ###

function thermoModel1!(F,x)
    """
    Equations for region 1: only S8 solid
    """

    global dod
    # variables to solve
    CLi = exp(x[1])
    CS8l = exp(x[2])
    CS8m = exp(x[3])
    CS6m = exp(x[4])
    CS4m = exp(x[5])
    CS2m = exp(x[6])
    CSm2 = exp(x[7])
    phi1 = exp(x[8])
    epcS8 = exp(x[9])

    # porosites
#     epcS8 = 0
    epcLi2S = 0
    epsS8 = 0
    epsLi2S = 0
    ec = 1-epcS8-epcLi2S-efc
    es = 1-epsS8-epsLi2S-efs

    # potential equations, setting overpotentials=0
    F[1] = phi1 - U2 + f *(sS8m_2*log(CS8m/1000) + sS8l_2*log(CS8l/1000) )
    F[2] = phi1 - U3 + f *(sS6m_3*log(CS6m/1000) + sS8m_3*log(CS8m/1000) )
    F[3] = phi1 - U4 + f *(sS4m_4*log(CS4m/1000) + sS6m_4*log(CS6m/1000) )
    F[4] = phi1 - U5 + f *(sS2m_5*log(CS2m/1000) + sS4m_5*log(CS4m/1000) )
    F[5] = phi1 - U6 + f *(sSm2_6*log(CSm2/1000) + sS2m_6*log(CS2m/1000) )

    # charge balance
    F[6] = (ec*Lcat + es*Lsep) * 2*(CS8m + CS6m + CS4m + CS2m + CSm2) +
            epcLi2S/VLi2S*Lcat*2 + epsLi2S/VLi2S*Lsep*2 - S_bal*2*dod/100 - inichg

    # electroneutrality
    F[7] = CAref*(1-efc-e_S8init)/ec + 2*(CS8m + CS6m + CS4m + CS2m + CSm2) - CLi

    # S mass balance
    F[8] = (ec*Lcat + es*Lsep) * (CS8l*8 + CS8m*8 + CS6m*6 + CS4m*4 + CS2m*2 + CSm2*1) +
            epcS8/VS8*Lcat*8 + epsS8/VS8*Lsep*8 + epcLi2S/VLi2S*Lcat*1 + epsLi2S/VLi2S*Lsep*1 -S_bal

    # Precipitation
    F[9] = CS8l - KspS8

end

function thermoModel2!(F,x)
    """
    Equations for region 2: no solids
    """

    global dod
    # variables to solve
    CLi = exp(x[1])
    CS8l = exp(x[2])
    CS8m = exp(x[3])
    CS6m = exp(x[4])
    CS4m = exp(x[5])
    CS2m = exp(x[6])
    CSm2 = exp(x[7])
    phi1 = exp(x[8])

    # porosites
    epcS8 = 0
    epcLi2S = 0
    epsS8 = 0
    epsLi2S = 0
    ec = 1-epcS8-epcLi2S-efc
    es = 1-epsS8-epsLi2S-efs

    # potential equations, setting overpotentials=0
    F[1] = phi1 - U2 + f *(sS8m_2*log(CS8m/1000) + sS8l_2*log(CS8l/1000) )
    F[2] = phi1 - U3 + f *(sS6m_3*log(CS6m/1000) + sS8m_3*log(CS8m/1000) )
    F[3] = phi1 - U4 + f *(sS4m_4*log(CS4m/1000) + sS6m_4*log(CS6m/1000) )
    F[4] = phi1 - U5 + f *(sS2m_5*log(CS2m/1000) + sS4m_5*log(CS4m/1000) )
    F[5] = phi1 - U6 + f *(sSm2_6*log(CSm2/1000) + sS2m_6*log(CS2m/1000) )

    # charge balance
    F[6] = (ec*Lcat + es*Lsep) * 2*(CS8m + CS6m + CS4m + CS2m + CSm2) +
            epcLi2S/VLi2S*Lcat*2 + epsLi2S/VLi2S*Lsep*2 - S_bal*2*dod/100 - inichg

    # electroneutrality
    F[7] = CAref*(1-efc-e_S8init)/ec + 2*(CS8m + CS6m + CS4m + CS2m + CSm2) - CLi

    # S mass balance
    F[8] = (ec*Lcat + es*Lsep) * (CS8l*8 + CS8m*8 + CS6m*6 + CS4m*4 + CS2m*2 + CSm2*1) +
            epcS8/VS8*Lcat*8 + epsS8/VS8*Lsep*8 + epcLi2S/VLi2S*Lcat*1 + epsLi2S/VLi2S*Lsep*1 - S_bal

end

function thermoModel3!(F,x)
    """
    Equations for region 3: only Li2S solid
    """

    global dod
    # variables to solve
    CLi = exp(x[1])
    CS8l = exp(x[2])
    CS8m = exp(x[3])
    CS6m = exp(x[4])
    CS4m = exp(x[5])
    CS2m = exp(x[6])
    CSm2 = exp(x[7])
    phi1 = exp(x[8])
    epcLi2S = exp(x[9])

    # porosites
    epcS8 = 0
#     epcLi2S = 0
    epsS8 = 0
    epsLi2S = 0
    ec = 1-epcS8-epcLi2S-efc
    es = 1-epsS8-epsLi2S-efs

    # potential equations, setting overpotentials=0
    F[1] = phi1 - U2 + f *(sS8m_2*log(CS8m/1000) + sS8l_2*log(CS8l/1000) )
    F[2] = phi1 - U3 + f *(sS6m_3*log(CS6m/1000) + sS8m_3*log(CS8m/1000) )
    F[3] = phi1 - U4 + f *(sS4m_4*log(CS4m/1000) + sS6m_4*log(CS6m/1000) )
    F[4] = phi1 - U5 + f *(sS2m_5*log(CS2m/1000) + sS4m_5*log(CS4m/1000) )
    F[5] = phi1 - U6 + f *(sSm2_6*log(CSm2/1000) + sS2m_6*log(CS2m/1000) )

    # charge balance
    F[6] = (ec*Lcat + es*Lsep) * 2*(CS8m + CS6m + CS4m + CS2m + CSm2) +
            epcLi2S/VLi2S*Lcat*2 + epsLi2S/VLi2S*Lsep*2 - S_bal*2*dod/100 - inichg

    # electroneutrality
    F[7] = CAref*(1-efc-e_S8init)/ec + 2*(CS8m + CS6m + CS4m + CS2m + CSm2) - CLi

    # S mass balance
    F[8] = (ec*Lcat + es*Lsep) * (CS8l*8 + CS8m*8 + CS6m*6 + CS4m*4 + CS2m*2 + CSm2*1) +
            epcS8/VS8*Lcat*8 + epsS8/VS8*Lsep*8 + epcLi2S/VLi2S*Lcat*1 + epsLi2S/VLi2S*Lsep*1 - S_bal

    # Precipitation
    F[9] = CLi^2 * CSm2 - KspLi2S

end

function findIdx(dod,step)
    """Finds index of dod in range 1:step:99"""
    if step < 1
        idx = ceil(Int, dod/step) -1
    else
        idx = ceil(Int, dod/step)
    end

    return idx
end

function runThermoRegion!(dod_arr, x0, thermoModel, sol_mat, tol, step)
    """
    Solves thermoModel for a specific region over a range dod_arr.
    dod_arr has step size of Int step.
    Starts from initial guess x0.
    Store voltage solutions in place in V_arr.
    Store all solutions in place in sol_mat (initialized as 10 columns)
    Solver used is nlsolve, newton method.
    Breaks if solution does not converge OR
    if we hit a tolerance limit (i.e. vol fracs < tol)
    """

    for i in 1:length(dod_arr)
        global dod = dod_arr[i]

        sol = nlsolve(thermoModel, x0, autodiff = :forward, method=:newton,xtol=1e-10,ftol=1e-10)

        # update initial guess to previous DOD solution
        global x0 = sol.zero

        if converged(sol)
            if ((thermoModel == thermoModel1!) && (exp.(x0[9])>tol)) ||
                thermoModel == thermoModel2! ||
                ((thermoModel == thermoModel3!) && (exp.(x0[9])>tol)) ||

                ((thermoModel == thermoModel1_S3!) && (exp.(x0[9])>tol)) ||
                    thermoModel == thermoModel2_S3! ||
                    ((thermoModel == thermoModel3_S3!) && (exp.(x0[9])>tol))

                # store solutions in sol_mat
                dod_ind  = findIdx(dod,step)   # index of dod
                store_sol!(thermoModel, sol_mat, sol, dod_ind)
#                 print(dod, ", ") #'\n')
            else
#                 print("\n hits tol=", exp.(x0[9]), "\n")
                break
            end
        else
#             print("\n not converged for dod=", dod, "\n")
            break
        end
    end
end

function store_sol!(thermoModel, sol_mat, sol, dod_ind)
    """
    Stores solutions into sol_mat depending on which model is used.
    Stored at the index Int dod_ind.
    For indexing purposes as different models have variables in different indices.
    Variable order: CLi, CS8l, CS8m, CS6m, CS4m, CS2m, CSm2, phi1, epcS8, epcLi2S, CS3
    """
    ## models without S3 concentration as a variable
    if size(sol_mat)[2]==10
        if thermoModel == thermoModel1!
            sol_mat[dod_ind,:] = append!(exp.(sol.zero),0)
        elseif thermoModel == thermoModel2!
            sol_mat[dod_ind,:] = append!(exp.(sol.zero),[0,0])
        elseif thermoModel == thermoModel3!
            sol_mat[dod_ind,:] = insert!(exp.(sol.zero),9,0)
        else
            print("\n solution array is > Nvar (number of variables)")
        end
    end

    ## with S3 radical concentration in index 11
    if size(sol_mat)[2]==11
        if thermoModel == thermoModel1_S3!
            sol_mat[dod_ind,:] = insert!(exp.(sol.zero),10,0)
        elseif thermoModel == thermoModel2_S3!
            sol_mat[dod_ind,:] = insert!(insert!(exp.(sol.zero),9,0),9,0)
        elseif thermoModel == thermoModel3_S3!
            sol_mat[dod_ind,:] = insert!(exp.(sol.zero),9,0)
        else
            print("\n solution array is > Nvar (number of variables)")
        end
    end

end

function solveThermoAll!(sol_mat, x0_99, x0_1; tol=1e-3, step=1)

    """
    Algorithm:
    Solve Region 3 in descending order, starting from 99% DOD till no more solution found.
    Solve Region 1 in ascending order, starting from 1% DOD till no more solution found.
    Solve Region 2 in ascending order, starting from where Region 1 ended till where Region 3 begun.

    Stores all solutions in sol_mat, a len(dod_arr) x Nvar=10 matrix
    x0_99: initial guess for 99% dod,
    x0_1: initial guess for 1% dod

    also return dod_r2LB and dod_r2UB which are the lower and upper bounds of region 2, inclusive

    Default arguments:
    tol: tolerance limit for solver for vol fracs (see runThermoRegion!)                             step: step size for dod_arr

    """
    #### RUNNING REGION 3 #####
    #print("\nREGION 3\n")
    global x0 = x0_99
    dod_arr = 99:-step:1
    thermoModel = thermoModel3!
    runThermoRegion!(dod_arr, x0, thermoModel, sol_mat, tol, step)
    dod_r2UB = dod   # dod where region 3 begins and region 2 ends

    #### RUNNING REGION 1 #####
    #print("\nREGION 1\n")
    global x0 = x0_1
    dod_arr = 1:step:99
    thermoModel = thermoModel1!
    runThermoRegion!(dod_arr, x0, thermoModel, sol_mat, tol, step)
    dod_r2LB = dod

    ### RUNNING REGION 2 #####
    #print("\nREGION 2\n")
    dod_arr = dod_r2LB:step:dod_r2UB
    if dod_r2LB==1
        global x0 = x0_1[1:8]
    else         # look at element before  the LB for initial guess
        idx  = findIdx(dod,step)
        global x0 = log.(sol_mat[idx-1,1:8])
    end
    thermoModel = thermoModel2!
    runThermoRegion!(dod_arr, x0, thermoModel, sol_mat, tol, step)

    return dod_r2LB, dod_r2UB
end

######################### S3 RADICAL MODEL ############################

function thermoModel1_S3!(F,x)
    """
    Equations for region 1: only S8 solid
    With S3 radical chemistry
    """

    global dod
    # variables to solve
    CLi = exp(x[1])
    CS8l = exp(x[2])
    CS8m = exp(x[3])
    CS6m = exp(x[4])
    CS4m = exp(x[5])
    CS2m = exp(x[6])
    CSm2 = exp(x[7])
    phi1 = exp(x[8])
    epcS8 = exp(x[9])
    CS3 = exp(x[10])

    # porosites
#     epcS8 = 0
    epcLi2S = 0
    epsS8 = 0
    epsLi2S = 0
    ec = 1-epcS8-epcLi2S-efc
    es = 1-epsS8-epsLi2S-efs

    # potential equations, setting overpotentials=0
    F[1] = phi1 - U2 + f *(sS8m_2*log(CS8m/1000) + sS8l_2*log(CS8l/1000) )
    F[2] = phi1 - U3 + f *(sS6m_3*log(CS6m/1000) + sS8m_3*log(CS8m/1000) )
    F[3] = phi1 - U4 + f *(sS4m_4*log(CS4m/1000) + sS6m_4*log(CS6m/1000) )
    F[4] = phi1 - U5 + f *(sS2m_5*log(CS2m/1000) + sS4m_5*log(CS4m/1000) )
    F[5] = phi1 - U6 + f *(sSm2_6*log(CSm2/1000) + sS2m_6*log(CS2m/1000) )

    # charge balance
    F[6] = (ec*Lcat + es*Lsep) * 2*(CS8m + CS6m + CS4m + CS2m + CSm2 + CS3/2) +
            epcLi2S/VLi2S*Lcat*2 + epsLi2S/VLi2S*Lsep*2 - S_bal*2*dod/100 - inichg

    # electroneutrality
    F[7] = CAref*(1-efc-e_S8init)/ec + 2*(CS8m + CS6m + CS4m + CS2m + CSm2) + CS3 - CLi

    # S mass balance
    F[8] = (ec*Lcat + es*Lsep) * (CS8l*8 + CS8m*8 + CS6m*6 + CS4m*4 + CS2m*2 + CSm2*1 + CS3*3) + epcS8/VS8*Lcat*8 + epsS8/VS8*Lsep*8 + epcLi2S/VLi2S*Lcat*1 + epsLi2S/VLi2S*Lsep*1 -S_bal

    # Precipitation
    F[9] = CS8l - KspS8

    # Radical chemistry
    F[10] = CS3^2/CS6m - KS3
end

function thermoModel2_S3!(F,x)
    """
    Equations for region 2: no solids
    With S3 radical chemistry
    """

    global dod
    # variables to solve
    CLi = exp(x[1])
    CS8l = exp(x[2])
    CS8m = exp(x[3])
    CS6m = exp(x[4])
    CS4m = exp(x[5])
    CS2m = exp(x[6])
    CSm2 = exp(x[7])
    phi1 = exp(x[8])
    CS3 = exp(x[9])

    # porosites
    epcS8 = 0
    epcLi2S = 0
    epsS8 = 0
    epsLi2S = 0
    ec = 1-epcS8-epcLi2S-efc
    es = 1-epsS8-epsLi2S-efs

    # potential equations, setting overpotentials=0
    F[1] = phi1 - U2 + f *(sS8m_2*log(CS8m/1000) + sS8l_2*log(CS8l/1000) )
    F[2] = phi1 - U3 + f *(sS6m_3*log(CS6m/1000) + sS8m_3*log(CS8m/1000) )
    F[3] = phi1 - U4 + f *(sS4m_4*log(CS4m/1000) + sS6m_4*log(CS6m/1000) )
    F[4] = phi1 - U5 + f *(sS2m_5*log(CS2m/1000) + sS4m_5*log(CS4m/1000) )
    F[5] = phi1 - U6 + f *(sSm2_6*log(CSm2/1000) + sS2m_6*log(CS2m/1000) )

    # charge balance
    F[6] = (ec*Lcat + es*Lsep) * 2*(CS8m + CS6m + CS4m + CS2m + CSm2 + CS3/2) +
            epcLi2S/VLi2S*Lcat*2 + epsLi2S/VLi2S*Lsep*2 - S_bal*2*dod/100 - inichg

    # electroneutrality
    F[7] = CAref*(1-efc-e_S8init)/ec + 2*(CS8m + CS6m + CS4m + CS2m + CSm2) + CS3 - CLi

    # S mass balance
    F[8] = (ec*Lcat + es*Lsep) * (CS8l*8 + CS8m*8 + CS6m*6 + CS4m*4 + CS2m*2 + CSm2*1 + CS3*3) +  epcS8/VS8*Lcat*8 + epsS8/VS8*Lsep*8 + epcLi2S/VLi2S*Lcat*1 + epsLi2S/VLi2S*Lsep*1 - S_bal

    # Radical chemistry
    F[9] = CS3^2/CS6m - KS3

end

function thermoModel3_S3!(F,x)
    """
    Equations for region 3: only Li2S solid
    With S3 radical chemistry
    """

    global dod
    # variables to solve
    CLi = exp(x[1])
    CS8l = exp(x[2])
    CS8m = exp(x[3])
    CS6m = exp(x[4])
    CS4m = exp(x[5])
    CS2m = exp(x[6])
    CSm2 = exp(x[7])
    phi1 = exp(x[8])
    epcLi2S = exp(x[9])
    CS3 = exp(x[10])

    # porosites
    epcS8 = 0
#     epcLi2S = 0
    epsS8 = 0
    epsLi2S = 0
    ec = 1-epcS8-epcLi2S-efc
    es = 1-epsS8-epsLi2S-efs

    # potential equations, setting overpotentials=0
    F[1] = phi1 - U2 + f *(sS8m_2*log(CS8m/1000) + sS8l_2*log(CS8l/1000) )
    F[2] = phi1 - U3 + f *(sS6m_3*log(CS6m/1000) + sS8m_3*log(CS8m/1000) )
    F[3] = phi1 - U4 + f *(sS4m_4*log(CS4m/1000) + sS6m_4*log(CS6m/1000) )
    F[4] = phi1 - U5 + f *(sS2m_5*log(CS2m/1000) + sS4m_5*log(CS4m/1000) )
    F[5] = phi1 - U6 + f *(sSm2_6*log(CSm2/1000) + sS2m_6*log(CS2m/1000) )

    # charge balance
    F[6] = (ec*Lcat + es*Lsep) * 2*(CS8m + CS6m + CS4m + CS2m + CSm2 + CS3/2) +
            epcLi2S/VLi2S*Lcat*2 + epsLi2S/VLi2S*Lsep*2 - S_bal*2*dod/100 - inichg

    # electroneutrality
    F[7] = CAref*(1-efc-e_S8init)/ec + 2*(CS8m + CS6m + CS4m + CS2m + CSm2) + CS3 - CLi

    # S mass balance
    F[8] = (ec*Lcat + es*Lsep) * (CS8l*8 + CS8m*8 + CS6m*6 + CS4m*4 + CS2m*2 + CSm2*1 + CS3*3) + epcS8/VS8*Lcat*8 + epsS8/VS8*Lsep*8 + epcLi2S/VLi2S*Lcat*1 + epsLi2S/VLi2S*Lsep*1 - S_bal

    # Precipitation
    F[9] = CLi^2 * CSm2 - KspLi2S

    # Radical chemistry
    F[10] = CS3^2/CS6m - KS3

end

function solveThermoAll_S3!(sol_mat, x0_99, x0_1; tol=1e-3, step=1)

    """
    Algorithm:
    Solve Region 3 in descending order, starting from 99% DOD till no more solution found.
    Solve Region 1 in ascending order, starting from 1% DOD till no more solution found.                   Solve Region 2 in ascending order, starting from where Region 1 ended till where Region 3 begun.

    Stores all solutions in sol_mat, a len(dod_arr) x Nvar=11 matrix
    x0_99: initial guess for 99% dod,
    x0_1: initial guess for 1% dod

    also return dod_r2LB and dod_r2UB which are the lower and upper bounds of region 2, inclusive
    """

    #### RUNNING REGION 3 #####
    #print("\nREGION 3\n")
    global x0 = copy(x0_99)
    append!(x0,CS3_0)
    dod_arr = 99:-step:1
    thermoModel = thermoModel3_S3!
    runThermoRegion!(dod_arr, x0, thermoModel,sol_mat, tol,step)
    dod_r2UB = dod   # dod where region 3 begins and region 2 ends

    # #### RUNNING REGION 1 #####
    #print("\nREGION 1\n")
    global x0 = copy(x0_1)
    append!(x0,CS3_0)
    dod_arr = 1:step:99
    thermoModel = thermoModel1_S3!
    runThermoRegion!(dod_arr, x0, thermoModel, sol_mat, tol, step)
    dod_r2LB = dod

    # ### RUNNING REGION 2 #####
    #print("\nREGION 2\n")
    dod_arr = dod_r2LB:step:dod_r2UB
    if dod_r2LB==1
        global x0 = x0_1[1:8]
    else         # look at element before  the LB for initial guess
        idx  = findIdx(dod,step)
        global x0 = log.(sol_mat[idx-1,1:8])
    end
    append!(x0,CS3_0)
    thermoModel = thermoModel2_S3!
    runThermoRegion!(dod_arr, x0, thermoModel, sol_mat, tol, step)

    return dod_r2LB, dod_r2UB
end


text:
using Test
using TestSetExtensions

using JSON
using LinearAlgebra: SingularException

using Nemo

include("../src/ExactODEReduction.jl")
using .ExactODEReduction: from_dense, to_dense, from_COO, zero_sparsik, random_sparsik,
                          construct_jacobians, check_invariance!, @sparse, linear_span!,
                          check_inclusion!, density, is_thorough, transpose!, tr, dim, scale,
                          to_cartesian, to_plain, load_COO_set, load_MTX, find_radical,
                          rational_reconstruction, modular_reduction, Sparsik, inner,
                          scale!, reduce!, first_nonzero, find_basis_1_β, find_basis_1,
                          minimal_polynomial_wiedemann, unit_sparsik, evaluate, reconstruct!,
                          square_nonsingular_deterministic_wiedemann, square_nonsingular_randomized_wiedemann

@info "Testing started"

@testset "All the tests" begin

    @includetests ARGS

end

@info "All tests OK"


text:
function Base.show(io::IO, mime::MIME, plt::Plot2D)
    show(io, mime, Plots.plot(plt))
    return
end

Base.showable(::MIME, ::Plot2D) = false
Base.showable(::MIME"text/plain", ::Plot2D) = true
Base.showable(::MIME"image/png", ::Plot2D) = true
Base.showable(::MIME"image/jpeg", ::Plot2D) = true