text:
abstract type DeterministicSolution{dType, tType, N} <: Solution{dType, tType, N} end

"""
```julia
Solution(equation, Δt, ntime; kwargs...)
```

Create the appropriate `Solution` for the given `equation` type for a
simulation with `ntime` time steps of step size `Δt`.

"""
function Solution end

# Create solution for ODE.
function Solution(equation::AbstractEquationODE, Δt, ntime::Int; kwargs...)
    SSolutionODE(equation, Δt, ntime; kwargs...)
end

# Create solution for partitioned ODE.
function Solution(equation::AbstractEquationPODE, Δt, ntime::Int; kwargs...)
    SSolutionPODE(equation, Δt, ntime; kwargs...)
end

# Create solution for DAE.
function Solution(equation::AbstractEquationDAE, Δt, ntime::Int; kwargs...)
    SSolutionDAE(equation, Δt, ntime; kwargs...)
end

# Create solution for partitioned DAE.
function Solution(equation::AbstractEquationPDAE, Δt, ntime::Int; kwargs...)
    SSolutionPDAE(equation, Δt, ntime; kwargs...)
end

# Print error for solutions of equations not implemented, yet.
function Solution(equation::Equation, Δt, ntime::Int; kwargs...)
    error("No solution found for equation ", equation)
end

"""
```julia
ParallelSolution(equation, Δt, ntime; kwargs...)
```

Create the appropriate `ParallelSolution` for the given `equation` type for a
simulation with `ntime` time steps of step size `Δt`.

"""
function ParallelSolution end

# Create parallel solution for ODE.
function ParallelSolution(equation::AbstractEquationODE, Δt, ntime::Int; kwargs...)
    PSolutionODE(equation, Δt, ntime; kwargs...)
end

# Create parallel solution for partitioned ODE.
function ParallelSolution(equation::AbstractEquationPODE, Δt, ntime::Int; kwargs...)
    PSolutionPODE(equation, Δt, ntime; kwargs...)
end

# Create parallel solution for DAE.
function ParallelSolution(equation::AbstractEquationDAE, Δt, ntime::Int; kwargs...)
    PSolutionDAE(equation, Δt, ntime; kwargs...)
end

# Create parallel solution for partitioned DAE.
function ParallelSolution(equation::AbstractEquationPDAE, Δt, ntime::Int; kwargs...)
    PSolutionPDAE(equation, Δt, ntime; kwargs...)
end

# Print error for parallel solutions of equations not implemented, yet.
function ParallelSolution(equation::Equation, Δt, ntime::Int; kwargs...)
    error("No parallel solution found for equation ", equation)
end


text:
mutable struct Pressures
    pressures::Vector{Float}

    Pressures(
        mass_component::MassComponent,
        energy_component::EnergyComponent,
        eos::EquationOfState,
    ) = new([
        gas_pressure_from(eos, ρ, u)
        for
        (ρ, u) in zip(
            get_mass_densities(mass_component),
            get_specific_energies(energy_component),
        )
    ])
end

"""
    get_pressures(pressures::Pressures)
        -> Vector{Float}

Return the pressures of the particles.
"""
get_pressures(pressures::Pressures) = pressures.pressures

"""
    get_pressures(pressures::Pressures, (i, j)::Tuple{Unsigned,Unsigned})
        -> (Float, Float)

Return the pressures of the particles with the given indices.
"""
function get_pressures(pressures::Pressures, (i, j)::Tuple{Unsigned,Unsigned})
    pressures.pressures[i], pressures.pressures[j]
end

"""
    get_pressure(pressures::Pressures, i::Unsigned) -> (Float, Float)

Return the pressure of the particle with the given index.
"""
function get_pressure(pressures::Pressures, i::Unsigned)
    pressures.pressures[i]
end

"""
    updatepressures!(
        pressures::Pressures,
        mass_component::MassComponent,
        energy_component::EnergyComponent,
        eos::EquationOfState,
    )

Update the state of the pressure component according to the given mass
and energy component.
"""
updatepressures!(
    pressures::Pressures,
    mass_component::MassComponent,
    energy_component::EnergyComponent,
    eos::EquationOfState,
) = map!(
    (ρ, u) -> gas_pressure_from(eos, ρ, u),
    pressures.pressures,
    get_mass_densities(mass_component),
    get_specific_energies(energy_component),
)

"""
    compute_Δa_pressure(
        (ρᵢ, ρⱼ)::Tuple{Number,Number},
        (Pᵢ, Pⱼ)::Tuple{Number,Number},
        m_ddr_W_r_hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
    ) -> (Float, Float)

Compute pairwise particle accelerations due to pressure.
"""
function compute_Δa_pressure(
    (ρᵢ, ρⱼ)::Tuple{Number,Number},
    (Pᵢ, Pⱼ)::Tuple{Number,Number},
    m_ddr_W_r_hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
)
    compute_Δa_pressure(m_ddr_W_r_hᵢ_hⱼ, (Pᵢ / ρᵢ^2, Pⱼ / ρⱼ^2))
end

"""
    compute_Δa_pressure(
        m_ddr_W_r_h::Number,
        (Pᵢ_over_ρᵢ², Pⱼ_over_ρⱼ²)::Tuple{Number,Number},
    ) -> (Float, Float)

Compute pairwise particle accelerations due to pressure.
"""
function compute_Δa_pressure(
    m_ddr_W_r_h::Number,
    (Pᵢ_over_ρᵢ², Pⱼ_over_ρⱼ²)::Tuple{Number,Number},
)
    m_ddr_W_r_h * Pᵢ_over_ρᵢ², m_ddr_W_r_h * Pⱼ_over_ρⱼ²
end

"""
    compute_Δa_pressure(
        (m_ddr_W_r_hᵢ, m_ddr_W_r_hⱼ)::Tuple{Number,Number},
        (Pᵢ_over_ρᵢ², Pⱼ_over_ρⱼ²)::Tuple{Number,Number},
    ) -> (Float, Float)

Compute pairwise particle accelerations due to pressure.
"""
function compute_Δa_pressure(
    (m_ddr_W_r_hᵢ, m_ddr_W_r_hⱼ)::Tuple{Number,Number},
    (Pᵢ_over_ρᵢ², Pⱼ_over_ρⱼ²)::Tuple{Number,Number},
)
    m_ddr_W_r_hᵢ * Pᵢ_over_ρᵢ², m_ddr_W_r_hⱼ * Pⱼ_over_ρⱼ²
end

ext: jl

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# BUILD OPTIONS
const DEFAULT_SAVE_BUILD = false
const DEFAULT_OVERWRITE = false
const DEFAULT_VERSION = "1.0"
const DEFAULT_DATASET = "state_model"
const BUILD_STEPS = ["partition", "share", "calibrate", "disagg"];
const PARAM_DIR = "parameters"
const SET_DIR = "sets"

# CALIBRATION
const DEFAULT_PENALTY_NOKEY = 1e4
const DEFAULT_CALIBRATE_LOWER_BOUND = 0.1
const DEFAULT_CALIBRATE_UPPER_BOUND = 5

# NUMERICAL
const DEFAULT_TOL = 1e-6
const DEFAULT_SMALL = missing
const DEFAULT_ROUND_DIGITS = 10

# STATE MODEL
const SUB_ELAST = Dict()
SUB_ELAST[:va] = 1      # value-added nest
SUB_ELAST[:y] = 0       # top-level Y nest (VA,M)
SUB_ELAST[:m] = 0       # materials nest
SUB_ELAST[:a] = 0       # top-level A nest for aggregate demand (margins, goods)
SUB_ELAST[:mar] = 0     # margin supply
SUB_ELAST[:d] = 2       # domestic demand aggregation nest (intranational)
SUB_ELAST[:f] = 4       # domestic and foreign demand aggregation nest (international)

# Transportation elasticity: disposition, distribute regional supply to local, national, export
const TRANS_ELAST = Dict(:x => 4)

# Variable lower bound
const MODEL_LOWER_BOUND = 0.00

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
export link!

function link!(dst::Module, src::Module)
    status = convert(Core.Bool, API.LLVMLinkModules2(dst, src))
    @assert !status # caught by diagnostics handler

    return nothing
end


text:
"""
    @voiceargs args function

Declare some or all arguments of the `function` definition to be arguments that can be obtained by voice.

# Arguments
- `args`: a voicearg, a pair `voicearg=>kwargs` or a tuple with multiple voiceargs, which can each have kwargs or not (see examples below).
- `function`: the function definition.
!!! note "Keyword arguments definable for each voice argument in `args`"
    - `model::String=DEFAULT_MODEL_NAME`: the name of the model to be used for the function argument (the name must be one of the keys of the modeldirs dictionary passed to `init_jsi`).
    - `valid_input::AbstractArray{String}`: the valid speech input (e.g. `["up", "down"]`).
    - `valid_input_auto::Bool`: whether the valid speech input can automatically be derived from the type of the function argument.
    - `interpret_function::Function`: a function to interpret the token (mapping a String to a different String).
    - `use_max_speed::Bool=false`: whether to use maxium speed for the recognition of the next token (rather than maximum accuracy). It is generally only recommended to set `use_max_speed=true` for single word commands or very specfic use cases that require immediate minimal latency action when a command is said.
    - `ignore_unknown::Bool=false`: whether to ignore unknown tokens in the speech (and consume the next instead). It is generally not recommended to set `ignore_unkown=true` - in particular not in combination with a limited valid input - as then the function will block until it receives a token it recognizes.
    - `vararg_end::String`: a token to signal the end of a vararg (only valid if the function argument is a vararg).
    - `vararg_max::Integer=∞`: the maximum number of arguments the vararg can contain (only valid if the function argument is a vararg).
    - `vararg_timeout::AbstractFloat`: timeout after which to abort waiting for a next token to be spoken (only valid if the function argument is a vararg).

# Examples
```jldoctest
@voiceargs (b, c) function f(a, b::String, c::String, d)
    #(...)
    return
end

@voiceargs (b=>(use_max_speed=true), c) function f(a, b::String, c::String, d)
    #(...)
    return
end

@enum TypeMode words formula
@voiceargs (mode=>(valid_input_auto=true), token=>(model=TYPE_MODEL_NAME, vararg_timeout=2.0)) function type_tokens(mode::TypeMode, tokens::String...)
    #(...)
    return
end
```
"""
macro voiceargs(args...) checkargs(args...); handle_voiceargs(__module__, args...); end

## CONSTANTS

const USE_PARTIAL_RECOGNITIONS_DEFAULT = false
const USE_IGNORE_UNKNOWN_DEFAULT       = false
const VARARG_TIMEOUT_DEFAULT           = 60.0

## ARGUMENT CHECKS

function checkargs(args...)
    if (length(args) != 2) @ArgumentError("wrong number of arguments.") end
    if !(is_symbol(args[1]) || is_pair(args[1]) || is_tuple(args[1])) @ArgumentError("the first argument must be a voicearg, a pair `voicearg=>kwargs` or a tuple with multiple voiceargs, which can each have kwargs or not (obtained: $(args[1])).") end
    if !is_function(args[end]) @ArgumentError("the last argument must be a function definition (obtained: $(args[end])).") end
end

## VOICEARGS FUNCTIONS
function handle_voiceargs(caller::Module, voiceargs_arg::Union{Symbol,Expr}, f_expr::Expr)
    voiceargs_vect = extract_voiceargs(voiceargs_arg)
    simpleargs, complexargs_expr = split_args(voiceargs_vect)
    complexargs = split_complexargs(complexargs_expr)
    voiceargs = create_voiceargs_dict(simpleargs, complexargs)
    f_name, f_args = parse_signature(f_expr)
    eval_kwargs!(caller, voiceargs)
    validate_voiceargs(voiceargs, f_args)
    handle_valid_inputs!(caller, voiceargs, f_args)
    f_wrapper = wrap_f(f_name, f_args, f_expr, voiceargs)
    set_voiceargs(f_name, voiceargs)
    return quote
        $f_wrapper
        $(esc(f_expr))
    end
end

let
    global voiceargs, voicearg_f_names, recognizer, set_voiceargs, set_recognizer
    _f_voiceargs::Dict{Symbol, Dict{Symbol, Dict{Symbol, Any}}}                            = Dict{Symbol, Dict{Symbol, Dict{Symbol, Any}}}()
    voiceargs(f_name::Symbol)::Dict{Symbol, Dict{Symbol, Any}}                             = _f_voiceargs[f_name]
    voicearg_f_names()::Base.KeySet{Symbol, Dict{Symbol, Dict{Symbol, Dict{Symbol, Any}}}} = keys(_f_voiceargs)
    recognizer(f_name::Symbol, voicearg::Symbol)::PyObject                                 = _f_voiceargs[f_name][voicearg][:recognizer]
    set_voiceargs(f_name::Symbol, v::Dict{Symbol, Dict{Symbol, Any}})                      = (_f_voiceargs[f_name] = v; return)
    set_recognizer(f_name::Symbol, voicearg::Symbol, r::PyObject)                          = (_f_voiceargs[f_name][voicearg][:recognizer] = r; return)
end

## FUNCTIONS TO DEAL WITH VOICEARGS PARSING

is_function(arg)                 = isdef(arg)
is_symbol(arg)                   = isa(arg, Symbol)
is_pair(arg)                     = isa(arg, Expr) && (arg.head==:call) && (arg.args[1]==:(=>))
is_kwarg(arg)                    = isa(arg, Expr) && (arg.head==:(=))
is_tuple(arg)                    = isa(arg, Expr) && (arg.head==:tuple)
is_kwargexpr(arg)                = is_kwarg(arg) || ( (arg.head==:tuple) && all([is_kwarg(x) for x in arg.args]) )
complexarg_key(complexarg)       = complexarg.args[2]
complexarg_kwargexpr(complexarg) = complexarg.args[3]

function extract_voiceargs(arg)
    if (is_symbol(arg) || is_pair(arg)) return [arg]
    elseif is_tuple(arg)                return arg.args
    else                                @ArgumentError("the argument does not contain voiceargs.")
    end
end

function extract_kwargs(arg)
    if is_kwarg(arg)     return Dict{Symbol, Any}(split_kwarg(arg))
    elseif is_tuple(arg) return Dict{Symbol, Any}(split_kwarg(x) for x in arg.args)
    else                 @KeywordArgumentError("the argument does not contain keyword arguments.")
    end
end

function split_kwarg(kwarg)
    if !is_kwarg(kwarg) @KeywordArgumentError("the argument is not a keyword arguments.") end
    @capture(kwarg, key_ = value_)
    return key => value
end

function split_args(args)
    simpleargs  = Symbol[x for x in args if is_symbol(x)]
    complexargs = Expr[x for x in args if is_pair(x)]
    invalidargs = [x for x in args if !is_symbol(x) && !is_pair(x)]
    if !isempty(invalidargs) @ArgumentError("one or multiple arguments are not valid:\nInvalid: $(join(invalidargs, "\nInvalid: "))") end
    return simpleargs, complexargs
end

function split_complexargs(complexargs_expr)
    if !all(is_pair.(complexargs_expr)) @ArgumentError("not all of complexargs_expr are `voicearg=>kwargs` pairs.") end
    if !all([is_kwargexpr(complexarg_kwargexpr(x)) for x in complexargs_expr]) @ArgumentError("not all keyword arguments are parsable") end
    return Dict{Symbol, Dict{Symbol, Any}}(complexarg_key(x) => extract_kwargs(complexarg_kwargexpr(x)) for x in complexargs_expr)
end

function create_voiceargs_dict(simpleargs, complexargs::Dict{Symbol, Dict{Symbol, Any}})
    return merge(complexargs, Dict{Symbol, Dict{Symbol, Any}}(x => Dict{Symbol, Any}() for x in simpleargs))
end

function parse_signature(f_expr)
    f_def  = splitdef(f_expr)
    f_name = f_def[:name]
    f_args = Dict(x[1] => Dict(:arg_type=>x[2], :slurp=>x[3], :default=>x[4]) for x in splitarg.(f_def[:args]))
    if !is_symbol(f_name) @ArgumentError("the function name $f_name in the function definition is not valid.") end
    for arg_name in keys(f_args)
        if !is_symbol(arg_name) @ArgumentError("the argument $arg_name in the function signature is not valid.") end
    end
    return f_name, f_args
end

function eval_kwargs!(caller::Module, voiceargs)
    for voicearg in keys(voiceargs)
        for kwarg in keys(voiceargs[voicearg])
            if isa(voiceargs[voicearg][kwarg], Union{Symbol,Expr})
                voiceargs[voicearg][kwarg] = eval_arg(caller, voiceargs[voicearg][kwarg])
            end
        end
    end
end

function eval_arg(caller::Module, arg)
    try
        return @eval(caller, $arg)
    catch e
        @ArgumentEvaluationError("argument $arg could not be evaluated at parse time (in module $caller).")
    end
end

function validate_voiceargs(voiceargs, f_args)
    for voicearg in keys(voiceargs)
        if !haskey(f_args, voicearg) @ArgumentError("voicearg $(voicearg) is not part of the (positional) function arguments (keyword arguments cannot be voiceargs).") end
        for kwarg in keys(voiceargs[voicearg])
            if !haskey(VALID_VOICEARGS_KWARGS, kwarg) @KeywordArgumentError("keyword $(kwarg) of voicearg $(voicearg) is not valid.") end
            if !isa(voiceargs[voicearg][kwarg], VALID_VOICEARGS_KWARGS[kwarg]) @KeywordArgumentError("keyword argument $(kwarg)=$(voiceargs[voicearg][kwarg]) of voicearg $(voicearg) is of the wrong type (obtained type: $(typeof(voiceargs[voicearg][kwarg])); expected type: $(VALID_VOICEARGS_KWARGS[kwarg])).") end
        end
        if haskey(voiceargs[voicearg], :valid_input) && haskey(voiceargs[voicearg], :valid_input_auto) @IncoherentArgumentError("the keywords valid_input and valid_input_auto are incompatible. Set can set only one of them, not both.") end
        if  (haskey(voiceargs[voicearg], :vararg_end) || haskey(voiceargs[voicearg], :vararg_max) || haskey(voiceargs[voicearg], :vararg_timeout)) && !f_args[voicearg][:slurp] @KeywordArgumentError("the keywords vararg_end, vararg_max and vararg_timeout are only valid for vararg arguments (as e.g. `args...`)") end
        if !(haskey(voiceargs[voicearg], :vararg_end) || haskey(voiceargs[voicearg], :vararg_max) || haskey(voiceargs[voicearg], :vararg_timeout)) &&  f_args[voicearg][:slurp] @ArgumentError("at least one of the keywords vararg_end, vararg_max and vararg_timeout must be set for vararg arguments (as e.g. `args...`)") end
    end
end

# NOTE: This function requires that the voiceargs as given by the user have been validated (else it would randomly fail).
function handle_valid_inputs!(caller::Module, voiceargs, f_args)
    for voicearg in keys(voiceargs)
        if haskey(voiceargs[voicearg], :valid_input_auto) && (voiceargs[voicearg][:valid_input_auto] == true)
            type = eval_arg(caller, f_args[voicearg][:arg_type])
            voiceargs[voicearg][:valid_input] = generate_valid_input(type)
            if !isa(voiceargs[voicearg][:valid_input], VALID_VOICEARGS_KWARGS[:valid_input]) @ArgumentError("generation of valid input by type inspection failed for voicearg $voicearg (generated: $(voiceargs[voicearg][:valid_input])).") end
        end
    end
end

generate_valid_input(type::Type{<:Enum}) = [string.(instances(type))...]
generate_valid_input(type)               = @ArgumentError("type $(esc(type)) is not supported to define the valid input of a voicearg; currently supported are: Enum types. Other types can be supported by overloading `generate_valid_input`. Else, you can define the valid input explicitly witht the `valid_input` keyword argument.")

# Construct the wrapper based on the definition of the original function, as it will only differ in the arguments and body.
function wrap_f(f_name, f_args, f_expr, voiceargs)
    use_dynamic_recognizers = true
    if haskey(ENV, "JSI_USE_DYNAMIC_RECOGNIZERS") use_dynamic_recognizers = (parse(Int64, ENV["JSI_USE_DYNAMIC_RECOGNIZERS"]) > 0); end

    # Generate the code for the recognition of the voiceargs.
    recognitions = fill(:(begin end), length(voiceargs))
    token = gensym("token")
    i = 1
    for voicearg in keys(voiceargs)
        kwargs                   = voiceargs[voicearg]
        modelname                = haskey(kwargs,:model) ? kwargs[:model] : DEFAULT_MODEL_NAME
        use_partial_recognitions = haskey(kwargs,:use_max_speed) ? kwargs[:use_max_speed] : USE_PARTIAL_RECOGNITIONS_DEFAULT
        ignore_unknown           = haskey(kwargs,:ignore_unknown) ? kwargs[:ignore_unknown] : USE_IGNORE_UNKNOWN_DEFAULT
        f_arg                    = f_args[voicearg]
        f_name_sym               = :(Symbol($(string(f_name))))
        voicearg_sym             = :(Symbol($(string(voicearg))))
        voicearg_esc             = esc(voicearg)
        is_vararg                = f_args[voicearg][:slurp]
        if haskey(kwargs, :valid_input)
            valid_input = kwargs[:valid_input]
            if (use_dynamic_recognizers) recognizer_or_info = :(($f_name_sym, $voicearg_sym, $valid_input, $modelname))
            else                         recognizer_or_info = :(recognizer($f_name_sym, $voicearg_sym))
            end
        else
            recognizer_or_info = :(recognizer($modelname))
        end
        if is_vararg
            tic_call   = :(begin end)
            conditions = :()
            vararg_timeout = VARARG_TIMEOUT_DEFAULT
            if haskey(kwargs, :vararg_timeout) # NOTE: the timeout check shall always remain the first condition in order not to wait e.g. for an end keyword that is never coming...
                vararg_timeout = kwargs[:vararg_timeout]
                has_timed_out = gensym("has_timed_out")
                conditions = :(!$has_timed_out)
            end
            if haskey(kwargs, :vararg_max)
                vararg_max = kwargs[:vararg_max]
                reached_max = :(length($voicearg_esc) == $vararg_max)
                conditions = (conditions!=:()) ? :($conditions && !$reached_max) : :(!$reached_max)
            end
            if haskey(kwargs, :vararg_end) # NOTE: the end keyword check shall always remain the last condition in order not to wait for an end keyword that is never coming...
                vararg_end = kwargs[:vararg_end]
                is_end = :(_is_next($vararg_end, $recognizer_or_info, noises($modelname); consume_if_match=true, use_partial_recognitions=$use_partial_recognitions))
                conditions = (conditions!=:()) ? :($conditions && !$is_end) : :(!$is_end)
            end
            recognition = quote
                $voicearg_esc = []
                $has_timed_out = false
                while $conditions
                    $token = next_token($recognizer_or_info, noises($modelname); use_partial_recognitions=$use_partial_recognitions, timeout=$vararg_timeout, ignore_unknown=$ignore_unknown)
                    if ($token == UNKNOWN_TOKEN) @InsecureRecognitionException("@voiceargs: argument not recognised.") end
                    push!($voicearg_esc, $(interpret_and_parse_calls(token, kwargs, f_arg)))
                    $has_timed_out = ($token == "")
                end
            end
        else
            recognition = quote
                $token = next_token($recognizer_or_info, noises($modelname); use_partial_recognitions=$use_partial_recognitions, ignore_unknown=$ignore_unknown)
                if ($token == UNKNOWN_TOKEN) @InsecureRecognitionException("@voiceargs: argument not recognised.") end
                if ($token == "") @InsecureRecognitionException($("time out waiting for voice argument $voicearg in function $f_name.")) end
                $voicearg_esc = $(interpret_and_parse_calls(token, kwargs, f_arg))
            end
        end
        recognitions[i] = recognition
        i += 1
    end
    voiceargs_recognition = quote $(recognitions...) end

    # Generate the call of the original function.
    f_def = splitdef(f_expr)
    f_call_args = [splitarg(x)[3] ? :($(splitarg(x)[1])...) : splitarg(x)[1] for x in f_def[:args]] # NOTE: the f_args dict is not used here as it does not preserve the argument order.
    f_call_kwargs = [:($(splitarg(x)[1]) = $(splitarg(x)[1])) for x in f_def[:kwargs]]
    f_call = :($(esc(f_name))($(esc.(f_call_args)...); $(esc.(f_call_kwargs)...)))

    # Assign new arguments and body, and escape the other parts of the function definition.
    f_def[:name]        = esc(f_def[:name])
    f_def[:kwargs]      = esc.(f_def[:kwargs])
    f_def[:whereparams] = esc.(f_def[:whereparams])
    f_def[:args]        = [esc(x) for x in f_def[:args] if !haskey(voiceargs, splitarg(x)[1])] # The voiceargs will not be part of the wrapper function signature.
    f_def[:body]        = quote
        $voiceargs_recognition
        $f_call
    end
    return combinedef(f_def)
end

function interpret_and_parse_calls(input, kwargs, f_arg)
    type = esc(f_arg[:arg_type])
    if haskey(kwargs, :interpret_function)
        interpret_function = esc(kwargs[:interpret_function])
        return :(parse($type, $interpret_function($input)))
    else
        return :(parse($type, $input))
    end
end

import Base.parse

parse(type::Type{String}, val::AbstractString) = string(val)

function parse(type::Type{<:Enum}, val::AbstractString)
    parser = Dict(string(x) => x for x in instances(type))
    if !haskey(parser, val)
        reset_all()
        @InsecureRecognitionException("Insecurity in recognition: token '$val' was not expected (it is not parsable to type '$type')")
    end
    return parser[val]
end

## TEMPORARY FUNCTION DEFINITIONS TO BE MERGED IN MACROTOOLS (https://github.com/FluxML/MacroTools.jl/pull/173)

isdef(ex)     = isshortdef(ex) || islongdef(ex)
islongdef(ex) = @capture(ex, function (fcall_ | fcall_) body_ end)
isshortdef(ex) = MacroTools.isshortdef(ex)


text:
# First, generate mock data
include("generate_mockdata.jl")

# Then, run tests on mock data
using Test
using GEOTRACES
using Unitful
using OceanographyCruises
using DataFrames

#@testset "Test dummy cruise" begin
#    ct = dummy_cruise_track()
#    @test ct isa CruiseTrack{Float64,DateTime}
#    @test ct.name isa String
#    @test ct.name == "Dummy Cruise Track"
#    @test ct.lon isa Vector{Float64}
#    @test ct.lat isa Vector{Float64}
#    @test ct.date isa Vector{DateTime}
#end

@testset "Cruises" begin
    IDP21_cruise_list = ["GA01", "GA02", "GA03", "GA04N", "GA04S", "GA06", "GA08", "GA10", "GA11", "GA13", "GAc01", "GAc02", "GAc03", "GApr08", "GApr09", "GI01", "GI02", "GI03", "GI04", "GI05", "GI06", "GIPY01", "GIPY02", "GIPY04", "GIPY05", "GIPY06", "GIPY11", "GIPY13", "GIpr01", "GIpr05", "GIpr06", "GN01", "GN02", "GN03", "GN04", "GN05", "GP02", "GP06", "GP12", "GP13", "GP15", "GP16", "GP18", "GP19", "GPc01", "GPc02", "GPc03", "GPc05", "GPc06", "GPpr01", "GPpr02", "GPpr04", "GPpr05", "GPpr07", "GPpr08", "GPpr10", "GPpr11", "GS01", "GSc01", "GSc02"]
    @test unique(GEOTRACES.list_of_cruises()) == IDP21_cruise_list
    @test GEOTRACES.list_of_stations() isa Vector{String}
    @test GEOTRACES.cruisetrack("GA03") isa CruiseTrack
end

@testset "observations" begin
    @testset "$var" for var in ["Cd", "PO₄", "δCd", "Fe"]
        x = GEOTRACES.observations(var)
        @test x isa DataFrame
        @test getproperty(x, var) isa Vector{<:Quantity}
        x_QC = GEOTRACES.qualitycontrols(var)
        @test x_QC isa Vector{Int}
    end
    @testset "($var1, $var2)" for (var1, var2) in zip(["Cd", "PO₄"], ["δCd", "Fe"])
        x = GEOTRACES.observations(var1, var2)
        @test x isa DataFrame
        @test getproperty(x, var1) isa Vector{<:Quantity}
        @test getproperty(x, var2) isa Vector{<:Quantity}
        x_QC, y_QC = GEOTRACES.qualitycontrols(var1, var2)
        @test x_QC isa Vector{Int}
        @test y_QC isa Vector{Int}
        @test length(x_QC) == length(y_QC)
    end
end

@testset "transects" begin
    @testset "$var" for var in ["Cd", "PO₄", "δCd", "Fe"]
        x = GEOTRACES.transects(var)
        @test x isa Transects
        x = GEOTRACES.transect(var; cruise="GA02")
        @test x isa Transect
    end
    @testset "($var1, $var2)" for (var1, var2) in zip(["Cd", "PO₄"], ["δCd", "Fe"])
        obs = GEOTRACES.observations(var1, var2)
        x = GEOTRACES.transects(obs, var1)
        y = GEOTRACES.transects(obs, var2)
        @test x isa Transects
        @test y isa Transects
    end
end

@testset "helper functions" begin
    @test GEOTRACES.matchingvariables("cd") isa Vector{Pair{String, String}}
    @testset "variable names" begin
        @test GEOTRACES.varname("He") == "var22"
        @test GEOTRACES.varname("O2") == "var35"
        @test GEOTRACES.varname("P") == "var37"
        @test GEOTRACES.varname("NO3") == "var40"
        @test GEOTRACES.varname("Si") == "var39"
        @test GEOTRACES.varname("Cd") == "var85"
        @test GEOTRACES.varname("δCd") == "var158"
        @test GEOTRACES.varname("Fe") == "var88"
        @test GEOTRACES.varname("δFe") == "var160"
        @test GEOTRACES.varname("Ni") == "var104"
        @test GEOTRACES.varname("Nd") == "var259"
        @test GEOTRACES.varname("eNd") == "var251"
        @test GEOTRACES.varname("εNd") == "var251"
        @test GEOTRACES.varname("V") == "var109"
    end
end


text:
using DiffEqProblemLibrary
using Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems
using DiffEqProblemLibrary.DDEProblemLibrary: importddeproblems
using DiffEqProblemLibrary.DAEProblemLibrary: importdaeproblems
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems
using DiffEqProblemLibrary.JumpProblemLibrary: importjumpproblems

importodeproblems()
importddeproblems()
importjumpproblems()
importsdeproblems()
importdaeproblems()

# The test is simply that all of the examples build!


text:
using Documenter, Microbiome

makedocs(
    format = :html,
    sitename = "Microbiome.jl",
    pages = [
        "Home" => "index.md",
        "Microbial Abundances" => "abundances.md",
        "Distances & Dissimilarity" => "distances.md",
        "Contributing" => "contributing.md"
    ],
    authors = "Kevin Bonham, PhD"
)

deploydocs(
    repo = "github.com/BioJulia/Microbiome.jl.git",
    julia = "0.6",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
export qft_circuit

"""
    cphase(i, j)

Control phase gate.
"""
cphase(i, j) = control(i, j=> shift(2π/(2^(i-j+1))));

hcphases(n, i) = chain(n, i==j ? put(i=>H) : cphase(j, i) for j in i:n);

"""
    qft_circuit(n)

Create a Quantum Fourer Transform circuit. See also [`QFT`](@ref).
"""
qft_circuit(n::Int) = chain(n, hcphases(n, i) for i = 1:n)

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
struct Model
    T::DataType
    params::Params
    constants::Constants
    geometry::Geometry
    current_datetime::DateTime
    end_datetime::DateTime
    spectral_trans::SpectralTrans
    boundaries::Boundaries
    prognostics::Prognostics
    horizontal_diffusion::HorizontalDiffusion
    implicit::Implicit
end

function Model(;
    # Model grid point resolution
    nlon, nlat,
    # Number of model levels
    nlev,
    # Model spectral resolution
    trunc,
    # Number of time steps per day
    n_steps_day,
    # Real number type
    real_type = Float64,
    # Radius of Earth
    Rₑ = 6.371e+6,
    # Angular frequency of Earth's rotation
    Ω = 7.292e-05,
    # Gravitational acceleration
    g = 9.81,
    # Start and end datetimes
    start_datetime = DateTime(1982,1,1),
    end_datetime = DateTime(1982,1,2),
    # Name of restart file
    restart_file = nothing,
    # Number of tracers (including humidity)
    n_tr = 1,
    # Frequency of diagnostic checks
    n_diag = 36*5
    )

    # Default gas parameters
    akap = 2.0/7.0
    R    = akap*1000.4

    # Default dynamical constant parameters
    γ      = 6.0       # Reference temperature lapse rate (-dT/dz in deg/km)
    hscale = 7.5       # Reference scale height for pressure (in km)
    hshum  = 2.5       # Reference scale height for specific humidity (in km)
    refrh1 = 0.7       # Reference relative humidity of near-surface air

    current_datetime = start_datetime

    params = Params(n_diag, n_steps_day, 86400.0/real_type(n_steps_day), real_type(0.5),)
    constants = Constants(real_type, Rₑ, Ω, g, akap, R, γ, hscale, hshum, refrh1)
    geometry = Geometry(real_type, constants, nlon, nlat, nlev, trunc)
    spectral_trans = SpectralTrans(real_type, geometry, constants.Rₑ)
    boundaries = Boundaries(real_type, geometry, spectral_trans, g)

    if restart_file == nothing
        prognostics = initialize_from_rest(real_type, geometry, constants, boundaries.ϕ₀ₛ,
                                           spectral_trans, n_tr)
       # Write initial data
       output(geometry, params, spectral_trans, current_datetime, 1, prognostics)
    else
        throw("Restart functionality not yet implemented")
    end

    horizontal_diffusion = HorizontalDiffusion(real_type, geometry, constants)
    implicit = Implicit(real_type, geometry, constants, params, horizontal_diffusion)

    Model(real_type, params, constants, geometry, current_datetime, end_datetime, spectral_trans, boundaries,
          prognostics, horizontal_diffusion, implicit)
end


text:
using RLinearAlgebra
using LinearAlgebra
using UnicodePlots

n = 10
d = 10
κ = 8.0

A = RLinearAlgebra.generate_matrix(n, d, d, κ)
b = randn(n)

sol = LinearSolver(TypeRPM())
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Default RPM", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerKaczmarzCYC()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Cyclic Kaczmarz", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerKaczmarzCYC(), ProjectionLowCore()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Cyclic with Low Core", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerKaczmarzCYC(), ProjectionFullCore()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Cyclic with Full Core", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerMotzkin()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Motzkin", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerMotzkin(), ProjectionLowCore()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Motzkin low core", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerMotzkin(), ProjectionFullCore()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Motzkin full core", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerGaussSketch()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Gauss", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerGaussSketch(), ProjectionLowCore()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Gauss Low Core", xlabel = "iteration", ylabel = "residual")
println(plt)

sol = LinearSolver(TypeRPM(SamplerGaussSketch(), ProjectionFullCore()))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Gauss Full Core", xlabel = "iteration", ylabel = "residual")
println(plt)

samp = SamplerMotzkin()
samp.sampled = true
sol = LinearSolver(TypeRPM(samp))
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "RPM Convergence", xlabel = "iteration", ylabel = "residual")
println(plt)

println("Plot collection of samplers")
plt = lineplot(0:500, zeros(501), width = 100, height = 20, xlabel = "iteration", ylabel = "residual",
               ylim = [0, norm(b)])
samplers = RPMSamplers()
for s in samplers
    solt = LinearSolver(TypeRPM(s))
    xx = rsolve(solt, A, b)
    lineplot!(plt, solt.log.residual_hist)
end
println(plt)

println("Randomized Gauss seidel")
sol = LinearSolver(TypeRGS())
x = rsolve(sol, A, b)
plt = lineplot(sol.log.residual_hist, title = "Randomized Gauss Seidel", xlabel = "iteration", ylabel = "residual")
println(plt)


text:
using OrdinaryDiffEq
using Flux
using Flux: logitcrossentropy
using MLDatasets: MNIST
using MLDataUtils
using DiffEqFlux
using CuArrays; CuArrays.allowscalar(false)
using NNlib
using IterTools: ncycle

function loadmnist(batchsize=bs)
	# Use MLDataUtils LabelEnc for natural onehot conversion
    onehot(labels_raw) = convertlabel(LabelEnc.OneOfK, labels_raw, LabelEnc.NativeLabels(collect(0:9)))
	# Load MNIST
	imgs, labels_raw = MNIST.traindata();
	# Process images into (H,W,C,BS) batches
	x_train = reshape(imgs,size(imgs,1),size(imgs,2),1,size(imgs,3))|>gpu
	x_train = batchview(x_train,batchsize);
	# Onehot and batch the labels
	y_train = onehot(labels_raw)|>gpu
	y_train = batchview(y_train,batchsize)
	return x_train, y_train
end

# Main
const bs = 128
x_train, y_train = loadmnist(bs);

down = Chain(
             Conv((3,3),1=>64,relu,stride=1), GroupNorm(64,64),
             Conv((4,4),64=>64,relu,stride=2,pad=1), GroupNorm(64,64),
             Conv((4,4),64=>64,stride=2,pad=1),
            )|>gpu;
dudt = Chain(
           Conv((3,3),64=>64,relu,stride=1,pad=1),
           Conv((3,3),64=>64,relu,stride=1,pad=1)
          ) |>gpu;
fc = Chain(GroupNorm(64,64),
           x->relu.(x),
           MeanPool((6,6)),
           x -> reshape(x,(64,bs)),
           Dense(64,10)
          )|>gpu;

solver_kwargs = Dict(:save_start=>false,:save_everystep=>false, :reltol=>1e-3, :abstol=>1e-3)
neural_ode_layer = DiffEqFlux.NeuralODE(dudt, (0.f0,1.f0), Tsit5(), solver_kwargs)

model = Chain(
              down,             #(28,28,1,BS) -> (6,6,64,BS)
              neural_ode_layer, #(6,6,64,BS) -> (6,6,64,BS)
              fc                #(6,6,64,BS) -> (10, BS)
             )

# Showing this works with forward solve
x_m = model(x_train[1])

# Define loss and test with forward solve
function loss(x,y)
    y_hat = model(x)
    return logitcrossentropy(y_hat,y)
end

loss(x_train[1],y_train[1])

# Define accuracy, diagnostic but not objective
classify(x) = argmax.(eachcol(x))

function accuracy(model,data; n_batches=100)
    total_correct = 0
    total = 0
    for (x,y) in collect(data)[1:n_batches]
        target_class = classify(cpu(y))
        predicted_class = classify(cpu(model(x)))
        total_correct += sum(target_class .== predicted_class)
        total += length(target_class)
    end
    return total_correct/total
end

accuracy(model, zip(x_train,y_train))

# Logging and accuracy in callback
using TensorBoardLogger, Logging

lg=TBLogger("tensorboard_logs/run", min_level=Logging.Info)
iter = 0
cb() = begin
    global iter += 1
    li = Flux.data(loss(x_train[1],y_train[1]))
    with_logger(lg) do
        @info "loss" li
    end
    if iter%100 == 0
        ai = accuracy(model, zip(x_train,y_train))
        with_logger(lg) do
            @info "accuracy" ai
        end
        ai > 0.97 && Flux.stop()
    end
end

#Main Train Loop
opt=ADAM()
Flux.train!(loss,params(model),ncycle(zip(x_train,y_train),100),opt, cb=cb)

# Save and load trained models for plotting
using BSON: @save, @load

#model_cpu = cpu(model)
#@save "saved_models/mnist_node.bson" model_cpu

@load "saved_models/mnist_node.bson" model_cpu

dwn_loaded = gpu(model_cpu[1])
node_loaded= gpu(model_cpu[2])

dense_solver_kwargs = Dict(:saveat=> collect(0.f0:0.01:1.f0),:reltol=>1e-3, :abstol=>1e-3)
dense_node = DiffEqFlux.NeuralODE(node_loaded.model,node_loaded.tspan,node_loaded.solver,node_loaded.args,dense_solver_kwargs)

# Animation code for JuliaCon talk:
slv = cpu(Flux.data(dense_node(dwn_loaded(x_train[1]))))

using Plots; pyplot();

cs = palette(:default)

clfy(x) = argmax.(eachcol(x))

for ti in 1:size(slv,5)
    plot(
         collect(eachrow(slv[2,3,7,:,1:ti])),
         collect(eachrow(slv[3,3,7,:,1:ti])),
         legend=false,
         ticks=false,
         c=cs[clfy(y_train[1])]',
        opacity=0.5,
        xlim= extrema(slv[2,3,7,:,1:end]),
        ylim = slv[3,3,7,:,1:end])
    scatter!(
         collect(eachrow(slv[2,3,7,:,ti])),
         collect(eachrow(slv[3,3,7,:,ti])),
         legend=false,
         ticks=false,
         c=cs[clfy(y_train[1])]')
    savefig("plots/mnist_dynamics_anim/$(lpad(ti,4,"0")).png")
end

max_stars_repo_name: jessebett/model-zoo

max_issues_repo_name: jessebett/model-zoo

max_forks_repo_name: jessebett/model-zoo


text:
include("help.jl")
function count_rectangular_borders(r::Robot)::Int
    moves!(r, Sud)
    moves!(r, Ost)

    count = 0
    side = West

    while !isborder(r, Nord)
        possible, is_rectangular_border = move_if_possible!(r, side)

        while possible
            if is_rectangular_border
                count += 1
            end

            possible, is_rectangular_border = move_if_possible!(r, side)
        end

        side = inverse(side)
        move!(r, Nord)
    end

    return count
end

function move_if_possible!(r::Robot, direct_side::HorizonSide)::Tuple{Bool, Bool}
    if direct_side == Ost
        orthogonal_side = Sud
    else
        orthogonal_side = next(direct_side)
    end

    reverse_side = inverse(orthogonal_side)
    num_steps=0

    while isborder(r, direct_side) == true
        if isborder(r, orthogonal_side) == false
            move!(r, orthogonal_side)
            num_steps += 1
        else
            break
        end
    end

    if isborder(r,direct_side) == false
        move!(r,direct_side)
        if num_steps > 0
            while isborder(r,reverse_side) && !isborder(r,direct_side)
                move!(r,direct_side)
            end
        end
        result = true
    else
        result = false
    end

    for _ in 1:num_steps
        move!(r, reverse_side)
    end

    return result, num_steps == 1 && result
end

function moves!(r::Robot,side::HorizonSide)
    possible, _ = move_if_possible!(r,side)
    while possible
        possible, _ = move_if_possible!(r,side)
    end
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
function lin_interp(grid, vals)
    function func(x::Real)
        if x <= grid[1]
            index_1 = 1
            index_2 = 2
        elseif x >= grid[length(grid)]
            index_1 = length(grid) - 1
            index_2 = length(grid)
        else
            index_1 = searchsortedfirst(grid, x)
            index_2 = searchsortedlast(grid, x)
        end
        x_1 = grid[index_1]
        x_2 = grid[index_2]
        y_1 = vals[index_1]
        y_2 = vals[index_2]

        y = ((y_2- y_1)/(x_2 - x_1))*(x - x_1) + y_1
        return y
    end

    function func{T<:Real}(x::AbstractVector{T})
        y = zeros(length(x))
        for i in 1:length(x)
            y[i] = func(x[i])
        end
        return y
    end

    return func
end


text:
using NeuralQuantum
using LinearAlgebra, SparseArrays
using Test

function test_cache(sc_1)
    @testset "test zero" begin
        zero!(sc_1)
        test_zero(sc_1)
    end

    @testset "test addition" begin
        sc_2 = deepcopy(sc_1)
        random_data!(sc_1)
        random_data!(sc_2)
        sc_3 = deepcopy(sc_1)

        NeuralQuantum.add!(sc_3, sc_2)
        test_sum!(sc_1, sc_2, sc_3)
    end
end

Nsites = 4
lind = quantum_ising_lind(SquareLattice([Nsites],PBC=true), g=1.3, V=2.0, γ=1.0)
net  = cached(NDM(Nsites, 2, 3))

@testset "Test sampling caches" begin
    include("MCMC_SR_cache.jl")
    @testset "Test sampling cache MCMCSREvaluationCache" begin
        prob = SteadyStateProblem(lind, operators=false, variance=false)
        sc   = SamplingCache(SR(), prob, net)
        @test sc isa MCMCSREvaluationCache
        test_cache(sc)
    end

    include("MCMC_SRL_cache.jl")
    @testset "Test sampling cache MCMCSRLEvaluationCache" begin
        prob = SteadyStateProblem(lind, operators=false, variance=true)
        sc   = SamplingCache(SR(), prob, net)
        @test sc isa MCMCSRLEvaluationCache
        test_cache(sc)
    end

    include("MCMC_Grad_cache.jl")
    @testset "Test sampling cache MCMCGradientEvaluationCache" begin
        prob = SteadyStateProblem(lind, operators=false, variance=false)
        sc   = SamplingCache(Gradient(), prob, net)
        @test sc isa MCMCGradientEvaluationCache
        test_cache(sc)
    end

    include("MCMC_GradL_cache.jl")
    @testset "Test sampling cache MCMCGradientLEvaluationCache" begin
        prob = SteadyStateProblem(lind, operators=false, variance=true)
        sc   = SamplingCache(Gradient(), prob, net)
        @test sc isa MCMCGradientLEvaluationCache
        test_cache(sc)
    end
end


text:
function LUDecomposer(arg0::jdouble)
    return LUDecomposer((jdouble,), arg0)
end

function decompose(obj::LUDecomposer, arg0::RealMatrix)
    return jcall(obj, "decompose", DecompositionSolver, (RealMatrix,), arg0)
end


text:
"""
    read_gmsh_mesh(filename) -> mesh::Mesh
"""
function read_gmsh_mesh(fn::AbstractString; kwargs...)
    open(fn, "r") do io
        read_gmsh_mesh(io; kwargs...)
    end
end

"""
    read_gmsh_mesh(iostream) -> mesh::Mesh

Reads the mesh nodes and elements stored in the input .msh file (`io`, output by GMSH)
    into arrays of node vectors and vertex vectors respectively.

Returns an object `mesh` of type `Mesh`, comprising both vector arrays.
"""
function read_gmsh_mesh(io; physical=nothing, dimension=2, sort=true, T=Float64)

    entity_tag = 0
    if physical != nothing
        @assert physical isa String
        while true
            thisLine = io |> readline |> strip
            thisLine == "\$PhysicalNames" && break
            thisLine == "\$Nodes" && error("Mesh file does not define physical entities")
        end
        thisLine = io |> readline |> strip
        s = split(thisLine)
        NP = parse(Int,s[1])
        println("Mesh file defined $(NP) physical entities:")
        for i in 1:NP
            thisLine = io |> readline |> strip
            s = split(thisLine)
            entity_name = s[3]
            @show entity_name
            if entity_name == "\"$(physical)\""
                entity_tag = parse(Int, s[2])
                println("Target entity \"$(physical)\" is $i out of $(NP).")
                break
            end
        end
        entity_tag == 0 && error("Specified physical entitiy not found")
    end

    thisLine = io |> readline |> strip
    while thisLine != "\$Nodes"
        thisLine = io |> readline |> strip
    end

    thisLine = io |> readline |> strip
    s = split(thisLine)
    NV = parse(Int, s[1])

    P = SVector{3,T}
    v = Vector{P}(undef,NV)
    for i in 1:NV
        thisLine = io |> readline |>  strip
        d = readdlm(IOBuffer(thisLine), T)
        v[i] = P(d[2], d[3], d[4])
    end

    while thisLine != "\$Elements"
        thisLine = io |> readline |> strip
    end

    thisLine = io |> readline |> strip
    s = split(thisLine)
    NF = parse(Int, s[1])

    dim2type = [1,2,4]
    type = dim2type[dimension]

    thisLine = io |> readline |> strip
    I = SVector{dimension+1,Int}
    f = Vector{I}(undef,NF)
    i = 0
    while thisLine != "\$EndElements"
        d = readdlm(IOBuffer(thisLine), Int)
        thisLine = io |> readline |> strip
        d[2] == type || continue
        entity_tag == d[4] || entity_tag == 0 || continue
        # f[i +=1] = SVector(d[end-2], d[end-1], d[end-0])
        f[i+=1] = d[end-dimension:end]
    end
    resize!(f,i)

    @show length(v)
    @show length(f)

    if sort
        @info "sorting..."
        f = sort_sfc(v,f)
    end
    # Q = Pt{3,Float64}
    # ctrs = Q[]
    # for tr in f
    #         ctr = sum(v[tr])/3
    #         push!(ctrs,ctr)
    # end
    # sorted = sort_sfc(ctrs)
    # f = f[sorted]

    return Mesh(v, f)

end

function sort_sfc(v, f)
    T = eltype(eltype(v))
    Q = Pt{3,T}
    ctrs = Q[]
    for tr in f
            ctr = sum(v[tr])/3
            push!(ctrs,ctr)
    end
    sorted = sort_sfc(ctrs)
    return f[sorted]
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
a = eye(100)
b = eye(100)

dm_a = data(a)
dm_b = data(b)

dm_a_na = copy(dm_a)
dm_a_na[:, :] = NA
dm_b_na = copy(dm_b)
dm_b_na[:, :] = NA

f1() = *(a, b)
f2() = *(dm_a, dm_b)
f3() = *(dm_a_na, dm_b_na)

df1 = benchmark(f1,
	            "Linear Algebra",
                "Matrix Multiplication w/ No NA's",
                1_000)
df2 = benchmark(f2,
	            "Linear Algebra",
                "DataMatrix Multiplication w/ No NA's",
                1_000)
df3 = benchmark(f3,
	            "Linear Algebra",
                "DataMatrix Multiplication w/ NA's",
                1_000)

# TODO: Keep permanent record
printtable(vcat(df1, df2, df3), header=false)

# Compare with R
# We're 10x as fast!
# a <- diag(100)
# b <- diag(100)
# a %*% b
# s <- Sys.time(); a %*% b; e <- Sys.time(); e - s


text:
function update_feed!(
    state::Tuple{AbstractGraph, AbstractArray}, agent_idx::Integer,
    config::Config
)
    graph, agent_list = state
    this_agent = agent_list[agent_idx]
    unique!(this_agent.feed)
    deleted_tweets = Integer[]
    for (index, tweet) in enumerate(this_agent.feed)
        if tweet.weight == -1 || !(tweet.source_agent in inneighbors(graph, agent_idx))
            push!(deleted_tweets, index)
        else
            tweet.weight = config.feed_props.tweet_decay * tweet.weight
        end
    end
    deleteat!(this_agent.feed, deleted_tweets)
    sort!(this_agent.feed, lt=<, rev=true)
    if length(this_agent.feed) > config.feed_props.feed_size
        this_agent.feed = this_agent.feed[1:config.feed_props.feed_size]
    end
    return state
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
using JLD2, Test

function writeloop(f, sz)
    for i = 1:sz
        write(f, UInt8(6))
    end
end

# Force growing an MmapIO
name, io = mktemp()
f = JLD2.MmapIO(name, true, true, true)
close(io)
sz = JLD2.MMAP_GROW_SIZE+5
writeloop(f, sz)
JLD2.truncate_and_close(f, position(f))

@test filesize(name) == sz
f = open(name)
rd = read!(f, Array{UInt8}(undef, sz))
@test all(rd .== UInt8(6))
close(f)

f = JLD2.MmapIO(name, true, true, true)
write(f, rd)
JLD2.truncate_and_close(f, position(f))

@test filesize(name) == sz
f = open(name)
rd = read!(f, Array{UInt8}(undef, sz))
@test all(rd .== UInt8(6))
close(f)

@test JLD2.size_flag(1) === UInt8(0)
@test JLD2.size_flag(256) === UInt8(1)
@test JLD2.size_flag(65536) === UInt8(2)
@test JLD2.size_flag(UInt64(4294967296)) === UInt8(3)

buf = IOBuffer()
JLD2.write_size(buf, 175)
seek(buf, 0)
@test read(buf, UInt8) == 175
seek(buf, 0)
@test JLD2.read_size(buf, UInt8(0)) == 175
seek(buf, 0)

JLD2.write_size(buf, 42580)
seek(buf, 0)
@test read(buf, UInt16) == 42580
seek(buf, 0)
@test JLD2.read_size(buf, UInt8(1)) == 42580
seek(buf, 0)

JLD2.write_size(buf, 1902153053)
seek(buf, 0)
@test read(buf, UInt32) == 1902153053
seek(buf, 0)
@test JLD2.read_size(buf, UInt8(2)) == 1902153053
seek(buf, 0)

if Int == Int64
    # Only test this on 64-bit platforms. We can't read data of size >4GB on 32-bit anyway.
    JLD2.write_size(buf, 3804306107)
    seek(buf, 0)
    @test read(buf, UInt32) == 3804306107
    seek(buf, 0)
    @test JLD2.read_size(buf, UInt8(2)) == 3804306107
    seek(buf, 0)

    JLD2.write_size(buf, 7832227080891617460)
    seek(buf, 0)
    @test read(buf, UInt64) == 7832227080891617460
    seek(buf, 0)
    @test JLD2.read_size(buf, UInt8(3)) == 7832227080891617460
    seek(buf, 0)
end

@test JLD2.size_size(1) === 1
@test JLD2.size_size(256) === 2
@test JLD2.size_size(65536) === 4
@test JLD2.size_size(UInt64(4294967296)) === 8

# test with one of its own types
@test JLD2.typename(JLD2.BasicDatatype) == "JLD2.BasicDatatype"
@test JLD2.typename(String) == "Core.String"

module Foo
    module Bar
        struct Baz end
    end
end

@test endswith(JLD2.typename(Foo.Bar.Baz), "Foo.Bar.Baz")


text:
"""
    rekursif(f,a,b,n)

berisi program untuk mencari integral numerik dari fungsi
`f(x)` pada interval `[a,b]` menggunakan Aturan Rekursif Trapesium.
Program ini secara default berisi 4 masukan, yaitu fungsi `f(x)`,
titik ujung interval `[a,b]` dan banyaknya iterasi `n`

# Example
```jl
julia> f(x) = 1/x;

julia> T = rekursif(f,1,5,5)
6×1 Array{Float64,2}:
 2.4
 1.8666666666666667
 1.6833333333333333
 1.628968253968254
 1.6144063238103485
 1.6106858960792332
```
return solusi integral numerik pada setiap iterasi `T`
"""
function rekursif(f,a,b,n)
  M = 1;
  h = b-a;
  T = Array{Float64}(undef,n+1,1)
  T[1] = h/2*(f(a)+f(b))
  for j = 1:n
    h = h/2;
    s = 0;
    for k = 1:M
      x = a+h*(2*k-1);
      s = s+f(x);
    end
    M = 2*M;
    T[j+1] = T[j]/2+h*s;
  end
  return T
end


text:
import AutoGrad: addto!, addtoindex!, matches
using CUDA: CuArray, CuPtr
using Knet.KnetArrays: KnetArray, DevArray, Cptr
using Knet.LibKnet8: @knet8
using AutoGrad: Sparse, recording

function addto!(a::DevArray{T},b::DevArray{T}) where {T}
    if recording(); a = copy(a); end  # support highorder gradients
    axpy!(1,b,a) # (a+b)
end

function addto!(a::DevArray,b::Sparse)
    @assert matches(a, b.container)
    if recording(); a = copy(a); end  # support highorder gradients
    for (idx,val) in zip(b.indices, b.values)
        addtoindex!(a, val, idx...)
    end
    return a
end

addto!(a::Sparse, b::DevArray) = addto!(b, a)

matches(a::DevArray,b::DevArray)=(typeof(a)===typeof(b) && size(a)===size(b))

# Setindex! works for addtoindex! when there are no repeated indices. This is true for index types:
# Real, (Real...), CartesianIndex, Colon, AbstractArray{Bool}, Range, EmptyArray
# and pairs of Union{Real,AbstractUnitRange,Colon} and (Colon,Range)

# AutoGrad defines the following 3 methods for addtoindex!
# This is most specific, catches cases of I::Union{Real, AbstractArray}... that are safe.
# [1] addtoindex!(A::AbstractArray, X, I::Union{Colon, Real, AbstractRange, AbstractArray{Bool,N} where N}...) in AutoGrad at /dev/shm/dyuret/.julia/packages/AutoGrad/VFrAv/src/addto.jl:133
# This specifically catches cases where indices can repeat (i.e. arrays of ints or cartesianindices) and handles them
# [2] addtoindex!(A::AbstractArray, x, I::Union{Real, AbstractArray}...) in AutoGrad at /dev/shm/dyuret/.julia/packages/AutoGrad/VFrAv/src/addto.jl:107
# This is most general and handles cases that are not caught by #2 that are safe.
# [3] addtoindex!(A::AbstractArray, X, I...) in AutoGrad at /dev/shm/dyuret/.julia/packages/AutoGrad/VFrAv/src/addto.jl:135

# This extends #1, do KnetArray and CuArray separately because CuArray <: AbstractArray which may cause ambiguity
addtoindex!(A::KnetArray, X, I::Union{Real,Colon,AbstractRange,AbstractArray{Bool}}...) = setindex!(A, addto!(getindex(A,I...), X), I...)
addtoindex!(A::CuArray, X, I::Union{Real,Colon,AbstractRange,AbstractArray{Bool}}...) = setindex!(A, addto!(getindex(A,I...), X), I...)

# These extend #2, however only a small number of cases have been implemented, the rest should throw method errors
addtoindex!(A::KnetArray, X, I::Union{Real, AbstractArray}...) = _addtoindex!(A, X, I...)
addtoindex!(A::CuArray, X, I::Union{Real, AbstractArray}...) = _addtoindex!(A, X, I...)
addtoindex!(A::KnetArray, X, I::Union{Real, Colon, AbstractArray}...) = _addtoindex!(A, X, I...)
addtoindex!(A::CuArray, X, I::Union{Real, Colon, AbstractArray}...) = _addtoindex!(A, X, I...)
addtoindex!(A::KnetArray, X, I::AbstractArray{<:CartesianIndex}) = _addtoindex!(A, X, c2i(size(A),I))
addtoindex!(A::CuArray, X, I::AbstractArray{<:CartesianIndex}) = _addtoindex!(A, X, c2i(size(A),I))
c2i(d::Dims,i::AbstractArray{<:CartesianIndex}) = Int32[(LinearIndices(d))[c.I...] for c in i]

# This extends #3, which is the remaining cases safe to setindex!
addtoindex!(A::KnetArray, X, I...) = setindex!(A, addto!(getindex(A,I...), X), I...)
addtoindex!(A::CuArray, X, I...) = setindex!(A, addto!(getindex(A,I...), X), I...)

# These are the actual implementations for #2 that can handle repeat indices:
for (R,P) in ((CuArray,CuPtr),(KnetArray,Ptr)), T in (Float16, Float32, Float64); F = sizeof(T)*8
    @eval begin

        function _addtoindex!(A::$R{$T}, X, I::AbstractArray{<:Real})
            I = convert($R{Int32},I)
            X = convert($R{$T},X)
            @knet8($("addents_$F"),(Cint,$P{Int},$P{$T},$P{$T}),
                   length(I), I, A, X)
            return A
        end

        function _addtoindex!(A::$R{$T}, X, ::Colon, I::AbstractArray{<:Real})
            I = convert($R{Int32},I)
            X = convert($R{$T},X)
            @knet8($("addcols_$F"),(Cint,Cint,Cint,$P{Int},$P{$T},$P{$T}),
                   size(A,1), size(A,2), length(I), I, A, X)
            return A
        end

        function _addtoindex!(A::$R{$T}, X, I::AbstractArray{<:Real}, ::Colon)
            I = convert($R{Int32},I)
            X = convert($R{$T},X)
            @knet8($("addrows_$F"),(Cint,Cint,Cint,$P{Int},$P{$T},$P{$T}),
                   size(A,1), size(A,2), length(I), I, A, X)
            return A
        end

    end
end


text:
struct Margins
    top::Integer
    bottom::Integer
    left::Integer
    right::Integer
    function Margins(; top, bottom, left, right)
        top = isnothing(top) ? 0 : top
        bottom = isnothing(bottom) ? 0 : bottom
        left = isnothing(left) ? 0 : left
        right = isnothing(right) ? 0 : right
        new(top, bottom, left, right)
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Line <: ACBranch
        name::String
        available::Bool
        activepower_flow::Float64
        reactivepower_flow::Float64
        arc::Arc
        r::Float64
        x::Float64
        b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}
        rate::Float64
        anglelimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        services::Vector{Service}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `activepower_flow::Float64`
- `reactivepower_flow::Float64`
- `arc::Arc`
- `r::Float64`: System per-unit value, validation range: (0, 4), action if invalid: error
- `x::Float64`: System per-unit value, validation range: (0, 4), action if invalid: error
- `b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}`: System per-unit value, validation range: (0, 100), action if invalid: error
- `rate::Float64`
- `anglelimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`, validation range: (-1.571, 1.571), action if invalid: error
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Line <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    "System per-unit value"
    r::Float64
    "System per-unit value"
    x::Float64
    "System per-unit value"
    b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}
    rate::Float64
    anglelimits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Line(name, available, activepower_flow, reactivepower_flow, arc, r, x, b, rate, anglelimits, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    Line(name, available, activepower_flow, reactivepower_flow, arc, r, x, b, rate, anglelimits, services, ext, forecasts, InfrastructureSystemsInternal(), )
end

function Line(; name, available, activepower_flow, reactivepower_flow, arc, r, x, b, rate, anglelimits, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    Line(name, available, activepower_flow, reactivepower_flow, arc, r, x, b, rate, anglelimits, services, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function Line(::Nothing)
    Line(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        b=(from=0.0, to=0.0),
        rate=0.0,
        anglelimits=(min=-1.571, max=1.571),
        services=Device[],
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::Line) = value.name
"""Get Line available."""
get_available(value::Line) = value.available
"""Get Line activepower_flow."""
get_activepower_flow(value::Line) = value.activepower_flow
"""Get Line reactivepower_flow."""
get_reactivepower_flow(value::Line) = value.reactivepower_flow
"""Get Line arc."""
get_arc(value::Line) = value.arc
"""Get Line r."""
get_r(value::Line) = value.r
"""Get Line x."""
get_x(value::Line) = value.x
"""Get Line b."""
get_b(value::Line) = value.b
"""Get Line rate."""
get_rate(value::Line) = value.rate
"""Get Line anglelimits."""
get_anglelimits(value::Line) = value.anglelimits
"""Get Line services."""
get_services(value::Line) = value.services
"""Get Line ext."""
get_ext(value::Line) = value.ext

InfrastructureSystems.get_forecasts(value::Line) = value.forecasts
"""Get Line internal."""
get_internal(value::Line) = value.internal


text:
using HybridSystems, MathematicalSystems

function unstatedep(s::AbstractSwitchedSystem)
    if typeof(s) <: StateDepDiscreteSwitchedLinearSystem
        n = nstates(s.automaton)
        oldmodes = s.modes
        modes = Vector{ContinuousIdentitySystem}(undef, n)
        for (i, mode) in enumerate(oldmodes)
            modes[i] = ContinuousIdentitySystem(mode.statedim)
        end
        G = s.automaton
        oldrmaps = s.resetmaps
        rmaps = Vector{LinearMap}(undef, n)
        for (i, map) in enumerate(oldrmaps)
            rmaps[i] = LinearMap(map.A)
        end
        sw = s.switchings
        s = HybridSystem(G, modes, rmaps, sw)
    end
    s
end

function unconstrained(s::AbstractSwitchedSystem)
    if typeof(s) <: ConstrainedDiscreteSwitchedLinearSystem
        s = HybridSystem(OneStateAutomaton(nstates(s.automaton)), s.modes, s.resetmaps, s.switchings)
    elseif typeof(s) <: StateDepDiscreteSwitchedLinearSystem
        n = nstates(s.automaton)
        G = LightAutomaton(n)
        for σ in 1:n
            for σnext in 1:n
                add_transition!(G, σ, σnext, σ)
            end
        end
        s = HybridSystem(G, s.modes, s.resetmaps, s.switchings)
    end
    s
end


text:
struct Covariate{T<:AbstractFloat} <: AbstractCovariate{T}
    name::String
    data::Vector{T}
end

Base.length(covariate::Covariate{T}) where {T<:AbstractFloat} = length(covariate.data)

function Covariate(data::AbstractVector{T}) where {T<:AbstractFloat}
    Covariate{T}("", data)
end

function Covariate(name::String, data::AbstractVector{T}) where {T<:Integer}
    Covariate{Float32}(name, convert(Vector{Float32}, data))
end

function Covariate(data::AbstractVector{T}) where {T<:Integer}
    Covariate{Float32}("", convert(Vector{Float32}, data))
end

function Covariate{T}(name::String, length::Integer, datpath::String) where {T<:AbstractFloat}
    data = Vector{T}(undef, length)
    open(datpath) do f
        read!(f, data)
    end
    Covariate{T}(name, data)
end

function slice(covariate::Covariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slice(covariate.data, fromobs, toobs, slicelength)
end

function Base.map(covariate::Covariate, dataframe::AbstractDataFrame)
    dataframe[getname(covariate)]
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/libim_image_/include/im_image.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static is_windows()? (const libim_image_ = "im") : (const libim_image_ = "libim")  # Name of IM shared lib.

function imImageCreate(width::Integer, height::Integer, color_space::Integer, data_type::Integer)
  ccall( (:imImageCreate, libim_image_), Ptr{imImage}, (Cint, Cint, Cint, Cint), width, height, color_space, data_type)
end
function imImageInit(width::Integer, height::Integer, color_mode::Integer, data_type::Integer, data_buffer::Ptr{Void}, palette::Ptr{Clong}, palette_count::Integer)
  ccall( (:imImageInit, libim_image_), Ptr{imImage}, (Cint, Cint, Cint, Cint, Ptr{Void}, Ptr{Clong}, Cint), width, height, color_mode, data_type, data_buffer, palette, palette_count)
end
function imImageCreateBased(image::Ptr{imImage}, width::Integer, height::Integer, color_space::Integer, data_type::Integer)
  ccall( (:imImageCreateBased, libim_image_), Ptr{imImage}, (Ptr{imImage}, Cint, Cint, Cint, Cint), image, width, height, color_space, data_type)
end
function imImageDestroy(image::Ptr{imImage})
	ccall( (:imImageDestroy, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageAddAlpha(image::Ptr{imImage})
  ccall( (:imImageAddAlpha, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageSetAlpha(image::Ptr{imImage}, alpha::Cfloat)
  ccall( (:imImageSetAlpha, libim_image_), Void, (Ptr{imImage}, Cfloat), image, alpha)
end
function imImageRemoveAlpha(image::Ptr{imImage})
  ccall( (:imImageRemoveAlpha, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageReshape(image::Ptr{imImage}, width::Integer, height::Integer)
  ccall( (:imImageReshape, libim_image_), Void, (Ptr{imImage}, Cint, Cint), image, width, height)
end
function imImageCopy(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopy, libim_image_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageCopyData(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopyData, libim_image_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageCopyAttributes(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopyAttributes, libim_image_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageMergeAttributes(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageMergeAttributes, libim_image_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageCopyPlane(src_image::Ptr{imImage}, src_plane::Integer, dst_image::Ptr{imImage}, dst_plane::Integer)
  ccall( (:imImageCopyPlane, libim_image_), Void, (Ptr{imImage}, Cint, Ptr{imImage}, Cint), src_image, src_plane, dst_image, dst_plane)
end
function imImageDuplicate(image::Ptr{imImage})
  ccall( (:imImageDuplicate, libim_image_), Ptr{imImage}, (Ptr{imImage},), image)
end
function imImageClone(image::Ptr{imImage})
  ccall( (:imImageClone, libim_image_), Ptr{imImage}, (Ptr{imImage},), image)
end
function imImageSetAttribute(image::Ptr{imImage}, attrib::String, data_type::Integer, count::Integer, data::Ptr{Void})
  ccall( (:imImageSetAttribute, libim_image_), Void, (Ptr{imImage}, Ptr{UInt8}, Cint, Cint, Ptr{Void}), image, attrib, data_type, count, data)
end
function imImageGetAttribute(image::Ptr{imImage}, attrib::String, data_type::Ptr{Int}, count::Ptr{Int})
  ccall( (:imImageGetAttribute, libim_image_), Ptr{Void}, (Ptr{imImage}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}), image, attrib, data_type, count)
end
function imImageGetAttribInteger(image::Ptr{imImage}, attrib::String, index::Integer)
  ccall( (:imImageGetAttribInteger, libim_image_), Ptr{Void}, (Ptr{imImage}, Ptr{UInt8}, Ptr{Cint}), image, attrib, index)
end
function imImageGetAttribReal(image::Ptr{imImage}, attrib::String, index::Integer)
  ccall( (:imImageGetAttribReal, libim_image_), Ptr{Void}, (Ptr{imImage}, Ptr{UInt8}, Ptr{Cint}), image, attrib, index)
end
function imImageGetAttribString(image::Ptr{imImage}, attrib::Ptr{UInt8})
  ccall( (:imImageGetAttribString, libim_image_), Ptr{UInt8}, (Ptr{imImage}, Ptr{UInt8}), image, attrib)
end
function imImageGetAttributeList(image::Ptr{imImage}, attrib::Ptr{Ptr{UInt8}}, attrib_count::Ptr{Int})
  ccall( (:imImageGetAttributeList, libim_image_), Void, (Ptr{imImage}, Ptr{Ptr{UInt8}}, Cint), image, attrib, attrib_count)
end
function imImageClear(image::Ptr{imImage})
  ccall( (:imImageClear, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageIsBitmap(image::Ptr{imImage})
  ccall( (:imImageIsBitmap, libim_image_), Cint, (Ptr{imImage},), image)
end
function imImageSetPalette(image::Ptr{imImage}, palette::Ptr{Int}, palette_count::Integer)
  ccall( (:imImageSetPalette, libim_image_), Void, (Ptr{imImage}, Ptr{Clong}, Cint), image, palette, palette_count)
end
function imImageMatchSize(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchSize, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatchColor(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchColor, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatchDataType(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchDataType, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatchColorSpace(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchColorSpace, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatch(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatch, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageSetMap(image::Ptr{imImage})
  ccall( (:imImageSetMap, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageSetBinary(image::Ptr{imImage})
  ccall( (:imImageSetBinary, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageSetGray(image::Ptr{imImage})
  ccall( (:imImageSetGray, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageMakeBinary(image::Ptr{imImage})
  ccall( (:imImageMakeBinary, libim_image_), Void, (Ptr{imImage},), image)
end
function imImageMakeGray(image::Ptr{imImage})
  ccall( (:imImageMakeGray, libim_image_), Void, (Ptr{imImage},), image)
end
function imFileLoadImage(ifile::Ptr{Int}, index::Integer, error::Ptr{Int})
  ccall( (:imFileLoadImage, libim_image_), Ptr{imImage}, (Ptr{Cint}, Cint, Ptr{Cint}), ifile, index, error)
end
function imFileLoadImageFrame(ifile::Ptr{Int}, index::Integer, image::Ptr{imImage}, error::Ptr{Int})
    ccall( (:imFileLoadImageFrame, libim_image_), Void, (Ptr{Cint}, Cint, Ptr{imImage}, Ptr{Cint}), ifile,
            index, image, error)
end
function imFileLoadBitmap(ifile::Ptr{Int}, index::Integer, error::Ptr{Int})
  ccall( (:imFileLoadBitmap, libim_image_), Ptr{imImage}, (Ptr{Cint}, Cint, Ptr{Cint}), ifile, index, error)
end
function imFileLoadImageRegion(ifile::Ptr{Int}, index::Integer, bitmap::Integer, error::Ptr{Cint}, xmin::Integer,
                               xmax::Integer, ymin::Integer, ymax::Integer, width::Integer, height::Integer)
    ccall( (:imFileLoadImageRegion, libim_image_), Ptr{imImage}, (Ptr{Cint}, Cint, Cint, Ptr{Cint}, Cint, Cint,
          Cint, Cint, Cint, Cint), ifile, index, bitmap, error, xmin, xmax, ymin, ymax, width, height)
end
function imFileLoadBitmapFrame(ifile::Ptr{Int}, index::Integer, image::Ptr{imImage}, error::Ptr{Int})
    ccall( (:imFileLoadBitmapFrame, libim_image_), Void, (Ptr{Cint}, Cint, Ptr{imImage}, Ptr{Cint}), ifile,
            index, image, error)
end
function imFileSaveImage(ifile::Ptr{Int}, image::Ptr{imImage})
  ccall( (:imFileSaveImage, libim_image_), Cint, (Ptr{Cint}, Ptr{imImage}), ifile, image)
end
function imFileImageLoad(file_name::String, index::Integer, error::Ptr{Cint})
  ccall( (:imFileImageLoad, libim_image_), Ptr{imImage}, (Ptr{UInt8}, Cint, Ptr{Cint}), file_name, index, error)
end
function imFileImageLoadBitmap(file_name::String, index::Integer, error::Ptr{Int})
  ccall( (:imFileImageLoadBitmap, libim_image_), Ptr{imImage}, (Ptr{UInt8}, Cint, Ptr{Cint}), file_name, index, error)
end
function imFileImageLoadRegion(file_name::String, index::Integer, bitmap::Integer, error::Ptr{Int}, xmin::Integer,
                               xmax::Integer, ymin::Integer, ymax::Integer, width::Integer, height::Integer)
    ccall( (:imFileImageLoadRegion, libim_image_), Ptr{imImage}, (Ptr{UInt8}, Cint, Cint, Ptr{Cint}, Cint, Cint,
             Cint, Cint, Cint, Cint), file_name, index, bitmap, error, xmin, xmax, ymin, ymax, width, height)
end
function imFileImageSave(file_name::String, format::String, image::Ptr{imImage})
  ccall( (:imFileImageSave, libim_image_), Cint, (Ptr{UInt8}, Ptr{UInt8}, Ptr{imImage}), file_name, format, image)
end


text:
#
println("Cg) Test of the reduced 1- and 2-particle density matrices & natural orbitals.")
#
setDefaults("print summary: open", "zzz-ReducedDensityMatrix.sum")
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=Radial.Grid(true), nuclearModel=Nuclear.Model(26.),
                        configs=[Configuration("[Ne] 3s^2 3p^5")],
                        propertySettings=[ ReducedDensityMatrix.Settings() ] )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
module DiscretisedFluidQueues

import Jacobi, LinearAlgebra, SparseArrays, StaticArrays
import Base: *, +, size, show, getindex, setindex!, sum

# Types
export Model, DiscretisedFluidQueue, BoundedFluidQueue, Phase, PhaseSet # Queues are <:Model
export AbstractMatrixExponential, ConcentratedMatrixExponential, MatrixExponential
export Mesh, DGMesh, FRAPMesh, FVMesh # are <:Mesh
export Generator, FullGenerator, LazyGenerator # are <:Generator
export SFMDistribution, AbstractIntegrationMethod, AutoQuadrature, Quadrature, TrapezoidRule
export ExplicitRungeKuttaScheme, ForwardEuler, Heuns, StableRK3, StableRK4 # Euler, RungeKutta4 <: TimeIntegrationScheme
export Simulation

# Functions
export augment_model, membership, N₋, N₊, n_phases, phases, rates # Model methods
export cell_nodes, Δ, n_bases_per_cell, n_bases_per_phase, n_intervals, total_n_bases # Mesh methods
export interior_point_mass, left_point_mass, right_point_mass, integrate_time # SFMDistribution methods
export simulate, fixed_time, n_jumps, first_exit_x # Simulation methods
export build_lazy_generator, build_full_generator, static_generator
export cme_params, pdf, ccdf, cdf, build_me, cell_probs
export normalised_closing_operator_cdf, normalised_closing_operator_pdf
export naive_normalised_closing_operator_cdf, naive_normalised_closing_operator_pdf
export unnormalised_closing_operator_cdf, unnormalised_closing_operator_pdf
export limit, Limiter, NoLimiter, GeneralisedMUSCL

# model
include("1_SFM.jl")

include("2_abstract_mesh.jl") # things which apply to all meshes
include("2a_discretised_fluid_queue.jl")
include("3_lazy_generators.jl")
include("4_full_generators.jl")

# auxillary functions
include("6_ME_tools.jl") # used in FRAPApproximation.jl
include("7_polynomials.jl") # used in discontinuous_Galerkin.jl

include("8_discontinuous_Galerkin.jl")
include("9_finite_volume_method.jl")
include("10_FRAP_approximation.jl")

include("11_distributions.jl")
include("11c_extend_Base_operations.jl")
include("11d_limiters.jl")

include("12_time_integration.jl")

include("13_simulate.jl")

include("14_error_metrics.jl")

end

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

# SUITE = BenchmarkGroup()
model = "ROBE21"
cases = ["1", "2", "3"]
# SUITE[model] = BenchmarkGroup()

include("robertson.jl")
num_timesteps = Int[]
final_width = Float64[]
elapsed_time = Float64[]

# ----------------------------------------
#  Case 1: alpha = 0.4, beta = 1.0e2, gamma = 1.0e3
# ----------------------------------------
println("Case 1")
const param1 = [0.4, 1.0e2, 1.0e3]

prob = robertson()
alg = TMJets21a(abstol=1e-10, orderT=5, orderQ=1, maxsteps=3500)

# warm-up run
@elapsed sol1 = solve(prob, T=Tf, alg=alg, params=param1);
eltime1 = @elapsed sol1 = solve(prob, T=Tf, alg=alg, params=param1);

# For later use
dom1 = @. getindex(getfield(sol1,:Δt),1);
# Compute the sum by summing the Taylor models, and then evaluating
sum1 = @. evaluate( evaluate( sum(getfield(sol1,:X)),
    domain(getindex(getfield(sol1,:X),1)) ), (normalized_box,))
# Evaluate the Taylor models first, and then sum; yields a much wider result!
# sum1 = [sum(evaluate.(evaluate.(getfield(sol1[ind], :X),
#     domain(getindex(getfield(sol1[ind],:X),1)) ),
#     (normalized_box,))) for ind in eachindex(sol1)];

# Store results
push!(num_timesteps, length(sol1)-1)
push!(final_width, diam(sum1[end]))
push!(elapsed_time, eltime1)
println("Time steps   $(cases[1]) : $(num_timesteps[1])")
println("Final width  $(cases[1]) : $(final_width[1])")
println("Elapsed time $(cases[1]) : $(elapsed_time[1])")

# # benchmark
# SUITE[model][cases[1]] = @benchmarkable solve($prob, T=$Tf, alg=$alg, params=$param1)

# ----------------------------------------
#  Case 2: alpha = 0.4, beta = 1.0e3, gamma = 1.0e5
# ----------------------------------------
println("Case 2")
const param2 = [0.4, 1.0e3, 1.0e5]

# prob = robertson()
alg = TMJets21a(abstol=1e-10, orderT=7, orderQ=1, maxsteps=32000)

# warm-up run
sol2 = solve(prob, T=Tf, alg=alg, params=param2);
eltime2 = @elapsed sol2 = solve(prob, T=Tf, alg=alg, params=param2);

# For later use
dom2 = @. getindex(getfield(sol2,:Δt),1);
# Compute the sum by summing the Taylor models, and then evaluating
sum2 = @. evaluate( evaluate( sum(getfield(sol2,:X)),
    domain(getindex(getfield(sol2,:X),1)) ), (normalized_box,))
# Evaluate the Taylor models first, and then sum; yields a much wider result!
# sum2 = [sum(evaluate.(evaluate.(getfield(sol2[ind], :X),
#     domain(getindex(getfield(sol2[ind],:X),1)) ),
#     (normalized_box,))) for ind in eachindex(sol2)];

# Store results
push!(num_timesteps, length(sol2)-1)
push!(final_width, diam(sum2[end]))
push!(elapsed_time, eltime2)
println("Time steps   $(cases[2]) : $(num_timesteps[2])")
println("Final width  $(cases[2]) : $(final_width[2])")
println("Elapsed time $(cases[2]) : $(elapsed_time[2])")

# # benchmark
# SUITE[model][cases[2]] = @benchmarkable solve($prob, T=$Tf, alg=$alg, params=$param2)

# ----------------------------------------
#  Case 3: alpha = 0.4, beta = 1.0e3, gamma = 1.0e7
# ----------------------------------------
println("Case 3")
const param3 = [0.4, 1.0e3, 1.0e7]

# prob = robertson()
alg = TMJets21a(abstol=1e-10, orderT=10, orderQ=1, maxsteps=75000)

# warm-up run
eltime3 = @elapsed sol3 = solve(prob, T=Tf, alg=alg, params=param3);

# For later use
dom3 = @. getindex(getfield(sol3,:Δt),1);
# Compute the sum by summing the Taylor models, and then evaluating
sum3 = @. evaluate( evaluate( sum(getfield(sol3,:X)),
    domain(getindex(getfield(sol3,:X),1)) ), (normalized_box,))
# Evaluate the Taylor models first, and then sum; yields a much wider result!
# sum3 = [sum(evaluate.(evaluate.(getfield(sol3[ind], :X),
#     domain(getindex(getfield(sol3[ind],:X),1)) ),
#     (normalized_box,))) for ind in eachindex(sol3)];

# Store results
push!(num_timesteps, length(sol3)-1)
push!(final_width, diam(sum3[end]))
push!(elapsed_time, eltime3)
println("Time steps   $(cases[3]) : $(num_timesteps[3])")
println("Final width  $(cases[3]) : $(final_width[3])")
println("Elapsed time $(cases[3]) : $(elapsed_time[3])")

# # benchmark
# SUITE[model][cases[3]] = @benchmarkable solve($prob, T=$Tf, alg=$alg, params=$param3)

# # ==============================================================================
# # Execute benchmarks and save benchmark results
# # ==============================================================================

# # tune parameters
# tune!(SUITE)

# # run the benchmarks
# results = run(SUITE, verbose=true)

# # return the sample with the smallest time value in each test
# println("minimum time for each benchmark:\n", minimum(results))

# # return the median for each test
# println("median time for each benchmark:\n", median(results))

# # export runtimes
# runtimes = Dict()
# for (i, c) in enumerate(cases)
#      t = median(results[model][c]).time * 1e-9
#      runtimes[c] = t
# end

validation = ""
for (i, c) in enumerate(cases)
    print(io, "JuliaReach, $model, $c, $validation, $(elapsed_time[i]), " *
              "$(num_timesteps[i]), $(final_width[i])\n")
end

# ==============================================================================
# Plot
# ==============================================================================

fig = Plots.plot(legend=:topright, legendfontsize=18,
    tickfont=font(30, "Times"), guidefontsize=45,
    xlab=L"t", # \raisebox{-0.5mm}{\textcolor{white}{.}}",
    ylab=L"s", # \raisebox{2mm}{\textcolor{white}{.}}",
    xtick=0:10:40, ytick=0.999:0.0002:1.001,
    xlims=(0., 42.), ylims=(0.999, 1.00105),
    # bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
    size=(1000, 1000))

Plots.plot!(fig, IntervalBox.(dom1,sum1),
    color=:blue,  linecolor=:match, alpha=0.8, label="Case 1")
Plots.plot!(fig, IntervalBox.(dom2,sum2),
    color=:green, linecolor=:match, alpha=0.8, label="Case 2")
Plots.plot!(fig, IntervalBox.(dom3,sum3),
    color=:red,   linecolor=:match, alpha=0.8, label="Case 3")

savefig(fig, "ARCH-COMP21-JuliaReach-Robertson.png")


text:
f = open("test.txt")
input = [5.0;3.0;-1.0;2.5;-0.5;0.3;0.7;1.9;3.5;-0.3]
hiddenSizes=[25,9];
action=1;
qtarget = 0.5;

weights = readlines(f);
counter = 1;
W = Any[]
x = length(input)
for h in hiddenSizes
    params = split(weights[counter],",")
    w = parse.(Float64,params)
    w = reshape(w,h,x)
    push!(W,w)
    counter += 1
    params = split(weights[counter],",")
    b = parse.(Float64,params)
    b = reshape(b,h,1)
    push!(W,b)
    counter += 1
    x = h
end

function forward(w,x)
    for i=1:2:length(w)
        x = w[i]*x .+ w[i+1]
        if i<length(w)-1
            x = tanh.(x)
        end
    end
    return x
end

function loss(w,x,action,target)
    x3 = forward(w,x);
    loss = x3[action]-target
    return loss^2
end

using Knet

gradloss = grad(loss);

gradloss(W,input,1,qtarget)

max_stars_repo_name: ekinakyurek/Deep-RL-Cars

max_issues_repo_name: ekinakyurek/Deep-RL-Cars

max_forks_repo_name: ekinakyurek/Deep-RL-Cars


text:
export findsitemap
export findorbitalmap

"""
    findsitemap(unitcell, translation_operation)

Find which site gets mapped to which site with what lattice displacement.
Return a vector of tuples `(j, Rj)` at index `i`, representing that the site `i`
gets mapped to site `j` with lattice displacement `Rj`.

Since periodic boundary condition is always assumed, the return value always exists.

This function is used to find the mapping between sites under translation symmetry in a
given lattice. The overall transformation will be a combiation of (1) the result here, and
(2) Bravais lattice transformation. Since all the integer translations will be
taken care of by the Bravais lattice transformation, `Rj` here is all zero vectors.
"""
function findsitemap(
    unitcell::UnitCell,
    tsym_op::TranslationOperation{<:Integer},
)::Vector{Tuple{Int, Vector{Int}}}
    norb = numsite(unitcell)
    dim = dimension(unitcell)
    return [(i, zeros(Int, dim)) for i in 1:norb]
end

"""
    findsitemap(unitcell, translation_symmetry)

Find the sitemap for every element of the translation symmetry.
Return `[findsitemap(unitcell, m) for m in elements(translation_symmetry)]`.
"""
function findsitemap(unitcell::UnitCell, tsym::FiniteTranslationSymmetry)
    return [findsitemap(unitcell, m) for m in elements(tsym)]
end

"""
    findsitemap(unitcell, point_operation)

Find which site gets mapped to which site with what lattice displacement.
Return a vector of tuples `(j, Rj)` at index `i`, representing that the site `i`
gets mapped to site `j` with lattice displacement `Rj`.
If the unit cell is not compatible with the point operation, return `nothing`.
"""
function findsitemap(
    unitcell::UnitCell,
    psym_op::PointOperation{<:Integer},
)::Union{Nothing, Vector{Tuple{Int, Vector{Int}}}}
    norb = numsite(unitcell)
    map = Tuple{Int, Vector{Int}}[]
    for (orbname, orbfc) in unitcell.sites
        j, Rj = findsiteindex(unitcell, psym_op.matrix * orbfc)
        j <= 0 && return nothing  # throw(ArgumentError("site map not found with $unitcell and $psym_op"))
        push!(map, (j, Rj))
    end
    return map
end

"""
    findsitemap(unitcell, point_symmetry)

Find the sitemap for every element of the point symmetry.
Return `nothing` if the unit cell is incompatible with any element of the point symmetry.
"""
function findsitemap(
    unitcell::UnitCell,
    psym::PointSymmetry,
)::Union{Nothing, Vector{Vector{Tuple{Int, Vector{Int}}}}}
    out = Vector{Tuple{Int, Vector{Int}}}[]
    sizehint!(out, length(elements(psym)))
    for el in elements(psym)
        m = findsitemap(unitcell, el)
        isnothing(m) && return nothing
        push!(out, m)
    end
    return out
end


text:
# module Lupin

function Base.reverse(nt::NamedTuple{names, T}) where {names, T <: Tuple}
    NamedTuple{reverse(names)}((reverse ∘ values)(nt))
end

function Base.reverse(dict::Dict{K, V})::Dict{V, K} where {K, V}
    d = Dict{V, K}()
    for (k, v) in dict
        d[v] = k
    end
    d
end

# module Lupin


text:
"""
    vector_array_flatten(x::Vector{Array})

Convert a vector of `N`-dimensional arrays to an `N+1`-dimensional array by concatenating along a new leading dimension.
"""
function vector_array_flatten(x::Vector{Array{T,N}}) where {T,N}
    cat(itemize.(x)..., dims=1)
end

"""
    itemize(x::Array{T,N})

Convert an `N`-dimensional array to an `N+1`-dimension array of the same size with a singleton leading dimension.
"""
function itemize(x::Array{T,N}) where {T,N}
    y = Array{T,N+1}(undef, 1, size(x)...)
    for i in 1:length(x)
        y[i] = x[i]
    end
    return y
end

"""
    vector_array_transpose(x, dim=1)

"""
function vector_array_transpose(x::Vector{Array{T,N}}, dim=1) where {T,N}

end


text:
module MT1D
export calc_response, calc_rms

"""
Description
===========

Forward model the apparent resistivity and phase response for a model.

Arguments
=========

- `model::Matrix`: Matrix representing the model. Each row is a layer. Columns are:
    1. Depth to top of layer (m)
    2. Resistivity of layer (Ωm)
- `fs::Vector`: Frequencies to forward model the response at.

Returns
=======

- `Tuple{ρ::Vector,Φ::Vector}`: Forward modelled apparent resistivity and phase.
"""
function calc_response{T<:AbstractFloat}(model::Matrix{T}, fs::Vector{T})
    const μ = 4E-7 * π
    nf = length(fs)
    nlayers = length(model[:,1])

    ρa = zeros(nf)
    Φ = zeros(nf)

    for fidx in 1:nf
        k = sqrt(-im * 2π * fs[fidx] * μ / model[end, 2])
        C = 1 / k
        for n in nlayers-1:-1:1
            k = sqrt(-im * 2π * fs[fidx] * μ / model[n, 2])
            Δz = model[n+1, 1] - model[n, 1]
            C = (C * k + tanh(k * Δz)) / (C * k^2 * tanh(k * Δz) + k)
        end
        Z = 2π * fs[fidx] * μ * C
        ρa[fidx] = abs(Z)^2 / (2π * fs[fidx] * μ)
        Φ[fidx] = 90 - angle(Z) * 180 / π
    end

    return (ρa, Φ);
end

"""
Description
===========

Calculate the root mean squared misfit between a set of observations and
modelled responses.

Arguments
=========

- `obs::Vector`: Observation points
- `model::Vector`: Modelled responses
- `ϵ::Vector`: Uncertainties for `obs`

Returns
=======

- `r::Real`: Calculated RMS misfit
"""
function calc_rms(obs::Vector, model::Vector, ϵ::Vector)
    sqrt(sum(((obs - model) ./ ϵ).^2) / (2 * length(obs)))
end
end


text:
"""
    AbstractOperator

An abstract type for linear operators on Banach spaces.
"""
abstract type AbstractOperator end

"""
    KernelOperator{K,S} <: AbstractOperator

A type for representing kernel integral operators over Hilbert spaces.
"""
struct KernelOperator{K,S} <: AbstractOperator
    kernel::K
    solver::S
end

KernelOperator(k) = KernelOperator(k,HCubatureJL())

function eval_(ko::KernelOperator, f; solver=ko.solver)
    k = ko.kernel
    fun = function(u,p)
        k(p,u) * f(u)
    end
    return function(x)
        prob = QuadratureProblem{false}(fun,0,1,x)
        solve(prob, solver).u
    end
end

function (ko::KernelOperator)(x)
    eval_(ko, x)
end

Base.:*(ko::KernelOperator, x) = ko(x)

function ChainRulesCore.rrule(ko::KernelOperator, x)
    y = ko(x)
    function ko_pullback(δy)
        δko = NoTangent()
        δx = t -> kernel_transpose_integral(ko.kernel, δy, t, ko.solver)
        return δko, δx
    end
    return y, ko_pullback
end

"""
    FourierKernel{T<:Real}

An object representing a Fourier kernel function ``k(x,y) = sin (ω  π xy + θ)``
"""
struct FourierKernel{T<:Real}
    ω::T
    θ::T
end

FourierKernel(; a = 0.0, b = 1.0) = FourierKernel((b-a), a)

(ff::FourierKernel)(x,y) = sin(ff.ω*π*x*y + ff.θ)

"""
    PolynomialKernel{M<:AbstractMatrix}

A type representing a polynomial kernel function.
"""
struct PolynomialKernel{M<:AbstractMatrix}
    m::M
end

function (pk::PolynomialKernel)(x,y)
    (m,n) = size(pk.m)
    X = [x^(k-1) for k in 1:m]
    Y = [y^(k-1) for k in 1:n]
    X * pk.m * Y
end


text:
using Seaborn
using Base.Test


text:
@inline execute(ctx::Any, f, args...) = overdub_recurse(Val(false), ctx, f, args...)
@inline execute_inbounds(ctx::Any, f, args...) = overdub_recurse(Val(true), ctx, f, args...)

#@generated function overdub_recurse(inbounds, ctx, f, args...)
function overdub_recurse_gen(self, inbounds, ctx, f, args)
    # configure the global logger to use plain stderr so that we can log without task switches
    old_logger = global_logger()
    global_logger(Logging.ConsoleLogger(Core.stderr))
    @debug "Overdubbing function call" func=f types=args inbounds=inbounds

    # don't recurse into Core
    if parentmodule(f) == Core || f == typeof(Base.getindex)
        @debug "Can't overdub Core function"
        global_logger(old_logger)
        if inbounds <: Val{true}
            return :(@inbounds f(args...))
        else
            return :(f(args...))
        end
    end

    # refuse to overdub already-overdubbed functions
    if f <: typeof(execute)
        error("can't double-overdub")
    end

    # slot number flag indicating that it is already set
    # yes this is hacky :o
    mark_fixed = Int64(1 << 62)

    # retrieve code
    # NOTE: this could use code_lowered if it weren't for F being a function type
    world = typemax(UInt)
    ## initial Method
    matched_methods = Base._methods_by_ftype(Tuple{f,args...}, -1, world)
    length(matched_methods) == 1 || error("did not uniquely match method")
    type_signature, raw_static_params, method = first(matched_methods)
    if method.pure
        @debug "Won't overdub pure function"
        global_logger(old_logger)
        return :(f(args...))
    end

    ## initial CodeInfo
    method_instance = Core.Compiler.code_for_method(method, type_signature, raw_static_params, world, false)
    method_signature = method.sig
    static_params = Any[raw_static_params...]
    code_info = Core.Compiler.retrieve_code_info(method_instance)
    isa(code_info, Core.CodeInfo) || error("could not retrieve original code")

    # prepare for rewriting
    body = Expr(:block)
    original_code_info = code_info
    code_info = Core.Compiler.copy_code_info(code_info)
    body.args = code_info.code
    insert_point = findfirst(item->isa(item,Expr), body.args) # past meta, lineno, ...

    # NOTE: most of the complexity here comes from the disconnect between (the signature of)
    #       the called generated function, and the function whose code we're injecting.
    #       arguably, the Julia compiler should handle this (see Cassette.jl#7).

    # substitute static parameters (the called generator doesn't have any)
    Core.Compiler.substitute!(body, 0, Any[], method_signature, static_params, 0, :propagate)

    # rewrite function calls
    # these slotnumbers will not be bumped
    inbounds_slot = Core.SlotNumber(2 | mark_fixed)
    context_slot = Core.SlotNumber(3 | mark_fixed)

    worklist = Any[map(item->(item,item), body.args)...] # item & pos to insert before
    inbounds_stack = Bool[]
    global_ref_ssa_values = Dict{Core.SSAValue, GlobalRef}()
    while !isempty(worklist)
        item, paren = popfirst!(worklist)

        if Meta.isexpr(item, :inbounds)
            if item.args[1] == true
                push!(inbounds_stack, true)
            elseif item.args[1] == :pop
                pop!(inbounds_stack)
            end
            continue
        end

        if isa(item, Expr)
            # DFS expr arguments
            prepend!(worklist, map(item->(item,paren), item.args))
        end

        if Meta.isexpr(item, :(=))
            lhs, rhs = item.args
            if isa(lhs, Core.SSAValue) && isa(rhs, GlobalRef)
                global_ref_ssa_values[lhs] = rhs
            end
        end

        if Meta.isexpr(item, :call)
            orig_func = item.args[1]

            # TODO: Fix me, this is hacky
            check_func = orig_func
            if isa(check_func, Core.SSAValue) && haskey(global_ref_ssa_values, check_func)
                check_func = global_ref_ssa_values[check_func]
            end

            if (isa(check_func, GlobalRef) && check_func.mod == Core) ||
               (isa(check_func, GlobalRef) && check_func.mod == Core.Intrinsics) ||
               (isa(check_func, GlobalRef) && check_func.name == :Val) ||
               (isa(check_func, GlobalRef) && check_func.name == :getfield) ||
               (isa(check_func, GlobalRef) && check_func.name == :getproperty) ||
               (isa(check_func, GlobalRef) && check_func.name == :not_int) ||
               check_func == GlobalRef(Main, :eltype) ||
               check_func == GlobalRef(Main, :(:)) ||
               check_func == GlobalRef(Base, :iterate) ||
               check_func == GlobalRef(Base, :to_indices)
                continue
            end

            if !isempty(inbounds_stack) || inbounds <: Val{true}
                item.args[1] = GlobalRef(@__MODULE__, :execute_inbounds)
            else
                item.args[1] = GlobalRef(@__MODULE__, :execute)
            end
            insert!(item.args, 2, context_slot)
            insert!(item.args, 3, orig_func)
        end
    end

    # destructure the splatted argument tuple
    argc = length(args)
    paramc = method.nargs - 1
    splat = Core.SlotNumber(5) # this won't be bumped later on
    ## fix up codeinfo arrays
    code_info.slotnames = Any[code_info.slotnames[1], Symbol("#inbounds#"), Symbol("#ctx#"), Symbol("#f#"), Symbol("#args#"), code_info.slotnames[2:end]...]
    code_info.slotflags = Any[code_info.slotflags[1], 0x00                , 0x00           , 0x00         , 0x00            , code_info.slotflags[2:end]...]
    ## generate new slots
    prelude = Expr[]
    for i in 1:paramc
        # insert new slot
        slotnum = i+5
        slot = Core.SlotNumber(slotnum)
        code_info.slotflags[slotnum] |= 0x01 << 0x01    # mark the slot as assigned to

        # populate it with the actual argument
        arg = Expr(:call, GlobalRef(Core, :getfield), splat, i)
        push!(prelude, :($slot = $arg))
    end
    ## fix uses of slots
    function replace_nodes!(f, code)
        for (i,node) in enumerate(code)
            replacement = f(node)
            if replacement !== nothing
                code[i] = replacement
            elseif isa(node, Expr)
                # visit expr arguments
                replace_nodes!(f, node.args)
            end
        end
    end
    replace_nodes!(body.args) do node
        if isa(node, Core.SlotNumber) && node.id == 1
            return Core.SlotNumber(4)
        elseif isa(node, Core.SlotNumber) && (node.id & mark_fixed) == mark_fixed
            return Core.SlotNumber(node.id & ~mark_fixed)
        elseif isa(node, Core.SlotNumber)
            return Core.SlotNumber(node.id+4)
        elseif isa(node, Core.NewvarNode) && node.slot.id > 1
            return Core.NewvarNode(Core.SlotNumber(node.slot.id+4))
        end
    end
    ## special handling for vararg parameters
    if method.isva
        # the previous final slot assignment is wrong
        isempty(prelude) || pop!(prelude)
        # instead create and assign a tuple containing all trailing arguments
        vararg = Expr(:call, GlobalRef(Core, :tuple))
        for i in paramc:argc
            ssa = Core.SSAValue(code_info.ssavaluetypes)
            arg = Expr(:call, GlobalRef(Core, :getfield), splat, i)
            push!(prelude, :($ssa = $arg))
            push!(vararg.args, ssa)
            code_info.ssavaluetypes += 1
        end
        push!(prelude, :($(Core.SlotNumber(paramc + 5)) = $vararg))
    end
    ## insert slot definitions
    for expr in reverse(prelude)
        insert!(body.args, insert_point, expr)
    end

    # fix labels and references to them
    changes = Dict{Int,Int}()
    for (i, stmnt) in enumerate(code_info.code)
        if isa(stmnt, Core.LabelNode)
            code_info.code[i] = Core.LabelNode(i)
            changes[stmnt.label] = i
        end
    end
    for (i, stmnt) in enumerate(code_info.code)
        if isa(stmnt, Core.GotoNode)
            code_info.code[i] = Core.GotoNode(get(changes, stmnt.label, stmnt.label))
        elseif Meta.isexpr(stmnt, :enter)
            stmnt.args[1] = get(changes, stmnt.args[1], stmnt.args[1])
        elseif Meta.isexpr(stmnt, :gotoifnot)
            stmnt.args[2] = get(changes, stmnt.args[2], stmnt.args[2])
        end
    end

    # set inference limit heursitics
    code_info.method_for_inference_limit_heuristics = method
    code_info.inlineable = true

    # validate
    errors = Core.Compiler.validate_code(method_instance, code_info)
    for e in errors
        @error "Encountered invalid code" code=body.args error=e
    end

    @debug "Rewriting code" original=original_code_info overdubbed=code_info
    global_logger(old_logger)
    return code_info
end

@eval function overdub_recurse(inbounds, ctx, f, args...)
    # manual construction of the generated function in order to control the expand_early arg
    $(begin
          stub = Expr(:new,
                      Core.GeneratedFunctionStub,
                      :overdub_recurse_gen,
                      Any[:overdub_recurse, :inbounds, :ctx, :f, :args],
                      Any[],
                      @__LINE__,
                      QuoteNode(Symbol(@__FILE__)),
                      true)
          Expr(:meta, :generated, stub)
      end)
end