text:
function concat(multiple_src::TupleOrVector{<:Memory{T,N}}, dim::Integer) where {T,N}
    nargs = length(multiple_src)
    outputdim = sum(i -> size(i, dim), multiple_src)
    dims = size(first(multiple_src))
    dst_dims = ntuple(i -> (i == dim) ? outputdim : dims[i], Val(N))
    dst_md = memorydesc(T, dst_dims, dnnl_format_any())

    return temp_primitive(
        Lib.dnnl_concat_primitive_desc_create,
        dst_md,
        nargs,
        N - dim,
        map(memorydesc, multiple_src),
        noattributes(),
        global_engine(),
    ) do p, pd
        dst_md = query_md(pd, @query(dst))
        dst = similar(first(multiple_src), T, dst_dims, dst_md)
        execute!(p, @dnnl_args dst multiple_src)
        return dst
    end
end

#####
##### Slicing
#####

# TODO: See if the #15276 style problems still exist in Julia 1.6.
#
# `map` seems to be having Julia issue #15276 is problems when keeping track of where
# we are indexing to create views.
#
# As such, we have to build this `Slicer` struct below in order to give inference
# some help.
mutable struct Slicer{T,N,A<:AbstractArray{T,N}}
    current_index::Int
    concat_dim::Int
    captured_array::A
end

function (S::Slicer{T,N})(sz) where {T,N}
    current_index = S.current_index
    range = current_index:(current_index + sz - 1)
    inds = ntuple(i -> i == S.concat_dim ? range : 1:size(S.captured_array, i), Val(N))
    S.current_index += sz
    return view(S.captured_array, inds...)
end

#####
##### concat backprop
#####

function ChainRulesCore.rrule(
    ::typeof(concat), multiple_src::TupleOrVector{T}, dim::Integer
) where {T<:Memory}
    # Capture sizes for reconstruction.
    lengths = size.(multiple_src, dim)
    dst = concat(multiple_src, dim)

    function concat_pullback(Δ)
        f = Slicer(1, dim, materialize(Δ))
        δA = map(f, lengths)
        return (ChainRulesCore.NoTangent(), δA, ChainRulesCore.NoTangent())
    end
    return dst, concat_pullback
end


text:
include("ParamEstim.jl")
using .ParamEstim

function optimize(nth_param_set::Int64)
    if !isdir(strip(MODEL_PATH, '/') * "/fitparam")
        mkdir(strip(MODEL_PATH, '/') * "/fitparam")
    end

    try
        files = readdir(
            strip(MODEL_PATH, '/') * "/fitparam/$nth_param_set"
        )
        for file in files
            if occursin(".dat",file)
                rm(
                    strip(MODEL_PATH, '/') * "/fitparam/$nth_param_set/$file"
                )
            end
        end
    catch
        mkdir(strip(MODEL_PATH, '/') * "/fitparam/$nth_param_set")
    end

    search_rgn::Matrix{Float64} = get_search_region()

    max_generation::Int64 = 10000
    n_population::Int64 = 15*size(search_rgn, 2)
    n_children::Int64 = 50
    n_gene::Int64 = size(search_rgn, 2)
    allowable_error::Float64 = 0.0

    (best_indiv, best_fitness) = ga_v2(
        nth_param_set,
        max_generation,
        n_population,
        n_children,
        n_gene,
        allowable_error
    )
end

if abspath(PROGRAM_FILE) == @__FILE__
    optimize(parse(Int64,ARGS[1]))
end

ext: jl

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim


text:
# Demands like BEVs, where there is some flexbility, but also hard constraints.

# Import modules.
using DataFrames, Dates, MAT, TimeSeries
import Base: convert

"""
Segment of a charging plan for a BEV.

# Type parameters
- `L`            : network location

# Fields:
- `location        :: L`       : network location
-  duration        :: Float64` : duration [time] of the time interval
-  load            :: Float64` : average load [energy/time] on the power system
-  chargerate      :: Float64` : average charging rate [energy/time] of the battery, accounting for the charging efficiency
- `maxchargerate   :: Float64` : maximum allowable charging rate [energy/time] for the battery [energy/time]
- `consumptionrate :: Float64` : average discharge rate [energy/time] of the battery due to driving
- `batterylevel`   :: Float64` : the battery level [energy] at the beginning of the time interval
"""
const ChargingSegment{L} = NamedTuple{(:location, :duration, :load, :chargerate, :maxchargerate, :consumptionrate, :batterylevel), Tuple{L, Float64, Float64, Float64, Float64, Float64, Float64}}

"""
Charging plan for a BEV.

# Type parameters
- `T <: TimeType`: timestamp
- `L`            : network location
"""
const ChargingPlan{T <: TimeType, L} = TimeArray{ChargingSegment{L},1,T,Vector{ChargingSegment{L}}}

"""
Convert a charging plan to a data frame.
"""
function convert(::Type{DataFrame}, plan :: ChargingPlan{T, L}) where T where L
    insertcols!(DataFrame(values(plan)), 1, :timestamp => timestamp(plan))
end

"""
Extract the locations from a charging plan.
"""
function locations(plan :: ChargingPlan{T,L}) :: TimeArray{L,1,T,Vector{L}} where T where L
    map(x -> x.location, plan)
end

"""
Extract the duration [time] of the time interval from a charging plan.
"""
function durations(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.duration, plan)
end

"""
Extract the located demand from a charging plan.
"""
function locateddemand(plan :: ChargingPlan{T,L}) :: LocatedDemand{T,L} where T where L
    map(x -> (x.location, x.load), plan)
end

"""
Extract the average system load [energy/time] from a charging plan.
"""
function loads(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.load, plan)
end

"""
Extract the average charge rate [energy/time] of the battery, accounting for the charging efficiency, from a charging plan.
"""
function chargerates(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.chargerate, plan)
end

"""
Extract the charge amount [energy] during the time interval from a charging plan.
"""
function chargeamounts(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.duration * x.chargerate, plan)
end

"""
Extract the maximum charge rate [energy/time] from a charging plan.
"""
function maxchargerates(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.maxchargerate, plan)
end

"""
Extract the average consumption rate [energy/time], the average discharge rate of the battery due to driving, from a charging plan.
"""
function consumptionrates(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.consumptionrate, plan)
end

"""
Extract the consumption amount [energy], the dicharge of the battery due to driving during the time interval, from a charging plan.
"""
function consumptionamounts(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.duration * x.consumptionrate, plan)
end

"""
Extract the battery level [energy] at the beginning of the time interval from a charging plan.
"""
function batterylevels(plan :: ChargingPlan{T,L}) :: TimeArray{Float64,1,T,Vector{Float64}} where T where L
    map(x -> x.batterylevel, plan)
end

"""
Detailed representation of a flexible, mobile demand.  Typically, this may represent single or multiple battery-electric vehicles.

This records the location of the demand and its consumption of its stored energy (i.e., its battery) as a function of time, along with physical parameters and constraints: storage capacity, charge/discharge effiency, maximum charge/discharge rate, and temporal boundary conditions.  Charging represents grid-to-vehicle transfer of energy and discharging represents vehicle-to-grid transfer.

# Type parameters
- `T <: TimeType`: timestamp
- `L`            : network location

# Fields
- `locations`    : time/event series of locations of the vehicle, and their maximum AC and DC charging rates [energy/time]
- `power`        : time/event series of consumption rate of the vehicle [energy/time]
- `capacity`     : constraint on minimum and maximum battery level [energy]
- `rate`         : constraint on minimum and maximum AC and DC charging rates of the battery [energy/time]
- `efficiency`   : efficiency of charging the battery [energy/energy], applied to computations before the charging rate
- `timeboundary` : `nothing` for cyclic boundary conditions in time, or a tuple of minimum and maximum battery level allowed at the time boundaries [energy]

# Example
```
example = BevDemand(
    TimeArray(
        [Time(0)          , Time(8)         , Time(9)              , Time(17)       , Time(18)         , Time(23,59,59)   ], # [h]
        [("Home #23", (ac=1.4, dc=0.)), ("Road #14", (ac=0., dc=0.)), ("Workplace #3", (ac=7.7, dc=0.)), ("Road #9", (ac=0., dc=0.)), ("Home #23", (ac=1.4, dc=0.)), ("Home #23", (ac=1.4, dc=0.))]  # [kW]
    ),
    TimeArray(
        [Time(0), Time(8), Time(9), Time(17), Time(18), Time(23,59,59)], # [h]
        [     0.,     10.,      0.,      11.,       0.,             0.]  # [kW]
    ),
    (min=0., max=40.), # [kWh]
    (ac=(min=0., max=20.), dc=(min=0., max=50.)), # [kW]
    (in=0.90, out=0.), # [kWh/kWh]
    nothing,
)
```
"""
struct BevDemand{T,L} <: FlexibleDemand{T,L}
    locations    :: MobileDemand{T,Tuple{L,NamedTuple{(:ac, :dc),Tuple{Float64,Float64}}}}
    power        :: TemporalDemand{T}
    capacity     :: NamedTuple{(:min, :max),Tuple{Float64,Float64}}
    rate         :: NamedTuple{(:ac, :dc),Tuple{NamedTuple{(:min, :max),Tuple{Float64,Float64}},NamedTuple{(:min, :max),Tuple{Float64,Float64}}}}
    efficiency   :: NamedTuple{(:in, :out),Tuple{Float64,Float64}}
    timeboundary :: Union{Tuple{Float64,Float64},Nothing}
end

"""
The "envelope" of minimum and maximum allowable demands at each time point, with a location for the demand.

# Arguments
- `demand :: BevDemand{T,L}`: the demand
"""
function envelope(demand :: BevDemand{T,L}) :: LocatedEnvelope{T,L} where L where T <: TimeType
    map(
        d -> (d[1][1], d[2][1][2], d[2][2][2]),
        aligntimes(
            demand.locations,
            aligntimes(
                earliestdemands(demand),
                latestdemands(demand)
            )
        )
    )
end

"""
Apply efficiency factors to relate energy at the vehicle to energy at the charger.

# Arguments
- `demand :: BevDemand{T,L}`: the demand

# Returns
- a function that converts energy at the vehicle to energy at the charger
"""
function applyefficiencies(demand :: BevDemand{T,L}) where L where T <: TimeType
    function f(x)
        if x == 0
            0.
        elseif x > 0
            x / demand.efficiency.in
        else
            x * demand.efficiency.out
        end
    end
    f
end

"""
Apply efficiency factors to relate energy at the vehicle to energy at the charger.

# Arguments
- `demand :: BevDemand{T,L}`: the demand

# Returns
- a function that converts energy at the charger to energy at the vehicle
"""
function applyefficienciesinverse(demand :: BevDemand{T,L}) where L where T <: TimeType
    function f(x)
        if x == 0
            0.
        elseif x > 0
            x * demand.efficiency.in
        else
            x / demand.efficiency.out
        end
    end
    f
end

"""
The demand if charging takes place as early as possible.

# Arguments
- `demand :: BevDemand{T,L}`: the demand

# Returns
the demand over time
"""
function earliestdemands(demand :: BevDemand{T,L}) :: LocatedDemand{T,L} where L where T <: TimeType
    if demand.timeboundary == nothing
        b = earliestdemands(demand, demand.capacity.max)[2]
        earliestdemands(demand, b)[1]
    else
        earliestdemands(demand, demand.timeboundary[2])[1]
    end
end

"""
The demand if charging takes place as early as possible.

# Arguments
- `demand  :: BevDemand{T,L}`: the demand
- `initial :: Float64`       : the battery level at the first time

# Returns
- the demand over time
- the battery level at the last time
"""
function earliestdemands(demand :: BevDemand{T,L}, initial :: Float64) :: Tuple{LocatedDemand{T,L},Float64} where L where T <: TimeType
    resolution = 1 / 60 / 60 / 1000 # Resolve to the nearest millisecond.
    eff = applyefficiencies(demand)
    onehour = Time(1) - Time(0)
    x = aligntimes(demand.locations, demand.power)
    xt = timestamp(x)
    xv = values(x)
    zt = Array{T,1}()
    zv = Array{Float64,1}()
    b = initial
    ix = 2
    tnow = xt[ix-1]
    while ix <= length(x)
        tnext = (xt[ix] - tnow) / onehour
        ((_, (chargingac, chargingdc)), consumption) = xv[ix-1]
        chargingac1 = min(b < demand.capacity.max ? chargingac : 0, demand.rate.ac.max)
        chargingdc1 = min(b < demand.capacity.max ? chargingdc : 0, demand.rate.dc.max)
        charging1 = max(chargingac1, chargingdc1)
        net = charging1 - consumption
        tcrit = net > 0 ? (demand.capacity.max - b) / net : Inf
        if tcrit < resolution
            b = demand.capacity.max
        else
            push!(zt, tnow)
            push!(zv, eff(charging1))
            if tcrit < tnext
                b = b + tcrit * net
                tnow = addhours(tnow, tcrit)
            else
                b = b + tnext * net
                tnow = xt[ix]
                ix = ix + 1
            end
        end
    end
    push!(zt, tnow)
    push!(zv, NaN)
    (
        aligntimes(map(v -> v[1], demand.locations), TimeArray(zt, zv)),
        b
    )
end

"""
The demand if charging takes place as late as possible.

# Arguments
- `demand :: BevDemand{T,L}`: the demand

# Returns
- the demand over time
"""
function latestdemands(demand :: BevDemand{T,L}) :: LocatedDemand{T,L} where L where T <: TimeType
    if demand.timeboundary == nothing
        b = latestdemands(demand, demand.capacity.min)[2]
        latestdemands(demand, b)[1]
    else
        latestdemands(demand, demand.timeboundary[2])[1]
    end
end

"""
The demand if charging takes place as late as possible.

# Arguments
- `demand :: BevDemand{T,L}`: the demand
- `final  :: Float64`       : the battery level at the last time

# Returns
- the demand over time
- the battery level at the first time
"""
function latestdemands(demand :: BevDemand{T,L}, final :: Float64) :: Tuple{LocatedDemand{T,L},Float64} where L where T <: TimeType
    resolution = 1 / 60 / 60 / 1000 # Resolve to the nearest millisecond.
    eff = applyefficiencies(demand)
    onehour = Time(1) - Time(0)
    x = aligntimes(demand.locations, demand.power)
    xt = timestamp(x)
    xv = values(x)
    zt = Array{T,1}()
    zv = Array{Float64,1}()
    b = final
    ix = length(x) - 1
    tnow = xt[ix+1]
    push!(zt, tnow)
    push!(zv, NaN)
    while ix >= 1
        tnext = (tnow - xt[ix]) / onehour
        ((_, (chargingac, chargingdc)), consumption) = xv[ix]
        chargingac1 = min(b > demand.capacity.min ? chargingac : 0, demand.rate.ac.max)
        chargingdc1 = min(b > demand.capacity.min ? chargingdc : 0, demand.rate.dc.max)
        charging1 = max(chargingac1, chargingdc1)
        net = charging1 - consumption
        tcrit = net > 0 ? (b - demand.capacity.min) / net : Inf
        if tcrit < resolution
            b = demand.capacity.min
        else
            if tcrit < tnext
                b = b - tcrit * net
                tnow = addhours(tnow, - tcrit)
            else
                b = b - tnext * net
                tnow = xt[ix]
                ix = ix - 1
            end
            push!(zt, tnow)
            push!(zv, eff(charging1))
        end
    end
    zt = reverse(zt)
    zv = reverse(zv)
    (
        aligntimes(map(v -> v[1], demand.locations), TimeArray(zt, zv)),
        b
    )
end

"""
Procedural approach for BEV Greedy charging scenario, maximizing the BEV's battery level.

# Arguments
- `demand :: BevDemand{T,L}`: the BEV demand
- `prices :: TimeArray{T}`  : the electricity prices

# Returns
- `ChargingPlan{T,L}` : the charging plan
"""
function greedydemands(demand :: BevDemand{T,L}) where L where T <: TimeType
    pricing = map(v -> 1., demand.power)
    greedydemands(demand, pricing)
end

"""
Procedural approach for BEV Greedy charging scenario, maximizing the BEV's battery level.

# Arguments
- `demand :: BevDemand{T,L}`: the BEV demand
- `prices :: TimeArray{T}`  : the electricity prices

# Returns
- `ChargingPlan{T,L}` : the charging plan
"""
function greedydemands(demand :: BevDemand{T,L}, prices :: TimeArray{Float64,1,T,Vector{Float64}}) :: ChargingPlan{T,L} where L where T <: TimeType
    @warn "`PowerSystems.greedydemands` fails some tests."
    eff = applyefficiencies(demand)

    onehour = Time(1) - Time(0)
    eff = applyefficiencies(demand)
    x = aligntimes(aligntimes(demand.locations, demand.power), prices)
    xt = timestamp(x)
    xv = values(x)

    NT = length(x)
    NP = NT - 1

    hour = map(t -> t.instant / onehour, xt)
    location = map(v -> v[1][1][1], xv)
    duration = (xt[2:NT] - xt[1:NP]) / onehour
    chargemin = duration .* map(v -> min(v[1][1][2].ac, demand.rate.ac.min), xv[1:NP])
    chargemax = duration .* map(v -> min(v[1][1][2].ac, demand.rate.ac.max), xv[1:NP])
    consumption = duration .* map(v -> v[1][2], xv[1:NP])

    #Find value of greatest consumption in back-to-back time interval
    consumpMax = [0.0]
    consumpCompare = [0.0]
    for i in 1:NP
        if i == 1
            if consumption[1] > 0
                consumpCompare[1] += consumption[1]
            end
        else
            if consumption[i] > 0 && consumption[i-1] == 0
                consumpCompare[1] += consumption[i]
            elseif consumption[i] > 0 && consumption[i-1] > 0
                consumpCompare[1] += consumption[i]
            elseif consumption[i] == 0
                consumpMax[1] = max(consumpMax[1], consumpCompare[1])
                consumpCompare[1] = 0
            end
        end
    end

    ##Attempt greedy charging procedure with inital battery level at maximum capacity

    #Initialize battery levels
    battery = [demand.capacity.max] #Initial and final battery values?
    charge = []

    for i in 1:NP
        if (battery[i] < demand.capacity.max) && (chargemax[i] > 0)
            if (battery[i] + chargemax[i]) > demand.capacity.max
                append!(battery, demand.capacity.max)
                append!(charge, (demand.capacity.max - battery[i]))
            else
                append!(battery, battery[i] + chargemax[i])
                append!(charge, chargemax[i])
            end
        elseif (battery[i] == demand.capacity.max) && (consumption[i] == 0)
            append!(battery, battery[i])
            append!(charge, 0.0)
        elseif (chargemax[i] == 0.0) && (consumption[i] == 0)
            append!(battery, battery[i])
            append!(charge, 0.0)
        elseif consumption[i] > 0
            append!(battery, battery[i] - consumption[i])
            append!(charge, 0.0)
        end
    end

    if length(battery[battery .< 0]) > 0
        battery = []
        for i in 1:NT
            append!(battery, NaN)
        end

        charge = []
        for i in 1:NP
            append!(charge, NaN)
        end

    elseif (sum(consumption) > sum(chargemax)) && (length(battery[battery .< 0]) == 0)
        battery = []
        for i in 1:NT
            append!(battery, NaN)
        end

        charge = []
        for i in 1:NP
            append!(charge, NaN)
        end

#   elseif (battery[96] < demand.capacity.max) && length(battery[battery .< 0]) == 0 && sum(consumption) <= sum(chargemax)
    elseif (battery[NT] < demand.capacity.max) && length(battery[battery .< 0]) == 0 && sum(consumption) <= sum(chargemax)
        #Create curve ignoring maximum battery capacity
        battery = [demand.capacity.max]
        charge = []
        for i in 1:NP
            if chargemax[i] > 0
                append!(battery, battery[i] + chargemax[i])
                append!(charge, chargemax[i])
            elseif consumption[i] > 0
                append!(battery, battery[i] - consumption[i])
                append!(charge, 0.0)
            elseif (chargemax[i] == 0.0) && (consumption[i] == 0)
                append!(battery, battery[i])
                append!(charge, 0.0)
            end
        end

        #Check if there is a starting charge value for which the ending charge value can be equal.
        if consumpMax[1] > (demand.capacity.max - demand.capacity.min)
            battery = []
            for i in 1:NT
                append!(battery, NaN)
            end

            charge = []
            for i in 1:NP
                append!(charge, NaN)
            end

        elseif battery[NT] >= battery[1] && consumpMax[1] <= (demand.capacity.max - demand.capacity.min)
            #Shift data down such that the maximum battery level lies along the
            #maximum battery capacity line.
            shift1 = minimum(battery) - demand.capacity.min
            for i in 1:NT
                battery[i] -= shift1
            end
            if(maximum(battery) > demand.capacity.max)
                shift = maximum(battery) - demand.capacity.max
            else
                shift = 0
            end
            battery = [battery[NT] - shift]
            charge = []

            #Apply greedy charging procedure with new inital battery level
            for i in 1:NP
                if (battery[i] < demand.capacity.max) && (chargemax[i] > 0)
                    if (battery[i] + chargemax[i]) > demand.capacity.max
                        append!(battery, demand.capacity.max)
                        append!(charge, (demand.capacity.max - battery[i]))
                    else
                        append!(battery, battery[i] + chargemax[i])
                        append!(charge, chargemax[i])
                    end
                elseif (battery[i] == demand.capacity.max) && (consumption[i] == 0)
                    append!(battery, battery[i])
                    append!(charge, 0.0)
                elseif (chargemax[i] == 0.0) && (consumption[i] == 0)
                    append!(battery, battery[i])
                    append!(charge, 0.0)
                elseif consumption[i] > 0
                    append!(battery, battery[i] - consumption[i])
                    append!(charge, 0.0)
                end
            end
        end
    end

    convert(Array{Float64}, battery)
    convert(Array{Float64}, charge)
    #Contains optimization charging results with charging rate from charger during each time interval

    TimeArray(
        xt,
        vcat(
            [(
                location        = location[i]                 ,
                duration        = duration[i]                 ,
                load            = eff(charge[i] / duration[i]),
                chargerate      = charge[i] / duration[i]     ,
                maxchargerate   = chargemax[i]                ,
                consumptionrate = consumption[i]              ,
                batterylevel    = battery[i]                  ,
            ) for i in 1:NP],
            (
                location        = location[NT]                ,
                duration        = NaN                         ,
                load            = NaN                         ,
                chargerate      = NaN                         ,
                maxchargerate   = NaN                         ,
                consumptionrate = NaN                         ,
                batterylevel    = battery[NT]                 ,
            )
        )
    )

end

"""
Check the validity of a charging plan.

# Arguments
- `demand    :: BevDemand{T,L}`        : the BEV demand demand
- `charging  :: LocatedDemand{T,L}`    : the charging plan
- `tolerance :: Float64`               : tolerance for testing limits
- `message   :: Union{String,Nothing}` : the prefix of the warning messages, or `nothing` if no messages are to be logged

# Returns
- whether charging does not differ from consumption
- whether charging rates do not exceed limits
- whether battery levels do not exceed limits
"""
function verify(demand :: BevDemand{T,L}, charging :: LocatedDemand{T,L}; tolerance :: Float64 = 1e-5, message :: Union{String,Nothing} = nothing) :: NamedTuple{(:balance, :rates, :battery),Tuple{Bool,Bool,Bool}} where L where T <: TimeType
    verbose = !isnothing(message)
    the_shortfall = shortfall(demand, charging)
    check_shortfall = abs(the_shortfall) <= tolerance
    if verbose && !check_shortfall
        @warn string(message, " has charging shortfall of ", the_shortfall, " kWh.")
    end
    check_rates     = verifylimits(demand, charging, tolerance=tolerance)
    if verbose && !check_rates
        @warn string(message, " violates charging limits.")
    end
    check_battery   = verifybattery(demand, charging, tolerance=tolerance)
    if verbose && !check_battery
        @warn string(message, " violates battery limits.")
    end
    (
        balance = check_shortfall,
        rates   = check_rates    ,
        battery = check_battery  ,
    )
end

"""
Compute the discrepancy between the BEV constraints and the charging plan.

# Arguments
- `demand   :: BevDemand{T,L}`     : the BEV demand demand
- `charging :: LocatedDemand{T,L}` : the charging plan

# Returns
- the shortfall of the charging plan's meeting of the demand demand
"""
function shortfall(demand :: BevDemand{T,L}, charging :: LocatedDemand{T,L}) :: Float64 where L where T <: TimeType
    onehour = Time(1) - Time(0)
    function durations(x)
        xt = timestamp(x)
        (xt[2:end] - xt[1:end-1]) / onehour
    end
    function rates(x)
        xv = map(y -> y[end], values(x))
        xv[1:end-1]
    end
    function levels(x)
        durations(x) .* rates(x)
    end
    function total(x)
        sum(levels(x))
    end
    total(demand.power) / demand.efficiency.in - total(charging)
end

"""
Verify charging limits.

# Arguments
- `demand    :: BevDemand{T,L}`     : the BEV demand demand
- `charging  :: LocatedDemand{T,L}` : the charging plan
- `tolerance :: Float64`            : tolerance for testing limits

# Returns
- whether the charging limits are satisfied
"""
function verifylimits(demand :: BevDemand{T,L}, charging :: LocatedDemand{T,L}; tolerance :: Float64 = 1e-5) :: Bool where L where T <: TimeType
    onehour = Time(1) - Time(0)
    eff = applyefficienciesinverse(demand)
    x = aligntimes(demand.locations, charging)
    xv = values(x[1:end-1])
    powers = eff.(map(v -> v[2][2], xv))
    limits = map(v -> max(min(v[1][2].ac, demand.rate.ac.max), min(v[1][2].dc, demand.rate.dc.max)), xv)
    all(powers .<= limits .+ tolerance)
end

"""
Compute battery levels.

# Arguments
- `demand    :: BevDemand{T,L}`     : the BEV demand demand
- `charging  :: LocatedDemand{T,L}` : the charging plan

# Returns
- the battery levels
"""
function chargelevels(demand :: BevDemand{T,L}, charging :: LocatedDemand{T,L}) :: TemporalDemand{T} where L where T <: TimeType
    onehour = Time(1) - Time(0)
    eff = applyefficienciesinverse(demand)
    x = aligntimes(demand.power, charging)
    xt = timestamp(x)
    durations = (xt[2:end] - xt[1:end-1]) ./ onehour
    xv = values(x[1:end-1])
    nets = (eff.(map(v -> v[2][2], xv)) .- map(v -> v[1], xv)) .* durations
    nets = prepend!(nets, [0.])
    nets = cumsum(nets)
    nets = nets .- minimum(nets) .+ demand.capacity.min
    TimeArray(xt, nets)
end

"""
Verify battery levels.

# Arguments
- `demand    :: BevDemand{T,L}`     : the BEV demand demand
- `charging  :: LocatedDemand{T,L}` : the charging plan
- `tolerance :: Float64`            : tolerance for testing limits

# Returns
- whether the battery level constraints are satisfied
"""
function verifybattery(demand :: BevDemand{T,L}, charging :: LocatedDemand{T,L}; tolerance :: Float64 = 1e-5) :: Bool where L where T <: TimeType
    levels = values(chargelevels(demand, charging))
    maxokay = maximum(levels) <= demand.capacity.max + tolerance
    if isnothing(demand.timeboundary)
        boundaryokay = abs(levels[1] - levels[end]) <= tolerance
    else
        boundaryokay  = demand.timeboundary[1] - tolerance <= levels[1  ] <= demand.timeboundary[2] + tolerance
        boundaryokay &= demand.timeboundary[1] - tolerance <= levels[end] <= demand.timeboundary[2] + tolerance
    end
    maxokay && boundaryokay
end

"""
Function to populate a BevDemand struct. Receives location of ".mat" file as string input
and index of Bev entry to read in ".mat".
The ".mat" file is output from EVI-Pro tool.

# Arguments
- `data_location :: String`: the path of the file containing the BEV data
- 'index :: Int64': the Bev entry to read in ".mat" file
"""
function populate_BEV_demand(file :: MAT.MAT_v5.Matlabv5File, index :: Int64)
    workaround1 = true

    if index < 1
        print("index out of range")
        return []
    else
        @time entry = read(file, "bev$index")
        num = length(entry["locations"]["bus_id"]);

        if workaround1
            num = parse(Int64, string(num))
        end

        #             Locations
        loc_tuples = Array{Tuple{String,NamedTuple{(:ac, :dc),Tuple{Float64,Float64}}}}(undef,num);

        for j in range(1,num,step=1)
            loc = string("bus id #",entry["locations"]["bus_id"][j])
            a_ch = parse(Float64, string(entry["locations"]["max_ac_charging_rate"][j]))
            d_ch = parse(Float64, string(entry["locations"]["max_dc_charging_rate"][j]))
            loc_tuples[j] = (loc,(ac=a_ch,dc=d_ch))
        end

        num = length(entry["locations"]["time_stamp"])

        temp_time_loc = [fldmod(Int(ceil(entry["consumptions"]["time_stamp"][j]*24*60)),60) for j in
                range(1,num, step = 1)];
        time_stamp_loc = [Time(temp_time_loc[j][1],temp_time_loc[j][2]) for j in range(1,num, step = 1)];
        locations = TimeArray(time_stamp_loc, loc_tuples);

        #Processing the data for consumption
        num = length(entry["consumptions"]["time_stamp"]);
        temp_time = [fldmod(Int(ceil(entry["consumptions"]["time_stamp"][j]*24*60)),60) for j in
                range(1,num, step = 1)];
        time_stamp = [Time(temp_time[j][1],temp_time[j][2]) for j in range(1,num, step = 1)];
        cons_val = [entry["consumptions"]["consumptions_rates"][k] for k in (1:num)];
        consumptions = TimeArray(time_stamp,cons_val);

        # Adding remaining attributes
        capacity = (min=entry["storagemin"],max =entry["storagemax"])
        rate = (ac=(min= 0.0, max = entry["AC_chargeratemax"]),
            dc=(min=0.0,max=entry["DC_chargeratemax"]))
        efficiency = (in=entry["chargeEfficiency"],out=entry["dischargeEfficiency"])
        timeboundary = nothing

#             Adding created BevDemand struct in array
        populated_BEV_demand = BevDemand(locations, consumptions, capacity, rate,efficiency, timeboundary)
    end

    return populated_BEV_demand;
end

"""
Function to populate a BevDemand struct. Receives location of ".mat" file as string input
and range of Bev entries to read in ".mat".
The ".mat" file is output from EVI-Pro tool.

# Arguments
- `data_location :: String`: the path of the file containing the BEV data
- 'start_i :: Int64': index of first Bev entry to read in ".mat" file
- 'end_i :: Int64': index of last Bev entry to read in ".mat" file
"""
function populate_BEV_demand(data_location :: String, start_i:: Int64, end_i:: Int64) :: Array{BevDemand{Time,String}}
    file = matopen(data_location)

    #Getting file indices from file name
    pieces = split(data_location, '_', keepempty=false)
    num_pieces = length(pieces)
    start_index = parse(Int64,pieces[num_pieces-1])
    end_pieces = split(pieces[num_pieces],'.', keepempty=false)
    end_index = parse(Int64,end_pieces[1])

    num_el = end_index - start_index + 1

    workaround1 = true

    if start_i < 1 || end_i > (start_i + num_el-1) || start_i < start_index || end_i > end_index
        print("indices out of range")
        return []
    else
        size_array = end_i - start_i + 1;
        populated_BEV_demand = Array{BevDemand}(undef,size_array)

        b_index = 1;
        for i in range(start_i,end_i, step = 1)
            populated_BEV_demand[b_index] = populate_BEV_demand(file, i)
            b_index = b_index + 1;
        end
    end

    return populated_BEV_demand;
end

"""
Function to populate an array of BevDemand structs. Receives location of ".mat" file as string input.
The ".mat" file is output from EVI-Pro tool.

# Arguments
- `data_location :: String`: the path of the file containing the BEV data
"""
function populate_BEV_demand(data_location :: String) :: Array{BevDemand{Time,String}}
    file = matopen(data_location)

    #Getting file indices from file name
    pieces = split(data_location, '_', keepempty=false)
    num_pieces = length(pieces)
    start_index = parse(Int64,pieces[num_pieces-1])
    end_pieces = split(pieces[num_pieces],'.', keepempty=false)
    end_index = parse(Int64,end_pieces[1])

    num_el = end_index - start_index + 1

    populated_BEV_demand = Array{BevDemand}(undef,num_el)

    i = 1
    for index in range(start_index, end_index, step = 1)
#         index = parse(Int64,(replace(var, "bev" => "")))
        populated_BEV_demand[i] = populate_BEV_demand(file, index)
        i = i+1
    end

    return populated_BEV_demand
end


text:
default = LocalCSRMatrix{Float32, UInt32}()
@test 0 == numRows(default)
@test 0 == numCols(default)
@test_throws BoundsError getRowView(default, 1)
@test_throws BoundsError getRowView(default, 5)

rawVals = Float32[5, 8, 6, 2, 1, 6]
rawCols =  UInt32[2, 4, 5, 2, 3, 1]
mat = LocalCSRMatrix(4, 5, rawVals, UInt32[1, 2, 4, 6, 7], rawCols)
@test isa(mat, LocalCSRMatrix{Float32, UInt32})
@test 4 == numRows(mat)
@test 5 == numCols(mat)
@test Float32[5] == vals(getRowView(mat, 1))
@test  UInt32[2] == cols(getRowView(mat, 1))
@test Float32[8, 6] == vals(getRowView(mat, 2))
@test  UInt32[4, 5] == cols(getRowView(mat, 2))
@test Float32[2, 1] == vals(getRowView(mat, 3))
@test  UInt32[2, 3] == cols(getRowView(mat, 3))
@test Float32[6] == vals(getRowView(mat, 4))
@test  UInt32[1] == cols(getRowView(mat, 4))


text:
function build_model_generation_report(problem_dictionary::Dict{String,Any})::VLProgramComponent

    # initialize -
    buffer = Array{String,1}()
    program_component = VLProgramComponent()

    # get metadata -
    metadata_dictionary = problem_dictionary["author_metadata_dictionary"]

    # get metadata stuff -
    author_first_name = metadata_dictionary["author_first_name"]
    author_last_name = metadata_dictionary["author_last_name"]
    author_email_address = metadata_dictionary["author_email_address"]
    author_address=metadata_dictionary["author_address"]
    job_description=metadata_dictionary["job_description"]
    job_version=metadata_dictionary["job_version"]

    # get some path information -
    path_to_model_file = problem_dictionary["path_to_model_file"]
    path_to_output_dir = problem_dictionary["path_to_output_dir"]
    elapsed_time = problem_dictionary["elapsed_time"]

    # add a metadata section -
    +(buffer, "// --------------------------------------------------------------------------- //\n")
    +(buffer, "// Metadata - \n")
    +(buffer, "// author: ")
    +(buffer, "$(author_first_name) $(author_last_name)\n")
    +(buffer, "// email: $(author_email_address)\n")
    +(buffer, "// address: $(author_address)\n")
    +(buffer, "// description: $(job_description)\n")
    +(buffer, "// version: $(job_version)\n")
    +(buffer, "//\n");
    +(buffer, "// path to model file: $(path_to_model_file)\n")
    +(buffer, "// path to output: $(path_to_output_dir)\n")
    +(buffer, "// elapsed time: $(elapsed_time) s\n")
    +(buffer, "// --------------------------------------------------------------------------- //\n")
    +(buffer, "\n");

    # add reaction header -
    +(buffer, "// --------------------------------------------------------------------------- //\n")
    +(buffer, "// Reaction list - \n")
    +(buffer, "// record: index name [ec number(s)] left<=>right is_reversible {true|false}\n")
    +(buffer, "// --------------------------------------------------------------------------- //\n")

    # grab the reaction ordering array -
    reaction_order_array = problem_dictionary["metabolic_reaction_order_array"]
    reaction_dictionary = problem_dictionary["metabolic_reaction_dictionary"]
    for (index, reaction_key) in enumerate(reaction_order_array)

        # check do we have this key?
        if (haskey(reaction_dictionary,reaction_key) == true)

            # VLMetabolicReaction -
            # record::String
            # reaction_name::String
            # ec_number::String
            # left_phrase::String
            # right_phrase::String
            # reverse::String
            # forward::String

            # grab this reaction object -
            reaction_object = reaction_dictionary[reaction_key]
            +(buffer, "$(index) ")
            +(buffer, "$(reaction_object.reaction_name) ")
            +(buffer, "$(reaction_object.ec_number) ")
            +(buffer, "$(reaction_object.left_phrase)")
            +(buffer, "<=>")
            +(buffer, "$(reaction_object.right_phrase)")

            # direction -
            reverse_flag = reaction_object.reversible
            +(buffer," ")
            +(buffer, reverse_flag)
            +(buffer, "\n")

        else
            error("$(reaction_key) not found in the problem dictionary")
        end
    end

    # add metabolite header -
    +(buffer,"\n")
    +(buffer, "// --------------------------------------------------------------------------- //\n")
    +(buffer, "// Metabolite list - \n")
    +(buffer, "// record: index name type \n")
    +(buffer, "// --------------------------------------------------------------------------- //\n")

    # metabolites
    list_of_species_objects = problem_dictionary["list_of_species"]
    for (index, species_object) in enumerate(list_of_species_objects)

        # grab the symbol -
        species_symbol = species_object.species_symbol
        species_type = string(species_object.species_type)

        # write -
        +(buffer, "$(index) ")
        +(buffer, "$(species_symbol) ")
        +(buffer, "$(species_type)")
        +(buffer, "\n")

    end

    # collapse -
    flat_buffer = ""
    [flat_buffer *= line for line in buffer]

    # add stuff to component -
    program_component.filename = "Report.log"
    program_component.buffer = flat_buffer
    program_component.type = :buffer

    # return -
    return program_component
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
# test the three functions in utilities
using Test
include("utilities.jl")
@test checkmem() > 0
@test checkcpu() > 0
@test success(callout("process finished"))


text:
module p2

# Return boundary conditions values and indices, ϕb and β, for problem 2,
# with N subdivisions
function boundaryconditions2(N::Integer)
    ϕb = zeros(N+1,N+1)                 # boundary conditions
    β  = zeros(Int,N+1,N+1)             # boundary condition locations

    # Conditions on the perimeter
    ϕb[:,1] = ϕb[:,N+1] = ϕb[1,:] = ϕb[N+1,:] = 0.0
    β[:,1]  = β[:,N+1]  = β[1,:]  = β[N+1,:]  = 1

    # Conditions on the interior
    # xmin < x < xmax  →  N*xmin + 1 < i < N*xmax + 1
    xmin = 0.5
    xmax = 0.75
    ymin = 0.625
    ymax = 0.875

    xs = int([ceil(N*xmin+1):floor(N*xmax+1)])
    ys = int([ceil(N*ymin+1):floor(N*ymax+1)])

    ϕb[xs,ys] = 0.0
    β[xs,ys]  = 1

    return ϕb, β
end

# Create the inverse iteration matrix, A = Δ - μI, of size N
function invit(Δ::Array, μ::Real, β::Array)
    return Δ - μ .* eye(size(Δ)[1]) .* ((1.-β)*(1.-β)')
end

end


text:
function qsort(n, from, to)
    if from >= to
        return
    end

    pivot = n[from]

    j = from
    for i = (from + 1):to
        if n[i] < pivot
            t = n[j]
            n[j] = n[i]
            n[i] = t

            j += 1
        end

        i += 1
    end

    n[j] = pivot

    qsort(n, from, j-1)
    qsort(n, j+1, to)

    return n
end

function quicksort(n)
    return qsort(n, 1, length(n))
end

println(quicksort([5, 6, 3, 2, 8]))


text:
using HSA
using FactCheck

facts("AQL Packets") do
    # Little Endian!
    dispatch_bytes = UInt8[
        HSA.PacketTypeKernelDispatch,
        (0x01 << 0) | # Barrier Bit
        (HSA.FenceScopeAgent << 1) | # Acquire Scope
        (HSA.FenceScopeSystem << 3), # Release Scope
        # Dispatch Packet
        0x03, 0x00, # setup (dimensions)
        0x0A, 0x00, # wg_s_x
        0x09, 0x00, # wg_s_y
        0x08, 0x00, # wg_s_z
        0x00, 0x00, # reserved
        0x14, 0x00, 0x00, 0x00, # gr_s_x
        0x12, 0x00, 0x00, 0x00, # gr_s_y
        0x10, 0x00, 0x00, 0x00, # gr_s_z
        0x20, 0x00, 0x00, 0x00, # pr_seg_size
        0x30, 0x00, 0x00, 0x00, # gr_seg_size
        0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # kern_obj_addr
        0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # kern_arg_addr
        0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # reserved
        0x10, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # completion_signal
    ]
    dispatch_ptr = convert(Ptr{Void}, pointer(dispatch_bytes))

    dispatch_1d_bytes = copy(dispatch_bytes)
    dispatch_1d_bytes[3] = 0x01
    dispatch_1d_bytes[7] = 0x01
    dispatch_1d_bytes[9] = 0x01

    dispatch_1d_ptr = convert(Ptr{Void}, pointer(dispatch_1d_bytes))

    agent_bytes = UInt8[
        HSA.PacketTypeAgentDispatch,
        (0x00 << 0) | # Barrier Bit
        (HSA.FenceScopeNone << 1) | # Acquire Scope
        (HSA.FenceScopeNone << 3), # Release Scope
        # Agent Dispatch Packet
        0xFF, 0xFF, # type
        0x00, 0x00, 0x00, 0x00, # reserved
        0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # return_address
        0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # arg1
        0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # arg2
        0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # arg3
        0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # arg4
        0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # reserved
        0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # completion_signal
    ]
    agent_ptr = convert(Ptr{Void}, pointer(agent_bytes))

    barrier_bytes = UInt8[
        HSA.PacketTypeBarrierAnd,
        (0x00 << 0) | # Barrier Bit
        (HSA.FenceScopeNone << 1) | # Acquire Scope
        (HSA.FenceScopeNone << 3), # Release Scope
        # Barrier Packet
        0x00, 0x00, # UInt16 reserved
        0x00, 0x00, 0x00, 0x00, # UInt32 reserved
        0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # dep1
        0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # dep2
        0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # dep3
        0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # dep4
        0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # dep5
        0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # reserved
        0x00, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, # completion_signal
    ]
    barrier_ptr = convert(Ptr{Void}, pointer(barrier_bytes))

    context("PacketHeader") do
        context("can be loaded") do
            header = HSA.unsafe_convert(HSA.PacketHeader, dispatch_ptr)
            @fact header.typ --> HSA.PacketTypeKernelDispatch
            @fact header.barrier --> true
            @fact header.acquire_fence_scope --> HSA.FenceScopeAgent
            @fact header.release_fence_scope --> HSA.FenceScopeSystem
        end

        context("can be stored") do
            hdr_bytes = Array(UInt8, 2)
            hdr_ptr = convert(Ptr{Void}, pointer(hdr_bytes))

            hdr = HSA.unsafe_convert(HSA.PacketHeader, dispatch_ptr)
            HSA.unsafe_store!(hdr_ptr, hdr)

            @fact hdr_bytes --> dispatch_bytes[1:2]
        end

        context("can be created in various ways") do
            ph1 = HSA.PacketHeader(
                1, # type
                barrier = true,
                acquire = HSA.FenceScopeSystem,
                release = HSA.FenceScopeNone
            )

            @fact ph1.typ --> 1
            @fact ph1.barrier --> true
            @fact ph1.acquire_fence_scope --> HSA.FenceScopeSystem
            @fact ph1.release_fence_scope --> HSA.FenceScopeNone

            ph2 = HSA.PacketHeader(
                HSA.PacketTypeInvalid
            )

            @fact ph2.typ --> HSA.PacketTypeInvalid
            @fact ph2.barrier --> false
            @fact ph2.acquire_fence_scope --> HSA.FenceScopeNone
            @fact ph2.release_fence_scope --> HSA.FenceScopeNone
        end

        context("Can be copied") do
            hdr1 = PacketHeader(HSA.PacketTypeInvalid)

            hdr2 = copy(hdr1)

            @fact hdr1 --> hdr2
            @fact object_id(hdr1) --> not(object_id(hdr2))
        end

        context("Can be compared for equality") do
            h1 = PacketHeader(HSA.PacketTypeKernelDispatch)
            h2 = PacketHeader(HSA.PacketTypeKernelDispatch)

            @fact isequal(h1, h2) --> true
        end
    end

    context("VendorPackets") do

    end

    context("KernelDispatchPackets") do
        context("can be loaded") do
            dp = HSA.unsafe_convert(HSA.AQLPacket, dispatch_ptr)

            @fact isa(dp, HSA.KernelDispatchPacket) --> true

            @fact dp.header.typ --> HSA.PacketTypeKernelDispatch
            @fact dp.header.barrier --> true
            @fact dp.header.acquire_fence_scope --> HSA.FenceScopeAgent
            @fact dp.header.release_fence_scope --> HSA.FenceScopeSystem

            @fact dp.dimensions --> 0x0003
            @fact dp.workgroup_size_x --> 0x000A
            @fact dp.workgroup_size_y --> 0x0009
            @fact dp.workgroup_size_z --> 0x0008
            @fact dp.grid_size_x --> 0x00000014
            @fact dp.grid_size_y --> 0x00000012
            @fact dp.grid_size_z --> 0x00000010
            @fact dp.private_segment_size --> 0x00000020
            @fact dp.group_segment_size --> 0x00000030
            @fact dp.kernel_object --> 0x0000000000000001
            @fact dp.kernarg_address --> 0x0000000000000002
            @fact dp.completion_signal.handle --> 0x0000000000000010

            dp1d = HSA.unsafe_convert(HSA.AQLPacket, dispatch_1d_ptr)
            @fact dp1d.dimensions --> 0x0001
        end

        context("can be stored") do
            pkt_bytes = Array(UInt8, 64)
            pkt_ptr = convert(Ptr{Void}, pointer(pkt_bytes))

            dp = HSA.unsafe_convert(HSA.AQLPacket, dispatch_ptr)
            HSA.unsafe_store!(pkt_ptr, dp)

            @fact pkt_bytes --> dispatch_bytes
        end

        context("can be created in various ways") do
            @fact_throws HSA.KernelDispatchPacket(0)
            @fact_throws HSA.KernelDispatchPacket(4, 1,2,3,4,5,6,7,8)

            dp1 = HSA.KernelDispatchPacket(0x0, (2,))
            @fact dp1.dimensions --> 1
            @fact dp1.workgroup_size_x --> 1
            @fact dp1.workgroup_size_y --> 1
            @fact dp1.workgroup_size_z --> 1
            @fact dp1.grid_size_x --> 2
            @fact dp1.grid_size_y --> 1
            @fact dp1.grid_size_z --> 1
            @fact dp1.private_segment_size --> 0
            @fact dp1.group_segment_size --> 0
            @fact dp1.kernel_object --> 0
            @fact dp1.kernarg_address --> 0
            @fact dp1.completion_signal.handle --> 0
        end
    end

    context("InvalidPackets") do
        in_disp_bytes = copy(dispatch_bytes)
        in_disp_bytes[1] = convert(UInt8, HSA.PacketTypeInvalid)
        in_disp_ptr = convert(Ptr{Void}, pointer(in_disp_bytes))

        context("Can be loaded") do
            p = HSA.unsafe_convert(HSA.AQLPacket, in_disp_ptr)

            @fact isa(p, InvalidPacket) --> true
            @fact p.bytes --> in_disp_bytes
        end

        context("Can be reinterpreted as a KernelDispatchPacket") do
            p = HSA.unsafe_convert(HSA.AQLPacket, in_disp_ptr)
            id = convert(KernelDispatchPacket, p)

            @fact id.header.typ --> HSA.PacketTypeInvalid
        end
    end

    context("AgentDispatchPackets") do
        context("can be loaded") do
            ad = HSA.unsafe_convert(HSA.AQLPacket, agent_ptr)

            @fact isa(ad, HSA.AgentDispatchPacket) --> true

            @fact ad.header.typ --> HSA.PacketTypeAgentDispatch
            @fact ad.header.barrier --> false
            @fact ad.header.acquire_fence_scope --> HSA.FenceScopeNone
            @fact ad.header.release_fence_scope --> HSA.FenceScopeNone

            @fact ad.typ --> 0xFFFF
            @fact ad.return_address --> 0x0000000000000002
            @fact ad.arg[1] --> 0x0000000000000003
            @fact ad.arg[2] --> 0x0000000000000004
            @fact ad.arg[3] --> 0x0000000000000001
            @fact ad.arg[4] --> 0x0000000000000002
            @fact ad.completion_signal.handle --> 0x0000000000000100
        end

        context("can be stored") do
            pkt_bytes = Array(UInt8, 64)
            pkt_ptr = convert(Ptr{Void}, pointer(pkt_bytes))

            ad = HSA.unsafe_convert(HSA.AQLPacket, agent_ptr)
            HSA.unsafe_store!(pkt_ptr, ad)

            @fact pkt_bytes --> agent_bytes
        end
    end

    context("BarrierPacket") do
        context("can be loaded") do
            bp = HSA.unsafe_convert(HSA.AQLPacket, barrier_ptr)

            @fact isa(bp, HSA.BarrierPacket) --> true

            @fact bp.header.typ --> HSA.PacketTypeBarrierAnd

            @fact bp.dep_signal[1] --> 0x0000000000000001
            @fact bp.dep_signal[2] --> 0x0000000000000002
            @fact bp.dep_signal[3] --> 0x0000000000000003
            @fact bp.dep_signal[4] --> 0x0000000000000004
            @fact bp.dep_signal[5] --> 0x0000000000000005
            @fact bp.completion_signal.handle --> 0x0000000000000100
        end

        context("can be stored") do
            pkt_bytes = Array(UInt8, 64)
            pkt_ptr = convert(Ptr{Void}, pointer(pkt_bytes))

            b = HSA.unsafe_convert(HSA.AQLPacket, barrier_ptr)
            HSA.unsafe_store!(pkt_ptr, b)

            @fact pkt_bytes --> barrier_bytes
        end
    end

    context("Can be compared for equality") do
        p1 = HSA.unsafe_convert(HSA.AQLPacket, dispatch_ptr)
        p2 = HSA.unsafe_convert(HSA.AQLPacket, dispatch_ptr)

        @fact p1 == p2 --> true
    end

end


text:
using LinearAlgebra
using Plots

function get_mass2_mat(nu1_mass::Float16, nu2_mass::Float16)::Array
    """Create a 2D mass-squared diagonal matrix"""
    ret_array::Array{Float16} = Diagonal([nu1_mass^2, nu2_mass^2])
    return ret_array
end

function get_mixing_matrix(mixing_angle::Float16, isdeg::Bool = true)::Array
    """Create a 2D mixing/rotation matrix with a given angle"""
    if isdeg
        mixing_angle = deg2rad(mixing_angle)
    end
    mixing_mat::Array{Float16} = [
        +cos(mixing_angle) -sin(mixing_angle)
        +sin(mixing_angle) +cos(mixing_angle)
    ]
    return mixing_mat
end

function time_propagator_mat(
        energy_nu::Float16,
        time_step::Float16,
        delta_mat::Array{Float16}
)::Array{ComplexF16}
    """Calculate the time propagator. Due to a bug in Julia, complex
    matrix arguments for the exponential function are not supported.
    So instead, I use the Euler relation instead.
    """
    exp_arg::Array{ComplexF16} = time_step / (2.0 * energy_nu) * delta_mat
    exp_mat::Array{ComplexF16} = cos(exp_arg) + (0 + 1im) * sin(exp_arg)
    return exp_mat
end

function schrodingerODE()
    """Run a 2D toy example of neutrino oscillations"""

    # Define scalars
    nu1_mass::Float16 = 1.0
    nu2_mass::Float16 = 2.0
    mixing_angle_deg::Float16 = 40
    time_steps::Int = 100
    time_step_size::Float16 = 0.01
    energy_nu::Float16 = 10.0

    # Define physics matrices
    mass2_mat::Array{Float16} = get_mass2_mat(nu1_mass, nu2_mass)
    mixing_mat::Array{Float16} = get_mixing_matrix(mixing_angle_deg, true)
    delta_mat::Array{Float16} = transpose(conj(mixing_mat)) * mass2_mat * mixing_mat

    # Store results in arrays
    # initial state is nu1 = 100%, nu2 = 0%
    initial_state_vector::Vector{ComplexF16} = Vector([1 + 0im, 0 + 0im])
    state_vectors::Array{ComplexF16} = zeros(ComplexF16, (2, time_steps))
    state_fractions::Array{Float16} = zeros(Float16, (2, time_steps))
    state_vectors[:, 1] = initial_state_vector
    state_fractions[:, 1] = Array([1.0; 0.0])

    # Evaluate the wave-function over time and save the state over each step
    for step_index::Int = 2:time_steps
        previous_step_index::Int = step_index - 1
        previous_state_vector::Array{ComplexF16} = state_vectors[:, previous_step_index]
        time_step::Float16 = (step_index - 1) * time_step_size
        time_propagator::Array{ComplexF16} =
            time_propagator_mat(energy_nu, time_step, delta_mat)
        new_state::Vector{ComplexF16} = time_propagator * previous_state_vector
        new_state_fraction::Vector{Float16} =
            Array([abs2(new_state[1]); abs2(new_state[2])])
        state_vectors[:, step_index] = new_state
        state_fractions[:, step_index] = new_state_fraction
    end

    # Plot the results
    Plots.pyplot()
    fig = plot(1:time_steps, state_fractions[1, :], label="\$\\nu_1\$")
    plot!(1:time_steps, state_fractions[2, :], label="\$\\nu_2\$")
    xlabel!("Time step")
    ylabel!("Wave Function Component")
    title!("Flavor Oscillations")
    @show fig
end

schrodingerODE()


text:
module PTimersTests

include("../test_p_timers.jl")

nparts = 5
prun(test_p_timers,sequential,nparts)

end # module


text:
using BVP1D
using PyPlot

const L = 2.0
const γ0 = -1.0
const γL = 2.5
const m = 5.0

const A = γ0 + m
const B = ( γL - A + m*exp(-L)) / L

u(x) = A + B*x - m*exp(-x)
f(x) = m * exp(-x)

P = 6
U, x = solve_bvp(L, f, γ0, γL, P)

xx = range(0; stop=L, length=201)
plot(xx, u.(xx), "-", x[0:P], U[0:P], "o")
xticks(x, [ latexstring("\$x_$p\$") for p = 0:6 ])
legend(("Exact solution", "Finite difference approx"))
grid(true)
#xlabel(L"$x$")
savefig("bvp1d_example.pdf")
~

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
using Test
using Dash
using Dash:DevTools
@testset "DevTools creation" begin
    @test true
    tools = DevTools()
    @test tools.ui == false
    @test tools.props_check == false
    @test tools.serve_dev_bundles == false
    @test tools.hot_reload == false
    @test tools.silence_routes_logging == false
    @test tools.prune_errors == false

    @test tools.hot_reload_interval ≈ 3.
    @test tools.hot_reload_watch_interval ≈ 0.5
    @test tools.hot_reload_max_retry ≈ 8

    tools = DevTools(false)
    @test tools.ui == false
    @test tools.props_check == false
    @test tools.serve_dev_bundles == false
    @test tools.hot_reload == false
    @test tools.silence_routes_logging == false
    @test tools.prune_errors == false

    @test tools.hot_reload_interval ≈ 3.
    @test tools.hot_reload_watch_interval ≈ 0.5
    @test tools.hot_reload_max_retry ≈ 8

    tools = DevTools(true)
    @test tools.ui == true
    @test tools.props_check == true
    @test tools.serve_dev_bundles == true
    @test tools.hot_reload == true
    @test tools.silence_routes_logging == true
    @test tools.prune_errors == true

    @test tools.hot_reload_interval ≈ 3.
    @test tools.hot_reload_watch_interval ≈ 0.5
    @test tools.hot_reload_max_retry ≈ 8

    tools = DevTools(true,
        props_check = false,
        serve_dev_bundles = true,
        hot_reload = false,
        silence_routes_logging = true,
        prune_errors = false,
        hot_reload_interval = 2.2,
        hot_reload_watch_interval = 0.1,
        hot_reload_max_retry = 2
    )
    @test tools.ui == true
    @test tools.props_check == false
    @test tools.serve_dev_bundles == true
    @test tools.hot_reload == false
    @test tools.silence_routes_logging == true
    @test tools.prune_errors == false

    @test tools.hot_reload_interval ≈ 2.2
    @test tools.hot_reload_watch_interval ≈ 0.1
    @test tools.hot_reload_max_retry ≈ 2
end

@testset "DevTools end" begin
        UI = false
        ENV["DASH_UI"] = "False"
        ENV["DASH_PROPS_CHECK"] = "False"
        ENV["DASH_SERVE_DEV_BUNDLES"] = "False"
        ENV["DASH_HOT_RELOAD"] = "False"
        ENV["DASH_SILENCE_ROUTES_LOGGING"] = "false"
        ENV["DASH_PRUNE_ERRORS"] = "0"
        ENV["DASH_HOT_RELOAD_INTERVAL"] = "2.8"
        ENV["DASH_HOT_RELOAD_WATCH_INTERVAL"] = "0.34"
        ENV["DASH_HOT_RELOAD_MAX_RETRY"] = "7"
        tools = DevTools(true)

        @test tools.ui == false
        @test tools.props_check == false
        @test tools.serve_dev_bundles == false
        @test tools.hot_reload == false
        @test tools.silence_routes_logging == false
        @test tools.prune_errors == false

        @test tools.hot_reload_interval ≈ 2.8
        @test tools.hot_reload_watch_interval ≈ 0.34
        @test tools.hot_reload_max_retry ≈ 7

        tools = DevTools(true, hot_reload_interval = 1.6)
        @test tools.hot_reload_interval ≈ 1.6

        delete!(ENV, "DASH_UI")
        delete!(ENV, "DASH_PROPS_CHECK")
        delete!(ENV, "DASH_SERVE_DEV_BUNDLES")
        delete!(ENV, "DASH_HOT_RELOAD")
        delete!(ENV, "DASH_SILENCE_ROUTES_LOGGING")
        delete!(ENV, "DASH_PRUNE_ERRORS")
        delete!(ENV, "DASH_HOT_RELOAD_INTERVAL")
        delete!(ENV, "DASH_HOT_RELOAD_WATCH_INTERVAL")
        delete!(ENV, "DASH_HOT_RELOAD_MAX_RETRY")
end


text:
#using Queryverse
using XLSX
using DataFrames
using Dates
using DelimitedFiles
using CSV

# Parameters

const N_Gens =  132 # number of conventional generators
const N_StorgUs = 2 # number of storage units
const N_Zones = 2
const M_Zones = 2
const N_Blocks =7
const INITIAL_DAY = 1
const FINAL_DAY = 3
const INITIAL_HR_FUCR = 6 # represents the running time for the first WA unit commitment run. INITIAL_HR_FUCR=0 means the FUCR's optimal outcomes are ready at 00:00
const INITIAL_HR_SUCR = 17 # represents the running time for the second WA unit commitment run. INITIAL_HR_SUCR=17 means the SUCR's optimal outcomes are ready at 17:00
const N_Hrs_FUCR = 162 # N_Hrs_FUCR = 168-INITIAL_HR_FUCR, and INITIAL_HR_FUCR runs from 0 to 23; INITIAL_HR_FUCR=6 ==> N_Hrs_FUCR =162
const N_Hrs_SUCR = 151  # N_Hrs_SUCR = 168-INITIAL_HR_SUCR, and INITIAL_HR_SUCR runs from 17 to 23; INITIAL_HR_SUCR=17 ==> N_Hrs_FUCR =168

const FILE_ACCESS_OVER = "w+"
const FILE_ACCESS_APPEND = "a+"

##
########################### Importing input data from the input spreadsheets
# Generators' specification
# Logging file
io_log = open(
    string(
        ".//outputs//log//ReadInputs_",
        Dates.format(now(), "yyyy-mm-dd_HH-MM-SS"),
        ".txt",
    ),
    FILE_ACCESS_APPEND,
)

t1 = time_ns()

DF_Generators = CSV.read(".//inputs//csv//data_generators.csv", DataFrame);
Map_Gens = readdlm(".//inputs//csv//location_generators.csv", ','; header = true);
DF_Storage = CSV.read(".//inputs//csv//data_storage.csv", DataFrame);
Map_Storage = readdlm(".//inputs//csv//location_storage.csv", ','; header = true);

FUCR_Demands = readdlm(".//inputs//csv//data_demand.csv", ','; header = true);
SUCR_Demands = readdlm(".//inputs//csv//data_demand_updated.csv", ','; header = true);
BUCR_Demands = readdlm(".//inputs//csv//data_demand_actual.csv", ','; header = true);

# solar generation data at each location
FUCR_SolarGs = readdlm(".//inputs//csv//data_solar.csv", ','; header = true);
SUCR_SolarGs = readdlm(".//inputs//csv//data_solar_updated.csv", ','; header = true);
BUCR_SolarGs = readdlm(".//inputs//csv//data_solar_actual.csv", ','; header = true);

# wind energy data for each location
FUCR_WindGs = readdlm(".//inputs//csv//data_wind.csv", ','; header = true);
SUCR_WindGs = readdlm(".//inputs//csv//data_wind_updated.csv", ','; header = true);
BUCR_WindGs = readdlm(".//inputs//csv//data_wind_actual.csv", ','; header = true);

#hydro generation data for each location
FUCR_HydroGs = readdlm(".//inputs//csv//data_hydro.csv", ','; header = true);
SUCR_HydroGs = readdlm(".//inputs//csv//data_hydro_updated.csv", ','; header = true);
BUCR_HydroGs = readdlm(".//inputs//csv//data_hydro_actual.csv", ','; header = true);

tranC = readdlm(".//inputs//csv//LineCapacity.csv", ','; header = true);
tranS = readdlm(".//inputs//csv//LineSusceptance.csv", ','; header = true);
reserve_reqs = readdlm(".//inputs//csv//data_reserve_reqs.csv", ','; header = true);

fuelPrice = readdlm(".//inputs//csv//data_fuel_price.csv", ','; header = true);

fuelPrice_head = fuelPrice[2];
fuelPrice = fuelPrice[1];
fuelPrice = fuelPrice[2:N_Gens+1, 4:368];
#FuelPrice = XLSX.readdata(".\\inputs\\data_generators.XLSX", "data_fuel_price", "D3:ND134")

# Reorganize data
FUCR_Demands_head = FUCR_Demands[2];
FUCR_Demands = FUCR_Demands[1];
FUCR_Demands = FUCR_Demands[1:53999, 4:(4+N_Zones-1)];

SUCR_Demands_head = SUCR_Demands[2];
SUCR_Demands = SUCR_Demands[1];
SUCR_Demands = SUCR_Demands[1:53999, 4:(4+N_Zones-1)];

BUCR_Demands_head = BUCR_Demands[2];
BUCR_Demands = BUCR_Demands[1];
BUCR_Demands = BUCR_Demands[1:7728, 3:(3+N_Zones-1)];

# Solar
FUCR_SolarGs_head = FUCR_SolarGs[2];
FUCR_SolarGs = FUCR_SolarGs[1];
FUCR_SolarGs = FUCR_SolarGs[1:7728, 3:(3+N_Zones-1)];

SUCR_SolarGs_head = SUCR_SolarGs[2];
SUCR_SolarGs = SUCR_SolarGs[1];
SUCR_SolarGs = SUCR_SolarGs[1:7728, 3:(3+N_Zones-1)];

BUCR_SolarGs_head = BUCR_SolarGs[2];
BUCR_SolarGs = BUCR_SolarGs[1];
BUCR_SolarGs = BUCR_SolarGs[1:7728, 3:(3+N_Zones-1)];

#Wind
FUCR_WindGs_head = FUCR_WindGs[2];
FUCR_WindGs = FUCR_WindGs[1];
FUCR_WindGs = FUCR_WindGs[1:7728, 3:(3+N_Zones-1)];

SUCR_WindGs_head = SUCR_WindGs[2];
SUCR_WindGs = SUCR_WindGs[1];
SUCR_WindGs = SUCR_WindGs[1:7728, 3:(3+N_Zones-1)];

BUCR_WindGs_head = BUCR_WindGs[2];
BUCR_WindGs = BUCR_WindGs[1];
BUCR_WindGs = BUCR_WindGs[1:7728, 3:(3+N_Zones-1)];

#Hydro
FUCR_HydroGs_head = FUCR_HydroGs[2];
FUCR_HydroGs = FUCR_HydroGs[1];
FUCR_HydroGs = FUCR_HydroGs[1:7728, 2:(2+N_Zones-1)];

SUCR_HydroGs_head = SUCR_HydroGs[2];
SUCR_HydroGs = SUCR_HydroGs[1];
SUCR_HydroGs = SUCR_HydroGs[1:7728, 2:(2+N_Zones-1)];

BUCR_HydroGs_head = BUCR_HydroGs[2];
BUCR_HydroGs = BUCR_HydroGs[1];
BUCR_HydroGs = BUCR_HydroGs[1:7728, 2:(2+N_Zones-1)];

Map_Gens_head = Map_Gens[2];
Map_Gens = Map_Gens[1];
Map_Gens = Map_Gens[:,3:N_Zones+2]

Map_Storage_head = map_storage[2];
Map_Storage = Map_Storage[1];
Map_Storage = Map_Storage[:,2:N_Zones+1]
#map_storage = convert(Array{Int32,2}, map_storage[:,2:N_Zones+1]);

tranC_head = tranC[2];
tranC = tranC[1];
tranC = tranC[1:N_Zones,2:(2+N_Zones-1)];

tranS_head = tranS[2];
tranS = tranS[1];
tranS = tranS[1:N_Zones,2:(2+N_Zones-1)];

reserve_reqs_head = reserve_reqs[2];
reserve_reqs = reserve_reqs[1];
reserve_req_up = reserve_reqs
#reserve_req_dn = reserve_reqs[:,3];

##
# Importing input data from the input spreadsheets
# Generators' specification
#DF_Generators = DataFrame(XLSX.readtable(".\\inputs\\data_generators.XLSX", "data_generators")...)

# Generators location map: if a generator g is located in zone z Map_Gens[g,z]=1; and 0 otherwise
#DF_Map_Gens = DataFrame(XLSX.readtable(".\\inputs\\data_generators.XLSX", "location_generators")...)
#Map_Gens = convert(Matrix, DF_Map_Gens[:,3:N_Zones+2])

# Storage Units' specification and location
#DF_Storage = DataFrame(XLSX.readtable(".\\inputs\\data_storage.XLSX", "data_storage")...) # storage specs
#DF_Map_Storage = DataFrame(XLSX.readtable(".\\inputs\\data_storage.XLSX", "location_storage")...) # storage location as a dataframe
#Map_Storage = convert(Matrix, DF_Map_Storage[:,2:N_Zones+1]) # convert storage location data to  a matrix

# energy demand at each location
#DF_Dems = DataFrame(XLSX.readtable(".\\Inputs\\data_demand.XLSX", "data_demand")...)
#FUCR_Demands = XLSX.readdata(".\\inputs\\data_demand.XLSX", "data_demand", "D2:E54000")
#SUCR_Demands = XLSX.readdata(".\\inputs\\data_demand.XLSX", "data_demand_updated", "D2:E54000")
#BUCR_Demands = XLSX.readdata(".\\inputs\\data_demand.XLSX", "data_demand_actual", "C2:D7729")
# There is no map for the demand data. Instead we take the input demand data for each zone. In other words, Demand[t,z] represents demand at zone z and time t

# solar generation data at each location
#FUCR_SolarGs = XLSX.readdata(".\\inputs\\data_solar.XLSX", "data_solar", "C2:D7729")
#SUCR_SolarGs = XLSX.readdata(".\\inputs\\data_solar.XLSX", "data_solar_updated", "C2:D7729")
#BUCR_SolarGs = XLSX.readdata(".\\inputs\\data_solar.XLSX", "data_solar_actual", "C2:D7729")

# wind energy data for each location
#FUCR_WindGs = XLSX.readdata(".\\inputs\\data_wind.XLSX", "data_wind", "C2:D7729")
#SUCR_WindGs = XLSX.readdata(".\\inputs\\data_wind.XLSX", "data_wind_updated", "C2:D7729")
#BUCR_WindGs = XLSX.readdata(".\\inputs\\data_wind.XLSX", "data_wind_actual", "C2:D7729")

#hydro generation data for each location
#FUCR_HydroGs = XLSX.readdata(".\\inputs\\data_hydro.XLSX", "data_hydro", "B2:C7729")
#SUCR_HydroGs = XLSX.readdata(".\\inputs\\data_hydro.XLSX", "data_hydro_updated", "B2:C7729")
#BUCR_HydroGs = XLSX.readdata(".\\inputs\\data_hydro.XLSX", "data_hydro_actual", "B2:C7729")

#Zonal Reserve Targets
#Reserve_Req_Up = XLSX.readdata(".\\inputs\\data_reserve_reqs.XLSX", "data_reserve_reqs", "A2:B2")

# Transmission system data (Capacity and Susceptance)
#TranC = XLSX.readdata(".\\inputs\\data_tranSmission.XLSX", "LineCapacity","B2:C3") # Transmission line capacity
#TranS = XLSX.readdata(".\\inputs\\data_tranSmission.XLSX", "LineSusceptance","B2:C3")# Transmission line susceptance

# Daily Fuel Price data
#FuelPrice = XLSX.readdata(".\\inputs\\data_generators.XLSX", "data_fuel_price", "D3:ND134")

##
t2 = time_ns()
elapsedTime = (t2 -t1)/1.0e9;

write(io_log, "Whole program time execution (s):\t $elapsedTime\n")
close(io_log);

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
module ShallowWater

export SWModel, SWProblem

using StaticArrays
using ..VariableTemplates
using LinearAlgebra: I, dot
using ..PlanetParameters: grav

import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                        vars_diffusive, vars_integrals, flux_nondiffusive!,
                        flux_diffusive!, source!, wavespeed,
                        boundary_state!,
                        gradvariables!, init_aux!, init_state!,
                        LocalGeometry, diffusive!

using ..DGmethods.NumericalFluxes

×(a::SVector, b::SVector) = StaticArrays.cross(a, b)
⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b)
⊗(a::SVector, b::SVector) = a * b'

abstract type SWProblem end

abstract type TurbulenceClosure end
struct LinearDrag{L} <: TurbulenceClosure
  λ::L
end
struct ConstantViscosity{L} <: TurbulenceClosure
  ν::L
end

abstract type AdvectionTerm end
struct NonLinearAdvection <: AdvectionTerm end

struct SWModel{P, T, A, S} <: BalanceLaw
  problem::P
  turbulence::T
  advection::A
  c::S
end

function vars_state(m::SWModel, T)
  @vars begin
    U::SVector{3, T}
    η::T
  end
end

function vars_aux(m::SWModel, T)
  @vars begin
    f::SVector{3, T}
    τ::SVector{3, T}  # value includes τₒ, g, and ρ
  end
end

function vars_gradient(m::SWModel, T)
  @vars begin
    U::SVector{3, T}
  end
end

function vars_diffusive(m::SWModel, T)
  @vars begin
    ν∇U::SMatrix{3, 3, T, 9}
  end
end

@inline function flux_nondiffusive!(m::SWModel, F::Grad, q::Vars,
                                    α::Vars, t::Real)
  U = q.U
  η = q.η
  H = m.problem.H

  F.η += U
  F.U += grav * H * η * I

  advective_flux!(m, m.advection, F, q, α, t)

  return nothing
end

advective_flux!(::SWModel, ::Nothing, _...) = nothing

@inline function advective_flux!(m::SWModel, A::NonLinearAdvection, F::Grad,
                                 q::Vars, α::Vars, t::Real)
  U = q.U
  H = m.problem.H

  F.U += 1 / H * U ⊗ U

  return nothing
end

function gradvariables!(m::SWModel, f::Vars, q::Vars, α::Vars, t::Real)
  gradvariables!(m.turbulence, f, q, α, t)
end

gradvariables!(::LinearDrag, _...) = nothing

@inline function gradvariables!(T::ConstantViscosity, f::Vars, q::Vars,
                                α::Vars, t::Real)
  f.U = q.U

  return nothing
end

function diffusive!(m::SWModel, σ::Vars, δ::Grad, q::Vars, α::Vars, t::Real)
  diffusive!(m.turbulence, σ, δ, q, α, t)
end

diffusive!(::LinearDrag, _...) = nothing

@inline function diffusive!(T::ConstantViscosity, σ::Vars, δ::Grad, q::Vars,
                            α::Vars, t::Real)
  ν  = T.ν
  ∇U = δ.U

  σ.ν∇U = ν * ∇U

  return nothing
end

function flux_diffusive!(m::SWModel, G::Grad, q::Vars, σ::Vars, ::Vars,
                         α::Vars, t::Real)
  flux_diffusive!(m.turbulence, G, q, σ, α, t)
end

flux_diffusive!(::LinearDrag, _...) = nothing

@inline function flux_diffusive!(::ConstantViscosity, G::Grad, q::Vars,
                                 σ::Vars, α::Vars, t::Real)
  G.U -= σ.ν∇U

  return nothing
end

@inline wavespeed(m::SWModel, n⁻, q::Vars, α::Vars, t::Real) = m.c

@inline function source!(m::SWModel{P}, S::Vars, q::Vars,
                         diffusive::Vars, α::Vars, t::Real) where P
  τ = α.τ
  f = α.f
  U = q.U
  S.U += τ - f × U

  linear_drag!(m.turbulence, S, q, α, t)

  return nothing
end

linear_drag!(::ConstantViscosity, _...) = nothing

@inline function linear_drag!(T::LinearDrag, S::Vars, q::Vars,
                              α::Vars, t::Real)
  λ = T.λ
  U = q.U

  S.U -= λ * U

  return nothing
end

function shallow_init_aux! end
function init_aux!(m::SWModel, aux::Vars, geom::LocalGeometry)
  shallow_init_aux!(m.problem, aux, geom)
end

function shallow_init_state! end
function init_state!(m::SWModel, state::Vars, aux::Vars, coords, t)
  shallow_init_state!(m.problem, m.turbulence, state, aux, coords, t)
end

function boundary_state!(nf, m::SWModel, state⁺::Vars, aux⁺::Vars, n⁻,
                         state⁻::Vars, aux⁻::Vars, bctype, t, _...)
    shallow_boundary_state!(nf, m, m.turbulence, state⁺, aux⁺, n⁻, state⁻,
                          aux⁻, t)
end

@inline function shallow_boundary_state!(::Rusanov, m::SWModel, ::LinearDrag, state⁺,
                               aux⁺, n⁻, state⁻, aux⁻, t)
  U⁻ = state⁻.U
  n⁻ = SVector(n⁻)

  state⁺.η = state⁻.η
  state⁺.U = U⁻ - 2 * (n⁻⋅U⁻) * n⁻

  return nothing
end

shallow_boundary_state!(::CentralNumericalFluxGradient, m::SWModel,
                        ::LinearDrag, _...) = nothing

shallow_boundary_state!(::CentralNumericalFluxDiffusive, m::SWModel,
                        ::LinearDrag, _...) = nothing

function boundary_state!(nf, m::SWModel, q⁺::Vars, σ⁺::Vars, α⁺::Vars,
                         n⁻, q⁻::Vars, σ⁻::Vars, α⁻::Vars, bctype, t, _...)
  shallow_boundary_state!(nf, m, m.turbulence, q⁺, σ⁺, α⁺, n⁻, q⁻, σ⁻, α⁻, t)
end

@inline function shallow_boundary_state!(::Rusanov, m::SWModel,
                                         ::ConstantViscosity,
                                         q⁺, α⁺, n⁻, q⁻, α⁻, t)
  q⁺.η =  q⁻.η
  q⁺.U = -q⁻.U

  return nothing
end

@inline function shallow_boundary_state!(::CentralNumericalFluxGradient, m::SWModel,
                                         ::ConstantViscosity, q⁺, α⁺, n⁻, q⁻, α⁻, t)
  q⁺.U = 0

  return nothing
end

@inline function shallow_boundary_state!(::CentralNumericalFluxDiffusive,
                                         m::SWModel,
                                         ::ConstantViscosity, q⁺, σ⁺, α⁺,
                                         n⁻, q⁻, σ⁻, α⁻, t)
  q⁺.U   = -q⁻.U
  σ⁺.ν∇U =  σ⁻.ν∇U

  return nothing
end

end

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
module AutoDiff

using SparseArrays

using CUDA
import CUDA.CUSPARSE
import ForwardDiff
import SparseDiffTools
using KernelAbstractions

using ..ExaPF: Spmat, xzeros, State, Control

import Base: show

"""
    AbstractJacobian

Automatic differentiation for the compressed Jacobian of
any nonlinear constraint ``h(x)``.

"""
abstract type AbstractJacobian end

"""
    AbstractHessian

Automatic differentiation for the adjoint-Hessian-vector product ``λ^⊤ H v`` of
any nonlinear constraint ``h(x)``.

"""
abstract type AbstractHessian end

abstract type AbstractAdjointStack{VT} end

"""
    jacobian!(form::AbstractFormulation, jac::AutoDiff.AbstractJacobian, x)

Update inplace the Jacobian ``J`` stored inside `jac` at a new point `x`.
"""
function jacobian! end

"""
    adj_hessian_prod!(form::AbstractFormulation, H::AutoDiff.AbstractHessian, hv, x, λ, v)

Compute the adjoint-Hessian-vector product ``λ^⊤ H v`` at a given
point `x`, and store the result inplace in vector `hv`.

"""
function adj_hessian_prod! end

"""
    AutoDiff.Jacobian <: AbstractJacobian

Creates an object to compute the Jacobian with ForwardDiff.

### Attributes

* `func::Func`: base function to differentiate
* `var::Union{State,Control}`: specify whether we are differentiating w.r.t. the state or the control.
* `J::SMT`: Sparse uncompressed Jacobian to be used by linear solver. This is either of type `SparseMatrixCSC` or `CuSparseMatrixCSR`.
* `compressedJ::MT`: Dense compressed Jacobian used for updating values through AD either of type `Matrix` or `CuMatrix`.
* `coloring::VI`: Row coloring of the Jacobian.
* `t1sseeds::VP`: The seeding vector for AD built based on the coloring.
* `t1sF::VD`: Output array of active (AD) type.
* `x::VT`: Input array of passive type. This includes both state and control.
* `t1sx::VD`: Input array of active type.
* `map::VI`: State and control mapping to array `x`
* `varx::SubT`: View of `map` on `x`
* `t1svarx::SubD`: Active (AD) view of `map` on `x`

"""
struct Jacobian{Func, VI, VT, MT, SMT, VP, VD, SubT, SubD} <: AbstractJacobian
    func::Func
    var::Union{State,Control}
    J::SMT
    compressedJ::MT
    coloring::VI
    t1sseeds::VP
    t1sF::VD
    x::VT
    t1sx::VD
    map::VI
    # Cache views on x and its dual vector to avoid reallocating on the GPU
    varx::SubT
    t1svarx::SubD
end

"""
    AutoDiff.ConstantJacobian <: AbstractJacobian

Creates a constant Jacobian object for a linear function ``h(x)``.
Using a `ConstantJacobian` object allows to avoid computing
the full Jacobian with AutoDiff when it is not necessary.
"""
struct ConstantJacobian{SMT} <: AbstractJacobian
    J::SMT
end

"""
    AutoDiff.Hessian

Creates an object for computing Hessian adjoint tangent projections.

* `func::Func`: base function to differentiate.
* `host_t1sseeds::VHP`: Seeding vector for seeding on the host.
* `t1sseeds::VP`: The seeding vector for AD built based on the coloring.
* `x::VT`: Input array of passive type. This includes both state and control.
* `t1sF::VD`: Output array of active (AD) type.
* `∂t1sF::VD`: Adjoint of the output array.
* `t1sx::VD`: Input array of active type.
* `∂t1sx::VD`: Adjoint of the input array.
* `map::VI`: State and control mapping to array `x`
* `varx::SubT`: View of `map` on `x`
* `t1svarx::SubD`: Active (AD) view of `map` on `x`
* `buffer::Buff`: cache for computing the adjoint (could be `Nothing`)
"""
struct Hessian{Func, VI, VT, VHP, VP, VD, SubT, SubD, Buff} <: AbstractHessian
    func::Func
    host_t1sseeds::VHP # Needed because seeds have to be created on the host
    t1sseeds::VP
    x::VT
    t1sF::VD
    ∂t1sF::VD
    t1sx::VD
    ∂t1sx::VD
    map::VI
    varx::SubT
    t1svarx::SubD
    buffer::Buff
end

# Cache for adjoint
"""
    TapeMemory{F, S, I}

This object is used as a buffer to compute the adjoint of a given function
``h(x)``. It stores internally all intermediate values necessary
to compute the adjoint, and cache the stack used in the backward pass.

## Note
This structure is largely inspired from [ChainRulesCore.jl](https://juliadiff.org/ChainRulesCore.jl/stable/design/changing_the_primal.html#The-Journey-to-rrule).
"""
struct TapeMemory{F, S, I}
    func::F
    stack::S
    intermediate::I
end

# Seeding
# TODO
@kernel function _init_seed!(t1sseeds, t1sseedvecs, coloring, ncolor)
    i = @index(Global, Linear)
    t1sseedvecs[:,i] .= 0
    @inbounds for j in 1:ncolor
        if coloring[i] == j
            t1sseedvecs[j,i] = 1.0
        end
    end
    t1sseeds[i] = ForwardDiff.Partials{ncolor, Float64}(NTuple{ncolor, Float64}(t1sseedvecs[:,i]))
end

function init_seed(coloring, ncolor, nmap)
    t1sseeds = Vector{ForwardDiff.Partials{ncolor, Float64}}(undef, nmap)
    t1sseedvecs = zeros(Float64, ncolor, nmap)
    # The seeding is always done on the CPU since it's faster
    ev = _init_seed!(CPU())(t1sseeds, t1sseedvecs, Array(coloring), ncolor, ndrange=nmap, dependencies=Event(CPU()))
    wait(ev)
    return t1sseeds
end

@kernel function seed_kernel!(
    duals::AbstractArray{ForwardDiff.Dual{T, V, N}}, x,
    seeds::AbstractArray{ForwardDiff.Partials{N, V}}
) where {T,V,N}
    i = @index(Global, Linear)
    duals[i] = ForwardDiff.Dual{T,V,N}(x[i], seeds[i])
end

"""
    seed!

Calling the seeding kernel.
Seeding is parallelized over the `ncolor` number of duals.

"""
function seed!(t1sseeds, varx, t1svarx)
    if isa(t1sseeds, Vector)
        device = CPU()
        kernel! = seed_kernel!(CPU())
    else
        device = CUDADevice()
        kernel! = seed_kernel!(CUDADevice())
    end
    ev = kernel!(t1svarx, varx, t1sseeds, ndrange=length(t1svarx), dependencies=Event(device))
    wait(ev)
end

# Get partials
@kernel function getpartials_kernel_cpu!(compressedJ, t1sF)
    i = @index(Global, Linear)
    compressedJ[:, i] .= ForwardDiff.partials.(t1sF[i]).values
end

@kernel function getpartials_kernel_gpu!(compressedJ, t1sF)
    i = @index(Global, Linear)
    for j in eachindex(ForwardDiff.partials.(t1sF[i]).values)
        @inbounds compressedJ[j, i] = ForwardDiff.partials.(t1sF[i]).values[j]
    end
end

# Get partials for Hessian projection
@kernel function getpartials_hv_kernel!(hv, adj_t1sx, map)
    i = @index(Global, Linear)
    hv[i] = ForwardDiff.partials(adj_t1sx[map[i]]).values[1]
end

"""
    getpartials_kernel!(compressedJ, t1sF)

Calling the partial extraction kernel.
Extract the partials from the AutoDiff dual type on the target
device and put it in the compressed Jacobian `compressedJ`.

"""
function getpartials_kernel!(hv::AbstractVector, adj_t1sx, map)
    if isa(hv, Array)
        device = CPU()
        kernel! = getpartials_hv_kernel!(CPU())
    else
        device = CUDADevice()
        kernel! = getpartials_hv_kernel!(CUDADevice())
    end
    ev = kernel!(hv, adj_t1sx, map, ndrange=length(hv), dependencies=Event(device))
    wait(ev)
end

function getpartials_kernel!(compressedJ::AbstractMatrix, t1sF)
    if isa(compressedJ, Array)
        device = CPU()
        kernel! = getpartials_kernel_cpu!(CPU())
    else
        device = CUDADevice()
        kernel! = getpartials_kernel_gpu!(CUDADevice())
    end
    ev = kernel!(compressedJ, t1sF, ndrange=length(t1sF), dependencies=Event(device))
    wait(ev)
end

# Uncompress kernels
@kernel function uncompress_kernel_gpu!(J_rowPtr, J_colVal, J_nzVal, compressedJ, coloring)
    i = @index(Global, Linear)
    @inbounds for j in J_rowPtr[i]:J_rowPtr[i+1]-1
        @inbounds J_nzVal[j] = compressedJ[coloring[J_colVal[j]], i]
    end
end

@kernel function uncompress_kernel_cpu!(J_colptr, J_rowval, J_nzval, compressedJ, coloring)
    # CSC is column oriented: nmap is equal to number of columns
    i = @index(Global, Linear)
    @inbounds for j in J_colptr[i]:J_colptr[i+1]-1
        @inbounds J_nzval[j] = compressedJ[coloring[i], J_rowval[j]]
    end
end

"""
    uncompress_kernel!(J, compressedJ, coloring)

Uncompress the compressed Jacobian matrix from `compressedJ`
to sparse CSC (on the CPU) or CSR (on the GPU).
"""
function uncompress_kernel!(J, compressedJ, coloring)
    if isa(J, SparseArrays.SparseMatrixCSC)
        kernel! = uncompress_kernel_cpu!(CPU())
        ev = kernel!(J.colptr, J.rowval, J.nzval, compressedJ, coloring, ndrange=size(J,2), dependencies=Event(CPU()))
    else
        kernel! = uncompress_kernel_gpu!(CUDADevice())
        ev = kernel!(J.rowPtr, J.colVal, J.nzVal, compressedJ, coloring, ndrange=size(J,1), dependencies=Event(CUDADevice()))
    end
    wait(ev)
end

function Base.show(io::IO, jacobian::Jacobian)
    println(io, "A AutoDiff Jacobian for $(jacobian.func) (w.r.t. $(jacobian.var))")
    ncolor = size(jacobian.compressedJ, 1)
    print(io, "Number of Jacobian colors: ", ncolor)
end

end


text:
println("---> START <---")
#cd(raw"C:\Users\Nicky\Downloads\COVID_Simulator_Toronto\final_codes")
#pwd()
# Set up environment
using Pkg
Pkg.activate(".")

using JSON
using Glob
using CSV
using Random
using DataFrames
using Conda
using PyCall
using Serialization
using SparseArrays

using Statistics
using StatsBase

using Distributions
using DataStructures
using Parameters

using Geodesy
using SimpleWeightedGraphs
using OpenStreetMapX
using LightGraphs

using Plots
using ColorSchemes
using Colors
using Dates
using GR

println("---> PACKAGES DONE <---")

# Include codes
include("TTC_data_read_and_prep_functions.jl")
include("graph_from_TTC_data_functions.jl")
include("read_and_prep_osm_map_functions.jl")
include("complete_map_with_TTC_functions.jl")
include("create_simulation_objects_functions.jl")
include("simulation_dynamic_functions.jl")
include("init_and_run_sim_functions.jl")
include("gather_statistics_from_sim_functions.jl")
include("post_simulation_plots_functions.jl")
println("---> INCLUDES DONE <---")

# Read TTC data and create needed structures
gg,secs_between_stops,routes_dict,route_types,routes_filtered,routes_types_both_dir = read_and_prep_TTC_data(data_folder = "./data_input/");

# Define simulation parameters
sim_params_comb_df = create_df_with_sim_params(TTC_car_freq = 3*60
                            , max_load_gov_restr = 1.0
                            , when_to_run_people = 60.0*60.0*1.0
                            , when_to_run_wagons = 0.0
                            , sim_run_time_in_secs = 3*24*60*60
                            , N_agents = 1000
                            , agents_speed_in_sec_per_m = 0.9
                            , max_load_subway = 10
                            , max_load_streetcar = 10
                            , p0 = 1/1000);
CSV.write("./data_output/sim_params_comb_df.csv",sim_params_comb_df);

# Run one time
sim_params = Simulation_Parameters(copy(sim_params_comb_df[1,2:end])...);
freq_to_check = deepcopy(sim_params.TTC_car_freq)#600
speed_to_check = deepcopy(sim_params.agents_speed_in_sec_per_m)#0.8
if isfile("./data_output/map_prepared_$freq_to_check.....$speed_to_check.bin") == false
    m, orig_map_nodes_num, routes_path, routes_distances  = read_and_prep_osm_map(map_path = "./data_input/central_toronto.osm"
                                                                                , routes_dict = routes_dict
                                                                                , secs_between_stops = secs_between_stops
                                                                                , parset_id = 1
                                                                                , serialize_results = true

                                                                                , TTC_car_freq = sim_params.TTC_car_freq#600
                                                                                , agents_speed_in_sec_per_m = sim_params.agents_speed_in_sec_per_m)
else
    println("Objects with suffix $freq_to_check.....$speed_to_check.bin already exist and will be used.")
    m = deserialize("./data_output/map_prepared_$freq_to_check.....$speed_to_check.bin")
    orig_map_nodes_num = deserialize("./data_output/orig_map_nodes_num_prepared_$freq_to_check.....$speed_to_check.bin")
    routes_path = deserialize("./data_output/routes_path_prepared_$freq_to_check.....$speed_to_check.bin")
    routes_distances = deserialize("./data_output/routes_distances_prepared_$freq_to_check.....$speed_to_check.bin")
end

@time simul_res = initialize_sim_and_run(m = m
        , routes_path = routes_path
        , routes_distances = routes_distances
        , routes_types_both_dir = routes_types_both_dir
        , orig_map_nodes_num = orig_map_nodes_num
        , serialize_finished_sim_object = false
        , parset_id = 1

        , when_to_run_people = sim_params.when_to_run_people
        , when_to_run_wagons = sim_params.when_to_run_wagons
        , TTC_freq = sim_params.TTC_car_freq
        , max_load_subway= sim_params.max_load_subway
        , max_load_streetcar = sim_params.max_load_streetcar
        , max_load_gov_restr=  sim_params.max_load_gov_restr
        , N_agents =  sim_params.N_agents
        , sim_run_time_in_secs = sim_params.sim_run_time_in_secs
        , p0 =  sim_params.p0);

simul_res

# Select one parameters set
for par_set_id in 1:nrow(sim_params_comb_df)
    println("CURRENT PARSET ID: ", par_set_id)
    sim_params = Simulation_Parameters(copy(sim_params_comb_df[par_set_id,2:end])...);
    println("Simulation parameters selected:\n", sim_params)
    freq_to_check = deepcopy(sim_params.TTC_car_freq)#600
    speed_to_check = deepcopy(sim_params.agents_speed_in_sec_per_m)#0.8
    # Create map and add TTC
    if isfile("./data_output/map_prepared_$freq_to_check.....$speed_to_check.bin") == false
        m, orig_map_nodes_num, routes_path, routes_distances  = read_and_prep_osm_map(map_path = "./data_input/central_toronto.osm"
                                                                                , routes_dict = routes_dict
                                                                                , secs_between_stops = secs_between_stops
                                                                                , parset_id = par_set_id
                                                                                , serialize_results = true

                                                                                , TTC_car_freq = sim_params.TTC_car_freq#600
								                                                        , agents_speed_in_sec_per_m = sim_params.agents_speed_in_sec_per_m)
    else
        println("Objects with suffix $freq_to_check.....$speed_to_check.bin already exist and will be used.")
        m = deserialize("./data_output/map_prepared_$freq_to_check.....$speed_to_check.bin")
        orig_map_nodes_num = deserialize("./data_output/orig_map_nodes_num_prepared_$freq_to_check.....$speed_to_check.bin")
        routes_path = deserialize("./data_output/routes_path_prepared_$freq_to_check.....$speed_to_check.bin")
        routes_distances = deserialize("./data_output/routes_distances_prepared_$freq_to_check.....$speed_to_check.bin")
    end

    # Run simulation
    simul_res = initialize_sim_and_run_N_times_and_gather_results(m = m
            , routes_path = routes_path
            , routes_distances = routes_distances
            , routes_types_both_dir = routes_types_both_dir
            , orig_map_nodes_num = orig_map_nodes_num
            , serialize_finished_sim_object = false
            , parset_id = par_set_id
            , N_times_to_run = 2

            , when_to_run_people = sim_params.when_to_run_people
            , when_to_run_wagons = sim_params.when_to_run_wagons
            , TTC_freq = sim_params.TTC_car_freq
            , max_load_subway= sim_params.max_load_subway
            , max_load_streetcar = sim_params.max_load_streetcar
            , max_load_gov_restr=  sim_params.max_load_gov_restr
            , N_agents =  sim_params.N_agents
            , sim_run_time_in_secs = sim_params.sim_run_time_in_secs
            , p0 =  sim_params.p0);

    # Save results to JSON
    simul_res_to_json = JSON.json(simul_res)
    open("./data_output/simul_res_parset_id_$par_set_id.json", "w") do f
            write(f, simul_res_to_json)
         end
end

# Load json to Julia
my_json=join(readlines("./data_output/simul_res_parset_id_1.json"))
my_json_dict=JSON.parse(my_json)

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
using Base.Test
using StatsBase
import OnlineMoments: Mean, Variance, Moments

const N = 10000
const n = 100

for r = 1:1000
  x = rand(N)
  m = let
    m = zeros(Mean, n)
    for j = 1:n
      for i = j:n:N
        m[j] += x[i]
      end
    end
    sum(m)
  end
  @test nobs(m) == length(x)
  @test_approx_eq_eps mean(m) mean(x) 1e-14
end

for r = 1:1000
  x = rand(N)
  m = let
    m = zeros(Variance, n)
    for j = 1:n
      for i = j:n:N
        m[j] += x[i]
      end
    end
    sum(m)
  end
  @test nobs(m) == length(x)
  @test_approx_eq_eps mean(m) mean(x) 1e-14
  @test_approx_eq_eps var(m) var(x) 1e-14
  @test_approx_eq_eps std(m) std(x) 1e-14
  @test_approx_eq_eps sem(m) sem(x) 1e-14
end

for r = 1:1000
  x = rand(N)
  m = let
    m = zeros(Moments, n)
    for j = 1:n
      for i = j:n:N
        m[j] += x[i]
      end
    end
    sum(m)
  end
  @test nobs(m) == length(x)
  @test_approx_eq_eps mean(m) mean(x) 1e-14
  @test_approx_eq_eps var(m) var(x) 1e-14
  @test_approx_eq_eps std(m) std(x) 1e-14
  @test_approx_eq_eps sem(m) sem(x) 1e-14
  @test_approx_eq_eps skewness(m) skewness(x) 1e-13
  @test_approx_eq_eps kurtosis(m) kurtosis(x) 1e-13
end


text:
function generate_network(rng::AbstractRNG, n_nodes::Int, n_comms::Int, p_intra_community::Float64, p_inter_community::Float64)
    network = CommunityNetwork()
    network.graph = Graph(n_nodes)
    network.n_comms = n_comms
    network.labels = Array(Int, n_nodes)
    for i in 1:n_nodes
        network.labels[i] = ceil(rand(rng)*n_comms)
    end
    # generate with SBM
    for i in 1:n_nodes
        for j in i+1:n_nodes
            if network.labels[i] == network.labels[j]
                if rand(rng) < p_intra_community
                    add_edge!(network.graph, i, j)
                end
            else
                if rand(rng) < p_inter_community
                    add_edge!(network.graph, i, j)
                end
            end
        end
    end
    return network
end

function generate_problem(rng::AbstractRNG, net::CommunityNetwork, budget::Int, node::Int, reward::Float64, n_known_labels::Int, n_known_edges::Int, p_intra::Float64, p_inter::Float64)
    problem = CommunityAffiliationPOMDP(budget,
                                        node,
                                        reward,
                                        net,
                                        Set{Int}(),
                                        Set{Edge}(),
                                        Dict{Int,Int}(),
                                        p_intra,
                                        p_inter)

    problem.true_network.budget = budget

    edge_vector = collect(net.graph.edges)
    available = ones(Bool, length(edge_vector))
    for i in 1:n_known_edges
        r = ceil(rand(rng)*sum(available))
        ed = edge_vector[available][r]
        push!(problem.known_edges, ed)
        push!(problem.known_nodes, ed.first)
        push!(problem.known_nodes, ed.second)
        available[available][r] = false
    end

    node_vector = collect(problem.known_nodes)
    available = ones(Bool, length(node_vector))
    available[find(n->(n==node),node_vector)] = false
    for i in 1:n_known_labels
        r = ceil(rand(rng)*sum(available))
        n = node_vector[available][r]
        problem.known_labels[n] = problem.true_network.labels[n]
        available[available][r] = false
    end

    return problem
end