text:
@testset "Initialise" begin
    model = initialise()
    @test sort(collect(keys(model.properties))) == [
        :critical_nutrients,
        :identifier,
        :init_nutrients,
        :init_pike_mortality,
        :knowledge,
        :lake,
        :max_sewage_water,
        :nutrient_change,
        :nutrient_series,
        :outcomes_upgraded_households_sum,
        :outcomes_year_of_full_upgrade,
        :outcomes_year_when_desired_level_is_back,
        :outcomes_year_when_desired_pike_is_back,
        :outcomes_year_when_informed,
        :outcomes_year_when_nutrients_became_critical,
        :outcomes_year_when_pike_became_critical,
        :pike_expectation,
        :policy,
        :recycling_rate,
        :target_nutrients,
        :year,
    ]

    @test sort(collect(keys(LimnoSES.type2dict(model.policy)))) == [
        :current_term_only,
        :every,
        :max_steps,
        :max_time,
        :objectives,
        :opt_replicates,
        :start,
        :target,
        :trace_mode,
    ]

    @test sort(collect(keys(LimnoSES.type2dict(model.policy; prefix = "x")))) == [
        :x_current_term_only,
        :x_every,
        :x_max_steps,
        :x_max_time,
        :x_objectives,
        :x_opt_replicates,
        :x_start,
        :x_target,
        :x_trace_mode,
    ]

    @testset "plan" begin
        schedule = plan(Angling)
        @test collect(keys(schedule)) == [-1]
        @test length(schedule) == 1
        @test all(i isa Angling for i in Iterators.flatten(values(schedule)))
        @test count(_ -> true, Iterators.flatten(values(schedule))) == 1
        @test all(i.rate ≈ 0.000225 for i in Iterators.flatten(values(schedule)))
        schedule = plan(Angling; rate = 2.5e-3)
        @test collect(keys(schedule)) == [-1]
        @test length(schedule) == 1
        @test all(i isa Angling for i in Iterators.flatten(values(schedule)))
        @test count(_ -> true, Iterators.flatten(values(schedule))) == 1
        @test all(i.rate ≈ 0.0025 for i in Iterators.flatten(values(schedule)))
        schedule = plan(Angling, 7; rate = 3.2e-3)
        @test collect(keys(schedule)) == [7]
        @test length(schedule) == 1
        @test all(i isa Angling for i in Iterators.flatten(values(schedule)))
        @test count(_ -> true, Iterators.flatten(values(schedule))) == 1
        @test all(i.rate ≈ 0.0032 for i in Iterators.flatten(values(schedule)))
        schedule = plan(Angling, 3:5)
        @test sort(collect(keys(schedule))) == [3, 4, 5]
        @test length(schedule) == 3
        @test all(i isa Angling for i in Iterators.flatten(values(schedule)))
        @test count(_ -> true, Iterators.flatten(values(schedule))) == 3
        @test all(i.rate ≈ 0.000225 for i in Iterators.flatten(values(schedule)))
        schedule =
            plan(Angling, [(period = 1:4,), (year = 5, rate = 7.2e-3), (period = 7:9,)])
        @test sort(collect(keys(schedule))) == [1, 2, 3, 4, 5, 7, 8, 9]
        @test length(schedule) == 8
        @test all(i isa Angling for i in Iterators.flatten(values(schedule)))
        @test count(_ -> true, Iterators.flatten(values(schedule))) == 8
        @test all(
            i.rate ≈ 0.000225
            for i in Iterators.flatten(v for (k, v) in pairs(schedule) if k != 5)
        )
        @test all(
            i.rate ≈ 0.0072
            for i in Iterators.flatten(v for (k, v) in pairs(schedule) if k == 5)
        )
        # An edge case:
        @test length(plan(Angling, [(year = 1, rate = 1), (year = 2,)])) == 2
        @test length(plan(Angling, [(year = 1,), (year = 2,)])) == 2
        @test length(plan(Angling, [(year = 1,)])) == 1

        @test_throws AssertionError plan(Angling, [(rate = 1,)])
        @test_throws MethodError plan(Angling, [(year = 1,)]; rate = 1e-3)
    end
    @testset "planner" begin
        interventions = planner(plan(Angling))
        schedule = plan(Angling)
        @test collect(keys(schedule)) == collect(keys(interventions))
        @test length(schedule) == length(interventions)
        @test collect(i isa Angling for i in Iterators.flatten(values(schedule))) ==
              collect(i isa Angling for i in Iterators.flatten(values(interventions)))
        @test count(_ -> true, Iterators.flatten(values(schedule))) ==
              count(_ -> true, Iterators.flatten(values(interventions)))
        @test all(i.rate ≈ 0.000225 for i in Iterators.flatten(values(interventions)))

        interventions = planner(plan(Planting, 2; rate = 5e-3), plan(Trawling, 1:3))
        @test sort(collect(keys(interventions))) == [1, 2, 3]
        @test length(interventions) == 3
        @test count(_ -> true, Iterators.flatten(values(interventions))) == 4
        @test all(
            i isa Planting || i isa Trawling
            for i in Iterators.flatten(values(interventions))
        )
        @test count(i -> i isa Planting, Iterators.flatten(values(interventions))) == 1

        @test_throws AssertionError planner(plan(Angling), plan(Angling, 7))
    end
    @testset "scan" begin
        default_planting = (rate = (0.0, 0.01),)
        @test scan(Planting) == Dict(Planting => default_planting)
        @test LimnoSES.default_policy(Planting) == default_planting
        @test scan(Planting; rate = (0.01, 1.0)) == Dict(Planting => (rate = (0.01, 1.0),))
        default_trawling = (rate = (0.0, 0.01),)
        @test scan(Trawling) == Dict(Trawling => default_trawling)
        @test LimnoSES.default_policy(Trawling) == default_trawling
        @test scan(Trawling; rate = (0.001, 1.0)) ==
              Dict(Trawling => (rate = (0.001, 1.0),))
        default_angling = (rate = (2.25e-3, 2.7e-3),)
        @test scan(Angling) == Dict(Angling => default_angling)
        @test LimnoSES.default_policy(Angling) == default_angling
        @test scan(Angling; rate = (0.0, 1.0)) == Dict(Angling => (rate = (0.0, 1.0),))
    end
    @testset "policy" begin
        plant = scan(Planting)
        pol = policy(scan(Planting))
        @test pol == plant

        pol = policy(scan(Planting), scan(Angling))
        @test all(map(k -> k in [Angling, Planting], collect(keys(pol))))
    end
    @testset "objective" begin
        obj = objective(min_time)
        @test typeof(obj) <: Tuple{<:Function,Float64}
        @test last(obj) == 1.0

        obj = objective(min_time, 5)
        @test typeof(obj) <: Tuple{<:Function,Int}
        @test last(obj) == 5

        obj = objective(min_time, 1.6)
        @test typeof(obj) <: Tuple{<:Function,Float64}
        @test last(obj) == 1.6
    end
    @testset "objectives" begin
        obj = objective(min_time)
        objs = objectives(objective(min_time))
        @test objs[1] == obj

        objs = objectives(objective(min_time), objective(min_cost))
        @test typeof(objs) <: NTuple{N,Tuple{<:Function,Float64}} where {N}
        @test length(objs) == 2
        @test sum(last.(objs)) == 1.0

        objs = objectives(objective(min_time, 57), objective(min_cost, 1.5))
        @test typeof(objs) <: NTuple{N,Tuple{<:Function,Float64}} where {N}
        @test sum(last.(objs)) == 1.0
    end
end


text:
# Assigning sources to threads based on which sources overlap others.
#

function union_find!(i, components_tree)
    root = i
    while components_tree[i] != i
        i = components_tree[i]
    end

    while root != i
        root2 = components_tree[root]
        components_tree[root] = i
        root = root2
    end
    return i
end

"""
Computes connected components given a conflict graph, and the range of sources
to process. Writes to the components argument.

- source_start - start source
- end_source - end source (inclusively processed)
- sources - the actual source ids to find connected components for
    (indexed by source_start and source_end)
- neighbor_map - conflict graph of sources
- components_tree - the components array to do union find on.
- components - Dict{Int64, Vector{Int64}} dictionary from cc_id to
              target_source index in that component
- src_to_target_index - a mapping from light source id -> index
            within the target_sources array. This is required
            since we need to write the index within target sources
            to the output, rather than the source id itself.
Returns:
- Nothing
"""
function compute_connected_components(source_start, source_end, sources, neighbor_map,
                                      components_tree, components, src_to_target_index)
    # Construct reverse map from source_id -> index in source. Only do this for
    # the specified portion we are computing connected components.
    source_to_index = Dict{Int64, Int64}()
    index_to_source = Dict{Int64, Int64}()
    for i = source_start:source_end
        source_to_index[sources[i]] = i-source_start+1
        index_to_source[i-source_start+1] = sources[i]
    end

    for i = source_start:source_end
        components_tree[i-source_start+1] = i-source_start+1
    end

    # Run union find algorithm to find connected components.
    for i = source_start:source_end
        target = union_find!(i-source_start+1, components_tree)
        for neighbor in neighbor_map[index_to_source[i-source_start+1]]
            if haskey(source_to_index, neighbor)
                neighbor_idx = source_to_index[neighbor]
                component_source = union_find!(neighbor_idx, components_tree)
                components_tree[component_source] = target
            end

        end
    end

    # Write to components dictionary.
    for i = source_start:source_end
        index = union_find!(i-source_start+1, components_tree)
        if !haskey(components, index)
            components[index] = Vector{Int64}()
        end
        push!(components[index], src_to_target_index[index_to_source[i-source_start+1]])
    end
end

"""
Partitions sources via the cyclades algorithm.

- n_threads - number of threads to which to distribute sources
- target_sources - array of target sources. Elements should match keys of neighbor_map.
- neighbor_map - graph of connections of sources

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources(indices)])
"""
function partition_cyclades(n_threads, target_sources, neighbor_map; batch_size=60)
    #Log.info("Starting Cyclades partitioning...")
    tic()

    n_sources = length(target_sources)
    n_total_batches = convert(Int64, ceil(n_sources / batch_size))

    # Construct src_to_indx map
    src_to_indx = Dict{Int64, Int64}()
    for i=1:n_sources
        src_to_indx[target_sources[i]] = i
    end

    # The final workload distribution.
    thread_sources_assignment = Vector{Vector{Vector{Int64}}}(n_threads)
    for thread = 1:n_threads
        thread_sources_assignment[thread] = Vector{Vector{Int64}}(n_total_batches)
        for batch = 1:n_total_batches
            thread_sources_assignment[thread][batch] = Vector{Int64}()
        end
    end

    # First shuffle the sources. Note Cyclades is serially equivalent
    # to this permutation of sources.
    sources = [x for x in keys(neighbor_map)]
    shuffle!(sources)

    # Process batch_size number of sources at a time.
    # TODO max - parallelize everything below (particularly CC computation)

    # The components tree is for union-find.
    components_tree = [Vector{Int64}(batch_size) for i=1:n_threads]

    # We have n_total_batches components, where each component is a dictionary
    # from the component id, to a list of sources in that component
    components = [Dict{Int64, Vector{Int64}}() for i=1:n_total_batches]
    sources_in_components = 0
    for (cur_batch, source_index) in enumerate(collect(1:batch_size:n_sources))
        source_start = source_index
        source_end = min(source_index+batch_size-1, n_sources)
        @assert source_end - source_start + 1 <= batch_size

        # TODO max - parallelize this
        compute_connected_components(source_start, source_end, sources, neighbor_map,
                                     components_tree[1],
                                     components[cur_batch],
                                     src_to_indx)
    end

    assigned_sources = 0
    # Load balance the connected components within each batch into thread_sources_assignment.
    for (cur_batch, cur_batch_component) in enumerate(components)
        # Priority queue for load balancing.
        pqueue = PriorityQueue(Base.Order.Forward, (i, 0) for i=1:n_threads)

        # Assign non-conflicting group of sources to different threads
        for (component_group_id, sources_of_component) in cur_batch_component
            least_loaded_thread = peek(pqueue)[1]
            for source_to_assign in sources_of_component
                push!(thread_sources_assignment[least_loaded_thread][cur_batch], source_to_assign)
                assigned_sources += 1
            end
            pqueue[least_loaded_thread] += length(sources_of_component)
        end
    end

    #Log.info("Cyclades - Assigned sources: $(assigned_sources) vs correct number of sources: $(n_sources)")
    @assert assigned_sources == n_sources
    #Log.info("Cyclades - Number of batches: $(n_total_batches)")
    #Log.info("Finished Cyclades partitioning.  Elapsed time: $(toq()) seconds")

    for cur_batch = 1:length(collect(1:batch_size:n_sources))
        load_balance_for_batch = [length(thread_sources_assignment[t][cur_batch]) for t=1:n_threads]
        #Log.info("Load balance for batch $(cur_batch) - $(load_balance_for_batch)")
    end

    thread_sources_assignment
end

"""
Partitions sources via the cyclades algorithm.

- target_sources - array of target sources. Elements should match keys of neighbor_map.
- neighbor_map - graph of connections of sources

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources(indices)])
"""
function partition_cyclades_dynamic(target_sources, neighbor_map; batch_size=60)
    #Log.info("Starting Cyclades partitioning...")
    tic()

    n_sources = length(target_sources)
    n_total_batches = convert(Int64, ceil(n_sources / batch_size))

    # Construct src_to_indx map
    src_to_indx = Dict{Int64, Int64}()
    for i=1:n_sources
        src_to_indx[target_sources[i]] = i
    end

    # The final workload distribution.
    thread_sources_assignment = Vector{Vector{Vector{Int64}}}(n_total_batches)
    for batch = 1:n_total_batches
        thread_sources_assignment[batch] = Vector{Vector{Int64}}()
    end

    # First shuffle the sources. Note Cyclades is serially equivalent
    # to this permutation of sources.
    sources = [x for x in keys(neighbor_map)]
    shuffle!(sources)

    # Process batch_size number of sources at a time.
    # TODO max - parallelize everything below (particularly CC computation)

    # The components tree is for union-find.
    components_tree = Vector{Int64}(batch_size)

    # We have n_total_batches components, where each component is a dictionary
    # from the component id, to a list of sources in that component
    components = [Dict{Int64, Vector{Int64}}() for i=1:n_total_batches]
    sources_in_components = 0
    for (cur_batch, source_index) in enumerate(collect(1:batch_size:n_sources))
        source_start = source_index
        source_end = min(source_index+batch_size-1, n_sources)
        @assert source_end - source_start + 1 <= batch_size

        # TODO max - parallelize this
        compute_connected_components(source_start, source_end, sources, neighbor_map,
                                     components_tree,
                                     components[cur_batch],
                                     src_to_indx)
    end

    assigned_sources = 0
    # Add sources to sources assignment
    for (cur_batch, cur_batch_component) in enumerate(components)
        for (component_group_id, sources_of_component) in cur_batch_component
            push!(thread_sources_assignment[cur_batch], copy(sources_of_component))
        assigned_sources += length(sources_of_component)
        end
    end

    #Log.info("Cyclades - Assigned sources: $(assigned_sources) vs correct number of sources: $(n_sources)")
    @assert assigned_sources == n_sources
    #Log.info("Cyclades - Number of batches: $(n_total_batches)")
    #Log.info("Finished Cyclades partitioning.  Elapsed time: $(toq()) seconds")

    #Log.info("Assigned sources: $(thread_sources_assignment)")

    thread_sources_assignment
end

"""
Partitions sources across threads equally.

- n_threads - number of threads to which to distribute processing the sources.
- n_sources - the number of total sources to process.

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources]).
  Note for partition equally, there is only 1 batch.

"""
function partition_equally(n_threads, n_sources)
    #Log.info("Starting basic source partitioning...")
    tic()
    n_sources_per_thread = floor(Int64, n_sources / n_threads)
    thread_sources_assignment = Vector{Vector{Vector{Int64}}}(n_threads)
    n_sources_assigned = 0
    for thread = 1:n_threads
        start_source = (thread-1) * n_sources_per_thread
        end_source = thread * n_sources_per_thread
        if thread == n_threads
            end_source = n_sources
        end
        # Only 1 batch for partitioning equally
        thread_sources_assignment[thread] = Vector{Vector{Int64}}(1)
        thread_sources_assignment[thread][1] = Vector{Int64}(end_source-start_source)
        for source_assignment = start_source:end_source-1
            thread_sources_assignment[thread][1][source_assignment-start_source+1] = source_assignment+1
            n_sources_assigned += 1
        end
    end
    @assert n_sources_assigned == n_sources
    #Log.info("Finished basic source partitioning. Elapsed time: $(toq()) seconds")
    thread_sources_assignment
end

ext: jl


text:
struct SNAPParams{S<:Integer, T<:AbstractFloat}
    n_atoms    :: S
    twojmax    :: S
    elements   :: Vector{Symbol}

    rcutfac    :: T
    rmin0      :: T
    rfac0      :: T
    radii      :: Vector{T}
    weight     :: Vector{T}
    chem_flag  :: Bool
    bzero_flag :: Bool
    bnorm_flag :: Bool
    switch_flag :: Bool
    wselfall_flag :: Bool
    prebuilt_flag :: Bool
    prebuilt_arrays :: PrebuiltArrays{T}
end

function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T},
    chem_flag::Bool, bzero_flag::Bool, bnorm_flag::Bool,
    switch_flag::Bool, prebuilt_flag::Bool, prebuilt_arrays::PrebuiltArrays{T} ) where {S<:Integer,T<:AbstractFloat}
    if chem_flag != bnorm_flag
        error("chem_flag and bnorm_flag should be the same, setting bnorm_flag == chem_flag")
        bnorm_flag = chem_flag
    end
    return SNAPParams{S, T}(n_atoms, twojmax, elements,
        rcutfac, rmin0, rfac0, radii, weight, chem_flag, bzero_flag, bnorm_flag,
        switch_flag, wselfall_flag, true, prebuilt_arrays)

end

function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T},
    chem_flag::Bool, bzero_flag::Bool, bnorm_flag::Bool,
    switch_flag::Bool, wselfall_flag::Bool, prebuilt_flag::Bool ) where {S<:Integer,T<:AbstractFloat}
    if chem_flag != bnorm_flag
        error("chem_flag and bnorm_flag should be the same, setting bnorm_flag == chem_flag")
        bnorm_flag = chem_flag
    end
    if prebuilt_flag
        error("Cannot pass prebuilt_flag = true without also passing prebuilt arrays.")
    else
        return SNAPParams{S, T}(n_atoms, twojmax, elements,
        rcutfac, rmin0, rfac0, radii, weight, chem_flag, bzero_flag,
        bnorm_flag, switch_flag, wselfall_flag, true,
        initialize_prebuilt_arrays(twojmax, length(elements), chem_flag, bzero_flag, bnorm_flag))
    end
end

function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T},
        chem_flag::Bool, bzero_flag::Bool, bnorm_flag::Bool,
        switch_flag::Bool, wselfall_flag::Bool ) where {S<:Integer,T<:AbstractFloat}
        if chem_flag != bnorm_flag
            error("chem_flag and bnorm_flag should be the same, setting bnorm_flag == chem_flag")
            bnorm_flag = chem_flag
        end
        return SNAPParams{S, T}(n_atoms, twojmax, elements,
        rcutfac, rmin0, rfac0, radii, weight, chem_flag, bzero_flag, bnorm_flag,
        switch_flag, wselfall_flag, false)
end
function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T}, chem_flag::Bool,
    bzero_flag::Bool, bnorm_flag::Bool, switch_flag::Bool) where {S<:Integer,T<:AbstractFloat}
    if chem_flag != bnorm_flag
        error("chem_flag and bnorm_flag should be the same, setting bnorm_flag == chem_flag")
        bnorm_flag = chem_flag
    end
    return SNAPParams(n_atoms, twojmax, elements,
    rcutfac, rmin0, rfac0, radii, weight, chem_flag, bzero_flag, bnorm_flag,
    false, false, false)
end

function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T},
    chem_flag::Bool, bzero_flag::Bool, bnorm_flag::Bool) where {S<:Integer,T<:AbstractFloat}
    if chem_flag != bnorm_flag
        error("chem_flag and bnorm_flag should be the same, setting bnorm_flag == chem_flag")
        bnorm_flag = chem_flag
    end
    return SNAPParams(n_atoms, twojmax, elements,
    rcutfac, rmin0, rfac0, radii, weight, chem_flag, bzero_flag, bnorm_flag,
    false, false, false)
end
function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T},
    chem_flag::Bool, bzero_flag::Bool) where {S<:Integer,T<:AbstractFloat}

    # setting bzero_flag == chem_flag
    bzero_flag = chem_flag
    return SNAPParams(n_atoms, twojmax, elements,
    rcutfac, rmin0, rfac0, radii, weight, chem_flag, bzero_flag, bzero_flag,
    false, false, false)
end

function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T}, chem_flag::Bool) where {S<:Integer,T<:AbstractFloat}
    return SNAPParams(n_atoms, twojmax, elements,
    rcutfac, rmin0, rfac0, radii, weight, chem_flag, false, bnorm_flag,
    false, false, false)
end

function SNAPParams(n_atoms::S, twojmax::S, elements::Vector{Symbol},
    rcutfac::T, rmin0::T, rfac0::T, radii::Vector{T}, weight ::Vector{T},) where {S<:Integer,T<:AbstractFloat}
    return SNAPParams(n_atoms, twojmax, elements,
    rcutfac, rmin0, rfac0, radii, weight, false, false, false,
    false, false, false)
end


text:
######## GAS DEFINITIONS #######

export AbstractGas, PerfectGas

abstract type AbstractGas <: AbstractFluid end

# default values
const DefaultGasRefTemperature = Temperature(15u"°C")
const DefaultGasRefPressure = Pressure(1u"atm")
const DefaultGasConstant = GasConstant(287)
const DefaultSpecificHeatRatio = SpecificHeatRatio(1.4)
const DefaultGasViscosity = Viscosity(1.79e-5)

SpecificHeatPressure(;γ::SpecificHeatRatio=DefaultSpecificHeatRatio,
                      R::GasConstant=DefaultGasConstant) = SpecificHeatPressure(γ*R/(γ-1))
SpecificHeatVolume(;γ::SpecificHeatRatio=DefaultSpecificHeatRatio,
                    R::GasConstant=DefaultGasConstant) = SpecificHeatVolume(R/(γ-1))

struct PerfectGas <: AbstractGas
    Tref :: Temperature
    pref :: Pressure
    ρref :: Density
    specwt :: SpecificWeight
    μ :: Viscosity
    ν :: KinematicViscosity
    γ :: SpecificHeatRatio
    R :: GasConstant
    cp :: SpecificHeatPressure
    cv :: SpecificHeatVolume
end

function PerfectGas(;
            Tref::Temperature=DefaultGasRefTemperature,
            pref::Pressure=DefaultGasRefPressure,
            μ::Viscosity=DefaultGasViscosity,
            γ::SpecificHeatRatio=DefaultSpecificHeatRatio,
            R::GasConstant=DefaultGasConstant)

    ρref = Density(pref/Tref/R)
    return PerfectGas(Tref,pref,ρref,SpecificWeight(ρref*1u"ge"),μ,KinematicViscosity(μ/ρref),γ,R,SpecificHeatPressure(γ=γ,R=R),SpecificHeatVolume(γ=γ,R=R))
end

Density(g::PerfectGas) = g.ρref
Viscosity(g::PerfectGas) = g.μ
KinematicViscosity(g::PerfectGas) = g.ν
SpecificWeight(g::PerfectGas) = g.specwt
SpecificHeatRatio(g::PerfectGas) = g.γ
GasConstant(g::PerfectGas) = g.R
SpecificHeatPressure(g::PerfectGas) = g.cp
SpecificHeatVolume(g::PerfectGas) = g.cv

function Base.show(io::IO, m::MIME"text/plain", g::PerfectGas)
    println(io,"Perfect gas with")
    println(io,"   Density = $(value(Density(g)))")
    println(io,"   Viscosity = $(value(Viscosity(g)))")
    println(io,"   Specific heat ratio = $(value(SpecificHeatRatio(g)))")
    println(io,"   Gas constant = $(value(GasConstant(g)))")
    println(io,"   cp = $(value(SpecificHeatPressure(g)))")
    println(io,"   cv = $(value(SpecificHeatVolume(g)))")
    println(io,"   at reference temp $(value(g.Tref))")
end

max_stars_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_issues_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_forks_repo_name: UCLAMAEThreads/ThermofluidQuantities


text:
using myJuliaUtils
using Test
using LinearAlgebra

N = rand(1:100)
A = genPDmat(N)
cholA = cholPSD(A)
@test cholA == cholesky(A).U
@test cholA'cholA ≈ A

# test1=0
# test2=0
# NoTimes=20
# for i=1:NoTimes
# N = round(Int, 1 + rand() * 100)
# A = randn(N, N)
# A = A'A
# cholA = cholPSD(A)
# test1+=cholA == cholesky(A).U
# test2+=cholA'cholA ≈ A
# end
# @test test1==NoTimes
# @test test2==NoTimes

N = rand(2:100)
K = rand(1:N - 1)
A = genPDmat(N)
S, Q = eigen(A)
# S[1:K] = ones(K) * eps(0.0)
S[1:K] = zeros(K)
# A = Q * Diagonal(S) * inv(Q)
A = Q * Diagonal(S)
A[:,1:K] = zeros(N, K)
A .= A / Q
cholA = cholPSD(A)
if K > 1
    @test size(cholA) == (N - K, N)
else
    @test all(size(cholA) .>= (N - K, N)) # this is done to avoid approximation error
end

@test cholA'cholA ≈ A

# Nseed=round(Int, 1 + rand() * 5000)
# srand(Nseed)
# a=randn(1, N) * cholx(A)
# srand(Nseed)
# a=randn(1, N) * cholx(A)

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
# global state management
#
# The functionality in this file serves to create a coherent environment to perform
# computations in, with support for Julia constructs like tasks (and executing those on
# multiple threads), using a GPU of your choice (with the ability to reset that device, or
# use different devices on different threads).

export context, context!, device, device!, device_reset!, deviceid

## hooks

"""
    CUDA.attaskswitch(f::Function)

Register a function to be called after switching to or initializing a task on a thread.

Use this hook to invalidate thread-local state that depends on the current task.
"""
attaskswitch(f::Function) = (pushfirst!(task_hooks, f); nothing)
const task_hooks = []
_attaskswitch() = foreach(f->Base.invokelatest(f), task_hooks)

"""
    CUDA.atdeviceswitch(f::Function)

Register a function to be called after switching to or initializing a device on a thread.

Use this hook to invalidate thread-local state that depends on the current device. If that
state is also context dependent, be sure to query the context in your callback.
"""
atdeviceswitch(f::Function) = (pushfirst!(device_switch_hooks, f); nothing)
const device_switch_hooks = []
_atdeviceswitch() = foreach(f->Base.invokelatest(f), device_switch_hooks)

"""
    CUDA.atdevicereset(f::Function)

Register a function to be called after resetting devices. The function is passed one
argument: the device which has been reset.

Use this hook to invalidate global state that depends on the current device.
"""
atdevicereset(f::Function) = (pushfirst!(device_reset_hooks, f); nothing)
const device_reset_hooks = []
_atdevicereset(dev) = foreach(f->Base.invokelatest(f, dev), device_reset_hooks)

## initialization

"""
    CUDA.prepare_cuda_call()

Prepare state for calling CUDA API functions.

Many CUDA APIs, like the CUDA driver API used by CUDA.jl, use implicit thread-local state
to determine, e.g., which device to use. With Julia however, code is grouped in tasks.
Execution can switch between them, and tasks can be executing on (and in the future migrate
between) different threads. To synchronize these two worlds, call this function before any
CUDA API call to update thread-local state based on the current task and its context.

If you need to maintain your own task-local state, subscribe to device and task switch
events using [`CUDA.atdeviceswitch`](@ref) and [`CUDA.attaskswitch`](@ref) for
proper invalidation. If your state is device-specific, but global (i.e. not task-bound), it
suffices to index your state with the current [`deviceid()`](@ref) and invalidate that state
when the device is reset by subscribing to [`CUDA.atdevicereset()`](@ref).
"""
@inline function prepare_cuda_call()
    tid = Threads.threadid()

    # detect when a different task is now executing on a thread
    task = @inbounds thread_tasks[tid]
    if task === nothing || task.value === nothing || task.value::Task !== current_task()
        switched_tasks(tid, current_task())
    end

    # initialize CUDA state when first executing on a thread
    state = @inbounds thread_state[tid]
    if state === nothing
        initialize_thread(tid)
    end

    # FIXME: this is expensive. Maybe kernels should return a `wait`able object, a la KA.jl,
    #        which then performs the necessary checks. Or only check when launching kernels.
    check_exceptions()

    return
end

# Julia executes with tasks, so we need to keep track of the active task for each thread
# in order to detect task switches and update the thread-local state accordingly.
# doing so using task_local_storage is too expensive.
const thread_tasks = Union{Nothing,WeakRef}[]
@noinline function switched_tasks(tid::Int, task::Task)
    thread_tasks[tid] = WeakRef(task)
    _attaskswitch()

    # switch contexts if task switched to was already bound to one
    ctx = get(task_local_storage(), :CuContext, nothing)
    if ctx !== nothing && isvalid(ctx)
        # NOTE: the context may be invalid if another task reset it (which we detect here
        #       since we can't touch other tasks' local state from `device_reset!`)
        context!(ctx)
    else
        thread_state[tid] = nothing  # trigger `initialize_thread` down the line
        # NOTE: actually deactivating the CUDA context would be more correct,
        #       but that confuses CUDA and leads to invalid contexts later on.
    end
end

# the default device unitialized tasks will use, set when switching devices.
# this behavior differs from the CUDA Runtime, where device 0 is always used.
# this setting won't be used when switching tasks on a pre-initialized thread.
const default_device = Ref{Union{Nothing,CuDevice}}(nothing)

# CUDA uses thread-bound state, but calling CuCurrent* all the time is expensive,
# so we maintain our own thread-local copy keeping track of the current CUDA state.
const CuCurrentState = NamedTuple{(:ctx, :dev), Tuple{CuContext,CuDevice}}
const thread_state = Union{Nothing,CuCurrentState}[]
@noinline function initialize_thread(tid::Int)
    dev = something(default_device[], CuDevice(0))
    device!(dev)

    # NOTE: we can't be compatible with externally initialize contexts here (i.e., reuse
    #       the CuCurrentContext we don't know about) because of how device_reset! works:
    #       contexts that got reset remain bound, and are not discernible from regular ones.
end

## context-based API

"""
    context()::CuContext

Get or create a CUDA context for the current thread (as opposed to
`CuCurrentContext` which may return `nothing` if there is no context bound to the
current thread).
"""
@inline function context()
    tid = Threads.threadid()

    prepare_cuda_call()
    state = @inbounds thread_state[tid]::CuCurrentState

    if Base.JLOptions().debug_level >= 2
        @assert state.ctx == CuCurrentContext()
    end
    state.ctx
end

"""
    context!(ctx::CuContext)

Bind the current host thread to the context `ctx`.

Note that the contexts used with this call should be previously acquired by calling
[`context`](@ref), and not arbitrary contexts created by calling the `CuContext` constructor.
"""
function context!(ctx::CuContext)
    # update the thread-local state
    tid = Threads.threadid()
    state = @inbounds thread_state[tid]
    if state === nothing || state.ctx != ctx
        activate(ctx)
        dev = CuCurrentDevice()
        thread_state[tid] = (;ctx=ctx, dev=dev)
        _atdeviceswitch()
    end

    # update the task-local state
    task_local_storage(:CuContext, ctx)

    return
end

"""
    context!(f, ctx)

Sets the active context for the duration of `f`.
"""
function context!(f::Function, ctx::CuContext)
    @assert isvalid(ctx)
    old_ctx = CuCurrentContext()
    if ctx != old_ctx
        context!(ctx)
    end
    try
        f()
    finally
        if ctx !== old_ctx && old_ctx !== nothing
            context!(old_ctx)
        end
    end
end

# macro version for maximal performance (avoiding closures)
macro context!(ctx_expr, expr)
    quote
        ctx = $(esc(ctx_expr))
        @assert isvalid(ctx)
        old_ctx = CuCurrentContext()
        if ctx != old_ctx
            context!(ctx)
        end
        try
            $(esc(expr))
        finally
            if ctx !== old_ctx && old_ctx !== nothing
                context!(old_ctx)
            end
        end
    end
end

## device-based API

"""
    device()::CuDevice

Get the CUDA device for the current thread, similar to how [`context()`](@ref) works
compared to [`CuCurrentContext()`](@ref).
"""
@inline function device()
    tid = Threads.threadid()

    prepare_cuda_call()
    state = @inbounds thread_state[tid]::CuCurrentState

    if Base.JLOptions().debug_level >= 2
        @assert state.dev == CuCurrentDevice()
    end
    state.dev
end

const device_contexts = Union{Nothing,CuContext}[]
function context(dev::CuDevice)
    tid = Threads.threadid()
    devidx = deviceid(dev)+1

    # querying the primary context for a device is expensive, so cache it
    ctx = @inbounds device_contexts[devidx]
    if ctx !== nothing
        return ctx
    end

    # configure the primary context
    pctx = CuPrimaryContext(dev)
    ctx = CuContext(pctx)
    @inbounds device_contexts[devidx] = ctx
    return ctx
end

"""
    device!(dev::Integer)
    device!(dev::CuDevice)

Sets `dev` as the current active device for the calling host thread. Devices can be
specified by integer id, or as a `CuDevice` (slightly faster).

If your library or code needs to perform an action when the active device changes,
add a hook using [`CUDA.atdeviceswitch`](@ref).
"""
function device!(dev::CuDevice, flags=nothing)
    tid = Threads.threadid()
    devidx = deviceid(dev)+1

    # configure the primary context flags
    if flags !== nothing
        if @inbounds device_contexts[devidx] !== nothing
            error("Cannot set flags for an active device. Do so before calling any CUDA function, or reset the device first.")
        end
        cuDevicePrimaryCtxSetFlags(dev, flags)
    end

    # make this device the new default
    default_device[] = dev

    # bail out if switching to the current device
    state = @inbounds thread_state[tid]
    if state !== nothing && state.dev == dev
        return
    end

    # actually switch contexts
    ctx = context(dev)
    context!(ctx)
end

"""
    device!(f, dev)

Sets the active device for the duration of `f`.

Note that this call is intended for temporarily switching devices, and does not change the
default device used to initialize new threads or tasks.
"""
function device!(f::Function, dev::CuDevice)
    ctx = context(dev)
    context!(f, ctx)
end

macro device!(dev_expr, expr)
    quote
        dev = $(esc(dev_expr))
        ctx = context(dev)
        @context! ctx $(esc(expr))
    end
end

"""
    device_reset!(dev::CuDevice=device())

Reset the CUDA state associated with a device. This call with release the underlying
context, at which point any objects allocated in that context will be invalidated.

If your library or code needs to perform an action when the active context changes,
add a hook using [`CUDA.atdevicereset`](@ref). Resetting the device will also cause
subsequent API calls to fire the [`CUDA.atdeviceswitch`](@ref) hook.
"""
function device_reset!(dev::CuDevice=device())
    # unconditionally reset the primary context (don't just release it),
    # as there might be users outside of CUDA.jl
    pctx = CuPrimaryContext(dev)
    unsafe_reset!(pctx)

    # wipe the thread-local state for all threads using this device
    for (tid, state) in enumerate(thread_state)
        if state !== nothing && state.dev == dev
            thread_state[tid] = nothing
        end
    end

    # wipe the device-specific state
    devidx = deviceid(dev)+1
    device_contexts[devidx] = nothing

    _atdevicereset(dev)

    return
end

## integer device-based API

device!(dev::Integer, flags=nothing) = device!(CuDevice(dev), flags)
device!(f::Function, dev::Integer) = device!(f, CuDevice(dev))

"""
    deviceid()::Int
    deviceid(dev::CuDevice)::Int

Get the ID number of the current device of execution. This is a 0-indexed number,
corresponding to the device ID as known to CUDA.
"""
deviceid(dev::CuDevice=device()) = Int(convert(CUdevice, dev))

## helpers

# helper struct to maintain state per device
# - make it possible to index directly with CuDevice (without converting to integer index)
# - initialize function to fill state based on a constructor function
# - automatically wiping state on device reset
struct PerDevice{T,F} <: AbstractVector{T}
    inner::Vector{T}
    ctor::F

    PerDevice{T,F}(ctor::F) where {T,F<:Function} = new(Vector{T}(), ctor)
end

PerDevice{T}(ctor::F) where {T,F} = PerDevice{T,F}(ctor)

function initialize!(x::PerDevice, n::Integer)
    @assert isempty(x.inner)

    resize!(x.inner, n)
    for i in 0:n-1
        x[i] = x.ctor(CuDevice(i))
    end

    atdevicereset() do dev
        x[dev] = x.ctor(dev)
    end

    return
end

# 0-based indexing for using CUDA device identifiers
Base.getindex(x::PerDevice, devidx::Integer) = x.inner[devidx+1]
Base.setindex!(x::PerDevice, val, devidx::Integer) = (x.inner[devidx+1] = val; )

# indexing using CuDevice objecs
Base.getindex(x::PerDevice, dev::CuDevice) = x[deviceid(dev)]
Base.setindex!(x::PerDevice, val, dev::CuDevice) = (x[deviceid(dev)] = val; )

Base.length(x::PerDevice) = length(x.inner)
Base.size(x::PerDevice) = size(x.inner)

function Base.show(io::IO, mime::MIME"text/plain", x::PerDevice{T}) where {T}
    print(io, "PerDevice{$T} with $(length(x)) entries")
end


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

JLLWrappers.@generate_main_file_header("Clipper")
JLLWrappers.@generate_main_file("Clipper", UUID("1721f0f4-5627-55cb-8b31-c466f04189fe"))
end  # module Clipper_jll


text:
"""
    Schedule(ts, values, interpolations; unit = :epoch; kwargs...)

Describes how the values of a hyperparameter change over the
training.

## Examples

```julia
# one-cycle schedule over 10 epochs from 0.01 over 0.1 to 0.001
Schedule([0, 2, 10], [0.01, 0.1, 0.001], [sineio(), sineio()])
```

`Schedule` is a wrapper around `Animations.Animation`, see the
[documentation](https://jkrumbiegel.github.io/Animations.jl/dev)
for more detailed information.
"""
struct Schedule{T}
    animation::Animation{T}
    unit::Symbol
end

Schedule(args...; unit = :epoch, kwargs...) =
    Schedule(Animation(args...; kwargs...), unit)

struct Scheduler <: Callback
    schedules::Dict{Type{<:HyperParameter}, Schedule}
end

Base.show(io::IO, scheduler::Scheduler) = print(
    io, "Scheduler(", join(keys(scheduler.schedules), ", "), ")")

function stateaccess(scheduler::Scheduler)
    # TODO: implement proper merging of permissions
    if length(keys(scheduler.schedules)) == 0
        hpstateaccess = (;)
    else
        hpstateaccess = merge(stateaccess.(keys(scheduler.schedules))...)
    end
    return (
        data = Read(), cbstate = (history = Read(), hyperparams = Write()),
        hpstateaccess...
    )
end

function on(::Init, phase, scheduler::Scheduler, learner)
    learner.cbstate.hyperparams = MVHistory()
    # convert schedules to use step as a unit
    epochlength = length(learner.data[1])
    for (H, schedule) in scheduler.schedules
        if schedule.unit == :epoch
            scheduler.schedules[H] = Schedule(schedule.animation * epochlength, :step)
        end
    end
end

function on(::BatchBegin, ::AbstractTrainingPhase, scheduler::Scheduler, learner)
    step = learner.cbstate.history.steps
    for (H, schedule) in scheduler.schedules
        value = Animations.at(schedule.animation, step)
        sethyperparameter!(learner, H, Animations.at(schedule.animation, step))
        push!(learner.cbstate.hyperparams, Symbol(H), step, value)
    end
end

function onecycle(
        epochs, max_lr,
        start_lr = max_lr / 10,
        end_lr = max_lr / 30,
        start_pctg = 0.1)
    return Schedule(
        [0, epochs * start_pctg, epochs],
        [start_lr, max_lr, end_lr],
        [Animations.sineio(), Animations.sineio()]
    )

end


text:
# routines that implement different preconditioners

export ic02!, ic02, ilu02!, ilu02

"""
    ic02!(A::CuSparseMatrix, index::SparseChar)

Incomplete Cholesky factorization with no pivoting.
Preserves the sparse layout of matrix `A`.
"""
ic02!(A::CuSparseMatrix, index::SparseChar)
for (bname,aname,sname,elty) in ((:cusparseScsric02_bufferSize, :cusparseScsric02_analysis, :cusparseScsric02, :Float32),
                                 (:cusparseDcsric02_bufferSize, :cusparseDcsric02_analysis, :cusparseDcsric02, :Float64),
                                 (:cusparseCcsric02_bufferSize, :cusparseCcsric02_analysis, :cusparseCcsric02, :ComplexF32),
                                 (:cusparseZcsric02_bufferSize, :cusparseZcsric02_analysis, :cusparseZcsric02, :ComplexF64))
    @eval begin
        function ic02!(A::CuSparseMatrixCSR{$elty},
                       index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            info = csric02Info_t[0]
            cusparseCreateCsric02Info(info)
            @workspace size=@argout(
                    $bname(handle(), m, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), m, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXcsric02_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural/numerical zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), m, A.nnz,
                           desc, A.nzVal, A.rowPtr, A.colVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyCsric02Info(info[1])
            A
        end
    end
end

# cscic02
for (bname,aname,sname,elty) in ((:cusparseScsric02_bufferSize, :cusparseScsric02_analysis, :cusparseScsric02, :Float32),
                                 (:cusparseDcsric02_bufferSize, :cusparseDcsric02_analysis, :cusparseDcsric02, :Float64),
                                 (:cusparseCcsric02_bufferSize, :cusparseCcsric02_analysis, :cusparseCcsric02, :ComplexF32),
                                 (:cusparseZcsric02_bufferSize, :cusparseZcsric02_analysis, :cusparseZcsric02, :ComplexF64))
    @eval begin
        function ic02!(A::CuSparseMatrixCSC{$elty},
                       index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            info = csric02Info_t[0]
            cusparseCreateCsric02Info(info)
            @workspace size=@argout(
                    $bname(handle(), m, A.nnz, desc,
                           A.nzVal, A.colPtr, A.rowVal, info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), m, A.nnz, desc,
                           A.nzVal, A.colPtr, A.rowVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXcsric02_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural/numerical zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), m, A.nnz,
                           desc, A.nzVal, A.colPtr, A.rowVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyCsric02Info(info[1])
            A
        end
    end
end

"""
    ilu02!(A::CuSparseMatrix, index::SparseChar)

Incomplete LU factorization with no pivoting.
Preserves the sparse layout of matrix `A`.
"""
ilu02!(A::CuSparseMatrix, index::SparseChar)
for (bname,aname,sname,elty) in ((:cusparseScsrilu02_bufferSize, :cusparseScsrilu02_analysis, :cusparseScsrilu02, :Float32),
                                 (:cusparseDcsrilu02_bufferSize, :cusparseDcsrilu02_analysis, :cusparseDcsrilu02, :Float64),
                                 (:cusparseCcsrilu02_bufferSize, :cusparseCcsrilu02_analysis, :cusparseCcsrilu02, :ComplexF32),
                                 (:cusparseZcsrilu02_bufferSize, :cusparseZcsrilu02_analysis, :cusparseZcsrilu02, :ComplexF64))
    @eval begin
        function ilu02!(A::CuSparseMatrixCSR{$elty},
                        index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            info = csrilu02Info_t[0]
            cusparseCreateCsrilu02Info(info)
            @workspace size=@argout(
                    $bname(handle(), m, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), m, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXcsrilu02_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), m, A.nnz,
                           desc, A.nzVal, A.rowPtr, A.colVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyCsrilu02Info(info[1])
            A
        end
    end
end

# cscilu02
for (bname,aname,sname,elty) in ((:cusparseScsrilu02_bufferSize, :cusparseScsrilu02_analysis, :cusparseScsrilu02, :Float32),
                                 (:cusparseDcsrilu02_bufferSize, :cusparseDcsrilu02_analysis, :cusparseDcsrilu02, :Float64),
                                 (:cusparseCcsrilu02_bufferSize, :cusparseCcsrilu02_analysis, :cusparseCcsrilu02, :ComplexF32),
                                 (:cusparseZcsrilu02_bufferSize, :cusparseZcsrilu02_analysis, :cusparseZcsrilu02, :ComplexF64))
    @eval begin
        function ilu02!(A::CuSparseMatrixCSC{$elty},
                        index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            info = csrilu02Info_t[0]
            cusparseCreateCsrilu02Info(info)
            @workspace size=@argout(
                    $bname(handle(), m, A.nnz, desc,
                           A.nzVal, A.colPtr, A.rowVal, info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), m, A.nnz, desc,
                           A.nzVal, A.colPtr, A.rowVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXcsrilu02_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), m, A.nnz,
                           desc, A.nzVal, A.colPtr, A.rowVal, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyCsrilu02Info(info[1])
            A
        end
    end
end

# bsric02
for (bname,aname,sname,elty) in ((:cusparseSbsric02_bufferSize, :cusparseSbsric02_analysis, :cusparseSbsric02, :Float32),
                                 (:cusparseDbsric02_bufferSize, :cusparseDbsric02_analysis, :cusparseDbsric02, :Float64),
                                 (:cusparseCbsric02_bufferSize, :cusparseCbsric02_analysis, :cusparseCbsric02, :ComplexF32),
                                 (:cusparseZbsric02_bufferSize, :cusparseZbsric02_analysis, :cusparseZbsric02, :ComplexF64))
    @eval begin
        function ic02!(A::CuSparseMatrixBSR{$elty},
                       index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_UPPER, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            mb = div(m,A.blockDim)
            info = bsric02Info_t[0]
            cusparseCreateBsric02Info(info)
            @workspace size=@argout(
                    $bname(handle(), A.dir, mb, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, A.blockDim, info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), A.dir, mb, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, A.blockDim, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXbsric02_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural/numerical zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), A.dir, mb, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, A.blockDim, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyBsric02Info(info[1])
            A
        end
    end
end

# bsrilu02
for (bname,aname,sname,elty) in ((:cusparseSbsrilu02_bufferSize, :cusparseSbsrilu02_analysis, :cusparseSbsrilu02, :Float32),
                                 (:cusparseDbsrilu02_bufferSize, :cusparseDbsrilu02_analysis, :cusparseDbsrilu02, :Float64),
                                 (:cusparseCbsrilu02_bufferSize, :cusparseCbsrilu02_analysis, :cusparseCbsrilu02, :ComplexF32),
                                 (:cusparseZbsrilu02_bufferSize, :cusparseZbsrilu02_analysis, :cusparseZbsrilu02, :ComplexF64))
    @eval begin
        function ilu02!(A::CuSparseMatrixBSR{$elty},
                        index::SparseChar)
            desc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_UPPER, CUSPARSE_DIAG_TYPE_NON_UNIT, index)
            m,n = A.dims
            if m != n
                throw(DimensionMismatch("A must be square, but has dimensions ($m,$n)!"))
            end
            mb = div(m,A.blockDim)
            info = bsrilu02Info_t[0]
            cusparseCreateBsrilu02Info(info)
            @workspace size=@argout(
                    $bname(handle(), A.dir, mb, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, A.blockDim, info[1],
                           out(Ref{Cint}(1)))
                )[] buffer->begin
                    $aname(handle(), A.dir, mb, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, A.blockDim, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                    posit = Ref{Cint}(1)
                    cusparseXbsrilu02_zeroPivot(handle(), info[1], posit)
                    if posit[] >= 0
                        error("Structural/numerical zero in A at ($(posit[]),$(posit[])))")
                    end
                    $sname(handle(), A.dir, mb, A.nnz, desc,
                           A.nzVal, A.rowPtr, A.colVal, A.blockDim, info[1],
                           CUSPARSE_SOLVE_POLICY_USE_LEVEL, buffer)
                end
            cusparseDestroyBsrilu02Info(info[1])
            A
        end
    end
end

for elty in (:Float32, :Float64, :ComplexF32, :ComplexF64)
    @eval begin
        function ilu02(A::CuSparseMatrix{$elty},
                       index::SparseChar)
            ilu02!(copy(A),index)
        end
        function ic02(A::CuSparseMatrix{$elty},
                      index::SparseChar)
            ic02!(copy(A),index)
        end
        function ilu02(A::HermOrSym{$elty,CuSparseMatrix{$elty}},
                       index::SparseChar)
            ilu02!(copy(A.data),index)
        end
        function ic02(A::HermOrSym{$elty,CuSparseMatrix{$elty}},
                      index::SparseChar)
            ic02!(copy(A.data),index)
        end
    end
end


text:
"""

$(TYPEDEF)

Structure to contain the density values obtained from the calculation.

$(TYPEDFIELDS)

"""
@with_kw mutable struct Density
  solute::Float64 = 0.
  solvent::Float64 = 0.
  solvent_bulk::Float64 = 0.
end

function reset!(d::Density)
  d.solute = 0.
  d.solvent = 0.
  d.solvent_bulk = 0.
  return nothing
end

#function Base.show(io::IO, d::Density )
#  println(" Mean solute density: $(d.solute) ")
#  println(" Mean solvent density: $(d.solvent) ")
#  println(" Mean solvent bulk density: $(d.solvent_bulk) ")
#end


text:
function mdl_filter_ode!(device_states,
        output_ode,
        current_r,
        current_i,
        sys_Sbase,
        f0,
        device::DynInverter{C,O,VC,DC,P,LCLFilter}) where {C <: Converter,
                                                   O <: OuterControl,
                                                   VC<: VSControl,
                                                   DC<: DCSource,
                                                   P <: FrequencyEstimator}

    #Obtain external states inputs for component
    #TODO: If converter has dynamics, need to reference states:
    #external_ix = device.input_port_mapping[device.converter]
    #vcvd = device_states[external_ix[1]]
    #vcvq = device_states[external_ix[2]]
    external_ix = device.input_port_mapping[device.filter]
    δ = device_states[external_ix[1]]

    #Obtain inner variables for component
    V_tR = device.inner_vars[VR_inv_var]
    V_tI = device.inner_vars[VI_inv_var]
    vcvd = device.inner_vars[Vdcnv_var]
    vcvq = device.inner_vars[Vqcnv_var]

    #Get parameters
    ωb = 2*pi*f0
    lf = device.filter.lf
    rf = device.filter.rf
    cf = device.filter.cf
    lg = device.filter.lg
    rg = device.filter.rg
    MVABase = get_inverter_Sbase(device)
    ωg = 1.0 #TODO: create getter later

    #RI to dq transformation
    V_dq = ri_dq(δ)*[V_tR; V_tI]
    V_g = sqrt(V_tR^2 + V_tI^2)

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.filter]
    #@show local_ix

    #Define internal states for filter
    internal_states = @view device_states[local_ix]
    icvd = internal_states[1]
    icvq = internal_states[2]
     vod = internal_states[3]
     voq = internal_states[4]
     iod = internal_states[5]
     ioq = internal_states[6]

    #Inputs (control signals) - N/A

    #Compute 6 states ODEs (D'Arco EPSR122 Model)
    #Inverter Output Inductor (internal state)
    #𝜕id_c/𝜕t
    output_ode[local_ix[1]] = ( ωb/lf*vcvd
                              - ωb/lf*vod
                              - ωb*rf/lf*icvd
                              + ωb*ωg*icvq )
    #𝜕iq_c/𝜕t
    output_ode[local_ix[2]] = ( ωb/lf*vcvq
                              - ωb/lf*voq
                              - ωb*rf/lf*icvq
                              - ωb*ωg*icvd )
    #LCL Capacitor (internal state)
    #𝜕vd_o/𝜕t
    output_ode[local_ix[3]] = ( ωb/cf*icvd
                              - ωb/cf*iod  #i_gd was specified; use equivalent i_od
                              + ωb*ωg*voq )
    #𝜕vq_o/𝜕t
    output_ode[local_ix[4]] = ( ωb/cf*icvq
                              - ωb/cf*ioq  #i_gq was specified; use equivalent i_oq
                              - ωb*ωg*vod )
    #Grid Inductance (internal state)
    #𝜕id_o/𝜕t
    output_ode[local_ix[5]] = ( ωb/lg*vod
                              - ωb/lg*V_dq[2] #vgd
                              - ωb*rg/lg*iod
                              + ωb*ωg*ioq )
    #𝜕iq_o/𝜕t
    output_ode[local_ix[6]] = ( ωb/lg*voq
                              + ωb/lg*V_dq[1] #vgq
                              - ωb*rg/lg*ioq
                              - ωb*ωg*iod )

    #Update inner_vars
    device.inner_vars[Vdo_var] = vod
    device.inner_vars[Vqo_var] = voq
    #TODO: If PLL models at PCC, need to update inner vars:
    #device.inner_vars[Vdo_var] = V_dq[q::dq_ref]
    #device.inner_vars[Vqo_var] = V_dq[q::dq_ref]

    #Compute current from the generator to the grid
    I_RI = (MVABase/sys_Sbase)*dq_ri(δ)*[iod; ioq]
    #@show MVABase
    #@show sys_Sbase
    #Update current
    current_r[1] += I_RI[1]
    current_i[1] += I_RI[2]
end


text:
const uint64_t = UInt64
const UINT64_C = UInt64


text:
module KnetModules

using Knet, AutoGrad

include("core.jl")
export
    # Ctx
    ParamCtx, active_ctx, switch_ctx!, reset_ctx!, actx,
    # Param
    Param, aval, val, setval!,
    # KnetModule
    KnetModule, getgrad, params, modules,
    convert_params!, convert_buffers!,
    gpu!, cpu!, training!, testing!,
    clean_ctx!, is_ctx_clean,
    save_module, load_module

include("linear.jl")
export Linear

include("conv.jl")
export Conv, kaiming, conv_mode!

include("norm.jl")
export Norm, BatchNorm

include("func.jl")
export
    FnModule,
    Activation, ReLU, Sigmoid, Tanh,
    Pool, MaxPool, AvgPool,
    Dropout

include("container.jl")
export KnetContainer, Sequential, add!

include("rnn.jl")
export AbstractRNN, RNNTanh, RNNRelu, LSTM, GRU

include("embedding.jl")
export Embedding, EmbeddingMul, EmbeddingLookup

include("data.jl")

end #module


text:
using Documenter
using SVG
using Dates

makedocs(
    modules = [SVG],
    sitename="SVG.jl",
    authors = "Matt Lawless",
    format = Documenter.HTML(),
)

deploydocs(
    repo = "github.com/lawless-m/SVG.jl.git",
    devbranch = "master",
    push_preview = true,
)


text:
"""
A list of jsrenderables, that allows fast updating & rendering!
"""
struct DiffList
    children::Vector
    append::Observable{Vector{Any}}
    setindex::Observable{Tuple{Vector{Int}, Any}}
    delete::Observable{Vector{Int}}
    insert::Observable{Tuple{Int, Any}}
    empty::Observable{Bool}
    attributes::Dict{Symbol, Any}
    update_lock::Any
end

function DiffList(children::Vector; attributes...)
    return DiffList(
        children,
        Observable(Any[]),
        Observable{Tuple{Vector{Int}, Any}}((Int[], Any[])),
        Observable(Int[]),
        Observable{Tuple{Int64, Any}}((0, nothing)),
        Observable(false),
        Dict{Symbol, Any}(attributes),
        ReentrantLock()
    )
end

function synchronize_update(f, difflist::DiffList)
    Base.lock(difflist.update_lock)
    try
        return f()
    catch e
        rethrow(e)
    finally
        Base.unlock(difflist.update_lock)
    end
end

function Hyperscript.children(difflist::DiffList)
    synchronize_update(difflist) do
        return difflist.children
    end
end

function Base.length(difflist::DiffList)
    synchronize_update(difflist) do
        return length(difflist.children)
    end
end

function Base.isempty(difflist::DiffList)
    synchronize_update(difflist) do
        return isempty(difflist.children)
    end
end

function Base.getindex(difflist::DiffList, i)
    synchronize_update(difflist) do
        return difflist.children[i]
    end
end

"""
    connect!(observable_list::Observable{<:AbstractVector}, difflist::DiffList)
Connects an Observable list with the events (empty, append, setindex, delete, insert) from
`difflist`.
"""
function Observables.connect!(observable_list::Observable{<:AbstractVector}, difflist::DiffList)
    if observable_list[] != difflist.children
        empty!(observable_list[])
        append!(observable_list[], difflist.children)
    end
    on(difflist.empty) do bool
        empty!(observable_list[])
        observable_list[] = observable_list[]
    end

    on(difflist.append) do values
        append!(observable_list[], values)
        observable_list[] = observable_list[]
    end

    on(difflist.setindex) do (indices, values)
        observable_list[][indices] = values
        observable_list[] = observable_list[]
    end

    on(difflist.setindex) do (indices, values)
        observable_list[][indices] = values
        observable_list[] = observable_list[]
    end

    on(difflist.delete) do indices
        i = 0
        filter!(x->(i+=1; !(i in indices)), observable_list[])
        observable_list[] = observable_list[]
    end

    on(difflist.insert) do (index, item)
        insert!(observable_list[], index, item)
        observable_list[] = observable_list[]
    end
end

function Base.empty!(difflist::DiffList)
    synchronize_update(difflist) do
        difflist.empty[] = true
        empty!(children(difflist))
    end
end

Base.push!(difflist::DiffList, value) = append!(difflist, [value])

function Base.append!(difflist::DiffList, values::Vector)
    synchronize_update(difflist) do
        difflist.append[] = values
        append!(children(difflist), values)
        return values
    end
end

Base.setindex!(difflist::DiffList, value, idx::Integer) = setindex!(difflist, [value], Int[idx])

function Base.setindex!(difflist::DiffList, values::Vector, idx::Union{Integer, AbstractVector{<: Integer}, AbstractRange})
    synchronize_update(difflist) do
        indices = convert(Vector{Int}, idx)
        length(values) != length(indices) && error("Dimensions must match!")
        difflist.setindex[] = (indices, values)
        children(difflist)[indices] = values
        return values
    end
end

function Base.delete!(difflist::DiffList, idx::Union{Integer, AbstractVector{<: Integer}, AbstractRange})
    synchronize_update(difflist) do
        indices = idx isa Integer ? Int[idx] : convert(Vector{Int}, idx)
        difflist.delete[] = indices
        i = 0
        filter!(x->(i+=1; !(i in idx)), children(difflist))
        return idx
    end
end

function Base.insert!(difflist::DiffList, index::Integer, item)
    synchronize_update(difflist) do
        difflist.insert[] = (Int(index), item)
        insert!(children(difflist), index, item)
    end
end

function replace_children(difflist::DiffList, list::Vector; batch = 100)
    empty!(difflist)
    isempty(list) && return
    task = @async begin
        synchronize_update(difflist) do
            try
                for i in 1:batch:length(list)
                    append!(difflist, list[i:min(i - 1 + batch, length(list))])
                    yield()
                end
            catch e
                @warn "error in list updates" exception=CapturedException(e, Base.catch_backtrace())
            end
        end
    end
    # Wait untill task started running & locking
    while !istaskstarted(task)
        yield()
    end
    return
end

function JSServe.jsrender(session::JSServe.Session, difflist::DiffList)
    # We start with an ampty node, to not stress rendering too much!
    # We will fill the div async after creation.
    node = DOM.div(; difflist.attributes...)
    append = map(difflist.append) do values
        return JSServe.jsrender.((session,), values)
    end

    onjs(session, append, js"""function (nodes){
        var nodes_array = materialize(nodes);
        var node = $(node);
        for(var idx in nodes_array){
            node.appendChild(nodes_array[idx]);
        }
    }""")

    setindex = map(difflist.setindex) do (indices, values)
        return (indices, JSServe.jsrender.((session,), values))
    end

    onjs(session, setindex, js"""function (indices_children){
        var indices = deserialize_js(indices_children[0]); // 1 based indices
        var children = materialize(indices_children[1]);
        // children.length == indices.length should be checked in julia
        var node = $(node);
        for(var idx in children){
            var replace = node.children[indices[idx] - 1];
            node.replaceChild(children[idx], replace);
        }
    }""")

    onjs(session, difflist.delete, js"""function (indices){
        var indices = deserialize_js(indices);
        var node = $(node);
        var children2remove = indices.map(x=> node.children[x - 1]);
        for(var idx in children2remove){
            node.removeChild(children2remove[idx]);
        }
    }""")

    insert = map(difflist.insert) do (index, element)
        if element !== nothing
            return (index, JSServe.jsrender(session, element))
        else
            (index, nothing)
        end
    end

    onjs(session, insert, js"""function (index_element){
        var node = $(node);
        var element = materialize(index_element[1]);
        node.insertBefore(element, node.children[index_element[0] - 1]);
    }""")

    onjs(session, difflist.empty, js"""function (empty){
        var node = $(node);
        while (node.firstChild) {
            node.removeChild(node.firstChild);
        }
    }""")
    # Schedule to fill in nodes async & batched
    # copy == otherwise it will call empty! on children(difflist)
    # before replacing!
    replace_children(difflist, copy(children(difflist)))

    return node
end


text:
"""
Case 9:
This case study a three bus system with 1 machine (One d- One q-: 4th order model), a VSM of 19 states and an infinite source.
The perturbation increase the reference power (analogy for mechanical power) of the inverter from 1.0 to 1.2.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(dirname(@__FILE__), "data_tests/test09.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

#time span
tspan = (0.0, 20.0);
case_inv = collect(PSY.get_components(PSY.DynamicInverter, threebus_sys))[1]

#Define Fault using Callbacks
Pref_change = LITS.ControlReferenceChange(1.0, case_inv, LITS.P_ref_index, 1.2)

path = (joinpath(pwd(), "test-09"))
!isdir(path) && mkdir(path)
try
    #Define Simulation Problem
    sim = LITS.Simulation(path, threebus_sys, tspan, Pref_change)

    small_sig = small_signal_analysis(sim)

    #Solve problem in equilibrium
    run_simulation!(sim, IDA(), dtmax = 0.02)

    #Obtain data for angles
    series = get_state_series(sim, ("generator-103-1", :θ_oc))

    diff = [0.0]
    res = get_init_values_for_comparison(sim)
    for (k, v) in test09_x0_init
        diff[1] += LinearAlgebra.norm(res[k] - v)
    end
    @test (diff[1] < 1e-3)
    @test sim.solution.retcode == :Success
    @test small_sig.stable
finally
    @info("removing test files")
    rm(path, force = true, recursive = true)
end


text:
function getmedia()
  """
    function getMediaElement(e, t) {
      function E() {
          s.style.marginLeft = i.clientWidth - s.clientWidth - 2 + "px";
          if (d) {
              d.style.width = i.clientWidth / 3 + "px";
              l.style.marginLeft = i.clientWidth / 3 - 30 + "px";
              if (v) v.style["border-top-right-radius"] = "5px"
          } else {
              l.style.marginLeft = i.clientWidth - l.clientWidth - 2 + "px"
          }
          l.style.marginTop = i.clientHeight - l.clientHeight - 2 + "px";
          if (w < 10) {
              w++;
              setTimeout(E, 1e3)
          } else w = 0
      }
      t = t || {};
      if (!e.nodeName || e.nodeName.toLowerCase() != "audio" && e.nodeName.toLowerCase() != "video") {
          if (!e.getVideoTracks().length) {
              return getAudioElement(e, t)
          }
          var n = e;
          e = document.createElement(n.getVideoTracks().length ? "video" : "audio");
          e[!!navigator.mozGetUserMedia ? "mozSrcObject" : "src"] = !!navigator.mozGetUserMedia ? n : window.webkitURL.createObjectURL(n)
      }
      if (e.nodeName && e.nodeName.toLowerCase() == "audio") {
          return getAudioElement(e, t)
      }
      e.controls = false;
      var r = t.buttons || ["mute-audio", "mute-video", "full-screen", "volume-slider", "stop"];
      r.has = function(e) {
          return r.indexOf(e) !== -1
      };
      t.toggle = t.toggle || [];
      t.toggle.has = function(e) {
          return t.toggle.indexOf(e) !== -1
      };
      var i = document.createElement("div");
      i.className = "media-container";
      var s = document.createElement("div");
      s.className = "media-controls";
      i.appendChild(s);
      if (r.has("mute-audio")) {
          var o = document.createElement("div");
          o.className = "control " + (t.toggle.has("mute-audio") ? "unmute-audio selected" : "mute-audio");
          s.appendChild(o);
          o.onclick = function() {
              if (o.className.indexOf("unmute-audio") != -1) {
                  o.className = o.className.replace("unmute-audio selected", "mute-audio");
                  e.muted = false;
                  e.volume = 1;
                  if (t.onUnMuted) t.onUnMuted("audio")
              } else {
                  o.className = o.className.replace("mute-audio", "unmute-audio selected");
                  e.muted = true;
                  e.volume = 0;
                  if (t.onMuted) t.onMuted("audio")
              }
          }
      }
      if (r.has("mute-video")) {
          var u = document.createElement("div");
          u.className = "control " + (t.toggle.has("mute-video") ? "unmute-video selected" : "mute-video");
          s.appendChild(u);
          u.onclick = function() {
              if (u.className.indexOf("unmute-video") != -1) {
                  u.className = u.className.replace("unmute-video selected", "mute-video");
                  e.muted = false;
                  e.volume = 1;
                  e.play();
                  if (t.onUnMuted) t.onUnMuted("video")
              } else {
                  u.className = u.className.replace("mute-video", "unmute-video selected");
                  e.muted = true;
                  e.volume = 0;
                  e.pause();
                  if (t.onMuted) t.onMuted("video")
              }
          }
      }
      if (r.has("take-snapshot")) {
          var a = document.createElement("div");
          a.className = "control take-snapshot";
          s.appendChild(a);
          a.onclick = function() {
              if (t.onTakeSnapshot) t.onTakeSnapshot()
          }
      }
      if (r.has("stop")) {
          var f = document.createElement("div");
          f.className = "control stop";
          s.appendChild(f);
          f.onclick = function() {
              i.style.opacity = 0;
              setTimeout(function() {
                  if (i.parentNode) {
                      i.parentNode.removeChild(i)
                  }
              }, 800);
              if (t.onStopped) t.onStopped()
          }
      }
      var l = document.createElement("div");
      l.className = "volume-control";
      if (r.has("record-audio")) {
          var c = document.createElement("div");
          c.className = "control " + (t.toggle.has("record-audio") ? "stop-recording-audio selected" : "record-audio");
          l.appendChild(c);
          c.onclick = function() {
              if (c.className.indexOf("stop-recording-audio") != -1) {
                  c.className = c.className.replace("stop-recording-audio selected", "record-audio");
                  if (t.onRecordingStopped) t.onRecordingStopped("audio")
              } else {
                  c.className = c.className.replace("record-audio", "stop-recording-audio selected");
                  if (t.onRecordingStarted) t.onRecordingStarted("audio")
              }
          }
      }
      if (r.has("record-video")) {
          var h = document.createElement("div");
          h.className = "control " + (t.toggle.has("record-video") ? "stop-recording-video selected" : "record-video");
          l.appendChild(h);
          h.onclick = function() {
              if (h.className.indexOf("stop-recording-video") != -1) {
                  h.className = h.className.replace("stop-recording-video selected", "record-video");
                  if (t.onRecordingStopped) t.onRecordingStopped("video")
              } else {
                  h.className = h.className.replace("record-video", "stop-recording-video selected");
                  if (t.onRecordingStarted) t.onRecordingStarted("video")
              }
          }
      }
      if (r.has("volume-slider")) {
          var p = document.createElement("div");
          p.className = "control volume-slider";
          l.appendChild(p);
          var d = document.createElement("input");
          d.type = "range";
          d.min = 0;
          d.max = 100;
          d.value = 100;
          d.onchange = function() {
              e.volume = "." + d.value.toString().substr(0, 1)
          };
          p.appendChild(d)
      }
      if (r.has("full-screen")) {
          var v = document.createElement("div");
          v.className = "control " + (t.toggle.has("zoom-in") ? "zoom-out selected" : "zoom-in");
          if (!d && !c && !h && v) {
              s.insertBefore(v, s.firstChild)
          } else l.appendChild(v);
          v.onclick = function() {
              if (v.className.indexOf("zoom-out") != -1) {
                  v.className = v.className.replace("zoom-out selected", "zoom-in");
                  g()
              } else {
                  v.className = v.className.replace("zoom-in", "zoom-out selected");
                  m(i)
              }
          };

          function m(e) {
              if (e.requestFullscreen) {
                  e.requestFullscreen()
              } else if (e.mozRequestFullScreen) {
                  e.mozRequestFullScreen()
              } else if (e.webkitRequestFullscreen) {
                  e.webkitRequestFullscreen(Element.ALLOW_KEYBOARD_INPUT)
              }
          }

          function g() {
              if (document.fullscreen) {
                  document.cancelFullScreen()
              }
              if (document.mozFullScreen) {
                  document.mozCancelFullScreen()
              }
              if (document.webkitIsFullScreen) {
                  document.webkitCancelFullScreen()
              }
          }

          function y(e) {
              if (e.srcElement != i) return;
              var n = document.webkitIsFullScreen || document.mozFullScreen || document.fullscreen;
              i.style.width = (n ? window.innerWidth - 20 : t.width) + "px";
              i.style.display = n ? "block" : "inline-block";
              if (t.height) {
                  b.style.height = (n ? window.innerHeight - 20 : t.height) + "px"
              }
              if (!n && t.onZoomout) t.onZoomout();
              if (n && t.onZoomin) t.onZoomin();
              if (!n && v.className.indexOf("zoom-out") != -1) {
                  v.className = v.className.replace("zoom-out selected", "zoom-in");
                  if (t.onZoomout) t.onZoomout()
              }
              setTimeout(E, 1e3)
          }
          document.addEventListener("fullscreenchange", y, false);
          document.addEventListener("mozfullscreenchange", y, false);
          document.addEventListener("webkitfullscreenchange", y, false)
      }
      if (r.has("volume-slider") || r.has("full-screen") || r.has("record-audio") || r.has("record-video")) {
          i.appendChild(l)
      }
      var b = document.createElement("div");
      b.className = "media-box";
      i.appendChild(b);
      b.appendChild(e);
      if (!t.width) t.width = innerWidth / 2 - 50;
      i.style.width = t.width + "px";
      if (t.height) {
          b.style.height = t.height + "px"
      }
      b.querySelector("video").style.maxHeight = innerHeight + "px";
      var w = 0;
      if (t.showOnMouseEnter || typeof t.showOnMouseEnter === "undefined") {
          i.onmouseenter = i.onmousedown = function() {
              E();
              s.style.opacity = 1;
              l.style.opacity = 1
          };
          i.onmouseleave = function() {
              s.style.opacity = 0;
              l.style.opacity = 0
          }
      } else {
          setTimeout(function() {
              E();
              setTimeout(function() {
                  s.style.opacity = 1;
                  l.style.opacity = 1
              }, 300)
          }, 700)
      }
      E();
      i.toggle = function(e) {
          if (typeof e != "string") {
              for (var t = 0; t < e.length; t++) {
                  i.toggle(e[t])
              }
              return
          }
          if (e == "mute-audio" && o) o.onclick();
          if (e == "mute-video" && u) u.onclick();
          if (e == "record-audio" && c) c.onclick();
          if (e == "record-video" && h) h.onclick();
          if (e == "stop" && f) f.onclick();
          return this
      };
      i.media = e;
      return i
  }

  function getAudioElement(e, t) {
      function p() {
          i.style.marginLeft = r.clientWidth - i.clientWidth - 7 + "px";
          i.style.marginTop = r.clientHeight - i.clientHeight - 6 + "px";
          if (h < 10) {
              h++;
              setTimeout(p, 1e3)
          } else h = 0
      }
      t = t || {};
      if (!e.nodeName || e.nodeName.toLowerCase() != "audio" && e.nodeName.toLowerCase() != "video") {
          var n = e;
          e = document.createElement("audio");
          e[!!navigator.mozGetUserMedia ? "mozSrcObject" : "src"] = !!navigator.mozGetUserMedia ? n : window.webkitURL.createObjectURL(n)
      }
      t.toggle = t.toggle || [];
      t.toggle.has = function(e) {
          return t.toggle.indexOf(e) !== -1
      };
      e.controls = false;
      e.play();
      var r = document.createElement("div");
      r.className = "media-container";
      var i = document.createElement("div");
      i.className = "media-controls";
      r.appendChild(i);
      var s = document.createElement("div");
      s.className = "control " + (t.toggle.has("mute-audio") ? "unmute-audio selected" : "mute-audio");
      i.appendChild(s);
      s.style["border-top-left-radius"] = "5px";
      s.onclick = function() {
          if (s.className.indexOf("unmute-audio") != -1) {
              s.className = s.className.replace("unmute-audio selected", "mute-audio");
              e.muted = false;
              if (t.onUnMuted) t.onUnMuted("audio")
          } else {
              s.className = s.className.replace("mute-audio", "unmute-audio selected");
              e.muted = true;
              if (t.onMuted) t.onMuted("audio")
          }
      };
      if (!t.buttons || t.buttons && t.buttons.indexOf("record-audio") != -1) {
          var o = document.createElement("div");
          o.className = "control " + (t.toggle.has("record-audio") ? "stop-recording-audio selected" : "record-audio");
          i.appendChild(o);
          o.onclick = function() {
              if (o.className.indexOf("stop-recording-audio") != -1) {
                  o.className = o.className.replace("stop-recording-audio selected", "record-audio");
                  if (t.onRecordingStopped) t.onRecordingStopped("audio")
              } else {
                  o.className = o.className.replace("record-audio", "stop-recording-audio selected");
                  if (t.onRecordingStarted) t.onRecordingStarted("audio")
              }
          }
      }
      var u = document.createElement("div");
      u.className = "control volume-slider";
      u.style.width = "auto";
      i.appendChild(u);
      var a = document.createElement("input");
      a.style.marginTop = "11px";
      a.style.width = " 200px";
      if (t.buttons && t.buttons.indexOf("record-audio") == -1) {
          a.style.width = " 241px"
      }
      a.type = "range";
      a.min = 0;
      a.max = 100;
      a.value = 100;
      a.onchange = function() {
          e.volume = "." + a.value.toString().substr(0, 1)
      };
      u.appendChild(a);
      var f = document.createElement("div");
      f.className = "control stop";
      i.appendChild(f);
      f.onclick = function() {
          r.style.opacity = 0;
          setTimeout(function() {
              if (r.parentNode) {
                  r.parentNode.removeChild(r)
              }
          }, 800);
          if (t.onStopped) t.onStopped()
      };
      f.style["border-top-right-radius"] = "5px";
      f.style["border-bottom-right-radius"] = "5px";
      var l = document.createElement("div");
      l.className = "media-box";
      r.appendChild(l);
      var c = document.createElement("h2");
      c.innerHTML = t.title || "Audio Element";
      c.setAttribute("style", "position: absolute;color: rgb(160, 160, 160);font-size: 20px;text-shadow: 1px 1px rgb(255, 255, 255);padding:0;margin:0;");
      l.appendChild(c);
      l.appendChild(e);
      r.style.width = "329px";
      l.style.height = "90px";
      c.style.width = r.style.width;
      c.style.height = "50px";
      c.style.overflow = "hidden";
      var h = 0;
      if (t.showOnMouseEnter || typeof t.showOnMouseEnter === "undefined") {
          r.onmouseenter = r.onmousedown = function() {
              p();
              i.style.opacity = 1
          };
          r.onmouseleave = function() {
              i.style.opacity = 0
          }
      } else {
          setTimeout(function() {
              p();
              setTimeout(function() {
                  i.style.opacity = 1
              }, 300)
          }, 700)
      }
      p();
      r.toggle = function(e) {
          if (typeof e != "string") {
              for (var t = 0; t < e.length; t++) {
                  r.toggle(e[t])
              }
              return
          }
          if (e == "mute-audio" && s) s.onclick();
          if (e == "record-audio" && o) o.onclick();
          if (e == "stop" && f) f.onclick();
          return this
      };
      r.media = e;
      return r
  }
  document.write('<link rel="stylesheet" href="https://cdn.WebRTC-Experiment.com/getMediaElement.css">')
  """
end


text:
const NMAX = 100

for TYP in (:BigFloat,:Float16,:Float32,:Float64)
    A = Meta.parse("A"*string(TYP))
    @eval begin
        const $A = Array{Complex{$TYP}}(undef,NMAX,NMAX)
        riemann_array(::Type{$TYP}) = $A
    end
end

lhelm_riemanncomplex(z::Complex{T},z0::Complex{T},ζ::Complex{T},ζ0::Complex{T},E::T) where {T<:AbstractFloat} = lhelm_riemanncomplex(z,z0,ζ,ζ0,E,riemann_array(T))

function lhelm_riemanncomplex(z::Complex{T},z0::Complex{T},ζ::Complex{T},ζ0::Complex{T},E::T,A::Array{Complex{T}}) where T<:AbstractFloat
    cst = -(E/4+(z0-ζ0)/(8im))
    executerecurrence!(A,cst)

    zz0,ζζ0 = z-z0,ζ-ζ0
    ret = A[min(NMAX,2NMAX),NMAX]
    for i=min(NMAX,2NMAX)-1:-1:div(NMAX+1,2)
        ret = muladd(zz0,ret,A[i,NMAX])
    end
    ret = ret*zz0^div(NMAX+3,2)

    for j=NMAX-1:-1:1
        temp = A[min(NMAX-1,2j),j]
        for i=min(NMAX-1,2j)-1:-1:div(j+1,2)
            temp = muladd(zz0,temp,A[i,j])
        end
        temp = temp*zz0^div(j+3,2)
        ret = muladd(ret,ζζ0,temp)
    end
    1+ret
end

lhelm_riemann(x::AbstractFloat,x0::AbstractFloat,y::AbstractFloat,y0::AbstractFloat,E::AbstractFloat) = lhelm_riemanncomplex(complex(x,y),complex(x0,y0),complex(x,-y),complex(x0,-y0),E)
lhelm_riemann(x::Union{T1,Complex{T1}},y::Union{T2,Complex{T2}},E::AbstractFloat) where {T1<:AbstractFloat,T2<:AbstractFloat} = lhelm_riemann(real(x),real(y),imag(x),imag(y),E)

lhelm_riemann(x::Vector,y::Vector,E::AbstractFloat) = promote_type(eltype(x),eltype(y),typeof(E))[lhelm_riemann(x[k],y[k],E) for k in 1:length(x)]
lhelm_riemann(x::Matrix,y::Matrix,E::AbstractFloat) = reshape(lhelm_riemann(vec(x),vec(y),E),size(x))

function executerecurrence!(A::Array{Complex{T}},cst::Complex{T}) where T<:AbstractFloat
    eighthim = one(T)/(8im)
    sixteenthim = eighthim/2

    A[1,1] = cst
    A[2,1] = -sixteenthim

    A[1,2] = sixteenthim
    A[2,2] = cst*A[1,1]/4
    A[3,2] = (cst*A[2,1] - eighthim*A[1,1])/6
    A[4,2] = -eighthim*A[2,1]/8

    for j=3:NMAX
        if isodd(j)
            i = div(j+1,2)
            A[i,j] = (cst*A[i-1,j-1] + eighthim*A[i-1,j-2])/(i*j)
        else
            i = div(j+1,2)
            A[i,j] = (eighthim*A[i-1,j-2])/(i*j)
            i+=1
            A[i,j] = (cst*A[i-1,j-1] + eighthim*A[i-1,j-2])/(i*j)
        end
        for i = div(j+4,2):min(2j-3,NMAX)
            A[i,j] = (cst*A[i-1,j-1] + eighthim*A[i-1,j-2] - eighthim*A[i-2,j-1])/(i*j)
        end
        if 2j-2 ≤ NMAX
            i = 2j-2
            A[i,j] = (cst*A[i-1,j-1] - eighthim*A[i-2,j-1])/(i*j)
            if 2j-1 ≤ NMAX
                i = 2j-1
                A[i,j] = (cst*A[i-1,j-1] - eighthim*A[i-2,j-1])/(i*j)
                if 2j ≤ NMAX
                    i = 2j
                    A[i,j] = (-eighthim*A[i-2,j-1])/(i*j)
                end
            end
        end
    end
end


text:
using Logging, MiniLoggers
global_logger(MiniLogger(level=Logging.Debug))

module Configuration

using ODBC

const conn1 = ODBC.Connection("Driver=/usr/lib/x86_64-linux-gnu/odbc/psqlodbcw.so;Server=127.0.0.1;Port=<db port>;Database=<db name>;Uid=<user name>;Pwd=<password>;")

precompile(ODBC.Connection, (String, ))

end # module

const conn = ZulipReminderBot.Adapter(Configuration.conn1, ZulipReminderBot.PostgresODBC)


text:
function rad(mtime)

  global SW,LW,albedo

  ########################################################################
  # radiation
  ########################################################################
  AOI=angle_of_incidence(lat,lon,jd,t);
  AOI<0 ? AOI=0 : AOI

  #define albedo
  if theta_l[end]<0.1
    albedo = 0.25;
  else
    if theta_l[end]>0.25
      albedo = 0.1;
    else
      albedo = 0.35-theta_l[end];
    end
  end

  SW=albedo.*(SunConstant.*AOI);
  # longwave radiation at the surface from Stensrud eq. 8.24-8.25
  ϵa=0.725+0.17*log10(sum(qA/1000*1.2*ΔZA)*100)
  LW=emissivity*sb*((TS[end])^4 - ϵa*(TAup)^4);

end


text:
#####
## Layer and Network Architectures for TintiNet.jl
#

# 1. 1D Inception V2 layer with 3 window sizes (Inception4)

"""
Structure for the 1D Inception V2 layer with 3 window sizes (singleton, trigram, pentagram) and a pooling layer
"""
struct Inception4{L1, L2, L3, L4}
    Solo_One_Layer::L1
    One_Three_Layer::L2
    One_Five_Layer::L3
    Pool_One_Layer::L4
end

@functor Inception4

function Inception4(input_channels)
    output_channels = fld(input_channels, 4)
    return Inception4(
        Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 3), output_channels => output_channels, relu; pad = (0,1))
        ),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 5), output_channels => output_channels, relu; pad = (0,2))
        ),
        Chain(
            MaxPool((1,3); pad=(0,1), stride=1),
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0))
        ),
    )
end

function (m::Inception4)(x)

    res1 = m.Solo_One_Layer(x)
    res2 = m.One_Three_Layer(x)
    res3 = m.One_Five_Layer(x)
    res4 = m.Pool_One_Layer(x)

    concatenated_result = cat(res1, res2, res3, res4; dims = (3))

end

# 2. 1D Inception V2 layer with 7 window sizes (Inception8 or Inceptigor)

"""
Structure for the 1D Inception V2 layer with 7 window sizes and a pooling layer
"""
struct Inception8{L1, L2, L3, L4, L5, L6, L7, L8}
    Solo_One_Layer::L1
    One_Three_Layer::L2
    One_Five_Layer::L3
    One_Seven_Layer::L4
    One_Nine_Layer::L5
    One_Eleven_Layer::L6
    One_Thirteen_Layer::L7
    Pool_One_Layer::L8
end

@functor Inception8

function Inception8(input_channels)
    output_channels = fld(input_channels, 8)
    return Inception8(
        Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 3), output_channels => output_channels, relu; pad = (0,1))
        ),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 5), output_channels => output_channels, relu; pad = (0,2))
        ),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 7), output_channels => output_channels, relu; pad = (0,3))
        ),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 9), output_channels => output_channels, relu; pad = (0,4))
        ),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 11), output_channels => output_channels, relu; pad = (0,5))
        ),
        Chain(
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0)),
            Conv((1, 13), output_channels => output_channels, relu; pad = (0,6))
        ),
        Chain(
            MaxPool((1,3); pad=(0,1), stride=1),
            Conv((1, 1), input_channels => output_channels, relu; pad = (0,0))
        ),
    )
end

function (m::Inception8)(x)

    res1 = m.Solo_One_Layer(x)
    res2 = m.One_Three_Layer(x)
    res3 = m.One_Five_Layer(x)
    res4 = m.One_Seven_Layer(x)
    res5 = m.One_Nine_Layer(x)
    res6 = m.One_Eleven_Layer(x)
    res7 = m.One_Thirteen_Layer(x)
    res8 = m.Pool_One_Layer(x)

    concatenated_result = cat(res1, res2, res3, res4, res5, res6, res7, res8; dims = (3))

end


text:
# Dantzig-Wolfe Reformulation and Column Generation
# Data for General Assignment Problem 1
# Edward J. Xu
# 2019.5.1

function getData()
    vec_c = [6 4 6 1 3 4 1 2 8]
    mat_a =[
    1 0 0 1 0 0 1 0 0 ;
    0 1 0 0 1 0 0 1 0 ;
    0 0 1 0 0 1 0 0 1 ;
    7 2 8 0 0 0 0 0 0 ;
    0 0 0 8 7 6 0 0 0 ;
    0 0 0 0 0 0 9 1 9]
    vec_b = hcat([1;1;1;9;7;10])
    index_sub = [[1,2,3], [4,5,6], [7,8,9]]
    num_sub = 3
    return (vec_c, mat_a, vec_b, index_sub, num_sub)
end

(vec_c, mat_a, vec_b, index_sub, num_sub) = getData()


text:
module epidemic
export Cell, Parameters, create_map, to_colors, interact!, update, count_infections, count_deaths

"Enumerate possible states of a single cell"
@enum InfectionStatus uninfected infected dead recovered immune

"Data structure containing the infection status of a cell"
mutable struct Cell
    status::InfectionStatus
    infection_time::Int8
end

"Parameters for a simulation"
mutable struct Parameters
    infection_rate::Float32
    death_probability::Float32
    recovery_time::Int8
    immunity_rate::Float32
end

"Create a map of immune and uninfected cells with 1 infected in the middle."
function make_cells(width::Integer=7, height::Integer=7, immunity_rate::Float32)
    cells = Matrix{Cell}(undef, width, height)
    for i in 1:size(cells)[1]
        for j in 1:size(cells)[2]
            if rand(1)[1] < immunity_rate
                cells[i,j] = Cell(immune, 0)
            else
                cells[i,j] = Cell(uninfected, 0)
            end
        end
    end
    cells[width÷2+1,height÷2+1] = 1
end

"Map the cells to colors for plotting"
function to_colors(cell::Cell)
    if cell.status == uninfected
        return RGB(0.0,0.8,0.0)
    end
    if cell.status == infected
        return RGB(0.8,0.0,0.0)
    end
    if cell.status == dead
        return RGB(0.1,0.1,0.1)
    end
    if cell.status == recovered
        return RGB(0.0,0.0,0.8)
    end
    if cell.status == immune
        return RGB(0.2,0.2,1.0)
    end
end

"""
Simulate the interaction between one cell and a neighbour.

If the neighbour is infected, it infect this cell with the propability parameters.infection_rate.
"""
function interact!(new_cell::Cell, other_cell::Cell, parameters::Parameters)
    if new_cell.status == uninfected && other_cell.status == infected
        if rand(1)[1] < parameters.infection_rate
            new_cell.status = infected
            new_cell.infection_time = 0
        end
    end
end

"""
Update a single cell, not accounting for it's interactions with the neighbours.
"""
function update!(new_cell::Cell, parameters::Parameters)
    if new_cell.status == infected
        # Increase the number of time steps since infection
        new_cell.infection_time += 1

        # Check if the cell recovers
        if new_cell.infection_time > parameters.recovery_time
            new_cell.status = recovered
        end

        # Check if it dies
        if rand(1)[1] < parameters.death_probability
            new_cell.status = dead
        end
    end
end

"""
Update the cells in the 2D array of Cells, using given parameters.
"""
function update!(cells::Matrix{Cell}, parameters::Parameters)
    old_cells = deepcopy(cells)
    for i in 1:size(cells)[1]
        for j in 1:size(cells)[2]
            update!(cells[i,j], parameters)
        end
    end
    for i in 1:size(cells)[1]-1
        for j in 1:size(cells)[2]
            interact!(cells[i,j], old_cells[i+1,j], parameters)
            interact!(cells[i+1,j], old_cells[i,j], parameters)
        end
    end
    for i in 1:size(cells)[1]
        for j in 1:size(cells)[2]-1
            interact!(cells[i,j], old_cells[i,j+1], parameters)
            interact!(cells[i,j+1], old_cells[i,j], parameters)
        end
    end
end

"Count the current number of infections"
function count_infections(cells::Matrix{Cell})
    infections = 0
    for i in 1:size(cells)[1]
        for j in 1:size(cells)[2]
            if cells[i,j].status == infected
                infections += 1
            end
        end
    end
    return infections
end

"Count the number of dead cells"
function count_deaths(cells::Matrix{Cell})
    deaths = 0
    for i in 1:size(cells)[1]
        for j in 1:size(cells)[2]
            if cells[i,j].status == dead
                deaths += 1
            end
        end
    end
    return deaths
end

# End module
end

max_stars_repo_name: tpfau/julia-introduction

max_issues_repo_name: tpfau/julia-introduction

max_forks_repo_name: tpfau/julia-introduction