text:
@testset "infer generated function" begin
    function code_llvm_lines(g, t)
        buf = IOBuffer()
        code_llvm(buf, g, t)
        return filter(x -> !startswith(x, ";"), split(String(buf), "\n"))
    end

    function llvm_eq(a, b)
        if length(a) != length(b)
            return false
        end

        for (i, j) in zip(a, b)
            r = r"define(.*)@julia"
            mi = match(r, i)
            mj = match(r, j)
            if mi != nothing && mj != nothing
                if mi[1] != mj[1]
                    return false
                end
            elseif mi != nothing || mj != nothing
                return false
            elseif i != j
                return false
            end
        end
        return true
    end

    @testset "simple" begin
        @eval @generated function f(t)
            if t <: Int
                :(t + 3)
            else
                println(t)
            end
        end

        @eval function g(x)
            f(x)
        end

        orig = code_llvm_lines(g, Tuple{Int64})

        dub_tt = Tuple{Nothing, typeof(g), Int64}
        over = code_llvm_lines(TinyCassette.execute, dub_tt)

        @test llvm_eq(orig, over)
    end

    @testset "with Val" begin
        @eval @generated function f(::Val{name}) where {name}
            i = "llvm.nvvm.read.ptx.sreg.tid.$name"
            :(Int(ccall($i, llvmcall, UInt32, ())))
        end

        @eval @inline x() = f(Val(:x))

        @eval function g()
            x()
        end

        orig = code_llvm_lines(g, Tuple{})

        dub_tt = Tuple{Nothing, typeof(g)}
        over = code_llvm_lines(TinyCassette.execute, dub_tt)

        @test llvm_eq(orig, over)
    end
end


text:
# Code for augmented Lagrangian evaluator. Inspired by the excellent NLPModels.jl:
# Ref: https://github.com/JuliaSmoothOptimizers/Percival.jl/blob/master/src/AugLagModel.jl
# Two-sided Lagrangian

const CONSTRAINTS_TYPE = Union{Val{:inequality}, Val{:equality}, Val{:mixed}}

@doc raw"""
    AugLagEvaluator{Evaluator<:AbstractNLPEvaluator, T, VT} <: AbstractPenaltyEvaluator

Augmented-Lagrangian evaluator.

### Description

Takes as input any [`AbstractNLPEvaluator`](@ref) encoding a non-linear problem
```math
\begin{aligned}
       \min_u \quad & f(u)\\
\mathrm{s.t.} \quad & h^♭ ≤ h(u) ≤ h^♯,\\
                    & u^♭ ≤  u   ≤ u^♯,
\end{aligned}
```
and return a new evaluator reformulating the original problem
by moving the $m$ constraints $h^♭ ≤ h(u) ≤ h^♯$ into the objective
using a set of penalties $ϕ_1, ⋯, ϕ_m$ and multiplier estimates
$λ_1, ⋯, λ_m$:
```math
\begin{aligned}
    \min_u \quad & f(u) + \sum_{i=1}^m ϕ_i(h_i, λ_i)   \\
\mathrm{s.t.} \quad &  u^♭ ≤  u   ≤  u^♯,
\end{aligned}
```

This evaluator considers explicitly the inequality constraints,
without reformulating them by introducing slack variables. Each
penalty $ϕ_i$ is defined as
```math
ϕ_i(h_i, λ_i) = λ_i^⊤ φ_i(h_i) + \frac \rho2 \| φ_i(h_i) \|_2^2
```
with $φ_i$ a function to compute the current infeasibility
```math
φ_i(h_i, λ_i) = \max\{0 , λ_i + ρ (h_i - h_i^♯)   \} + \min\{0 , λ_i + ρ (h_i - h_i^♭)   \}
```

### Attributes

* `inner::Evaluator`: original problem.
* `cons_type`: type of the constraints of the original problem (equalities or inequalities).
* `cons::VT`: a buffer storing the current evaluation of the constraints for the inner evaluator.
* `rho::T`: current penalty.
* `λ::VT`: current multiplier.
* `scaler::MaxScaler{T,VT}`: a scaler to rescale the range of the constraints in the original problem.

"""
mutable struct AugLagEvaluator{Evaluator<:AbstractNLPEvaluator, T, VT} <: AbstractPenaltyEvaluator
    inner::Evaluator
    # Type
    cons_type::CONSTRAINTS_TYPE
    cons::VT
    ρ::T
    λ::VT
    λc::VT
    # Scaling
    scaler::MaxScaler{T, VT}
    # Stats
    counter::NLPCounter
end
function AugLagEvaluator(
    nlp::AbstractNLPEvaluator, u0;
    scale=false, c₀=0.1,
)
    if !is_constrained(nlp)
        error("Model specified in `nlp` is unconstrained. Instead of using" *
              " an Augmented Lagrangian algorithm, you could use any "*
              "bound constrained solver instead.")
    end
    cons_type = Val(constraints_type(nlp))
    if !isa(cons_type, CONSTRAINTS_TYPE)
        error("Constraints $(constraints_type(nlp)) is not supported by" *
              " AugLagEvaluator.")
    end

    g_min, g_max = bounds(nlp, Constraints())
    cx = similar(g_min) ; fill!(cx, 0)
    λc = similar(g_min) ; fill!(λc, 0)
    λ = similar(g_min) ; fill!(λ, 0)

    scaler = scale ?  MaxScaler(nlp, u0) : MaxScaler(g_min, g_max)
    return AugLagEvaluator(nlp, cons_type, cx, c₀, λ, λc, scaler, NLPCounter())
end
function AugLagEvaluator(
    datafile::String; options...
)
    nlp = ReducedSpaceEvaluator(datafile)
    u0 = initial(nlp)
    return AugLagEvaluator(nlp, u0; options...)
end

has_hessian(nlp::AugLagEvaluator) = has_hessian(nlp.inner)

# Default fallback
function _update_internal!(ag::AugLagEvaluator, ::CONSTRAINTS_TYPE)
    # Update (shifted) infeasibility error
    g♭ = ag.scaler.g_min
    g♯ = ag.scaler.g_max
    ag.λc .= max.(0, ag.λ .+ ag.ρ .* (ag.cons .- g♯)) .+
             min.(0, ag.λ .+ ag.ρ .* (ag.cons .- g♭))
    ag.cons .= max.(0, ag.cons .- g♯) .+ min.(0, ag.cons .- g♭)
end

# Specialization for equality constraints
function _update_internal!(ag::AugLagEvaluator, ::Val{:equality})
    ag.λc .= ag.λ .+ ag.ρ .* ag.cons
end

function update!(ag::AugLagEvaluator, u)
    conv = update!(ag.inner, u)
    # Update constraints
    constraint!(ag.inner, ag.cons, u)
    # Rescale
    ag.cons .*= ag.scaler.scale_cons
    _update_internal!(ag, ag.cons_type)
    return conv
end

function update_penalty!(ag::AugLagEvaluator; η=10.0)
    ag.ρ = min(η * ag.ρ, 10e12)
end

function update_multipliers!(ag::AugLagEvaluator)
    ag.λ .= ag.λc
    return
end

function objective(ag::AugLagEvaluator, u)
    ag.counter.objective += 1
    base_nlp = ag.inner
    cx = ag.cons
    obj = ag.scaler.scale_obj * objective(base_nlp, u) +
        0.5 * ag.ρ * dot(cx, cx) + dot(ag.λ, cx)
    return obj
end

function inner_objective(ag::AugLagEvaluator, u)
    return ag.scaler.scale_obj * objective(ag.inner, u)
end

function gradient!(ag::AugLagEvaluator, grad, u)
    ag.counter.gradient += 1
    base_nlp = ag.inner
    scaler = ag.scaler
    σ = scaler.scale_obj
    mask = abs.(ag.cons) .> 0
    v = scaler.scale_cons .* ag.λc .* mask
    ojtprod!(base_nlp, grad, u, σ, v)
    return
end

function hessprod!(ag::AugLagEvaluator, hessvec, u, w)
    ag.counter.hprod += 1
    scaler = ag.scaler
    cx = ag.cons
    mask = abs.(cx) .> 0

    σ = scaler.scale_obj
    y = (scaler.scale_cons .* ag.λc .* mask)
    ρ = ag.ρ .* (scaler.scale_cons .* scaler.scale_cons .* mask)

    hessian_lagrangian_penalty_prod!(ag.inner, hessvec, u, y, σ, w, ρ)::Nothing
    return
end

function hessian!(ag::AugLagEvaluator, H, u)
    ag.counter.hessian += 1
    scaler = ag.scaler
    cx = ag.cons
    mask = abs.(cx) .> 0

    σ = scaler.scale_obj
    y = (scaler.scale_cons .* ag.λc .* mask)
    ρ = ag.ρ .* (scaler.scale_cons .* scaler.scale_cons .* mask)
    hessian_lagrangian_penalty!(ag.inner, H, u, y, σ, ρ)
    return
end

function estimate_multipliers(ag::AugLagEvaluator, u)
    J = Diagonal(ag.scaler.scale_cons) * jacobian(ag.inner, u)
    ∇f = gradient(ag.inner, u)
    ∇f .*= ag.scaler.scale_obj
    λ = - (J * J') \ (J * ∇f)
    return λ
end

function reset!(ag::AugLagEvaluator)
    reset!(ag.inner)
    empty!(ag.counter)
    fill!(ag.cons, 0)
    fill!(ag.λ, 0)
    fill!(ag.λc, 0)
end

ext: jl


text:
# included in runtests.jl

# Test sending / receiving messages correctly,
# closing from within websocket handlers,
# symmetry of client and server side websockets,
# stress tests opening and closing a sequence of servers.
# At this time, we unfortunately get irritating messages
# 'Workqueue inconsistency detected:...'
import Sockets: IPAddr,
                InetAddr,
                IPv4

include("logformat.jl")
if !@isdefined SUBPROTOCOL
    const SUBPROTOCOL = "Server start the conversation"
    const SUBPROTOCOL_CLOSE = "Server start the conversation and close it from within websocket handler"
end
addsubproto(SUBPROTOCOL)
addsubproto(SUBPROTOCOL_CLOSE)
if !@isdefined(PORT)
    const PORT = 8000
    const SURL = "127.0.0.1"
    const EXTERNALWSURI = "ws://echo.websocket.org"
    const EXTERNALHTTP = "http://httpbin.org/ip"
    const MSGLENGTHS = [0 , 125, 126, 127, 2000]
end

@info "External server http request"
@test 200 == WebSockets.HTTP.request("GET", EXTERNALHTTP).status

@info "String to IPAddr conversion"
let
    ip = parse(IPAddr, SURL)
    @test WebSockets.to_ipaddr(SURL) == WebSockets.to_ipaddr(ip) == ip
end

@info "ServerWS: Open, http response, close. Repeat three times. Takes a while."
for i = 1:3
    let
        ip = parse(IPAddr, SURL)
        server = startserver(url=ip)
        @test 200 == WebSockets.HTTP.request("GET", "http://$SURL:$PORT").status
        close(server)
    end
end

@info "ServerWS: Client side initates message exchange."
let
    server = startserver()
    WebSockets.open(initiatingws, "ws://$SURL:$PORT")
    close(server)
end

@info "ServerWS: Server side initates message exchange."
let
    server = startserver()
    WebSockets.open(echows, "ws://$SURL:$PORT", subprotocol = SUBPROTOCOL)
    close(server)
end

@info "ServerWS: Server side initates message exchange. Close from within server side handler."
let
    server = startserver()
    WebSockets.open(echows, "ws://$SURL:$PORT", subprotocol = SUBPROTOCOL_CLOSE)
    close(server)
end
nothing


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using Statistics
using MAT

mmesh = Mesh(joinpath(PDATA, "twoholes_large.stl"))
mmesh = Mesh(mmesh.nodes * 10, mmesh.elems, -1, 2)

ν = Variable(0.01)
K = ν*constant(compute_fem_laplace_matrix(mmesh))
B = constant(compute_interaction_matrix(mmesh))
Z = [K -B'
    -B spzero(size(B,1))]

bd1 = bcnode((x,y)->y>0.23-1e-5, mmesh)
bd2 = bcnode((x,y)->y<=0.23-1e-5, mmesh)
bd = [bd1;bd2]
bd = [bd; bd .+ mmesh.ndof; 2mmesh.ndof + 1]

rhs = zeros(2mmesh.ndof + mmesh.nelem)
bdval = zeros(length(bd))
bdval[1:length(bd1)] .= 1.0
Z, rhs = impose_Dirichlet_boundary_conditions(Z, rhs, bd, bdval)
sol = Z\rhs
U = matread("fenics/fwd1.mat")["u"]

loss = sum((sol[2mmesh.ndof+1:end] - U[2mmesh.ndof+1:end])^2)
sess = Session(); init(sess)

BFGS!(sess, loss)


text:
@testset "EquivalenceNot" begin
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Int}()))
    config = MOIT.Config()

    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [1, 2])),
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [1, 2])),
    )
    COIDT.equivalencenottest(mock, config)
end


text:
# This file demonstrates an implementation of FLAC support for Onda.jl. Note
# that it's a naive implementation - it just shells out and assumes you have
# the `flac` command line utility installed and available on your system.

using Onda, Test, Random, Dates

#####
##### FLAC
#####

"""
     FLAC(lpcm::LPCM; sample_rate, level=5)
     FLAC(signal::Signal; level=5)

Return a `FLAC<:AbstractLPCMSerializer` instance that represents the
FLAC format assumed for signal files with the ".flac" extension.

The `sample_rate` keyword argument corresponds to `flac`'s `--sample-rate` flag,
while `level` corresponds to `flac`'s `--compression-level` flag.

Note that FLAC is only applicable for signals where `1 <= signal.channel_count <= 8`
and `sizeof(signal.sample_type) in (1, 2)`. The corresponding `signal.file_options`
value may be either `nothing` or `Dict(:level => i)` where `0 <= i <= 8`.

See https://xiph.org/flac/ for details about FLAC.

See also: [`Zstd`](@ref)
"""
struct FLAC{S} <: Onda.AbstractLPCMSerializer
    lpcm::LPCM{S}
    sample_rate::Int
    level::Int
    function FLAC(lpcm::LPCM{S}; sample_rate, level=5) where {S}
        sizeof(S) in (1, 2) || throw(ArgumentError("bit depth must be 8 or 16"))
        1 <= lpcm.channel_count <= 8 || throw(ArgumentError("channel count must be between 1 and 8"))
        return new{S}(lpcm, sample_rate, level)
    end
end

FLAC(signal::Signal; kwargs...) = FLAC(LPCM(signal); sample_rate=signal.sample_rate,
                                       kwargs...)

Onda.serializer_constructor_for_file_extension(::Val{:flac}) = FLAC

function flac_raw_specification_flags(serializer::FLAC{S}) where {S}
    return (level="--compression-level-$(serializer.level)",
            endian="--endian=little",
            channels="--channels=$(serializer.lpcm.channel_count)",
            bps="--bps=$(sizeof(S) * 8)",
            sample_rate="--sample-rate=$(serializer.sample_rate)",
            is_signed=string("--sign=", S <: Signed ? "signed" : "unsigned"))
end

function Onda.deserialize_lpcm(io::IO, serializer::FLAC, args...)
    flags = flac_raw_specification_flags(serializer)
    command = pipeline(`flac - --totally-silent -d --force-raw-format $(flags.endian) $(flags.is_signed)`; stdin=io)
    return open(lpcm_io -> deserialize_lpcm(lpcm_io, serializer.lpcm, args...), command, "r")
end

function Onda.serialize_lpcm(io::IO, samples::AbstractMatrix, serializer::FLAC)
    flags = flac_raw_specification_flags(serializer)
    command = pipeline(`flac --totally-silent $(flags) -`; stdout=io)
    return open(lpcm_io -> serialize_lpcm(lpcm_io, samples, serializer.lpcm), command, "w")
end

#####
##### tests
#####

if VERSION >= v"1.1.0"
    @testset "FLAC example" begin
        signal = Signal([:a, :b, :c], Nanosecond(0), Nanosecond(0), :unit, 0.25, 0.0, Int16, 50.0, :flac, Dict(:level => 2))
        samples = encode(Samples(signal, false, rand(MersenneTwister(1), 3, 50 * 10))).data
        s = serializer(signal)
        bytes = serialize_lpcm(samples, s)
        @test deserialize_lpcm(bytes, s) == samples
        @test deserialize_lpcm(IOBuffer(bytes), s) == samples
        io = IOBuffer(); serialize_lpcm(io, samples, s); seekstart(io)
        @test take!(io) == bytes
    end
else
    @warn "This example may be broken on Julia versions lower than v1.1 due to https://github.com/JuliaLang/julia/issues/33117"
end


text:
using Test
using MPI
using project2

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
size = MPI.Comm_size(comm)

if rank == 0
    g = Grid(6, 6)
    set_boundary_condition!(g, (x, y) -> 10, sides=("top", "bottom"))
else
    g = nothing
end

my_g = partition_grid(g, comm)
solve!(my_g, comm, quiet=true)

if size == 2
    if rank == 0
        ans = [10.0 10.0 10.0 10.0;
               0.0 4.999999999999994 6.666666666666657 6.666666666666653;
               0.0 3.333333333333324 4.999999999999988 4.999999999999982;
               0.0 3.3333333333333255 4.999999999999989 4.999999999999986;
               0.0 4.999999999999996 6.666666666666661 6.666666666666658;
               10.0 10.0 10.0 10.0]
        @test all(isapprox.(my_g.u, ans, atol=1.0e-3))
    elseif rank == 1
        ans = [10.0 10.0 10.0 10.0;
               6.666666666666653 6.666666666666657 4.999999999999996 0.0;
               4.999999999999982 4.999999999999988 3.3333333333333273 0.0;
               4.999999999999986 4.999999999999989 3.3333333333333277 0.0;
               6.666666666666658 6.666666666666661 4.999999999999997 0.0;
               10.0 10.0 10.0 10.0]
        @test all(isapprox.(my_g.u, ans, atol=1.0e-3))
    end
end

if size == 4
    if rank == 0
        ans = [10.0 10.0 10.0;
               0.0 4.999999999999973 6.666666666666611;
               0.0 3.333333333333297 4.999999999999923;
               0.0 3.333333333333302 4.999999999999935;
               0.0 4.999999999999984 6.666666666666633;
               10.0 10.0 10.0]
        @test all(isapprox.(my_g.u, ans, atol=1.0e-3))
    elseif rank == 1
        ans = [10.0 10.0 10.0 10.0;
               4.999999999999964 6.666666666666626 6.666666666666633 4.999999999999969;
               3.333333333333283 4.999999999999945 4.999999999999954 3.3333333333332904;
               3.3333333333332904 4.999999999999954 4.99999999999996 3.333333333333297;
               4.999999999999977 6.666666666666642 6.666666666666646 4.999999999999981;
               10.0 10.0 10.0 10.0]
        @test all(isapprox.(my_g.u, ans, atol=1.0e-3))
    elseif rank == 2
        ans = [10.0 10.0 10.0;
               6.66666666666662 4.999999999999977 0.0;
               4.999999999999935 3.333333333333302 0.0;
               4.999999999999945 3.333333333333307 0.0;
               6.666666666666638 4.999999999999987 0.0;
               10.0 10.0 10.0]
        @test all(isapprox.(my_g.u, ans, atol=1.0e-3))
    elseif rank == 3
        ans = [10.0 10.0;
               4.999999999999969 0.0;
               3.3333333333332904 0.0;
               3.333333333333297 0.0;
               4.999999999999981 0.0;
               10.0 10.0]
        @test all(isapprox.(my_g.u, ans, atol=1.0e-3))
    end
end

get_solution!(g, my_g, comm)

if rank == 0
    ans = [10.0 10.0 10.0 10.0 10.0 10.0;
           0.0 4.999999999999994 6.666666666666657 6.666666666666657 4.999999999999996 0.0;
           0.0 3.333333333333324 4.999999999999988 4.999999999999988 3.3333333333333273 0.0;
           0.0 3.3333333333333255 4.999999999999989 4.999999999999989 3.3333333333333277 0.0;
           0.0 4.999999999999996 6.666666666666661 6.666666666666661 4.999999999999997 0.0;
           10.0 10.0 10.0 10.0 10.0 10.0]
    @test all(isapprox.(g.u, ans, atol=1.0e-3))
end

MPI.Finalize()


text:
#!/usr/bin/env julia

#Start Test Script
using TinySparse
using Base.Test

function test_pack()
    println("Testing sparse matrix pack...")
    pack(sprand(100,100,0.1))
    return true
end

function test_spmv()
    println("Testing out of place multiplication.")
    A = sprand(100,100,0.1)
    spA = pack(A)
    x = rand(100)
    if A*x == spA*x
        return true
    else
        return false
    end
end

function test_spmv_ip()
    println("Testing in place multiplication.")
    A = sprand(100,100,0.1)
    spA = pack(A)
    x = rand(100)
    b = Vector{Float64}(100)
    if A_mul_B!(b, A, x) == A_mul_B!(b, A, x)
        return true
    else
        return false
    end
end

function test_spmv_t_ip()
    println("Testing in place transposed multiplication.")
    A = sprand(100,100,0.1)
    spA = pack(A)
    x = rand(100)
    b = Vector{Float64}(100)
    if Ac_mul_B!(b, A, x) == Ac_mul_B!(b, A, x)
        return true
    else
        return false
    end
end

function run_tests()
    @test test_pack()
    @test test_spmv()
    @test test_spmv_ip()
    @test test_spmv_t_ip()
end

run_tests()


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

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion <: SwaggerModel
    conversionReviewVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: conversionReviewVersions
    strategy::Any # spec type: Union{ Nothing, String } # spec name: strategy
    webhookClientConfig::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1WebhookClientConfig } # spec name: webhookClientConfig

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion(;conversionReviewVersions=nothing, strategy=nothing, webhookClientConfig=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion, Symbol("conversionReviewVersions"), conversionReviewVersions)
        setfield!(o, Symbol("conversionReviewVersions"), conversionReviewVersions)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion, Symbol("strategy"), strategy)
        setfield!(o, Symbol("strategy"), strategy)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion, Symbol("webhookClientConfig"), webhookClientConfig)
        setfield!(o, Symbol("webhookClientConfig"), webhookClientConfig)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion = Dict{Symbol,Symbol}(Symbol("conversionReviewVersions")=>Symbol("conversionReviewVersions"), Symbol("strategy")=>Symbol("strategy"), Symbol("webhookClientConfig")=>Symbol("webhookClientConfig"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion = Dict{Symbol,String}(Symbol("conversionReviewVersions")=>"Vector{String}", Symbol("strategy")=>"String", Symbol("webhookClientConfig")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1WebhookClientConfig")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion)
    (getproperty(o, Symbol("strategy")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1CustomResourceConversion }, name::Symbol, val)
end


text:
using Pathfinder
using Documenter

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

makedocs(;
    modules=[Pathfinder],
    authors="Seth Axen <seth.axen@gmail.com> and contributors",
    repo="https://github.com/sethaxen/Pathfinder.jl/blob/{commit}{path}#{line}",
    sitename="Pathfinder.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://sethaxen.github.io/Pathfinder.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Library" => [
            "Public" => "lib/public.md",
            "Internals" => "lib/internals.md",
        ],
        "Examples" => [
            "Quickstart" => "examples/quickstart.md",
            "Initializing HMC" => "examples/initializing-hmc.md",
            "Turing usage" => "examples/turing.md",
        ]
    ],
)

deploydocs(; repo="github.com/sethaxen/Pathfinder.jl", devbranch="main")


text:
using Ensemble
using CARMA
using HDF5

BLAS.set_num_threads(1)

usage = "julia run_multi_carma.jl P Q [NLIVE] FILE1 FILE2 ..."

if size(ARGS, 1) < 4
    println(usage)
    exit(1)
end

p = parse(Int, ARGS[1])
q = parse(Int, ARGS[2])

nlive = 1024
nlive_given = false
try
    nlive = parse(Int, ARGS[3])
    nlive_given = true
catch
    # Do nothing
end

files = ARGS[3:end]
if nlive_given
    files = ARGS[4:end]
end

nmcmc = 128

outfile = "state-$(p)-$(q).hdf5"
ckpt_file = "state-$(p)-$(q).ckpt"

ts = Array{Float64, 1}[]
ys = Array{Float64, 1}[]
dys = Array{Float64, 1}[]
for f in files
    data = readdlm(f)

    t = data[:,1]
    y = data[:,2]
    dy = data[:,3]

    # Sorted in time
    inds = sortperm(t)
    t = t[inds]
    y = y[inds]
    dy = dy[inds]

    # Get rid of NaNs
    sel = ~isnan(y)
    t = t[sel]
    y = y[sel]
    dy = dy[sel]

    # Set dy == 0 to minimum positive dy
    @assert any(dy .> 0)
    sel = dy .== 0
    dy[sel] = minimum(dy[~sel])

    push!(ts, t)
    push!(ys, y)
    push!(dys, dy)
end

post = Kalman.MultiSegmentCARMAKalmanPosterior(ts, ys, dys, p, q)
function logl(x)
    Kalman.log_likelihood(post, x)
end
function logp(x)
    Kalman.log_prior(post, x)
end

if ispath(ckpt_file)
    nest_state = h5open(ckpt_file, "r") do f
        EnsembleNest.NestState(f, logl=logl, logp=logp)
    end
else
    nest_state = EnsembleNest.NestState(logl, logp, Kalman.init(post, nlive), nmcmc)
end

EnsembleNest.run!(nest_state, 0.01, verbose=true, ckpt_file=ckpt_file)

h5open(outfile, "w") do f
    ng = g_create(f, "nest_state")
    write(ng, nest_state)

    pg = g_create(f, "carma_data")
    pg["p"] = p
    pg["q"] = q
    pg["n"] = size(ts, 1)
    for i in eachindex(ts)
        gi = g_create(pg, "$(i)")
        gi["ts", "compress", 3, "shuffle", ()] = ts[i]
        gi["ys", "compress", 3, "shuffle", ()] = ys[i]
        gi["dys", "compress", 3, "shuffle", ()] = dys[i]
    end
end

if ispath(ckpt_file)
    rm(ckpt_file)
end


text:
###############################################################################
# Damper Force
###############################################################################
function damper_force(relative::Symbol, joint::Rotational{T}, qa::Quaternion, ϕa::AbstractVector,
        qb::Quaternion, ϕb::AbstractVector, timestep;
        rotate::Bool=true,
        unitary::Bool=false) where T

    damper = unitary ? 1.0 : joint.damper
    A = nullspace_mask(joint)
    Aᵀ = zerodimstaticadjoint(A)
    axis_offset = joint.axis_offset

    velocity = minimal_velocities(joint, szeros(3), szeros(3), qa, ϕa, szeros(3), szeros(3), qb, ϕb, timestep)
    if relative == :parent
        force = 1.0000 * damper * Aᵀ * velocity # currently assumes same damper constant in all directions
        rotate && (force = vector_rotate(force, axis_offset)) # rotate back to frame a
        return [szeros(T, 3); force]
    elseif relative == :child
        force = - 1.0000 * damper * Aᵀ * velocity # currently assumes same damper constant in all directions
        rotate && (force = vector_rotate(force, inv(qb) * qa * axis_offset)) # rotate back to frame b
        return [szeros(T, 3); force]
    end
end

damper_impulses(relative::Symbol, joint::Rotational, pbody::Node, cbody::Node, timestep; unitary::Bool=false) =
    timestep * damper_force(relative, joint, current_configuration(pbody.state)[2], pbody.state.ϕsol[2],
    current_configuration(cbody.state)[2], cbody.state.ϕsol[2], timestep, unitary=unitary)

damper_impulses(relative::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep; unitary::Bool=false) where T = szeros(T, 6)

###############################################################################
# Damper Jacobians
################################################################################

function damper_jacobian_configuration(relative::Symbol, jacobian::Symbol,
        joint::Rotational, pbody::Node, cbody::Node,
        timestep::T;
        ) where T

    Aᵀ = zerodimstaticadjoint(nullspace_mask(joint))
    xa, va, qa, ϕa = current_configuration_velocity(pbody.state)
    xb, vb, qb, ϕb = current_configuration_velocity(cbody.state)
    axis_offset = joint.axis_offset
    X = szeros(T, 3, 3)
    Z = szeros(T, 3, 6)

    force = damper_force(relative, joint, qa, ϕa, qb, ϕb, timestep; rotate = false)[SUnitRange(4,6)]
    ∂vel = minimal_velocities_jacobian_configuration(jacobian, joint, xa, va, qa, ϕa, xb, vb, qb, ϕb, timestep)[:, SUnitRange(4,6)]

    if relative == :parent
        Q = rotation_matrix(axis_offset) * 1.0000 * joint.damper * Aᵀ * ∂vel
    elseif relative == :child
        Q = rotation_matrix(inv(qb) * qa * axis_offset) * -1.0000 * joint.damper * Aᵀ * ∂vel
        if jacobian == :parent
            Q += rotation_matrix(inv(qb)) * ∂rotation_matrix∂q(qa, vector_rotate(force, axis_offset), attjac=true)
        elseif jacobian == :child
            Q += ∂rotation_matrix_inv∂q(qb, vector_rotate(force, qa * axis_offset), attjac=true)
        end
    end
    return timestep * [Z; X Q]
end

function damper_jacobian_velocity(relative::Symbol, jacobian::Symbol,
    joint::Rotational, pbody::Node, cbody::Node,
    timestep::T) where T

    Aᵀ = zerodimstaticadjoint(nullspace_mask(joint))
    xa, va, qa, ϕa = current_configuration_velocity(pbody.state)
    xb, vb, qb, ϕb = current_configuration_velocity(cbody.state)
    axis_offset = joint.axis_offset

    force = damper_force(relative, joint, qa, ϕa, qb, ϕb, timestep; rotate = false)[SUnitRange(4,6)]
    ∂vel = minimal_velocities_jacobian_velocity(jacobian, joint, xa, va, qa, ϕa, xb, vb, qb, ϕb, timestep)

    if relative == :parent
        VΩ = rotation_matrix(axis_offset) * 1.0000 * joint.damper * Aᵀ * ∂vel
    elseif relative == :child
        VΩ = rotation_matrix(inv(qb) * qa * axis_offset) * -1.0000 * joint.damper * Aᵀ * ∂vel
    end
    return timestep * [szeros(T, 3, 6); VΩ]
end

damper_jacobian_configuration(relative::Symbol, jacobian::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep::T; attjac::Bool = true) where T = attjac ? szeros(T, 6, 6) : szeros(T, 6, 7)
damper_jacobian_velocity(relative::Symbol, jacobian::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep::T) where T = szeros(T, 6, 6)

#
#
#
# using Test
# using Plots
# using FiniteDiff
# vis = Visualizer()
# open(vis)
#
#
#
# mech = Dojo.get_snake(gravity=0.00, num_bodies=2, damper=0.3, joint_type=:Revolute)
# Dojo.initialize_snake!(mech)
# function ctrl!(m,k)
#     set_input!(m, 0.01*m.timestep*ones(minimal_dimension(m)))
# end
# storage = Dojo.simulate!(mech, 1.0, ctrl!)
# Dojo.visualize(mech, storage, vis=vis)
#
#
# rot0 = mech.joints[2].rotational
# timestep0 = mech.timestep
# pbody0 = mech.bodies[1]
# cbody0 = mech.bodies[2]
# xa0, va0, qa0, ϕa0 = Dojo.current_configuration_velocity(pbody0.state)
# xb0, vb0, qb0, ϕb0 = Dojo.current_configuration_velocity(cbody0.state)
#
# # Configuration
# J0 = Dojo.damper_jacobian_configuration(:parent, :parent, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:parent, rot0,
#         Dojo.Quaternion(xq[4:7]...,true), ϕa0, qb0, ϕb0, timestep0; rotate=true, unitary=false),
#     [xa0; Dojo.vector(qa0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qa0), dims=(1,2))
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
# J0 = Dojo.damper_jacobian_configuration(:parent, :child, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:parent, rot0, qa0, ϕa0,
#         Dojo.Quaternion(xq[4:7]...,true), ϕb0, timestep0; rotate=true, unitary=false),
#     [xb0; Dojo.vector(qb0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qb0), dims=(1,2))
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
# J0 = Dojo.damper_jacobian_configuration(:child, :parent, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:child, rot0,
#         Dojo.Quaternion(xq[4:7]...,true), ϕa0, qb0, ϕb0, timestep0; rotate=true, unitary=false),
#     [xa0; Dojo.vector(qa0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qa0), dims=(1,2))
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
#
# J0 = Dojo.damper_jacobian_configuration(:child, :child, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:child, rot0, qa0, ϕa0,
#         Dojo.Quaternion(xq[4:7]...,true), ϕb0, timestep0; rotate=true, unitary=false),
#     [xb0; Dojo.vector(qb0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qb0), dims=(1,2))
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
# mech = Dojo.get_snake(gravity=0.00, num_bodies=2, damper=0.3, joint_type=:Revolute)
# Dojo.initialize_snake!(mech)
# function ctrl!(m,k)
#     set_input!(m, 0.01*m.timestep*ones(minimal_dimension(m)))
# end
# storage = Dojo.simulate!(mech, 1.0, ctrl!)
# Dojo.visualize(mech, storage, vis=vis)
#
# mech.joints[2]
# rot0 = mech.joints[2].rotational
# timestep0 = mech.timestep
# pbody0 = mech.bodies[1]
# cbody0 = mech.bodies[2]
# xa0, va0, qa0, ϕa0 = Dojo.current_configuration_velocity(pbody0.state)
# xb0, vb0, qb0, ϕb0 = Dojo.current_configuration_velocity(cbody0.state)
#
# # Configuration
# J0 = Dojo.damper_jacobian_configuration(:parent, :parent, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:parent, rot0, Dojo.Quaternion(xq[4:7]...,true), ϕa0, qb0, ϕb0, timestep0; rotate=true, unitary=false),
#     [xa0; Dojo.vector(qa0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qa0), dims=(1,2))
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
# J0 = Dojo.damper_jacobian_configuration(:parent, :child, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:parent, rot0, qa0, ϕa0, Dojo.Quaternion(xq[4:7]...,true), ϕb0, timestep0; rotate=true, unitary=false),
#     [xb0; Dojo.vector(qb0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qb0), dims=(1,2))
# norm(J0 - J1, Inf)
# J0 + J1
# @test norm(J0 - J1, Inf) < 1e-6
#
# J0 = Dojo.damper_jacobian_configuration(:child, :parent, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:child, rot0, Dojo.Quaternion(xq[4:7]...,true), ϕa0, qb0, ϕb0, timestep0; rotate=true, unitary=false),
#     [xa0; Dojo.vector(qa0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qa0), dims=(1,2))
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
# J0
# J1
#
# J0 = Dojo.damper_jacobian_configuration(:child, :child, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> timestep0 * Dojo.damper_force(:child, rot0, qa0, ϕa0, Dojo.Quaternion(xq[4:7]...,true), ϕb0, timestep0; rotate=true, unitary=false),
#     [xb0; Dojo.vector(qb0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qb0), dims=(1,2))
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
# J0
# J1
#
#
#
# timestep0 = 0.01
# xa0 = szeros(3)
# qa0 = Quaternion(1,0,0,0.0,true)
# va0 = szeros(3)
# ϕa0 = szeros(3)
# xb0 = szeros(3)
# qb0 = Quaternion(1,0,0,0.0,true)
# vb0 = szeros(3)
# ϕb0 = szeros(3)
# J0 = minimal_velocities_jacobian_velocity(:parent,
#     rot0, xa0, va0, qa0, ϕa0, xb0, vb0, qb0, ϕb0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     xq -> Dojo.minimal_velocities(rot0, xq[SUnitRange(1,3)], va0, Dojo.Quaternion(xq[4:7]...,true), ϕa0,
#         xb0, vb0, qb0, ϕb0, timestep0),
#     [xa0; Dojo.vector(qa0)]) * Dojo.cat(I(3), Dojo.LVᵀmat(qa0), dims=(1,2))
# norm(J0 - J1, Inf)
#
# using Rotations
#
#
#
#
#
#
#
#
#
#
#
#
# q0 = rand(QuatRotation).q
# J0 = drotation_vectordq(q0)
# J1 = FiniteDiff.finite_difference_jacobian(q -> rotation_vector(q), vector(q0))
# norm(J0 - J1, Inf)
#
# timestep0 = 0.01
# q0 = rand(QuatRotation).q
# ϕ0 = srand(3)
# J0 = rotational_integrator_jacobian_velocity(q0, -ϕ0, timestep0)
# J1 = -FiniteDiff.finite_difference_jacobian(ϕ -> vector(next_orientation(q0, -ϕ, timestep0)), ϕ0)
# norm(J0 - J1, Inf)
#
#
#
#
#
#
# function minvel(joint::Rotational, qb::Quaternion, ϕb::AbstractVector, timestep)
# 	A = nullspace_mask(joint)
# 	qb1 = next_orientation(qb, -ϕb, timestep)
#     return A * rotation_vector(qb1)
# end
#
# function minvel_jacobian_velocity(joint::Rotational{T}, qb::Quaternion,
#         ϕb::AbstractVector, timestep) where T
#
# 	A = nullspace_mask(joint)
#     qb1 = next_orientation(qb, -ϕb, timestep)
#     Ω = A * drotation_vectordq(qb1) * -rotational_integrator_jacobian_velocity(qb, -ϕb, timestep)
# 	return Ω
# end
#
# LVᵀmat(qb0) * VLᵀmat(qb0)
#
# timestep0 = 1.0
# xa0 = srand(3)
# # qa0 = rand(QuatRotation).q
# qa0 = Quaternion(1,0,0,0.0,true)
# va0 = srand(3)
# ϕa0 = srand(3)
# xb0 = srand(3)
# qb0 = rand(QuatRotation).q
# # qb0 = Quaternion(1,0,0,0.0,true)
# vb0 = srand(3)
# ϕb0 = srand(3)
# rot0
# J0 = minvel_jacobian_velocity(rot0, qb0, ϕb0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     ϕ -> Dojo.minvel(rot0, qb0, ϕ, timestep0),
#     ϕb0)
# norm(J0 - J1, Inf)
#
#
#
#
#
#
#
#
#
#
# # Velocity
# J0 = Dojo.damper_jacobian_velocity(:parent, :parent, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     vϕ -> timestep0 * Dojo.damper_force(:parent, rot0, qa0, vϕ[Dojo.SUnitRange(4,6)], qb0, ϕb0, timestep0; rotate=true, unitary=false),
#     [va0; ϕa0])
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
# J0 = Dojo.damper_jacobian_velocity(:parent, :child, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     vϕ -> timestep0 * Dojo.damper_force(:parent, rot0, qa0, ϕa0, qb0, vϕ[Dojo.SUnitRange(4,6)], timestep0; rotate=true, unitary=false),
#     [vb0; ϕb0])
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
# J0 = Dojo.damper_jacobian_velocity(:child, :parent, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     vϕ -> timestep0 * Dojo.damper_force(:child, rot0, qa0, vϕ[Dojo.SUnitRange(4,6)], qb0, ϕb0, timestep0; rotate=true, unitary=false),
#     [va0; ϕa0])
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
# J0 = Dojo.damper_jacobian_velocity(:child, :child, rot0, pbody0, cbody0, timestep0)
# J1 = FiniteDiff.finite_difference_jacobian(
#     vϕ -> timestep0 * Dojo.damper_force(:child, rot0, qa0, ϕa0, qb0, vϕ[Dojo.SUnitRange(4,6)], timestep0; rotate=true, unitary=false),
#     [vb0; ϕb0])
# norm(J0 - J1, Inf)
# @test norm(J0 - J1, Inf) < 1e-6
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#


text:
using Documenter, Connectomes

makedocs(
	sitename="Connectomes.jl",
	modules = [Connectomes],
	pages = [
	"Home" => "index.md",
	"Connectomes" => "connectomes.md"]
)

deploydocs(;
    repo = "github.com/PavanChaggar/Connectomes.jl.git",
	devbranch="main",
)


text:
"""
    is_option(x)

Check if `x` is an option type or not.
"""
is_option(x) = false

"""
    compare_options(a, b, xs...)::Bool

Compare option types check if they are the same.
"""
function compare_options(a, b, xs...)::Bool
    compare_options(a, b) || return false
    compare_options(b, xs...)
end

compare_options(a, b) = false

function compare_options(a::A, b::A) where {A}
    is_option(A) || return a == b # fall through to normal compare
    for idx in 1:nfields(a)
        compare_options(getfield(a, idx), getfield(b, idx)) || return false
    end
    return true
end

function tryparse_jltype(s)
    s isa String || throw(ArgumentError("expect type String, got: $(typeof(s))"))
    type_ex = Meta.parse(s)
    is_datatype_expr(type_ex) || throw(ArgumentError("expect type expression got: $type_ex"))
    try
        return eval(type_ex)
    catch
        return
    end
end

function parse_jltype(s)
    type = tryparse_jltype(s)
    type === nothing && throw(ArgumentError("cannot parse $s to a Julia type"))
    return type
end

# NOTE: copied from JLD
# https://github.com/JuliaIO/JLD.jl/blob/83ea0c5ef7293c78d7d9c8ffdf9ede599b54dc4c/src/JLD00.jl#L991
# we only have DataType to serialize
function full_typename(jltype::DataType)
    tname = string(jltype.name.module, ".", jltype.name.name)
    if isempty(jltype.parameters)
        return tname
    else
        params_str = join([full_typename(x) for x in jltype.parameters], ",")
        return string(tname, "{", params_str, "}")
    end
end

function contains_reflect_type(::Type{T}) where T
    for idx in 1:fieldcount(T)
        Reflect === fieldtype(T, idx) && return true
    end
    return false
end

function is_union_of_multiple_options(::Type{T}) where T
    T isa Union || return false
    T.a === Nothing && return is_union_of_multiple_options(T.b)
    T.b === Nothing && return is_union_of_multiple_options(T.a)

    # not option type
    return is_option_maybe(T.a) && is_option_maybe(T.b)
end


text:
T_FWHM = 0.05ps # pulse time paramater passed as FWHM
betacoeff = [0.0/m, 0.0ps/m, -0.12ps^2/m, 0.0ps^3/m, 0.005e-3ps^4/m]
γnl = 1.0W/m
ω0 = 2π*c0/(1550nm) |> THz
P_pulse = 10.0W
fiber_length = 0.020m
nlo = ToyModel(;ω0=ω0, betacoeff=betacoeff,
                   nonlinearcoeff=γnl,
                   linearloss=0.0/m,
                   power_in=P_pulse,
                   pulsetime=T_FWHM,
                   length=fiber_length,
                   has_shock=true,
                   has_raman=true)

sol_nlo = solve(nlo, DynamicNLSE();
                tpoints=2^13, time_window=10.0ps)


text:
@testset "loglikelihoods.jl" begin
    for m in DynamicPPL.TestUtils.DEMO_MODELS
        vi = VarInfo(m)

        vns = vi.metadata.m.vns
        if length(vns) == 1 && length(vi[vns[1]]) == 1
            # Only have one latent variable.
            DynamicPPL.setval!(vi, [1.0], ["m"])
        else
            DynamicPPL.setval!(vi, [1.0, 1.0], ["m[1]", "m[2]"])
        end

        lls = pointwise_loglikelihoods(m, vi)

        if isempty(lls)
            # One of the models with literal observations, so we just skip.
            continue
        end

        loglikelihood = if length(keys(lls)) == 1 && length(m.args.x) == 1
            # Only have one observation, so we need to double it
            # for comparison with other models.
            2 * sum(lls[first(keys(lls))])
        else
            sum(sum, values(lls))
        end

        @test loglikelihood ≈ -324.45158270528947
    end
end

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
using Statistics

using Oceananigans.Fields: CenterField, ZFaceField, compute_at!
using Oceananigans.Grids: halo_size

@testset "Averaged fields" begin
    @info "Testing averaged fields..."

    for arch in archs
        @testset "Averaged fields [$(typeof(arch))]" begin
            @info "  Testing AveragedFields [$(typeof(arch))]"
            for FT in float_types

                grid = RegularRectilinearGrid(topology = (Periodic, Periodic, Bounded),
                                                size = (2, 2, 2),
                                                   x = (0, 2), y = (0, 2), z = (0, 2))

                w = ZFaceField(arch, grid)
                T = CenterField(arch, grid)

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

                set!(T, trilinear)
                set!(w, trilinear)

                @compute T̃ = AveragedField(T, dims=(1, 2, 3))

                # Note: halo regions must be *filled* prior to computing an average
                # if the average within halo regions is to be correct.
                fill_halo_regions!(T, arch)
                @compute T̅ = AveragedField(T, dims=(1, 2))

                fill_halo_regions!(T, arch)
                @compute T̂ = AveragedField(T, dims=1)

                @compute w̃ = AveragedField(w, dims=(1, 2, 3))
                @compute w̅ = AveragedField(w, dims=(1, 2))
                @compute ŵ = AveragedField(w, dims=1)

                Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz

                @test T̃[1, 1, 1] ≈ 3

                @test Array(interior(T̅))[1, 1, :] ≈ [2.5, 3.5]
                @test Array(interior(T̂))[1, :, :] ≈ [[2, 3] [3, 4]]

                @test w̃[1, 1, 1] ≈ 3

                @test Array(interior(w̅))[1, 1, :] ≈ [2, 3, 4]
                @test Array(interior(ŵ))[1, :, :] ≈ [[1.5, 2.5] [2.5, 3.5] [3.5, 4.5]]
            end
        end

        @testset "Conditional computation of AveragedFields [$(typeof(arch))]" begin
            @info "  Testing conditional computation of AveragedFields [$(typeof(arch))]"
            for FT in float_types
                grid = RegularRectilinearGrid(FT, size=(2, 2, 2), extent=(1, 1, 1))
                c = CenterField(arch, grid)

                for dims in (1, 2, 3, (1, 2), (2, 3), (1, 3), (1, 2, 3))
                    C = AveragedField(c, dims=dims)

                    @test !isnothing(C.status)

                    # Test conditional computation
                    set!(c, 1)
                    compute_at!(C, FT(1)) # will compute
                    @test all(interior(C) .== 1)
                    @test C.status.time == FT(1)

                    set!(c, 2)
                    compute_at!(C, FT(1)) # will not compute because status == 1
                    @test C.status.time == FT(1)
                    @test all(interior(C) .== 1)

                    compute_at!(C, FT(2))
                    @test C.status.time == FT(2)
                    @test all(interior(C) .== 2)
                end
            end
        end
    end
end


text:
type FdmSolverDesc{F <: FdmMesher, C <: FdmInnerValueCalculator}
  mesher::F
  bcSet::FdmBoundaryConditionSet
  condition::FdmStepConditionComposite
  calculator::C
  maturity::Float64
  timeSteps::Int
  dampingSteps::Int
end

## Solvers ##
type FdmBackwardSolver
  map::FdmLinearOpComposite
  bcSet::FdmBoundaryConditionSet
  condition::FdmStepConditionComposite
  schemeDesc::FdmSchemeDesc
end

function FdmBackwardSolver(map::FdmLinearOpComposite, bcSet::FdmBoundaryConditionSet, schemeDesc::FdmSchemeDesc)
  condition = FdmStepConditionComposite(Vector{Float64}(), Vector{StepCondition}())

  return FdmBackwardSolver(map, bcSet, condition, schemeDesc)
end

function rollback!(bsolv::FdmBackwardSolver, schemeType::Hundsdorfer, rhs::Vector{Float64}, from::Float64, to::Float64, steps::Int, dampingSteps::Int)
  deltaT = from - to
  allSteps = steps + dampingSteps
  dampingTo = from - (deltaT * dampingSteps) / allSteps

  # if dampingSteps > 0
  #   implicitEvolver = ImplicitEulerScheme(bsolv.map, bsolv.bcSet)
  # end

  dampingSteps == 0 || error("damping steps shoudl be 0")

  hsEvolver = HundsdorferScheme(bsolv.schemeDesc.theta, bsolv.schemeDesc.mu, bsolv.map, bsolv.bcSet)
  hsModel = FiniteDifferenceModel{HundsdorferScheme}(hsEvolver, bsolv.condition.stoppingTimes)
  rollback!(hsModel, rhs, dampingTo, to, steps, bsolv.condition)

  return bsolv, rhs
end

function rollback!(bsolv::FdmBackwardSolver, schemeType::Douglas, rhs::Vector{Float64}, from::Float64, to::Float64, steps::Int, dampingSteps::Int)
  deltaT = from - to
  allSteps = steps + dampingSteps
  dampingTo = from - (deltaT * dampingSteps) / allSteps

  dampingSteps == 0 || error("damping steps shoudl be 0")

  dsEvolver = DouglasScheme(bsolv.schemeDesc.theta, bsolv.map, bsolv.bcSet)
  dsModel = FiniteDifferenceModel{DouglasScheme}(dsEvolver, bsolv.condition.stoppingTimes)
  rollback!(dsModel, rhs, dampingTo, to, steps, bsolv.condition)

  return bsolv, rhs
end

type Fdm1DimSolver{FD <: FdmLinearOpComposite} <: LazyObject
  lazyMixin::LazyMixin
  solverDesc::FdmSolverDesc
  schemeDesc::FdmSchemeDesc
  op::FD
  thetaCondition::FdmSnapshotCondition
  conditions::FdmStepConditionComposite
  initialValues::Vector{Float64}
  resultValues::Vector{Float64}
  x::Vector{Float64}
  interpolation::NaturalCubicSpline

  function Fdm1DimSolver{FD}(solverDesc::FdmSolverDesc, schemeDesc::FdmSchemeDesc, op::FD)
    thetaCondition = FdmSnapshotCondition(0.99 * min(1.0 / 365.0, length(solverDesc.condition.stoppingTimes) == 0 ? solverDesc.maturity : solverDesc.condition.stoppingTimes[1]))
    conditions = join_conditions_FdmStepConditionComposite(thetaCondition, solverDesc.condition)

    layout = solverDesc.mesher.layout

    x = zeros(layout.size)
    initialValues = zeros(layout.size)
    resultValues = zeros(layout.size)

    coords = ones(Int, length(layout.dim))

    for i = 1:layout.size
      initialValues[i], solverDesc.calculator = avg_inner_value(solverDesc.calculator, coords, i, solverDesc.maturity)
      x[i] = get_location(solverDesc.mesher, coords, 1)

      iter_coords!(coords, layout.dim)
    end

    new(LazyMixin(), solverDesc, schemeDesc, op, thetaCondition, conditions, initialValues, resultValues, x)
  end
end

type Fdm2DimSolver{FD <: FdmLinearOpComposite} <: LazyObject
  lazyMixin::LazyMixin
  solverDesc::FdmSolverDesc
  schemeDesc::FdmSchemeDesc
  op::FD
  thetaCondition::FdmSnapshotCondition
  conditions::FdmStepConditionComposite
  initialValues::Vector{Float64}
  resultValues::Matrix{Float64}
  x::Vector{Float64}
  y::Vector{Float64}
  interpolation::BicubicSpline

  function Fdm2DimSolver{FD}(solverDesc::FdmSolverDesc, schemeDesc::FdmSchemeDesc, op::FD)
    thetaCondition = FdmSnapshotCondition(0.99 * min(1.0 / 365.0, length(solverDesc.condition.stoppingTimes) == 0 ? solverDesc.maturity : solverDesc.condition.stoppingTimes[1]))
    conditions = join_conditions_FdmStepConditionComposite(thetaCondition, solverDesc.condition)

    layout = solverDesc.mesher.layout

    initialValues = zeros(layout.size)
    resultValues = zeros(layout.dim[1], layout.dim[2])

    x = zeros(layout.dim[1])
    y = zeros(layout.dim[2])
    x_count = 1
    y_count = 1

    coords = ones(Int, length(layout.dim))

    for i = 1:layout.size
      initialValues[i], solverDesc.calculator = avg_inner_value(solverDesc.calculator, coords, i, solverDesc.maturity)

      if coords[2] == 1
        x[x_count] = get_location(solverDesc.mesher, coords, 1)
        x_count += 1
      end

      if coords[1] == 1
        y[y_count] =  get_location(solverDesc.mesher, coords, 2)
        y_count += 1
      end

      iter_coords!(coords, layout.dim)
    end

    return new{FD}(LazyMixin(), solverDesc, schemeDesc, op, thetaCondition, conditions, initialValues, resultValues, x, y)
  end
end

get_interpolation(solv::Fdm1DimSolver, x::Float64) = solv.interpolation(x)
get_interpolation(solv::Fdm2DimSolver, x::Float64, y::Float64) = solv.interpolation.spline(x, y)

function interpolate_at(solv::Fdm1DimSolver, x::Float64)
  calculate!(solv)
  return get_interpolation(solv, x)
end

function interpolate_at(solv::Fdm2DimSolver, x::Float64, y::Float64)
  calculate!(solv)
  return get_interpolation(solv, x, y)
end

function perform_calculations!(solv::Fdm1DimSolver)
  rhs = copy(solv.initialValues)

  bsolver = FdmBackwardSolver(solv.op, solv.solverDesc.bcSet, solv.conditions, solv.schemeDesc)
  rollback!(bsolver, solv.schemeDesc.schemeType, rhs, solv.solverDesc.maturity, 0.0, solv.solverDesc.timeSteps, solv.solverDesc.dampingSteps)
  solv.resultValues = copy(rhs)
  solv.interpolation = NaturalCubicSpline(solv.x, solv.resultValues)

  return solv
end

function perform_calculations!(solv::Fdm2DimSolver)
  rhs = copy(solv.initialValues)

  bsolver = FdmBackwardSolver(solv.op, solv.solverDesc.bcSet, solv.conditions, solv.schemeDesc)
  rollback!(bsolver, solv.schemeDesc.schemeType, rhs, solv.solverDesc.maturity, 0.0, solv.solverDesc.timeSteps, solv.solverDesc.dampingSteps)
  solv.resultValues = reshape(rhs, length(solv.x), length(solv.y))
  solv.interpolation = BicubicSpline(solv.x, solv.y, solv.resultValues)

  return solv
end

type FdmG2Solver <: LazyObject
  lazyMixin::LazyMixin
  model::G2
  solverDesc::FdmSolverDesc
  schemeDesc::FdmSchemeDesc
  solver::Fdm2DimSolver

  function FdmG2Solver(model::G2, solverDesc::FdmSolverDesc, schemeDesc::FdmSchemeDesc)
    op = FdmG2Op(solverDesc.mesher, model, 1, 2)
    solver = Fdm2DimSolver{FdmG2Op}(solverDesc, schemeDesc, op)
    new(LazyMixin(), model, solverDesc, schemeDesc, solver)
  end
end

type FdmHullWhiteSolver <: LazyObject
  lazyMixin::LazyMixin
  model::HullWhite
  solverDesc::FdmSolverDesc
  schemeDesc::FdmSchemeDesc
  solver::Fdm1DimSolver

  function FdmHullWhiteSolver(model::HullWhite, solverDesc::FdmSolverDesc, schemeDesc::FdmSchemeDesc)
    op = FdmHullWhiteOp(solverDesc.mesher, model, 1)
    solver = Fdm1DimSolver{FdmHullWhiteOp}(solverDesc, schemeDesc, op)
    new(LazyMixin(), model, solverDesc, schemeDesc, solver)
  end
end

function value_at(solv::FdmG2Solver, x::Float64, y::Float64)
  return interpolate_at(solv.solver, x, y)
end

function value_at(solv::FdmHullWhiteSolver, x::Float64)
  return interpolate_at(solv.solver, x)
end


text:
""" Interoperation with mesh files.

This module enables delta sets to be imported from mesh files supported by
MeshIO.jl and for delta sets to be converted to meshes, mainly for the purposes
of plotting. Meshes are represented by the `GeometryBasics.Mesh` type.
"""
module MeshInterop

using FileIO, MeshIO, GeometryBasics

using Catlab.CategoricalAlgebra: copy_parts!
using ..SimplicialSets, ..DiscreteExteriorCalculus
import ..SimplicialSets: EmbeddedDeltaSet2D

# Export meshes
###############

""" Construct a Mesh object from an embedded delta set.
"""
function GeometryBasics.Mesh(ds::EmbeddedDeltaSet2D)
  points = Point{3, Float64}[point(ds)...]
  tris = TriangleFace{Int}[zip(triangle_vertices(ds)...)...]
  Mesh(points, tris)
end

""" Construct a Mesh object from a dual embedded delta set.
"""
function GeometryBasics.Mesh(ds::EmbeddedDeltaDualComplex2D{O,R,P};
                             primal=false) where {O,R,P}
  if primal
    ds′ = EmbeddedDeltaSet2D{O,P}()
    copy_parts!(ds′, ds)
    Mesh(ds′)
  else
    points = Point{3, Float64}[dual_point(ds)...]
    tris = TriangleFace{Int}[zip(dual_triangle_vertices(ds)...)...]
    Mesh(points, tris)
  end
end

# Import meshes
###############

""" Construct EmbeddedDeltaSet2D from Mesh object

This operator should work for any triangular mesh object. Note that it will
not preserve any normal, texture, or other attributes from the Mesh object.

The `force_unique` flag merges all points which are at the same location. This
is necessary to process `.stl` files which references points by location.
"""
function EmbeddedDeltaSet2D(m::GeometryBasics.Mesh; force_unique=false)
  coords = metafree.(coordinates(m))
  ind_map = 1:length(coords)
  if(force_unique)
    indices = unique(i->coords[i], 1:length(coords))
    val2ind = Dict(coords[indices[i]]=>i for i in 1:length(indices))
    ind_map = map(c->val2ind[c], coords)
    coords = coords[indices]
  end
  tris = faces(m)
  s = EmbeddedDeltaSet2D{Bool, eltype(coords)}()
  add_vertices!(s, length(coords), point=coords)
  for tri in tris
    tri = ind_map[convert.(Int64, tri)]
    glue_sorted_triangle!(s, tri...)
  end
  # Properly orient the delta set
  orient!(s)
  s
end

""" Construct EmbeddedDeltaSet2D from mesh file

This operator should work for any file support for import from MeshIO. Note
that it will not preserve any normal, texture, or other data beyond points and
triangles from the mesh file.
"""
function EmbeddedDeltaSet2D(fn::String)
  if(splitext(fn)[2] == ".stl")
    # The .stl format references points by location so we apply the
    # force_unique flag
    EmbeddedDeltaSet2D(FileIO.load(fn); force_unique=true)
  else
    EmbeddedDeltaSet2D(FileIO.load(fn))
  end
end

end


text:
tests = [
    "types",
    "utils",
    "samplers",
    "categorical",
    "univariates",
    "fit",
    "multinomial",
    "dirichlet",
    "mvnormal",
    "mvtdist",
    "kolmogorov",
    "edgeworth",
    "matrix",
    "noncentralhypergeometric",
    "vonmisesfisher",
    "compoundvariate",
    "conjugates",
    "conjugates_normal",
    "conjugates_mvnormal",
    "mixture",
    "gradlogpdf"]

print_with_color(:blue, "Running tests:\n")

srand(345678)

for t in tests
    test_fn = "$t.jl"
    print_with_color(:green, "* $test_fn\n")
    include(test_fn)
end


text:
using OrdinaryDiffEq, Test

@testset "Differentiation Trait Tests" begin
  jac_called = Ref(false)
  tgrad_called = Ref(false)

  function Lotka(du,u,p,t)
    du[1] = u[1] - u[1] * u[2] # REPL[7], line 3:
    du[2] = -3 * u[2] + 1 * u[1] * u[2]
    nothing
  end

  function Lotka_jac(J,u,p,t)
    jac_called.x = true
    J[1,1] = 1.0 - u[2]
    J[1,2] = -u[1]
    J[2,1] = 1 * u[2]
    J[2,2] = -3 + u[1]
    nothing
  end

  function Lotka_tgrad(grad,u,p,t)
    tgrad_called.x = true
    grad[1] = 1 * 0
    grad[2] = 1 * 0
    nothing
  end

  Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,tgrad=Lotka_tgrad)
  prob = ODEProblem(Lotka_f,ones(2),(0.0,10.0))

  good_sol = solve(prob,Rosenbrock23())

  @test jac_called[]
  @test tgrad_called[]

  prob2 = ODEProblem(Lotka,ones(2),(0.0,10.0))

  sol = solve(prob2,Rosenbrock23(autodiff=true))
  @test ≈(good_sol[end],sol[end],rtol=1e-2)

  sol = solve(prob2,Rosenbrock23(autodiff=true,chunk_size=1))
  @test ≈(good_sol[end],sol[end],rtol=1e-2)

  sol = solve(prob2,Rosenbrock23(autodiff=false))
  @test ≈(good_sol[end],sol[end],rtol=1e-2)

  invW_called = Ref(false)

  function Lotka_invW(iW,u,p,γ,t)
    invW_called.x = true
    iW[1,1] = (1 - (1 * u[1] * u[2] * γ ^ 2) / (((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) * (1 - (1 - 1 * u[2]) * γ))) / (1 - (1 - 1 * u[2]) * γ)
    iW[1,2] = (-(1) * u[1] * γ) / (((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) * (1 - (1 - 1 * u[2]) * γ))
    iW[2,1] = (u[2] * γ) / (((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) * (1 - (1 - 1 * u[2]) * γ))
    iW[2,2] = ((1 - (-3 + u[1]) * γ) + (1 * u[1] * u[2] * γ ^ 2) / (1 - (1 - 1 * u[2]) * γ)) ^ -1
    nothing
  end

  Lotka_f2 = ODEFunction(Lotka,jac=Lotka_jac,tgrad=Lotka_tgrad,invW = Lotka_invW)
  prob3 = ODEProblem(Lotka_f2,ones(2),(0.0,10.0))
  inv_sol = solve(prob3,Rosenbrock23())

  @test invW_called[]
  @test ≈(good_sol[end],inv_sol[end],rtol=1e-2)
end


text:
module XML

using LabelMe
using LightXML

export load, annotation, source, image_size, object, polygon, element

function load(filename::AbstractString)
    e = root(parse_file(filename))
    return annotation(e)
end

function annotation(e::XMLElement)
    a = Annotation()
    a.filename = element_content(e, "filename")
    a.folder = element_content(e, "folder")
    a.source = source(find_element(e, "source"))
    a.image_size = image_size(find_element(e, "imagesize"))
    for object_element in get_elements_by_tagname(e, "object")
        push!(a.objects, object(object_element))
    end
    return a
end

function source(e::Void)
    return Source()
end

function source(e::XMLElement)
    s = Source()
    s.sourceImage = element_content(e, "sourceImage")
    s.sourceAnnotation = element_content(e, "sourceAnnotation")
    return s
end

function image_size(e::XMLElement)
    nrows = tryparse(Int64, element_content(e, "nrows"))
    ncols = tryparse(Int64, element_content(e, "ncols"))
    if isnull(nrows) || isnull(ncols)
        return nothing
    end
    return (get(nrows), get(ncols))
end

function object(e::XMLElement)
    o = Object()
    o.id = tryparse(Int64, element_content(e, "id"))
    o.name = element_content(e, "name")
    o.deleted = tryparse(UInt8, element_content(e, "deleted"))
    o.verified = tryparse(UInt8, element_content(e, "verified"))
    o.occluded = tryparse(UInt8, element_content(e, "occluded"))
    o.attributes = element_content(e, "attributes")
    o.parts_parent = tryparse(Int64, element_content(find_element(e, "parts"), "ispartof"))
    hasparts = element_content(find_element(e, "parts"), "hasparts")
    o.parts_children = [parse(Int64, p) for p in split(hasparts, ','; keep=false)]
    o.date = DateTime(element_content(e, "date"), "d-u-yyyy H:M:S")
    o.polygon = polygon(find_element(e, "polygon"))
    return o
end

function polygon(e::XMLElement)
    p = Polygon()
    p.username = element_content(e, "username")
    for pt_element in get_elements_by_tagname(e, "pt")
        push!(p.points, (
            parse(Float32, element_content(pt_element, "x")),
            parse(Float32, element_content(pt_element, "y"))
            ))
    end
    return p
end

function element_content(e::XMLElement, tag::AbstractString)
    first = find_element(e, tag)
    if first == nothing
        return ""
    end
    return content(first)
end

function element(a::Annotation)
    e = new_element("annotation")
    add_text(new_child(e, "filename"), a.filename)
    add_text(new_child(e, "folder"), a.folder)
    add_child(e, element(a.source))
    i = new_child(e, "imagesize")
    add_text(new_child(i, "nrows"), isnull(a.image_size) ? "" : string(get(a.image_size)[1]))
    add_text(new_child(i, "ncols"), isnull(a.image_size) ? "" : string(get(a.image_size)[2]))
    for o in a.objects
        add_child(e, element(o))
    end
    return e
end

function element(s::Source)
    e = new_element("source")
    add_text(new_child(e, "sourceImage"), s.sourceImage)
    add_text(new_child(e, "sourceAnnotation"), s.sourceAnnotation)
    return e
end

function element(o::Object)
    e = new_element("object")
    if !isnull(o.id)
        add_text(new_child(e, "id"), string(get(o.id)))
    end
    add_text(new_child(e, "name"), o.name)
    if !isnull(o.deleted)
        add_text(new_child(e, "deleted"), string(UInt8(get(o.deleted))))
    end
    if !isnull(o.verified)
        add_text(new_child(e, "verified"), string(UInt8(get(o.verified))))
    end
    if !isnull(o.occluded)
        add_text(new_child(e, "occluded"), string(UInt8(get(o.occluded))))
    end
    add_text(new_child(e, "attributes"), o.attributes)
    p = new_child(e, "parts")
    if !isnull(o.parts_parent)
        add_text(new_child(p, "ispartof"), string(get(o.parts_parent)))
    end
    add_text(new_child(p, "hasparts"), join(o.parts_children, ","))
    add_text(new_child(e, "date"), Dates.format(o.date, "d-u-yyyy H:M:S"))
    add_child(e, element(o.polygon))
    return e
end

function element(p::Polygon)
    e = new_element("polygon")
    add_text(new_child(e, "username"), p.username)
    for pt in p.points
        pt_e = new_child(e, "pt")
        add_text(new_child(pt_e, "x"), string(pt[1]))
        add_text(new_child(pt_e, "y"), string(pt[2]))
    end
    return e
end

end # module


text:
using Base.Test
using LibFoo

@test call_fooifier(2.2, 1.1) ≈ 2*2.2^2 - 1.1
@test call_libfoo(2.2, 1.1) ≈ 2*2.2^2 - 1.1


text:
# Compare two digits
print("enter digit: a = ")
a = read(stdin, Char)
print("\n")
#x = readline()
print("enter digit: b = ")
b = read(stdin, Char)
print("\n")
# x = readline()
# Two branches conditional
if a > b
   println("a > b")
else
   println("a <= b")
end

print("done.")

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
for N in [Float64, Rational{Int}, Float32]
    # ==============================
    # Testing box approximation
    # ==============================
    # Box approximation of a 2D square
    b = BallInf(N[1, 1], N(1//10))
    h = box_approximation(b)
    hexp = Hyperrectangle(N[1, 1], N[1//10, 1//10])
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Box approximation of a 2D unit ball in the 1-norm
    b = Ball1(N[1, -2], N(0.2))
    h = box_approximation(b)
    hexp = Hyperrectangle(N[1, -2], N[0.2, 0.2])
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Box approximation of a 3D unit ball in the 1-norm
    b = Ball1(N[1, 2, 0], N(1))
    h = box_approximation(b)
    hexp = Hyperrectangle(N[1, 2, 0], N[1, 1, 1])
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Box approximation of a Hyperrectangle
    hexp = Hyperrectangle(N[0, 0], N[1, 1])
    h = box_approximation(hexp)
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Box approximation of a VPolytope
    P = VPolytope([N[0, 0], N[2, 1], N[1, 2]])
    H = box_approximation(P)
    @test H == Hyperrectangle(N[1, 1], N[1, 1])

    # empty set
    E = EmptySet{N}(2)
    @test box_approximation(E) == E

    # ===================================================================
    # Testing box_approximation_symmetric (= symmetric interval hull)
    # ===================================================================
    # Box approximation of a 2D square
    b = BallInf(N[1, 1], N(1//10))
    h = box_approximation_symmetric(b)
    hexp = Hyperrectangle(N[0, 0], N[11//10, 11//10])
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Box approximation of a 2D unit ball in the 1-norm
    b = Ball1(N[1, -2], N(2//10))
    h = box_approximation_symmetric(b)
    hexp = Hyperrectangle(N[0, 0], N[12//10, 22//10])
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Box approximation of a 3D unit ball in the 1-norm
    b = Ball1(N[1, 2, 0], N(1//10))
    h = box_approximation_symmetric(b)
    hexp = Hyperrectangle(N[0, 0, 0], N[11//10, 21//10, 1//10])
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Box approximation of a 4D hyperrectangle
    b = Hyperrectangle(N[-15//10, -25//10, 24//10, -4//10],
                       N[1//10, 2//10, 3//10, 4//10])
    h = box_approximation_symmetric(b)
    hexp = Hyperrectangle(zeros(N, 4), N[16//10, 27//10, 27//10, 8//10])
    @test h.center ≈ hexp.center
    @test h.radius ≈ hexp.radius

    # Testing alias symmetric_interval_hull
    h = symmetric_interval_hull(b)

    # empty set
    E = EmptySet{N}(2)
    @test box_approximation_symmetric(E) == E
end

for N in [Float64]
    # empty intersection due to line search
    X = HalfSpace(N[-1], N(0)) ∩ HalfSpace(N[1], N(-1e-15))
    @test box_approximation(X) isa EmptySet{N}
end


text:
using ATools
using DataFrames
using Glob
using CSV
using Plots

"""
    readdf(dir)

Reads all csv files found in dir and returns a single DataFrame.
"""
function readdf(dir)
    files = glob("*.csv",dir)
    dfs =[DataFrame(CSV.File(file)) for file in files]
    evtdf=vcat(dfs...)
end

function writemdf(dir, file, df)
    path = string(dir,"/", file)
    CSV.write(path, df)
end

function pathfinder(rootdir, conf)
    mdf = string(conf, "-mdf")
    lor = string(conf, "-lor")
    phif = string("mdf-phistd-",conf,".csv")
    zf = string("mdf-zstd-",conf,".csv")

    mdfdir = joinpath(rootdir, mdf)
    phipath = joinpath(mdfdir, phif)
    zpath = joinpath(mdfdir, zf)

    lordir = joinpath(rootdir, lor)

    if isdir(lordir) == false
        mkdir(lordir)
    end

    philp = string("lor-phistd-",conf)
    zlp = string("lor-zstd-",conf)
    return (phipath= phipath, zpath = zpath, lordir=lordir, philp = philp, zlp = zlp)

end

function getlorpath(lordir, prefix, posfix="first-true")
    fullpath = string(prefix,"-",posfix,".h5")
    joinpath(lordir, fullpath)
end

"""
    selectinterval(df, column1, column2, xmin, xmax)

Select interval over 2 hemisphere columns (e.g, q1, q2) for df
"""
function selectinterval(df, column1, column2, xmin, xmax)
    e1 = select_by_column_value_interval(df, column1, xmin,xmax)
    select_by_column_value_interval(e1, column2, xmin,xmax)
end

function pout(p,filename)
    if filename != ""
        png(p,filename)
    end
end

function plot_and_save(p1::Plots.Plot, tit::String="", filename::String="", legend=true)
    if tit != ""
        p = plot(p1, layout= (1, 1), title=tit, legend=legend, fmt = :png, size = (1000, 400),
            left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    else
        p = plot(p1, layout= (1, 1), legend=legend, fmt = :png, size = (1000, 400),
            left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    end
    pout(p,filename)
    return p
end

function plot_and_save(p1::Plots.Plot,p2::Plots.Plot,
                       tit::String="", filename::String="", legend=true)
    if tit != ""
        p = plot(p1, p2, layout= (1, 2), title=tit, legend=legend, fmt = :png, size = (1000, 400),
            left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    else
        p = plot(p1, p2, layout= (1, 2), legend=legend, fmt = :png, size = (1000, 400),
            left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    end
    pout(p,filename)
    return p
end

function plot_and_save(p1::Plots.Plot,p2::Plots.Plot, p3::Plots.Plot, tit::String="",
    filename::String="", legend=true)
    if tit != ""
        p = plot(p1, p2, p3, layout= (1, 3), title=tit, legend=legend, fmt = :png, size = (1000, 400),
            left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    else
        p = plot(p1, p2, p3, layout= (1, 3), legend=legend, fmt = :png, size = (1000, 400),
            left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    end
    pout(p,filename)
    return p
end

# function plot_and_save(p1,p2,p3,p4, tit="", filename="")
#     if tit != ""
#         p = plot(p1, p2, p3, p4, layout= (2, 2), title=tit, legend=false, fmt = :png, size = (1000, 400),
#             left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
#     else
#         p = plot(p1, p2, p3, p4, layout= (2, 2), legend=false, fmt = :png, size = (1000, 400),
#             left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
#     end
#     pout(p,filename)
#     return p
# end

# function condplot(p1, tit, filename, save)
#     if save
#         plot_and_save(p1, tit, filename)
#     else
#         plot(p1, layout= (1, 1), legend=false, fmt = :png, size = (1000, 400),
#             left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
#     end
# end

function condplot(p1, p2, tit, filename, save)
    if save
        plot_and_save(p1,p2, tit, filename)
    else
        plot(p1, p2, layout= (1, 2), legend=false, fmt = :png, size = (1000, 400),
        left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    end
end

# function condplot(p1, p2, p3, tit, filename, save)
#     if save
#         plot_and_save(p1,p2, p3, tit, filename)
#     else
#         plot(p1, p2, p3, layout= (1, 3), legend=false, fmt = :png, size = (1000, 400),
#         left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
#     end
# end

function condplot(p1, p2, p3, p4, tit, filename, save)
    if save
        plot_and_save(p1,p2, p3, p4, tit, filename)
    else
        plot(p1, p2, p3, p4, layout= (2, 2), legend=false, fmt = :png, size = (1000, 400),
        left_margin=5Plots.mm, right_margin=1Plots.mm, bottom_margin=5Plots.mm)
    end
end

function q1vsq2(df, nbins=100, qmin=50.0, qmax=5000.0; tit="", filename="", save=false)
    h1,p1 = hist2d(df.q1, df.q2, nbins, "q1 (pes)","q2 (pes)", qmin, qmax, qmin, qmax)
    h2,p2 = hist1d(df.q1, "q1 (pes)", nbins, qmin, qmax);
    condplot(p1, p2, tit, filename, save)

end

function r1q1(df, nbins=100, qmin=50.0, qmax=5000.0, rmin=300.0, rmax=400.0; tit="", filename="", save=false)
    h1,p1 = hist2d(df.q1, df.r1, nbins, "q1 (pes) ","r1 (mm)",qmin, qmax, rmin,rmax)
    h2,p2 = hist1d(df.r1, "r1", nbins, rmin,rmax)

    condplot(p1, p2, tit, filename, save)
end

function zstd(df, nbins=100, zmin=0.0, zmax=40.0,rmin=300.0, rmax=400.0; tit="", filename="", save=false)
    h1,p1 = hist2d(df.zstd1, df.r1, nbins, "σz (mm) ","1 (mm)",zmin, zmax, rmin, rmax)
    h2,p2 = hist1d(df.zstd1, "σz (mm)", nbins, zmin, zmax)
    condplot(p1, p2, tit, filename, save)

end

function phistd(df, nbins=100,phimin=0.0, phimax=0.1,rmin=300.0, rmax=400.0; tit="", filename="", save=false)
    h1,p1 = hist2d(df.phistd1, df.r1, nbins, "σϕ (mm) ","r1 (mm)",phimin, phimax,rmin, rmax)
    h2,p2 = hist1d(df.phistd1, "σϕ (mm)", 100, phimin, phimax)
    condplot(p1, p2, tit, filename, save)
end

function plotreso(r1t, r1x, tx1, ty1, xmin, xmax, nbins=150; tit="", filename="", save=false)
    h1,p1 = hist2d(r1x, r1t, nbins, tx1, ty1)
    h2,p2 = hist1d(r1t - r1x, tit, nbins, xmin, xmax)
    condplot(p1, p2, tit, filename, save)
end

max_stars_repo_name: jjgomezcadenas/mytest

max_issues_repo_name: jjgomezcadenas/mytest

max_forks_repo_name: jjgomezcadenas/mytest


text:
# hide!(vis::Visualizer) = settransform!(vis, LinearMap(0I))

struct CentroidalTrajectoryVisualizer
    vis::Visualizer
    com_visualizer::Visualizer
    region_data::Vector{ContactRegion{Float64}}
    force_visualizers::Vector{ArrowVisualizer{Visualizer}}
    cone_visualizers::Vector{Visualizer}
    contact_position_visualizers::Vector{Visualizer}
    gravity_mag::Float64
end

function CentroidalTrajectoryVisualizer(vis::Visualizer,
        region_data::Vector{ContactRegion{Float64}}, gravity_mag, num_contacts::Number)

    # CoM
    com_visualizer = vis[:com]
    img = PngImage(joinpath(@__DIR__, "..", "assets", "checkerboard.png"))
    com_material = MeshBasicMaterial(map=Texture(image=img))#, wrap=(1, 1), repeat=(1, 1)))
    setobject!(com_visualizer, HyperSphere(Point(0., 0, 0), 0.05), com_material)
    setvisible!(com_visualizer, false)

    # Forces / contact force cones
    force_visualizers = [ArrowVisualizer(vis["f_$i"]) for i = 1 : num_contacts]
    cone_visualizers = [vis["cone_$i"] for i = 1 : num_contacts]
    force_orange = RGB(243 / 255, 118 / 255, 32 / 255)
    force_orange_transparent = RGBA(force_orange, 0.2)
    for (force_visualizer, cone_visualizer) in zip(force_visualizers, cone_visualizers)
        setobject!(force_visualizer, MeshLambertMaterial(; color=force_orange))
        cone_height = 0.2
        cone = Cone(Point(0., 0., cone_height), Point(0., 0., 0.), cone_height)
        setobject!(cone_visualizer, cone, MeshLambertMaterial(color=force_orange_transparent))
        setvisible!(cone_visualizer, false)
    end

    # Contact positions
    contact_position_visualizers = [vis["p_$i"] for i = 1 : num_contacts]
    for contact_position_visualizer in contact_position_visualizers
        setobject!(contact_position_visualizer, HyperSphere(Point(0., 0, 0), 0.02), MeshLambertMaterial(color=RGB(0.1, 0.1, 0.1)))
        setvisible!(contact_position_visualizer, false)
    end

    # Regions
    for (i, region) in enumerate(region_data)
        extruded = extrude(hrep(region.A, region.b), 1e-4; zmax=1e-4)
        world_to_body = inv(region.transform)
        geometry = GLNormalMesh(polyhedron(hrep(extruded.A * world_to_body.linear, extruded.b - extruded.A * world_to_body.translation)))
        setobject!(vis["region_$i"], geometry, MeshLambertMaterial(color=RGB(0.5, 0.5, 0.5)))
    end

    CentroidalTrajectoryVisualizer(vis, com_visualizer, region_data,
        force_visualizers, cone_visualizers, contact_position_visualizers, gravity_mag)
end

function set_com_trajectory!(vis::CentroidalTrajectoryVisualizer, result::CentroidalTrajectoryResult)
    ts = range(first(result.break_times), last(result.break_times); length=100)
    geometry = PointCloud([Point(result.center_of_mass(t)) for t in ts])
    setobject!(vis.vis[:com_trajectory], LineSegments(geometry, LineBasicMaterial(linewidth=3)))
    vis
end

function set_state!(vis::CentroidalTrajectoryVisualizer, result::CentroidalTrajectoryResult, t::Number; show_cones::Bool=false)
    # CoM
    c = result.center_of_mass(t)
    setvisible!(vis.com_visualizer, true)
    settransform!(vis.com_visualizer, Translation(c))

    # Forces
    for (force_vis, cone_vis, cop, force, indicator) in zip(vis.force_visualizers, vis.cone_visualizers, result.centers_of_pressure, result.contact_forces, result.contact_indicators)
        r = cop(t)
        f = force(t)
        settransform!(force_vis, Point(r), Vec(f / vis.gravity_mag))
        zs = indicator(t)
        if sum(zs) > 0.5 && show_cones
            setvisible!(cone_vis, true)
            region_index = argmax(zs)
            region = vis.region_data[region_index]
            n = region.transform.linear[:, 3]
            μ = region.μ
            R = rotation_between(SVector(0, 0, 1), n)
            settransform!(cone_vis, Translation(r) ∘ LinearMap(R) ∘ LinearMap(Diagonal(SVector(1, 1, 1 / μ))))
        else
            setvisible!(cone_vis, false)
            settransform!(cone_vis, LinearMap(0I))
        end
    end

    # Contact positions
    for (position_vis, position, indicator) in zip(vis.contact_position_visualizers, result.contact_positions, result.contact_indicators)
        p = position(t)
        zs = indicator(t)
        if sum(zs) > 0.5
            setvisible!(position_vis, true)
            settransform!(position_vis, Translation(p))
        else
            setvisible!(position_vis, false)
            settransform!(position_vis, LinearMap(0I))
        end
    end
    vis
end

function MeshCat.Animation(
        cvis::CentroidalTrajectoryVisualizer,
        result::CentroidalTrajectoryResult;
        fps::Number=30, kwargs...)
    t0 = first(result.break_times)
    tf = last(result.break_times)
    animation = Animation(fps)
    num_frames = floor(Int, (tf - t0) * fps)
    for frame in 0 : num_frames
        t = t0 + frame / fps
        atframe(animation, frame) do
            set_state!(cvis, result, t; kwargs...)
        end
    end
    return animation
end


text:
# Autogenerated wrapper script for GoogleTest_jll for x86_64-apple-darwin
export libgtest, libgtest_main

JLLWrappers.@generate_wrapper_header("GoogleTest")
JLLWrappers.@declare_library_product(libgtest, "@rpath/libgtest.dylib")
JLLWrappers.@declare_library_product(libgtest_main, "@rpath/libgtest_main.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgtest,
        "lib/libgtest.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgtest_main,
        "lib/libgtest_main.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("kruskal.jl")

"""
    count_spanning_trees(g::AGraph)

Returns the number of spanning trees of `g`, computed through
[Kirchhoff's theorem](https://en.wikipedia.org/wiki/Kirchhoff%27s_theorem).
The return type is a float, since the number can be very large.
"""
function count_spanning_trees(g::AGraph)
    @assert is_connected(g) "The graph has to be connected"
    nv(g) <= 1 && return 1.
    Δ = laplacian_matrix(g)[2:nv(g),2:nv(g)]
    return round(det(Δ))
end


text:
using Metal

@show devices()
dev = MtlDevice(1)

bufferSize = 128
bufferA = MtlArray{Float32,1}(undef, tuple(bufferSize), storage=Shared)
bufferB = MtlArray{Float32,1}(undef, tuple(bufferSize), storage=Shared)
bufferC = MtlArray{Float32,1}(undef, tuple(bufferSize), storage=Shared)

vecA = unsafe_wrap(Vector{Float32}, bufferA.buffer, tuple(bufferSize))
vecB = unsafe_wrap(Vector{Float32}, bufferB.buffer, tuple(bufferSize))
vecC = unsafe_wrap(Vector{Float32}, bufferC.buffer, tuple(bufferSize))

using Random
rand!.([vecA, vecB])
vecC .= 0.0

## Setup
src = read(dirname(pathof(Metal))*"/Metal/kernels/vadd.metal", String)
opts = MtlCompileOptions()
lib = MtlLibrary(dev, src, opts)

fun = MtlFunction(lib, "add_vectors")
pip_addfun = MtlComputePipelineState(dev, fun)
queue = global_queue(dev) #MtlCommandQueue(dev)

##
argbufenc = MtlArgumentEncoder(fun, 1)
argbuf1 = alloc(Cchar, dev, sizeof(argbufenc), storage=Shared)
MTL.assign_argument_buffer!(argbufenc, argbuf1, 1)
MTL.set_field!(argbufenc, 128, 1)
set_buffer!(argbufenc, bufferA.buffer, 0, 2)
content1 = (128, Base.unsafe_load(Base.convert(Ptr{Int}, MTL.content(argbuf1)),2))

argbufenc = MtlArgumentEncoder(fun, 2)
argbuf2 = alloc(Cchar, dev, sizeof(argbufenc), storage=Shared)
MTL.assign_argument_buffer!(argbufenc, argbuf2, 1)
MTL.set_field!(argbufenc, 128, 1)
set_buffer!(argbufenc, bufferB.buffer, 0, 2)
content2 = (128, Base.unsafe_load(Base.convert(Ptr{Int}, MTL.content(argbuf2)),2))

argbufenc = MtlArgumentEncoder(fun, 3)
argbuf3 = alloc(Cchar, dev, sizeof(argbufenc), storage=Shared)
MTL.assign_argument_buffer!(argbufenc, argbuf3, 1)
MTL.set_field!(argbufenc, 128, 1)
set_buffer!(argbufenc, bufferC.buffer, 0, 2)
content3 = (128, Base.unsafe_load(Base.convert(Ptr{Int}, MTL.content(argbuf3)),2))

rc1 = Base.RefValue(content1)
rc2 = Base.RefValue(content2)
rc3 = Base.RefValue(content3)
pt1 = Base.unsafe_convert(Ptr{typeof(content1)}, rc1)
pt2 = Base.unsafe_convert(Ptr{typeof(content1)}, rc2)
pt3 = Base.unsafe_convert(Ptr{typeof(content1)}, rc3)

cmd = MTL.commit!(queue) do cmdbuf
    MtlComputeCommandEncoder(cmdbuf) do enc
        MTL.use!(enc, bufferA.buffer, MTL.ReadWriteUsage)
        MTL.use!(enc, bufferB.buffer, MTL.ReadWriteUsage)
        MTL.use!(enc, bufferC.buffer, MTL.ReadWriteUsage)
        #MTL.set_buffers!(enc,
        #                        [argbuf1, argbuf2, argbuf3],
        #                        [0,0,0], 1:3)
        MTL.set_bytes!(enc, pt1, 16, 1)
        MTL.set_bytes!(enc, pt2, 16, 2)
        MTL.set_bytes!(enc, pt3, 16, 3)
        gridSize = MtSize(length(vecA), 1, 1)
        MTL.set_function!(enc, pip_addfun)
        threadGroupSize = min(length(vecA), pip_addfun.maxTotalThreadsPerThreadgroup)
        threadGroupSize = MTL.MtSize(threadGroupSize, 1, 1)
        MTL.append_current_function!(enc, gridSize, threadGroupSize)
    end
end

####

# Execute
wait(cmd)

@show vecC


text:
function spline!(n, x, f, cf)
# !INPUT/OUTPUT PARAMETERS:
#   n  : number of points (in,integer)
#   x  : abscissa array (in,real(n))
#   f  : input data array (in,real(n))
#   cf : cubic spline coefficients (out,real(3,n))
# !DESCRIPTION:
#   Calculates the coefficients of a cubic spline fitted to input data. In other
#   words, given a set of data points $f_i$ defined at $x_i$, where
#   $i=1\ldots n$, the coefficients $c_j^i$ are determined such that
#   $$ y_i(x)=f_i+c_1^i(x-x_i)+c_2^i(x-x_i)^2+c_3^i(x-x_i)^3, $$
#   is the interpolating function for $x\in[x_i,x_{i+1})$. The coefficients are
#   determined piecewise by fitting a cubic polynomial to adjacent points.
#
# integer, intent(in) :: n
# real(8), intent(in) :: x(n),f(n)
# real(8), intent(out) :: cf(3,n)

    if n <= 0
        error("Error(spline) n is invalid")
    end

    if n == 1
        for i in 1:3
            cf[i,1] = 0.0
        end
        return
    end

    if n == 2
        cf[1,1] = ( f[2] - f[1])/( x[2] - x[1] )
        cf[2,1] = 0.0
        cf[3,1] = 0.0
        cf[1,2] = cf[1,1]
        cf[2,2] = 0.0
        ct[3,2] = 0.0
        return
    end

    if n == 3
        x0 = x[1]
        x1 = x[2] - x0
        x2 = x[3] - x0
        y0 = f[1]
        y1 = f[2] - y0
        y2 = f[3] - y0
        t0 = 1.0/(x1*x2*(x2-x1))
        t3 = x1*y2
        t4 = x2*y1
        t1 = t0*(x2*t4 - x1*t3)
        t2 = t0*(t3 - t4)
        cf[1,1] = t1
        cf[2,1] = t2
        cf[3,1] = 0.0
        t3 = 2.0*t2
        cf[1,2] = t1 + t3*x1
        cf[2,2] = t2
        cf[3,2] = 0.0
        cf[1,3] = t1 + t3*x2
        cf[2,3] = t2
        cf[3,3] = 0.0
        return
    end

    x0 = x[1]
    x1 = x[2] - x0
    x2 = x[3] - x0
    x3 = x[4] - x0
    t4 = x1 - x2
    t5 = x1 - x3
    t6 = x2 - x3
    y0 = f[1]
    y1 = f[2] - y0
    y2 = f[3] - y0
    y3 = f[4] - y0
    t1 = x1*x2*y3
    t2 = x2*x3*y1
    t3 = x1*x3
    t0 = 1.0/(x2*t3*t4*t5*t6)
    t3 = t3*y2
    t7 = t0*(t1*t4 + t2*t6 - t3*t5)
    t4 = x1^2
    t5 = x2^2
    t6 = x3^2
    y1 = t3*t6 - t1*t5
    y3 = t2*t5 - t3*t4
    y2 = t1*t4 - t2*t6
    t1 = t0*(x1*y1 + x2*y2 + x3*y3)
    t2 = -t0*(y1 + y2 + y3)

    cf[1,1] = t1
    cf[2,1] = t2
    cf[3,1] = t7
    cf[1,2] = t1 + 2.0*t2*x1 + 3.0*t7*t4
    cf[2,2] = t2 + 3.0*t7*x1
    cf[3,2] = t7

    if n == 4
        cf[1,3] = t1 + 2.0*t2*x2 + 3.0*t7*t5
        cf[2,3] = t2 + 3.0*t7*x2
        cf[3,3] = t7
        cf[1,4] = t1 + 2.0*t2*x3 + 3.0*t7*t6
        cf[2,4] = t2 + 3.0*t7*x3
        cf[3,4] = t7
        return
    end

    for i in 3:n-2
        x0 = x[i]
        x1 = x[i-1] - x0
        x2 = x[i+1] - x0
        x3 = x[i+2] - x0
        t4 = x1 - x2
        t5 = x1 - x3
        t6 = x2 - x3
        y0 = f[i]
        y1 = f[i-1] - y0
        y2 = f[i+1] - y0
        y3 = f[i+2] - y0
        t1 = x1*x2*y3
        t2 = x2*x3*y1
        t3 = x1*x3
        t0 = 1.0/(x2*t3*t4*t5*t6)
        t3 = t3*y2
        t7 = t0*(t1*t4+t2*t6-t3*t5)
        t4 = x1^2
        t5 = x2^2
        t6 = x3^2
        y1 = t3*t6 - t1*t5
        y2 = t1*t4 - t2*t6
        y3 = t2*t5 - t3*t4
        t1 = t0*(x1*y1 + x2*y2 + x3*y3)
        t2 = -t0*(y1 + y2 + y3)
        cf[1,i] = t1
        cf[2,i] = t2
        cf[3,i] = t7
    end

    cf[1,n-1] = t1 + 2.0*t2*x2 + 3.0*t7*t5
    cf[2,n-1] = t2 + 3.0*t7*x2
    cf[3,n-1] = t7
    cf[1,n] = t1 + 2.0*t2*x3 + 3.0*t7*t6
    cf[2,n] = t2 + 3.0*t7*x3
    cf[3,n] = t7
    return
end


text:
module MooreDDQN

using Flux
using Flux.NNlib
using Juno
using DataStructures: CircularBuffer
using BSON: @save, @load
import BSON
using ProgressMeter
using LearnBase: DiscreteSet
using Reinforce
import Reinforce: action, reset!, finished
import StatsBase

include("sequences.jl")
export ConstantSequence,ExponentialSequence,HeavisideSequence,ArithmeticSequence,LinearSequence

include("memory.jl")
export ReplayMemoryBuffer, PriorityReplayMemoryBuffer

include("policies.jl")
export QPolicy, ϵGreedyPolicy, DeepQPolicy, Double_DeepQPolicy, DuelingDouble_DeepQPolicy, reset!, save_policy, load_policy

include("utils.jl")

include("learn.jl")
export learn!

end # module


text:
abstract type AbstractKvec{T<:AbstractFloat} <: AbstractVector{T} end

Base.IndexStyle(::Type{<:AbstractKvec}) = IndexLinear()
@inline Base.length(a::AbstractKvec) = a.n
@inline Base.size(a::AbstractKvec) = (a.n,)

struct RKvec{T<:AbstractFloat} <: AbstractKvec{T}
    n::Int
    l::T
    @inline function RKvec{T}(nx::Integer,lx::Real) where {T}
        n = div(nx,2) + 1
        return new{T}(n,T(lx))
    end
end

struct SRKvec{T<:AbstractFloat,n,l} <: AbstractKvec{T}
    @inline function SRKvec(nx::Integer,lx::T) where {T}
        n = div(nx,2) + 1
        return new{T,n,lx}()
    end
end

@inline function Base.getproperty(a::SRKvec{T,n,l},s::Symbol) where {T,n,l}
    if s === :n
        return n
    elseif s === :l
        return l
    else
        return getfield(a,s)
    end
end

@inline function Base.getindex(a::Union{<:SRKvec,<:RKvec},i::Integer)
    l = a.l
    n = a.n
    @boundscheck checkbounds(a,i)
    return (i-1)/l
end

struct Kvec{T<:AbstractFloat} <: AbstractKvec{T}
    n::Int
    l::T
    p::Int
    @inline function Kvec{T}(n::Integer,l::Real) where {T}
        p = (n+1)÷2
        return new{T}(n,T(l),p)
    end
end

struct SKvec{T<:AbstractFloat,n,l,p} <: AbstractKvec{T}
    @inline function SKvec(n::Integer,l::T) where {T}
        p = (n+1)÷2
        return new{T,n,l,p}()
    end
end

@inline function Base.getproperty(a::SKvec{T,n,l,p},s::Symbol) where {T,n,l,p}
    if s === :n
        return n
    elseif s === :l
        return l
    elseif s === :p
        return p
    else
        return getfield(a,s)
    end
end

@inline function Base.getindex(a::Union{<:SKvec,<:Kvec},i::Integer)
    l = a.l
    n = a.n
    p = a.p
    @boundscheck checkbounds(a,i)
    no2 = (i <= p)
    return ifelse(no2, (i-1)/l, (i-n-1)/l)
end


text:
using Dates
using GitCommand
using GitHub
using JSON
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

const timestamp_regex = r"integration\/(\d\d\d\d-\d\d-\d\d-\d\d-\d\d-\d\d-\d\d\d)\/"

function wait_pr_compute_mergeability(repo::GitHub.Repo, pr::GitHub.PullRequest; auth::GitHub.Authorization)
    while !(pr.mergeable isa Bool)
        sleep(5)
        pr = GitHub.pull_request(repo, pr.number; auth = auth)
    end
    return pr
end

function close_all_pull_requests(repo::GitHub.Repo;
                                 auth::GitHub.Authorization,
                                 state::String)
    all_pull_requests = AutoMerge.get_all_pull_requests(repo,
                                                        state;
                                                        auth = auth)
    for pr in all_pull_requests
        try
            GitHub.close_pull_request(repo, pr; auth = auth)
        catch
        end
    end
    return nothing
end

function templates(parts...)
    this_filename = @__FILE__
    test_directory = dirname(this_filename)
    templates_directory = joinpath(test_directory, "templates")
    result = joinpath(templates_directory, parts...)
    return result
end

function get_random_number_from_system()
    result = parse(Int, strip(read(pipeline(`cat /dev/random`, `od -vAn -N4 -D`), String)))::Int
    return result
end

function my_two_datetimes_now()
    _now_zdt = now_localzone()
    _now_utc = astimezone(_now_zdt, tz"UTC")
    _now_et = astimezone(_now_zdt, tz"America/New_York")
    return _now_utc, _now_et
end

function my_two_times_now()
    _now_utc, _now_et = my_two_datetimes_now()
    utc_string = @sprintf "%02d:%02d UTC" Hour(_now_utc).value Minute(_now_utc).value
    et_string = @sprintf "%02d:%02d %s" Hour(_now_et).value Minute(_now_et).value _now_et.zone.name
    result = "$(et_string) ($(utc_string))"
end

@inline now_localzone() = now(localzone())

function utc_to_string(zdt::ZonedDateTime)
    zdt_as_utc = astimezone(zdt, tz"UTC")
    year = Year(zdt_as_utc.utc_datetime).value
    month = Month(zdt_as_utc.utc_datetime).value
    day = Day(zdt_as_utc.utc_datetime).value
    hour = Hour(zdt_as_utc.utc_datetime).value
    minute = Minute(zdt_as_utc.utc_datetime).value
    second = Second(zdt_as_utc.utc_datetime).value
    millisecond = Millisecond(zdt_as_utc.utc_datetime).value
    result = @sprintf "%04d-%02d-%02d-%02d-%02d-%02d-%03d" year month day hour minute second millisecond
    return result
end

function string_to_utc(s::AbstractString)
    df = DateFormat("yyyy-mm-dd-HH-MM-SS-sss")
    dt = DateTime(s, df)
    zdt = ZonedDateTime(dt, tz"UTC")
    return zdt
end

function list_all_origin_branches(git_repo_dir; GIT)
    result = Vector{String}(undef, 0)
    original_working_directory = pwd()
    cd(git_repo_dir)
    a = try
        read(`$(GIT) branch -a`, String)
    catch
        ""
    end
    b = split(strip(a), '\n')
    b_length = length(b)
    c = Vector{String}(undef, b_length)
    for i = 1:b_length
        c[i] = strip(strip(strip(b[i]), '*'))
        c[i] = first(split(c[i], "->"))
        c[i] = strip(c[i])
    end
    my_regex = r"^remotes\/origin\/(.*)$"
    for i = 1:b_length
        if occursin(my_regex, c[i])
            m = match(my_regex, c[i])
            if m[1] != "HEAD"
                push!(result, m[1])
            end
        end
    end
    cd(original_working_directory)
    return result
end

function delete_stale_branches(AUTOMERGE_INTEGRATION_TEST_REPO; GIT)
    with_cloned_repo(AUTOMERGE_INTEGRATION_TEST_REPO; GIT = GIT) do git_repo_dir
        cd(git_repo_dir)
        all_origin_branches = list_all_origin_branches(git_repo_dir; GIT=GIT)::Vector{String}
        for b in all_origin_branches
            if occursin(timestamp_regex, b)
                try
                    run(`$(GIT) push origin --delete $(b)`)
                catch
                end
            end
        end
    end
    return nothing
end

function empty_git_repo(git_repo_dir::AbstractString)
    original_working_directory = pwd()
    cd(git_repo_dir)
    for x in readdir(git_repo_dir)
        if x != ".git"
            path = joinpath(git_repo_dir, x)
            rm(path; force = true, recursive = true)
        end
    end
    cd(original_working_directory)
    return nothing
end

function with_temp_dir(f)
    original_working_directory = pwd()
    tmp_dir = mktempdir()
    atexit(() -> rm(tmp_dir; force = true, recursive = true))
    cd(tmp_dir)
    result = f(tmp_dir)
    cd(original_working_directory)
    rm(tmp_dir; force = true, recursive = true)
    return result
end

function with_cloned_repo(f, repo_url; GIT)
    original_working_directory = pwd()
    result = with_temp_dir() do dir
        git_repo_dir = joinpath(dir, "REPO")
        cd(dir)
        try
            run(`$(GIT) clone $(repo_url) REPO`)
        catch
        end
        cd(git_repo_dir)
        return f(git_repo_dir)
    end
    cd(original_working_directory)
    return result
end

function get_git_current_head(dir)
    original_working_directory = pwd()
    cd(dir)
    result = convert(String, strip(read(`git rev-parse HEAD`, String)))::String
    cd(original_working_directory)
    return result
end

function with_pr_merge_commit(f::Function,
                              pr::GitHub.PullRequest,
                              repo_url::AbstractString;
                              GIT)
    original_working_directory = pwd()
    result = with_cloned_repo(repo_url; GIT = GIT) do git_repo_dir
        cd(git_repo_dir)
        number = pr.number
        run(`$(GIT) fetch origin +refs/pull/$(number)/merge`)
        run(`$(GIT) checkout -qf FETCH_HEAD`)
        head = get_git_current_head(git_repo_dir)
        merge_commit_sha = pr.merge_commit_sha
        @test strip(head) == strip(merge_commit_sha)
        return f(git_repo_dir)
    end
    cd(original_working_directory)
    return result
end

function _generate_branch_name(name::AbstractString)
    sleep(0.1)
    _now = now_localzone()
    _now_utc_string = utc_to_string(_now)
    b = "integration/$(_now_utc_string)/$(rand(UInt32))/$(name)"
    sleep(0.1)
    return b
end

function generate_branch(name::AbstractString,
                         path_to_content::AbstractString,
                         parent_branch::AbstractString = "master";
                         GIT,
                         repo_url)
    original_working_directory = pwd()
    b = _generate_branch_name(name)
    with_cloned_repo(repo_url; GIT = GIT) do git_repo_dir
        cd(git_repo_dir)
        run(`$(GIT) checkout $(parent_branch)`)
        run(`$(GIT) branch $(b)`)
        run(`$(GIT) checkout $(b)`)
        empty_git_repo(git_repo_dir)
        for x in readdir(path_to_content)
            src = joinpath(path_to_content, x)
            dst = joinpath(git_repo_dir, x)
            rm(dst; force = true, recursive = true)
            cp(src, dst; force = true)
        end
        cd(git_repo_dir)
        try
            run(`$(GIT) add -A`)
        catch
        end
        try
            run(`$(GIT) commit -m "Automatic commit - AutoMerge integration tests"`)
        catch
        end
        try
            run(`$(GIT) push origin $(b)`)
        catch
        end
        cd(original_working_directory)
        rm(git_repo_dir; force = true, recursive = true)
    end
    return b
end

function generate_master_branch(path_to_content::AbstractString,
                                parent_branch::AbstractString = "master";
                                GIT,
                                repo_url)
    name = "master"
    b = generate_branch(name, path_to_content, parent_branch; GIT = GIT, repo_url = repo_url)
    return b
end

function generate_feature_branch(path_to_content::AbstractString,
                                 parent_branch::AbstractString;
                                 GIT,
                                 repo_url)
    name = "feature"
    b = generate_branch(name,
                        path_to_content,
                        parent_branch;
                        GIT = GIT,
                        repo_url = repo_url)
    return b
end

function with_master_branch(f::Function,
                            path_to_content::AbstractString,
                            parent_branch::AbstractString;
                            GIT,
                            repo_url)
    b = generate_master_branch(path_to_content,
                               parent_branch;
                               GIT = GIT,
                               repo_url = repo_url)
    result = f(b)
    return result
end

function with_feature_branch(f::Function,
                             path_to_content::AbstractString,
                             parent_branch::AbstractString;
                             GIT,
                             repo_url)
    b = generate_feature_branch(path_to_content,
                                parent_branch;
                                GIT = GIT,
                                repo_url = repo_url)
    result = f(b)
    return result
end


text:
"""

alert(message="Done!")

Display a cross-platform notification.

On MacOS, displays a notification window. On linux, tries to use notify-send,
zenity, kdialog or xmessage, in that order. On Windows or WSL2, uses a toast
notification.

Other platforms are not yet supported.

You can customize this function using [`set_alert_backend!`](@ref).

"""
alert(message="Done!") = alert_backend[](message)

"""

set_alert_backend!(fn)

Defines a custom backend for how `alert` sends messages to the user.
The argument should be a function of one argument (a string) which displays
the message to the user via some native UX api call.

If you wish to revert to the default backend, call this method with no arguments.

"""
set_alert_backend!(fn::Function) = alert_backend[] = fn
set_alert_backend!() = alert_backend[] = default_backend[]

const default_backend = Ref{Function}()
const alert_backend = Ref{Function}()

function init_alert_backends()
# select an appropriate implementation for sending alert notifications
default_backend[] = @static if Sys.isapple()
    message -> run(`osascript -e 'display notification "'$message'" with title "Julia"'`)
elseif Sys.islinux()
    if iswsl()
        win_toast
    elseif !isnothing(Sys.which("notify-send"))
        message -> run(`notify-send $message`)
    elseif !isnothing(Sys.which("zenity"))
        message -> run(pipeline(`echo $message`,`zenity --notification --listen`))
    elseif !isnothing(Sys.which("kdialog"))
        message -> run(`kdialog --title "Julia" --passivepopup $message 10`,wait=false)
    elseif !isnothing(Sys.which("xmessage"))
        message -> run(`xmessage $message`)
    else
        @warn("The `alert` method has no available local backend (call `alert()` for details).")
        message -> error("No viable messaging program available, install"*
        " 'notify-send', 'zenity', 'kdialog' or 'xmessage'.")
    end
elseif Sys.iswindows()
    win_toast
else
    @error "Unsupported operating system, no messages sent by `alert` will be "*
        "displayed. Consider using `AlertPushover`."
    message -> error("No alert backend installed! "*
        "Reload julia and read your error messages.")
end

alert_backend[] = default_backend[]
end

# Windows specific methods: `win_toast` and `iswsl`
# -----------------------------------------------------------------

@static if Sys.iswindows() || Sys.islinux()

# determine if the current system is Linux under WSL
function iswsl()
    try
        v = read("/proc/version", String)
        return occursin(r"microsoft",v)
    catch
        return false
    end
end

# display windows notification via powershell script
function win_toast(content)
    # format script input
    io = IOBuffer()
    show(io, content)
    disp_content = String(take!(io))

    # create a powershell script
    posh_script = """
    \$ErrorActionPreference = "Stop"

    \$notificationTitle = $disp_content

    [Windows.UI.Notifications.ToastNotificationManager, Windows.UI.Notifications, ContentType = WindowsRuntime] > \$null
    \$template = [Windows.UI.Notifications.ToastNotificationManager]::GetTemplateContent([Windows.UI.Notifications.ToastTemplateType]::ToastText01)

    #Convert to .NET type for XML manipuration
    \$toastXml = [xml] \$template.GetXml()
    \$toastXml.GetElementsByTagName("text").AppendChild(\$toastXml.CreateTextNode(\$notificationTitle)) > \$null

    #Convert back to WinRT type
    \$xml = New-Object Windows.Data.Xml.Dom.XmlDocument
    \$xml.LoadXml(\$toastXml.OuterXml)

    \$toast = [Windows.UI.Notifications.ToastNotification]::new(\$xml)
    \$toast.Tag = "PowerShell"
    \$toast.Group = "PowerShell"
    \$toast.ExpirationTime = [DateTimeOffset]::Now.AddMinutes(5)
    #\$toast.SuppressPopup = \$true

    \$notifier = [Windows.UI.Notifications.ToastNotificationManager]::CreateToastNotifier("PowerShell")
    \$notifier.Show(\$toast);
    """

    # encode as Base64 of Unicode-16
    io = IOBuffer()
    encode = Base64EncodePipe(io)
    write(encode, transcode(UInt16, posh_script))
    close(encode)
    str = String(take!(io))

    # send to powershell (works on windows and WSL2)
    run(`powershell.exe -Enc $str`)
end

end # if Sys.iswindows()