text:
module PWscf

using Compat: only
# using Dates: DateTime, DateFormat
using DataFrames: AbstractDataFrame, DataFrame, groupby
using Parameters: @with_kw
using QuantumESPRESSOBase.Inputs.PWscf
using VersionParsing: vparse

using ..Outputs: SubroutineError

export Diagonalization,
    Preamble,
    Davidson,
    ConjugateGradient,
    ProjectedPreconditionedConjugateGradient,
    parse_fft_base_info,
    parse_symmetries,
    parse_ibz,
    parse_stress,
    parse_iteration_time,
    parse_bands,
    parse_all_electron_energy,
    parse_energy_decomposition,
    parse_paw_contribution,
    parse_smearing_energy,
    parse_version,
    parse_parallel_info,
    parse_fft_dimensions,
    parse_iteration_head,
    parse_electrons_energies,
    parse_clock,
    parse_input_name,
    isoptimized,
    isjobdone,
    tryparsefirst,
    parsefirst,
    tryparseall,
    parseall,
    tryparselast,
    parselast,
    tryparsenext,
    parsenext,
    tryparsefinal,
    parsefinal

include("regexes.jl")

# From https://discourse.julialang.org/t/aliases-for-union-t-nothing-and-union-t-missing/15402/4
const Maybe{T} = Union{T,Nothing}  # Should not be exported

abstract type Diagonalization end
struct Davidson <: Diagonalization end
struct ConjugateGradient <: Diagonalization end
struct ProjectedPreconditionedConjugateGradient <: Diagonalization end

@with_kw struct Preamble
    ibrav::Int
    alat::Float64
    omega::Float64
    nat::Int
    ntyp::Int
    nelec::Float64
    nelup::Maybe{Float64} = nothing
    neldw::Maybe{Float64} = nothing
    nbnd::Int
    ecutwfc::Float64
    ecutrho::Float64
    ecutfock::Maybe{Float64} = nothing
    conv_thr::Maybe{Float64} = nothing
    mixing_beta::Maybe{Float64} = nothing
    mixing_ndim::Maybe{Int} = nothing
    mixing_mode::Maybe{String} = nothing
    xc::String
    nstep::Maybe{Int} = nothing
end

"""
    parse_fft_base_info(str::AbstractString)

Parse the FFT base information from `pw.x`'s output and return a `DataFrame`.

If there are more than one processors, the title is "Parallelization info" and three
rows, i.e., "Min", "Max", and "Sum" are printed. If not, the title is
"G-vector sticks info" and only the "Sum" row is printed. If no information is found,
return `nothing`. The `DataFrame` is grouped by "sticks" and "gvecs".
"""
function parse_fft_base_info(str::AbstractString)::Maybe{AbstractDataFrame}
    df =
        DataFrame(kind = String[], stats = String[], dense = Int[], smooth = Int[], PW = [])
    m = match(FFT_BASE_INFO, str)
    if m === nothing
        @info("The FFT base info is not found!")
        return
    end
    body = m[:body]
    for line in split(body, r"\R+")  # Don’t want empty lines
        # "Min",4X,2I8,I7,12X,2I9,I8
        sp = split(line, " ", keepempty = false)  # Don't want empty strings
        numbers = map(x -> parse(Int, x), sp[2:7])
        push!(df, ["sticks" sp[1] numbers[1:3]...])
        push!(df, ["gvecs" sp[1] numbers[4:6]...])
    end
    return df
end # function parse_fft_base_info

function parse_symmetries(str::AbstractString)
    m = match(SYM_OPS, str)
    m === nothing && return
    num_sym_ops = isempty(m[:n]) ? 0 : parse(Int, m[:n])
end # function parse_symmetries

# Return `nothing`, `(cartesian_coordinates, nothing)`, `(nothing, crystal_coordinates)`, `(cartesian_coordinates, crystal_coordinates)`
function parse_ibz(str::AbstractString)::Maybe{Tuple}
    m = match(K_POINTS_BLOCK, str)
    if m === nothing
        @info("The k-points info is not found!")
        return
    end
    nk = parse(Int, m[:nk])
    result = []
    kinds = (:cart => "tpiba", :cryst => "crystal")
    for (k, v) in kinds
        if m[k] !== nothing
            x = Matrix{Float64}(undef, nk, 4)
            for (i, m) in enumerate(eachmatch(K_POINTS_ITEM, m[k]))
                x[i, :] = map(x -> parse(Float64, x), m.captures[1:end])
            end
            push!(result, SpecialPointsCard(x, v))
        else
            push!(result, nothing)
        end
    end
    return Tuple(result)
end # function parse_ibz

function parse_stress(str::AbstractString)
    pressures = Float64[]
    atomic_stresses, kbar_stresses = Matrix{Float64}[], Matrix{Float64}[]
    for m in eachmatch(STRESS_BLOCK, str)
        pressure, content = m.captures[1], m.captures[3]
        push!(pressures, parse(Float64, pressure))

        stress_atomic, stress_kbar = ntuple(_ -> Matrix{Float64}(undef, 3, 3), 2)
        for (i, line) in enumerate(split(content, r"\R+"))
            tmp = map(x -> parse(Float64, x), split(line, " ", keepempty = false))
            stress_atomic[i, :], stress_kbar[i, :] = tmp[1:3], tmp[4:6]
        end
        push!(atomic_stresses, stress_atomic)
        push!(kbar_stresses, stress_kbar)
    end
    return pressures, atomic_stresses, kbar_stresses
end # function parse_stress

function _iterationwise!(f::Function, df::AbstractDataFrame, str::AbstractString)
    # Loop relax steps
    for (i, scf) in enumerate(eachmatch(SELF_CONSISTENT_CALCULATION_BLOCK, str))
        # Loop scf iterations
        for (j, iter) in enumerate(eachmatch(ITERATION_BLOCK, scf[1]))
            push!(df, [i j f(iter[1])...])
        end
    end
    return df
end # function _iterationwise

function parse_iteration_head(str::AbstractString)
    df = DataFrame(
        step = Int[],  # Step number
        iteration = Int[],  # Iteration number
        ecut = Float64[],  # Cutoff energy
        β = Float64[],  # Mixing beta
    )
    return _iterationwise!(_parse_iteration_head, df, str)
end # function parse_iteration_head
# This is a helper function and should not be exported.
function _parse_iteration_head(str::AbstractString)
    head = match(ITERATION_HEAD, str)
    return map(x -> parse(Float64, x), head.captures[2:3])
end # function _parse_iteration_head

function parse_iteration_time(str::AbstractString)
    df = DataFrame(step = Int[], iteration = Int[], time = Float64[])
    return _iterationwise!(_parse_iteration_time, df, str)
end # function parse_iteration_time
# This is a helper function and should not be exported.
function _parse_iteration_time(str::AbstractString)
    return parse(Float64, match(TOTAL_CPU_TIME, str)[1])
end # function _parse_iteration_time

function parse_diagonalization(str::AbstractString)
    df = DataFrame(
        step = Int[],
        iteration = Int[],
        diag = Diagonalization[],  # Diagonalization style
        ethr = Float64[],  # Energy threshold
        avg = Float64[],  # Average # of iterations
    )
    return _iterationwise!(_parse_diagonalization, df, str)
end # function parse_diagonalization
# This is a helper function and should not be exported.
function _parse_diagonalization(str::AbstractString)
    solver, ethr, avg_iter = nothing, nothing, nothing  # Initialization
    m = match(C_BANDS, str)
    if m !== nothing
        solver = if m[:diag] == "Davidson diagonalization with overlap"
            Davidson()
        elseif m[:diag] == "CG style diagonalization"
            ConjugateGradient()
        elseif m[:diag] == "PPCG style diagonalization"
            ProjectedPreconditionedConjugateGradient()
        else
            error("unknown diagonalization style!")
        end
        ethr, avg_iter = map(x -> parse(Float64, x), m.captures[2:end])
    end  # Keep them `nothing` if `m` is `nothing`
    return solver, ethr, avg_iter
end # function _parse_diagonalization

function _parse_nonconverged_energy(str::AbstractString)
    ɛ, hf, δ = nothing, nothing, nothing  # Initialization
    m = match(UNCONVERGED_ELECTRONS_ENERGY, str)
    if m !== nothing
        ɛ, hf, δ = map(x -> x === nothing ? x : parse(Float64, x), m.captures)
    end  # Keep them `nothing` if `m` is `nothing`
    return ɛ, hf, δ
end # function _parse_nonconverged_energy
function _parse_electrons_energies(str::AbstractString, ::Val{:nonconverged})
    df = DataFrame(
        step = Int[],
        iteration = Int[],
        ɛ = Maybe{Float64}[],  # Total energy
        hf = Maybe{Float64}[],  # Harris-Foulkes estimate
        δ = Maybe{Float64}[],  # Estimated scf accuracy
    )
    return _iterationwise!(_parse_nonconverged_energy, df, str)
end # function _parse_electrons_energies
function _parse_electrons_energies(str::AbstractString, ::Val{:converged})
    df = DataFrame(
        step = Int[],
        ɛ = Maybe{Float64}[],  # Total energy
        hf = Maybe{Float64}[],  # Harris-Foulkes estimate
        δ = Maybe{Float64}[],  # Estimated scf accuracy
    )
    for (i, m) in enumerate(eachmatch(CONVERGED_ELECTRONS_ENERGY, str))
        data = if m !== nothing
            map(x -> x === nothing ? x : parse(Float64, x), m.captures[1:3])
        else
            ntuple(_ -> nothing, 3)
        end  # Keep them `nothing` if `m` is `nothing`
        push!(df, [i data...])
    end
    return df
end # function _parse_electrons_energies
function _parse_electrons_energies(str::AbstractString, ::Val{:combined})
    converged = parse_electrons_energies(str, :converged)
    nonconverged = parse_electrons_energies(str, :nonconverged)
    # TODO: Very ugly hack
    m = 1  # Initial step number
    for (i, n) in enumerate(nonconverged.step)
        if n != m
            @assert(all(==(nothing), nonconverged[i-1, 3:5]))
            # nonconverged[i - 1, 3:5] = converged[n, 2:4]  # Converged energies do not have `iteration` column
            nonconverged[i-1, 3] = converged[n, 2]
            nonconverged[i-1, 4] = converged[n, 3]
            nonconverged[i-1, 5] = converged[n, 4]
        end
        m = n  # Save the last step number
    end
    if m == 1
        nonconverged[end, 3] = converged[end, 2]
        nonconverged[end, 4] = converged[end, 3]
        nonconverged[end, 5] = converged[end, 4]
    end
    return nonconverged
end # function _parse_electrons_energies
function parse_electrons_energies(str::AbstractString, option::Symbol)
    @assert(option ∈ (:combined, :converged, :nonconverged))
    return _parse_electrons_energies(str, Val(option))
end # function parse_electrons_energies

# See https://github.com/QEF/q-e/blob/4132a64/PW/src/print_ks_energies.f90#L10.
function parse_bands(str::AbstractString)
    str == "Number of k-points >= 100: set verbosity='high' to print the bands." && return
    kpts, bands = nothing, nothing  # Initialization
    m = match(KS_ENERGIES_BLOCK, str)
    if m !== nothing
        kpts, bands = Vector{Float64}[], Vector{Float64}[]
        regex =
            match(KS_ENERGIES_BANDS, str) === nothing ? KS_ENERGIES_BAND_ENERGIES :
            KS_ENERGIES_BANDS
        for m in eachmatch(regex, str)
            push!(
                kpts,
                map(x -> parse(Float64, x[1]), eachmatch(Regex(GENERAL_REAL), m[:k])),
            )
            push!(
                bands,
                map(x -> parse(Float64, x[1]), eachmatch(Regex(GENERAL_REAL), m[:band])),
            )
        end
        len, nbnd = length(kpts), length(bands[1])
        kpts, bands = reshape(Iterators.flatten(kpts) |> collect, len, 3),
        reshape(Iterators.flatten(bands) |> collect, len, nbnd)
    end  # Keep them `nothing` if `m` is `nothing`
    return kpts, bands
end # function parse_bands

function parse_all_electron_energy(str::AbstractString)
    df = DataFrame(step = Int[], ae = Maybe{Float64}[])
    for (i, m) in enumerate(eachmatch(CONVERGED_ELECTRONS_ENERGY, str))
        ae = if any(==(nothing), (m, m[:ae]))
            nothing
        else
            parse(Float64, match(Regex(FIXED_POINT_REAL), m[:ae])[1])
        end
        push!(df, [i ae])
    end
    return df
end # function parse_all_electron_energy

function parse_energy_decomposition(str::AbstractString)
    df = DataFrame(
        step = Int[],
        onelectron = Maybe{Float64}[],
        hartree = Maybe{Float64}[],
        xc = Maybe{Float64}[],
        ewald = Maybe{Float64}[],
    )
    for (i, m) in enumerate(eachmatch(CONVERGED_ELECTRONS_ENERGY, str))
        data = if any(==(nothing), (m, m[:decomp]))
            ntuple(_ -> nothing, 4)
        else
            map(x -> parse(Float64, x[1]), eachmatch(Regex(FIXED_POINT_REAL), m[:decomp]))
        end
        push!(df, [i data...])
    end
    return df
end # function parse_energy_decomposition

function parse_paw_contribution(str::AbstractString)
    df = DataFrame(
        step = Int[],
        one_electron = Maybe{Float64}[],
        hartree_ae = Maybe{Float64}[],
        hartree_ps = Maybe{Float64}[],
        xc_ae = Maybe{Float64}[],
        xc_ps = Maybe{Float64}[],
        eh = Maybe{Float64}[],
        exc = Maybe{Float64}[],
    )
    for (i, m) in enumerate(eachmatch(CONVERGED_ELECTRONS_ENERGY, str))
        data = if any(==(nothing), (m, m[:one]))
            ntuple(_ -> nothing, 6)
        else
            map(x -> parse(Float64, x[1]), eachmatch(Regex(FIXED_POINT_REAL), m[:one]))
        end
        push!(df, [i data...])
    end
    return df
end # function parse_paw_contribution

function parse_smearing_energy(str::AbstractString)
    df = DataFrame(step = Int[], smearing = Maybe{Float64}[])
    for (i, m) in enumerate(eachmatch(CONVERGED_ELECTRONS_ENERGY, str))
        smearing = if any(==(nothing), (m, m[:smearing]))
            nothing
        else
            parse(Float64, match(Regex(FIXED_POINT_REAL), m[:smearing])[1])
        end
        push!(df, [i smearing])
    end
    return df
end # function parse_smearing_energy

function parse_version(str::AbstractString)::Maybe{VersionNumber}
    m = match(PWSCF_VERSION, str)
    m !== nothing ? vparse(m[:version]) : return
end # function parse_version

function parse_parallel_info(str::AbstractString)::Maybe{Tuple{String,Int}}
    m = match(PARALLEL_INFO, str)
    m === nothing && return
    return m[:kind], m[:num] === nothing ? 1 : parse(Int, m[:num])
end # function parse_parallel_info

function parse_fft_dimensions(str::AbstractString)::Maybe{NamedTuple}
    m = match(FFT_DIMENSIONS, str)
    m === nothing && return
    parsed = map(x -> parse(Int, x), m.captures)
    return (; zip((:ng, :nr1, :nr2, :nr3), parsed)...)
end # function parse_fft_dimensions

function parse_clock(str::AbstractString)::Maybe{AbstractDataFrame}
    m = match(TIME_BLOCK, str)
    m === nothing && return
    content = only(m.captures)

    info = DataFrame(
        subroutine = String[],
        item = String[],
        CPU = Float64[],
        wall = Float64[],
        calls = Int[],
    )
    for regex in [
        SUMMARY_TIME_BLOCK
        INIT_RUN_TIME_BLOCK
        ELECTRONS_TIME_BLOCK
        C_BANDS_TIME_BLOCK
        SUM_BAND_TIME_BLOCK
        EGTERG_TIME_BLOCK
        H_PSI_TIME_BLOCK
        GENERAL_ROUTINES_TIME_BLOCK
        PARALLEL_ROUTINES_TIME_BLOCK
    ]
        block = match(regex, content)
        if block !== nothing
            for m in eachmatch(TIME_ITEM, block[:body])
                push!(
                    info,
                    [block[:head] m[1] map(x -> parse(Float64, x), m.captures[2:4])...],
                )
            end
        end
    end
    # m = match(TERMINATED_DATE, content)
    # info["terminated date"] = parse(DateTime, m.captures[1], DateFormat("H:M:S"))
    return info
end # function parse_clock

function parse_input_name(str::AbstractString)
    m = match(READING_INPUT_FROM, str)
    return m === nothing ? nothing : only(m)
end # function parse_input_name

isoptimized(str::AbstractString) =
    match(FINAL_COORDINATES_BLOCK, str) === nothing ? false : true

isjobdone(str::AbstractString) = match(JOB_DONE, str) !== nothing

# This is an internal function and should not be exported.
function Base.tryparse(::Type{Preamble}, str::AbstractString)
    dict = Dict{Symbol,Any}()
    m = match(SUMMARY_BLOCK, str)
    return if m !== nothing
        body = only(m.captures)
        f = (T, x) -> T == String ? string(x) : parse(T, x)
        for (field, regex, T) in zip(
            (
                :ibrav,
                :alat,
                :omega,
                :nat,
                :ntyp,
                :nelec,
                :nbnd,
                :ecutwfc,
                :ecutrho,
                :ecutfock,
                :conv_thr,
                :mixing_beta,
                :mixing_ndim,
                :xc,
                :nstep,
            ),
            (
                NUMBER_OF_ATOMS_PER_CELL,
                LATTICE_PARAMETER,
                UNIT_CELL_VOLUME,
                NUMBER_OF_ATOMS_PER_CELL,
                NUMBER_OF_ATOMIC_TYPES,
                NUMBER_OF_ELECTRONS,
                NUMBER_OF_KOHN_SHAM_STATES,
                KINETIC_ENERGY_CUTOFF,
                CHARGE_DENSITY_CUTOFF,
                CUTOFF_FOR_FOCK_OPERATOR,
                CONVERGENCE_THRESHOLD,
                MIXING_BETA,
                NUMBER_OF_ITERATIONS_USED,
                EXCHANGE_CORRELATION,
                NSTEP,
            ),
            (
                Int,
                Float64,
                Float64,
                Int,
                Int,
                Float64,
                Int,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Int,
                String,
                Int,
            ),
        )
            m = match(regex, body)
            if m !== nothing
                dict[field] = f(T, m[1])
            end
        end
        # 2 special cases
        let x = match(NUMBER_OF_ELECTRONS, body), y = match(NUMBER_OF_ITERATIONS_USED, body)
            if all(!=(nothing), x.captures[2:end])
                dict[:nelup], dict[:neldw] = parse.(Float64, x.captures[2:end])
            end
            if y !== nothing
                dict[:mixing_mode] = y[2]
            end
        end
        Preamble(; dict...)
    end
end # function Base.tryparse
function Base.tryparse(::Type{SubroutineError}, str::AbstractString)
    # According to my observation, a QE output can have at most one type of
    # `SubroutineError`. Warn me if there can be multiple types of errors.
    m = match(ERROR_BLOCK, str)
    return if m !== nothing
        # `tryparse` returns nothing if the string does not contain what we want,
        # while `parse` raises an error.
        body = strip(m[:body])
        # Referenced from https://stackoverflow.com/a/454919/3260253
        e, msg = map(strip, split(body, r"\R+"))
        m = match(ERROR_IN_ROUTINE, e)
        SubroutineError(m[1], m[2], msg)
    end
end # function Base.tryparse

function tryparse_internal(::Type{CellParametersCard}, str::AbstractString)
    m = match(CELL_PARAMETERS_BLOCK, str)
    return if m !== nothing
        body, data = m[:data], Matrix{Float64}(undef, 3, 3)  # Initialization
        for (i, matched) in enumerate(eachmatch(CELL_PARAMETERS_ITEM, body))
            data[i, :] = map(x -> parse(Float64, x), matched.captures)
        end
        if m[:option] == "alat"
            alat = parse(Float64, m[:alat])
            CellParametersCard(alat * data, "bohr")
        else
            CellParametersCard(data, m[:option])
        end
    end
end # function tryparse_internal
function tryparse_internal(::Type{AtomicPositionsCard}, str::AbstractString)
    atomic_positions = AtomicPositionsCard[]
    m = match(ATOMIC_POSITIONS_BLOCK, str)
    return if m !== nothing
        option = string(m[1])
        body = m[2]
        data = AtomicPosition[]
        for matched in eachmatch(ATOMIC_POSITIONS_ITEM, body)
            captured = matched.captures
            if_pos = map(x -> x === nothing ? 1 : parse(Int, x), captured[5:7])
            atom, pos = captured[1], map(x -> parse(Float64, x), captured[2:4])
            push!(data, AtomicPosition(atom, pos, if_pos))
        end
        AtomicPositionsCard(data, option)
    end
end # function tryparse_internal

const AtomicStructure = Union{CellParametersCard,AtomicPositionsCard}

function Base.parse(
    ::Type{T},
    str::AbstractString,
) where {T<:Union{Preamble,SubroutineError}}
    x = tryparse(T, str)
    x === nothing ? throw(Meta.ParseError("cannot find `$(T)`!")) : x
end # function Base.parse

function _parse(::Type{T}, str::AbstractString) where {T<:AtomicStructure}
    x = tryparse(T, str)
    x === nothing ? throw(Meta.ParseError("cannot find `$(T)`!")) : x
end # function _parse

const REGEXOF = (
    CellParametersCard = CELL_PARAMETERS_BLOCK,
    AtomicPositionsCard = ATOMIC_POSITIONS_BLOCK,
)

tryparsefirst(::Type{T}, str::AbstractString) where {T<:AtomicStructure} =
    tryparse_internal(T, str)
parsefirst(::Type{T}, str::AbstractString) where {T<:AtomicStructure} = _parse(T, str)

function tryparseall(::Type{T}, str::AbstractString) where {T<:AtomicStructure}
    return map(eachmatch(REGEXOF[nameof(T)], str)) do x
        try
            tryparse_internal(T, x.match)
        catch
            nothing
        end
    end
end # function parseall
function parseall(::Type{T}, str::AbstractString) where {T<:AtomicStructure}
    return map(eachmatch(REGEXOF(T), str)) do x
        try
            tryparse_internal(T, x.match)
        catch
            Meta.ParseError("Pass failed!")
        end
    end
end # function parseall

tryparselast(::Type{T}, str::AbstractString) where {T<:AtomicStructure} =
    tryparseall(T, str)[end]
parselast(::Type{T}, str::AbstractString) where {T<:AtomicStructure} = parseall(T, str)[end]

function _parsenext_internal(
    ::Type{T},
    str::AbstractString,
    start::Integer,
    raise::Bool,
) where {T}
    x = findnext(REGEXOF(T), str, start)
    if x === nothing
        raise ? throw(Meta.ParseError("Nothing found for next!")) : return
    end
    return tryparse_internal(T, str[x])
end # function parsenext
tryparsenext(::Type{T}, str::AbstractString, start::Integer) where {T} =
    _parsenext_internal(T, str, start, false)
parsenext(::Type{T}, str::AbstractString, start::Integer) where {T} =
    _parsenext_internal(T, str, start, true)

function tryparsefinal(::Type{T}, str::AbstractString) where {T<:AtomicStructure}
    m = match(FINAL_COORDINATES_BLOCK, str)
    m === nothing && return
    m = match(REGEXOF[nameof(T)], m.match)
    m === nothing && return
    return tryparse_internal(T, m.match)
end # function parsefinal
function parsefinal(::Type{T}, str::AbstractString) where {T<:AtomicStructure}
    m = match(FINAL_COORDINATES_BLOCK, str)
    m === nothing && throw(Meta.ParseError("No final coordinates found!"))
    m = match(REGEXOF[nameof(T)], m.match)
    m === nothing && throw(Meta.ParseError("No `CELL_PARAMETERS` found!"))
    return tryparse_internal(T, m.match)
end # function parsefinal

end

ext: jl


text:
using CometLogger
using Documenter

DocMeta.setdocmeta!(CometLogger, :DocTestSetup, :(using CometLogger); recursive=true)

makedocs(;
    modules=[CometLogger],
    authors="rejuvyesh <mail@rejuvyesh.com> and contributors",
    repo="https://github.com/rejuvyesh/CometLogger.jl/blob/{commit}{path}#{line}",
    sitename="CometLogger.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rejuvyesh.github.io/CometLogger.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/rejuvyesh/CometLogger.jl",
)


text:
@with_kw mutable struct OnPolicySolver <: Solver
    agent::PolicyParams # Policy
    S::AbstractSpace # State space
    N::Int = 1000 # Number of environment interactions
    ΔN::Int = 200 # Number of interactions between updates
    max_steps::Int = 100 # Maximum number of steps per episode
    log::Union{Nothing, LoggerParams} = nothing # The logging parameters
    i::Int = 0 # The current number of environment interactions
    param_optimizers::Dict{Any, TrainingParams} = Dict() # Training parameters for the parameters
    a_opt::TrainingParams # Training parameters for the actor
    c_opt::Union{Nothing, TrainingParams} = nothing # Training parameters for the critic
    𝒫::NamedTuple = (;) # Parameters of the algorithm

    # On-policy-specific parameters
    λ_gae::Float32 = 0.95 # Generalized advantage estimation parameter
    required_columns = isnothing(c_opt) ? [:return, :logprob] : [:return, :advantage, :logprob] # Extra data columns to store
    post_batch_callback = (𝒟; kwargs...) -> nothing # Callback that that happens after sampling a batch
    loop_start_callback = (𝒮) -> nothing # Callback that happens at the beginning of each experience gathering iteration

    # Parameters specific to cost constraints (a separete value network)
    Vc::Union{ContinuousNetwork, Nothing} = nothing # Cost value approximator
    cost_opt::Union{Nothing, TrainingParams} = nothing # Training parameters for the cost value
end

function POMDPs.solve(𝒮::OnPolicySolver, mdp)
    # Construct the training buffer, constants, and sampler
    𝒟 = ExperienceBuffer(𝒮.S, 𝒮.agent.space, 𝒮.ΔN, 𝒮.required_columns, device=device(𝒮.agent.π))
    γ, λ = Float32(discount(mdp)), 𝒮.λ_gae
    s = Sampler(mdp, 𝒮.agent, S=𝒮.S, required_columns=𝒮.required_columns, λ=λ, max_steps=𝒮.max_steps, Vc=𝒮.Vc)
    isnothing(𝒮.log.sampler) && (𝒮.log.sampler = s)

    # Log the pre-train performance
    log(𝒮.log, 𝒮.i)

    # Loop over the desired number of environment interactions
    for 𝒮.i = range(𝒮.i, stop=𝒮.i + 𝒮.N - 𝒮.ΔN, step=𝒮.ΔN)
        # Call the loop start callback function
        𝒮.loop_start_callback(𝒮)

        # Sample transitions into the batch buffer
        push!(𝒟, steps!(s, Nsteps=𝒮.ΔN, reset=true, explore=true, i=𝒮.i))

        # Info to collect during training
        info = Dict()

        # Call the post-batch callback function
        𝒮.post_batch_callback(𝒟, info=info)

        # Train parameters
        for (θs, p_opt) in 𝒮.param_optimizers
            batch_train!(θs, p_opt, 𝒮.𝒫, 𝒟, info=info, π_loss=𝒮.agent.π)
        end

        # Train the actor
        batch_train!(actor(𝒮.agent.π), 𝒮.a_opt, 𝒮.𝒫, 𝒟, info=info)

        # Train the critic (if applicable)
        if !isnothing(𝒮.c_opt)
            batch_train!(critic(𝒮.agent.π), 𝒮.c_opt, 𝒮.𝒫, 𝒟, info=info)
        end

        if !isnothing(𝒮.cost_opt)
            batch_train!(𝒮.Vc, 𝒮.cost_opt, 𝒮.𝒫, 𝒟, info=info)
        end

        # Log the results
        log(𝒮.log, 𝒮.i + 1:𝒮.i + 𝒮.ΔN, info)
    end
    𝒮.i += 𝒮.ΔN
    𝒮.agent.π
end


text:
using ArrayFire
using Images
using LinearAlgebra

const blinker = [0 0 0; 1 1 1; 0 0 0]
const glider = [0 0 1; 1 0 1; 0 1 1]
const lwss = [0 1 1 1 1; 1 0 0 0 1; 0 0 0 0 1; 1 0 0 1 0]
const glidergun = [
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0;
  0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0;
  0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
  0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

function lifegame(fname, initializer, mapsize=30, imgsteps=50, upscaleratio=20)
    kernel = convert(Array{Float32}, [1 1 1; 1 0 1; 1 1 1]) |> AFArray
    initialstate = zeros(Bool, mapsize, mapsize)
    mid = div(mapsize, 2)
    (xlen, ylen), (halfx, halfy) = size(initializer), div.(size(initializer), 2)
    x1, x2 = mid - halfx, isodd(xlen) ? mid + halfx : mid + halfx - 1
    y1, y2 = mid - halfy, isodd(ylen) ? mid + halfy : mid + halfy - 1
    initialstate[x1:x2, y1:y2] .= initializer
    state = initialstate .+ Float32(.0) |> AFArray
    img = zeros(Float32, mapsize * upscaleratio, mapsize * upscaleratio, imgsteps)
    for i in 1:imgsteps
        nb = convolve2(state, kernel, UInt32(0), UInt32(0))
        a = (nb == 2)
        b = (nb == 3)
        state = ((state .* a .+ b) > 0) + Float32(0)
        frame = imresize(state, ratio=upscaleratio)
        img[:, :, i] .= frame
    end
    save(fname, img)
end

lifegame("blinker.gif", blinker)
lifegame("glider.gif", glider)
lifegame("lwss.gif", lwss)
lifegame("glidergun.gif", glidergun, 90, 200)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module HomotopyContinuation

import DoubleFloats
import DynamicPolynomials
import ElasticArrays
import FixedPolynomials
import LinearAlgebra
import MixedSubdivisions
import MultivariatePolynomials
import PrettyTables
import Printf
import ProjectiveVectors
import Random
import StaticArrays
import StaticPolynomials
import TreeViews

using Base: @propagate_inbounds
using LinearAlgebra: cond
using Parameters: @pack!, @unpack
using DoubleFloats: Double64, ComplexDF64
using DynamicPolynomials: @polyvar, subs, differentiate
using ProjectiveVectors: PVector,
                         dims,
                         dimension_indices,
                         fubini_study,
                         affine_chart,
                         ×,
                         component,
                         components,
                         combine
using StaticArrays: SVector, @SVector
using Test: @test
using MixedSubdivisions: mixed_volume

const FP = FixedPolynomials
const MP = MultivariatePolynomials
const SP = StaticPolynomials

export @polyvar, subs, differentiate
export mixed_volume
export cond
export PVector,
       dims,
       dimension_indices,
       affine_chart,
       fubini_study,
       ×,
       component,
       components,
       combine

include("progress_meter.jl")
import .ProgressMeter

include("norms.jl")
include("utilities.jl")
include("linear_algebra.jl")
include("affine_patches.jl")

include("systems_and_homotopies.jl")
include("input.jl")
include("problems.jl")
include("totaldegree.jl")

include("predictors.jl")
include("newton_corrector.jl")
include("core_tracker.jl")

include("valuation.jl")
include("cauchy_endgame.jl")
include("path_tracker.jl")
include("result.jl")
include("polyhedral.jl")
include("overdetermined.jl")
include("solver.jl")
include("monodromy.jl")
include("path_info.jl")

end


text:
#network
From=[1
1
1
2
2
3
3
4
5
6
7
8
8
9
9
10
10
11
11
12
12
13
14
15
15
15
16
16
17
17
18
19
20
21
2
3
5
4
6
9
24
9
10
10
8
9
10
11
12
11
12
13
14
13
23
23
16
16
21
24
17
19
18
22
21
20
23
22]

To=[2
3
5
4
6
9
24
9
10
10
8
9
10
11
12
11
12
13
14
13
23
23
16
16
21
24
17
19
18
22
21
20
23
22
1
1
1
2
2
3
3
4
5
6
7
8
8
9
9
10
10
11
11
12
12
13
14
15
15
15
16
16
17
17
18
19
20
21]

#Line 1-2: 175 (uncongested)
LineCapacity=[1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1500
1000
1000
1000
1500
1000
1000
1000
1000
1000
1000
1000
1000
1500
1500
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1000
1500
1000
1000
1000
1500
1000
1000
1000
1000
1000
1000
1000
1000
1500
1500
1000]
Susceptance=[68.49
4.44
11.03
7.37
4.88
7.87
11.90
9.01
10.64
15.58
15.34
5.68
5.68
11.90
11.90
11.90
11.90
20.49
23.47
20.49
10.15
11.31
16.84
58.14
40.16
18.90
38.02
42.74
69.93
9.35
75.76
49.26
89.29
14.45
68.49
4.44
11.03
7.37
4.88
7.87
11.90
9.01
10.64
15.58
15.34
5.68
5.68
11.90
11.90
11.90
11.90
20.49
23.47
20.49
10.15
11.31
16.84
58.14
40.16
18.90
38.02
42.74
69.93
9.35
75.76
49.26
89.29
14.45]


text:
module Radiation

export AbstractRadiation

    abstract type AbstractRadiation end

    using GSL
    using Interpolations
    # const GSL_CONST = GSL.Const

    # const kpc =             1e3*GSL_CONST.CGSM_PARSEC
    # const cee =             GSL_CONST.CGSM_SPEED_OF_LIGHT
    # const emgm =            GSL_CONST.CGSM_MASS_ELECTRON
    # const pmgm =            GSL_CONST.CGSM_MASS_PROTON
    # const kboltz =	        GSL_CONST.CGSM_BOLTZMANN
    # const kboltz_kev2erg =  1.6022e-9 		  # Boltzman constant in keV/erg
    # const emerg =           (GSL_CONST.CGSM_MASS_ELECTRON*GSL_CONST.CGSM_SPEED_OF_LIGHT^2)
    # const charg =           4.8e-10
    # const sigtom =          GSL_CONST.CGSM_THOMSON_CROSS_SECTION
    # const herg =            GSL_CONST.CGSM_PLANCKS_CONSTANT_H
    # const hkev =            (GSL_CONST.CGSM_PLANCKS_CONSTANT_H*6.2415e8)
    # const mjy =             1.e-26
    # const re0 =             2.81794e-13
    # const gconst =          GSL_CONST.CGSM_GRAVITATIONAL_CONSTANT
    # const sbconst =         GSL_CONST.CGSM_STEFAN_BOLTZMANN_CONSTANT
    # const aconst =          7.56e-15
    # const msun =            GSL_CONST.CGSM_SOLAR_MASS

    struct cyclosyn_emis_params
        nu::Float64
        b::Float64
        syn::Interpolations.Extrapolation
        derivs::Interpolations.Extrapolation
        # acc_syn::GSL.gsl_interp_accel
        # acc_derivs::GSL.gsl_interp_accel
    end

    struct cyclosyn_abs_params
        nu::Float64
        b::Float64
        syn::Interpolations.Extrapolation
        derivs::Interpolations.Extrapolation
        # acc_syn::gsl_interp_accel
        # acc_derivs::gsl_interp_accel
    end

    struct comint_params
        eph::Float64
        ephmin::Float64
        ephmax::Float64
        eldis::Interpolations.Extrapolation
        phodis::Interpolations.Extrapolation
        # acc_eldis::gsl_interp_accel
        # acc_phodis::gsl_interp_accel
    end

    struct comfnc_params
        game::Float64
        e1::Float64
        phodis::Interpolations.Extrapolation
        # acc_phodis::gsl_interp_accel
    end

    struct disk_obs_params
        tin::Float64
        rin::Float64
        nu::Float64
    end

    struct disk_ic_params
        gamma::Float64
        beta::Float64
        tin::Float64
        rin::Float64
        rout::Float64
        h::Float64
        z::Float64
        nu::Float64
    end

    function get_energy(rad::AbstractRadiation)     return rad.en_phot end
    function get_nphot(rad::AbstractRadiation)      return num_phot end
    function get_energy_obs(rad::AbstractRadiation) return en_phot_obs end
    function get_nphot_obs(rad::AbstractRadiation)  return num_phot_obs end
    function get_size(rad::AbstractRadiation)       return size end

    function set_beaming!(rad::AbstractRadiation, θ, speed, doppler)
        rad.angle = deg2rad(theta)
        rad.beta = speed
        rad.dopfac = doppler
    end

    function set_inclination!(rad::AbstractRadiation, θ)
        rad.angle = deg2rad(θ)
    end

    function set_geometry!(rad::AbstractRadiation, geom, l1, l2)

        if isequal(geom, "cylinder")
            r = l1
            z = l2
            vol = π*r^2*z
            dopnum = 2.0
        elseif isequal(geom, "sphere")
            r = l1
            z = r
            vol = (4/3)*π*r^3
            dopnum = 3.0
        else
            @warn "Input wrong, assuming sphere of radius 1 cm. \n Choose either sphere or cylinder"
            r = 1.0
            z = 1.0
            vol = 4π/3
            dopnum = 3.0
        end

        rad.r = r
        rad.z = z
        rad.vol = vol
        rad.dopnum = dopnum
        rad.geometry = geom
    end

    function integrated_lumunosity!(rad::AbstractRadiation, numin, numax)
        temp = 0.0
        for i ∈ 1:rad.size - 1
            if ((rad.en_phot_obs[i]  /rad.herg > numin) &&
                (rad.en_phot_obs[i+1]/rad.herg < numax))

                temp = temp + 0.5*(rad.en_phot_obs[i+1]/rad.herg - rad.en_phot_obs[i]/rad.herg)*
                                (rad.num_phot_obs[i+1] + rad.num_phot_obs[i])
            end
        end
        return temp
    end

    function set_counterjet!(rad::AbstractRadiation, flag)
        rad.counterjet = flag
    end

end


text:
using DataStructures

include("types.jl")

"""
Depth First Search.
"""
function dfs(val, target, s_fns, max_depth)
    path = []
    nodes = Stack{TreeNode}()
    push!(nodes, TreeNode(val, 1))

    while !isempty(nodes)
        node = pop!(nodes)

        if node.depth > max_depth
            continue
        end

        push!(path, node.value)

        if node.value == target
            return node
        end

        for fn in s_fns
            new_val = fn(node)
            new_node = TreeNode(new_val, node.depth+1)
            push!(nodes, new_node)
        end
    end
end

"""
Breadth First Search.
"""
function bfs(val, target, s_fns)
    path = []
    nodes = Queue{TreeNode}()
    enqueue!(nodes, TreeNode(val, 1))

    while !isempty(nodes)
        node = dequeue!(nodes)

        push!(path, node.value)

        if node.value == target
            return node
        end

        for fn in s_fns
            new_val = fn(node)
            new_node = TreeNode(new_val, node.depth+1)
            enqueue!(nodes, new_node)
        end
    end
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
import Base: done, next, start

#### Iterator ####
# An n-dimensional grid iterator. This allows you to do things in
# arbitrary dimensions that in fixed dimensions you might do with a
# comprehension.
#
# For example, for two dimensions the equivalent of
#    p = [ i*j for i = 1:m, j = 1:n ]
# is
#    sz = (m, n)
#    p = Array(Int, sz)
#    index = 1
#    for c = Counter(sz)
#      p[index] = prod(c)
#      index += 1
#    end
# However, with sz appropriately defined, this version works for
# arbitrary dimensions.

type Counter
    max::Vector{Int}
end
Counter(sz::Tuple) = Counter(Int[sz...])

function start(c::Counter)
    N = length(c.max)
    state = ones(Int,N)
    if N > 0
        state[1] = 0 # because of start/done/next sequence, start out one behind
    end
    return state
end
function done(c::Counter, state)
    if isempty(state)
        return true
    end
    # we do the increment as part of "done" to make exit-testing more efficient
    state[1] += 1
    i = 1
    max = c.max
    while state[i] > max[i] && i < length(state)
        state[i] = 1
        i += 1
        state[i] += 1
    end
    state[end] > max[end]
end
next(c::Counter, state) = state, state


text:
function moon_position(time::Float64)
    moon_angle = MOON_INITIAL_ANGLE + 2.0pi * time / MOON_PERIOD
    x::Float64 = DISTANCE_TO_MOON * cos(moon_angle)
    y::Float64 = DISTANCE_TO_MOON * sin(moon_angle)

    return [x, y]
end

function update(me::Moon, time)
    me.position = moon_position(time)

    me
end

max_stars_repo_name: rodrigolece/grupo-fisica-compu

max_issues_repo_name: rodrigolece/grupo-fisica-compu

max_forks_repo_name: rodrigolece/grupo-fisica-compu


text:
using Cameras
using Test

@testset "Acquisition" begin
    image_size = (2, 2)

    function produce_images!(image_source::Channel)
        while true
            put!(image_source, zeros(image_size))
        end
    end
    image_source = Channel(produce_images!)

    @testset "Manually triggered" begin
        camera = SimulatedCamera(image_source)

        @test !isopen(camera)
        open!(camera)
        @test isopen(camera)

        @test !isrunning(camera)
        start!(camera)
        @test isrunning(camera)

        @testset "Synchronous" begin
            trigger!(camera)
            img = take!(camera)
            @test size(img) == image_size
        end

        @testset "Asynchronous" begin
            let img = zeros(0, 0)
                @sync begin
                    @async img = take!(camera)
                    @async trigger!(camera)
                end
                @test size(img) == image_size
            end
        end

        @assert isrunning(camera)
        stop!(camera)
        @test !isrunning(camera)

        @assert isopen(camera)
        close!(camera)
        @test !isopen(camera)
    end

    @testset "Continuously triggered" begin
        period = 0.04 # seconds

        function produce_triggers!(trigger_source::Channel{UInt64})
            while true
                put!(trigger_source, time_ns())
                sleep(period)
            end
        end
        trigger_source = Channel(produce_triggers!; ctype = UInt64)

        continuous_camera = SimulatedCamera(trigger_source, image_source)

        start!(continuous_camera)

        let i = 0
            while i < 5
                i += 1
                img = take!(continuous_camera)
                @assert size(img) == image_size
            end
            @test i == 5
        end
    end
end


text:
####    GridOperators

# Local finite-difference operator functions.

export advect
export grad, ∇_pos, ∇_neg

export ∇_z_flux, ∇_z_centered, ∇_z_dual, ∇_z_upwind
export Δ_z, Δ_z_dual

export advect_old
export UpwindAdvective,
       UpwindCollocated,
       OneSidedUp,
       OneSidedDn,
       CenteredUnstable

function advect(f, w, grid::Grid)
  @assert length(f)==3
  @assert length(w)==3
  if w[2] < 0
    return (f[3]-f[2])*grid.Δzi
  else
    return (f[2]-f[1])*grid.Δzi
  end
end

"""
    grad(f, grid::Grid{FT}) where FT

Computes local finite-difference gradient of `f`: `∇f`
"""
function grad(f, grid::Grid{FT}) where FT
  @assert length(f)==3 || length(f)==2
  if length(f)==2
    return (f[2]-f[1])*grid.Δzi
  else
    return (f[3]-f[1])*FT(0.5)*grid.Δzi
  end
end

"""
    ∇_pos(f, grid::Grid)

Computes a one-sided (up) local finite-difference gradient of `f`: `∇f`
"""
function ∇_pos(f, grid::Grid)
  @assert length(f)==3
  return (f[2]-f[1])*grid.Δzi
end

"""
    ∇_neg(f, grid::Grid)

Computes a one-sided (down) local finite-difference gradient of `f`: `∇f`
"""
function ∇_neg(f, grid::Grid)
  @assert length(f)==3
  return (f[3]-f[2])*grid.Δzi
end

"""
    ∇_z_flux(f, grid::Grid)

Computes the local derivative of field ``f``, which is
assumed to live on the dual grid: ``∇f``
"""
function ∇_z_flux(f, grid::Grid)
  @assert length(f)==2
  (f[2] - f[1])*grid.Δzi
end

"""
    ∇_z_centered(f, grid::Grid)

Computes the local derivative of field ``f``, which is
assumed to live on the dual grid: ``∇f``
"""
function ∇_z_centered(f, grid::Grid)
  @assert length(f)==3
  (f[3] - f[1])*grid.Δzi/2
end

"""
    ∇_z_dual(f, grid::Grid)

Computes the local derivative of field ``f``, which is
assumed to live on the dual grid: ``∇f``
"""
function ∇_z_dual(f, grid::Grid)
  @assert length(f)==3
  return [∇_z_flux(f[1:2], grid), ∇_z_flux(f[2:3], grid)]
end

"""
    ∇_z_upwind(ϕ, w, grid::Grid)

Local upwind gradient operator ``∇ϕ``. This
operator is stable but numerically diffusive.
"""
function ∇_z_upwind(ϕ, w, grid::Grid)
  @assert length(ϕ)==3
  @assert length(w)==3
  return w[2]>0 ? ∇_z_flux(ϕ[1:2], grid) : ∇_z_flux(ϕ[2:3], grid)
end

"""
    Δ_z(f, grid::Grid)

Computes the local Laplacian of field ``f``:
``∇ • (∇f)``
"""
function Δ_z(f, grid::Grid)
  @assert length(f)==3
  (f[3] - 2*f[2] + f[1])*grid.Δzi2
end

"""
    Δ_z(f, grid::Grid, K)

Computes the local Laplacian of field ``f``
with a variable coefficient:
``∇ • (K ∇f)``
"""
function Δ_z(f, grid::Grid, K)
  @assert length(f)==3
  K∇f_M = (K[1]+K[2])/2*(f[2]-f[1])*grid.Δzi
  K∇f_P = (K[2]+K[3])/2*(f[3]-f[2])*grid.Δzi
  return (K∇f_P - K∇f_M)*grid.Δzi
end

"""
    Δ_z_dual(f, grid::Grid, K)

Computes the local Laplacian of field ``f``
with a variable coefficient:
``∇ • (K ∇f)``
"""
function Δ_z_dual(f, grid::Grid, K)
  @assert length(f)==3
  @assert length(K)==2
  K∇f_M = K[1]*(f[2]-f[1])*grid.Δzi
  K∇f_P = K[2]*(f[3]-f[2])*grid.Δzi
  return (K∇f_P - K∇f_M)*grid.Δzi
end

abstract type AdvectForm end
struct UpwindAdvective <: AdvectForm end

abstract type ConservativeForm end
struct UpwindCollocated <: ConservativeForm end
struct OneSidedUp <: ConservativeForm end
struct OneSidedDn <: ConservativeForm end
struct CenteredUnstable <: ConservativeForm end

"""
    advect(ϕ, w, grid::Grid, ::AdvectForm, Δt)

Advection operator, which computes

 - ``w • ∇ϕ`` given `AdvectiveForm`
 - ``∇ • (wϕ)`` given `ConservativeForm`

which is dispatched by `AdvectForm`.
"""
function advect_old(ϕ, ϕ_dual, w, w_dual, grid::Grid, ::UpwindAdvective, Δt)
  @assert length(ϕ)==3
  @assert length(ϕ_dual)==2
  @assert length(w)==3
  @assert length(w_dual)==2
  if w[2] > 0
    return w[ 2 ]*(ϕ[2] - ϕ[1]) * grid.Δzi
  else
    return w[ 2 ]*(ϕ[3] - ϕ[2]) * grid.Δzi
  end
end
function advect_old(ϕ, ϕ_dual, w, w_dual, grid::Grid, ::UpwindCollocated, Δt)
  @assert length(ϕ)==3
  @assert length(ϕ_dual)==2
  @assert length(w)==3
  @assert length(w_dual)==2
  if sum(w)/length(w) > 0
    return (w[2]*ϕ[2] - w[1]*ϕ[1]) * grid.Δzi
  else
    return (w[3]*ϕ[3] - w[2]*ϕ[2]) * grid.Δzi
  end
end
function advect_old(ϕ, ϕ_dual, w, w_dual, grid::Grid, ::OneSidedUp, Δt)
  @assert length(ϕ)==3
  @assert length(ϕ_dual)==2
  @assert length(w)==3
  @assert length(w_dual)==2
  return (w[2]*ϕ[2] - w[1]*ϕ[1]) * grid.Δzi
end
function advect_old(ϕ, ϕ_dual, w, w_dual, grid::Grid, ::OneSidedDn, Δt)
  @assert length(ϕ)==3
  @assert length(ϕ_dual)==2
  @assert length(w)==3
  @assert length(w_dual)==2
  return (w[3]*ϕ[3] - w[2]*ϕ[2]) * grid.Δzi
end
function advect_old(ϕ, ϕ_dual, w, w_dual, grid::Grid, ::CenteredUnstable, Δt)
  @assert length(ϕ)==3
  @assert length(ϕ_dual)==2
  @assert length(w)==3
  @assert length(w_dual)==2
  return (w_dual[2]*ϕ_dual[2]-w_dual[1]*ϕ_dual[1])*grid.Δzi
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Multivector.jl defines the Multivector type and interface
"""

# --- Exports

# ------ Types

export Multivector

# --- Type definitions

using DataStructures: SortedDict

"""
    struct Multivector{T<:AbstractFloat} <: AbstractMultivector

Multivector represented with the floating-point precision of type `T`.
"""
struct Multivector{T<:AbstractFloat} <: AbstractMultivector{T}
    #=
      Fields
      ------
      * `parts`: collection of k-vectors that sum to the multivector

      * `norm`: norm of the multivector
    =#
    dim::Int
    parts::SortedDict{Int, Vector{AbstractBlade}}
    norm::T

    """
    Construct a Multivector from a of vector of blades. For each grade ``k``,
    the blades used to represent the ``k``-vector part of the multivector form
    an orthogonal basis for the subspace of ``k``-vectors.
    """
    function Multivector{T}(
            blades::Vector{<:AbstractBlade}) where {T<:AbstractFloat}

        # --- Check arguments

        # Check that all non-scalar blades have the same dimension
        dim_ = nothing
        for B in blades
            if !(B isa AbstractScalar)
                if isnothing(dim_)
                    dim_ = dim(B)
                else
                    if dim_ != dim(B)
                        message = "Non-scalar blades in `blades` have " *
                                  "differing dimensions."
                        throw(DimensionMismatch(message))
                    end
                end
            end
        end

        # --- Handle edge cases

        # Return Scalar(0) if number of blades is 0.
        if length(blades) == 0
            return zero(Blade{T})
        end

        # Return a Scalar if all of the blades are scalars
        blades_all_scalar = true
        for B in blades
            if !(B isa AbstractScalar)
                blades_all_scalar = false
                break
            end
        end
        if blades_all_scalar
            return Scalar{T}(value(reduce(+, blades)))
        end

        # --- Construct Multivector

        # Construct `parts`. Sort blades by grade.
        parts = SortedDict{Int, Vector{AbstractBlade}}()
        for B in blades
            # Ignore Scalar(0)
            if volume(B) == 0
                continue
            end

            k_vectors = get!(parts, grade(B), [])
            push!(k_vectors, convert(AbstractBlade{T}, B))
        end

        # --- Reduce multivector

        # Reduce scalar, vector, and pseudoscalar parts
        for k in [0, 1, dim_]
            if k in keys(parts)
                parts[k] = [reduce(+, parts[k])]
            end
        end

        # TODO: reduce remaining k-vectors

        # --- Compute norm

        norm_ = 0
        for k_vectors in parts
            norm_ += reduce(+, map(B -> norm(B)^2, k_vectors))
        end
        norm_ = sqrt(norm_)

        # --- Construct result

        # Return AbstractBlade if there is only one term
        if length(parts) == 1
            k_vectors = first(values(parts))
            if length(k_vectors) == 1
                return first(k_vectors)
            end
        end

        # Return new Multivector
        new(dim_, parts, norm_)
    end

    """
    Type conversion constructor. Construct a copy of the Multivector converted
    to the specified precision.
    """
    Multivector{T}(M::Multivector) where {T<:AbstractFloat} =
        T == typeof(norm(M)) ?
            M :
            Multivector{T}(reduce(vcat, blades(M)))
end

"""
    Multivector(blades::Vector{<:AbstractBlade})

Construct a Multivector from a of vector of blades. For each grade ``k``, the
blades used to represent the ``k``-vector part of the multivector form an
orthogonal basis for the subspace of ``k``-vectors.

The precision of the Multivector is inferred from the precision of the first
element of `blades`.
"""
Multivector(blades::Vector{<:AbstractBlade}) =
    length(blades) > 0 ?
        Multivector{typeof(volume(blades[1]))}(blades) :
        Multivector{Float64}(blades)

Multivector(multivectors::Vector{<:AbstractMultivector}) =
    length(multivectors) > 0 ?
        Multivector(reduce(vcat, map(blades, multivectors))) :
        Multivector{Float64}(Vector{AbstractBlade}())

# --- Method definitions for AbstractMultivector interface functions

dim(M::Multivector) = M.dim

grades(M::Multivector; collect=true) =
    collect ?
        Base.collect(keys(M.parts)) :
        keys(M.parts)

blades(M::Multivector) = reduce(vcat, values(M.parts))

norm(M::Multivector) = M.norm

Base.getindex(M::Multivector, k::Int) =
    k in grades(M) ?  M.parts[k] : Vector{AbstractBlade}()

# --- Utility methods

convert(::Type{T}, M::Multivector) where {T<:AbstractFloat} =
    T == typeof(norm(M)) ? M : Multivector{T}(M)


text:
function as_sub(x::AbstractVector)
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = sub(y, 2:2:length(y))
    y[:] = x[:]
    y
end
function as_sub(x::AbstractMatrix)
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = sub(y, 2:2:size(y,1), 2:2:size(y,2))
    for j=1:size(x,2)
        for i=1:size(x,1)
            y[i,j] = x[i,j]
        end
    end
    y
end
function as_sub{T}(x::AbstractArray{T,3})
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = sub(y, 2:2:size(y,1), 2:2:size(y,2), 2:2:size(y,3))
    for k=1:size(x,3)
        for j=1:size(x,2)
            for i=1:size(x,1)
                y[i,j,k] = x[i,j,k]
            end
        end
    end
    y
end

bittest(f::Function, ewf::Function, a...) = (@test ewf(a...) == bitpack(broadcast(f, a...)))
n1 = 21
n2 = 32
n3 = 17
rb = 1:5

for arr in (identity, as_sub)
    @test broadcast(+, arr(eye(2)), arr([1, 4])) == [2 1; 4 5]
    @test broadcast(+, arr(eye(2)), arr([1  4])) == [2 4; 1 5]
    @test broadcast(+, arr([1  0]), arr([1, 4])) == [2 1; 5 4]
    @test broadcast(+, arr([1, 0]), arr([1  4])) == [2 5; 1 4]
    @test broadcast(+, arr([1, 0]), arr([1, 4])) == [2, 4]
    @test broadcast(+, arr([1, 0]), 2) == [3, 2]

    @test @inferred(arr(eye(2)) .+ arr([1, 4])) == arr([2 1; 4 5])
    @test arr(eye(2)) .+ arr([1  4]) == arr([2 4; 1 5])
    @test arr([1  0]) .+ arr([1, 4]) == arr([2 1; 5 4])
    @test arr([1, 0]) .+ arr([1  4]) == arr([2 5; 1 4])
    @test arr([1, 0]) .+ arr([1, 4]) == arr([2, 4])
    @test arr([1]) .+ arr([]) == arr([])

    A = arr(eye(2)); @test broadcast!(+, A, A, arr([1, 4])) == arr([2 1; 4 5])
    A = arr(eye(2)); @test broadcast!(+, A, A, arr([1  4])) == arr([2 4; 1 5])
    A = arr([1  0]); @test_throws ErrorException broadcast!(+, A, A, arr([1, 4]))
    A = arr([1  0]); @test broadcast!(+, A, A, arr([1  4])) == arr([2 4])
    A = arr([1  0]); @test broadcast!(+, A, A, 2) == arr([3 2])

    @test arr([ 1    2])   .* arr([3,   4])   == [ 3 6; 4 8]
    @test arr([24.0 12.0]) ./ arr([2.0, 3.0]) == [12 6; 8 4]
    @test arr([1 2]) ./ arr([3, 4]) == [1/3 2/3; 1/4 2/4]
    @test arr([1 2]) .\ arr([3, 4]) == [3 1.5; 4 2]
    @test arr([3 4]) .^ arr([1, 2]) == [3 4; 9 16]
    @test arr(bitpack([true false])) .* arr(bitpack([true, true])) == [true false; true false]
    @test arr(bitpack([true false])) .^ arr(bitpack([false, true])) == [true true; true false]
    @test arr(bitpack([true false])) .^ arr([0, 3]) == [true true; true false]

    M = arr([11 12; 21 22])
    @test broadcast_getindex(M, eye(Int, 2).+1,arr([1, 2])) == [21 11; 12 22]
    @test_throws BoundsError broadcast_getindex(M, eye(Int, 2).+1,arr([1, -1]))
    @test_throws BoundsError broadcast_getindex(M, eye(Int, 2).+1,arr([1, 2]), [2])
    @test broadcast_getindex(M, eye(Int, 2).+1,arr([2, 1]), [1]) == [22 12; 11 21]

    A = arr(zeros(2,2))
    broadcast_setindex!(A, arr([21 11; 12 22]), eye(Int, 2).+1,arr([1, 2]))
    @test A == M
    broadcast_setindex!(A, 5, [1,2], [2 2])
    @test A == [11 5; 21 5]
    broadcast_setindex!(A, 7, [1,2], [1 2])
    @test A == fill(7, 2, 2)
    A = arr(zeros(3,3))
    broadcast_setindex!(A, 10:12, 1:3, 1:3)
    @test A == diagm(10:12)
    @test_throws BoundsError broadcast_setindex!(A, 7, [1,-1], [1 2])

    for (f, ewf) in (((==), (.==)),
                     ((<) , (.<) ),
                     ((!=), (.!=)),
                     ((<=), (.<=)))
        bittest(f, ewf, arr(eye(2)), arr([1, 4]))
        bittest(f, ewf, arr(eye(2)), arr([1  4]))
        bittest(f, ewf, arr([0, 1]), arr([1  4]))
        bittest(f, ewf, arr([0  1]), arr([1, 4]))
        bittest(f, ewf, arr([1, 0]), arr([1, 4]))
        bittest(f, ewf, arr(rand(rb, n1, n2, n3)), arr(rand(rb, n1, n2, n3)))
        bittest(f, ewf, arr(rand(rb,  1, n2, n3)), arr(rand(rb, n1,  1, n3)))
        bittest(f, ewf, arr(rand(rb,  1, n2,  1)), arr(rand(rb, n1,  1, n3)))
        bittest(f, ewf, arr(randbool(n1, n2, n3)), arr(randbool(n1, n2, n3)))
    end
end

r1 = 1:1
r2 = 1:5
ratio = [1,1/2,1/3,1/4,1/5]
@test r1.*r2 == [1:5]
@test r1./r2 == ratio
m = [1:2]'
@test m.*r2 == [1:5 2:2:10]
@test_approx_eq m./r2 [ratio 2ratio]
@test_approx_eq m./[r2] [ratio 2ratio]

@test @inferred([0,1.2].+reshape([0,-2],1,1,2)) == reshape([0 -2; 1.2 -0.8],2,1,2)
rt = Base.return_types(.+, (Array{Float64, 3}, Array{Int, 1}))
@test length(rt) == 1 && rt[1] == Array{Float64, 3}
rt = Base.return_types(broadcast, (Function, Array{Float64, 3}, Array{Int, 1}))
@test length(rt) == 1 && rt[1] == Array{Float64, 3}
rt = Base.return_types(broadcast!, (Function, Array{Float64, 3}, Array{Float64, 3}, Array{Int, 1}))
@test length(rt) == 1 && rt[1] == Array{Float64, 3}

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# pitch and z case

using Caesar
using RoME
using TransformUtils

# modify ROV init orientation
using CoordinateTransformations
using Rotations

# start visualizer with some scene
vc = startdefaultvisualization(draworigin=true)
sleep(3.0)

# load necessary data
rovt = loadmodel(:rov)
rovt(vc)

am = Translation(0,0,0.0) ∘ LinearMap(Rotations.AngleAxis(pi/2,0,0,1.0))
rovt(vc, am )

# defaultscene01!(vc)
sc01 = loadmodel(:scene01)
sc01(vc)

# ground truth positions
gt = Dict{Symbol, Vector{Float64}}()
gt[:x1] = [0.0;0.0;0;0;0;pi/2]
gt[:x2] = [0.0;4.0-4.0/sqrt(2);4.0/sqrt(2);0;0;pi/2]
gt[:l1] = [0.0;4.0;0.7]

# prepare for animation
as = ArcPointsRangeSolve(
[0.0;0;0.0],
[0.0;4.0-4.0/sqrt(2);4.0/sqrt(2)],
[0.0;4.0;4.0], 4.0 )
findaxiscenter!(as)

# covariance parameters
initCov = 0.01*eye(6); [initCov[i,i] = 0.001 for i in 4:6];
odoCov = 0.001*eye(6); [odoCov[i,i] = 0.001 for i in 4:6];
rangecov, bearingcov = 3e-4, 2e-3
tfx1 = SE3(gt[:x1][1:3],Euler(gt[:x1][4:6]) )

sleep(3)
# factor graph
N=300
fg = identitypose6fg(initpose=tfx1, initCov=initCov, N=N)

visualizeallposes!(vc, fg, drawlandms=false)

rho1 = norm(tfx1.t-gt[:l1])
addLinearArrayConstraint(fg, (rho1, 0.0), :x1, :l1, rangecov=rangecov,bearingcov=bearingcov)

visualizeDensityMesh!(vc, fg, :l1)

# animate visualization
# @async begin

initrot=Rotations.AngleAxis(pi/2,0,0,1.0)
animatearc(vc, rovt, as, initrot=initrot, delaytime=0.02,to=0.5)
# end

am = parameterizeArcAffineMap(0.5, as, initrot=initrot )
am = LinearMap(Rotations.AngleAxis(-pi/2,0,0,1.0)) ∘ am
px2 = SE3(am.v[1:3] ,Quaternion(am.m.w, [am.m.x;am.m.y;am.m.z] ))
@show px2.t

addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(px2), odoCov) ) )

@show rho2 = norm((tfx1*px2).t-gt[:l1])
addLinearArrayConstraint(fg, (rho2, 0.0), :x2, :l1, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (rho2, 0.0), :x2, :l2, rangecov=rangecov,bearingcov=bearingcov)
visualizeDensityMesh!(vc, fg, :l2)
sleep(3)

batchSolve(fg)
visualize(fg, vc, densitymeshes=[:l1, :l2], N=N, drawlandms=false)

# should just be the point where landmark likelihoods sightings intersect
deletemeshes!(vc)
visualizeDensityMesh!(vc, fg, :l1)

#


text:
function read_stdin_all()
  text = ""

  for line in readlines()
    text *= line
    text *= "\n"
  end

  text
end

function find_index(xs, x)
  for i = firstindex(xs):lastindex(xs)
    try
      if xs[i] == x
        return i
      end
      i += 1
    catch
      # ignore
    end
  end

  0
end

function is_digit(c)
  '0' <= c && c <= '9'
end

function match_int(rest)
  bi = 1

  if rest[1] == '-'
    bi = 2
  end

  while bi <= lastindex(rest)
    c = rest[bi]
    if is_digit(c)
      bi += 1
    else
      break
    end
  end

  bi - 1
end

function match_str(str)
  if str[1] != '"'
    return -1
  end

  find_index(str[2:end], '"') - 1
end


text:
####################################################################
function gpmodelVI(x, flux, σ2, idx; maxiter = 10, seed = 1, show_trace = true)
####################################################################

    #--------------------------------------------------------
    # Auxiliary
    #--------------------------------------------------------

    rg = MersenneTwister(seed)

    JITTER = 1e-8

    Σdiag = Diagonal(σ2)

    numF  = length(unique(idx))

    @assert(all(sort(unique(idx)) .== collect(1:numF)))

    N = length(x)

    @assert(N == length(flux) == length(idx) == length(σ2))

    @printf("Running gpmodelVIexact with %d and %d filters\n", N, numF)

    #--------------------------------------------------------
    # Initialise GP hyperparameters with ML estimation
    #--------------------------------------------------------

    a, b, logs = gpmodel(x, flux, σ2, idx; maxiter = 10_000, seed = 1, maxrandom=1000,show_trace = show_trace)

    #--------------------------------------------------------
    # Calculate kernel matrix and keep fixed!
    #--------------------------------------------------------

    dist2 = pairwise(SqEuclidean(), reshape(x, 1, N), dims=2)

    K1 = zeros(N, N)

    calculatekernelmatrix!(K1, dist2, rbf, [0.0; logs])

    K1 = K1 + JITTER * I

    K = PDMat((K1 + K1') / 2) # This helps us speed calculations

    #--------------------------------------------------------------------
    # Define functions for calling lower bound.
    # Actual calculation takes place in function "complete_lower_bound"
    #--------------------------------------------------------------------

    # Auxiliary function that calls the complete lower bound

    function lb(μα, Σα, μβ, Σβ, μf, Σf)

        local v²ₐ = v²ᵦ = 1000.0

        complete_lower_bound(; y=flux, idx=idx, S=Σdiag, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    end

    # Unpacks vectorised parameters

    function unpack(param)

        @assert(length(param) == numF*4)

        local μα =               param[1+0*numF:1*numF]

        local Σα = Diagonal(exp.(param[1+1*numF:2*numF]))

        local μβ =               param[1+2*numF:3*numF]

        local Σβ = Diagonal(exp.(param[1+3*numF:4*numF]))

        return μα, Σα, μβ, Σβ

    end

    # Constructs necessary matrices

    function negative_lb_helper(param)

        local μα, Σα, μβ, Σβ = unpack(param)

        local Ā  = Diagonal(μα[idx])

        local Σa = Diagonal(Σα.diag[idx])

        local b̄  = μβ[idx]

        # Calculate optimal posterior q(f)

        # local Σf = ((Ā*inv(Σdiag)*Ā + inv(Σdiag)*Σa) + inv(K))\I

        local invDiagPart = (Ā*(Σdiag\Ā) + Σdiag\Σa)\I #
                                                       # This is faster/more stable
        local Σf = invDiagPart*((invDiagPart + K)\K)   #

        local μf = Σf * Ā * (Σdiag \ (flux - b̄))

        -1.0 * lb(μα, Σα, μβ, Σβ, μf, Σf)

    end

    # Initial solution given by point estimates

    initsol = [a; log.(1.0*ones(numF)); b; log.(1.0*ones(numF))]

    # Use the LBFGS optimiser

    # opt = Optim.Options(show_trace = show_trace, iterations = maxiter, show_every = 1)
    #
    # result = Optim.optimize(negative_lb_helper, initsol, LBFGS(), opt, autodiff=:forward)

    # result = bboptimize(negative_lb_helper; SearchRange=(-10.0,10.0), NumDimensions = numF*4, Method=:xnes)

    opt = NLopt.Opt(:LN_BOBYQA, numF*4)

    opt.xtol_rel = 1e-8

    COUNTER, BESTF = 0, Inf

    function myfunc(x, dummy)

        COUNTER +=1

        F = negative_lb_helper(x)

        BESTF = min(F, BESTF)

        mod(COUNTER, 100)==0 ? @printf("Iter %d f=%e\n", COUNTER, BESTF) : nothing

        return F

    end

    opt.min_objective = myfunc

    (_minf, minx, _ret) = NLopt.optimize(opt, initsol)

    # Retrieve result and return posterior distribution of type Distributions.MvNormal

    μα, Σα, μβ, Σβ = unpack(minx)#result.minimizer)

    return MvNormal([μα; μβ], Diagonal([Σα.diag; Σβ.diag]))

end


text:
function welcome_message()

    println("___________________________________________________________")
    println(" ")
    println(" Welcome to the Zilindroa Code")
    println("___________________________________________________________")

end


text:
# Shock capturing using Prof. Hicken's (as yet unnamed) projection method,
#
function calcShockCapturing(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             sensor::AbstractShockSensor,
                             capture::ProjectionShockCapturing)

  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    coords_i = ro_sview(mesh.coords, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i)
    jac_i = ro_sview(mesh.jac, :, i)
    res_i = sview(eqn.res, :, :, i)

    projectionShockCapturing(eqn.params, sbp, sensor, capture, q_i, i, coords_i,
                             dxidx_i, jac_i, res_i)
  end

  return nothing
end

"""
Returns a nxn matrix that projects out modes lower than sbp.degree
"""
function getFilterOperator!(sbp::TriSBP{T}, diss::AbstractArray{T,2}) where {T}
  x = calcnodes(sbp)
  # loop over ortho polys up to degree d
  #println("TEMP: filter operator set to sbp.degree-1!!!!!")
  d = 0
#  d = sbp.degree
  V = zeros(T, (sbp.numnodes, convert(Int, (d+1)*(d+2)/2)) )
  ptr = 0
  for r = 0:d
    for j = 0:r
      i = r-j
      V[:,ptr+1] = SummationByParts.OrthoPoly.proriolpoly(vec(x[1,:]), vec(x[2,:]), i, j)
      ptr += 1
    end
  end
  #F = eye(sbp.numnodes) - (pinv(V.')*V.').'
  #Minv = diagm(1./diag(V.'*diagm(sbp.w)*V))
  # Minv is the idenity, since the V's are orthogonal in L2
  diss[:,:] = eye(sbp.numnodes) - (V*V.')*diagm(sbp.w)
end

#TODO: this doesn't work yet
function getFilterOperator!(sbp::TetSBP{T}, diss::AbstractArray{T,2}) where {T}

  fill!(diss, 0)
end

function projectionShockCapturing(params::ParamType, sbp::AbstractOperator,
                             sensor::AbstractShockSensor,
                             capture::ProjectionShockCapturing,
                             q::AbstractMatrix, elnum::Integer,
                             coords::AbstractMatrix,
                             dxidx::Abstract3DArray, jac::AbstractVector,
                             res::AbstractMatrix{Tres}) where {Tres}

  numDofPerNode, numNodesPerElement = size(q)
  dim = size(dxidx, 1)
  @unpack capture t1 t2 w

  #TODO: stash these somewhere
  Se = zeros(Tres, dim, numNodesPerElement)
  ee = zeros(Tres, dim, numNodesPerElement)
  is_nonconst = getShockSensor(params, sbp, sensor, q, elnum, coords, dxidx,
                               jac, Se, ee)
  if is_nonconst # if there is a shock
    # For scalar equations, the operator is applied -epsilon * P^T M P * u
    # For vector equations, P needs to be applied to all equations as once:
    # utmp^T = P*u^T
    # Instead, work with utmp = (P*u^T)^T = u*P^T

    # convert to entropy variables to make this term entropy-stable
    #w = zeros(eltype(q), numDofPerNode, numNodesPerElement)
    for i=1:numNodesPerElement
      w_i = sview(w, :, i)
      q_i = sview(q, :, i)
      convertToIR(params, q_i, w_i)
    end

    # apply P
    smallmatmatT!(w, capture.filt, t1)

    # apply mass matrix
    @simd for i=1:numNodesPerElement
      fac = sbp.w[i]/jac[i]
      @simd for j=1:numDofPerNode
        t1[j, i] *= fac
      end
    end

    # apply P^T
    smallmatmat!(t1, capture.filt, t2)
    @simd for i=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        # this uses only the x direction viscoscity
        res[j, i] -= ee[1, i]*t2[j, i]
      end
    end

  end  # end if

  return nothing
end

#------------------------------------------------------------------------------
# Differentiated version

function calcShockCapturing_diff(mesh::AbstractMesh, sbp::AbstractOperator,
                                 eqn::EulerData, opts,
                                 sensor::AbstractShockSensor,
                                 capture::ProjectionShockCapturing,
                                 assem::AssembleElementData)

  data = eqn.params.calc_volume_integrals_data
  res_jac = data.res_jac

  for i=1:mesh.numEl
    q_i = ro_sview(eqn.q, :, :, i)
    coords_i = ro_sview(mesh.coords, :, :, i)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i)
    jac_i = ro_sview(mesh.jac, :, i)

    fill!(res_jac, 0)
    nonzero_jac = projectionShockCapturing_diff(eqn.params, sbp, sensor,
                              capture, q_i, i, coords_i, dxidx_i, jac_i, res_jac)

    # assembling into a sparse matrix is non-trivially expensive, don't do
    # it unless this element has shock capturing active
    if nonzero_jac
      if eqn.params.use_Minv == 1
        applyMinvElement(jac_i, sbp.w, res_jac)
      end

      # assemble element level jacobian into the residual
      assembleElement(assem, mesh, i, res_jac)
    end  # if nonzero_jac
  end  # end i

  return nothing
end

"""
  Differentiated version of `projectionShockCapturing` for
  [`ProjectionShockCapturing`](@ref).
"""
function projectionShockCapturing_diff(params::ParamType, sbp::AbstractOperator,
                                  sensor::AbstractShockSensor,
                                  capture::ProjectionShockCapturing,
                                  u::AbstractMatrix, elnum::Integer,
                                  coords::AbstractMatrix,
                                  dxidx::Abstract3DArray,
                                  jac::AbstractVector{Tmsh},
                                  res_jac::AbstractArray{Tres, 4}) where {Tmsh, Tres}

  numDofPerNode, numNodesPerElement = size(u)
  dim = size(coords, 1)
  @unpack capture t1 t2 w Se_jac ee_jac A0inv

  #TODO: stash these somewhere
  Se = zeros(Tres, dim, numNodesPerElement)
  ee = zeros(Tres, dim, numNodesPerElement)
  #TODO: make shock capturing and shock sensing independent choices
  is_shock = isShockElement(params, sbp, sensor, u, elnum, coords, dxidx, jac)

  if is_shock
    fill!(Se_jac, 0); fill!(ee_jac, 0)
    # only compute the derivative if there is a shock
    getShockSensor(params, sbp, sensor, u, elnum, coords, dxidx, jac, Se, ee)
    ee_constant = getShockSensor_diff(params, sbp, sensor, u, elnum, coords,
                                              dxidx, jac, Se_jac, ee_jac)

    # the operator (for a scalar equation) is A = P^T * M * P * v, so
    # dR[p]/v[q] = (P^T * M * P)[p, q].  It then needs to be converted back
    # to conservative variables

    # compute derivative contribution from v
    @simd for p=1:numNodesPerElement
      @simd for q=1:numNodesPerElement
        getIRA0inv(params, sview(u, :, q), A0inv)
        # calculate the A[p, q]
        Apq = zero(Tres)
        @simd for k=1:numNodesPerElement
          Apq += capture.filt[k, p]*(sbp.w[k]/jac[k])*capture.filt[k, q]
        end
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            # this uses only the x direction viscoscity
            res_jac[i, j, p, q] = -ee[1, p]*Apq*A0inv[i, j]
          end
        end
      end
    end

    # compute derivative contribution from ee
    if !ee_constant

      @simd for i=1:numNodesPerElement
        w_i = sview(w, :, i)
        q_i = sview(u, :, i)
        convertToIR(params, q_i, w_i)
      end

      # apply P
      smallmatmatT!(w, capture.filt, t1)

      # apply mass matrix
      @simd for i=1:numNodesPerElement
        fac = sbp.w[i]/jac[i]
        @simd for j=1:numDofPerNode
          t1[j, i] *= fac
        end
      end

      # apply P^T
      smallmatmat!(t1, capture.filt, t2)

      @simd for p=1:numNodesPerElement
        @simd for q=1:numNodesPerElement
          @simd for j=1:numDofPerNode
            @simd for i=1:numDofPerNode
              res_jac[i, j, p, q] -= ee_jac[1, j, p, q]*t2[i, p]
            end
          end
        end
      end

    end  # if ee_constant

  end  # end if

  return is_shock
end


text:
module App

using Redux
using CImGui

include("../../todomvc/TodoMVC.jl")
using .TodoMVC
import .TodoMVC: Todo, TodoFilter

include("../../Renderer.jl")
using .Renderer

include("../../todomvc/gui.jl")

Base.@ccallable function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    store = create_store(TodoMVC.todo_mvc, TodoMVC.INITIAL_STATE)
    t = Renderer.render(()->todo_mvc_ui(store), width=600, height=300, title="App: TodoMVC")
    wait(t)
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end # module


text:
using CSV
using DataFrames

function loadData()
    startDate = "1/22/20"
    endDate = "3/30/20"
    roi = [(0,"Italy"), (0,"Spain"), ("Hubei","China"), (0,"Germany"), (0,"France"), (0,"Iran"), (0,"United Kingdom"), (0,"Switzerland"), (0,"Turkey"), (0,"Belgium"), (0,"Netherlands"), (0,"Austria"), (0,"Korea, South")] #(Province/State, Country/Region); Province/State=0 if data empty

    for i = 1:length(roi)
        buildDataFrame(roi[i],startDate,endDate)
    end

end

function buildDataFrame(roi,startDate,endDate)
    csvfile_confirmed = CSV.File("time_series_covid19_confirmed_global.csv")
    csvfile_deaths = CSV.File("time_series_covid19_deaths_global.csv")
    csvfile_recovered = CSV.File("time_series_covid19_recovered_global.csv")
    dfconfirmed = DataFrame(csvfile_confirmed)
    dfdeaths = DataFrame(csvfile_deaths)
    dfrecovered = DataFrame(csvfile_recovered)

    # region of interest
    if roi[1] == 0
        rowind_confirmed = .!(typeof.(dfconfirmed[:,Symbol("Province/State")]) .== String) .& (dfconfirmed[:,Symbol("Country/Region")] .== roi[2])
        rowind_deaths = .!(typeof.(dfdeaths[:,Symbol("Province/State")]) .== String) .& (dfdeaths[:,Symbol("Country/Region")] .== roi[2])
        rowind_recovered = .!(typeof.(dfrecovered[:,Symbol("Province/State")]) .== String) .& (dfrecovered[:,Symbol("Country/Region")] .== roi[2])
    else
        rowind_confirmed = (dfconfirmed[:,Symbol("Province/State")] .== roi[1]) .& (dfconfirmed[:,Symbol("Country/Region")] .== roi[2])
        rowind_deaths = (dfdeaths[:,Symbol("Province/State")] .== roi[1]) .& (dfdeaths[:,Symbol("Country/Region")] .== roi[2])
        rowind_recovered = (dfrecovered[:,Symbol("Province/State")] .== roi[1]) .& (dfrecovered[:,Symbol("Country/Region")] .== roi[2])
    end
    dfconfirmed_roi = dfconfirmed[rowind_confirmed,:];
    dfdeaths_roi = dfdeaths[rowind_deaths,:];
    dfrecovered_roi = dfrecovered[rowind_recovered,:];

    # time period of interest
    colind = collect(1:size(dfconfirmed_roi,2))
    startInd = colind[names(dfconfirmed_roi) .== Symbol(startDate)][1]
    endInd = colind[names(dfconfirmed_roi) .== Symbol(endDate)][1]

    # tailored to roi and toi
    nameInd = names(dfconfirmed_roi)[startInd:endInd]
    dfconfirmed_roi_toi = dfconfirmed_roi[startInd:endInd]
    dfdeaths_roi_toi = dfdeaths_roi[startInd:endInd]
    dfrecovered_roi_toi = dfrecovered_roi[startInd:endInd]

    # combine confirm, death, recover, active
    dflength = size(dfconfirmed_roi_toi,2)
    dfroi = DataFrame()
    dfroi.date = nameInd
    dfroi.cum_case = stack(dfconfirmed_roi_toi,1:dflength)[:,:value]
    dfroi.cum_death = stack(dfdeaths_roi_toi,1:dflength)[:,:value]
    dfroi.cum_recover = stack(dfrecovered_roi_toi,1:dflength)[:,:value]
    dfroi.active = dfroi.cum_case - dfroi.cum_death - dfroi.cum_recover

    CSV.write("covid_timeseries_" * roi[2] * ".csv",dfroi)

end

max_stars_repo_name: cannin/covid-sicr-test

max_issues_repo_name: cannin/covid-sicr-test

max_forks_repo_name: cannin/covid-sicr-test


text:
"""
Implement a variant of MGT using momentum accelaration.

*Experimental*
"""
module Momentum
import GraphIdx
import GraphIdx: Graph, Weights
import TreeLas: GapLas
import TreeLas.GapLas: gaplas, GapMem, gaplas!, Sol

include("broadcast.jl")

struct MomMem{N, WL<:Weights{Float64}}
    learn::Float64
    mass::Float64
    sol::Sol{N}
    sol1::Sol{N}        # solution, one step before
    sol2::Sol{N}        # solution, two steps before
    gmem::GapMem{N, WL}
    gamma::Vector{Float64}      # overload of `gmem.gamma` to be compliant in `gaplas!`
end

function MomMem(y::Array{Float64,N}, graph::Graph, lambda::Weights{Float64}) where {N}
    learn = 0.75
    mass = 0.95
    gmem = GapMem(y, graph, lambda)
    sol1 = similar(gmem.sol)
    sol1.x .= 0
    sol1.α .= 0
    sol2 = similar(gmem.sol)
    let η = learn, ζ = (1 - η) * mass, ι = (1 - η) * (1 - mass)
        @info "steps" η=round(η, digits=4) ζ=round(ζ, digits=4) ι=round(ι, digits=4)
    end
    MomMem(float(learn), float(mass), gmem.sol, sol1, sol2, gmem, gmem.gamma)
end

function gaplas!(
    mem::MomMem{N, W1},
    y::Array{Float64,N},
    graph::Graph,
    lambda::Weights{Float64},
    mu::Weights{Float64},
) where {N, W1<:Weights{Float64}}
    mem.sol2 .= mem.sol1
    mem.sol1 .= mem.sol
    GapLas.gaplas!(mem.gmem, y, graph, lambda, mu)
    let η = mem.learn, ζ = (1 - η) * mem.mass, ι = (1 - η) * (1 - mem.mass)
        @assert η + ζ + ι ≈ 1
        @. mem.sol.x = η * mem.sol.x + ζ * mem.sol1.x + ι * mem.sol2.x
        @. mem.sol.α = η * mem.sol.α + ζ * mem.sol1.α + ι * mem.sol2.α
    end
end

end # module Momentum

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
export TraceIdType,
    SpanIdType,
    INVALID_TRACE_ID,
    INVALID_SPAN_ID,
    current_span,
    TraceFlag,
    TraceState,
    SpanContext,
    INVALID_SPAN_CONTEXT,
    span_context,
    SPAN_KIND_UNSPECIFIED,
    SPAN_KIND_INTERNAL,
    SPAN_KIND_SERVER,
    SPAN_KIND_CLIENT,
    SPAN_KIND_PRODUCER,
    SPAN_KIND_CONSUMER,
    LimitInfo,
    Link,
    Event,
    SPAN_STATUS_UNSET,
    SPAN_STATUS_OK,
    SPAN_STATUS_ERROR,
    SpanStatus

const TraceIdType = UInt128
const SpanIdType = UInt64

const INVALID_TRACE_ID = zero(TraceIdType)
const INVALID_SPAN_ID = zero(SpanIdType)

const SPAN_KEY_IN_CONTEXT = create_key("current-span")

"""
    current_span([current_context])

Get the span in the current context.
"""
current_span() = current_span(current_context())
current_span(ctx::Context) = get(ctx, SPAN_KEY_IN_CONTEXT, nothing)

Base.@kwdef struct TraceFlag
    sampled::Bool = false
end

const MAXIMUM_TRACESTATE_KEYS = 32
const TRACESTATE_KEY_PATTERN =
    r"^[a-z][_0-9a-z\\-\\*\\/]{0,255}|[a-z0-9][_0-9a-z\\-\\*\\/]{0,240}@[a-z][_0-9a-z\\-\\*\\/]{0,13}$"
const TRACESTATE_VAL_PATTERN =
    r"^[\x20-\x2b\x2d-\x3c\x3e-\x7e]{0,255}[\x21-\x2b\x2d-\x3c\x3e-\x7e]$"

function is_valid_trace_state_kv(k, v)
    if isnothing(match(TRACESTATE_KEY_PATTERN, k))
        @warn "Invalid key $k found"
        false
    elseif isnothing(match(TRACESTATE_VAL_PATTERN, v))
        @warn "Invalid val $v found"
        false
    else
        true
    end
end

"""
    TraceState(entries::Pair{String,String}...)

`TraceState` carries vendor-specific trace identification data, represented as a list of key-value pairs. `TraceState` allows multiple tracing systems to participate in the same trace. It is fully described in the [W3C Trace Context specification](https://www.w3.org/TR/trace-context/#tracestate-header).
"""
struct TraceState{T<:NamedTuple}
    kv::T
    function TraceState(entries::Pair{String,String}...)
        if length(entries) > MAXIMUM_TRACESTATE_KEYS
            @warn "There can't be more than $MAXIMUM_TRACESTATE_KEYS key/value pairs"
        else
            kv = NamedTuple(
                Symbol(k) => v for (k, v) in entries if is_valid_trace_state_kv(k, v)
            )
            new{typeof(kv)}(kv)
        end
    end
end

Base.getindex(s::TraceState, key) = s.kv[key]
Base.haskey(s::TraceState, key) = haskey(s.kv, key)
Base.haskey(s::TraceState, key::String) = haskey(s, Symbol(key))
Base.length(s::TraceState) = length(s.kv)

function Base.show(io::IO, ts::TraceState)
    for (i, (k, v)) in enumerate(pairs(ts.kv))
        print(io, k)
        print(io, '=')
        print(io, v)
        if i != length(ts)
            print(io, ',')
        end
    end
end

"""
    SpanContext(;span_id, trace_id, is_remote, trace_flag=TraceFlag(), trace_state=TraceState())

A `SpanContext` represents the portion of a `Span` which must be serialized and
propagated along side of a distributed context. `SpanContext`s are immutable.

The OpenTelemetry `SpanContext` representation conforms to the [W3C TraceContext
specification](https://www.w3.org/TR/trace-context/). It contains two
identifiers - a `TraceId` and a `SpanId` - along with a set of common
`TraceFlags` and system-specific `TraceState` values.

`TraceId` A valid trace identifier is a 16-byte array with at least one
non-zero byte.

`SpanId` A valid span identifier is an 8-byte array with at least one non-zero
byte.

`TraceFlags` contain details about the trace. Unlike TraceState values,
TraceFlags are present in all traces. The current version of the specification
only supports a single flag called [sampled](https://www.w3.org/TR/trace-context/#sampled-flag).

`TraceState` carries vendor-specific trace identification data, represented as a list
of key-value pairs. TraceState allows multiple tracing
systems to participate in the same trace. It is fully described in the [W3C Trace Context
specification](https://www.w3.org/TR/trace-context/#tracestate-header).
"""
Base.@kwdef struct SpanContext{TS<:TraceState}
    span_id::SpanIdType
    trace_id::TraceIdType
    is_remote::Bool
    trace_flag::TraceFlag = TraceFlag()
    trace_state::TS = TraceState()
end

const INVALID_SPAN_CONTEXT =
    SpanContext(; trace_id = INVALID_TRACE_ID, span_id = INVALID_SPAN_ID, is_remote = false)

"""
    span_context([s::Span])

Get the `SpanContext` from a span `s`. If `s` is not specified, [`current_span()`](@ref) will be used. `nothing` is returned if no span context found.
"""
span_context(s) = s.span_context
span_context(::Nothing) = nothing
span_context() = span_context(current_span())

# !!! must be of the same order with https://github.com/open-telemetry/opentelemetry-proto/blob/main/opentelemetry/proto/trace/v1/trace.proto
@enum SpanKind begin
    SPAN_KIND_UNSPECIFIED
    SPAN_KIND_INTERNAL
    SPAN_KIND_SERVER
    SPAN_KIND_CLIENT
    SPAN_KIND_PRODUCER
    SPAN_KIND_CONSUMER
end

"""
    LimitInfo(;kw...)

Used in [`TracerProvider`](@ref) to configure generated [`Tracer`](@ref).

# Keyword arguments:

  - `span_attribute_count_limit::Int` = 128
  - `span_attribute_value_length_limit::Union{Nothing,Int}` = nothing
  - `span_event_count_limit::Int` = 128
  - `span_link_count_limit::Int` = 128
"""
Base.@kwdef struct LimitInfo
    span_attribute_count_limit::Int = 128
    span_attribute_value_length_limit::Union{Nothing,Int} = nothing
    span_event_count_limit::Int = 128
    span_link_count_limit::Int = 128
end

"""
    Link(span_context, attributes)

See more details at [links between spans](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/overview.md#links-between-spans).
"""
struct Link{C<:SpanContext,A<:StaticAttrs}
    context::C
    attributes::A
end

"""
    Event(name, timestamp=time()*10^9, attributes=StaticAttrs())

`timestamp` is the *nanoseconds*.
"""
Base.@kwdef struct Event{A<:StaticAttrs}
    name::String
    timestamp::UInt = UInt(time() * 10^9)
    attributes::A = StaticAttrs()
end

# !!! must be of the same order as https://github.com/open-telemetry/opentelemetry-proto/blob/main/opentelemetry/proto/trace/v1/trace.proto
@enum SpanStatusCode begin
    SPAN_STATUS_UNSET
    SPAN_STATUS_OK
    SPAN_STATUS_ERROR
end

"""
    SpanStatus(code, description=nothing)

Possible codes are:

  - `SPAN_STATUS_UNSET`
  - `SPAN_STATUS_ERROR`
  - `SPAN_STATUS_OK`

`description` is required when `code` is `SPAN_STATUS_ERROR`.
"""
struct SpanStatus
    code::SpanStatusCode
    description::Union{String,Nothing}
    function SpanStatus(code::SpanStatusCode, description = nothing)
        if code === SPAN_STATUS_ERROR
            isnothing(description) && @error "description not provided"
            new(code, description)
        else
            # description is ignored
            new(code, nothing)
        end
    end
end


text:
# # Zombie Outbreak
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/videos/agents/zombies.mp4?raw=true" type="video/mp4">
# </video>
# ```
#
# This model showcases an ABM running on a map, using [`OpenStreetMapSpace`](@ref).
#
# ## Constructing the end of days
using Agents

# We'll simulate a zombie outbreak in a city. To do so, we start with an agent which
# satisfies the OSMSpace conditions of having a `pos`ition of type
# `Tuple{Int,Int,Float64}`, a `route` vector and a `destination` with the same type as
# `pos`. For simplicity though we shall build this with the [`@agent`](@ref) macro.

@agent Zombie OSMAgent begin
    infected::Bool
end

# To be explicit, this macro builds the following type:
# ```julia
# mutable struct Zombie <: AbstractAgent
#     id::Int
#     pos::Tuple{Int,Int,Float64}
#     route::Vector{Int}
#     destination::Tuple{Int,Int,Float64}
#     infected::Bool
# end
# ```
# where a tuple `(i, j, x)::Tuple{Int,Int,Float64}` means a position
# on the road between nodes `i, j` of the map, having progressed `x` meters along the road.

# The model constructor we build consists of a map, and 100 agents scattered randomly
# around it. They have their own agenda and need to travel to some new destination.
# Unfortunately one of the population has turned and will begin infecting anyone who
# comes close.

function initialise(; map_path = OSM.TEST_MAP)
    model = ABM(Zombie, OpenStreetMapSpace(map_path))

    for id in 1:100
        start = random_position(model) # At an intersection
        finish = OSM.random_road_position(model) # Somewhere on a road
        route = OSM.plan_route(start, finish, model)
        human = Zombie(id, start, route, finish, false)
        add_agent_pos!(human, model)
    end
    ## We'll add patient zero at a specific (latitude, longitude)
    start = OSM.road((39.52320181536525, -119.78917553184259), model)
    finish = OSM.intersection((39.510773, -119.75916700000002), model)
    route = OSM.plan_route(start, finish, model)
    ## This function call creates & adds an agent, see `add_agent!`
    zombie = add_agent!(start, model, route, finish, true)
    return model
end

# In our model, zombies are seemingly oblivious to their state, since they keep going about their
# business, but start eating people along the way. Perhaps they can finally express their distaste
# for city commuting.

function agent_step!(agent, model)
    ## Each agent will progress 25 meters along their route
    move_along_route!(agent, model, 25)

    if is_stationary(agent, model) && rand(model.rng) < 0.1
        ## When stationary, give the agent a 10% chance of going somewhere else
        OSM.random_route!(agent, model)
        ## Start on new route
        move_along_route!(agent, model, 25)
    end

    if agent.infected
        ## Agents will be infected if they get within 50 meters of a zombie.
        map(i -> model[i].infected = true, nearby_ids(agent, model, 50))
    end
end

# ## Visualising the fall of humanity
#
# Plotting this space in a seamless manner is a work in progress. For now we
# use [OpenStreetMapXPlot](https://github.com/pszufe/OpenStreetMapXPlot.jl) and
# a custom routine.

# ```julia
# using OpenStreetMapXPlot
# using Plots
# gr()
# ```

ac(agent) = agent.infected ? :green : :black
as(agent) = agent.infected ? 6 : 5

function plotagents(model)
    ids = model.scheduler(model)
    colors = [ac(model[i]) for i in ids]
    sizes = [as(model[i]) for i in ids]
    markers = :circle
    pos = [OSM.map_coordinates(model[i], model) for i in ids]

    scatter!(
        pos;
        markercolor = colors,
        markersize = sizes,
        markershapes = markers,
        label = "",
        markerstrokewidth = 0.5,
        markerstrokecolor = :black,
        markeralpha = 0.7,
    )
end

# Let's see how this plays out!
# ```julia
# model = initialise()
#
# frames = @animate for i in 0:200
#     i > 0 && step!(model, agent_step!, 1)
#     plotmap(model.space.m)
#     plotagents(model)
# end
#
# gif(frames, "outbreak.gif", fps = 15)
# ```
#
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/videos/agents/zombies.mp4?raw=true" type="video/mp4">
# </video>
# ```


text:
### View types

# use ContiguousView when contiguousness can be determined statically
immutable ContiguousView{T,N,Arr<:Array} <: ArrayView{T,N,N}
    arr::Arr
    offset::Int
    len::Int
    shp::NTuple{N,Int}
end

ContiguousView{T,N}(arr::Array{T}, offset::Int, shp::NTuple{N,Int}) =
    ContiguousView{T,N,typeof(arr)}(arr, offset, prod(shp), shp)

ContiguousView(arr::Array, shp::Dims) = ContiguousView(arr, 0, shp)

immutable UnsafeContiguousView{T,N} <: UnsafeArrayView{T,N,N}
    ptr::Ptr{T}
    len::Int
    shp::NTuple{N,Int}
end

UnsafeContiguousView{T,N}(ptr::Ptr{T}, shp::NTuple{N,Int}) =
    UnsafeContiguousView{T,N}(ptr, prod(shp), shp)

UnsafeContiguousView{T,N}(ptr::Ptr{T}, offset::Int, shp::NTuple{N,Int}) =
    UnsafeContiguousView(ptr+offset*sizeof(T), shp)

UnsafeContiguousView(arr::Array, shp::Dims) = UnsafeContiguousView(pointer(arr), shp)

UnsafeContiguousView{T,N}(arr::Array{T}, offset::Int, shp::NTuple{N,Int}) =
    UnsafeContiguousView(pointer(arr, offset+1), shp)

# use StridedView when contiguousness can not be determined statically
# condition: M < N
immutable StridedView{T,N,M,Arr<:Array} <: ArrayView{T,N,M}
    arr::Arr
    offset::Int
    len::Int
    shp::NTuple{N,Int}
    strides::NTuple{N,Int}
end

function StridedView{T,N,M}(arr::Array{T}, offset::Int, shp::NTuple{N,Int},
                             ::Type{ContRank{M}}, strides::NTuple{N,Int})
    @assert M < N
    StridedView{T,N,M,typeof(arr)}(arr, offset, prod(shp), shp, strides)
end

function StridedView{T,N,M}(arr::Array{T}, shp::NTuple{N,Int},
                             ::Type{ContRank{M}}, strides::NTuple{N,Int})
    @assert M < N
    StridedView{T,N,M,typeof(arr)}(arr, 0, prod(shp), shp, strides)
end

immutable UnsafeStridedView{T,N,M} <: UnsafeArrayView{T,N,M}
    ptr::Ptr{T}
    len::Int
    shp::NTuple{N,Int}
    strides::NTuple{N,Int}
end

function UnsafeStridedView{T,N,M}(ptr::Ptr{T}, shp::NTuple{N,Int},
                                  ::Type{ContRank{M}}, strides::NTuple{N,Int})
    @assert M < N
    UnsafeStridedView{T,N,M}(ptr, prod(shp), shp, strides)
end

function UnsafeStridedView{T,N,M}(ptr::Ptr{T}, offset::Int, shp::NTuple{N,Int},
                                  ::Type{ContRank{M}}, strides::NTuple{N,Int})
    @assert M < N
    UnsafeStridedView(ptr+offset*sizeof(T), shp, ContRank{M}, strides)
end

function UnsafeStridedView{T,N,M}(arr::Array{T}, offset::Int, shp::NTuple{N,Int},
                                  ::Type{ContRank{M}}, strides::NTuple{N,Int})
    @assert M < N
    UnsafeStridedView(pointer(arr, offset+1), shp, ContRank{M}, strides)
end

function UnsafeStridedView{T,N,M}(arr::Array{T}, shp::NTuple{N,Int},
                                  ::Type{ContRank{M}}, strides::NTuple{N,Int})
    @assert M < N
    UnsafeStridedView(pointer(arr), shp, ContRank{M}, strides)
end

### basic methods

parent(a::ArrayView) = a.arr
parent(a::UnsafeArrayView) = error("Getting parent of an unsafe view is not allowed.")
parent_or_ptr(a) = parent(a)
parent_or_ptr(a::UnsafeArrayView) = a.ptr

uget(a::ArrayView, i::Int) = getindex(a.arr, a.offset + i)
uset!{T}(a::ArrayView{T}, v::T, i::Int) = setindex!(a.arr, v, a.offset + i)

uget(a::UnsafeArrayView, i::Int) = unsafe_load(a.ptr, i)
uset!{T}(a::UnsafeArrayView{T}, v::T, i::Int) = unsafe_store!(a.ptr, v, i)

offset(a::ArrayView) = a.offset
pointer(a::ArrayView) = pointer(parent(a), a.offset+1)
pointer(a::UnsafeArrayView) = a.ptr


text:
struct None <: AbstractProx end

function prox!(::None,step::Real,xprev::AbstractVector,gcurr::AbstractVector,xcurr::AbstractVector)
    xcurr .= xprev - step*gcurr
end


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 13, n}) where n
    io = openInput(puzzle)
    split(read(io, String), "\n", keepempty=false)
end

function solve(::Puzzle{2020, 13, 1}, input)
    nothing
end

function solve(::Puzzle{2020, 13, 2}, input)
    nothing
end

@testset "Day 13" begin
    TEST_INPUT_1 = """
"""

    @test_skip(solve(Puzzle(13, 1), TEST_INPUT_1) == nothing)
    @test_skip(solve(Puzzle(13, 1), input(Puzzle(13, 1))) == nothing)

    TEST_INPUT_2 = """
"""

    @test_skip(solve(Puzzle(13, 2), TEST_INPUT_2) == nothing)
    @test_skip(solve(Puzzle(13, 2), input(Puzzle(13, 2))) == nothing)
end


text:
## ExoplanetsSysSim/src/mission_constants.jl
## (c) 2015 Eric B. Ford

#module KeplerConstants
 const num_channels = 84
 const num_modules = 42
 const num_quarters = 17              # QUERY:  I'm favoring leaving out quarter 0, since that was engineering data.  Agree?
 const cdpp_durations = [1.5,2.,2.5,3.,3.5,4.5,5.,6.,7.5,9.,10.5,12.,12.5,15.]
 const duration_symbols = [:rrmscdpp01p5, :rrmscdpp02p0,:rrmscdpp02p5,:rrmscdpp03p0,:rrmscdpp03p5,:rrmscdpp04p5,:rrmscdpp05p0,:rrmscdpp06p0,:rrmscdpp07p5,:rrmscdpp09p0,:rrmscdpp10p5,:rrmscdpp12p0,:rrmscdpp12p5,:rrmscdpp15p0 ]
 const num_cdpp_timescales = 14
 @assert num_cdpp_timescales == length(cdpp_durations) == length(duration_symbols)
 const mission_data_span = 1459.789   # maximum(ExoplanetsSysSim.StellarTable.df[:dataspan])
 const mission_duty_cycle = 0.8751    # median(ExoplanetsSysSim.StellarTable.df[:dutycycle])

 const kepler_exp_time_internal  =  6.019802903/(24*60*60)    # https://archive.stsci.edu/kepler/manuals/archive_manual.pdf
 const kepler_read_time_internal = 0.5189485261/(24*60*60)    # https://archive.stsci.edu/kepler/manuals/archive_manual.pdf
 const num_exposures_per_LC = 270
 const num_exposures_per_SC = 9
 const LC_integration_time = kepler_exp_time_internal*num_exposures_per_LC
 const SC_integration_time = kepler_exp_time_internal*num_exposures_per_SC
 const LC_read_time = kepler_read_time_internal*num_exposures_per_LC
 const SC_read_time = kepler_read_time_internal*num_exposures_per_SC
 const LC_duration = LC_integration_time +  LC_read_time
 const SC_duration = SC_integration_time +  SC_read_time
 const LC_rate = 1.0/LC_duration
 const SC_rate = 1.0/SC_duration

# export num_channels,num_modules,num_quarters,num_cdpp_timescales,mission_data_span,mission_duty_cycle
# export kepler_exp_time_internal,kepler_read_time_internal,num_exposures_per_LC,num_exposures_per_SC
# export LC_integration_time,SC_integration_time,LC_read_time,SC_read_time,LC_duration,SC_duration,LC_rate,SC_rate
#end

#using KeplerConstants

# Standard conversion factors on which unit system is based
const global AU_in_m_IAU2012 = 149597870700.0
 const global G_in_mks_IAU2015 = 6.67384e-11
 const global G_mass_sun_in_mks = 1.3271244e20
 const global G_mass_earth_in_mks = 3.986004e14
 const global sun_radius_in_m_IAU2015 = 6.9566e8
 const global earth_radius_eq_in_m_IAU2015 = 6.3781e6
 const global sun_mass_in_kg_IAU2010 = 1.988547e30

# Constants used by this code
const global sun_mass = 1.0
 const global earth_mass = G_mass_earth_in_mks/G_mass_sun_in_mks  # about 3.0024584e-6
 const global earth_radius = earth_radius_eq_in_m_IAU2015 / sun_radius_in_m_IAU2015 # about 0.0091705248
 const global rsol_in_au = sun_radius_in_m_IAU2015 / AU_in_m_IAU2012  # about 0.00464913034
 const global sec_in_day = 24*60*60
 const global grav_const = G_in_mks_IAU2015 * sec_in_day^2 * sun_mass_in_kg_IAU2010 / AU_in_m_IAU2012^3 # about 2.9591220363e-4 in AU^3/(day^2 Msol)
 const global day_in_year = 365.25


text:
"""
    GeneralCache

Supertype for all kinds of cache types that are needed throughout the simulation.
"""
abstract type GeneralCache end

"""
    CellCache{C,A} <: GeneralCache

A special type of cache to handle `CellList` and auxiliary object.

# Fields
- `cell::C`: stores the cell list object, normally of type `CellListMap.CellList`.
- `aux::A`: stores the auxiliary lists for cell lists, useful when having to continously
    update the `cell` field.
"""
mutable struct CellCache{A} <: GeneralCache
    cell::CellListMap.CellList
    aux::A
end

function _build_cache(cell::CellListMap.CellList)
    aux = CellListMap.AuxThreaded(cell)

    return CellCache(cell, aux)
end


text:
# This contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

import Base.Order: Forward, Ordering, lt

# Heap operations on flat arrays
# ------------------------------

# Binary heap indexing
heapleft(i::Integer) = 2i
heapright(i::Integer) = 2i + 1
heapparent(i::Integer) = div(i, 2)

# Binary min-heap percolate down.
function percolate_down!(xs::AbstractArray, i::Integer, x=xs[i], o::Ordering=Forward, len::Integer=length(xs))
    @inbounds while (l = heapleft(i)) <= len
        r = heapright(i)
        j = r > len || lt(o, xs[l], xs[r]) ? l : r
        if lt(o, xs[j], x)
            xs[i] = xs[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
end

percolate_down!(xs::AbstractArray, i::Integer, o::Ordering, len::Integer=length(xs)) = percolate_down!(xs, i, xs[i], o, len)

# Binary min-heap percolate up.
function percolate_up!(xs::AbstractArray, i::Integer, x=xs[i], o::Ordering=Forward)
    @inbounds while (j = heapparent(i)) >= 1
        if lt(o, x, xs[j])
            xs[i] = xs[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
end

percolate_up!(xs::AbstractArray{T}, i::Integer, o::Ordering) where {T} = percolate_up!(xs, i, xs[i], o)

"""
    heappop!(v, [ord])

Given a binary heap-ordered array, remove and return the lowest ordered element.
For efficiency, this function does not check that the array is indeed heap-ordered.
"""
function heappop!(xs::AbstractArray, o::Ordering=Forward)
    x = xs[1]
    y = pop!(xs)
    if !isempty(xs)
        percolate_down!(xs, 1, y, o)
    end
    x
end

"""
    heappush!(v, x, [ord])

Given a binary heap-ordered array, push a new element `x`, preserving the heap property.
For efficiency, this function does not check that the array is indeed heap-ordered.
"""
function heappush!(xs::AbstractArray, x, o::Ordering=Forward)
    push!(xs, x)
    percolate_up!(xs, length(xs), x, o)
    xs
end

# Turn an arbitrary array into a binary min-heap in linear time.
"""
    heapify!(v, ord::Ordering=Forward)

In-place [`heapify`](@ref).
"""
function heapify!(xs::AbstractArray, o::Ordering=Forward)
    for i in heapparent(length(xs)):-1:1
        percolate_down!(xs, i, o)
    end
    xs
end

"""
    heapify(v, ord::Ordering=Forward)

Returns a new vector in binary heap order, optionally using the given ordering.
```jldoctest
julia> a = [1,3,4,5,2];

julia> heapify(a)
5-element Array{Int64,1}:
 1
 2
 4
 5
 3

julia> heapify(a, Base.Order.Reverse)
5-element Array{Int64,1}:
 5
 3
 4
 1
 2
```
"""
heapify(xs::AbstractArray, o::Ordering=Forward) = heapify!(copyto!(similar(xs), xs), o)

"""
    isheap(v, ord::Ordering=Forward)

Return `true` if an array is heap-ordered according to the given order.

```jldoctest
julia> a = [1,2,3]
3-element Array{Int64,1}:
 1
 2
 3

julia> isheap(a,Base.Order.Forward)
true

julia> isheap(a,Base.Order.Reverse)
false
```
"""
function isheap(xs::AbstractArray, o::Ordering=Forward)
    for i in 1:div(length(xs), 2)
        if lt(o, xs[heapleft(i)], xs[i]) ||
           (heapright(i) <= length(xs) && lt(o, xs[heapright(i)], xs[i]))
            return false
        end
    end
    true
end


text:
pop(α,i) = [α[1:i-1]...,α[i+1:end]...] # pop out ith entry
pad(α,i) = [α[1:i-1]...,zero(eltype(α)),α[i:end]...] # pad by 0

"
function recursive_nodes(α,node_family):
Note α is a barycentric coordinate here, not a multi-index.
"
function recursive_nodes(α,node_family)
    d = length(α)-1
    xn = node_family[sum(α)+1]
    if d==1
        b = [xn[α[1]+1],xn[α[2]+1]]
        return b
    else
        b = zeros(eltype(xn),d+1)
        weight = zero(eltype(xn))
        for i = 1:d+1
            br = pad(recursive_nodes(pop(α,i),node_family),i)
            x_i = xn[sum(pop(α,i))+1]
            @. b += x_i*br
            weight += x_i
        end
        return b./weight
    end
end

# recursive_nodes(d::Int,N::Int,node_family) =
#     recursive_nodes(MultiIndexSet(d,N),node_family)
#
# function recursive_nodes(set::MultiIndexSet,node_family)
#     X = ntuple(x->zeros(length(Ndset)),d)
#     for (index,α) in enumerate(Ndset)
#         b = recursive_nodes((α...,N-sum(α)),node_family)
#         for i = 1:d
#             X[i][index] = sum(b.*VXYZ[i])
#         end
#     end
#     return X
# end