text:
"""
```
dsgevar_likelihood(m::AbstractDSGEVARModel{S}, data::Matrix{S};
    apply_altpolicy::Bool = false) where {S<:Real}
```
evaluates the likelihood of a DSGE-VAR. It is assumed that
`get_λ(m)` retrieves the prior weight `λ` on the DSGE.
"""

# data is assumed an nobs x T+lags matrix,
# where the lags indicate we cut off the data for presampling
function dsgevar_likelihood(m::AbstractDSGEVARModel{S}, data::Matrix{S};
                            apply_altpolicy::Bool = false,
                            use_intercept::Bool = true) where {S<:Real}

    # Set up. Default to non-regime switching if these settings don't exist.
    lags = n_lags(m)
    n_regimes        = haskey(get_settings(m), :n_regimes) ?
        get_setting(m, :n_regimes) : 1
    regime_switching = haskey(get_settings(m), :regime_switching) && n_regimes > 1 ?
        get_setting(m, :regime_switching) : false # if n_regimes == 1, then no switching needed

    # Get population moments
    YYYY, XXYY, XXXX = compute_var_population_moments(data, lags; use_intercept = use_intercept)

    if regime_switching
        error("Regime switching has not been implemented yet.")
        yyyyd = Vector{Matrix{S}}(undef, n_regimes)
        xxyyd = Vector{Matrix{S}}(undef, n_regimes)
        xxxxd = Vector{Matrix{S}}(undef, n_regimes)

        out = compute_system(m; applya_altpolicy = apply_altpolicy,
                             get_population_moments = true, use_intercept = use_intercept)
        for reg = 1:n_regimes
            yyyyd[reg], xxyyd[reg], xxxxd[reg] = out[reg]
        end
    else
        yyyyd, xxyyd, xxxxd = compute_system(m; apply_altpolicy = apply_altpolicy,
                                             get_population_moments = true,
                                             use_intercept = use_intercept)
    end

    return dsgevar_likelihood(YYYY, XXYY, XXXX, yyyyd, xxyyd, xxxxd,
                      size(data, 2) - lags, get_λ(m))
end

"""
```
dsgevar_likelihood(YYYY::Matrix{S}, XXYY::Matrix{S},
    XXXX::Matrix{S}, YYYYD::Matrix{S}, XXYYD::Matrix{S},
    XXXXD::Matrix{S}, T::Int, λ::S,
    lags::Int, n_obs::Int) where {S<:Real}
```
evaluates the likelihood of a DSGE-VAR given the
population moments of the raw data (`YYYY`, `XXYY`, `XXXX`)
and the population moments implied by the DSGE prior (`YYYYD`, `XXYYD`, `XXXXD`).

### Other Inputs
* `T`: number of time periods in data
* `λ`: weight placed on the DSGE prior
* `lags`: number of lags in the VAR
* `n_obs`: number of observables (e.g. number of time series).
"""
function dsgevar_likelihood(YYYY::Matrix{S}, XXYY::Matrix{S},
                            XXXX::Matrix{S}, YYYYD::Matrix{S}, XXYYD::Matrix{S},
                            XXXXD::Matrix{S}, T::Int, λ::S) where {S<:Real}

    if isinf(λ)
        β = XXXXD \ XXYYD
        Σ = YYYYD - XXYYD' * β
        n_obs  = size(YYYY, 1)
        loglik = -(n_obs * T / 2) * log(2 * π) - (T / 2) * log(det(Σ)) -
            .5 * sum(diag(Σ \ (YYYY - 2 * β' * XXYY + β' * XXXX * β)))
        loglik = real(loglik)
    elseif λ < 0
        throw(DomainError("λ", "λ cannot be negative."))
    else
        # Dimensions
        λT     = λ * T
        n_obsc = T + λT
        k      = size(XXXX, 1) # 1 + lags * n_obs
        n_obs  = size(YYYY, 1)

        # Compute weighted population moments
        YYYYC = YYYY + λT * YYYYD
        XXYYC = XXYY + λT * XXYYD
        XXXXC = XXXX + λT * XXXXD

        # Compute OLS estimates of VAR
        Σ = (YYYYC - XXYYC' * (XXXXC \ XXYYC)) / n_obsc
        Σd = YYYYD - XXYYD' * (XXXXD \ XXYYD)

        # Add terms in normalization constant independent of θ
        loglik = (-n_obs * T / 2) * log(2 * π) + (T * n_obs / 2) * log(2)
        for i = 1:n_obs
            loglik += loggamma((n_obsc - k + 1 - i) / 2) - loggamma((λT - k + 1 - i) / 2)
        end

        # Calculate remainder of likelihood
        # p(Y | θ, λT) = ∫ p(Y | β, Σ) p(β, Σ | θ, λT) d(β, Σ), where
        # - p(Y | β, Σ) is the standard Gaussian likelihood of a VAR unconditional on θ,
        #               i.e. the likelihood ignoring the restrictions a DSGE places on a VAR,
        #                    or the likelihood when doing OLS on just the data.
        # - p(β, Σ | θ, λT) is the prior on (β, Σ) implied by a DSGE-VAR
        #                  with weighting λT.
        loglik += -((n_obsc - k) / 2) * log(det(n_obsc * Σ)) + ((λT - k) / 2) * log(det(λT * Σd)) -
            (n_obs / 2) * log(det(XXXXC) / det(λT * XXXXD))
        loglik = real(loglik)

        # Note that the posterior for DSGE parameters and weight λ is
        # p(θ | Y, λ) ∝ p(Y | θ, λ) p(θ | λ)
        # and joint posterior for VAR parameters is
        # p(β, Σ, θ | Y, λ) = p(β, Σ | θ, Y, λ) p(θ | Y, λ), where
        # p(β, Σ | θ, Y, λ) = p(Y, β, Σ | θ, λ) p(β, Σ | θ, λ).
        # Observe that p(β, Σ | θ, λ) is the prior on (β, Σ) implied by  DSGE-VAR
        #                             with weighting λ, and
        #              p(Y, β, Σ | θ, λ) = p(Y | β, Σ) * (p(β, Σ | θ, λ) * Jeffrey's prior),
        #                                and the right two terms form the Minnesota prior.
    end

    return loglik
end


text:
#!/usr/bin/env julia

using Base.Test
using Distributions
using ForwardDiff
using StaticArrays

using Celeste: Model, DeterministicVI

import Celeste: Infer, DeterministicVI, ParallelRun
import Celeste: PSF, SDSSIO, SensitiveFloats, Transform
import Celeste.SensitiveFloats.clear!
import Celeste.SDSSIO: RunCamcolField

include(joinpath(Pkg.dir("Celeste"), "test", "SampleData.jl"))

using SampleData

anyerrors = false

wd = pwd()
# Ensure that test images are available.
const datadir = joinpath(Pkg.dir("Celeste"), "test", "data")
cd(datadir)
run(`make`)
run(`make RUN=4263 CAMCOL=5 FIELD=119`)
# Ensure GalSim test images are available.
const galsim_benchmark_dir = joinpath(Pkg.dir("Celeste"), "benchmark", "galsim")
cd(galsim_benchmark_dir)
run(`make fetch`)
cd(wd)

# Check whether to run time-consuming tests.
long_running_flag = "--long-running"
test_long_running = long_running_flag in ARGS
test_files = setdiff(ARGS, [ long_running_flag ])

if length(test_files) > 0
    testfiles = ["test_$(arg).jl" for arg in test_files]
else
    testdir = joinpath(Pkg.dir("Celeste"), "test")
    testfiles = filter(r"test_.*\.jl", readdir(testdir))
end

if !test_long_running
    warn("Skipping long running tests.  ",
         "To test everything, run tests with the flag ", long_running_flag)
end

for testfile in testfiles
    try
        println("Running ", testfile)
        @time include(testfile)
        println("\t\033[1m\033[32mPASSED\033[0m: $(testfile)")
    catch e
        anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(testfile)")
        rethrow()  # Fail fast.
    end
end

ext: jl


text:
module Logger
# Follow approach taken by https://github.com/CliMA/Oceananigans.jl logger.jl

export QXLogger, QXLoggerMPIPerRank, QXLoggerMPIShared

using Logging
using Dates
using MPI
using TimerOutputs
import UUIDs: UUID, uuid4

import Logging: shouldlog, min_enabled_level, catch_exceptions, handle_message

const PerfLogger = Logging.LogLevel(-125)

struct QXLogger <: Logging.AbstractLogger
    stream::IO
    min_level::Logging.LogLevel
    message_limits::Dict{Any,Int}
    show_info_source::Bool
    session_id::UUID
end

struct QXLoggerMPIShared <: Logging.AbstractLogger
    stream::Union{MPI.FileHandle, Nothing}
    min_level::Logging.LogLevel
    message_limits::Dict{Any,Int}
    show_info_source::Bool
    session_id::UUID
    comm::MPI.Comm
end

struct QXLoggerMPIPerRank <: Logging.AbstractLogger
    stream::Union{MPI.FileHandle, Nothing}
    min_level::Logging.LogLevel
    message_limits::Dict{Any,Int}
    show_info_source::Bool
    session_id::UUID
    comm::MPI.Comm
    root_path::String
end

"""
    QXLogger(stream::IO=stdout, level=Logging.Info; show_info_source=false)

Single-process logger for QXContexts.
"""
function QXLogger(stream::IO=stdout, level=Logging.Info; show_info_source=false)
    return QXLogger(stream, level, Dict{Any,Int}(), show_info_source, uuid4())
end

"""
    QXLoggerMPIShared(stream=nothing,
                      level=Logging.Info;
                      show_info_source=false,
                      comm=MPI.COMM_WORLD,
                      path::String=".")

MPI-IO enabled logger that outputs to a single shared file for all ranks.
"""
function QXLoggerMPIShared(stream=nothing,
                           level=Logging.Info;
                           show_info_source=false,
                           comm=MPI.COMM_WORLD,
                           path::String=".")
    if MPI.Initialized()
        if MPI.Comm_rank(comm) == 0
            log_uid = uuid4()
        else
            log_uid = nothing
        end
        log_uid = MPI.bcast(log_uid, 0, comm)
        f_stream = MPI.File.open(comm, joinpath(path, "QXContexts_io_shared_$(log_uid).log"), read=true,  write=true, create=true)
    else
        error("""MPI is required for this logger. Pleasure ensure MPI is initialised. Use `QXLogger` for non-distributed logging""")
    end
    return QXLoggerMPIShared(f_stream, level, Dict{Any,Int}(), show_info_source, log_uid, comm)
end

"""
    QXLoggerMPIPerRank(stream=nothing,
                       level=Logging.Info;
                       show_info_source=false,
                       comm=MPI.COMM_WORLD,
                       path::String=".")

MPI-friendly logger that outputs to a new file per rank. Creates a UUIDs.uuid4 labelled directory and a per-rank log-file
"""
function QXLoggerMPIPerRank(stream=nothing,
                            level=Logging.Info;
                            show_info_source=false,
                            comm=MPI.COMM_WORLD,
                            path::String=".")
    if MPI.Initialized()
        if MPI.Comm_rank(comm) == 0
            log_uid = uuid4()
        else
            log_uid = nothing
        end
        log_uid = MPI.bcast(log_uid, 0, comm)
    else
        throw("""MPI is required for this logger. Pleasure ensure MPI is initialised. Use `QXLogger` for non-distributed logging""")
    end
    return QXLoggerMPIPerRank(stream, level, Dict{Any,Int}(), show_info_source, log_uid, comm, path)
end

shouldlog(logger::QXLogger, level, _module, group, id) = get(logger.message_limits, id, 1) > 0
shouldlog(logger::QXLoggerMPIShared, level, _module, group, id) = get(logger.message_limits, id, 1) > 0
shouldlog(logger::QXLoggerMPIPerRank, level, _module, group, id) = get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::QXLogger) = logger.min_level
min_enabled_level(logger::QXLoggerMPIShared) = logger.min_level
min_enabled_level(logger::QXLoggerMPIPerRank) = logger.min_level

catch_exceptions(logger::QXLogger) = false
catch_exceptions(logger::QXLoggerMPIShared) = false
catch_exceptions(logger::QXLoggerMPIPerRank) = false

function level_to_string(level)
    level == Logging.Error && return "ERROR"
    level == Logging.Warn  && return "WARN"
    level == Logging.Info  && return "INFO"
    level == Logging.Debug && return "DEBUG"
    level == PerfLogger && return "PERF"
    return string(level)
end

macro perf(expression)
    if haskey(ENV, "QXRUN_TIMER")
        ex = repr(expression)
        return :(t = @elapsed res = $expression; Logging.@logmsg(PerfLogger, (t, $ex)); res)
    else
        return esc(expression)
    end
end

"""
    stamp_builder(rank::Int)

Builds logger timestamp with rank and hostname capture.
Rank defaults to 0 for single-process evaluations
"""
function stamp_builder(rank::Int)
    io = IOBuffer()
    write(io, Dates.format(Dates.now(), "[yyyy/mm/dd-HH:MM:SS.sss]"))
    write(io, "[rank="*string(rank)*"]")
    write(io, "[host="*gethostname()*"]")
    s = String(take!(io))
    close(io)
    return s
end

function handle_message(logger::Union{QXLogger, QXLoggerMPIShared}, level, message, _module, group, id,
                                filepath, line; maxlog = nothing, kwargs...)
    if !isnothing(maxlog) && maxlog isa Int
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return nothing
    end

    buf = IOBuffer()
    level_name = level_to_string(level)
    if MPI.Initialized() && MPI.Comm_rank(MPI.COMM_WORLD) > 1
        rank = MPI.Comm_rank(logger.comm)
    else
        rank = 0
    end

    module_name = something(_module, "nothing")
    msg_timestamp = stamp_builder(rank)

    formatted_message = "$(msg_timestamp) $(level_name) $message"
    if logger.show_info_source || level != Logging.Info
        formatted_message *= " -@-> $(filepath):$(line)"
    end
    formatted_message *= "\n"

    if typeof(logger.stream) <: IO
        write(logger.stream,  formatted_message)
    else
        MPI.File.write_shared(logger.stream, formatted_message)
    end

    return nothing
end

function handle_message(logger::QXLoggerMPIPerRank, level, message, _module, group, id,
    filepath, line; maxlog = nothing, kwargs...)

    if !isnothing(maxlog) && maxlog isa Int
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return nothing
    end

    if !isdir(joinpath(logger.root_path, "QXContexts_io_" * string(logger.session_id)))
        mkdir(joinpath(logger.root_path, "QXContexts_io_" * string(logger.session_id)))
    end

    buf = IOBuffer()
    rank = MPI.Comm_rank(logger.comm)
    level_name = level_to_string(level)

    log_path = joinpath(logger.root_path, "QXContexts_io_" * string(global_logger().session_id), "rank_$(rank).log")
    file = open(log_path, read=true,  write=true, create=true, append=true)

    module_name = something(_module, "nothing")
    msg_timestamp = stamp_builder(rank)

    formatted_message = "$(msg_timestamp) $(level_name) $message"
    if logger.show_info_source || level != Logging.Info
        formatted_message *= " -@-> $(filepath):$(line)"
    end
    formatted_message *= "\n"

    write(file, formatted_message)
    close(file)
    return nothing
end

end


text:
abstract type AbstractNNPolicy  <: Policy end

## NN Policy interface

"""
    getnetwork(policy)
    return the  value network of the policy
"""
function getnetwork end

"""
    resetstate!(policy)
reset the hidden states of a policy
"""
function resetstate! end

struct NNPolicy{P <: Union{MDP, POMDP}, Q, A} <: AbstractNNPolicy
    problem::P
    qnetwork::Q
    action_map::Vector{A}
    n_input_dims::Int64
    sample_prob::Bool
    action_probability::Function
end

function getnetwork(policy::NNPolicy)
    return policy.qnetwork
end

function resetstate!(policy::NNPolicy)
    Flux.reset!(policy.qnetwork)
end

actionmap(p::NNPolicy) = p.action_map

function _action(policy::NNPolicy{P,Q,A}, o::AbstractArray{T, N}) where {P<:Union{MDP,POMDP},Q,A,T<:Real,N}
    if ndims(o) == policy.n_input_dims
        obatch = reshape(o, (size(o)...,1))
        vals = policy.qnetwork(obatch)
        if policy.sample_prob
            aprobs =  dropdims(policy.action_probability(obatch), dims = 2)
            probs = aprobs .* dropdims(vals, dims = 2)
            probs = (sum(probs) == 0) ? aprobs : probs ./ sum(probs)
            policy.action_map[rand(Categorical(probs))]
        else
            return policy.action_map[argmax(vals)]
        end
    else
        throw("NNPolicyError: was expecting an array with $(policy.n_input_dims) dimensions, got $(ndims(o))")
    end
end

function Distributions.logpdf(policy::NNPolicy{P,Q,A}, o, a) where {P<:Union{MDP,POMDP},Q,A}
    @assert policy.sample_prob
    obatch = reshape(o, (size(o)...,1))
    vals = policy.qnetwork(obatch)
    aprobs =  dropdims(policy.action_probability(obatch), dims = 2)
    probs = aprobs .* dropdims(vals, dims = 2)
    probs = (sum(probs) == 0) ? aprobs : probs ./ sum(probs)
    log(probs[findfirst(policy.action_map .== [a])])
end

function _actionvalues(policy::NNPolicy{P,Q,A}, o::AbstractArray{T,N}) where {P<:Union{MDP,POMDP},Q,A,T<:Real,N}
    if ndims(o) == policy.n_input_dims
        obatch = reshape(o, (size(o)...,1))
        return dropdims(policy.qnetwork(obatch), dims=2)
    else
        throw("NNPolicyError: was expecting an array with $(policy.n_input_dims) dimensions, got $(ndims(o))")
    end
end

function _value(policy::NNPolicy{P}, o::AbstractArray{T,N}) where {P<:Union{MDP,POMDP},T<:Real,N}
    if ndims(o) == policy.n_input_dims
        obatch = reshape(o, (size(o)...,1))
        vals = policy.qnetwork(obatch)
        if policy.sample_prob
            return sum(policy.action_probability(obatch) .* vals)
        else
            return maximum(vals)
        end
    else
        throw("NNPolicyError: was expecting an array with $(policy.n_input_dims) dimensions, got $(ndims(o))")
    end
end

function POMDPs.action(policy::NNPolicy{P}, s) where {P <: MDP}
    _action(policy, convert_s(Array{Float32}, s, policy.problem))
end

function POMDPs.action(policy::NNPolicy{P}, o) where {P <: POMDP}
    _action(policy, convert_o(Array{Float32}, o, policy.problem))
end

function POMDPPolicies.actionvalues(policy::NNPolicy{P}, s) where {P<:MDP}
    _actionvalues(policy, convert_s(Array{Float32}, s, policy.problem))
end

function POMDPPolicies.actionvalues(policy::NNPolicy{P}, o) where {P<:POMDP}
    _actionvalues(policy, convert_o(Array{Float32}, o, policy.problem))
end

function POMDPs.value(policy::NNPolicy{P}, s) where {P <: MDP}
    _value(policy, convert_s(Array{Float32}, s, policy.problem))
end

function POMDPs.value(policy::NNPolicy{P}, o) where {P <: POMDP}
    _value(policy, convert_o(Array{Float32}, o, policy.problem))
end


text:
using Distributions,StatsPlots, ForwardDiff

n_samples = 1000
n = 10_000

x = randn(n,50)
theta = [repeat([1],10); repeat([0],40)]
y = x*theta + randn(n)

prior = Laplace(0,2)

ll(theta) = sum(logpdf.(Normal.(0,1),y-x*theta))

post(theta) = ll(theta) + sum(logpdf.(prior,theta))

init = rand(prior,50)

grad(theta) = ForwardDiff.grad(post,theta)

### start here
#

old_par = init

pars = zeros(n_samples,50)

uni = Uniform(0,1)

accept = 0

var_trans = 0.20

for i = 1:n_samples
    old_transition = MvNormal(old_par,var_trans)
    new_par = rand(old_transition)
    new_transition = MvNormal(new_par,var_trans)
    r = post(new_par) + logpdf(new_transition,old_par) - (post(old_par)+ logpdf(old_transition,new_par))
    a = rand(uni)
    r = min(1,exp(r))
    if r > a
        pars[i,:] = new_par
        global old_par = new_par
        global accept = accept + 1
    else
        pars[i,:] = old_par
    end
    println("Accept ratio ", accept/i*100, "%")
end

StatsPlots.density(pars[:,21])

mm = mapslices(mode,pars,dims = 1)

scatter(mm')

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
# This file includes DAESystems

import DifferentialEquations: DAEProblem
import Sundials: IDA
import UUIDs: uuid4

"""
    @def_dae_system ex

where `ex` is the expression to define to define a new AbstractDAESystem component type. The usage is as follows:
```julia
@def_dae_system mutable struct MyDAESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractDAESystem
    param1::T1 = param1_default                 # optional field
    param2::T2 = param2_default                 # optional field
    param3::T3 = param3_default                 # optional field
        ⋮
    paramN::TN = paramN_default                 # optional field
    righthandside::RH = righthandside_function  # mandatory field
    readout::RO = readout_function              # mandatory field
    state::ST = state_default                   # mandatory field
    stateder::ST = stateder_default             # mandatory field
    diffvars::Vector{Bool} = diffvars_default   # mandatory field
    input::IP = input_default                   # mandatory field
    output::OP = output_default                 # mandatory field
end
```
Here, `MyDAESystem` has `N` parameters. `MyDAESystem` is represented by the `righthandside` and `readout` function. `state`,
'stateder`, `diffvars`, `input` and `output` is the initial state, initial value of differential variables, vector signifing
differetial variables, input port and output port of `MyDAESystem`.

!!! warning
    `righthandside` must have the signature
    ```julia
    function righthandside(out, dx, x, u, t, args...; kwargs...)
        out .= .... # update out
    end
    ```
    and `readout` must have the signature
    ```julia
    function readout(x, u, t)
        y = ...
        return y
    end
    ```

!!! warning
    New DAE system must be a subtype of `AbstractDAESystem` to function properly.

# Example
```julia
julia> @def_dae_system mutable struct MyDAESystem{RH, RO, ST, IP, OP} <: AbstractDAESystem
        righthandside::RH = function sfuncdae(out, dx, x, u, t)
                out[1] = x[1] + 1 - dx[1]
                out[2] = (x[1] + 1) * x[2] + 2
            end
        readout::RO = (x,u,t) -> x
        state::ST = [1., -1]
        stateder::ST = [2., 0]
        diffvars::Vector{Bool} = [true, false]
        input::IP = nothing
        output::OP = Outport(1)
        end

julia> ds = MyDAESystem();
```
"""
macro def_dae_system(ex)
    checksyntax(ex, :AbstractDAESystem)
    appendcommonex!(ex)
    foreach(nex -> appendex!(ex, nex), [
        :( alg::$ALG_TYPE_SYMBOL = Causal.IDA() ),
        :( integrator::$INTEGRATOR_TYPE_SYMBOL = Causal.construct_integrator(
            Causal.DAEProblem, input, righthandside, state, t, modelargs, solverargs; alg=alg,
            stateder=stateder,
            modelkwargs=(; zip((keys(modelkwargs)..., :differential_vars), (values(modelkwargs)..., diffvars))...),
            solverkwargs=solverkwargs, numtaps=3) )
        ])
    quote
        Base.@kwdef $ex
    end |> esc
end

##### Defien DAE system library

"""
    $TYPEDEF

Constructs a generic DAE system.

# Example
```jldoctest
julia> function sfuncdae(out, dx, x, u, t)
           out[1] = x[1] + 1 - dx[1]
           out[2] = (x[1] + 1) * x[2] + 2
       end;

julia> ofuncdae(x, u, t) = x;

julia> x0 = [1., -1];

julia> dx0 = [2., 0.];

julia> DAESystem(righthandside=sfuncdae, readout=ofuncdae, state=x0, input=nothing, output=Outport(1), diffvars=[true, false], stateder=dx0)
DAESystem(righthandside:sfuncdae, readout:ofuncdae, state:[1.0, -1.0], t:0.0, input:nothing, output:Outport(numpins:1, eltype:Outpin{Float64}))
```
"""
@def_dae_system mutable struct DAESystem{RH,
                                         RO,
                                         ST <: AbstractVector{<:Real},
                                         SD <: AbstractVector{<:Real},
                                         DV <: AbstractVector{<:Bool},
                                         IP <: Union{<:Inport, <:Nothing},
                                         OP <: Union{<:Outport,<:Nothing}} <: AbstractDAESystem
    "Right-hand-side function."
    righthandside::RH
    "Readout function"
    readout::RO
    "Initial state"
    state::ST
    "Initial state derivative"
    stateder::SD
    "The `true` elements of `diffvars` correspond to differential variables"
    diffvars::DV
    "Input. Expected to be an `Inport` of `Nothing`"
    input::IP
    "Output port"
    output::OP
end

"""
    $TYPEDEF

Constructs a Robertson systme with the dynamcis
```math
\\begin{array}{l}
    \\dot{x}_1 = -k_1 x_1 + k_3 x_2 x_3 \\\\[0.25cm]
    \\dot{x}_2 = k_1 x_1 - k_2 x_2^2 - k_3 x_2 x_3 \\\\[0.25cm]
    1 = x_1 + x_2 + x_3
\end{array}
```
"""
@def_dae_system mutable struct RobertsonSystem{T1 <: Real,
                                               T2 <: Real,
                                               T3 <: Real,
                                               RH,
                                               RO,
                                               IP <: Union{<:Inport, <:Nothing},
                                               OP <: Union{<:Outport,<:Nothing},
                                               T4 <: AbstractVector{<:Real},
                                               T5 <: AbstractVector{<:Real},
                                               T6 <: AbstractVector{<:Bool}} <: AbstractDAESystem
    "k1"
    k1::T1 = 0.04
    "k2"
    k2::T2 = 3e7
    "k3"
    k3::T3 = 1e4
    "Right-hand-side function"
    righthandside::RH = function robertsonrhs(out, dx, x, u, t)
        out[1] = -k1 * x[1] + k3 * x[2] * x[3] - dx[1]
        out[2] = k1 * x[1] - k2 * x[2]^2 - k3 * x[2] * x[3] - dx[2]
        out[3] = x[1] + x[2] + x[3] - 1
    end
    "Readout function"
    rightout::RO = (x, u, t) -> x[1:2]
    "Initial state"
    state::T4 = [1., 0., 0.]
    "Initial state derivative"
    stateder::T5 = [-k1, k1, 0.]
    "The `true` elements of `diffvars` correspond to differential variables"
    diffvars::T6 = [true, true, false]
    "Input. Expected to be an `Inport` or `Nothing`"
    input::IP = nothing
    "Output port"
    output::OP = Outport(2)
end

"""
    $TYPEDEF

Constructs a Pendulum systme with the dynamics
```math
\\begin{array}{l}
    \\dot{x}_1 = x_3 \\\\[0.25cm]
    \\dot{x}_2 = x_4 \\\\[0.25cm]
    \\dot{x}_3 = -\\dfrac{F}{m l} x_1 \\\\[0.25cm]
    \\dot{x}_4 = g \\dfrac{F}{l} x_2 \\\\[0.25cm]
    0 = x_1^2 + x_2^2 - l^2
\\end{array}
```
where ``F`` is the external force, ``l`` is the length, ``m`` is the mass and ``g`` is the accelaration of gravity.
"""
@def_dae_system mutable struct PendulumSystem{T1 <: Real,
                                              T2 <: Real,
                                              T3 <: Real,
                                              T4 <: Real,
                                              RH,
                                              RO,
                                              IP <: Union{<:Inport, <:Nothing},
                                              OP <: Union{<:Outport,<:Nothing},
                                              T5 <: AbstractVector{<:Real},
                                              T6 <: AbstractVector{<:Real},
                                              T7 <: AbstractVector{<:Bool}} <: AbstractDAESystem
    "Force"
    F::T1 = 1.
    "Length"
    l::T2 = 1.
    "Gravitational acceleration"
    g::T3 = 9.8
    "Mass"
    m::T4 = 1.
    "Right-hand-side function"
    righthandside::RH = function pendulumrhs(out, dx, x, u, t)
        out[1] = x[3] - dx[1]
        out[2] = x[4] - dx[2]
        out[3] = - F / (m * l) * x[1] - dx[3]
        out[4] = g * F / l  * x[2] - dx[4]
        out[5] = x[1]^2 + x[2]^2 - l^2
    end
    "Readout function"
    readout::RO = (x, u, t) -> x[1:4]
    "Initial state"
    state::T5 = [1., 0., 0., 0., 0.]
    "Initial state derivative"
    stateder::T6 = [0., 0., -1., 0., 0.]
    "`true` elements of `diffvars` correspond to diferential variables"
    diffvars::T7 = [true, true, true, true, false]
    "Input. Expected to be an `Inport` of `Nothing`"
    input::IP = nothing
    "Output port"
    output::OP = Outport(4)
end

"""
    $TYPEDEF

Construsts a RLC system with the dynamics
```math
\\begin{array}{l}
    \\dot{x}_1 = x_3 \\\\[0.25cm]
    \\dot{x}_2 = x_4 \\\\[0.25cm]
    \\dot{x}_3 = -\\dfrac{F}{m l} x_1 \\\\[0.25cm]
    \\dot{x}_4 = g \\dfrac{F}{l} x_2 \\\\[0.25cm]
    0 = x_1^2 + x_2^2 - l^2
\\end{array}
```
where ``F`` is the external force, ``l`` is the length, ``m`` is the mass and ``g`` is the accelaration of gravity.
"""
@def_dae_system mutable struct RLCSystem{T1 <: Real,
                                         T2 <: Real,
                                         T3 <: Real,
                                         RH,
                                         RO,
                                         IP <: Union{<:Inport, <:Nothing},
                                         OP <: Union{<:Outport,<:Nothing},
                                         T4 <: AbstractVector{<:Real},
                                         T5 <: AbstractVector{<:Real},
                                         T6 <: AbstractVector{<:Bool}} <: AbstractDAESystem
    "Resistance"
    R::T1 = 1.
    "Inductance"
    L::T2 = 1.
    "Capacitance"
    C::T3 = 1.
    "Right-hand-side function."
    righthandside::RH = function pendulumrhs(out, dx, x, u, t)
        out[1] = 1 / C * x[4] - dx[1]
        out[2] = 1 / L * x[4] - dx[2]
        out[3] = x[3] + R * x[5]
        out[4] = x[1] + x[2] + x[3] + u[1](t)
        out[5] = x[4] - x[5]
    end
    "Readout function"
    readout::RO = (x, u, t) -> x[1:2]
    "Initial state"
    state::T4 = [0., 0., 0., 0., 0.]
    "Initial state derivative"
    stateder::T5 = [0., 0., 0., 0., 0.]
    "`true` elements fo `diffvars` correspond to differential variables"
    diffvars::T6 = [true, true, false, false, false]
    "Input. Expected to an `Inport` or `Nothing`"
    input::IP = Inport(1)
    "Output port"
    output::OP = Outport(2)
end

##### Pretty printing

show(io::IO, ds::DAESystem) = print(io,
    "DAESystem(righthandside:$(ds.righthandside), readout:$(ds.readout), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::RobertsonSystem) = print(io,
    "RobersonSystem(k1:$(ds.k1), k2:$(ds.k2), k2:$(ds.k3), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::PendulumSystem) = print(io,
    "PendulumSystem(F:$(ds.F), m:$(ds.m), l:$(ds.l), g:$(ds.g), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::RLCSystem) = print(io,
    "RLCSystem(R:$(ds.R), L:$(ds.L), C:$(ds.C), state:$(ds.state), t:$(ds.t), ",
    "input:$(ds.input), output:$(ds.output))")


text:
@testset "NfRel" begin
  include("NfRel/NfRelOrd.jl")
  include("NfRel/Ideal.jl")
  include("NfRel/FracIdeal.jl")
  include("NfRel/NfRel.jl")
  include("NfRel/Elem.jl")
  include("NfRel/NEQ_Kirschmer.jl")
end


text:
using BlogPlots
using Documenter

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

makedocs(;
    modules=[BlogPlots],
    authors="Erik Engheim <erik.engheim@mac.com> and contributors",
    repo="https://github.com/ordovician/BlogPlots.jl/blob/{commit}{path}#{line}",
    sitename="BlogPlots.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
#FIXME update to manifolds
#FIXME add to tests
# test simulated PPE in reverse

using Test
using RoME
using DistributedFactorGraphs
using Manifolds: hat, exp
##

@testset "ensure solveParametricBinary is working" begin

##

fg = initfg()

addVariable!(fg, :x0, Pose2)
pp = PriorPose2( MvNormal(zeros(3), diagm(0.01*ones(3)) ))
addFactor!(fg, [:x0;], pp)

# reference ppe on :x0
refVal = zeros(3)
refKey = :simulated
ppe = DFG.MeanMaxPPE(refKey, refVal, refVal, refVal)
setPPE!(fg[:x0], refKey, DFG.MeanMaxPPE, ppe)

pp2 = Pose2Pose2( MvNormal([0;0;-pi+0.01], diagm(0.03*ones(3)) ))

##

M = getManifold(Pose2())
ϵ = getPointIdentity(Pose2())

X = hat(M, ϵ, [0;0;-pi]) #measurement
p = ϵ # variable from
q = ϵ  # variable to
@test isapprox( abs.(calcFactorResidualTemporary(pp2, (Pose2, Pose2), X, (p,  q))), [0;0;pi] )

q = exp(M, ϵ, hat(M, ϵ, [0;0;-pi])) # variable to
@test isapprox( calcFactorResidualTemporary(pp2, (Pose2, Pose2), X, (p,  q)), [0;0;0], atol=1e-14 )

q = exp(M, ϵ, hat(M, ϵ, [0;0;pi])) # variable to
@test isapprox( calcFactorResidualTemporary(pp2, (Pose2, Pose2), X, (p,  q)), [0;0;0], atol=1e-14 )

# wXjhat = SE2(zeros(3))*SE2([0;0;-pi])
# jXjhat = SE2(wxj) \ wXjhat
# return se2vee(jXjhat)

## do in factor graph

addVariable!(fg, :x1, Pose2)
addFactor!( fg, [:x0; :x1], pp2, graphinit=false )

##

refVal = accumulateFactorMeans(fg, [:x0f1; :x0x1f1])

@test isapprox(refVal[1:2], [0;0], atol=5e-4)
@test 0.9pi < abs(refVal[3])

## alternate API

isAlready, simPPE, genLabel = IIF._checkVariableByReference(fg, :x0, r"x\\d+", Pose2, pp2)

@test isapprox(simPPE.suggested[1:2], [0;0], atol=1e-2)
@test 0.9pi < abs(simPPE.suggested[3])

##

end

@testset "Secondary binary parametric solve test case" begin

##

fg = initfg()

addVariable!(fg, :x2, Pose2)
addFactor!(fg, [:x2],
            PriorPose2( MvNormal( [
                        15.000000000016204;
                        8.660254037814505;
                        2.0943951023931953; #+1e-5;
                      ],
                      diagm([0.01;0.01;0.01]) )) )
#

addVariable!(fg, :x3, Pose2)

pp = Pose2Pose2(MvNormal([10;0;pi/3], diagm([0.01;0.01;0.01])))
addFactor!(fg, [:x2;:x3], pp, graphinit=false)

##

# nlsolve failing ?
meas = [10.0, 0.0, 1.0471975511965976]
X1 = [15.000000000016204, 8.660254037814505, 2.0943951023931953]
# X2 = [0.0004014798555661571, 0.0006953833518543132, 1.206557444798529e-10]
X2 = [10.00004891350537; 17.320479835550103; 4.498439149584132e-6]

M = getManifold(Pose2())
ϵ = getPointIdentity(Pose2())

X = hat(M, ϵ, meas) #measurement
p = exp(M, ϵ, hat(M, ϵ, X1)) # variable from
q = exp(M, ϵ, hat(M, ϵ, X2)) # variable to

## NLsolve not moving the cost value and "terminating too early"
# res = [3.469264352442895e-5, -2.002964655985515e-5, -3.1415887889926792]
# x = [9.999965307450084, 17.32052810543953, -3.864597114787283e-6]
# res = [3.469264352798169e-5, -2.0029646563407863e-5, -3.141588788992679]
# x = [9.999965307450083, 17.32052810543953, -3.864597114889992e-6]
# res = [3.469264352798169e-5, -2.0029646563407863e-5, -3.141588788992679]
# x = [9.999965307450083, 17.32052810543953, -3.864597114941346e-6]

# @test isapprox(

resVal = calcFactorResidualTemporary( pp, (Pose2, Pose2), X, (p, q))

@test isapprox( resVal[1:2], [0;0], atol=1e-4)
@test isapprox( abs(resVal[3]), π, atol=1e-4)

# test positive angle
X2 = [10.00004891350537; 17.320479835550103; π]
q = exp(M, ϵ, hat(M, ϵ, X2)) # variable to
resVal = calcFactorResidualTemporary( pp, (Pose2, Pose2), X, (p, q))
@test isapprox( resVal[1:2], [0;0], atol=1e-4)
@test isapprox( abs(resVal[3]), 0, atol=1e-4)

# test negative angle
X2 = [10.00004891350537; 17.320479835550103; -π]
q = exp(M, ϵ, hat(M, ϵ, X2)) # variable to
resVal = calcFactorResidualTemporary( pp, (Pose2, Pose2), X, (p, q))
@test isapprox( resVal[1:2], [0;0], atol=1e-4)
@test isapprox( abs(resVal[3]), 0, atol=1e-4)

##

refVal = accumulateFactorMeans(fg, [:x2f1; :x2x3f1])

@test isapprox(refVal[1:2], [10;17.32], atol=1e-2)
@test isapprox( abs(refVal[3]), pi, atol = 1e-2)

##

end

@testset "Test canonical FG Honeycomb generation" begin

##

# build the graph
fg = RoME.generateCanonicalFG_Honeycomb!()

# check that pose :x3 has rotation near +-pi
t,m,gn = IIF._checkVariableByReference(fg, :x2, r"x\d+", Pose2, getFactorType(fg, :x2x3f1),destPrefix=:x, srcNumber=2)

@test isapprox( m.suggested[1], 10, atol = 1e-1)
@test isapprox( m.suggested[2], 17.32, atol = 1e-1)

@test isapprox( abs(m.suggested[3]), π, atol = 1e-2)

##

end

#


text:
using LibClang
using Documenter

makedocs(;
    modules=[LibClang],
    authors="Yupei Qi <qiyupei@gmail.com>",
    repo="https://github.com/Gnimuc/LibClang.jl/blob/{commit}{path}#L{line}",
    sitename="LibClang.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Gnimuc.github.io/LibClang.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Gnimuc/LibClang.jl",
)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualNetworkRulesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates the specified virtual network rule. During update, the virtual network rule with the specified name will be replaced with this new virtual network rule.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: virtualNetworkRuleName::String (required)
Param: parameters::CreateOrUpdateVirtualNetworkRuleParameters (required)
Param: api_version::String (required)
Return: VirtualNetworkRule
"""
function _swaggerinternal_virtualNetworkRulesCreateOrUpdate(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", VirtualNetworkRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/virtualNetworkRules/{virtualNetworkRuleName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkRuleName", virtualNetworkRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworkRulesCreateOrUpdate(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesCreateOrUpdate(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworkRulesCreateOrUpdate(_api::VirtualNetworkRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, parameters, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesCreateOrUpdate(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, parameters, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified virtual network rule from the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: virtualNetworkRuleName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_virtualNetworkRulesDelete(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/virtualNetworkRules/{virtualNetworkRuleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkRuleName", virtualNetworkRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworkRulesDelete(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesDelete(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworkRulesDelete(_api::VirtualNetworkRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesDelete(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified Data Lake Store virtual network rule.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: virtualNetworkRuleName::String (required)
Param: api_version::String (required)
Return: VirtualNetworkRule
"""
function _swaggerinternal_virtualNetworkRulesGet(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/virtualNetworkRules/{virtualNetworkRuleName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkRuleName", virtualNetworkRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworkRulesGet(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesGet(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworkRulesGet(_api::VirtualNetworkRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesGet(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists the Data Lake Store virtual network rules within the specified Data Lake Store account.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Return: VirtualNetworkRuleListResult
"""
function _swaggerinternal_virtualNetworkRulesListByAccount(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkRuleListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/virtualNetworkRules", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworkRulesListByAccount(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesListByAccount(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworkRulesListByAccount(_api::VirtualNetworkRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesListByAccount(_api, subscriptionId, resourceGroupName, accountName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates the specified virtual network rule.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: virtualNetworkRuleName::String (required)
Param: api_version::String (required)
Param: parameters::UpdateVirtualNetworkRuleParameters
Return: VirtualNetworkRule
"""
function _swaggerinternal_virtualNetworkRulesUpdate(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; parameters=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", VirtualNetworkRule, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.DataLakeStore/accounts/{accountName}/virtualNetworkRules/{virtualNetworkRuleName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkRuleName", virtualNetworkRuleName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function virtualNetworkRulesUpdate(_api::VirtualNetworkRulesApi, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; parameters=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesUpdate(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, api_version; parameters=parameters, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function virtualNetworkRulesUpdate(_api::VirtualNetworkRulesApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, accountName::String, virtualNetworkRuleName::String, api_version::String; parameters=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_virtualNetworkRulesUpdate(_api, subscriptionId, resourceGroupName, accountName, virtualNetworkRuleName, api_version; parameters=parameters, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export virtualNetworkRulesCreateOrUpdate, virtualNetworkRulesDelete, virtualNetworkRulesGet, virtualNetworkRulesListByAccount, virtualNetworkRulesUpdate


text:
abstract type AbstractContext end

mutable struct DefaultContext <: AbstractContext
    index::Index
    trans_units::Vector{TranslationUnit}
    options::Dict{String,Bool}
    libname::String
    force_name::String
    common_buffer::OrderedDict{Symbol,ExprUnit}
    api_buffer::Vector
    next_cursor::Dict{CLCursor,Union{CLCursor,Nothing}}
    queue::Queue{CLCursor}
    visited::Set{CLCursor}
    anonymous_counter::Int
    exclude_symbols::Set{String}
    only_select_symbols::Set{String}
    fields_align::Dict{Tuple{Symbol,Symbol},Int}
end
DefaultContext(index::Index, exclude_symbols, only_select_symbols, fields_align) = DefaultContext(index, TranslationUnit[], Dict{String,Bool}(),
                                              "libxxx", "", OrderedDict{Symbol,ExprUnit}(), [],
                                              Dict{CLCursor,Union{CLCursor,Nothing}}(),
                                              Queue{CLCursor}(), Set{CLCursor}(), 0,
                                              exclude_symbols, only_select_symbols, fields_align)
DefaultContext(diagnostic::Bool=true) = DefaultContext(Index(diagnostic))

parse_header!(ctx::AbstractContext, header::AbstractString; args...) = push!(ctx.trans_units, parse_header(header; args...))
parse_headers!(ctx::AbstractContext, headers::Vector{String}; args...) = (ctx.trans_units = parse_headers(headers; args...);)


text:
# ArgMacros
# argsparsing.jl

#=
The first two functions are used in "setting up" the environment for parsing
=#

"""Split up multiple flag arguments, e.g. "-zx-v2f" → ["-z", "-x", "-v", "2", "-f"]"""
function _split_multiflag(s::AbstractString)::Vector{String}
    splitflags = Vector{String}()

    max_i::Int = length(s) + 1
    i::Int = 2 # First character should always be '-'
    while i < max_i
        if s[i] == '-' && (i == max_i - 1 || (!isdigit(s[i+1]) && s[i+1] !== '.')) # Ignore
            i += 1
        elseif isdigit(s[i]) || s[i] === '-' || s[i] === '.' # Allow numbers inside group of flags with no space or =
            nextletter::Int = something(findnext((c -> isletter(c) || c == '-'), s, i + 1), max_i)
            push!(splitflags, s[i:nextletter - 1])
            i = nextletter
        elseif isletter(s[i]) # Split each flag off
            push!(splitflags, "-$(s[i])")
            i += 1
        else
            _quit_try_help("Error parsing argument(s) \"$s\", unrecognized character in flag sequence")
        end
    end
    return splitflags
end

"Split the arguments into a usable form with separated tokens"
function _split_arguments(args::Vector{String})::Vector{String}
    splitargs = Vector{String}()

    for arg in args
        # Flags/options
        if arg[1] == '-'
            arg_split = split(arg, '=', limit=2)::Vector{SubString{String}}

            # Handle chained flags or flags with Int values like "-xzfv" or "-O3"
            # But do not mess up negative numbers
            if length(arg_split[1]) > 2 && arg_split[1][2] !== '-' && arg_split[1][2] !== '.' && !isdigit(arg_split[1][2])
                append!(splitargs, _split_multiflag(arg_split[1]))

                if length(arg_split) == 2
                    push!(splitargs, arg_split[2])
                end
            # Handle other flags
            else
                append!(splitargs, arg_split)
            end
        # Values
        else
            push!(splitargs, arg)
        end
    end

    return splitargs
end

#=
The remaining functions are used in parsing individual arguments
=#

"Find the index of the (first) occurrence of one of the specified flags in args"
function _get_option_idx(args::Vector{String}, flags::Vector{String})::Union{Int, Nothing}
    return findfirst(in(flags), args)
end

"Remove the option and value for any flag in flags and return the given value"
function _pop_argval!(args::Vector{String}, flags::Vector{String})::Union{String, Nothing}
    option_idx = _get_option_idx(args, flags)
    if isnothing(option_idx)
        return nothing
    elseif option_idx == length(args) # Found flag but no value to go with it
        _quit_try_help("Argument $(args[option_idx]) requires a value")
    else
        val::String = args[option_idx + 1]
        deleteat!(args, (option_idx, option_idx + 1))
        return val
    end
end

"Remove a flag if present and return whether the flag was found"
function _pop_flag!(args::Vector{String}, flags::Vector{String})::Bool
    flag_idx = _get_option_idx(args, flags)

    if isnothing(flag_idx)
        return false
    else
        deleteat!(args, flag_idx)
        return true
    end
end

"Remove multiple occurrences of flag from args and return the count removed"
function _pop_count!(args::Vector{String}, flag::String)::Int
    initial_arg_length::Int = length(args)
    i::Int = 1
    while i <= length(args)
        if args[i] == flag
            deleteat!(args, i)
        else
            i += 1
        end
    end
    return initial_arg_length - length(args)
end

"Convert string to Number type with parse or other type with direct conversion."
function _converttype!(::Type{T}, s::String, name::String)::T where T
    try
        if T <: Number
            # Allow floating point value to be passed to Int argument
            return T(parse(Float64, s)::Float64)
        else
            return T(s)
        end
    catch
        _quit_try_help("Invalid type for argument $name")
    end
end

"Quit program if the value is nothing"
function _converttype!(::Type{T}, ::Nothing, name::String)::Nothing where T
    _quit_try_help("Argument $name missing")
end

"If x is not a string (i.e. literal passed as default) assume direct conversion is possible. Does NOT validate."
function _converttype!(::Type{T}, x::Any, ::String)::T where T
    return T(x)
end


text:
import Gen
import Distributions
import LinearAlgebra#: norm, cross, det
import Rotations
import StaticArrays

const Rot2 = Rotations.RotMatrix{2,Float64,4}

function Base.rand(::Type{Rot2})
    a = Gen.uniform(-pi, pi)
    return Rotations.RotMatrix{2,Float64}(a)
end

total_area(::Type{Rot2}) = 2 * pi

@inline function _to_array(x::Rot2)
    return Float64[x[1,1], x[2,1]]
end

export Rot2

################################################
# uniform distribution (Haar measure on SO(2)) #
################################################

struct UniformRot2 <: Gen.Distribution{Rot2} end

Gen.logpdf(::UniformRot2, x::Rot2) = -log(total_area(Rot2))

function Gen.logpdf_grad(::UniformRot2, x::Rot2)
    error("Not implemented") # TODO implement
    return (nothing,)
end

Gen.random(::UniformRot2) = rand(Rot2)

Gen.has_output_grad(::UniformRot2) = false

Gen.has_argument_grads(::UniformRot2) = ()

(::UniformRot2)() = Gen.random(UniformRot2())

const uniform_rot2 = UniformRot2()

export uniform_rot2

################################################
# von Mises Fisher based distribution on SO(2) #
################################################

struct VMFRot2 <: Gen.Distribution{Rot2} end

const von_mises_rot2 = VMFRot2()

function Gen.logpdf(::VMFRot2, x::Rot2, mu::Rot2, k::Real)
    dist = Distributions.VonMisesFisher(_to_array(mu), k)
    x_unit_vec = _to_array(x)
    return Distributions.logpdf(dist, x_unit_vec)
end

function Gen.logpdf_grad(::VMFRot2, x::Rot2, mu::Rot2, k::Real)
    error("Not implemented") # TODO implement
    return (nothing,nothing,nothing)
end

function Gen.random(::VMFRot2, mu::Rot2, k::Real)
    dist = Distributions.VonMisesFisher(_to_array(mu), k)
    v = rand(dist)
    mat = StaticArrays.SMatrix{2,2,Float64,4}(v[1], v[2], -v[2], v[1])
    return Rotations.RotMatrix{2,Float64,4}(mat)
end

Gen.has_output_grad(::VMFRot2) = false

Gen.has_argument_grads(::VMFRot2) = (false, false)

(::VMFRot2)(mu, k) = Gen.random(VMFRot2(), mu, k)

export von_mises_rot2


text:
module ReactiveMPNodeTest

using Test
using ReactiveMP

@testset "FactorNode" begin

    @testset "Common" begin

        @test ReactiveMP.as_node_functional_form(() -> nothing) === ReactiveMP.ValidNodeFunctionalForm()
        @test ReactiveMP.as_node_functional_form(2)             === ReactiveMP.UndefinedNodeFunctionalForm()

        @test isdeterministic(Deterministic()) === true
        @test isdeterministic(Deterministic)   === true
        @test isdeterministic(Stochastic())    === false
        @test isdeterministic(Stochastic)      === false
        @test isstochastic(Deterministic())    === false
        @test isstochastic(Deterministic)      === false
        @test isstochastic(Stochastic())       === true
        @test isstochastic(Stochastic)         === true

        @test sdtype(() -> nothing) === Deterministic()
        @test_throws MethodError sdtype(0)
    end

    @testset "@node macro" begin

        # Testing Stochastic node specification

        struct CustomStochasticNode end

        @node CustomStochasticNode Stochastic [ out, (x, aliases = [ xx ]), (y, aliases = [ yy ]), z ]

        @test ReactiveMP.interface_get_index(Val{ :CustomStochasticNode }, Val{ :out }) === 1
        @test ReactiveMP.interface_get_index(Val{ :CustomStochasticNode }, Val{ :x })   === 2
        @test ReactiveMP.interface_get_index(Val{ :CustomStochasticNode }, Val{ :y })   === 3
        @test ReactiveMP.interface_get_index(Val{ :CustomStochasticNode }, Val{ :z })   === 4

        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ :out }) === :out
        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ :x })   === :x
        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ :y })   === :y
        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ :z })   === :z

        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ :xx }) === :x
        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ :yy }) === :y

        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ 1 }) === :out
        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ 2 }) === :x
        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ 3 }) === :y
        @test ReactiveMP.interface_get_name(Val{ :CustomStochasticNode }, Val{ 4 }) === :z

        @test ReactiveMP.collect_factorisation(CustomStochasticNode, ((1, ), (2, ), (3, ), (4, ))) === ((1, ), (2, ), (3, ), (4, ))
        @test ReactiveMP.collect_factorisation(CustomStochasticNode, ((1, 2, ), (3, ), (4, )))     === ((1, 2, ), (3, ), (4, ))
        @test ReactiveMP.collect_factorisation(CustomStochasticNode, ((1, 2, 3, ), (4, )))         === ((1, 2, 3, ), (4, ))
        @test ReactiveMP.collect_factorisation(CustomStochasticNode, ((1, 2, 3, ), (4, )))         === ((1, 2, 3, ), (4, ))
        @test ReactiveMP.collect_factorisation(CustomStochasticNode, ((1, 2, 3, 4, ),))            === ((1, 2, 3, 4, ),)

        @test ReactiveMP.collect_factorisation(CustomStochasticNode, MeanField())         === ((1, ), (2, ), (3, ), (4, ))
        @test ReactiveMP.collect_factorisation(CustomStochasticNode, FullFactorisation()) === ((1, 2, 3, 4, ),)

        @test sdtype(CustomStochasticNode) === Stochastic()

        cx = constvar(:cx, 1.0)
        cy = constvar(:cy, 1.0)
        cz = constvar(:cz, 1.0)

        snode, svar = make_node(CustomStochasticNode, AutoVar(:cout), cx, cy, cz, factorisation = MeanField())

        @test snode        !== nothing
        @test typeof(svar) <: RandomVariable
        @test factorisation(snode) === ((1,), (2,), (3,), (4,))

        # Testing Deterministic node specification

        struct CustomDeterministicNode end

        CustomDeterministicNode(x, y, z) = x + y + z

        @node CustomDeterministicNode Deterministic [ out, (x, aliases = [ xx ]), (y, aliases = [ yy ]), z ]

        @test ReactiveMP.interface_get_index(Val{ :CustomDeterministicNode }, Val{ :out }) === 1
        @test ReactiveMP.interface_get_index(Val{ :CustomDeterministicNode }, Val{ :x })   === 2
        @test ReactiveMP.interface_get_index(Val{ :CustomDeterministicNode }, Val{ :y })   === 3
        @test ReactiveMP.interface_get_index(Val{ :CustomDeterministicNode }, Val{ :z })   === 4

        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ :out }) === :out
        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ :x })   === :x
        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ :y })   === :y
        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ :z })   === :z

        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ :xx }) === :x
        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ :yy }) === :y

        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ 1 }) === :out
        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ 2 }) === :x
        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ 3 }) === :y
        @test ReactiveMP.interface_get_name(Val{ :CustomDeterministicNode }, Val{ 4 }) === :z

        @test ReactiveMP.collect_factorisation(CustomDeterministicNode, ((1, ), (2, ), (3, ), (4, ))) === ((1, 2, 3, 4, ),)
        @test ReactiveMP.collect_factorisation(CustomDeterministicNode, ((1, 2, ), (3, ), (4, )))     === ((1, 2, 3, 4, ),)
        @test ReactiveMP.collect_factorisation(CustomDeterministicNode, ((1, 2, 3, ), (4, )))         === ((1, 2, 3, 4, ),)
        @test ReactiveMP.collect_factorisation(CustomDeterministicNode, ((1, 2, 3, ), (4, )))         === ((1, 2, 3, 4, ),)
        @test ReactiveMP.collect_factorisation(CustomDeterministicNode, ((1, 2, 3, 4, ),))            === ((1, 2, 3, 4, ),)

        @test ReactiveMP.collect_factorisation(CustomDeterministicNode, MeanField())         === ((1, 2, 3, 4, ),)
        @test ReactiveMP.collect_factorisation(CustomDeterministicNode, FullFactorisation()) === ((1, 2, 3, 4, ),)

        @test sdtype(CustomDeterministicNode) === Deterministic()

        cx = constvar(:cx, 1.0)
        cy = constvar(:cy, 1.0)
        cz = constvar(:cz, 1.0)

        snode, svar = make_node(CustomDeterministicNode, AutoVar(:cout), cx, cy, cz, factorisation = MeanField())

        @test snode        === nothing
        @test typeof(svar) <: ConstVariable

        # Testing expected exceptions

        struct DummyStruct end

        @test_throws Exception eval(:(@node DummyStruct NotStochasticAndNotDeterministic [ out, in, x ]))
        @test_throws Exception eval(:(@node DummyStruct Stochastic [ 1, in, x ]))
        @test_throws Exception eval(:(@node DummyStruct Stochastic [ (1, aliases = [ out ]), in, x ]))
        @test_throws Exception eval(:(@node DummyStruct Stochastic [ (out, aliases = [ out ]), in, x ]))
        @test_throws Exception eval(:(@node DummyStruct Stochastic [ (out, aliases = [ 1 ]), in, x ]))
        @test_throws Exception eval(:(@node DummyStruct Stochastic [  ]))

    end

end

end


text:
__precompile__()

module Bedgraph

include("record.jl")
include("helpers.jl")
include("header.jl")
include("reader.jl")
include("writer.jl")

end # module


text:
using Mocking
Mocking.activate()

using CloudWatchLogs
using CloudWatchLogs: MAX_EVENT_SIZE

import AWSCore
import AWSCore.Services: logs, sts
using AWSCore: AWSConfig, aws_config, AWSCredentials, AWSException
using Dates
using EzXML
using HTTP
using Printf
using Memento
using Memento.TestUtils
using Test
using TimeZones
using UUIDs

const LOGGER = getlogger(CloudWatchLogs)

function assume_role(config::AWSConfig, role_arn::AbstractString; kwargs...)
    response = sts(
        config,
        "AssumeRole";
        RoleArn=role_arn,
        RoleSessionName=session_name(),
        kwargs...
    )

    response_creds = response["Credentials"]
    response_user = response["AssumedRoleUser"]

    AWSCredentials(
        response_creds["AccessKeyId"],
        response_creds["SecretAccessKey"],
        response_creds["SessionToken"],
        response_user["Arn"],
    )
end

function session_name()
    user = get(ENV, "USER", "unknown")
    location = gethostname()

    name = "$user@$location"
    ts = string(round(Int64, time()))

    # RoleSessionName must be no more than 64 characters
    max_name_length = 64 - length(ts) - 1

    if length(name) > max_name_length
        name = name[1:max_name_length-3] * "..."
    end

    return "$name-$ts"
end

function stack_output(config::AWSConfig, stack_name::AbstractString)
    outputs = Dict{String, String}()

    # https://github.com/JuliaCloud/AWSCore.jl/issues/41
    response = AWSCore.service_query(
        config;
        service="cloudformation",
        version="2010-05-15",
        operation="DescribeStacks",
        args=[:StackName=>stack_name],
        return_raw=true,
    )

    xml = EzXML.root(EzXML.parsexml(response))
    ns = EzXML.namespace(xml)
    outputs_xml = findall("//ns:Stacks/ns:member[1]/ns:Outputs/ns:member", xml, ["ns"=>ns])
    for output_xml in outputs_xml
        key = string(findfirst("//ns:OutputKey/text()", output_xml, ["ns"=>ns]))
        val = string(findfirst("//ns:OutputValue/text()", output_xml, ["ns"=>ns]))
        outputs[key] = val
    end

    return outputs
end

@testset "CloudWatchLogs.jl" begin
    include("event.jl")
    include("mocked_aws.jl")
    include("online.jl")
end


text:
# function draw!(scene::Scene, h::OpticSim.ParaxialAnalysis.HeadEye.Head; draw_head=true, kwargs...) where {T<:Real}

#     # axes_mat = OpticSimVis.Material(color=RGBA(0.7, 0.2, 0.1, 1.0))
#     # axes = OpticSimVis.Axes(tr=Transform(Vec3(3.0, 3.0, 3.0)), material=axes_mat, name="My Axes")
#     # OpticSimVis.parent!(axes , OpticSimVis.root(scene))

#     name = "HeadEye" #-$(UUIDs.uuid1())"
#     mat = Material(;kwargs...)
#     # head_root = EmptySceneObject(material=mat, name=name)
#     head_root = Axes(
#         tr=transform(OpticSim.ParaxialAnalysis.HeadEye.tr(h)),
#         material=mat,
#         shaft_scale=0.4,
#         axes_scale=10.0,
#         name=name)
#     OpticSimVis.parent!(head_root , OpticSimVis.root(scene))

#     if (draw_head)
#         scale_factor = 60.0;
#         # face = ParaxialAnalysis.HeadEye.transform_mesh(ParaxialAnalysis.HeadEye.get_face_model(), ParaxialAnalysis.HeadEye.tr(h)*Geometry.scale(scale_factor, scale_factor, scale_factor))
#         face_mat = Material(color=RGBA(1.0,0.2,0.2, 0.3))
#         face = ParaxialAnalysis.HeadEye.transform_mesh(ParaxialAnalysis.HeadEye.get_face_model(), Geometry.scale(scale_factor, scale_factor, scale_factor))
#         face_mesh = Mesh(mesh=face, material=face_mat, name="Face")
#         OpticSimVis.parent!(face_mesh , head_root)
#     end

#     # ParaxialAnalysis.HeadEye.draw_local_frame(scene, ParaxialAnalysis.HeadEye.tr(h), "Head Transform")

#     # eyes
#     for eye in ParaxialAnalysis.HeadEye.eyes(h)
#         # Vis.draw!(scene, e, parent_transform = ParaxialAnalysis.HeadEye.tr(h))
#         local name = ParaxialAnalysis.HeadEye.text(eye)
#         eye_so = Axes(
#             tr=transform(OpticSim.ParaxialAnalysis.HeadEye.tr(eye)),
#             material=mat,
#             shaft_scale=0.2,
#             axes_scale=5.0,
#             name=name)
#         OpticSimVis.parent!(eye_so , head_root)

#         pupil = ParaxialAnalysis.HeadEye.pupil(eye)
#         name = "Pupil"
#         pupil_so = Axes(
#             tr=transform(OpticSim.ParaxialAnalysis.HeadEye.tr(pupil)),
#             material=mat,
#             shaft_scale=0.1,
#             axes_scale=2.0,
#             name=name)
#         OpticSimVis.parent!(pupil_so , eye_so)

#     end

#     return Dict{Symbol, Any}(
#         :scene => scene,
#         :scene_object => head_root,
#     )

# end


text:
"""
Implements the needed discharge for each purpose. (Purposes without discharge do not need to implement this function.)
The returned discharge may be an expression based on other variables.

For example, imposing a constraint on the reservoir level within the reservoir is not a discharge: the corresponding
function should return `0.`.
"""
hlp_JuMP_purposes_discharge(model::JuMP.Model, r::Reservoir, period::Period) =
  hlp_JuMP_purposes_discharge(model, getPurposes(r), period)

hlp_JuMP_purposes_discharge(model::JuMP.Model, purposes::Array{Purpose, 1}, period::Period) =
  sum([hlp_JuMP_purposes_discharge(model, p, period) for p in purposes])

# Generic purposes have no needed discharge.
hlp_JuMP_purposes_discharge(model::JuMP.Model, purpose::Purpose, period::Period) = 0.

hlp_JuMP_purposes_discharge(model::JuMP.Model, purpose::WaterWithdrawalPurpose, period::Period) =
  error("[hlp_JuMP_purposes_discharge] Helper not defined for " * string(typeof(purpose)) * ".")

hlp_JuMP_purposes_discharge(model::JuMP.Model, purpose::DeterministicWaterWithdrawalPurpose, period::Period) =
  getNeed(purpose, period)

# Hydropower has no imposed discharge.
hlp_JuMP_purposes_discharge(model::JuMP.Model, purpose::HydropowerDamOutput, period::Period) = 0.


text:
export quantization!, quantization

 # apply quantization and clipping
function quant(x, levels, minv, maxv)
    return clip_v(round((x - minv)/(maxv-minv) * levels) / levels * (maxv-minv) + minv, minv, maxv)
end

f_quant(levels, minv, maxv) = x -> quant(x, levels, minv, maxv)
comb_quant(x, n) = n

quantization(X::AbstractArray, levels; minv=0, maxv=1) =
    quantization!(copy(X), levels, minv=minv, maxv=maxv)

quantization!(X::AbstractArray, levels; minv=0, maxv=1) =
    apply_noise!(comb_quant, f_quant(levels-1, minv, maxv), X, false)

"""
    quantization(X, levels; minv=0, maxv=1)

Returns array `X` discretized to `levels` different values.
Therefore the array is discretized.
`levels` describes how many different value steps the resulting image has.
`minv=0` and `maxv` indicate the minimum and maximum possible values of the images.
In RGB and Gray images this is usually 0 and 1.
There is also `quantization!` available.
"""
quantization

max_issues_repo_name: roflmaostc/Noise

max_forks_repo_name: roflmaostc/Noise


text:
using Documenter
using FIRLSFilterDesign

makedocs(
    clean = false,
    doctest = false,
    sitename = "FIRLSFilterDesign.jl documentation",
    expandfirst = [],
    modules = [FIRLSFilterDesign],
    pages = [
        "Home" => "index.md",
        "Manual" => "manual.md"
    ],
    )


text:
module MultiCellMaps

using Gridap

export MultiCellMap
import Gridap.CellIntegration: integrate
import Base: +, -, *

struct MultiCellMap{N}
  blocks::Vector{<:CellMap}
  fieldids::Vector{NTuple{N,Int}}
end

function integrate(
  mcm::MultiCellMap,
  trian::Triangulation{Z},
  quad::CellQuadrature{Z}) where Z

  blocks = [ integrate(b,trian,quad) for b in mcm.blocks ]

  MultiCellArray(blocks,mcm.fieldids)

end

function (+)(a::MultiCellMap{N},b::MultiCellMap{N}) where N
  blocks = CellMap[]
  append!(blocks,a.blocks)
  append!(blocks,b.blocks)
  fieldids = NTuple{N,Int}[]
  append!(fieldids,a.fieldids)
  append!(fieldids,b.fieldids)
  MultiCellMap(blocks,fieldids)
end

function (-)(a::MultiCellMap{N},b::MultiCellMap{N}) where N
  blocks = CellMap[]
  append!(blocks,a.blocks)
  append!(blocks,[ -k for k in b.blocks])
  fieldids = NTuple{N,Int}[]
  append!(fieldids,a.fieldids)
  append!(fieldids,b.fieldids)
  MultiCellMap(blocks,fieldids)
end

function (+)(b::MultiCellMap{N}) where N
  b
end

function (-)(b::MultiCellMap{N}) where N
  blocks = CellMap[]
  append!(blocks,[ -k for k in b.blocks])
  MultiCellMap(blocks,b.fieldids)
end

function (*)(a::Real,b::MultiCellMap{N}) where N
  blocks = CellMap[]
  append!(blocks,[ a*k for k in b.blocks])
  MultiCellMap(blocks,b.fieldids)
end

end # module


text:
mqtt_read(s::IO, ::Type{UInt16}) = ntoh(read(s, UInt16))

function mqtt_read(s::IO, ::Type{String})
    len = mqtt_read(s, UInt16)
    return String(read(s, len))
end

function mqtt_write(stream::IO, x::Any)
    write(stream, x)
end

function mqtt_write(stream::IO, x::UInt16)
    write(stream, hton(x))
end

function mqtt_write(stream::IO, x::String)
    mqtt_write(stream, convert(UInt16, length(x)))
    write(stream, x)
end

function mqtt_write(stream::IO, x::Array{UInt8})
    mqtt_write(stream, convert(UInt16, length(x)))
    write(stream, x)
end

function write_len(s::IO, len::Int64)
    while true
        b = convert(UInt8, mod(len, 128))
        len = div(len, 128)
        if len > 0
            b = b | 0x80
        end
        write(s, b)
        if(len == 0)
            break
        end
    end
end

function read_len(s::IO)
    multiplier = 1
    value = 0
    while true
        b = read(s, UInt8)
        value += (b & 127) * multiplier
        multiplier *= 128
        if multiplier > 128 * 128 * 128
            throw(ErrorException("malformed remaining length"))
        end
        if (b & 128) == 0
            break
        end
    end
    return value
end

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
module rails
export RAILS

push!(LOAD_PATH, joinpath(pwd(), "models"))
push!(LOAD_PATH, joinpath(pwd(), "data"))

using Knet: Adam
using models: set_optim!, zero_grad_model, train!
using ego_model: EgoModel
using ego_dataset: EgoDataset
using metrics: l1_loss

using YAML

struct RAILS
    config::Dict
    device::AbstractString
    ego_model::EgoModel

    function RAILS(args)
        config = nothing
        open(args.config_path, 'r') do cfg_file
            config = YAML.load(cfg_file)
        end

        ego_model = EgoModel(1 / (args.fps * (args.num_repeat + 1)))
        # set Adam optimizer to ego_model
        set_optim!(ego_model, Adam, args.lr)

        new(config, args.device, ego_model)
    end
end

function train_ego(r::RAILS, locs, rots, spds, acts)
    locs = selectdim(locs, length(size(locs)), 1:2)
    yaws = selectdim(yaws, length(size(yaws)), w:size(yaws)[end]) .* pi ./ 180

    pred_locs = []
    pred_yaws = []

    pred_loc = locs[:, 1]
    pred_yaw = yaws[:, 1]
    pred_spd = spds[:, 1]

    for t = 1:size(locs, 2)-1
        act = acts[:, t]
        pred_loc, pred_yaw, pred_spd = r.ego_model(pred_loc, pred_yaw, pred_spd, act)
        append!(pred_locs, pred_loc)
        append!(pred_yaws, pred_yaw)
    end
    loc_loss = l1_loss(pred_locs, locs[:, 2:end])
    ori_loss = l1_loss(cos.(pred_yaws), cos.(yaws[:, 2:end])) + l1_loss(sin.(pred_yaws), sin.(yaws[:, 2:end]))

    loss = @diff loc_loss + ori_loss

    zero_grad_model(r.ego_model)
    train!(r.ego_model, loss)

    return Dict("loc_loss" => Float32(loc_loss), "ori_loss" => Float32(ori_loss),
        "pred_locs" => Array(pred_locs), "pred_yaws" => Array(pred_yaws),
        "locs" => Array(locs), "yaws" => Array(yaws))
end
end


text:
module TernaryPlots

export ternary_axes, ternary_axes!
export tern2cart, cart2tern

using RecipesBase
using LinearAlgebra: normalize

include("user_helpers.jl")
include("recipe_helpers.jl")
include("ternary_axes.jl")

end # module


text:
# This file was taken from https://github.com/bat/batjl-environment/blob/master/build_sysimage.jl
# Original Author: @oschulz

# This software is licensed under the MIT License (MIT).

import Pkg

replace_default_sysimg = "--replace-default" in ARGS

srcdir = @__DIR__

if replace_default_sysimg
    juliadir = dirname(Sys.BINDIR)

    # Store packages in system location, typically "/path/to/julia/local/share/julia":
    if any(s -> startswith(s, juliadir), DEPOT_PATH)
        filter!(s -> startswith(s, juliadir), DEPOT_PATH)
    end

    # Restrict load path to stdlib and current project:
    if any(s -> s == "@stdlib", LOAD_PATH)
        filter!(s -> s == "@stdlib", LOAD_PATH)
        pushfirst!(LOAD_PATH, "@")
    end

    sysprjdir = joinpath(first(DEPOT_PATH), "environments", "v$(VERSION.major).$(VERSION.minor)")
    mkpath(sysprjdir)
    cp(joinpath(srcdir, "Project.toml"), joinpath(sysprjdir, "Project.toml"), force = true)
    if ispath(joinpath(srcdir, "Manifest.toml"))
        cp(joinpath(srcdir, "Manifest.toml"), joinpath(sysprjdir, "Manifest.toml"), force = true)
    end

    Pkg.activate(sysprjdir)

    # IJulia should always be part of custom default system image
    Pkg.add("IJulia", preserve = Pkg.PRESERVE_ALL)

    # Need to have PackageCompiler installed
    Pkg.add("PackageCompiler", preserve = Pkg.PRESERVE_ALL)
else
    Pkg.activate(srcdir)
end

prj_file = Pkg.project().path
prj_dir = dirname(prj_file)

@info("Building custom system image in project \"$prj_dir\".")

@info("Instantiating project.")
Pkg.instantiate()

@info("Precompiling packages.")
Pkg.precompile()

prj_sysimage = get(Pkg.TOML.parsefile(prj_file), "sysimage", Dict{String,Any}())
sysimg_pkgs = sort(collect(keys(filter(x -> Bool(x[2]), prj_sysimage))))
@info("Package to include in system image: $(join(sysimg_pkgs, " "))")

@info("Building system image.")

import PackageCompiler, Libdl

sysimage_path = replace_default_sysimg ? nothing : joinpath(srcdir, "JuliaSysimage." * Libdl.dlext)

PackageCompiler.create_sysimage(
    Symbol.(sysimg_pkgs),
    sysimage_path = sysimage_path,
    precompile_execution_file = joinpath(srcdir, "precompile_exec.jl"),
    cpu_target = PackageCompiler.default_app_cpu_target(),
    replace_default = replace_default_sysimg
)

if replace_default_sysimg
    @info("Replaced default Julia system image.")
else
    @info("Created custom Julia system image \"$sysimage_path\".")
end

max_stars_repo_name: greimel/julia-binder

max_issues_repo_name: greimel/julia-binder

max_forks_repo_name: greimel/julia-binder


text:
using GmshTools
using Test

try
    gmsh.initialize()
    @test true
catch
    @test false
finally
    gmsh.finalize()
end

include("test_mesh.jl")


text:
using StaticArrays

struct BoundingBox{D,T}
    halflength::T
    center::SVector{D,T}
end

function BoundingBox(points::Vector{SVector{D,T}}) where {D,T}
    min_dim = Vector(points[1])
    max_dim = Vector(points[1])

    for i = 1:length(points)
        for j = 1:D
            min_dim[j] =  min_dim[j] < points[i][j] ? min_dim[j] : points[i][j]
            max_dim[j] =  max_dim[j] > points[i][j] ? max_dim[j] : points[i][j]
        end
    end

    center = @MVector zeros(D)

    length_dim = zeros(T,D)
    for j = 1:D
        length_dim[j] = max_dim[j] - min_dim[j]
        center[j] = (max_dim[j] + min_dim[j])/2.0
    end

    halflength = maximum(length_dim)/2.0

    return BoundingBox(halflength, SVector(center))
end

function getboxframe(bbox::BoundingBox{D,T}) where {D,T}
    if D == 2
        mat = zeros(5, 2)
        sign_x = [1 -1 -1 1 1]
        sign_y = [1 1 -1 -1 1]
        for i = 1:5
            mat[i,1] = bbox.center[1] + sign_x[i]*bbox.halflength
            mat[i,2] = bbox.center[2] + sign_y[i]*bbox.halflength
        end
    elseif D == 3
        sign_x = [1 -1 -1  1  1  1  1  1 NaN  1  1 -1 -1 NaN -1 -1 -1 NaN -1  1]
        sign_y = [1  1 -1 -1  1  1 -1 -1 NaN -1 -1 -1 -1 NaN -1  1  1 NaN  1  1]
        sign_z = [1  1  1  1  1 -1 -1  1 NaN -1 -1 -1  1 NaN -1 -1  1 NaN -1 -1]
        mat = zeros(length(sign_x), 3)
        for i = 1:length(sign_x)
            mat[i,1] = bbox.center[1] + sign_x[i]*bbox.halflength
            mat[i,2] = bbox.center[2] + sign_y[i]*bbox.halflength
            mat[i,3] = bbox.center[3] + sign_z[i]*bbox.halflength
        end
    else
        error("Dimension must be 2 or 3")
    end
    return mat
end

"""
    getchildbox(bbox, id)

Returns a bounding box that this is the `id`th quadrants (2D) /
octand (3D) of the bounding box `bbox`
"""
function getchildbox(bbox::BoundingBox{D,T}, id::Integer) where {D,T}

    hl = bbox.halflength / 2.0

    if D == 2
        sign_x = [1 -1 -1  1]
        sign_y = [1  1 -1 -1]

        center = bbox.center + hl .* SVector(sign_x[id], sign_y[id])
    elseif D == 3
        sign_x = [1 -1 -1  1  1 -1 -1  1]
        sign_y = [1  1 -1 -1  1  1 -1 -1]
        sign_z = [1  1  1  1 -1 -1 -1 -1]
        center = bbox.center + hl .* SVector(sign_x[id], sign_y[id],  sign_z[id])
    else
        error("We support only 2 or 3 dimensions")
    end
    return BoundingBox(hl, center)
end

function whichchildbox(bbox::BoundingBox{D,T}, point::SVector{D,T}) where {D,T}

    if point[1] >= bbox.center[1]
        if  point[2] >= bbox.center[2]
            id = 1
        else
            id = 4
        end
    else
        if point[2] >= bbox.center[2]
            id = 2
        else
            id = 3
        end
    end

    if D == 3 && point[3] < bbox.center[3]
        id += 4
    end

    return id
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
# ---
# title: 1018. Binary Prefix Divisible By 5
# id: problem1018
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/binary-prefix-divisible-by-5/description/>
# hidden: true
# ---
#
# Given an array `A` of `0`s and `1`s, consider `N_i`: the i-th subarray from
# `A[0]` to `A[i]` interpreted as a binary number (from most-significant-bit to
# least-significant-bit.)
#
# Return a list of booleans `answer`, where `answer[i]` is `true` if and only if
# `N_i` is divisible by 5.
#
# **Example 1:**
#
#
#
#     Input: [0,1,1]
#     Output: [true,false,false]
#     Explanation:
#     The input numbers in binary are 0, 01, 011; which are 0, 1, and 3 in base-10.  Only the first number is divisible by 5, so answer[0] is true.
#
#
# **Example 2:**
#
#
#
#     Input: [1,1,1]
#     Output: [false,false,false]
#
#
# **Example 3:**
#
#
#
#     Input: [0,1,1,1,1,1]
#     Output: [true,false,false,false,true,false]
#
#
# **Example 4:**
#
#
#
#     Input: [1,1,1,0,1]
#     Output: [false,false,false,false,false]
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 30000`
#   2. `A[i]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Inpaintings
# This module translates the `inpaint_nans` MATLAB code from John d'Erico

using SparseArrays, SuiteSparse, LinearAlgebra
using Match

include("inpaint.jl")

end # module