text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

# This example uses the social_network dataset. You can create the schema and upload data
# by following instructionns from https://docs.grakn.ai/docs/general/quickstart

using Revise

using TypeDBClient
import TypeDBClient.typedb.protocol as P

client = Client("127.0.0.1")
session = Session(client, "social_network", P.Session_Type.DATA)

# ------------------------------------------------------------------
# general query

function query(q::AbstractString)
    transaction = Transaction(
        session, session.sessionID, P.Transaction_Type.READ, TypeDBOptions()
    )
    result = TypeDBClient.match(transaction, q)
    close(transaction)
    return result
end

# Getting attributes
raw_result = query(raw"""
    match $p isa person, has full-name $n, has gender $g; get $n, $g;
    """)
result = map(r -> (fullname = r["n"].value, gender = r["g"].value),
    raw_result)

# Getting entities
raw_result = query(raw"""match $p isa person; get $p;""")
result = map(r -> (iid = r["p"].iid,), raw_result)

# Getting relations
raw_result = query(raw"""match $emp (employer: $x, employee: $y) isa employment; get $emp;""")

#=
typedb console returns the related iid's also!

Julia:
 ConceptMap:
    1. emp => TypeDBClient.Relation("aa8280058000000000000006",
        TypeDBClient.RelationType(:employment, false))

Console:
{ $emp iid 0xaa828005800000000000000a (
    employment:employer: iid 0x966e800a8000000000000003,
    employment:job: iid 0x966e800e8000000000000007,
    employment:employee: iid 0x966e8009800000000000000e) isa employment;
}
=#

# ------------------------------------------------------------------
# ConceptManager methods

function cm_call(f::Function, session::Session)
    transaction = Transaction(session, session.sessionID, P.Transaction_Type.READ, TypeDBOptions())
    concept_manager = ConceptManager(transaction)
    res = f(concept_manager)
    close(transaction)
    return res
end

cm_call(session) do concept_manager
    @show get(concept_manager, ThingType, "title")
    @show get(concept_manager, AttributeType, "title")  # same as ThingType
    @show get(concept_manager, EntityType, "title") # nothing
    @show get(concept_manager, RelationType, "title") # nothing
    nothing
end

cm_call(session) do concept_manager
    @show get(concept_manager, Entity, "966e800c8000000000000000")
    nothing
end

# ------------------------------------------------------------------
# ThingType methods

# Fetch a thing_type for testing
function thing_call(f::Function, session::Session, name::AbstractString)
    cm_call(session) do concept_manager
        thing_type_res = TypeDBClient.execute(
            concept_manager, TypeDBClient.ConceptManagerRequestBuilder.get_thing_type_req(name)
        )
        thing_type = TypeDBClient.instantiate(thing_type_res.get_thing_type_res.thing_type)
        remote_thing = TypeDBClient.as_remote(thing_type, concept_manager.transaction)
        f(remote_thing)
    end
end

res = thing_call(thing -> get_supertype(thing), session, "title")

res = thing_call(thing -> get_supertypes(thing), session, "approved-date")

res = thing_call(thing -> get_subtypes(thing), session, "event-date")

res = thing_call(thing -> get_instances(thing), session, "school") # entities
res = thing_call(thing -> get_instances(thing), session, "title") # attributes
res = thing_call(thing -> get_instances(thing), session, "ownership") # relations

res = thing_call(thing -> get_owns(thing), session, "person")
res = thing_call(thing -> get_owns(thing, VALUE_TYPE.BOOLEAN), session, "person")
res = thing_call(thing -> get_owns(thing, nothing, true), session, "person")

res = thing_call(thing -> get_plays(thing), session, "location")

# ------------------------------------------------------------------
# AttributeType methods

res = thing_call(thing -> get_regex(thing), session, "emotion")
res = thing_call(thing -> get_regex(thing), session, "title") # nothing

# TODO [INT03] Invalid Internal State: Illegal internal operation! This method should not have been called.
res = thing_call(thing -> set_regex(thing, raw".*"), session, "title")

# TODO [TYW03] Invalid Type Write: The type 'media' has instances, and cannot be set abstract.
# Should probably create my own attribute and set it to abstract
res = thing_call(thing -> set_abstract(thing), session, "media")

# TODO [INT03] Invalid Internal State: Illegal internal operation! This method should not have been called.
# Should probably create my own attribute and set it to abstract and then unset it
res = thing_call(thing -> unset_abstract(thing), session, "media")

# Isn't it a little dumb to get back what I passed?
res = thing_call(thing -> get(thing, "Biochemistry"), session, "title")
res = thing_call(thing -> get(thing, "Foo"), session, "title") # nothing

# TODO - what does it mean when calling Put method on an AttributeType
# https://discord.com/channels/762167454973296644/773324706924593192/851000661583986698
res = thing_call(thing -> TypeDBClient.put(thing, "Foo"), session, "title")

res = thing_call(thing -> get_subtypes(thing), session, "title")

# TODO - need an example to test the special get_subtypes filtering logic
# https://discord.com/channels/665254494820368395/837321963814912020/851164383714476074
# res = thing_call(thing -> get_subtypes(thing), session, "attribute")

close(session)
close(client)


text:
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_band.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_config.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_dense.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_direct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_iterative.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_nvector.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __N_VClone(w::N_Vector)
    ccall((:N_VClone,libsundials_sundials),N_Vector,(N_Vector,),w)
end

function N_VClone(w)
    __w = convert(NVector,w)
    __N_VClone(convert(N_Vector,__w))
end

function __N_VCloneEmpty(w::N_Vector)
    ccall((:N_VCloneEmpty,libsundials_sundials),N_Vector,(N_Vector,),w)
end

function N_VCloneEmpty(w)
    __w = convert(NVector,w)
    __N_VCloneEmpty(convert(N_Vector,__w))
end

function __N_VDestroy(v::N_Vector)
    ccall((:N_VDestroy,libsundials_sundials),Void,(N_Vector,),v)
end

function N_VDestroy(v)
    __v = convert(NVector,v)
    __N_VDestroy(convert(N_Vector,__v))
end

function __N_VSpace(v::N_Vector,lrw::Ptr{Clong},liw::Ptr{Clong})
    ccall((:N_VSpace,libsundials_sundials),Void,(N_Vector,Ptr{Clong},Ptr{Clong}),v,lrw,liw)
end

function N_VSpace(v,lrw,liw)
    __v = convert(NVector,v)
    __N_VSpace(convert(N_Vector,__v),pointer(lrw),pointer(liw))
end

function __N_VGetArrayPointer(v::N_Vector)
    ccall((:N_VGetArrayPointer,libsundials_sundials),Ptr{realtype},(N_Vector,),v)
end

function N_VGetArrayPointer(v)
    __v = convert(NVector,v)
    __N_VGetArrayPointer(convert(N_Vector,__v))
end

function __N_VSetArrayPointer(v_data::Ptr{realtype},v::N_Vector)
    ccall((:N_VSetArrayPointer,libsundials_sundials),Void,(Ptr{realtype},N_Vector),v_data,v)
end

function N_VSetArrayPointer(v_data,v)
    __v = convert(NVector,v)
    __N_VSetArrayPointer(pointer(v_data),convert(N_Vector,__v))
end

function __N_VLinearSum(a::realtype,x::N_Vector,b::realtype,y::N_Vector,z::N_Vector)
    ccall((:N_VLinearSum,libsundials_sundials),Void,(realtype,N_Vector,realtype,N_Vector,N_Vector),a,x,b,y,z)
end

function N_VLinearSum(a,x,b,y,z)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __z = convert(NVector,z)
    __N_VLinearSum(a,convert(N_Vector,__x),b,convert(N_Vector,__y),convert(N_Vector,__z))
end

function __N_VConst(c::realtype,z::N_Vector)
    ccall((:N_VConst,libsundials_sundials),Void,(realtype,N_Vector),c,z)
end

function N_VConst(c,z)
    __z = convert(NVector,z)
    __N_VConst(c,convert(N_Vector,__z))
end

function __N_VProd(x::N_Vector,y::N_Vector,z::N_Vector)
    ccall((:N_VProd,libsundials_sundials),Void,(N_Vector,N_Vector,N_Vector),x,y,z)
end

function N_VProd(x,y,z)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __z = convert(NVector,z)
    __N_VProd(convert(N_Vector,__x),convert(N_Vector,__y),convert(N_Vector,__z))
end

function __N_VDiv(x::N_Vector,y::N_Vector,z::N_Vector)
    ccall((:N_VDiv,libsundials_sundials),Void,(N_Vector,N_Vector,N_Vector),x,y,z)
end

function N_VDiv(x,y,z)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __z = convert(NVector,z)
    __N_VDiv(convert(N_Vector,__x),convert(N_Vector,__y),convert(N_Vector,__z))
end

function __N_VScale(c::realtype,x::N_Vector,z::N_Vector)
    ccall((:N_VScale,libsundials_sundials),Void,(realtype,N_Vector,N_Vector),c,x,z)
end

function N_VScale(c,x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VScale(c,convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VAbs(x::N_Vector,z::N_Vector)
    ccall((:N_VAbs,libsundials_sundials),Void,(N_Vector,N_Vector),x,z)
end

function N_VAbs(x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VAbs(convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VInv(x::N_Vector,z::N_Vector)
    ccall((:N_VInv,libsundials_sundials),Void,(N_Vector,N_Vector),x,z)
end

function N_VInv(x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VInv(convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VAddConst(x::N_Vector,b::realtype,z::N_Vector)
    ccall((:N_VAddConst,libsundials_sundials),Void,(N_Vector,realtype,N_Vector),x,b,z)
end

function N_VAddConst(x,b,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VAddConst(convert(N_Vector,__x),b,convert(N_Vector,__z))
end

function __N_VDotProd(x::N_Vector,y::N_Vector)
    ccall((:N_VDotProd,libsundials_sundials),realtype,(N_Vector,N_Vector),x,y)
end

function N_VDotProd(x,y)
    __x = convert(NVector,x)
    __y = convert(NVector,y)
    __N_VDotProd(convert(N_Vector,__x),convert(N_Vector,__y))
end

function __N_VMaxNorm(x::N_Vector)
    ccall((:N_VMaxNorm,libsundials_sundials),realtype,(N_Vector,),x)
end

function N_VMaxNorm(x)
    __x = convert(NVector,x)
    __N_VMaxNorm(convert(N_Vector,__x))
end

function __N_VWrmsNorm(x::N_Vector,w::N_Vector)
    ccall((:N_VWrmsNorm,libsundials_sundials),realtype,(N_Vector,N_Vector),x,w)
end

function N_VWrmsNorm(x,w)
    __x = convert(NVector,x)
    __w = convert(NVector,w)
    __N_VWrmsNorm(convert(N_Vector,__x),convert(N_Vector,__w))
end

function __N_VWrmsNormMask(x::N_Vector,w::N_Vector,id::N_Vector)
    ccall((:N_VWrmsNormMask,libsundials_sundials),realtype,(N_Vector,N_Vector,N_Vector),x,w,id)
end

function N_VWrmsNormMask(x,w,id)
    __x = convert(NVector,x)
    __w = convert(NVector,w)
    __id = convert(NVector,id)
    __N_VWrmsNormMask(convert(N_Vector,__x),convert(N_Vector,__w),convert(N_Vector,__id))
end

function __N_VMin(x::N_Vector)
    ccall((:N_VMin,libsundials_sundials),realtype,(N_Vector,),x)
end

function N_VMin(x)
    __x = convert(NVector,x)
    __N_VMin(convert(N_Vector,__x))
end

function __N_VWL2Norm(x::N_Vector,w::N_Vector)
    ccall((:N_VWL2Norm,libsundials_sundials),realtype,(N_Vector,N_Vector),x,w)
end

function N_VWL2Norm(x,w)
    __x = convert(NVector,x)
    __w = convert(NVector,w)
    __N_VWL2Norm(convert(N_Vector,__x),convert(N_Vector,__w))
end

function __N_VL1Norm(x::N_Vector)
    ccall((:N_VL1Norm,libsundials_sundials),realtype,(N_Vector,),x)
end

function N_VL1Norm(x)
    __x = convert(NVector,x)
    __N_VL1Norm(convert(N_Vector,__x))
end

function __N_VCompare(c::realtype,x::N_Vector,z::N_Vector)
    ccall((:N_VCompare,libsundials_sundials),Void,(realtype,N_Vector,N_Vector),c,x,z)
end

function N_VCompare(c,x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VCompare(c,convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VInvTest(x::N_Vector,z::N_Vector)
    ccall((:N_VInvTest,libsundials_sundials),Cint,(N_Vector,N_Vector),x,z)
end

function N_VInvTest(x,z)
    __x = convert(NVector,x)
    __z = convert(NVector,z)
    __N_VInvTest(convert(N_Vector,__x),convert(N_Vector,__z))
end

function __N_VConstrMask(c::N_Vector,x::N_Vector,m::N_Vector)
    ccall((:N_VConstrMask,libsundials_sundials),Cint,(N_Vector,N_Vector,N_Vector),c,x,m)
end

function N_VConstrMask(c,x,m)
    __c = convert(NVector,c)
    __x = convert(NVector,x)
    __m = convert(NVector,m)
    __N_VConstrMask(convert(N_Vector,__c),convert(N_Vector,__x),convert(N_Vector,__m))
end

function __N_VMinQuotient(num::N_Vector,denom::N_Vector)
    ccall((:N_VMinQuotient,libsundials_sundials),realtype,(N_Vector,N_Vector),num,denom)
end

function N_VMinQuotient(num,denom)
    __num = convert(NVector,num)
    __denom = convert(NVector,denom)
    __N_VMinQuotient(convert(N_Vector,__num),convert(N_Vector,__denom))
end

function __N_VCloneEmptyVectorArray(count::Cint,w::N_Vector)
    ccall((:N_VCloneEmptyVectorArray,libsundials_sundials),Ptr{N_Vector},(Cint,N_Vector),count,w)
end

function N_VCloneEmptyVectorArray(count,w)
    __w = convert(NVector,w)
    __N_VCloneEmptyVectorArray(convert(Cint,count),convert(N_Vector,__w))
end

function __N_VCloneVectorArray(count::Cint,w::N_Vector)
    ccall((:N_VCloneVectorArray,libsundials_sundials),Ptr{N_Vector},(Cint,N_Vector),count,w)
end

function N_VCloneVectorArray(count,w)
    __w = convert(NVector,w)
    __N_VCloneVectorArray(convert(Cint,count),convert(N_Vector,__w))
end

function __N_VDestroyVectorArray(vs::Ptr{N_Vector},count::Cint)
    ccall((:N_VDestroyVectorArray,libsundials_sundials),Void,(Ptr{N_Vector},Cint),vs,count)
end

function N_VDestroyVectorArray(vs,count)
    __N_VDestroyVectorArray(pointer(vs),convert(Cint,count))
end
# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_spbcgs.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_spgmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_sptfqmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_types.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_fnvector.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/astukalov/.julia/v0.5/Sundials/deps/usr/include/sundials/sundials_math.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

ext: jl


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

Create a REoptInputs for the Business As Usual scenario.
"""
function BAUInputs(p::REoptInputs)
    bau_scenario = BAUScenario(p.s)
    techs = Techs(p, bau_scenario)

    boiler_efficiency = Dict{String, Float64}()

    # REoptInputs indexed on techs.all:
    max_sizes = Dict(t => 0.0 for t in techs.all)
    min_sizes = Dict(t => 0.0 for t in techs.all)
    existing_sizes = Dict(t => 0.0 for t in techs.all)
    cap_cost_slope = Dict{String, Any}()
    om_cost_per_kw = Dict(t => 0.0 for t in techs.all)
    production_factor = DenseAxisArray{Float64}(undef, techs.all, p.time_steps)

    # export related inputs
    techs_by_exportbin = Dict(k => [] for k in p.s.electric_tariff.export_bins)
    export_bins_by_tech = Dict{String, Array{Symbol, 1}}()

    # REoptInputs indexed on techs.segmented
    n_segs_by_tech = Dict{String, Int}()
    seg_min_size = Dict{String, Any}()
    seg_max_size = Dict{String, Any}()
    seg_yint = Dict{String, Any}()

    # PV specific arrays
    pv_to_location = Dict(t => Dict{Symbol, Int}() for t in techs.pv)

    levelization_factor = Dict(t => 1.0 for t in techs.all)

    for pvname in techs.pv  # copy the optimal scenario inputs for existing PV systems
        production_factor[pvname, :] = p.production_factor[pvname, :]
        pv_to_location[pvname] = p.pv_to_location[pvname]
        existing_sizes[pvname] = p.existing_sizes[pvname]
        min_sizes[pvname] = p.existing_sizes[pvname]
        max_sizes[pvname] = p.existing_sizes[pvname]
        om_cost_per_kw[pvname] = p.om_cost_per_kw[pvname]
        levelization_factor[pvname] = p.levelization_factor[pvname]
        cap_cost_slope[pvname] = 0.0
        if pvname in p.techs.pbi
            push!(pbi_techs, pvname)
        end
        pv = get_pv_by_name(pvname, p.s.pvs)
        fillin_techs_by_exportbin(techs_by_exportbin, pv, pv.name)
        if !pv.can_curtail
            push!(techs.no_curtail, pv.name)
        end
    end

    if "Generator" in techs.all
        max_sizes["Generator"] = p.s.generator.existing_kw
        min_sizes["Generator"] = p.s.generator.existing_kw
        existing_sizes["Generator"] = p.s.generator.existing_kw
        cap_cost_slope["Generator"] = 0.0
        om_cost_per_kw["Generator"] = p.s.generator.om_cost_per_kw
        production_factor["Generator", :] = p.production_factor["Generator", :]
        fillin_techs_by_exportbin(techs_by_exportbin, p.s.generator, "Generator")
        if "Generator" in p.techs.pbi
            push!(pbi_techs, "Generator")
        end
        if !p.s.generator.can_curtail
            push!(techs.no_curtail, "Generator")
        end
    end

    if "ExistingBoiler" in techs.all
        setup_existing_boiler_inputs(bau_scenario, max_sizes, min_sizes, existing_sizes, cap_cost_slope, boiler_efficiency)
    end

    # filling export_bins_by_tech MUST be done after techs_by_exportbin has been filled in
    for t in techs.elec
        export_bins_by_tech[t] = [bin for (bin, ts) in techs_by_exportbin if t in ts]
    end

    t0, tf = p.s.electric_utility.outage_start_time_step, p.s.electric_utility.outage_end_time_step
    if tf > t0 && t0 > 0
        original_crit_lds = copy(p.s.electric_load.critical_loads_kw)
        update_bau_outage_outputs(bau_scenario, original_crit_lds, t0, tf, production_factor)

        if bau_scenario.outage_outputs.bau_critical_load_met_time_steps > 0
        # include critical load in bau load for the time that it can be met
            bau_scenario.electric_load.critical_loads_kw[
                t0 : t0 + bau_scenario.outage_outputs.bau_critical_load_met_time_steps
                ] = original_crit_lds[t0 : t0 + bau_scenario.outage_outputs.bau_critical_load_met_time_steps]
        end
    end

    REoptInputs(
        bau_scenario,
        techs,
        min_sizes,
        max_sizes,
        existing_sizes,
        cap_cost_slope,
        om_cost_per_kw,
        p.time_steps,
        p.time_steps_with_grid,
        p.time_steps_without_grid,
        p.hours_per_timestep,
        p.months,
        production_factor,
        levelization_factor,
        p.value_of_lost_load_per_kwh,
        p.pwf_e,
        p.pwf_om,
        p.pwf_fuel,
        p.third_party_factor,
        p.pvlocations,
        p.maxsize_pv_locations,
        pv_to_location,
        p.ratchets,
        techs_by_exportbin,
        export_bins_by_tech,
        n_segs_by_tech,
        seg_min_size,
        seg_max_size,
        seg_yint,
        p.pbi_pwf,
        p.pbi_max_benefit,
        p.pbi_max_kw,
        p.pbi_benefit_per_kwh,
        boiler_efficiency
    )
end

"""
    update_bau_outage_outputs(s::BAUScenario, crit_load, t0, tf, production_factors)

Update the `bau_critical_load_met` and `bau_critical_load_met_time_steps`  values.
"""
function update_bau_outage_outputs(s::BAUScenario, crit_load, t0, tf, production_factors)

    pv_kw_series = Float64[]  # actual output (not normalized)

    if any(pv.existing_kw > 0 for pv in s.pvs)  # fill in pv_kw_series
        for pv in s.pvs
            if pv.existing_kw > 0
                if length(pv_kw_series) == 0  # first non-zero existing_kw
                    pv_kw_series = pv.existing_kw * production_factors[pv.name, t0:tf]
                else
                    pv_kw_series += pv.existing_kw * production_factors[pv.name, t0:tf]
                end
            end
        end
    end

    s.outage_outputs.bau_critical_load_met, s.outage_outputs.bau_critical_load_met_time_steps =
        bau_outage_check(crit_load[t0:tf], pv_kw_series, s.generator, s.settings.time_steps_per_hour)
    nothing
end

"""
    bau_outage_check(critical_loads_kw::AbstractArray, pv_kw_series::AbstractArray, gen::Generator,
        time_steps_per_hour::Int)

Determine if existing generator and/or PV can meet critical load and for how long.

return: (Bool, Int) boolean for if the entire critical load is met and Int for number of time steps the existing
    generator and PV can meet the critical load
"""
function bau_outage_check(critical_loads_kw::AbstractArray, pv_kw_series::AbstractArray, gen::Generator,
    time_steps_per_hour::Int)

    fuel_gal = copy(gen.fuel_avail_gal)
    if gen.existing_kw == 0 && length(pv_kw_series) == 0
        return false, 0
    end

    if gen.existing_kw > 0
        if length(pv_kw_series) == 0
            pv_kw_series = zeros(length(critical_loads_kw))
        end

        for (i, (load, pv)) in enumerate(zip(critical_loads_kw, pv_kw_series))
            unmet = load - pv
            if unmet > 0
                fuel_kwh = (fuel_gal - gen.fuel_intercept_gal_per_hr) / gen.fuel_slope_gal_per_kwh
                gen_avail = minimum([fuel_kwh, gen.existing_kw * (1.0 / time_steps_per_hour)])
                gen_output = maximum([minimum([unmet, gen_avail]), gen.min_turn_down_pct * gen.existing_kw])
                fuel_needed = gen.fuel_intercept_gal_per_hr + gen.fuel_slope_gal_per_kwh * gen_output
                fuel_gal -= fuel_needed

                if gen_output < unmet
                    return false, i-1
                end
            end
        end

    else  # gen.existing_kw = 0 and pv.existing_kw > 0
        for (i, (load, pv)) in enumerate(zip(critical_loads_kw, pv_kw_series))
            unmet = load - pv
            if unmet > 0
                return false, i-1
            end
        end
    end

    return true, length(critical_loads_kw)
end


text:
using MPIParticleInCell
using Test
using MPI

@testset "MPIParticleInCell.jl" begin
    MPI.Init()

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    num_guard_cells = 1
    g = UniformGrid((0., 0.), (1., 1.), 2 .* (8, 5), num_guard_cells)
    rank == 0 && println(g)

    f = Field(g, "board", Int, true)
    rank == 0 && println(f)

    f.field_values .= rank
    gather_field_on_root!(f)
    rank == 0 && (display(f.field_values); println("\n"))

    f.field_values .= rank
    communicate_halo!(f, comm)
    rank == 0 && (display(f.field_values); println("\n"))

    MPI.Finalize()

    @test true
end


text:
using Poly
using Test
using StaticArrays
using BenchmarkTools

function poly_mul(A, B, out)
    n = size(out, 1)
    m = size(out, 2)
    r = size(A, 2)
    @poly_loop for i = 1:n
        for j = 1:m
            for k = 1:r
                out[i, j] += A[i, k] * B[k, j]
            end
        end
    end
end

function poly_mul_no_tile(A, B, out)
    n = size(out, 1)
    m = size(out, 2)
    r = size(A, 2)
    @poly_loop tile=0 for i = 1:n
        for j = 1:m
            for k = 1:r
                out[i, j] += A[i, k] * B[k, j]
            end
        end
    end
end

function mul(A, B, out)
    n = size(out, 1)
    m = size(out, 2)
    r = size(A, 2)
    @inbounds for i = 1:n
        for j = 1:m
            for k = 1:r
                out[i, j] += A[i, k] * B[k, j]
            end
        end
    end
end

function mul_right_order(A, B, out)
    n = size(out, 1)
    m = size(out, 2)
    r = size(A, 2)
    @inbounds for j = 1:m
        for k = 1:r
            for i = 1:n
                out[i, j] += A[i, k] * B[k, j]
            end
        end
    end
end

function tiled_mul(A, B, C)
    N = size(C)[1]
    TILE_DIM = 64
    NUM_TILES = Int(N/TILE_DIM)
    # tile1 = @MArray zeros(TILE_DIM, TILE_DIM)
    # tile2 = @MArray zeros(TILE_DIM, TILE_DIM)

    @inbounds for gj = 0:TILE_DIM:N-1
        for gk = 0:TILE_DIM:N-1
            for gi = 0:TILE_DIM:N-1
                # for j = 1:TILE_DIM
                #     for i = 1:TILE_DIM
                #         tile1[i, j] = A[gi + i, gk + j]
                #         tile2[i, j] = B[gk + i, gj + j]
                #     end
                # end
                for jj in 1:TILE_DIM
                    for k = 1:TILE_DIM
                        for ii = 1:TILE_DIM
                            # C[gi + ii, gj + jj] += tile1[ii, k] * tile2[k, jj]
                            C[gi + ii, gj + jj] += A[gi + ii, gk + k] * B[gk + k, gj + jj]
                        end
                    end
                end
            end
        end
    end
end

function tiled_mul_notsquare(A, B, C)
    N = size(C, 1)
    M = size(C, 2)
    R = size(A, 2)
    TILE_DIM = 64

    @inbounds for gj = 0:TILE_DIM:M-1
        for gk = 0:TILE_DIM:R-1
            for gi = 0:TILE_DIM:N-1
                for jj in 1:TILE_DIM
                    for k = 1:TILE_DIM
                        for ii = 1:TILE_DIM
                            C[gi + ii, gj + jj] += A[gi + ii, gk + k] * B[gk + k, gj + jj]
                        end
                    end
                end
            end
        end
    end
end

function tiled_mul_any(A, B, C)
    N = size(C, 1)
    M = size(C, 2)
    R = size(A, 2)
    TILE_DIM = 64

    @inbounds for gj = 1:TILE_DIM:M
        for gk = 1:TILE_DIM:R
            for gi = 1:TILE_DIM:N
                for jj in gj:min(gj + TILE_DIM, M)
                    for k = gk:min(gk + TILE_DIM, R)
                        for ii = gi:min(gi + TILE_DIM, N)
                            C[ii, jj] += A[ii, k] * B[k, jj]
                        end
                    end
                end
            end
        end
    end
end

function poly_array_double(A, out)
    n = size(out, 1)

    @poly_loop for i = 1:n
        for j = 1:n
            out[i, j] = A[i, j]*2
        end
    end
end

function array_double(A, out)
    n = size(out, 1)

    for i = 1:n
        for j = 1:n
            out[i, j] = A[i, j]*2
        end
    end
end

function poly_matmul_double(A, B, out)
    n = size(out, 1)
    m = size(out, 2)
    r = size(A, 2)

    @poly_loop for i = 1:n
        for j = 1:m
            for k = 1:r
                out[i, j] += A[i, k] * B[k, j]
            end
            out[i, j] *= 2
        end
    end
end

function matmul_double(A, B, out)
    n = size(out, 1)
    m = size(out, 2)
    r = size(A, 2)

    for i = 1:n
        for j = 1:m
            for k = 1:r
                out[i, j] += A[i, k] * B[k, j]
            end
            out[i, j] *= 2
        end
    end
end

function poly_dd_stride(arr)
    n = size(arr, 1)
    @poly_loop for i = 1:2:n
        for j = 1:2:n
            arr[i, j] += 1
        end
    end
end

function dd_stride(arr)
    n = size(arr, 1)
    for i = 1:2:n
        for j = 1:2:n
            arr[i, j] += 1
        end
    end
end

function poly_lu(PA, L, U)
    n = size(PA, 1)
    @poly_loop for j=1:n
        L[j, j] = 1.0
        for i=1:j
            s1 = 0.0
            for k=1:i
                s1 += U[k, j] * L[i, k]
            end
            U[i, j] = PA[i, j] - s1
        end

        for ii=j+1:n
            s2 = 0.0
            for kk=1:j
                s2 += U[kk, j] * L[ii, kk]
            end
            L[ii, j] = (PA[ii, j] - s2) / U[j, j]
        end
    end
end

function lu(PA, L, U)
    n = size(PA, 1)
    for j=1:n
        L[j, j] = 1.0
        for i=1:j
            s1 = 0.0
            for k=1:i
                s1 += U[k, j] * L[i, k]
            end
            U[i, j] = PA[i, j] - s1
        end

        for ii=j+1:n
            s2 = 0.0
            for kk=1:j
                s2 += U[kk, j] * L[ii, kk]
            end
            L[ii, j] = (PA[ii, j] - s2) / U[j, j]
        end
    end
end

@testset "performance tests" begin

    @testset "automatic loop reordering" begin
        # test matrix mul wrong loop order (check automatic reordering)
        DIM = 512
        A = rand(DIM, DIM)
        B = rand(DIM, DIM)
        out = zeros(DIM, DIM)

        poly_mul_no_tile(A, B, out) # allow for compiling once

        t_poly = @benchmark poly_mul_no_tile($A, $B, out) setup=(out = zeros($DIM, $DIM))
        t_orig = @benchmark mul($A, $B, out) setup=(out = zeros($DIM, $DIM))
        t_right_order = @benchmark mul_right_order($A, $B, out) setup=(out = zeros($DIM, $DIM))

        jorig = judge(minimum(t_poly), minimum(t_orig))
        jright = judge(minimum(t_poly), minimum(t_right_order))

        @test jorig.time == :improvement # want an improvement over bad order
        @test jright.time != :regression # want comparable or better than right order
    end

    @testset "automatic tiling square" begin
        # test performance nears tiled matrix multiplication (check automatic tiling)
        DIM = 1024
        A = rand(DIM, DIM)
        B = rand(DIM, DIM)
        out = zeros(DIM, DIM)

        poly_mul(A, B, out) # allow for compiling once

        t_poly = @benchmark poly_mul($A, $B, out) setup=(out = zeros($DIM, $DIM))
        t_tile = @benchmark tiled_mul($A, $B, out) setup=(out = zeros($DIM, $DIM))

        j = judge(minimum(t_poly), minimum(t_tile))

        @test j.time != :regression # want improvement or invariant (not regression)
    end

    @testset "automatic tiling non-square" begin
        # test performance nears tiled matrix multiplication (check automatic tiling), for non-square matricies that are multiples of tile dim
        N = 1024
        R = 2048
        M = 512
        A = rand(N, R)
        B = rand(R, M)
        out = zeros(N, M)

        poly_mul(A, B, out) # allow for compiling once

        t_poly = @benchmark poly_mul($A, $B, out) setup=(out = zeros($N, $M))
        t_tile = @benchmark tiled_mul_notsquare($A, $B, out) setup=(out = zeros($N, $M))

        j = judge(minimum(t_poly), minimum(t_tile))

        @test j.time != :regression # want improvement or invariant (not regression)
    end

    @testset "automatic tiling any dimensions" begin
        # test performance nears tiled matrix multiplication (check automatic tiling), for any size matrices
        N = 1020
        R = 900
        M = 981
        A = rand(N, R)
        B = rand(R, M)
        out = zeros(N, M)

        poly_mul(A, B, out) # allow for compiling once

        t_poly = @benchmark poly_mul($A, $B, out) setup=(out = zeros($N, $M))
        t_tile = @benchmark tiled_mul_any($A, $B, out) setup=(out = zeros($N, $M))

        j = judge(minimum(t_poly), minimum(t_tile))

        @test j.time != :regression # want improvement or invariant (not regression)
    end

    """
    All of the following tests check that code always gets faster than the original
    All of the code tested has misordered loops or can improve from tiling
    """

    @testset "array double" begin
        DIM = 512
        A = rand(DIM, DIM)
        out = zeros(DIM, DIM)

        poly_array_double(A, out) # allow for compiling once

        t_poly = @benchmark poly_array_double($A, out) setup=(out = zeros($DIM, $DIM))
        t_orig = @benchmark array_double($A, out) setup=(out = zeros($DIM, $DIM))

        j = judge(minimum(t_poly), minimum(t_orig))

        @test j.time == :improvement
    end

    @testset "matmul double" begin
        DIM = 512
        A = rand(DIM, DIM)
        B = rand(DIM, DIM)
        out = zeros(DIM, DIM)

        poly_matmul_double(A, B, out) # allow for compiling once

        t_poly = @benchmark poly_matmul_double($A, $B, out) setup=(out = zeros($DIM, $DIM))
        t_orig = @benchmark matmul_double($A, $B, out) setup=(out = zeros($DIM, $DIM))

        j = judge(minimum(t_poly), minimum(t_orig))

        @test j.time == :improvement
    end

    @testset "2D stride" begin
        DIM = 512
        arr = ones(DIM, DIM)

        poly_dd_stride(arr) # allow for compiling once

        t_poly = @benchmark poly_dd_stride(arr) setup=(arr = ones($DIM, $DIM))
        t_orig = @benchmark dd_stride(arr) setup=(arr = ones($DIM, $DIM))

        j = judge(minimum(t_poly), minimum(t_orig))

        @test j.time == :improvement

    end

    @testset "lu decomposition" begin
        DIM = 1024

        # this step is hard in plain julia (no functions), so we just do it ahead of time
        A = rand(DIM, DIM)*10
        P = zeros(DIM, DIM)
        for i = 1:DIM
            P[i, i] = 1.0
        end
        for j=1:DIM
            i = findmax(A[:, j])[2]
            if j != i
                P[j, :], P[i, :] = P[i, :], P[j, :]
            end
        end

        PA = P*A
        L = zeros(DIM, DIM)
        U = zeros(DIM, DIM)

        poly_lu(PA, L, U) # allow for compiling once

        t_poly = @benchmark poly_lu($PA, L, U) setup=(L = zeros($DIM, $DIM); U = zeros($DIM, $DIM))
        t_orig = @benchmark lu($PA, L, U) setup=(L = zeros($DIM, $DIM); U = zeros($DIM, $DIM))

        j = judge(minimum(t_poly), minimum(t_orig))

        @test j.time == :improvement
    end

end


text:
struct AlwaysEmpty{T} end
Base.eltype(::Type{AlwaysEmpty{T}}) where {T} = T
Base.isempty(::AlwaysEmpty) = true
putting(::AlwaysEmpty) = Return(Block())
taking(::AlwaysEmpty) = Return(Block())

abstract type ChannelEndpoint{T} end
Base.eltype(::Type{ChannelEndpoint{T}}) where {T} = T
ContextManagers.maybeenter(ch::ChannelEndpoint) = ch

struct SendEndpoint{T,Data,Send} <: ChannelEndpoint{T}
    eltype::Val{T}
    data::Data
    send::Send
    senderclosed::typeof(Promise{Nothing}())
    receiverclosed::typeof(Promise{Nothing}())
end

struct ReceiveEndpoint{T,Data,Receive} <: ChannelEndpoint{T}
    eltype::Val{T}
    data::Data
    receive::Receive
    senderclosed::typeof(Promise{Nothing}())
    receiverclosed::typeof(Promise{Nothing}())
end

function channel_endpoints(data)
    send, receive = Reagents.channel(eltype(data), Nothing)
    senderclosed = Promise{Nothing}()
    receiverclosed = Promise{Nothing}()
    return (
        SendEndpoint(Val(eltype(data)), data, send, senderclosed, receiverclosed),
        ReceiveEndpoint(Val(eltype(data)), data, receive, senderclosed, receiverclosed),
    )
end

Julio.channel(T::Type = Any) = channel_endpoints(AlwaysEmpty{T}())
Julio.stack(T::Type = Any) = channel_endpoints(TreiberStack{T}())
Julio.queue(T::Type = Any) = channel_endpoints(MSQueue{T}())

selfclosed(sept::SendEndpoint) = sept.senderclosed
dualclosed(sept::SendEndpoint) = sept.receiverclosed
selfclosed(rept::ReceiveEndpoint) = rept.receiverclosed
dualclosed(rept::ReceiveEndpoint) = rept.senderclosed

struct NoInput <: Exception end
struct NoReceivers <: Exception end

function tryputting(sept::SendEndpoint)
    (; send, data, receiverclosed) = sept
    if istracing()
        post_send = PostCommit(_ -> @trace(label = :putting_br_send))
        post_data = PostCommit(_ -> @trace(label = :putting_br_data))
    else
        post_send = Identity()
        post_data = Identity()
    end
    return (
        (send ⨟ Return(true) ⨟ post_send) |
        # TODO: check `receiverclosed` in a non-blocking manner before `putting(data)`?
        (putting(data) ⨟ post_data ⨟ Return(true)) |
        (fetching(receiverclosed) ⨟ Return(false))
    )
end

putting(sept::SendEndpoint) =
    tryputting(sept) ⨟ Map(x -> x ? nothing : Error(NoReceivers()))

block_if_nonempty(data) = Computed(x -> isempty(data) ? Return(x) : Block())

function maybetaking(rept::ReceiveEndpoint)
    (; receive, data, senderclosed) = rept
    item = Return(nothing) ⨟ (receive ⨟ block_if_nonempty(data) | taking(data)) ⨟ Map(Some)
    # recheck = item | Return(Error(NoInput()))
    # recheck = ReturnIfBlocked(Error(NoInput())) ⨟ item
    # recheck = Return(Error(NoInput()))
    recheck = taking(data) ⨟ Map(Some) | Return(nothing)
    return item | fetching(senderclosed) ⨟ Return(nothing) ⨟ recheck
end

taking(rept::ReceiveEndpoint) =
    maybetaking(rept) ⨟ Map(x -> x === nothing ? Error(NoInput()) : x)

function Base.close(ept::ChannelEndpoint)
    @trace(label = :close_start, issend = ept isa SendEndpoint)
    wasclosed = event(Julio.tryput!, selfclosed(ept), nothing)()
    @trace(label = :close_done, issend = ept isa SendEndpoint, wasclosed)
    return
end

event(::typeof(put!), sept::SendEndpoint, x) =
    Return(convert(eltype(sept), x)) ⨟ putting(sept)
event(::typeof(push!), sept::SendEndpoint) = event(put!, sept, x)

function Base.put!(sept::SendEndpoint, x)
    @trace label = :put_start x
    apply(put!, sept, x)
    @trace label = :put_done
    return
end

event(::typeof(Julio.maybetake!), rept::ReceiveEndpoint) = maybetaking(rept)
Julio.maybetake!(rept::ReceiveEndpoint) = apply(Julio.maybetake!, rept)

event(::typeof(take!), rept::ReceiveEndpoint) = taking(rept)

function Base.take!(rept::ReceiveEndpoint)
    @trace label = :take_start
    ans = apply(take!, rept)
    @trace label = :take_done x = ans
    return ans
end

function Base.push!(sept::SendEndpoint, x)
    put!(sept, x)
    return sept
end

Base.IteratorSize(::Type{<:ReceiveEndpoint}) = Base.SizeUnknown()

function Base.iterate(rept::ReceiveEndpoint, ::Nothing = nothing)
    @trace label = :iterate_start
    x = cancellable_react!(taking(rept))
    @trace label = :iterate_done ans = x
    if x === Error(NoInput())
        return nothing
    else
        return (something(x), nothing)
    end
end

function on_channel_handle_method_error(io::IO, exc::MethodError, argtypes, _kwargs)
    @nospecialize argtypes _kwargs
    putfns = (put!, Julio.tryput!)
    takefns = (take!, Julio.maybetake!)
    length(argtypes) > 0 || return
    EndpointType = argtypes[1]::Type
    if (
        (EndpointType <: ReceiveEndpoint && exc.f in putfns) ||
        (EndpointType <: SendEndpoint && exc.f in takefns)
    )
        required = exc.f in putfns ? "input" : "output"
        given = required == "input" ? "output" : "input"
        println(io)
        printstyled(io, string(exc.f); color = :cyan)
        print(io, " requires an ")
        printstyled(io, "$required endpoint"; color = :cyan, bold = true)
        print(io, " not an ")
        printstyled(io, "$given endpoint"; color = :cyan)
        print(io, ". ")
        print(
            io,
            """
            Use the other endpoint. For example:

                send_endpoint, receive_endpoint = Julio.channel()
                $(exc.f)($(required)_endpoint$(exc.f in putfns ? ", item" : ""))

            """,
        )
        println(io)
    end
end


text:
# Statistical Voice Conversion based on Spectrum Differential [Kobayashi 2014]

### basic setup

# Load source speaker's (`clb`) speech signal.
path = joinpath(Pkg.dir("VoiceConversion", "test", "data", "clb_a0028.txt"))
x = vec(readdlm(path))
fs = 16000
period = 5.0
order = 40
alpha = mcepalpha(fs)

# Mel-cepstrum extraction based on WORLD.
f0, timeaxis = dio(x, fs, DioOption(period=period))
f0 = stonemask(x, fs, timeaxis, f0)
spectrogram = cheaptrick(x, fs, timeaxis, f0)
src_clb28 = sp2mc(spectrogram, order, alpha)
@test all(isfinite.(src_clb28))

x_clb28 = copy(x)

function diffvc_base(src, mapper)
    # Perform parameter conversion
    converted = vc(mapper, src)
    @test all(isfinite.(converted))

    # remove power coef. in the converted signal
    converted[1,:] = 0.0

    # Waveform modification
    mf = MLSADF(order, alpha)
    hopsize = convert(Int, round(fs / (1000 / period)))
    synthesis!(mf, x_clb28, mc2b(converted, alpha), hopsize)
end

println("testing: voice conversion based on direct waveform modification")

let
    println("Female (`clb`) to female (`slt`) voice conversion")
    println("GMM-based frame-by-frame mapping")
    x = copy(src_clb28)

    # Load GMM to convert speech signal of `clb` to that of `slt`,
    # mixture: 32, order of mel-cepstrum: 40
    modelpath = joinpath(Pkg.dir("VoiceConversion"), "test", "models",
                         "clb_to_slt_gmm32_order40_diff.jld")
    gmm = load(modelpath)
    @assert gmm["diff"]

    # Construct GMM-based frame-by-frame mapping
    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])

    y = diffvc_base(x, mapper)
    @test all(isfinite.(y))
end

let
    println("Female (`clb`) to female (`slt`) voice conversion")
    println("GMM-based trajectory paramter mapping")
    x = copy(src_clb28)

    # add dynamic feature
    x = [reshape(x[1,:], 1, size(x,2)); push_delta(x[2:end,:])]

    # Load GMM to convert speech signal of `clb` to that of `slt`,
    # mixture: 32, order of mel-cepstrum: 40+40 (with delta feature)
    modelpath = joinpath(Pkg.dir("VoiceConversion"), "test", "models",
                         "clb_to_slt_gmm32_order40_diff_with_delta.jld")
    gmm = load(modelpath)
    @assert gmm["diff"]

    # Construct trajectory-based GMM parameter mapping
    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])
    mapper = TrajectoryGMMMap(mapper, 70)

    y = diffvc_base(x, mapper)
    @test all(isfinite.(y))
end


text:
# ================
# Homogeneous case
# ================

# homogeneous case; no invariant; recursive implementation
function reach_homog_INT!(F::Vector{ReachSet{N, Interval{N, IA.Interval{N}}}},
                          Ω0::Interval{N, IA.Interval{N}},
                          Φ::N,
                          NSTEPS::Integer,
                          δ::Float64,
                          X::Universe,
                          Δt0::TimeInterval) where {N}
    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    k = 2
    @inbounds while k <= NSTEPS
        Iₖ = Interval(Φ * set(F[k-1]).dat)
        Δt += δ
        F[k] = ReachSet(Iₖ, Δt)
        k += 1
    end
    return F
end

# homogeneous case; with invariant; recursive implementation
function reach_homog_INT!(F::Vector{ReachSet{N, Interval{N, IA.Interval{N}}}},
                          Ω0::Interval{N, IA.Interval{N}},
                          Φ::N,
                          NSTEPS::Integer,
                          δ::Float64,
                          X::LazySet,
                          Δt0::TimeInterval) where {N}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    k = 2
    @inbounds while k <= NSTEPS
        Iₖ = Interval(Φ * set(F[k-1]).dat)
        _is_intersection_empty(X, Iₖ) && break
        Δt += δ
        F[k] = ReachSet(Iₖ, Δt)
        k += 1
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end


text:
## Compile rand vars into callable functions
## =========================================

# There are two kinds of input a random variable may take as input
# 1. A event, a subset of the sample space
# 2. A psuedo-event which is a subset of a transformed sample space, e.g. U(0,10)

## Situation 1. Normal Sample Space
## ================================

dims(X::ExecutableRandVar) = X.dims
call(X::ExecutableRandVar, ω::Omega) = X.func(ω)
convert{T}(::Type{ExecutableRandVar{T}}, X::SymbolicRandVar{T}) =
  ExecutableRandVar{T}(lambda(X), dims(X))

all_functional_randvars = union(real_real_real, real_real, real_floating,
                                real_real_bool, ((:IfElseRandVar, :ifelse),),
                                bool_bool_bool, ((:NotRandVar, :!),))
for (name,op) in all_functional_randvars
  eval(
  quote
  function convert(::Type{Expr}, X::$name)
    Expr(:call,$op,[convert(Expr,arg) for arg in args(X)]...)
  end
  end)
end

convert(::Type{Expr}, X::ConstantRandVar) = X.val
#add one for julia/c++ indexing mismatch, basically a HACK
convert(::Type{Expr}, X::OmegaRandVar) = :(ω[$(X.dim)])
lambda_expr(X::SymbolicRandVar) = Expr(:(->),:ω,convert(Expr,X))
lambda(X::SymbolicRandVar) = eval(lambda_expr(X))

function convert{T <: ElementaryRandVar}(::Type{Expr}, X::T)
  Expr(:call, :quantile, distribution_type(T), :(ω[$X.dim]), [convert(Expr, arg) for arg in args(X)]...)
end

## Situation 2. Psuedo-Sample space
## ================================
# What kind of this is this new space, is it too  a sample space.
# - well it's certainly not bounded between 0 and 1
# - and it doesn't have a uniform measure
# - What happens if you do uniform(0,0,10) and uniform(0,0,20) we'll get two dimensions
# - sometings which are UNSAt could be deemed SAT,
# - So we should only allow one, in a way it really is constructing a new sample space.

convert_psuedo{T}(::Type{ExecutableRandVar{T}}, X::SymbolicRandVar{T}) =
  ExecutableRandVar{T}(lambda_psuedo(X), dims(X))

for (name,op) in all_functional_randvars
  eval(
  quote
  function convert_psuedo(::Type{Expr}, X::$name)
    Expr(:call,$op,[convert_psuedo(Expr,arg) for arg in args(X)]...)
  end
  end)
end

convert_psuedo(::Type{Expr}, X::ConstantRandVar) = X.val
#add one for julia/c++ indexing mismatch, basically a HACK
convert_psuedo(::Type{Expr}, X::OmegaRandVar) = :(ω[$(X.dim)])
lambda_psuedo_expr(X::SymbolicRandVar) = Expr(:(->),:ω,convert_psuedo(Expr,X))
lambda_psuedo(X::SymbolicRandVar) = eval(lambda_psuedo_expr(X))

# This is the only difference
convert_psuedo{T <: ElementaryRandVar}(::Type{Expr}, X::T) = :(ω[$(X.dim)])

## Rand Array
## ==========
function convert{T}(::Type{ExecutableRandArray{T}}, Xs::RandArray{T,1})
  ExecutableRandArray([convert(ExecutableRandVar{T}, Xs.array[i]) for i = 1:size(Xs.array,1)])
end

function convert{T}(::Type{ExecutableRandArray{T}}, Xs::RandArray{T,2})
  ExecutableRandArray([convert(ExecutableRandVar{T}, Xs.array[i,j]) for i = 1:size(Xs.array,1), j = 1:size(Xs.array,2)])
end

function convert_psuedo{T}(::Type{ExecutableRandArray{T}}, Xs::RandArray{T,1})
  ExecutableRandArray([convert_psuedo(ExecutableRandVar{T}, Xs.array[i]) for i = 1:size(Xs.array,1)])
end

function convert_psuedo{T}(::Type{ExecutableRandArray{T}}, Xs::RandArray{T,2})
  ExecutableRandArray([convert_psuedo(ExecutableRandVar{T}, Xs.array[i,j]) for i = 1:size(Xs.array,1), j = 1:size(Xs.array,2)])
end

function call{T}(Xs::ExecutableRandArray{T,1}, ω::Omega)
  T[call(Xs.array[i],ω) for i = 1:size(Xs.array,1)]
end

function call{T}(Xs::ExecutableRandArray{T,2}, ω::Omega)
  T[call(Xs.array[i,j],ω) for i = 1:size(Xs.array,1), j = 1:size(Xs.array,2)]
end

## Executionalize
executionalize{T,N}(Xs::RandArray{T,N}) = convert_psuedo(ExecutableRandArray{T}, Xs)

executionalize{T}(X::RandVar{T}) = convert_psuedo(ExecutableRandVar{T}, X)


text:
## In LinearAlgebra/src/qr.jl
## the QR factorization stores
## * Q in a compressed form using Householder transformations
## * R as an upper triangular matrix
##
## QRFactorization stores
## * Q in a compressed form using Givens rotations, either descending or twisted
## * R in a compressed form (RankOne, Pencil, or Identity) or as a full matrix that is upper triangular
## * this factorization gets updated by the AMRVW_algorithm!, and hence `eigvals`
##
abstract type AbstractQRRotatorFactorization  end

##  We  have  two types, though could consolidate
##
struct QRFactorization{QF, RF} <: AbstractQRRotatorFactorization
  QF::QF
  RF::RF
  end

# use identity R when not specified
function QRFactorization(QF::AbstractQFactorization{T,S}) where {T,S}
    RF = RFactorizationIdentity{T,S}()
    QRFactorization(QF, RF)
end

function Base.size(state::AbstractQRRotatorFactorization)
    m1, n1 = size(state.QF)
    m2, n2 = size(state.RF)
    return (max(m1,m2), max(n1, n2))
end

Base.copy(F::QRFactorization) = QRFactorization(copy(F.QF), copy(F.RF))
Base.length(state::QRFactorization) = length(state.QF)+1
Base.eltype(state::QRFactorization) = eltype(state.QF)

# return A = QR
function Base.Matrix(state::QRFactorization)

    Q = Matrix(state.QF)
    R = Matrix(state.RF)

    ## May have Q smaller than  R, so will pad in that case
    if R != I
        m,n =  size(Q)[1], size(R)[1]
        n  < m && error("R is  too  small")
        return Q * R[1:m, 1:m]
    else
        return Q
    end

end

"""

    qr_factorization(H::Matrix; unitary=false)

For a Hessenberg matrix `H` return a factorization,`Q⋅R`, where `Q` is a `QFactorization` object of descending rotators and `R` is either a `AMRVW.RFactorizationUnitaryDiagonal` object (if `unitary=true`) or a `RFactorizationUpperTriangular` object.

`H` may be a full matrix or the `.H` compontent  of a `hessenberg` factorization.

"""
function qr_factorization(H::AbstractMatrix; unitary=false)
    R = copy(H)
    qr_factorization!(R; unitary=unitary)
end

function qr_factorization!(R::AbstractMatrix{S}; unitary=false) where {S}

    n = size(R, 1)
    Qs = Vector{Rotator{real(S),S}}(undef, n-1)
    for i in 1:n-1
        c,s,r = givensrot(R[i,i], R[i+1,i])
        Ui =  Rotator(c,s,i)
        Qs[i] = Ui'
        j = i + 1
        R[i,i] = r
        R[j,i] = zero(S)
        for k in j:n
            rik, rjk =  R[i,k],  R[j,k]
            R[i,k] = c * rik + s * rjk
            R[j,k] = -conj(s) * rik + conj(c) * rjk
        end
    end

    _qr_factorization(DescendingChain(Qs), R, Val(unitary))
end

function _qr_factorization(Des::DescendingChain{T,S,V}, R, unitary::Val{true}) where {T, S, V}

    QF = QFactorization(Des)
    RF = RFactorizationUnitaryDiagonal(sign.(diag(R)))

    QRFactorization(QF, RF)
end

function _qr_factorization(Des, R, unitary::Val{false})

    QF = QFactorization(Des)
    RF = RFactorizationUpperTriangular(R)

    F = QRFactorization(QF, RF)
end


text:
@testset "density1D" begin
    df = (x=rand(1000), c=rand(["a", "b"], 1000))
    npoints = 500

    layer = data(df) * mapping(:x, color=:c) * AlgebraOfGraphics.density(; npoints)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    rgx1 = range(extrema(df.x)..., length=npoints)
    d1 = pdf(kde(x1), rgx1)

    x2 = df.x[df.c .== "b"]
    rgx2 = range(extrema(df.x)..., length=npoints)
    d2 = pdf(kde(x2), rgx2)

    rgx, d = processedlayer.positional

    @test rgx[1] ≈ rgx1
    @test d[1] ≈ d1

    @test rgx[2] ≈ rgx2
    @test d[2] ≈ d2

    layer = data(df) * mapping(:x, color=:c) * AlgebraOfGraphics.density(; npoints, datalimits=extrema)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    rgx1 = range(extrema(x1)..., length=npoints)
    d1 = pdf(kde(x1), rgx1)

    x2 = df.x[df.c .== "b"]
    rgx2 = range(extrema(x2)..., length=npoints)
    d2 = pdf(kde(x2), rgx2)

    rgx, d = processedlayer.positional

    @test rgx[1] ≈ rgx1
    @test d[1] ≈ d1

    @test rgx[2] ≈ rgx2
    @test d[2] ≈ d2

    @test processedlayer.primary == NamedArguments((color=["a", "b"],))
    @test isempty(processedlayer.named)
    @test processedlayer.attributes == NamedArguments()
    @test processedlayer.plottype == AlgebraOfGraphics.LinesFill

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, :color, "c")
    insert!(labels, 2, "pdf")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "density2d" begin
    df = (x=rand(1000), y=rand(1000), c=rand(["a", "b"], 1000))
    npoints = 500
    bandwidth = (0.01, 0.01)

    layer = data(df) * mapping(:x, :y, color=:c) * AlgebraOfGraphics.density(; npoints, bandwidth)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    rgx1 = range(extrema(df.x)..., length=npoints)
    rgy1 = range(extrema(df.y)..., length=npoints)
    d1 = pdf(kde((x1, y1); bandwidth), rgx1, rgy1)

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    rgx2 = range(extrema(df.x)..., length=npoints)
    rgy2 = range(extrema(df.y)..., length=npoints)
    d2 = pdf(kde((x2, y2); bandwidth), rgx2, rgy2)

    rgx, rgy, d = processedlayer.positional

    @test rgx[1] ≈ rgx1
    @test rgy[1] ≈ rgy1
    @test d[1] ≈ d1

    @test rgx[2] ≈ rgx2
    @test rgy[2] ≈ rgy2
    @test d[2] ≈ d2

    @test processedlayer.primary == NamedArguments((color=["a", "b"],))
    @test isempty(processedlayer.named)
    @test processedlayer.attributes == NamedArguments()
    @test processedlayer.plottype == Heatmap

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :color, "c")
    insert!(labels, 3, "pdf")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "expectation1d" begin
    df = (x=rand(["a", "b"], 1000), y=rand(1000), c=rand(["a", "b"], 1000))

    layer = data(df) * mapping(:x, :y, layout=:c) * expectation()
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    ma1 = mean(df.y[(df.x .== "a") .& (df.c .== "a")])
    mb1 = mean(df.y[(df.x .== "b") .& (df.c .== "a")])

    ma2 = mean(df.y[(df.x .== "a") .& (df.c .== "b")])
    mb2 = mean(df.y[(df.x .== "b") .& (df.c .== "b")])

    x, m = processedlayer.positional
    x1, m1 = x[1], m[1]
    x2, m2 = x[2], m[2]

    @test x1 == ["a", "b"]
    @test m1 ≈ [ma1, mb1]

    @test x2 == ["a", "b"]
    @test m2 ≈ [ma2, mb2]

    @test processedlayer.primary == NamedArguments((layout=["a", "b"],))
    @test isempty(processedlayer.named)
    @test processedlayer.attributes == NamedArguments()
    @test processedlayer.plottype == BarPlot

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :layout, "c")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "expectation2d" begin
    df = (x=rand(["a", "b"], 1000), y=rand(["a", "b"], 1000), z=rand(1000), c=rand(["a", "b"], 1000))

    layer = data(df) * mapping(:x, :y, :z, layout=:c) * expectation()
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    maa1 = mean(df.z[(df.x .== "a") .& (df.y .== "a") .& (df.c .== "a")])
    mab1 = mean(df.z[(df.x .== "a") .& (df.y .== "b") .& (df.c .== "a")])
    mba1 = mean(df.z[(df.x .== "b") .& (df.y .== "a") .& (df.c .== "a")])
    mbb1 = mean(df.z[(df.x .== "b") .& (df.y .== "b") .& (df.c .== "a")])

    maa2 = mean(df.z[(df.x .== "a") .& (df.y .== "a") .& (df.c .== "b")])
    mab2 = mean(df.z[(df.x .== "a") .& (df.y .== "b") .& (df.c .== "b")])
    mba2 = mean(df.z[(df.x .== "b") .& (df.y .== "a") .& (df.c .== "b")])
    mbb2 = mean(df.z[(df.x .== "b") .& (df.y .== "b") .& (df.c .== "b")])

    x, y, m = processedlayer.positional
    x1, y1, m1 = x[1], y[1], m[1]
    x2, y2, m2 = x[2], y[2], m[2]

    @test x1 == ["a", "b"]
    @test y1 == ["a", "b"]
    @test m1 ≈ [maa1 mab1; mba1 mbb1]

    @test x2 == ["a", "b"]
    @test y2 == ["a", "b"]
    @test m2 ≈ [maa2 mab2; mba2 mbb2]

    @test processedlayer.primary == NamedArguments((layout=["a", "b"],))
    @test isempty(processedlayer.named)
    @test processedlayer.attributes == NamedArguments()
    @test processedlayer.plottype == Heatmap

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, 3, "z")
    insert!(labels, :layout, "c")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "frequency1d" begin
    df = (x=rand(["a", "b"], 1000), c=rand(["a", "b"], 1000))

    layer = data(df) * mapping(:x, layout=:c) * frequency()
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    na1 = sum((df.x .== "a") .& (df.c .== "a"))
    nb1 = sum((df.x .== "b") .& (df.c .== "a"))

    na2 = sum((df.x .== "a") .& (df.c .== "b"))
    nb2 = sum((df.x .== "b") .& (df.c .== "b"))

    x, n = processedlayer.positional
    x1, n1 = x[1], n[1]
    x2, n2 = x[2], n[2]

    @test x1 == ["a", "b"]
    @test n1 ≈ [na1, nb1]

    @test x2 == ["a", "b"]
    @test n2 ≈ [na2, nb2]

    @test processedlayer.primary == NamedArguments((layout=["a", "b"],))
    @test isempty(processedlayer.named)
    @test processedlayer.attributes == NamedArguments()
    @test processedlayer.plottype == BarPlot

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, :layout, "c")
    insert!(labels, 2, "count")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "frequency2d" begin
    df = (x=rand(["a", "b"], 1000), y=rand(["a", "b"], 1000), c=rand(["a", "b"], 1000))

    layer = data(df) * mapping(:x, :y, layout=:c) * frequency()
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    naa1 = sum((df.x .== "a") .& (df.y .== "a") .& (df.c .== "a"))
    nab1 = sum((df.x .== "a") .& (df.y .== "b") .& (df.c .== "a"))
    nba1 = sum((df.x .== "b") .& (df.y .== "a") .& (df.c .== "a"))
    nbb1 = sum((df.x .== "b") .& (df.y .== "b") .& (df.c .== "a"))

    naa2 = sum((df.x .== "a") .& (df.y .== "a") .& (df.c .== "b"))
    nab2 = sum((df.x .== "a") .& (df.y .== "b") .& (df.c .== "b"))
    nba2 = sum((df.x .== "b") .& (df.y .== "a") .& (df.c .== "b"))
    nbb2 = sum((df.x .== "b") .& (df.y .== "b") .& (df.c .== "b"))

    x, y, n = processedlayer.positional
    x1, y1, n1 = x[1], y[1], n[1]
    x2, y2, n2 = x[2], y[2], n[2]

    @test x1 == ["a", "b"]
    @test y1 == ["a", "b"]
    @test n1 ≈ [naa1 nab1; nba1 nbb1]

    @test x2 == ["a", "b"]
    @test y2 == ["a", "b"]
    @test n2 ≈ [naa2 nab2; nba2 nbb2]

    @test processedlayer.primary == NamedArguments((layout=["a", "b"],))
    @test isempty(processedlayer.named)
    @test processedlayer.attributes == NamedArguments()
    @test processedlayer.plottype == Heatmap

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :layout, "c")
    insert!(labels, 3, "count")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "midpoints" begin
    edges = [1, 2, 10, 12]
    @test midpoints(edges) ≈ [1.5, 6, 11]

    edges_rg = 1:2:5
    edges_v = [1, 3, 5]
    @test midpoints(edges_v) ≈ midpoints(edges_rg) ≈ [2, 4]
end

@testset "histogram1D" begin
    df = (x=rand(1000), c=rand(["a", "b"], 1000))
    bins = 0:0.01:1

    layer = data(df) * mapping(:x, color=:c) * histogram(; bins)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    w1 = fit(Histogram, x1, bins).weights

    x2 = df.x[df.c .== "b"]
    w2 = fit(Histogram, x2, bins).weights

    rgx, w = processedlayer.positional
    width = processedlayer.named[:width]

    @test rgx[1] ≈ (bins[1:end-1] .+ bins[2:end]) ./ 2
    @test w[1] == w1
    @test width[1] ≈ diff(bins)

    @test rgx[2] ≈ (bins[1:end-1] .+ bins[2:end]) ./ 2
    @test w[2] == w2
    @test width[2] ≈ diff(bins)

    bins, closed = 12, :left
    layer = data(df) * mapping(:x, color=:c) * histogram(; bins, closed, datalimits=extrema)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    bins1 = histrange(extrema(x1)..., bins, closed)
    w1 = fit(Histogram, x1, bins1).weights

    x2 = df.x[df.c .== "b"]
    bins2 = histrange(extrema(x2)..., bins, closed)
    w2 = fit(Histogram, x2, bins2).weights

    rgx, w = processedlayer.positional
    width = processedlayer.named[:width]

    @test rgx[1] ≈ (bins1[1:end-1] .+ bins1[2:end]) ./ 2
    @test w[1] ≈ w1
    @test width[1] ≈ diff(bins1)

    @test rgx[2] ≈ (bins2[1:end-1] .+ bins2[2:end]) ./ 2
    @test w[2] ≈ w2
    @test width[2] ≈ diff(bins2)

    @test processedlayer.primary == NamedArguments((color=["a", "b"],))
    @test processedlayer.attributes == NamedArguments((gap=0, dodge_gap=0))
    @test keys(processedlayer.named) == Indices([:width])
    @test processedlayer.plottype == AlgebraOfGraphics.BarPlot

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, :color, "c")
    insert!(labels, 2, "count")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)

    bins = 12.3
    layer = data(df) * mapping(:x, color=:c) * histogram(; bins)
    @test_throws ArgumentError AlgebraOfGraphics.ProcessedLayer(layer)
end

@testset "weightedhistogram1d" begin
    df = (x=rand(1000), z=rand(1000), c=rand(["a", "b"], 1000))
    bins = collect(0:0.01:1) # test vector of bins

    layer = data(df) * mapping(:x, color=:c, weights=:z) * histogram(; bins)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    z1 = df.z[df.c .== "a"]
    w1 = fit(Histogram, x1, weights(z1), bins).weights

    x2 = df.x[df.c .== "b"]
    z2 = df.z[df.c .== "b"]
    w2 = fit(Histogram, x2, weights(z2), bins).weights

    rgx, w = processedlayer.positional
    width = processedlayer.named[:width]

    @test rgx[1] ≈ (bins[1:end-1] .+ bins[2:end]) ./ 2
    @test w[1] ≈ w1
    @test width[1] ≈ diff(bins)

    @test rgx[2] ≈ (bins[1:end-1] .+ bins[2:end]) ./ 2
    @test w[2] ≈ w2
    @test width[2] ≈ diff(bins)

    @test processedlayer.primary == NamedArguments((color=["a", "b"],))
    @test keys(processedlayer.named) == Indices([:width])
    @test processedlayer.attributes == NamedArguments((gap=0, dodge_gap=0))
    @test processedlayer.plottype == AlgebraOfGraphics.BarPlot

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, :color, "c")
    insert!(labels, :weights, "z")
    insert!(labels, 2, "count")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "histogram2d" begin
    df = (x=rand(1000), y=rand(1000), c=rand(["a", "b"], 1000))
    binsx, binsy = 0:0.01:1, 0:0.02:1
    bins = (binsx, binsy)

    layer = data(df) * mapping(:x, :y, layout=:c) * histogram(; bins)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    w1 = fit(Histogram, (x1, y1), bins).weights

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    w2 = fit(Histogram, (x2, y2), bins).weights

    rgx, rgy, w = processedlayer.positional

    @test rgx[1] ≈ (binsx[1:end-1] .+ binsx[2:end]) ./ 2
    @test rgy[1] ≈ (binsy[1:end-1] .+ binsy[2:end]) ./ 2
    @test w[1] == w1

    @test rgx[2] ≈ (binsx[1:end-1] .+ binsx[2:end]) ./ 2
    @test rgy[2] ≈ (binsy[1:end-1] .+ binsy[2:end]) ./ 2
    @test w[2] == w2

    bins, closed = 12, :left
    layer = data(df) * mapping(:x, :y, layout=:c) * histogram(; bins, closed, datalimits=extrema)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    binsx1 = histrange(extrema(x1)..., bins, closed)
    binsy1 = histrange(extrema(y1)..., bins, closed)
    w1 = fit(Histogram, (x1, y1), (binsx1, binsy1)).weights

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    binsx2 = histrange(extrema(x2)..., bins, closed)
    binsy2 = histrange(extrema(y2)..., bins, closed)
    w2 = fit(Histogram, (x2, y2), (binsx2, binsy2)).weights

    rgx, rgy, w = processedlayer.positional

    @test rgx[1] ≈ (binsx1[1:end-1] .+ binsx1[2:end]) ./ 2
    @test rgy[1] ≈ (binsy1[1:end-1] .+ binsy1[2:end]) ./ 2
    @test w[1] == w1

    @test rgx[2] ≈ (binsx2[1:end-1] .+ binsx2[2:end]) ./ 2
    @test rgy[2] ≈ (binsy2[1:end-1] .+ binsy2[2:end]) ./ 2
    @test w[2] == w2

    @test processedlayer.primary == NamedArguments((layout=["a", "b"],))
    @test isempty(processedlayer.named)
    @test isempty(processedlayer.attributes)
    @test processedlayer.plottype == AlgebraOfGraphics.Heatmap

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :layout, "c")
    insert!(labels, 3, "count")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)

    bins = rand(2, 2)
    layer = data(df) * mapping(:x, color=:c) * histogram(; bins)
    @test_throws ArgumentError AlgebraOfGraphics.ProcessedLayer(layer)
end

@testset "weightedhistogram2d" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), c=rand(["a", "b"], 1000))
    binsx, binsy = 0:0.01:1, 0:0.02:1
    bins = (binsx, binsy)

    layer = data(df) * mapping(:x, :y, layout=:c, weights=:z) * histogram(; bins)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    z1 = df.z[df.c .== "a"]
    w1 = fit(Histogram, (x1, y1), weights(z1), bins).weights

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    z2 = df.z[df.c .== "b"]
    w2 = fit(Histogram, (x2, y2), weights(z2), bins).weights

    rgx, rgy, w = processedlayer.positional

    @test rgx[1] ≈ (binsx[1:end-1] .+ binsx[2:end]) ./ 2
    @test rgy[1] ≈ (binsy[1:end-1] .+ binsy[2:end]) ./ 2
    @test w[1] == w1

    @test rgx[2] ≈ (binsx[1:end-1] .+ binsx[2:end]) ./ 2
    @test rgy[2] ≈ (binsy[1:end-1] .+ binsy[2:end]) ./ 2
    @test w[2] == w2

    @test processedlayer.primary == NamedArguments((layout=["a", "b"],))
    @test isempty(processedlayer.named)
    @test isempty(processedlayer.attributes)
    @test processedlayer.plottype == AlgebraOfGraphics.Heatmap

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :layout, "c")
    insert!(labels, :weights, "z")
    insert!(labels, 3, "count")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "intercept" begin
    x = rand(10)
    mat = AlgebraOfGraphics.add_intercept_column(x)
    @test mat == [ones(10) x]
    @test eltype(mat) == Float64

    x = rand(1:3, 10)
    mat = AlgebraOfGraphics.add_intercept_column(x)
    @test mat == [ones(10) x]
    @test eltype(mat) == Float64
end

@testset "linear" begin
    df = (x=rand(1000), y=rand(1000), c=rand(["a", "b"], 1000))
    npoints, dropcollinear = 150, false
    layer = data(df) * mapping(:x, :y, color=:c) * linear(; npoints, dropcollinear)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    lm1 = GLM.lm([fill(one(eltype(x1)), length(x1)) x1], y1; dropcollinear)
    x̂1 = range(extrema(x1)...; length=npoints)
    ŷ1, lower1, upper1 = map(vec, GLM.predict(lm1, [ones(length(x̂1)) x̂1]; interval=:confidence))

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    lm2 = GLM.lm([fill(one(eltype(x2)), length(x2)) x2], y2; dropcollinear)
    x̂2 = range(extrema(x2)...; length=npoints)
    ŷ2, lower2, upper2 = map(vec, GLM.predict(lm2, [ones(length(x̂2)) x̂2]; interval=:confidence))

    x̂, ŷ = processedlayer.positional
    lower, upper = processedlayer.named[:lower], processedlayer.named[:upper]

    @test x̂[1] ≈ x̂1
    @test ŷ[1] ≈ ŷ1
    @test lower[1] ≈ lower1
    @test upper[1] ≈ upper1

    @test x̂[2] ≈ x̂2
    @test ŷ[2] ≈ ŷ2
    @test lower[2] ≈ lower2
    @test upper[2] ≈ upper2

    @test processedlayer.primary == NamedArguments((color=["a", "b"],))
    @test isempty(processedlayer.attributes)

    @test processedlayer.plottype == LinesFill

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :color, "c")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)

    # Test `interval` and `level` custom values
    df = (x=rand(1000), y=rand(1000), c=rand(["a", "b"], 1000))
    npoints, dropcollinear = 150, false
    interval, level = :prediction, 0.9
    layer = data(df) * mapping(:x, :y, color=:c) * linear(; npoints, dropcollinear, interval, level)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    lm1 = GLM.lm([fill(one(eltype(x1)), length(x1)) x1], y1; dropcollinear)
    x̂1 = range(extrema(x1)...; length=npoints)
    ŷ1, lower1, upper1 = map(vec, GLM.predict(lm1, [ones(length(x̂1)) x̂1]; interval, level))

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    lm2 = GLM.lm([fill(one(eltype(x2)), length(x2)) x2], y2; dropcollinear)
    x̂2 = range(extrema(x2)...; length=npoints)
    ŷ2, lower2, upper2 = map(vec, GLM.predict(lm2, [ones(length(x̂2)) x̂2]; interval, level))

    x̂, ŷ = processedlayer.positional
    lower, upper = processedlayer.named[:lower], processedlayer.named[:upper]

    @test x̂[1] ≈ x̂1
    @test ŷ[1] ≈ ŷ1
    @test lower[1] ≈ lower1
    @test upper[1] ≈ upper1

    @test x̂[2] ≈ x̂2
    @test ŷ[2] ≈ ŷ2
    @test lower[2] ≈ lower2
    @test upper[2] ≈ upper2
end

@testset "weightedlinear" begin
    df = (x=rand(1000), y=rand(1000), z=rand(1000), c=rand(["a", "b"], 1000))
    npoints, dropcollinear = 150, false
    layer = data(df) * mapping(:x, :y, color=:c, weights=:z) * linear(; npoints, dropcollinear)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    z1 = df.z[df.c .== "a"]
    lm1 = GLM.lm([fill(one(eltype(x1)), length(x1)) x1], y1; dropcollinear, wts=z1)
    x̂1 = range(extrema(x1)...; length=npoints)
    ŷ1 = vec(GLM.predict(lm1, [ones(length(x̂1)) x̂1]; interval=nothing))

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    z2 = df.z[df.c .== "b"]
    lm2 = GLM.lm([fill(one(eltype(x2)), length(x2)) x2], y2; dropcollinear, wts=z2)
    x̂2 = range(extrema(x2)...; length=npoints)
    ŷ2 = vec(GLM.predict(lm2, [ones(length(x̂2)) x̂2]; interval=nothing))

    x̂, ŷ = processedlayer.positional

    @test x̂[1] ≈ x̂1
    @test ŷ[1] ≈ ŷ1

    @test x̂[2] ≈ x̂2
    @test ŷ[2] ≈ ŷ2

    @test processedlayer.primary == NamedArguments((color=["a", "b"],))
    @test isempty(processedlayer.named)
    @test isempty(processedlayer.attributes)

    @test processedlayer.plottype == Lines

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :color, "c")
    insert!(labels, :weights, "z")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
end

@testset "smooth" begin
    df = (x=rand(1000), y=rand(1000), c=rand(["a", "b"], 1000))
    npoints = 150
    layer = data(df) * mapping(:x, :y, color=:c) * smooth(; npoints)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    loess1 = Loess.loess(x1, y1)
    x̂1 = range(extrema(x1)...; length=npoints)
    ŷ1 = Loess.predict(loess1, x̂1)

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    loess2 = Loess.loess(x2, y2)
    x̂2 = range(extrema(x2)...; length=npoints)
    ŷ2 = Loess.predict(loess2, x̂2)

    x̂, ŷ = processedlayer.positional

    @test x̂[1] ≈ x̂1
    @test ŷ[1] ≈ ŷ1

    @test x̂[2] ≈ x̂2
    @test ŷ[2] ≈ ŷ2

    @test processedlayer.primary == NamedArguments((color=["a", "b"],))
    @test isempty(processedlayer.named)
    @test isempty(processedlayer.attributes)

    @test processedlayer.plottype == Lines

    labels = MixedArguments()
    insert!(labels, 1, "x")
    insert!(labels, 2, "y")
    insert!(labels, :color, "c")
    @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)

    # Also test integer input
    df = (x=rand(1:3, 1000), y=rand(1000), c=rand(["a", "b"], 1000))
    npoints = 150
    layer = data(df) * mapping(:x, :y, color=:c) * smooth(; npoints)
    processedlayer = AlgebraOfGraphics.ProcessedLayer(layer)

    x1 = df.x[df.c .== "a"]
    y1 = df.y[df.c .== "a"]
    loess1 = Loess.loess(x1, y1)
    x̂1 = range(extrema(x1)...; length=npoints)
    ŷ1 = Loess.predict(loess1, x̂1)

    x2 = df.x[df.c .== "b"]
    y2 = df.y[df.c .== "b"]
    loess2 = Loess.loess(x2, y2)
    x̂2 = range(extrema(x2)...; length=npoints)
    ŷ2 = Loess.predict(loess2, x̂2)

    x̂, ŷ = processedlayer.positional

    @test x̂[1] ≈ x̂1
    @test ŷ[1] ≈ ŷ1

    @test x̂[2] ≈ x̂2
    @test ŷ[2] ≈ ŷ2
end


text:
mutable struct Lattice{H<:Number} # H = HoppingEltype
  dim::Int
  sites::Int
  L::Int
  n_neighbors::Int
  n_bonds::Int
  t::Matrix{Float64} # colidx = flavor/band, rowidx = hor,ver
  tN::Matrix{Float64} # colidx = flavor/band, rowidx = just one (Yoni txy)
  tNN::Matrix{Float64} # colidx = flavor/band, rowidx = hor,ver
  time_neighbors::Matrix{Int} # colidx = slice, rowidx = up, down
  neighbors::Matrix{Int} # colidx = site, rowidx = up down right left
  Nneighbors::Matrix{Int} # ur, dr, dl, ul
  NNneighbors::Matrix{Int} # uu, rr, dd, ll

  bonds::Matrix{Int} # src, trg, type
  bond_vecs::Matrix{Float64}

  peirls::Matrix{Matrix{Float64}} # phis and NOT exp(im*phi); cols = flavor, row = spin; each inner matrix: trg, src (linidx)

  hopping_matrix::Matrix{H} # mu included
  hopping_matrix_exp::Matrix{H} # mu included
  hopping_matrix_exp_inv::Matrix{H}

  chkr_hop_half::Vector{SparseMatrixCSC{H, Int64}}
  chkr_hop_half_inv::Vector{SparseMatrixCSC{H, Int64}}
  chkr_hop_half_dagger::Vector{SparseMatrixCSC{H, Int64}}

  # -- chkr groups folding (regrouping)
  n_folded::Int
  folded::Vector{UnitRange{Int64}}

  chkr_hop_half_folded::Vector{SparseMatrixCSC{H, Int64}}
  chkr_hop_half_inv_folded::Vector{SparseMatrixCSC{H, Int64}}
  chkr_hop_half_dagger_folded::Vector{SparseMatrixCSC{H, Int64}}

  chkr_hop_half_folded_rev::Vector{SparseMatrixCSC{H, Int64}}
  chkr_hop_half_inv_folded_rev::Vector{SparseMatrixCSC{H, Int64}}
  chkr_hop_half_dagger_folded_rev::Vector{SparseMatrixCSC{H, Int64}}
  # --

  chkr_hop::Vector{SparseMatrixCSC{H, Int64}} # without prefactor 0.5 in matrix exponentials
  chkr_hop_inv::Vector{SparseMatrixCSC{H, Int64}}
  chkr_hop_dagger::Vector{SparseMatrixCSC{H, Int64}}
  chkr_mu_half::SparseMatrixCSC{H, Int64}
  chkr_mu_half_inv::SparseMatrixCSC{H, Int64}
  chkr_mu::SparseMatrixCSC{H, Int64}
  chkr_mu_inv::SparseMatrixCSC{H, Int64}

  # Currently NOT used in checkerboard.jl (but generic chkr)
  checkerboard::Matrix{Int} # src, trg, bondid
  groups::Vector{UnitRange}
  n_groups::Int

  Lattice{H}() where H = new{H}()
end

function load_lattice(mc::AbstractDQMC)
  mc.l.L = mc.p.L
  mc.p.hoppings != "none" && (mc.l.t = reshape([parse(Float64, f) for f in split(mc.p.hoppings, ',')],(2,2)))
  mc.p.Nhoppings != "none" && (mc.l.tN = reshape([parse(Float64, f) for f in split(mc.p.Nhoppings, ',')],(1,2)))
  mc.p.NNhoppings != "none" && (mc.l.tNN = reshape([parse(Float64, f) for f in split(mc.p.NNhoppings, ',')],(2,2)))

  init_lattice_from_filename(mc.p.lattice_file, mc.l)
  init_neighbors_table(mc)
  init_time_neighbors_table(mc)
  init_hopping_matrices(mc)
end

function init_lattice_from_filename(filename::String, l::Lattice)
  xdoc = parse_file(filename)
  graph = LightXML.root(xdoc)
  l.sites = 1

  l.sites = parse(Int, attribute(graph, "vertices"; required=true))
  l.dim = parse(Int, attribute(graph, "dimension"; required=true))

  edges = get_elements_by_tagname(graph, "EDGE")
  l.n_bonds = length(edges)

  # bonds & bond vectors
  println("\nLoading lattice with ", l.sites , " sites")
  l.bonds = zeros(l.n_bonds, 4)
  l.bond_vecs = zeros(l.n_bonds, l.dim)
  v = Vector{Float64}(undef, l.dim)
  for (i, edge) in enumerate(edges)
    src = 0
    trg = 0
    src = parse(Int, attribute(edge, "source"; required=true))
    trg = parse(Int, attribute(edge, "target"; required=true))
    typ = parse(Int, attribute(edge, "type"; required=true))
    id = parse(Int, attribute(edge, "id"; required=true))
    v = [parse(Float64, f) for f in split(attribute(edge, "vector"; required=true)," ")]

    if id != i error("Edges in lattice file must be sorted from 1 to N!") end

    l.bonds[i, 1] = src
    l.bonds[i, 2] = trg
    l.bonds[i, 3] = typ
    l.bond_vecs[i, :] = v
  end

  l.n_neighbors = count(x->x==1, l.bonds[:, 1]) + count(x->x==1, l.bonds[:, 2]) # neighbors of site 1
end

function init_neighbors_table(mc::AbstractDQMC)
  l = mc.l

  println("Initializing neighbor-tables")
  sql = reshape(1:l.sites, l.L, l.L)

  # Nearest neighbors
  if mc.p.hoppings != "none"
    up = circshift(sql,(-1,0))
    right = circshift(sql,(0,-1))
    down = circshift(sql,(1,0))
    left = circshift(sql,(0,1))
    l.neighbors = vcat(up[:]',right[:]',down[:]',left[:]')
  end

  # Next-nearest neighbors (xy, yx)
  if mc.p.Nhoppings != "none"
    ur = circshift(sql,(-1,-1))
    dr = circshift(sql,(1,-1))
    dl = circshift(sql,(1,1))
    ul = circshift(sql,(-1,1))
    l.Nneighbors = vcat(ur[:]',dr[:]',dl[:]',ul[:]')
  end

  # Next-next nearest neighbors (xx, yy)
  if mc.p.NNhoppings != "none"
    uu = circshift(sql,(-2,0))
    rr = circshift(sql,(0,-2))
    dd = circshift(sql,(2,0))
    ll = circshift(sql,(0,2))
    l.NNneighbors = vcat(uu[:]',rr[:]',dd[:]',ll[:]')
  end
  nothing
end

"""
Periodic boundary conditions in imaginary time
"""
function init_time_neighbors_table(mc::AbstractDQMC)
  l = mc.l
  p = mc.p

  l.time_neighbors = zeros(Int64, 2, p.slices)
  for s in 1:p.slices
    l.time_neighbors[1,s] = s == p.slices ? 1 : s+1
    l.time_neighbors[2,s] = s == 1 ? p.slices : s-1
  end
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
@testset "Matrix" begin
    @testset "Matrix: Determinant" begin
        M1 = [1 0; 2 2]
        M2 = rand(3, 3)
        M3 = rand(4, 4)

        @test determinant(M1) == det(M1)
        @test round(determinant(M2), digits = 4) == round(det(M2), digits = 4)
        @test round(determinant(M3), digits = 4) == round(det(M3), digits = 4)
    end

    @testset "Matrix: LU Decompose" begin
        mat = [
            2 -1 -2
            -4 6 3
            -4 -2 8
        ]

        L = [
            1 0 0
            -2 1 0
            -2 -1 1
        ]
        U = [
            2 -1 -2
            0 4 -1
            0 0 3
        ]

        @test lu_decompose(mat) == (L, U)
    end

    @testset "Matrix: rotation matrix" begin
        theta = pi / 6
        a = [1; 0]
        b = [0; 1]
        @test rotation_matrix(theta) * a == [cos(theta); sin(theta)]
        @test rotation_matrix(theta) * b == [-sin(theta); cos(theta)]
    end

    @testset "Matrix: Gauss-Jordan Elimination" begin
        M1 = [1 2 3; 4 5 6]
        M2 = [1 2 3; 4 8 12]
        M3 = [1 2; 4 8]
        M4 = [
            1 4 8 1 4
            4 5 6 8 11
            1 3 2 4 8
            4 5 67 23 0
        ]
        M5 = Float64[
            1 4 8 1 4
            4 5 6 8 11
            1 3 2 4 8
            4 5 67 23 0
        ]
        M6 = [
            1 4 8
            4 0 6
            1 3 2
        ]

        R3 = [
            1.0 0.0 0.0 0.0 -1.10637
            0.0 1.0 0.0 0.0 1.89743
            0.0 0.0 1.0 0.0 -0.444913
            0.0 0.0 0.0 1.0 1.07598
        ]

        @test gauss_jordan(M1) == Float64[1 0 -1; 0 1 2]
        @test_throws AssertionError gauss_jordan(M2)
        @test_throws AssertionError gauss_jordan(M3)
        @test isapprox(gauss_jordan(M4), R3, atol = 1e-5)
        @test isapprox(gauss_jordan(M5), R3, atol = 1e-5)
        @test isapprox(
            gauss_jordan(M6),
            Float64[1 0 0; 0 1 0; -0 -0 1],
            atol = 1e-5,
        )
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
# # Sugarscape
# **Growing Artificial Societies**

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../sugarvis.mp4" type="video/mp4">
# </video>
# ```

# (Descriptions below are from [this page](http://jasss.soc.surrey.ac.uk/12/1/6/appendixB/EpsteinAxtell1996.html))

# ---

# "Growing Artificial Societies" (Epstein & Axtell 1996) is a reference book for scientists
# interested in agent-based modelling and computer simulation. It represents one of the most
# paradigmatic and fascinating examples of the so-called generative approach to social science
# (Epstein 1999). In their book, Epstein & Axtell (1996) present a computational model where
# a heterogeneous population of autonomous agents compete for renewable resources that are
# unequally distributed over a 2-dimensional environment. Agents in the model are autonomous
# in that they are not governed by any central authority and they are heterogeneous in that
# they differ in their genetic attributes and their initial environmental endowments (e.g.
# their initial location and wealth). The model grows in complexity through the different
# chapters of the book as the agents are given the ability to engage in new activities such
# as sex, cultural exchange, trade, combat, disease transmission, etc. The core of Sugarscape
# has provided the basis for various extensions to study e.g. norm formation through cultural
# diffusion (Flentge et al. 2001) and the emergence of communication and cooperation in
# artificial societies (Buzing et al. 2005). Here we analyse the model described in the
# second chapter of Epstein & Axtell's (1996) book within the Markov chain framework.

# ## Rules of sugarscape

# The first model that Epstein & Axtell (1996) present comprises a finite population of
# agents who live in an environment. The environment is represented by a two-dimensional
# grid which contains sugar in some of its cells, hence the name Sugarscape. Agents' role
# in this first model consists in wandering around the Sugarscape harvesting the greatest
# amount of sugar they can find.

# ### Environment

# The environment is a 50×50 grid that wraps around forming a torus. Grid positions have
# both a sugar level and a sugar capacity c. A cell's sugar level is the number of units of
# sugar in the cell (potentially none), and its sugar capacity c is the maximum value the
# sugar level can take on that cell. Sugar capacity is fixed for each individual cell and
# may be different for different cells. The spatial distribution of sugar capacities depicts
# a sugar topography consisting of two peaks (with sugar capacity c = 4) separated by a
# valley, and surrounded by a desert region of sugarless cells (see Figure 1). Note,
# however, that the grid wraps around in both directions.

# The Sugarscape obbeys the following rule:

# Sugarscape growback rule G$\alpha$:
# At each position, sugar grows back at a rate of $\alpha$ units per time-step up to the cell's capacity c.

# ### Agents

# Every agent is endowed with individual (life-long) characteristics that condition her
# skills and capacities to survive in the Sugarscape. These individual attributes are:

# * A vision _v_, which is the maximum number of positions the agent can see in each of
# the four principal lattice directions: north, south, east and west.
# * A metabolic rate _m_, which represents the units of sugar the agent burns per time-step.
# * A maximum age _max-age_, which is the maximum number of time-steps the agent can live.

# Agents also have the capacity to accumulate sugar wealth _w_. An agent's sugar wealth
# is incremented at the end of each time-step by the sugar collected and decremented by the
# agent's metabolic rate. __Two agents are not allowed to occupy the same position in the grid.__

# The agents' behaviour is determined by the following two rules:

# #### Agent movement rule _M_:

# Consider the set of unoccupied positions within your vision (including the one you are
# standing on), identify the one(s) with the greatest amount of sugar, select the nearest one
# (randomly if there is more than one), move there and collect all the sugar in it.
# At this point, the agent's accumulated sugar wealth is incremented by the sugar collected
# and decremented by the agent's metabolic rate _m_. If at this moment the agent's sugar
# wealth is not greater than zero, then the agent dies.

# #### Agent replacement rule _R_:

# Whenever an agent dies it is replaced by a new agent of age 0 placed on a randomly chosen unoccupied position, having random attributes _v_, _m_ and _max-age_, and random initial wealth w0. All random numbers are drawn from uniform distributions with ranges specified in Table 1 below.

# ### Scheduling of events

# Scheduling is determined by the order in which the different rules _G_, _M_ and _R_ are fired in the model. Environmental rule _G_ comes first, followed by agent rule _M_ (which is executed by all agents in random order) and finally agent rule _R_ is executed (again, by all agents in random order).

# ### Parameterisation

# Our analysis corresponds to a model used by Epstein & Axtell (1996, pg. 33) to study the emergent wealth distribution in the agent population. This model is parameterised as indicated in Table 1 below (where U[a,b] denotes a uniform distribution with range [a,b]).

# Initially, each position of the Sugarscape contains a sugar level equal to its sugar capacity c, and the 250 agents are created at a random unoccupied initial location and with random attributes (using the uniform distributions indicated in Table 1).

# __Table 1__

# | Parameter                                | Value        |
# |------------------------------------------|--------------|
# | Lattice length L                         | 50           |
# | Number of sugar peaks                    | 2            |
# | Growth rate $\alpha$                     | 1            |
# | Number of agents N                       | 250          |
# | Agents' initial wealth w0 distribution   | U[5,25]      |
# | Agents' metabolic rate m distribution    | U[1,4]       |
# | Agents' vision v distribution            | U[1,6]       |
# | Agents' maximum age max-age distribution | U[60,100]    |

# ## Creating the ABM

using Agents, Random

mutable struct SugarSeeker <: AbstractAgent
    id::Int
    pos::Dims{2}
    vision::Int
    metabolic_rate::Int
    age::Int
    max_age::Int
    wealth::Int
end

# Functions `distances` and `sugar_caps` produce a matrix for the distribution of sugar capacities."

function distances(pos, sugar_peaks, max_sugar)
    all_dists = Array{Int,1}(undef, length(sugar_peaks))
    for (ind, peak) in enumerate(sugar_peaks)
        d = round(Int, sqrt(sum((pos .- peak) .^ 2)))
        all_dists[ind] = d
    end
    return minimum(all_dists)
end

function sugar_caps(dims, sugar_peaks, max_sugar, dia = 4)
    sugar_capacities = zeros(Int, dims)
    for i in 1:dims[1], j in 1:dims[2]
        sugar_capacities[i, j] = distances((i, j), sugar_peaks, max_sugar)
    end
    for i in 1:dims[1]
        for j in 1:dims[2]
            sugar_capacities[i, j] = max(0, max_sugar - (sugar_capacities[i, j] ÷ dia))
        end
    end
    return sugar_capacities
end

"Create a sugarscape ABM"
function sugarscape(;
    dims = (50, 50),
    sugar_peaks = ((10, 40), (40, 10)),
    growth_rate = 1,
    N = 250,
    w0_dist = (5, 25),
    metabolic_rate_dist = (1, 4),
    vision_dist = (1, 6),
    max_age_dist = (60, 100),
    max_sugar = 4,
    seed = 42
)
    sugar_capacities = sugar_caps(dims, sugar_peaks, max_sugar, 6)
    sugar_values = deepcopy(sugar_capacities)
    space = GridSpace(dims)
    properties = Dict(
        :growth_rate => growth_rate,
        :N => N,
        :w0_dist => w0_dist,
        :metabolic_rate_dist => metabolic_rate_dist,
        :vision_dist => vision_dist,
        :max_age_dist => max_age_dist,
        :sugar_values => sugar_values,
        :sugar_capacities => sugar_capacities,
    )
    model = AgentBasedModel(
        SugarSeeker,
        space,
        scheduler = Schedulers.randomly,
        properties = properties,
        rng = MersenneTwister(seed)
    )
    for ag in 1:N
        add_agent_single!(
            model,
            rand(model.rng, vision_dist[1]:vision_dist[2]),
            rand(model.rng, metabolic_rate_dist[1]:metabolic_rate_dist[2]),
            0,
            rand(model.rng, max_age_dist[1]:max_age_dist[2]),
            rand(model.rng, w0_dist[1]:w0_dist[2]),
        )
    end
    return model
end

model = sugarscape()

# Let's plot the spatial distribution of sugar capacities in the Sugarscape.
using CairoMakie
CairoMakie.activate!() # hide

fig = Figure(resolution = (600, 600))
ax, hm = heatmap(fig[1,1], model.sugar_capacities; colormap=cgrad(:thermal))
Colorbar(fig[1, 2], hm, width = 20)
fig

# ## Defining stepping functions

# Now we define the stepping functions that handle the time evolution of the model

function model_step!(model)
    ## At each position, sugar grows back at a rate of $\alpha$ units
    ## per time-step up to the cell's capacity c.
    togrow = findall(
        x -> model.sugar_values[x] < model.sugar_capacities[x],
        1:length(positions(model)),
    )
    model.sugar_values[togrow] .+= model.growth_rate
end

function movement!(agent, model)
    newsite = agent.pos
    ## find all unoccupied position within vision
    neighbors = nearby_positions(agent.pos, model, agent.vision)
    empty = collect(empty_positions(model))
    if length(empty) > 0
        ## identify the one(s) with greatest amount of sugar
        available_sugar = (model.sugar_values[x,y] for (x, y) in empty)
        maxsugar = maximum(available_sugar)
        if maxsugar > 0
            sugary_sites_inds = findall(x -> x == maxsugar, collect(available_sugar))
            sugary_sites = empty[sugary_sites_inds]
            ## select the nearest one (randomly if more than one)
            for dia in 1:(agent.vision)
                np = nearby_positions(agent.pos, model, dia)
                suitable = intersect(np, sugary_sites)
                if length(suitable) > 0
                    newsite = rand(model.rng, suitable)
                    break
                end
            end
            ## move there and collect all the sugar in it
            newsite != agent.pos && move_agent!(agent, newsite, model)
        end
    end
    ## update wealth (collected - consumed)
    agent.wealth += (model.sugar_values[newsite...] - agent.metabolic_rate)
    model.sugar_values[newsite...] = 0
    ## age
    agent.age += 1
end

function replacement!(agent, model)
    ## If the agent's sugar wealth become zero or less, it dies
    if agent.wealth <= 0 || agent.age >= agent.max_age
        kill_agent!(agent, model)
        ## Whenever an agent dies, a young one is added to a random pos.
        ## New agent has random attributes
        add_agent_single!(
            model,
            rand(model.rng, model.vision_dist[1]:model.vision_dist[2]),
            rand(model.rng, model.metabolic_rate_dist[1]:model.metabolic_rate_dist[2]),
            0,
            rand(model.rng, model.max_age_dist[1]:model.max_age_dist[2]),
            rand(model.rng, model.w0_dist[1]:model.w0_dist[2]),
        )
    end
end

function agent_step!(agent, model)
    movement!(agent, model)
    replacement!(agent, model)
end

# ## Plotting & Animating

# We can plot the ABM and the sugar distribution side by side using [`abm_plot`](@ref)
# and standard Makie.jl commands like so
using InteractiveDynamics

model = sugarscape()
fig, abmstepper = abm_plot(model; resolution = (800, 600))
ax, hm = heatmap(fig[1,2], model.sugar_values; colormap=cgrad(:thermal), colorrange=(0,4))
ax.aspect = AxisAspect(1) # equal aspect ratio for heatmap
Colorbar(fig[1, 3], hm, width = 15, tellheight=false)
rowsize!(fig.layout, 1, ax.scene.px_area[].widths[2]) # Colorbar height = axis height
fig

# To animate them both however, we will use the approach Makie.jl suggests for animations,
# which is based on `Observables`. We start similarly with a call to `abm_plot`,
# but now make the plotted heatmap an obsrvable
fig, abmstepper = abm_plot(model; resolution = (800, 600))
obs_heat = Observable(model.sugar_values)
ax, hm = heatmap(fig[1,2], obs_heat; colormap=cgrad(:thermal), colorrange=(0,4))
ax.aspect = AxisAspect(1) # equal aspect ratio for heatmap
Colorbar(fig[1, 3], hm, width = 15, tellheight=false)
rowsize!(fig.layout, 1, ax.scene.px_area[].widths[2]) # Colorbar height = axis height

# and also add a title for good measure
s = Observable(0) # counter of current step, also observable
t = lift(x -> "Sugarscape, step = $x", s)
supertitle = Label(fig[0, :], t, textsize = 24, halign = :left)
fig

# We animate the evolution of both the ABM and the sugar distribution using the following
# simple loop involving the abm stepper
record(fig, "sugarvis.mp4"; framerate = 3) do io
    for j in 0:50 # = total number of frames
        recordframe!(io) # save current state
        ## This updates the abm plot:
        Agents.step!(abmstepper, model, agent_step!, model_step!, 1)
        ## This updates the heatmap:
        obs_heat[] = model.sugar_values
        ## This updates the title:
        s[] = s[] + 1
    end
end

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../sugarvis.mp4" type="video/mp4">
# </video>
# ```

# ### Distribution of wealth across individuals
# First we produce some data that include the wealth
model2 = sugarscape()
adata, _ = run!(model2, agent_step!, model_step!, 20, adata = [:wealth])
adata[1:10,:]

# And now we animate the evolution of the distribution of wealth
figure = Figure(resolution = (600, 600))
title_text = Observable("Wealth distribution of individuals\nStep 1")
figure[1, 1] = Label(figure, title_text; textsize=30, tellwidth=false)
ax = figure[2, 1] = Axis(figure; xlabel="Wealth", ylabel="Number of agents")
histdata = Observable(adata[adata.step .== 20, :wealth])
hist!(ax, histdata; bar_position=:step)
record(figure, "sugarhist.mp4", 0:20; framerate=3) do i
    histdata[] = adata[adata.step .== i, :wealth]
    title_text[] = "Wealth distribution of individuals, step = $i"
    xlims!(ax, (0, max(histdata[]...)))
    ylims!(ax, (0, 50))
end
nothing # hide

# We see that the distribution of wealth shifts from a more or less uniform distribution to a skewed distribution.

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../sugarhist.mp4" type="video/mp4">
# </video>
# ```

# ## References

# BUZING P, Eiben A & Schut M (2005) Emerging communication and cooperation in evolving agent societies. Journal of Artificial Societies and Social Simulation 8(1)2. http://jasss.soc.surrey.ac.uk/8/1/2.html.

# EPSTEIN J M (1999) Agent-Based Computational Models And Generative Social Science. Complexity 4(5), pp. 41-60.

# EPSTEIN J M & Axtell R L (1996) Growing Artificial Societies: Social Science from the Bottom Up. The MIT Press.

# FLENTGE F, Polani D & Uthmann T (2001) Modelling the emergence of possession norms using memes. Journal of Artificial Societies and Social Simulation 4(4)3. http://jasss.soc.surrey.ac.uk/4/4/3.html.


text:
importall Bukdu

type HelloController <: ApplicationController
end

layout(::Layout, body) = "layout $body"

type CustomLayout <: ApplicationLayout
end
layout(::CustomLayout, body) = "custom $body"

custom(::HelloController) = render(Text/Layout/CustomLayout, "hello")
index(::HelloController) = render(Text/CustomLayout/Layout, "hello")

Router() do
    get("/", HelloController, index)
    get("/custom", HelloController, custom)
end

using Base.Test
conn = (Router)(get, "/")
@test 200 == conn.status
@test "layout custom hello" == conn.resp_body

conn = (Router)(get, "/custom")
@test 200 == conn.status
@test "custom layout hello" == conn.resp_body

@test "custom layout hello" == render(Text/Layout/CustomLayout, "hello").resp_body
@test "layout custom hello" == render(Text/CustomLayout/Layout, "hello").resp_body
@test "custom layout hello" == render(HTML/Layout/CustomLayout, "hello")
@test "layout custom hello" == render(HTML/CustomLayout/Layout, "hello")

@test "Text{T}" == string(Text)
@test "Text/Layout" == string(Text/Layout)
@test "Text/Layout/CustomLayout" == string(Text/Layout/CustomLayout)
@test "Text/CustomLayout/Layout" == string(Text/CustomLayout/Layout)
@test "Layout/CustomLayout" == string(Layout/CustomLayout)
@test "CustomLayout/Layout" == string(CustomLayout/Layout)