text:
module TypeInspections

import Crayons

_apply!(f, ::T, match_list, j, parent_type) where {T} = nothing # sometimes we need this...

function _apply!(f, ::Type{T}, match_list, j = nothing, parent_type = nothing) where {T}
    if f(T)
        push!(match_list, (T, j, parent_type))
    end
    for (i, p) in enumerate(T.parameters)
        _apply!(f, p, match_list, i, T)
    end
end

"""
    apply(::T) where {T <: Any}

Recursively traverse type `T` and apply
`f` to the types (and type parameters).

Returns a list of matches where `f(T)` is true.
"""
apply(f, ::T) where {T} = apply(f, T)
function apply(f, ::Type{T}) where {T}
    match_list = []
    _apply!(f, T, match_list)
    return match_list
end

function isconcretetype_params(::Type{T}) where {T}
    params = T.parameters
    isempty(params) && return true
    return all(p-> isconcretetype(p), params)
end

ith_name(i::Int) =
    if i == 1;     return "$(i)st"
    elseif i == 2; return "$(i)nd"
    elseif i == 3; return "$(i)rd"
    elseif 4 <= i <= 20; return "$(i)th"
    else; error("This is probably too many type parameters...")
    end

"""
    warntype_param(::T) where {T <: Any}

Print the type, warn (and highlight)
and non-concrete type parameters.
"""
warntype_param(::T) where {T} = warntype_param(T)
function warntype_param(::Type{T}) where {T}
    match_list = apply(x->!isconcretetype_params(x), T)
    if isempty(match_list)
        println(T)
    else
        sT_colored = string(T)
        for (ml, j, pt) in match_list
            sml = string(ml)
            ml_red = Crayons.Box.RED_FG(sml)
            sT_colored = replace(sT_colored, sml => ml_red)
        end
        println(sT_colored)
        s = "Type instance is not recursively concrete. "
        s *= "Non-concrete types are:\n"
        for (ml, j, pt) in match_list
            s *= "$(nameof(pt))'s ($(ith_name(j))) type parameter has a non-concrete eltype: "
            s *= "$(Crayons.Box.RED_FG(string(ml))) \n"
        end
        @warn s
    end
end

"""
    whose_type_parameter(::Type{T}, ::Type{FindType})

Finds where a type, `FindType` exists in the type-parameter
space of a complicated (and potentially nested) type `T`.

# Example
```julia
struct ComplicatedAndMaybeNestedFoo{A,B,C} end
struct Bar end
struct Baz end
struct Bingo end
f = ComplicatedAndMaybeNestedFoo{Bar, Baz, Bingo}()
whose_type_parameter(typeof(f), Bingo) # returns [(ComplicatedAndMaybeNestedFoo, 3)]
```
"""
function whose_type_parameter(::Type{T}, ::Type{FindType}) where {T, FindType}
    candidates = []
    whose_type_parameter!(candidates, T, FindType)
    return candidates
end

function whose_type_parameter!(candidates, ::Type{T}, ::Type{FindType}) where {T, FindType}
    if isempty(T.parameters)
        return
    else
        for (i, param) in enumerate(T.parameters)
            if param isa Type
                if param == FindType
                    push!(candidates, (nameof(T), i))
                else
                    whose_type_parameter!(candidates, param, FindType)
                end
            end
        end
    end
end

end # module


text:
# Assume that m ≫ n, so we want to throw out A
struct LassoADMMProblemData{T}
    ATA::AbstractMatrix{T}          # n x n matrix
    ATb::AbstractVector{T}          # n vector
    bTb::T                          # scalar
    m::Int
    n::Int
    function LassoADMMProblemData(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T <: Real}
        m, n = size(A)
        m != length(b) && error(DimensionMismatch("Invalid dimensions for A or b"))
        return new{T}(1/m * (A'*A), 1/m * (A'*b), 1/m * sum(x->x^2, b), m, n)
    end
end

#TODO: group metrics and params into vectors
struct LassoADMMProblem{T}
    data::LassoADMMProblemData{T}
    Ax::AbstractMatrix{T}       # matrix in x update
    bx::AbstractVector{T}       # vector in x update
    obj_val::MVector{1,T}       # objective (uses zk)
    sq_error::MVector{1,T}      # 0.5*||Ax - b||² (uses zk)
    dual_gap::MVector{1,T}      # duality gap
    xk::AbstractVector{T}       # primal var
    zk::AbstractVector{T}       # primal var
    uk::AbstractVector{T}       # dual var (ADMM)
    λ::MVector{1,T}
    ρ::MVector{1,T}
    α::MVector{1,T}
    function LassoADMMProblem(A::AbstractMatrix{T}, b::AbstractVector{T}, λ::T; ρ=1.0, α=1.5) where {T <: Real}
        n = size(A, 2)
        data = LassoADMMProblemData(A, b)
        return new{T}(
            data,
            similar(data.ATA),
            similar(data.ATb),
            SA[zero(T)],
            SA[zero(T)],
            SA[zero(T)],
            zeros(T, n),
            zeros(T, n),
            zeros(T, n),
            SA[λ],
            SA[ρ],
            SA[α]
        )
    end
end
#λmax = maximum(ATb) => x⋆ = 0
function reset!(prob::LassoADMMProblem{T}) where {T}
    n = length(prob.xk)
    prob.obj_val[1] = zero(T)
    prob.sq_error[1] = zero(T)
    prob.dual_gap[1] = zero(T)
    prob.xk .= zeros(T, n)
    prob.zk .= zeros(T, n)
    prob.uk .= zeros(T, n)
    return nothing
end

function update_λ!(prob::LassoADMMProblem{T}, λ::T) where {T}
    prob.λ[1] = λ
    return nothing
end

struct LassoADMMLog{T <: AbstractFloat}
    dual_gap::Union{AbstractVector{T}, Nothing}
    obj_val::Union{AbstractVector{T}, Nothing}
    iter_time::Union{AbstractVector{T}, Nothing}
    rp::Union{AbstractVector{T}, Nothing}
    rd::Union{AbstractVector{T}, Nothing}
    setup_time::T
    precond_time::T
    solve_time::T
end

struct LassoADMMResult{T}
    obj_val::T                 # primal objective
    sq_error::T                # 0.5*||Ax - b||²
    x::AbstractVector{T}       # primal soln
    dual_gap::T                       # duality gap
    # vT::AbstractVector{T}    # dual certificate
    log::LassoADMMLog
end

# --- substeps ---
# P is a preconditioner
function update_x!(prob::LassoADMMProblem{T}, solver::S, P) where {T <: Real, S <: CgSolver}
    # TODO: perhaps should decrease the ϵ here?
    # TODO: If A is fat, want to use matrix inversion lemma:
    #   (I + ρAᵀA)⁻¹ = I - ρAᵀ(I + ρAAᵀ)⁻¹A
    #   NOTE: should also adjust preconditioner to use AAᵀ
    @. prob.bx = prob.data.ATb + prob.ρ[1] * (prob.zk - prob.uk)
    cg!(solver, prob.Ax, prob.bx; M=P)
    !issolved(solver) && error("CG failed")
    prob.xk .= solver.x
end

function soft_threshold(x::T, y::T, κ::T) where {T <: Real}
    return max(zero(T), x + y - κ) - max(zero(T), -x - y - κ)
end

function update_z!(prob::LassoADMMProblem{T}; relax=true, xhat=nothing) where {T <: Real}
    if relax
        prob.zk .= soft_threshold.(xhat, prob.uk, prob.λ[1] / prob.ρ[1])
    else
        prob.zk .= soft_threshold.(prob.xk, prob.uk, prob.λ[1] / prob.ρ[1])
    end
end

function update_u!(prob::LassoADMMProblem{T}; relax=true, xhat=nothing) where {T <: Real}
    if relax
        @. prob.uk += xhat - prob.zk
    else
        @. prob.uk += prob.xk - prob.zk
    end
end

# updated_obj = true if
#   1) prob.sq_error has been updated
#   2) cache.v has been updated
function dual_gap!(prob::LassoADMMProblem{T}, cache; updated_obj=false) where {T}
    #TODO: benchmark norm(x, Inf) vs maximum(x)
    #TODO: benchmark norm(x, 1) vs sum(x->abs(x), x)
    x = prob.zk
    v = cache.v
    if !updated_obj
        mul!(v, prob.data.ATA, x)
    end
    u = prob.data.ATb
    w = prob.data.bTb
    @. cache.v_u = v - u
    z = norm(cache.v_u, Inf)
    uTx = dot(u, x)

    if !updated_obj
        prob.sq_error = 0.5*(dot(x, v) - 2uTx + w)
        prob.obj_val[1] = prob.sq_error[1] + prob.λ[1] * norm(x, 1)
    end

    prob.dual_gap[1] = prob.λ[1]^2/z^2*prob.sq_error[1] + prob.obj_val[1] + prob.λ[1]/z * (uTx - w)

    return prob.dual_gap[1]
end

function sq_error!(prob::LassoADMMProblem{T}, cache) where {T}
    x = prob.zk
    v = cache.v
    mul!(v, prob.data.ATA, x)

    uTx = dot(prob.data.ATb, x)

    prob.sq_error[1] = 0.5*(dot(x, v) - 2uTx + prob.data.bTb)
    prob.obj_val[1] = prob.sq_error[1] + prob.λ[1] * norm(x, 1)
    return prob.sq_error[1]
end

function update_rho!(prob::LassoADMMProblem, rp, rd, μ, τ_inc, τ_dec)
    if rp > μ * rd
        prob.ρ[1] = τ_inc * prob.ρ
        prob.Ax[diagind(prob.Ax)] .= prob.data.ATA[diagind(prob.data.ATA)] .+ prob.ρ[1]
        return prob.ρ[1]
    elseif rd > μ * rp
        prob.ρ[1] = prob.ρ / τ_dec
        prob.Ax[diagind(prob.Ax)] .= prob.data.ATA[diagind(prob.data.ATA)] .+ prob.ρ[1]
        return prob.ρ[1]
    else
        return prob.ρ[1]
    end
end

# --- main solver ---
function solve!(
    prob::LassoADMMProblem{T};
    relax::Bool=true,
    logging::Bool=false,
    precondition::Bool=true,
    tol=1e-6,
    max_iters::Int=100,
    print_iter::Int=25,
) where {T <: Real}
    setup_time_start = time_ns()
    @printf("Starting setup...")

    # --- parameters ---
    n = prob.data.n
    m = prob.data.m
    t = 1
    prob.dual_gap[1] = Inf
    r0 = 110
    μ = 10
    τ_inc = 2
    τ_dec = 2

    # --- enable multithreaded BLAS ---
    n_threads = Sys.CPU_THREADS
    BLAS.set_num_threads(n_threads ÷ 2)

    # TODO: scale problem data??

    # --- init cache and memory allocations ---
    ρ = prob.ρ[1]
    λ = prob.λ[1]
    α = prob.α[1]
    solver = CgSolver(n, n, typeof(prob.xk))
    if relax
        xhat = copy(prob.xk)                             # for relaxation
    end
    cache = (
        v=zeros(n),
        v_u=zeros(n)
    )
    prob.Ax .= prob.data.ATA
    prob.Ax[diagind(prob.Ax)] .+= ρ
    uk_old = similar(prob.uk)

    # --- Precondition ---
    if precondition
        @printf("\n\tPreconditioning...")
        precond_time_start = time_ns()
        ATA_nys = adaptive_nystrom_approx(prob.data.ATA, r0; q=10, tol=1e-7*n^2)
        P = RandomizedNystromPreconditionerInverse(ATA_nys, ρ)
        precond_time = (time_ns() - precond_time_start) / 1e9
        r = length(ATA_nys.Λ.diag)
        @printf("\n\tPreconditioned (rank %d) in %6.3fs", r, precond_time)
    else
        P = I
        precond_time = 0.0
    end

    # --- Logging ---
    if logging
        dual_gap_log = zeros(max_iters)
        obj_val_log = zeros(max_iters)
        iter_time_log = zeros(max_iters)
        rp_log = zeros(max_iters)
        rd_log = zeros(max_iters)
    else
        dual_gap_log = nothing
        obj_val_log = nothing
        iter_time_log = nothing
        rp_log = nothing
        rd_log = nothing
    end

    setup_time = (time_ns() - setup_time_start) / 1e9
    @printf("\nSetup in %6.3fs\n", setup_time)

    # --- Print Headers ---
    headers = ["Iteration", "Objective", "RMSE", "Dual Gap", "ρ", "Time"]
    print_header(headers)

    # --------------------------------------------------------------------------
    # --------------------- ITERATIONS -----------------------------------------
    # --------------------------------------------------------------------------
    solve_time_start = time_ns()
    #TODO: relative error for stopping criterion?? η/G(ν) ≤ ϵ
    while t <= max_iters && prob.dual_gap[1] > tol

        # --- Update Iterates ---
        update_x!(prob, solver, P)
        if relax
            @. xhat = α * prob.xk + (1-α) * prob.zk
        end
        update_z!(prob; relax=relax, xhat=xhat)
        uk_old .= prob.uk
        update_u!(prob; relax=relax, xhat=xhat)

        # Update ρ
        rp = norm(prob.xk - prob.zk)
        rd = norm(ρ*(prob.uk - uk_old))
        ρ = update_rho!(prob, rp, rd, μ, τ_inc, τ_dec)

        # --- Eval Termination Criterion
        sq_error!(prob, cache)
        dual_gap!(prob, cache; updated_obj=true)

        # --- Logging ---
        time_sec = (time_ns() - solve_time_start) / 1e9
        if logging
            dual_gap_log[t] = prob.dual_gap[1]
            obj_val_log[t] = prob.obj_val[1]
            iter_time_log[t] = time_sec
            rp_log[t] = rp
            rd_log[t] = rd
        end

        # --- Printing ---
        if t == 1 || t % print_iter == 0
            print_iter_func((
                string(t),
                prob.obj_val[1],
                sqrt(2*prob.sq_error[1]/m),
                prob.dual_gap[1],
                prob.ρ[1],
                time_sec
            ))
        end

        t += 1
    end

    # print final iteration if havent done so
    if (t-1) % print_iter != 0 && (t-1) != 1
        print_iter_func((
            string(t-1),
            prob.obj_val[1],
            sqrt(2*prob.sq_error[1]/m),
            prob.dual_gap[1],
            prob.ρ[1],
            (time_ns() - solve_time_start) / 1e9
        ))
    end
    solve_time = (time_ns() - solve_time_start) / 1e9
    @printf("\nSolved in %6.3fs, %d iterations\n", solve_time, t-1)
    @printf("Total time: %6.3fs\n", setup_time + solve_time)
    print_footer()

    # --- Construct Logs ---
    log = LassoADMMLog(
        dual_gap_log[1:t-1], obj_val_log[1:t-1], iter_time_log[1:t-1],
        rp_log[1:t-1], rd_log[1:t-1],
        setup_time, precond_time, solve_time
    )

    # --- Construct Solution ---
    res = LassoADMMResult(
        prob.obj_val[1],
        prob.sq_error[1],
        prob.xk,
        prob.dual_gap[1],
        log
    )

    return res, log

end

ext: jl


text:
using BeaconK8sUtilities
using Test, Logging, JSON3, Dates

function capture_logs(f)
    io = IOBuffer()
    Logging.with_logger(f, json_logger(; io))
    return JSON3.read(take!(io), Dict)
end

struct TestObj1
    f::Any
end

struct TestObj2
    f::Any
end

BeaconK8sUtilities.jsonable(obj::TestObj1) = (; obj.f)

@testset "BeaconK8sUtilities.jl" begin
    @testset "setup_tensorboard" begin
        testdir = mktempdir()
        test_args = (; app="test-app", logdir="test-logdir", ecr="test-ecr",
                     service_account="test-service-account", namespace="test-namespace")
        @test setup_tensorboard(testdir; test_args...) === nothing

        # Check we at least get the right files
        @test Set(readdir(testdir)) ==
              Set(["tensorboard.dockerfile", "tensorboard.sh", "tensorboard.yaml"])

        # No `overwrite=true`
        @test_throws ErrorException setup_tensorboard(testdir; test_args...)

        # Works with `overwrite=true`
        @test setup_tensorboard(testdir; test_args..., overwrite=true) === nothing
        @test Set(readdir(testdir)) ==
              Set(["tensorboard.dockerfile", "tensorboard.sh", "tensorboard.yaml"])
    end

    @testset "setup_follow" begin
        testdir = mktempdir()
        test_args = (; namespace="test-namespace")
        @test setup_follow(testdir; test_args...) === nothing

        # Check we at least get the right files
        @test readdir(testdir) == ["follow.sh"]

        # No `overwrite=true`
        @test_throws ErrorException setup_follow(testdir; test_args...)

        # Works with `overwrite=true`
        @test setup_follow(testdir; test_args..., overwrite=true) === nothing
        @test readdir(testdir) == ["follow.sh"]
    end

    @testset "json_logger" begin
        logs = capture_logs() do
            @info "Hi" x = Dict("a" => Dict("b" => 2))
        end
        @test logs["msg"] == "Hi"
        @test logs["kwargs"]["x"] == Dict("a" => Dict("b" => 2))
        @test logs["kwargs"]["worker_id"] == 1
        @test parse(DateTime, logs["kwargs"]["timestamp"],
                    dateformat"yyyy-mm-dd HH:MM:SS") isa DateTime

        logs = capture_logs() do
            @info "Hi" x = TestObj1("a")
        end
        @test logs["kwargs"]["x"] == Dict("f" => "a")

        logs = capture_logs() do
            @info "Hi" x = TestObj2("a")
        end
        @test logs["kwargs"]["x"] == "TestObj2(\"a\")"
        @test logs["kwargs"]["LoggingFormats.FormatError"] ==
              "ArgumentError: TestObj2 doesn't have a defined `StructTypes.StructType`"

        # test `maxlog`
        f() = @info "hi" maxlog = 1
        io = IOBuffer()
        Logging.with_logger(json_logger(; io)) do
            f()
            return f()
        end
        logs = String(take!(io))
        lines = split(logs, '\n')
        @test length(lines) == 2
        @test lines[2] == ""
        @test JSON3.read(lines[1])["msg"] == "hi"
    end
end


text:
using HalideCall
const libblur = "./libblur.so"
A = rand(Float32,2000,3000);
out = similar(A);
Abuf = HalideBuffer(A)
outbuf = HalideBuffer(out, (1,1))
fill!(out, 0);
ccall((:blur_twopasses,libblur), Void, (Ptr{Void}, Ptr{Void}), [Abuf], [outbuf])
@time (for i = 1:10; ccall((:blur_twopasses,libblur), Void, (Ptr{Void}, Ptr{Void}), [Abuf], [outbuf]); end)
@time (for i = 1:10; ccall((:blur_twopasses_vectorized,libblur), Void, (Ptr{Void}, Ptr{Void}), [Abuf], [outbuf]); end)
@time (for i = 1:10; ccall((:blur_tiled,libblur), Void, (Ptr{Void}, Ptr{Void}), [Abuf], [outbuf]); end)

using KernelTools

function blur_twopasses!(blur_xy, blur_y, A)
    interior_x, interior_y = 2:size(A,1)-1, 2:size(A,2)-1
    for y = interior_y, x = 1:size(A,1)
        @inbounds blur_y[x,y] = A[x,y-1]+A[x,y]+A[x,y+1]
    end
    for y = interior_y, x = interior_x
        @inbounds blur_xy[x,y] = (one(eltype(A))/9)*(blur_y[x-1,y] + blur_y[x,y] + blur_y[x+1,y])
    end
    blur_xy
end
blur_twopasses!(output, A) = blur_twopasses!(output, similar(A), A)
blur_twopasses(A) = blur_twopasses!(similar(A), similar(A), A)

function blur_tiled!(blur_xy, A)
    interior_x, interior_y = 2:size(A,1)-1, 2:size(A,2)-1
    KernelTools.@tile (x,1024,y,16) (@inbounds blur_x[x,y] = A[x-1,y]+A[x,y]+A[x+1,y]) begin
        for y = interior_y, x = interior_x
            @inbounds blur_xy[x,y] = (one(eltype(A))/9)*(blur_x[x,y-1] + blur_x[x,y] + blur_x[x,y+1])
        end
    end
    blur_xy
end
blur_tiled(A) = blur_tiled!(similar(A), A)

blur_twopasses!(out, A)
gc()
@time (for i = 1:10; blur_twopasses!(out, A); end)
tmp = similar(A)
gc()
@time (for i = 1:10; blur_twopasses!(out, tmp, A); end)
blur_tiled!(out, A)
@time (for i = 1:10; blur_tiled!(out, A); end)


text:
function mcP_Lv(Pk) #σᵥ² needed for P13 shortcut, no bias
    f = x -> exp(x)*Pk(exp(x))/(2π^2)
    quadgk(f,log(kmino), log(kmaxo))
end

function P13BiasTest(Pk)
    k, i1 = I1(Pk)
    _, i2 = I2(Pk)
    _, i3 = I3(Pk)
    p13 = @. Pk(k)*3*((2/63)*i1 + (1/42)*i2 - (1/18)*i3 - (1/18)*k^2*mcP_Lv(Pk)[1])
    return k , 2 .*p13
end

function P13BiasTest2(Pk)
    local13BiasDict = BiasDict13(b1=1, bη=0, bK²=0, bδη=0, bη²=0, bKKpara=0, btd=0,
                                 bδΠ2para=0, bηΠ2para=0, bΠ2Kpara=0, bΠ2para=0, bΠ3para=0)
    locallhdBiasDict = BiasDictlhd(b1=1, bη=0, P0_ϵ=0, b∇2δ=mcP_Lv(Pk)[1]/6, β∇2v=0,
                                   β∂2parav=0, P2_ϵ=0, P2_ϵε_η=0)
    k, p13EFT = P13ggs(Pk, [1.0], BiasOperatorDict13=local13BiasDict, f=0)
    pL13 = 2*p13EFT .+ Plhd(Pk, k, [1.0], BiasOperatorDictlhd=locallhdBiasDict, f=0)
    p13 = (pL13[:,1] .- Pk.(k))
    return k, p13
end

function P13BiasTest2f(Pk)
    local13BiasDict = BiasDict13(b1=1, bη=0, bK²=0, bδη=0, bη²=0, bKKpara=0, btd=0,
                                 bδΠ2para=0, bηΠ2para=0, bΠ2Kpara=0, bΠ2para=0, bΠ3para=0)
    locallhdBiasDict = BiasDictlhd(b1=1, bη=0, P0_ϵ=0, b∇2δ=mcP_Lv(Pk)[1]/6, β∇2v=0,
                                   β∂2parav=0, P2_ϵ=0, P2_ϵε_η=0)
    k, p13EFT = P13ggs(Pk, [1.0], BiasOperatorDict13=local13BiasDict, f=1)
    pL13 = 2*p13EFT .+ Plhd(Pk, k, [1.0], BiasOperatorDictlhd=locallhdBiasDict, f=1)
    p13 = (pL13[:,1] .- Pk.(k))
    return k, p13
end

function Itest(Pk)
    k, P1 = I1(Pk)
    _, P2 = I2(Pk)
    _, P3 = I3(Pk)
    _, P4 = I4(Pk)
    _, P5 = I5(Pk)
    plot(k,P1, label="I1", c="r")
    plot(k,P2, label="I2", c="b", ls="-.")
    plot(k,-P3, label="-I3", c="g", ls="--")
    plot(k, P4, label="I4", c="purple")
    plot(k, -P4, label="-I4", c="purple", ls="--")
    plot(k, P5, label="I5", c="orange", ls=":")
    xscale("log")
    yscale("log")
    legend()
    xlim(1e-3,1)
    ylim(1e-7,4)
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
nt = 1000
seed = 99
ssm = acqr(0.8,1.5,2,3,0.5,2)
m = acqr_simula(ssm,nt,seed)
y = m["y"]

#######################
start = acqr(0.5,0.5,1,1,0,0)
print("Testing acqr_em: ")
m1 = acqr_em(y,start,max_iter=10)
if sum(m1.fit.a^2) > 1e-6
    println("OK")
end


text:
# Computational cores of BBOB separable functions f1-f5

# f1 = Sphere function
sphere_function(x::AbstractVector{<:Number}) = sum(abs2, x, 1)


text:
#TODO: wrap this in a macro/function that takes priors,etc. as arguments
# Reminder that this might require using the esc function
# (See Julia manual metaprogramming section, hygiene subsection)
@model function fit_sam_model(data, prob, bkg, u0, k)
    σproc ~ Truncated(Normal(1,1), 0, 10)
    σobs ~ Truncated(Normal(2,2), 0, 10)
    α ~ Truncated(Normal(1.5,0.15), 1.3, 1.7)
    τ₊ ~ Truncated(Normal(900.0, 100.0), 0, 10000)
    τ₋ = 0.01
    # TODO: add initial condition uncertainty
    qPerUnitHZ = qCalc_powerLaw(1, τ₋, τ₊, α)
    p = [0.01, τ₊, α, k, qPerUnitHZ, 9.0, 18.0]
    h(p,t) = bkg
    prob = remake(prob, u0 = [u0], h = h, p = p)
    predicted = solve(prob, MethodOfSteps(RK4()))
    z = zeros(size(data)[2])
    for i = 1:size(data)[2]
        z[i] ~ Normal(predicted(data[1,i])[1],  σproc)
        data[2,i] ~ Normal( z[i],σobs)
    end
end

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
using MvrtCopula
using Test

include("./archmCopulaTest.jl")


text:
#creates a model using flux that uses levels (as opposed to differences)
#WARNING- DO NOT DELETE BELOW
function fluxlevelmlemodel(panel::AbstractDataFrame,  zs::ZSpec,
    ::Type{T}=PARAM[:fluxtype],
    ::Type{TV}=PARAM[:fluxgpu] ? CUDA.CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:fluxgpu] ? CUDA.CuMatrix{T} : Matrix{T},
    PredictionType::Val = Val(PARAM[:predictiontype])) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}

  PARAM[:predictiontype] ∈ [:level, :levellog] || error(
    "invalid predictiontype $PredictionType for model type $(PARAM[:fluxmodel])")

  #(PARAM[:fluxgpu] && (T≠Float32)) && error("use Float32 for fluxtype when using gpu")

  CUDA.allowscalar(false)
  Zms::MeasureSpec = zs.Zms

  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, Zms, T, TV, TM,
    weightdata = (Fw=Fw, FRLw=FRLw),
    Fvol=Zms.Fvol)
  dims = p.dims
  tidx::Dict = p.tidx
  cleanpanel::SubDataFrame = p.cleanpanel
  Ncleanpanel::Int = nrow(cleanpanel)
  ws::TM = p.dat[:Fw]
  RLws::TM = p.dat[:FRLw]
  ts::Vector{Int} = p.ts
  v::TV = p.dat[:Fvol]
  X = (ws=ws, RLws=RLws)
  #construct the parameter
  Θ = VolumePartsMLE(dims.T, dims.K, ts, T)
  Π=Flux.params(Θ)

  #define the loss function (in this case, log-likelihood)
  #Npanelones = ones(Npanel) |> fluxgpu
  iter::Int = 0
  function loss(X::NamedTuple, y::AbstractVector)
    θ = view(exp.(Θ.τ),Θ.expand.τidx)
    #keep count of each iteration
    iter += 1
    absμ = Θ(X, y, PredictionType)

    #θ = exp.(τlong) #here τ=logθ
    k = absμ ./ θ  #absμ ./ θ # μ = kθ,

    if (sum(θ .> 0) ≠ Npanel) || (sum(k .> 0) ≠ Npanel)
      return Inf #this will cause the program to stop
    end

    loglike = sum(vecgammalogpdf(k, θ, y))

    return -loglike
  end

  #NOTE:" cuda optimized version
  function loss(X::NamedTuple, y::CuVector)
    θ = CuVector(exp.(Θ.τ)[Θ.expand.τidx])
    #keep count of each iteration
    iter += 1
    absμ = Θ(X, y, PredictionType)

    #θ = CUDA.exp_fast.(τlong) #here τ=logθ
    k = absμ ./ θ # CUDA.div_fast  crashes for some reason

    if (CUDA.sum(θ .> 0) ≠ Ncleanpanel) || (CUDA.sum(k .> 0) ≠ Ncleanpanel)
      return Inf #this will cause the program to stop
    end

    loglike = CUDA.sum(vecgammalogpdf(k, θ, y))

    #println("minval: $minval, loglike: $loglike")
    return -loglike
  end

  #use this function to throw data from the main block over to
  #a seperate function to report error information
  debugcb() = debuginf(cleanpanel, tidx,  Θ, Matrix(ws), Matrix(RLws), Vector(v))

  @info "computing test gradient. Testing function: $(sum(Θ(X, v, Val{:level}())))"
  gs = gradient(()->loss(X, v),Π)
  #@info "gs[A]: $(gs[Θ.A])"

  test∇VolumeParts!(loss, X, v, Θ)

  @info "beginning training with $(length(Π)) parameter vectors"

  λ⁺=solvevol!(OPT_INDEX[PARAM[:fluxopt]], T;
    loss=loss,
    Θ=Θ,
    X=X, y=v,
    debugcb = debugcb)

  results::DataFrame = formresults(panel, Zms, Θ, PredictionType)

  return (Int(round(λ⁺)), results)
end

function cudadifsq(a,b)
  Δ = a - b
  return Δ*Δ
end

function fluxlevellsqmodel(panel::AbstractDataFrame,  zs::ZSpec,
    ::Type{T}=PARAM[:fluxtype],
    ::Type{TV}=PARAM[:fluxgpu] ? CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:fluxgpu] ? CuMatrix{T} : Matrix{T};
    #::Type{PredictionType} = Val(PARAM[:predictiontype])
    PredictionType::Val = Val(PARAM[:predictiontype])#Capacity.eval(:(Val(PARAM[:predictiontype])))
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}
  PARAM[:predictiontype] ∈ [:level, :levellog] || error(
    "invalid predictiontype $PredictionType for model type $(PARAM[:fluxmodel])")
  #(PARAM[:fluxgpu] && (T≠Float32)) && error("use Float32 for fluxtype when using gpu")

  CUDA.allowscalar(false)
  Zms::MeasureSpec = zs.Zms

  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, Zms,T, TV, TM,
    weightdata = (Fw=Fw, FRLw=FRLw),
    Fvol=Zms.Fvol)
  dims = p.dims
  tidx::Dict = p.tidx
  cleanpanel::SubDataFrame = p.cleanpanel
  Ncleanpanel::Int = nrow(cleanpanel)
  ws::TM = p.dat[:Fw]
  RLws::TM = p.dat[:FRLw]
  ts::Vector{Int} = p.ts
  X = (ws=ws, RLws=RLws)
  v::TV = p.dat[:Fvol]

  #construct the parameter
  Θ = VolumePartsLSq(dims.T, dims.K, ts, T)
  Π=Flux.params(Θ)
  length(Π) == 0 && error("No parameters noticed by Flux!")

  #define the loss function (in this case, log-likelihood)
  #Npanelones = ones(Npanel) |> fluxgpu
  #ptype2 = Val{:levellog}()

  function loss(X::NamedTuple,y::AbstractVector) where TX<:AbstractMatrix

    #keep count of each iteration
    absμ = Θ(X, y, PredictionType)

    sumsq = sum((absμ .- y) .^ 2)

    return sumsq
  end

  #NOTE:" cuda optimized version
  function loss(X::NamedTuple, y::CuVector)

    #keep count of each iteration
    absμ = Θ(X, y, PredictionType)

    sumsq = CUDA.sum(cudadifsq.(absμ,y))

    #println("minval: $minval, loglike: $loglike")
    return sumsq
  end

  #use this function to throw data from the main block over to
  #a seperate function to report error information
  debugcb() = debuginf(cleanpanel, tidx,  Θ, Matrix(ws), Matrix(RLws), Vector(v))

  @info "computing test gradient. Testing function: $(sum(Θ(X, v, PredictionType)))"

  test∇VolumeParts!(loss, X, v, Θ)

  @info "beginning training with $(length(Π)) parameter vectors"

  λ⁺=solvevol!(OPT_INDEX[PARAM[:fluxopt]], T;
    loss=loss,
    Θ=Θ,
    X=X,
    y=v,
    debugcb = debugcb)

  results::DataFrame = formresults(panel, Zms, Θ, PredictionType)

  return Int(round(λ⁺)), results
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using GAP, Hecke, Printf, ArgParse

include(joinpath(Hecke.pkgdir,"examples/FieldEnumeration/FieldEnumeration.jl"))

###############################################################################
#
#  Read-Write
#
###############################################################################

function ArgParse.parse_item(::Type{fmpz}, x::AbstractString)
  if in('^', x)
    l = split(x, '^')
    if length(l) != 2
      throw(error("Could not parse $x as fmpz"))
    end
    l = string.(l)
    return (parse(fmpz, l[1]))^parse(Int, l[2])
  else
    return parse(fmpz, string(x))
  end
end

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--order", "-o"
      help = "Order of the group"
      arg_type = Int
      default = -1
    "--id", "-i"
      help = "Id of the group"
      arg_type = Int
      default = -1
    "--number", "-n"
      help = "Number of the group"
      arg_type = Int
      default = -1
    "--tamely-ramified"
      help = "Only tamely ramified fields"
      action = :store_true
    "--only-real"
      help = "Only totally real fields"
      action = :store_true
    "--only-complex"
      help = "Only totally complex fields"
      action = :store_true
    "--disc-bound"
      help = "Discriminant bound"
      arg_type = fmpz
    "--rootdisc-bound"
      help = "Discriminant bound"
      arg_type = Float64
    "--only-cm"
      help = "Only CM fields"
      action = :store_true
    "--max-ab-subfields"
      help = "File containing maximal abelian subextensions"
    "--simplify"
      help = "Simplify the field"
      action = :store_true
    "--output"
      help = "Output file"
      arg_type = String
      default = ""
    "--output-type"
      help = "Type of output file. Options are default (uses _write_fields), polys, nfdb"
      default = "default"
    "--silent"
      help = "No verbose out"
      action = :store_true
  end

  return parse_args(s)
end

function main()
  parsed_args = parse_commandline()

  local grp_order::Int
  local grp_id::Int
  local dbound::Union{fmpz, Nothing}
  local only_real::Bool
  local only_complex::Bool
  local only_tame::Bool
  local grp_no::Int
  local only_cm::Bool
  local maxabsubfields::Union{String, Nothing}
  local simplify::Bool
  local output::String
  local rdbound::Union{Float64, Nothing}
  local silent::Bool

  for (arg, val) in parsed_args
    #println("$arg => $val")
    if arg == "order"
      grp_order = val
    elseif arg == "id"
      grp_id = val
    elseif arg == "disc-bound"
      dbound = val
    elseif arg == "rootdisc-bound"
      rdbound = val
    elseif arg == "only-real"
      only_real = val
    elseif arg == "only-complex"
      only_complex = val
    elseif arg == "tamely-ramified"
      only_tame = val
    elseif arg == "number"
      grp_no = val
    elseif arg == "only-cm"
      only_cm = val
    elseif arg == "max-ab-subfields"
      maxabsubfields = val
    elseif arg == "simplify"
      simplify = val
    elseif arg == "silent"
      silent = val
    elseif arg == "output"
      output = val
    end
  end

  sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

  if grp_id != -1
    n = grp_order
    i = grp_id
  end

  if !silent
    @info "Small group ($n, $i)"

    if dbound !== nothing
      @info "Discriminant bound $dbound"
    end

    if rdbound !== nothing
      @info "Root discriminant bound $rdbound"
    end

    @info "Simplifying fields: $simplify"

    if output == ""
      output = "fields_$(i)_$(n).nfdb"
    end

    if isfile(output)
      error("Output file $output already exists. You may specify the output file by setting --output")
    end

    @info "Maximal abelian fields: $maxabsubfields"
    @info "Only real fields: $only_real"
    @info "Only complex fields: $only_complex"
    @info "Only CM: $only_cm"

    @info "Output to: $output"

    set_verbose_level(:FieldsNonFancy, 3)
  end

  if dbound == nothing && rdbound == nothing
    error("One of --disc-bound or --rootdisc-bound must be specified")
  elseif dbound != nothing && rdbound != nothing
    error("Only one of --disc-bound or --rootdisc-bound can be set")
  end

  if only_complex && only_real
    error("Only one of --only-real or --only-complex can be set")
  end

  if dbound === nothing
    dbound = fmpz(BigInt(ceil(BigFloat(rdbound)^(n))))
    if !silent
      @info "Translated root discriminant bound to: $dbound"
    end
  end

  if grp_id != -1
    @assert grp_order != -1
    n = grp_order
    i = grp_id
    grp_no = findfirst(isequal((n, i)), small_solvable_groups)
  else
    @assert grp_no != -1
    n, i = small_solvable_groups[grp_no]
  end

  flush(stdout)

  if maxabsubfields isa String
    maxabsub = Hecke._read_from_file(maxabsubfields)
    l = fields(n, i, maxabsub, dbound, only_real = only_real)
  else
    l = fields(n, i, dbound, only_real = only_real)
  end

  flush(stdout)

  # Determine the automorphism groups
  if only_cm
    Hecke.assure_automorphisms.(l)
  end

  flush(stdout)

  ll = map(v -> v.field, l)

  flush(stdout)

  if only_cm
    ffields = [ (x, discriminant(maximal_order(x))) for x in ll if Hecke.iscm_field(x)[1]]
  else
    ffields = [ (x, discriminant(maximal_order(x))) for x in ll ]
  end

  flush(stdout)

  sort!(ffields, lt = (x, y) -> abs(x[2]) <= abs(y[2]))

  res = eltype(Hecke.NFDB{1})[]

  for (x, y) in ffields
    if simplify
      x, _ = Hecke.simplify(x)
    end
    r = Hecke._create_record(x)
    r[:discriminant] = y
    push!(res, r)
  end

  DBout = Hecke.NFDB(res)
  t = "fields($n, $i, $dbound, only_real = $only_real, simplify = $simplify)"
  if rdbound !== nothing
    t = t * " (rd bound $rdbound)"
  end
  if maxabsubfields isa String
    t = t * " with maximal abelian subfields from $maxabsubfields"
  end
  Hecke.add_meta!(DBout, :title => t)
  open(output, "w") do f
    Base.write(f, DBout)
  end
end

main()


text:
function ThomasLU(A)

  #A = the matrix of coefficients 'A' must be squared

    m,n =size(A)

    if m != n
        error("Matrix must be squared")
    end
    e=zeros(m-1)
    f=zeros(m)
    g=zeros(m-1)
    f[1]=A[1,1]
    for i=2:m
        e[i-1]=A[i,i-1]/f[i-1]
        g[i-1]=A[i-1,i]
        f[i]=A[i,i]-e[i-1]*g[i-1]

    end

    L=Tridiagonal(e,ones(m),zeros(m-1))
    U=Tridiagonal(zeros(m-1),f,g)

    return L, U

end


text:
@testset "CompoundRatingCurve type" begin

    G₁ = Gauging.([1,2],[3,4])
    rc₁ = RatingCurve(G₁, 1, 2, 3)
    G₂ = Gauging.([5,6],[7,8])
    rc₂ = RatingCurve(G₂, 4, 5, 6)
    crc = CompoundRatingCurve([4.5], [rc₁, rc₂])

    @test crc.threshold ≈ [4.5] atol = sqrt(eps())
    @test crc.component[1] == rc₁
    @test crc.component[2] == rc₂

    # Test the promotion of the threshold
    crc = CompoundRatingCurve([4], [rc₁, rc₂])
    @test typeof(crc.threshold[]) == Float64

    # Test for the wrong number of threholds
    @test_throws ErrorException CompoundRatingCurve([4.5, 5], [rc₁, rc₂])

end

@testset "bic CompoundRatingCurve" begin
    #TODO
end

@testset "cint CompoundRatingCurve" begin
    #TODO
end

@testset "discharge with compound rating curve" begin
    rc₁ = RatingCurve(Gauging[], 1, 0, 3)
    rc₂ = RatingCurve(Gauging[], 3, 1, 2)
    crc = CompoundRatingCurve([1], [rc₁, rc₂])
    @test discharge(crc,1/2) ≈ 1/8
    @test discharge(crc,1) ≈ 0.0
    @test discharge(crc,2) ≈ 3.0
end

@testset "level with compound rating curve" begin
    rc₁ = RatingCurve(Gauging[], 1, 0, 3)
    rc₂ = RatingCurve(Gauging[], 3, 1, 2)
    crc = CompoundRatingCurve([1], [rc₁, rc₂])
    @test level(crc, 3) ≈ 2.0
    @test level(crc,1/8) ≈ 1/2
    @test level(crc,12) ≈ 3
end

@testset "logdischarge with compound rating curve" begin
    rc₁ = RatingCurve(Gauging[], 1, 0, 3)
    rc₂ = RatingCurve(Gauging[], 3, 1, 2)
    crc = CompoundRatingCurve([1], [rc₁, rc₂])
    @test logdischarge(crc,1/2) ≈ log(1/8)
    @test logdischarge(crc,2) ≈ log(3)
end

@testset "pint CompoundRatingCurve" begin
    #TODO
end

@testset "pintlog CompoundRatingCurve" begin
    #TODO
end

@testset "sse of a CompoundRatingCruver" begin
    h₁ = exp.([1/5, 1/4, 1/2, 1/2])
    q₁ = exp.([8/5, 7/4, 5/2, 3])

    G₁ = Gauging.(h₁, q₁)

    rc₁ = RatingCurve(G₁, exp(1), 0, 3)

    h₂ = exp.([1, 3/2, 2, 2])
    q₂ = exp.([5/2, 7/2, 9/2, 5])

    G₂ = Gauging.(h₂, q₂)

    rc₂ = RatingCurve(G₂, exp(1/2), 0, 2)

    crc = CompoundRatingCurve([exp(3/4)], [rc₁, rc₂])

    SSE = RatingCurves.sse(crc)

    @test SSE[1] ≈ 1/4 atol=sqrt(eps())
    @test SSE[2] ≈ 1/4 atol=sqrt(eps())

end

@testset "var of RatingCurve" begin
    h₁ = exp.([1/5, 1/4, 1/2, 1/2])
    q₁ = exp.([8/5, 7/4, 5/2, 3])

    G₁ = Gauging.(h₁, q₁)

    rc₁ = RatingCurve(G₁, exp(1), 0, 3)

    h₂ = exp.([1, 3/2, 2, 2])
    q₂ = exp.([5/2, 7/2, 9/2, 5])

    G₂ = Gauging.(h₂, q₂)

    rc₂ = RatingCurve(G₂, exp(1/2), 0, 2)

    crc = CompoundRatingCurve([exp(3/4)], [rc₁, rc₂])

    σ̂² = RatingCurves.var(crc)

    @test σ̂²[1] ≈ 1/4 atol=sqrt(eps())
    @test σ̂²[2] ≈ 1/4 atol=sqrt(eps())

end


text:
n, p = 500, 5
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p; seed=525)

# NOTE: in these cases, if available, θ_newton is used as reference.
# so if the others are not far, it's assumed all have converged to the
# same (correct) minimizer. The reference J(θ) or J(θ1) is the objective
# at the generating θ; it should usually be >
# these are tests that should be considered as "sanity check" and are not
# stress tests testing corner cases.

@testset "HuberReg" begin
    # No intercept
    δ = 0.01
    λ = 3.0
    hr = HuberRegression(δ, λ; fit_intercept=false)
    J = objective(hr, X, y)
    o = RobustLoss(Huber(δ)) + λ * L2Penalty()
    @test J(θ) ≈ o(y, X*θ, θ)
    θ_newton   = fit(hr, X, y, solver=Newton())
    θ_newtoncg = fit(hr, X, y, solver=NewtonCG())
    θ_lbfgs    = fit(hr, X, y, solver=LBFGS())
    θ_iwls     = fit(hr, X, y, solver=IWLSCG())
    @test isapprox(J(θ),          5.456121, rtol=1e-5)
    @test isapprox(J(θ_newton),   4.87426,  rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 4.87426,  rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    4.87426,  rtol=1e-5)
    @test isapprox(J(θ_iwls),     4.87426,  rtol=1e-5)

    δ = 0.01
    λ = 3.0
    hr = HuberRegression(δ, λ; penalize_intercept=true)
    J = objective(hr, X, y1)
    o = RobustLoss(Huber(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1, X1*θ1, θ1)
    θ_newton   = fit(hr, X, y1, solver=Newton())
    θ_newtoncg = fit(hr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(hr, X, y1, solver=LBFGS())
    θ_iwls     = fit(hr, X, y1, solver=IWLSCG())
    @test isapprox(J(θ1),         7.601925, rtol=1e-5)
    @test isapprox(J(θ_newton),   6.259785, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 6.259785, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    6.259785, rtol=1e-5)
    @test isapprox(J(θ_iwls),     6.259785, rtol=1e-5)

    # don't penalize intercept
    δ = 0.01
    λ = 3.0
    hr = HuberRegression(δ, λ)
    J = objective(hr, X, y1)
    θ_newton   = fit(hr, X, y1, solver=Newton())
    θ_newtoncg = fit(hr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(hr, X, y1, solver=LBFGS())
    θ_iwls     = fit(hr, X, y1, solver=IWLSCG())
    @test isapprox(J(θ1),         7.536531, rtol=1e-5)
    @test isapprox(J(θ_newton),   6.200183, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 6.200183, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    6.200183, rtol=1e-5)
    @test isapprox(J(θ_iwls),     6.200183, rtol=1e-5)
end

# NOTE: small difference obtained with NCG and IWLS, might just be
# that the algorithm stops a bit earlier, negligible relative tol anyway.
@testset "AndrewsReg" begin
    δ = 0.1
    λ = 3.0
    rr = RobustRegression(rho=Andrews(δ), lambda=λ; penalize_intercept=true)
    J = objective(rr, X, y1)
    o = RobustLoss(AndrewsRho(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1, X1*θ1, θ1)
    θ_newton   = fit(rr, X, y1, solver=Newton())
    θ_newtoncg = fit(rr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(rr, X, y1, solver=LBFGS())
    θ_iwls     = fit(rr, X, y1, solver=IWLSCG())
    @test isapprox(J(θ1),         7.359476, rtol=1e-5)
    @test isapprox(J(θ_newton),   0.486388, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 0.486388, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    0.486388, rtol=1e-5)
    @test isapprox(J(θ_iwls),     0.486388, rtol=1e-5)
end

@testset "BisquareReg" begin
    δ = 0.1
    λ = 3.0
    rr = RobustRegression(rho=Bisquare(δ), lambda=λ; penalize_intercept=true)
    J = objective(rr, X, y1)
    o = RobustLoss(BisquareRho(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1, X1*θ1, θ1)
    θ_newton   = fit(rr, X, y1, solver=Newton())
    θ_newtoncg = fit(rr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(rr, X, y1, solver=LBFGS())
    θ_iwls     = fit(rr, X, y1, solver=IWLSCG())
    @test isapprox(J(θ1),         7.773512, rtol=1e-5)
    @test isapprox(J(θ_newton),   0.818677, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 0.818677, rtol=1e-4)
    @test isapprox(J(θ_lbfgs),    0.818677, rtol=1e-4)
    @test isapprox(J(θ_iwls),     0.818677, rtol=1e-4)
end

@testset "LogisticReg" begin
    δ = 1.5
    λ = 1.0
    rr = RobustRegression(rho=Logistic(delta=δ), lambda=λ, penalize_intercept=true)
    J = objective(rr, X, y1)
    o = RobustLoss(LogisticRho(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1, X1*θ1, θ1)
    θ_newton = fit(rr, X, y1, solver=Newton())
    θ_newtoncg = fit(rr, X, y1, solver=NewtonCG())
    θ_lbfgs = fit(rr, X, y1, solver=LBFGS())
    θ_iwls  = fit(rr, X, y1, solver=IWLSCG())
    @test isapprox(J(θ1),         5.00079, rtol=1e-5)
    @test isapprox(J(θ_newton),   4.98146, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 4.98146, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    4.98146, rtol=1e-5)
    @test isapprox(J(θ_iwls),     4.98146, rtol=1e-5)
end

@testset "FairReg" begin
    δ = 0.1
    λ = 3.0
    rr = RobustRegression(rho=Fair(δ), lambda=λ, penalize_intercept=true)
    J = objective(rr, X, y1)
    o = RobustLoss(Fair(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1, X1*θ1, θ1)
    θ_newton   = fit(rr, X, y1, solver=Newton())
    θ_newtoncg = fit(rr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(rr, X, y1, solver=LBFGS())
    θ_iwls     = fit(rr, X, y1, solver=IWLSCG())
    @test isapprox(J(θ1),         8.586636, rtol=1e-5)
    @test isapprox(J(θ_newton),   8.484584, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 8.484584, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    8.484584, rtol=1e-5)
    @test isapprox(J(θ_iwls),     8.484584, rtol=1e-5)
end

# NOTE: small difference obtained with NCG and IWLS, might just be
# that the algorithm stops a bit earlier, negligible relative tol anyway.
@testset "TalwarReg" begin
    δ = 0.1
    λ = 3.0
    rr = RobustRegression(rho=Talwar(δ), lambda=λ, penalize_intercept=true)
    J = objective(rr, X, y1)
    o = RobustLoss(Talwar(δ)) + λ * L2Penalty()
    @test J(θ1) ≈ o(y1, X1*θ1, θ1)
    θ_newton   = fit(rr, X, y1, solver=Newton())
    θ_newtoncg = fit(rr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(rr, X, y1, solver=LBFGS())
    θ_iwls     = fit(rr, X, y1, solver=IWLSCG())
    @test isapprox(J(θ1),          8.564334, rtol=1e-5)
    @test isapprox(J(θ_newton),    2.43470,  rtol=1e-5)
    @test isapprox(J(θ_newtoncg),  2.43803,  rtol=1e-5)
    @test isapprox(J(θ_lbfgs),     2.43202,  rtol=1e-5)
    @test isapprox(J(θ_iwls),      2.43803,  rtol=1e-5)
end

###########################
## With Sparsity penalty ##
###########################

n, p = 500, 100
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=51112, sparse=0.1)
# pepper with outliers
y1a  = outlify(y1, 0.1)

@testset "Robust+L1" begin
    δ  = 0.1
    λ  = 5.0
    γ  = 10.0

    # this is a nice case where the robust is actually smooth
    rr = RobustRegression(rho=Huber(δ), lambda=λ, gamma=γ, penalize_intercept=true)
    J  = objective(rr, X, y1a)
    θ_ls    = X1 \ y1a
    θ_fista = fit(rr, X, y1a, solver=FISTA())
    θ_ista  = fit(rr, X, y1a, solver=ISTA())
    @test isapprox(J(θ_ls),    453.12684, rtol=1e-5)
    @test isapprox(J(θ_fista), 124.20330, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_ista),  124.20330, rtol=1e-5) # ista stops a bit early?
    @test nnz(θ_ls)    == 101
    @test nnz(θ_fista) == 4
    @test nnz(θ_ista)  == 4
end


text:
# All or most of the code in this file is taken from:
# 1. https://github.com/NHDaly/DeepcopyModules.jl (license: MIT)
# 2. https://github.com/perrutquist/CodeTransformation.jl (license: MIT)

"""
    jl_method_def(argdata, ci, mod) - C function wrapper

This is a wrapper of the C function with the same name, found in the Julia
source tree at julia/src/method.c

Use `addmethod!` or `codetransform!` instead of calling this function directly.
"""
jl_method_def(argdata::Core.SimpleVector, ci::Core.CodeInfo, mod::Module) =
    ccall(:jl_method_def, Cvoid, (Core.SimpleVector, Any, Ptr{Module}), argdata, ci, pointer_from_objref(mod))
# `argdata` is `Core.svec(Core.svec(types...), Core.svec(typevars...))`

"Recursively get the typevars from a `UnionAll` type"
typevars(T::UnionAll) = (T.var, typevars(T.body)...)
typevars(T::DataType) = ()

@nospecialize # the below functions need not specialize on arguments

"Get the module of a function"
getmodule(F::Type{<:Function}) = F.name.mt.module
getmodule(f::Function) = getmodule(typeof(f))

"Create a call singature"
makesig(f::Function, args) = Tuple{typeof(f), args...}

"""
    argdata(sig[, f])

Turn a call signature into the 'argdata' `Core.svec` that `jl_method_def` uses
When a function is given in the second argument, it replaces the one in the
call signature.
"""
argdata(sig) = Core.svec(Base.unwrap_unionall(sig).parameters::Core.SimpleVector, Core.svec(typevars(sig)...))
argdata(sig, f::Function) = Core.svec(Core.svec(typeof(f), Base.unwrap_unionall(sig).parameters[2:end]...), Core.svec(typevars(sig)...))

"""
    addmethod!(f, argtypes, ci)

Add a method to a function.

The types of the arguments is given as a `Tuple`.

Example:
```
g(x) = x + 13
ci = code_lowered(g)[1]
function f end
addmethod!(f, (Any,), ci)
f(1) # returns 14
```
"""
addmethod!(f::Function, argtypes::Tuple, ci::Core.CodeInfo) = addmethod!(makesig(f, argtypes), ci)
"""
    addmethod(sig, ci)

Alternative syntax where the call signature is a `Tuple` type.

Example:
```
addmethod!(Tuple{typeof(f), Any}, ci)
```
"""
function addmethod!(sig::Type{<:Tuple{F, Vararg}}, ci::Core.CodeInfo) where {F<:Function}
    jl_method_def(argdata(sig), ci, getmodule(F))
end

@specialize # restore default

"""
    codetransform!(tr, dst, src)

Apply a code transformation function `tr` on the methods of a function `src`,
adding the transformed methods to another function `dst`.

Example: Search-and-replace a constant in a function.
```
g(x) = x + 13
function e end
codetransform!(g => e) do ci
    for ex in ci.code
        if ex isa Expr
            map!(x -> x === 13 ? 7 : x, ex.args, ex.args)
        end
    end
    ci
end
e(1) # returns 8
```
"""
function codetransform!(tr::Function, @nospecialize(dst::Function), @nospecialize(src::Function))
    mod = getmodule(dst)
    for m in methods(src).ms
        ci = Base.uncompressed_ast(m)
        ci = tr(ci)
        jl_method_def(argdata(m.sig, dst), ci, mod)
    end
end
"Alternative syntax: codetransform!(tr, src => dst)"
codetransform!(tr::Function, @nospecialize(p::Pair{<:Function, <:Function})) =
    codetransform!(tr, p.second, p.first)


text:
cd("/Users/federicoclaudi/Documents/Github/LocomotionControl/analysis/behavior")
import Pkg;
Pkg.activate(".")

include("/Users/federicoclaudi/Documents/Github/LocomotionControl/analysis/behavior/multiple_horizons/simulate.jl")

"""
Run a simulation in which the model can only plan for `planning_horizon` seconds ahead.
"""
function run_simulation_and_plot(; planning_horizon::Float64=.5, n_iter=480, Δt=.01)
    # run global solution
    track = Track(;start_waypoint=2, keep_n_waypoints=-1)

    _, bike, _, globalsolution = run_mtm(
        :dynamics,
        3;
        showtrials=nothing,
        track=track,
        n_iter=5000,
        fcond=State(; u=30, ω=0),
        timed=false,
        showplots=false,
    )

    # plot background & global trajectory
    colors = [
            range(HSL(colorant"red"), stop=HSL(colorant"green"), length=(Int ∘ floor)(n_iter/2))...,
            range(HSL(colorant"green"), stop=HSL(colorant"blue"), length=(Int ∘ ceil)(n_iter/2))...
    ]

    simtracker = SimTracker(Δt=Δt)
    p1 = draw(:arena)
    draw!(FULLTRACK; border_alpha=.25, alpha=0.0)

    for i in pbar(1:n_iter, redirectstdout=false, description="Running horizon: $planning_horizon seconds", expand=true)
        simtracker.iter = i

        # run simulation and store results
        initial_state, solution, shouldstop, track = step(simtracker, globalsolution, planning_horizon)
        shouldstop && break

        # plot stuff

        # draw!(track; color=colors[i], border_lw=5, alpha=0.0)
        i%5 == 0 && plot_bike_trajectory!(solution, bike; showbike=false, label=nothing, color=colors[i], alpha=.8, lw=4)

        if i%10 == 0
            draw!(initial_state; color=colors[i], alpha=1)
        end
        # plot!(; title="T: $(round(simtracker.t; digits=2))s | horizon: $planning_horizon s")

        simtracker.s > 260 && break
    end

    display(p1)
end

run_simulation_and_plot(planning_horizon=.75)
println("done")

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
function  fft2(fx, deltx::Real)
	d = 2.0
	c = complex( (deltx / √(2.0 * π))^d )
   fk = fft(fx)
   scale!(fk, c)
   fk
end
function  ifft2(fk, deltk::Real)
	d = 2.0
	c = (deltk / √(2.0 * π))^d
	fx = bfft(fk)
	scale!(fx, c)
	fx
end
function  ifft2r(fk, deltk::Real)
	d = 2.0
	c = (deltk / √(2.0 * π))^d
	fx = real(bfft(fk))
	scale!(fx, c)
	fx
end
fft2(fx, par::SpectrumGrids) = fft2(fx, par.grd.deltx)
ifft2(fk, par::SpectrumGrids) = ifft2(fk, par.grd.deltk)
ifft2r(fk, par::SpectrumGrids) = ifft2r(fk, par.grd.deltk)
function  rft(fx, par::SpectrumGrids)
  c = complex( (par.grd.deltx / √(2.0 * π))^2.0 )
  fk = rfft(fx)
  scale!(fk, c)
  fk
end
function  irft(fk, par::SpectrumGrids)
  c = (par.grd.deltk / √(2.0 * π))^2.0
  nint = int(par.grd.n)
  fx = brfft(fk, nint)
  scale!(fx, c)
  fx
end

max_stars_repo_name: EthanAnderes/BayesianCmbLensing

max_issues_repo_name: EthanAnderes/BayesianCmbLensing

max_forks_repo_name: EthanAnderes/BayesianCmbLensing


text:
using Test

using CompScienceMeshes

for T in [Float32, Float64]
    e0 = point(T,0,0,0)
    e1 = point(T,1,0,0)
    e2 = point(T,0,1,0)
    e3 = point(T,0,0,1)

    local p = simplex(e0,e1,e2)

    @test isinside(p, (e0+e1+e2)/3) == true
    @test isinside(p, point(T, -0.1, 0.3, 0.0)) == false
    @test isinside(p, point(T, 0.3, -0.3, 0.0)) == false
    @test isinside(p, point(T, 0.6, 0.6, 0.0)) == false
    @test isinside(p, point(T, 0.4, 0.6, 0.0)) == false

    @test isinclosure(p, (e0+e1+e2)/3) == true
    @test isinclosure(p, point(T, -0.1, 0.3, 0.0)) == false
    @test isinclosure(p, point(T, 0.3, -0.3, 0.0)) == false
    @test isinclosure(p, point(T, 0.6, 0.6, 0.0)) == false
    @test isinclosure(p, point(T, 0.4, 0.6, 0.0)) == true

    s = simplex(e0,e1,e2)
    line1 = point(T,1/3,1/3,-1), point(T,1/3,1/3,+1)
    @test CompScienceMeshes.touches(s, line1) == true
    line1 = point(T,1,1,-1), point(T,1,1,+1)
    @test CompScienceMeshes.touches(s, line1) == false
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
@testset "Landscape" begin

    import ComputationalHomology: scale, value

    p = 1.0 => 2.0
    q = 0.0 => 3.0

    @test scale(q) == (1.5 => 1.5)
    @test value(p, q, 1) == 2.0

    bar = [p,q]
    @testset for (i,j) in zip(scale(bar), [1.5=>1.5, 1.5=>0.5])
        @test i == j
    end

    l = Landscape([p,q], [0.0=>1.0, 2.0=>Inf])
    @test length(l) == 2
    @test l[1][1] === p

    bar = [1=>5, 2=>8, 3=>4, 6=>7, 5=>9]
    L = Landscape(
        [ (-Inf=>0.0), (1=>0.0), (3.0=>2.0), (3.5=>1.5), (5.0=>3.0), (6.5=>1.5), (7.0=>2.0), (9=>0.0), (Inf=>0.0) ],
        [ (-Inf=>0.0), (2=>0.0), (3.5=>1.5), (5=>0.0), (6.5=>1.5), (8=>0.0), (Inf=>0.0) ],
        [ (-Inf=>0.0), (3=>0.0), (3.5=>0.5), (4=>0.0), (6=>0.0), (6.5=>0.5), (7=>0.0), (Inf=>0.0) ]
    )
    Lk = landscape(bar)
    @testset for (l1, l2) in zip(L, Lk)
        for (i1, i2) in zip(l1, l2)
            @test i1.first == i2.first
            @test i1.second == i2.second
        end
    end

    # reductions on landscapes
    l1 = Landscape([(-Inf=>0.0), (1.0=>0.0), (2.0=>1.0), (3.0=>0.0), (Inf=>0.0)])
    s = reduce(+, l1, 1.0)
    @testset for (p,q) in zip(l1[1], s[1])
        @test p.first == q.first
        @test p.second+1 == q.second
    end

    s = reduce(+, l1, l1)
    @test s[1][3] == (2 => 2)
    s = reduce(-, l1, l1)
    @test s[1][3] == (2 => 0)

    l2 = Landscape([(-Inf=>0.0), (2.0=>0.0), (3.0=>1.0), (4.0=>1.0), (5=>0.0), (Inf=>0.0)])
    s = reduce(+, l1, l2)
    @testset for (i, v) in zip(2:4, [(1 => 0), (2 => 1), (3 => 1)])
        @test s[1][i] == v
    end

    l3 = Landscape(l2[1], [(-Inf=>0.0), (1.0=>0.0), (2.0=>2.0), (3.0=>0.0), (Inf=>0.0)])
    s = reduce(+, l1, l3)
    @testset for (i, v) in zip(2:4, [(1 => 0), (2 => 1), (3 => 1)])
        @test s[1][i] == v
    end
    @test s[2][3] == (2 => 2)

    s = reduce(*, l3, l1)
    @testset for (i, v) in zip(2:4, [(1 => 0), (2 => 0), (3 => 0)])
        @test s[1][i] == v
    end
    @test s[2][3] == (2 => 0)

    m = mean([l1, l2, l3])
    @testset for (i, v) in zip(3:5, [(2 => 1/3), (3 => 2/3), (4 => 2/3)])
        @test m[1][i] == v
    end
    @test m[2][3] == (2 => 2/3)

end


text:
#Tests on the normal reaction networks.
for rn in reaction_networks_standard
    p_vals = map(amp->amp*rand(length(rn.params)),[1.,10.,100.])
    for p in p_vals
        fps = steady_states(rn,p)
        stab = map(fp->stability(fp,rn,p), fps)
        for i=1:length(fps)
            stab[i] ? (tend = 10) : (tend = -10)
            end_point = OrdinaryDiffEq.solve(ODEProblem(rn,fps[i],(0.,tend),p),Rosenbrock23())[end]
            @test maximum(abs.(end_point-fps[i]))<0.0001
        end
        #start_points = map(amp->amp*rand(length(rn.syms)),[1.,10.,100.])
        #for sp in start_points
        #    end_point = OrdinaryDiffEq.solve(ODEProblem(rn,sp,(0.,1000.),p),Rosenbrock23())[end]
        #    @test any(map(fp->maximum(abs.(fp-end_point)),fps) .< 0.001)
        #end
    end
end


text:
using KroneckerProducts
using LinearAlgebra
using Test

a = rand(ComplexF64, 3, 2)
b = rand(ComplexF64, 4, 5)

k_big = kron(a,b)

k = a ⊗ b

@test k[1] == a[1]*b[1]
@test k[end] == a[end]*b[end]
@test k[1,1] == a[1,1]*b[1,1]
@test k[end,end] == a[end,end]*b[end,end]

@test k == k_big

@test k' == k_big'
@test transpose(k) == transpose(k_big)

@test size(k) == size(k_big)

@test eltype(k) == eltype(k_big)

y = rand(ComplexF64, size(k,2))
yr = rand(ComplexF64, size(k,1))

@test k*y ≈ k_big*y
@test yr'*k ≈ yr'*k_big

a = rand(ComplexF64, 3, 2)
b = rand(ComplexF64, 2, 3)
k_big = kron(a,b)
k = KroneckerProduct(a,b)

@test det(k) == 0
@test abs(det(k_big)) < eps()

@test tr(k) ≈ tr(k_big)

a = rand(ComplexF64, 3, 3)
b = rand(ComplexF64, 2, 2)
k_big = kron(a,b)
k = KroneckerProduct(a,b)

@test det(k) ≈ det(k_big)

@test tr(k) ≈ tr(k_big)

yr = rand(ComplexF64, size(k,1))
@test k\yr ≈ k_big\yr


text:
export compute_relative_feasibility

"""
Compute transform-domain relative feasibility w.r.t. a set
r_feas =  ||P(A*x)-A*x||/||A*x||
"""
function compute_relative_feasibility(m::Vector{TF},
                                                feasibility_initial::Vector{TF},
                                                TD_OP::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},
                                                P_sub
                                                ) where {TF<:Real,TI<:Integer}

feasibility_initial[1]=norm(P_sub[1](TD_OP[1]*m) .- TD_OP[1]*m) ./ (norm(TD_OP[1]*m)+(100*eps(TF)))
#gc() check if we need this in julia >0.7
end


text:
function MCTS.ranked_actions(policy::RobustMCTSPlanner, state)
    actions = keys(policy.tree[state].A)
    q_val(a) = policy.tree[state].A[a].Q
    return sort!(collect(actions), by=q_val, rev=true)
end


text:
"""
    CeedScalar

Scalar (floating point) type. Equivalent to `Float64`.
"""
const CeedScalar = C.CeedScalar

"""
    CeedInt

Integer type, used for indexing. Equivalent to `Int32`.
"""
const CeedInt = C.CeedInt

"""
    QuadMode

One of `GAUSS` or `GAUSS_LOBATTO`.
"""
const QuadMode = C.CeedQuadMode
const GAUSS = C.CEED_GAUSS
const GAUSS_LOBATTO = C.CEED_GAUSS_LOBATTO

"""
    MemType

One of `MEM_HOST` or `MEM_DEVICE`.
"""
const MemType = C.CeedMemType
const MEM_HOST = C.CEED_MEM_HOST
const MEM_DEVICE = C.CEED_MEM_DEVICE

"""
    CopyMode

One of `COPY_VALUES`, `USE_POINTER` or `OWN_POINTER`.

`OWN_POINTER` is not typically supported for objects created in Julia, because those must be
destroyed by the garbage collector, and cannot be freed from C.
"""
const CopyMode = C.CeedCopyMode
const COPY_VALUES = C.CEED_COPY_VALUES
const USE_POINTER = C.CEED_USE_POINTER
const OWN_POINTER = C.CEED_OWN_POINTER

"""
    EvalMode

Evaluation mode used in the specification of input and output fields for Q-functions, e.g.
in [`@interior_qf`](@ref).

One of:
- `EVAL_NONE`
- `EVAL_INTERP`
- `EVAL_GRAD`
- `EVAL_DIV`
- `EVAL_CURL`
- `EVAL_WEIGHT`

See the [libCEED
documentation](https://libceed.readthedocs.io/en/latest/api/CeedBasis/?highlight=EVAL_MODE#c.CeedEvalMode)
for further information.
"""
const EvalMode = C.CeedEvalMode
const EVAL_NONE = C.CEED_EVAL_NONE
const EVAL_INTERP = C.CEED_EVAL_INTERP
const EVAL_GRAD = C.CEED_EVAL_GRAD
const EVAL_DIV = C.CEED_EVAL_DIV
const EVAL_CURL = C.CEED_EVAL_CURL
const EVAL_WEIGHT = C.CEED_EVAL_WEIGHT

"""
    NormType

Denotes type of vector norm. One of `NORM_1`, `NORM_2`, or `NORM_MAX`.
"""
const NormType = C.CeedNormType
const NORM_1 = C.CEED_NORM_1
const NORM_2 = C.CEED_NORM_2
const NORM_MAX = C.CEED_NORM_MAX

"""
    TransposeMose

Denotes whether a linear transformation or its transpose should be applied. Either
`NOTRANSPOSE` or `TRANSPOSE`.
"""
const TransposeMode = C.CeedTransposeMode
const NOTRANSPOSE = C.CEED_NOTRANSPOSE
const TRANSPOSE = C.CEED_TRANSPOSE

"""
    Topology

Type of basis shape to create non-tensor H1 element basis. One of `LINE`, `TRIANGLE`,
`QUAD`, `TET`, `PYRAMID`, `PRISM`, or `HEX`.

The dimension can be extracted with bitshift:

    dim = Int(topology) >> 16
"""
const Topology = C.CeedElemTopology
const LINE = C.CEED_LINE
const TRIANGLE = C.CEED_TRIANGLE
const QUAD = C.CEED_QUAD
const TET = C.CEED_TET
const PYRAMID = C.CEED_PYRAMID
const PRISM = C.CEED_PRISM
const HEX = C.CEED_HEX

function set_globals()
    @doc """
        STRIDES_BACKEND

    Indicate that the stride is determined by the backend.
    """
    global STRIDES_BACKEND = C.CEED_STRIDES_BACKEND[]
end

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
module Katz
using LightGraphs
using GraphMatrices
using IterativeSolvers
export katz, stdbasis, ikatz, ρ

function katz(A, a, v)
    F = lufact(I-a*A)
    return A*(F\v)
end

function ikatz(A::KatzAdjacency, b::AbstractVector, tol::Real=1e-8, maxiter::Integer=100)
    x, convhist = cg(A, b, tol=tol, maxiter=maxiter)
    return x, convhist
end

function ρ(convhist::ConvergenceHistory)
    rs = @show convhist.residuals
    return mean([rs[i+1]/rs[i] for i in 1:length(rs)-1])
end

function stdbasis(n,i)
    x = ones(n)
    x[i] + 1
    return x
end

end # module


text:
"""
The threading macro is like a more flexible version of the `|>` operator.

    @> x f = f(x)
    @> x g f == f(g(x))
    @> x a b c d e == e(d(c(b(a(x)))))

Unlike |>, functions can have arguments - the value
preceding a function will be treated as its first argument

    @> x g(y, z) f == f(g(x, y, z))

    @> x g f(y, z) == f(g(x), y, z)

See also `@>>`, `@as`.
"""
macro >(exs...)
  thread(x) = isexpr(x, :block) ? thread(rmlines(x).args...) : x

  thread(x, ex) =
    isexpr(ex, :call, :macrocall) ? Expr(ex.head, ex.args[1], x, ex.args[2:end]...) :
    isexpr(ex, :block)            ? thread(x, rmlines(ex).args...) :
    Expr(:call, ex, x)

  thread(x, exs...) = reduce(thread, x, exs)

  esc(thread(exs...))
end

"""
Same as `@>`, but threads the last argument.

  @>> x g(y, z) f == f(g(y, z, x))

  @>> x g f(y, z) == f(y, z, g(x))
"""
macro >>(exs...)
  thread(x) = isexpr(x, :block) ? thread(rmlines(x).args...) : x

  thread(x, ex) =
    isexpr(ex, :call, :macrocall) ? Expr(ex.head, ex.args..., x) :
    isexpr(ex, :block)            ? thread(x, rmlines(ex).args...) :
    Expr(:call, ex, x)

  thread(x, exs...) = reduce(thread, x, exs)

  esc(thread(exs...))
end


text:
#API for calling/creating services. Syntax is practically identical to rospy.
export Service, ServiceProxy, wait_for_service

using Compat

type ServiceProxy{SrvType <: ServiceDefinition}
    o::PyObject

    function ServiceProxy(name::AbstractString; kwargs...)
        @debug("Creating <$SrvType> service proxy for '$name'")
        rospycls = _get_rospy_class(SrvType)
        new(__rospy__[:ServiceProxy](ascii(name), rospycls; kwargs...))
    end
end
function ServiceProxy{SrvType<:ServiceDefinition}(
    name::AbstractString,
    srv::Type{SrvType};
    kwargs...
)
    ServiceProxy{SrvType}(ascii(name); kwargs...)
end

@compat function (srv::ServiceProxy{SrvType}){SrvType <: ServiceDefinition}(
    req::SrvT
)
    if ! isa(req, _srv_reqtype(SrvType))
        throw(ArgumentError(
            string("Incorrect service request type: ",typeof(req))))
    end
    pyresp = pycall(srv.o, PyObject, convert(PyObject, req))
    resp = convert(_srv_resptype(SrvType), pyresp)
    resp
end

if VERSION >= v"0.5.0-dev+3692" #callbacks are broken

type Service{T}
end

Service(args...) = error(
"""Providing a service is currently broken on julia v0.5 and above. See
https://github.com/jdlangs/RobotOS.jl/issues/15 for ongoing efforts to fix this.""")

else #callbacks not broken

type Service{SrvType <: ServiceDefinition}
    o::PyObject
    jl_handler

    function Service(name::AbstractString, handler; kwargs...)
        @debug("Providing <$SrvType> service at '$name'")
        rospycls = _get_rospy_class(SrvType)
        ReqType = _srv_reqtype(SrvType)
        jl_hndl(req::PyObject) =
            convert(PyObject, handler(convert(ReqType,req)))
        try
            new(__rospy__[:Service](ascii(name), rospycls, jl_hndl; kwargs...),
                jl_hndl
            )
        catch err
            if isa(err, PyCall.PyError)
                error("Problem during service creation: $(err.val[:args][1])")
            else
                rethrow(err)
            end
        end
    end
end
function Service{SrvType<:ServiceDefinition}(
    name::AbstractString,
    srv::Type{SrvType},
    handler;
    kwargs...
)
    Service{SrvType}(ascii(name), handler; kwargs...)
end

end #version check

function wait_for_service(service::AbstractString; kwargs...)
    try
        __rospy__[:wait_for_service](ascii(service); kwargs...)
    catch ex
        error("Timeout exceeded waiting for service '$service'")
    end
end


text:
@safetestset "ishexadecimal" begin
    using Validator

    @test Validator.ishexadecimal("deadBEEF") === true
    @test Validator.ishexadecimal("ff0044") === true
    @test Validator.ishexadecimal("0xff0044") === true
    @test Validator.ishexadecimal("0XfF0044") === true
    @test Validator.ishexadecimal("0x0123456789abcDEF") === true
    @test Validator.ishexadecimal("0X0123456789abcDEF") === true
    @test Validator.ishexadecimal("0hfedCBA9876543210") === true
    @test Validator.ishexadecimal("0HfedCBA9876543210") === true
    @test Validator.ishexadecimal("0123456789abcDEF") === true

    @test Validator.ishexadecimal("abcdefg") === false
    @test Validator.ishexadecimal("") === false
    @test Validator.ishexadecimal("..") === false
    @test Validator.ishexadecimal("0xa2h") === false
    @test Validator.ishexadecimal("0xa20x") === false
    @test Validator.ishexadecimal("0x0123456789abcDEFq  ") === false
    @test Validator.ishexadecimal("0hfedCBA9876543210q") === false
    @test Validator.ishexadecimal("01234q56789abcDEF") === false
end


text:
# Semantics for abstract values #

"""
    lub(a, b)

Compute the least upper bound (i.e. join) of abstract values `a` and `b`.
"""
function lub end

"""
    glb(a, b)

Compute the greatest lower bound (i.e. meet) of abstract values `a` and `b`.
"""
function glb end

"""
    widen(a, b)

Compute the widening of abstract values `a` and `b`.
"""
function widen end

include("boolean.jl")
include("interval.jl")

const DEFAULT_ABSTRACTIONS = Dict(
    :boolean => BooleanAbs,
    :integer => IntervalAbs{Int},
    :numeric => IntervalAbs{Float64}
)


text:
using DataDeps

register(DataDep(
  "hybrid-price",
  "http://www.stat.ufl.edu/~winner/data/hybrid_reg.txt",
  "http://www.stat.ufl.edu/~winner/data/hybrid_reg.csv",
  "4c081c159b942b6c3e7b99f55047c03374a7c6833010b32a3508398cbe0b63d9",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=4,
        feature_indices=5:7,
        id_indices=2,
        skipto=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using Test
using RobotDynamics
using StaticArrays
using Rotations
using BenchmarkTools

r = @SVector rand(3)
q = rand(UnitQuaternion)
v = @SVector rand(3)
ω = @SVector rand(3)

x = RBState{Float64}(r, q, v, ω)
@test x[1:3] ≈ r ≈ position(x)
@test x[4:7] ≈ Rotations.params(q)
@test x[8:10] ≈ v ≈ RD.linear_velocity(x)
@test x[11:13] ≈ ω ≈ RD.angular_velocity(x)
@test length(x) == 13

x32 = RBState{Float32}(r, q, v, ω)
@test x32 isa RBState{Float32}
@test eltype(x32[1:3]) == Float32
@test RD.position(x32) isa SVector{3,Float32}
@test RD.orientation(x32) isa UnitQuaternion{Float32}

# Pass in other types of vectors
x = RBState{Float64}(Vector(r), q, v, ω)
@test RD.position(x) isa SVector{3,Float64}
x = RBState{Float64}(Vector(r), q, Vector(v), Vector(ω))
@test RD.position(x) isa SVector{3,Float64}
@test RD.linear_velocity(x) isa SVector{3,Float64}
@test RD.angular_velocity(x) isa SVector{3,Float64}

x = RBState{Float64}(view(r,:), q, view(Vector(v),1:3), Vector(ω))
@test RD.position(x) isa SVector{3,Float64}

@test RD.RBState{Float32}(x) isa RBState{Float32}
@test RD.RBState{Float64}(x32) isa RBState{Float64}
@test RD.RBState(x) isa RBState{Float64}
@test RD.RBState(x32) isa RBState{Float32}

# Pass in other rotations
x2 = RD.RBState{Float64}(r, RotMatrix(q), v, ω)
@test RD.orientation(x2) \ RD.orientation(x) ≈ one(UnitQuaternion)
@test RD.orientation(x2) isa UnitQuaternion

# Let the constructor figure out data type
@test RBState(r, q, v, ω) isa RBState{Float64}
@test RBState(Float32.(r), UnitQuaternion{Float32}(q), Float32.(v), Float32.(ω)) isa RBState{Float32}
@test RBState(Float32.(r), q, Float32.(v), Float32.(ω)) isa RBState{Float64}
@test RBState(r, UnitQuaternion{Float32}(q), Float32.(v), Float32.(ω)) isa RBState{Float64}

# Convert from another rotation type
x = [r; Rotations.params(MRP(q)); v; ω]
@test RBState(MRP, x) ≈ RBState(r,q,v,ω)
x = [r; Rotations.params(RodriguesParam(q)); v; ω]
@test !(RBState(MRP, x) ≈ RBState(r,q,v,ω))
@test RBState(RodriguesParam, x) ≈ RBState(r,q,v,ω)

# Pass in a vector for the quaternion
q_ = Rotations.params(q)
x = RBState(r, q_, v, ω)
@test RD.orientation(x) ≈ q
x = RBState(r, 2q_, v, ω)
@test RD.orientation(x) ≈ 2q  # shouldnt renormalize
@test RBState(Float32.(r), Float32.(q_), Float32.(v), Float32.(ω)) isa RBState{Float32}

# Pass in a vector for the entire state
using RobotDynamics: position, orientation, linear_velocity, angular_velocity
x_ = [r; 2q_; v; ω]
x = RBState(x_)
@test position(x) ≈ r
@test orientation(x) ≈ 2q  # shouldnt renormalize
@test linear_velocity(x) ≈ v
@test angular_velocity(x) ≈ ω

@test RBState{Float32}(x_) isa RBState{Float32}
@test RBState(Vector(x_)) isa RBState{Float64}

# Test comparison (with double-cover)
q = rand(UnitQuaternion)
x1 = RBState(r, q, v, ω)
x2 = RBState(r, -q, v, ω)
@test x1[4:7] ≈ -x2[4:7]
@test x1 ≈ x2
@test !(SVector(x1) ≈ SVector(x2))

# Test indexing
x = RBState(r, q, v, ω)
@test x[1:3] ≈ r
@test x[4:7] ≈ Rotations.params(q)
@test x[8:10] ≈ v
@test x[11:13] ≈ ω
@test x[4] ≈ q.w

@test getindex(x,2) == r[2]
@test getindex(x,8) == v[1]
@test getindex(x,13) == ω[3]

# Renorm
x2 = RBState(r,2*q,v,ω)
@test norm(orientation(x2)) ≈ 2
x = RobotDynamics.renorm(x2)
@test norm(orientation(x2)) ≈ 2
@test norm(orientation(x)) ≈ 1
@test isrotation(orientation(x))

# Test zero constructor
x0 = zero(RBState)
@test position(x0) ≈ zero(r)
@test orientation(x0) ≈ one(UnitQuaternion)
@test linear_velocity(x0) ≈ zero(v)
@test angular_velocity(x0) ≈ zero(v)
@test zero(RBState) isa RBState{Float64}
@test zero(RBState{Float32}) isa RBState{Float32}
x = rand(RBState)
@test zero(x) ≈ x0

# Test random generator
x = rand(RBState)
@test isrotation(orientation(x))
@test x isa RBState{Float64}
@test position(x) isa SVector{3,Float64}
x32 = rand(RBState{Float32})
@test isrotation(orientation(x))
@test x32 isa RBState{Float32}
@test position(x32) isa SVector{3,Float32}

# Addition and Subtraction
x1 = rand(RBState)
x2 = rand(RBState)
x = x1 + x2
@test position(x) ≈ position(x1) + position(x2)
@test orientation(x) ≈ orientation(x1) * orientation(x2)
@test linear_velocity(x) ≈ linear_velocity(x1) + linear_velocity(x2)
@test angular_velocity(x) ≈ angular_velocity(x1) + angular_velocity(x2)
x = x1 - x2
@test position(x) ≈ position(x1) - position(x2)
@test orientation(x) ≈ orientation(x2) \ orientation(x1)
@test linear_velocity(x) ≈ linear_velocity(x1) - linear_velocity(x2)
@test angular_velocity(x) ≈ angular_velocity(x1) - angular_velocity(x2)

dx = x1 ⊖ x2
@test dx[1:3] ≈ position(x1) - position(x2)
@test dx[4:6] ≈ orientation(x1) ⊖ orientation(x2)
@test dx[7:9] ≈ linear_velocity(x1) - linear_velocity(x2)
@test dx[10:12] ≈ angular_velocity(x1) - angular_velocity(x2)

@test dx isa SVector{12}
@test x2 ⊕ dx ≈ x1

q1 = orientation(x1)
q2 = orientation(x2)
@test Rotations.params(RodriguesParam(q2) \ RodriguesParam(q1)) ≈ dx[4:6]

# Test randbetween
xmin = RBState(fill(-1,3), one(UnitQuaternion), fill(-2,3), fill(-3,3))
xmax = RBState(fill(+1,3), one(UnitQuaternion), fill(+2,3), fill(+3,3))
@test randbetween(xmin, xmax) isa RBState{Float64}
for k = 1:10
    local x = randbetween(xmin, xmax)
    @test all(position(xmin) .<= position(x) .<= position(xmax))
    @test all(linear_velocity(xmin) .<= linear_velocity(x) .<= linear_velocity(xmax))
    @test all(angular_velocity(xmin) .<= angular_velocity(x) .<= angular_velocity(xmax))
end


text:
using RecipesBase, VisualRegressionTests, Plots, Colors
include("utils.jl")
# Note that the `false` in the test avoids to start Gtk
# I (kellertuer)  think it also avoids asking whether a new reference should be
# created but I am not 100% sure

@testset "Recipes Test" begin
    references_folder = joinpath(@__DIR__, "assets")
    function Hyp2PB_plot()
        M = Hyperbolic(2)
        p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
        p2 = convert.(Ref(PoincareBallPoint), p)
        return plot(M, p2)
    end
    @plottest Hyp2PB_plot joinpath(references_folder, "Hyp2PBPlot.png") false

    function Hyp2PB_plot_geo()
        M = Hyperbolic(2)
        p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
        p2 = convert.(Ref(PoincareBallPoint), p)
        return plot(M, p2; geodesic_interpolation = 80)
    end
    @plottest Hyp2PB_plot_geo joinpath(references_folder, "Hyp2PBPlotGeo.png") false

    function Hyp2PB_quiver()
        M = Hyperbolic(2)
        p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 0.0], [0.0, 1.0]])
        p2 = convert.(Ref(PoincareBallPoint), p)
        X = [log(M, p2[2], p2[1]), log(M, p2[1], p2[3])]
        return plot(M, [p2[2], p2[1]], X)
    end
    @plottest Hyp2PB_quiver joinpath(references_folder, "Hyp2PBQuiver.png") false

    function Hyp2PH_plot()
        M = Hyperbolic(2)
        p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
        p2 = convert.(Ref(PoincareHalfSpacePoint), p)
        return plot(M, p2)
    end
    @plottest Hyp2PH_plot joinpath(references_folder, "Hyp2PHPlot.png") false

    function Hyp2PH_plot_geo()
        M = Hyperbolic(2)
        p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
        p2 = convert.(Ref(PoincareHalfSpacePoint), p)
        return plot(M, p2; geodesic_interpolation = 80)
    end
    @plottest Hyp2PH_plot_geo joinpath(references_folder, "Hyp2PHPlotGeo.png")

    function Hyp2PH_quiver()
        M = Hyperbolic(2)
        p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 0.0], [0.0, 1.0]])
        p2 = convert.(Ref(PoincareHalfSpacePoint), p)
        X = [log(M, p2[2], p2[1]), log(M, p2[1], p2[3])]
        return plot(M, [p2[2], p2[1]], X)
    end
    @plottest Hyp2PH_quiver joinpath(references_folder, "Hyp2PHQuiver.png") false
end


text:
mutable struct LoadingsScaleGibbsOperator <: OperatorXMLElement
    el::XMLOrNothing
    likelihoods::Vector{MyXMLElement}
    prior::MyXMLElement
    weight::Float64

    function LoadingsScaleGibbsOperator(factors::LatentFactorModelXMLElement,
                                        prior::MyXMLElement)
        return new(nothing, [factors], prior, 1.0)
    end

    function LoadingsScaleGibbsOperator(factors::IntegratedFactorsXMLElement,
                                    tree_likelihood::TraitLikelihoodXMLElement)
        return new(nothing, [factors, tree_likelihood], factors.loadings_prior, 1.0)
    end
end

function make_xml(lsgo::LoadingsScaleGibbsOperator)
    el = new_element(bn.LOCADINGS_SCALE_GIBBS)
    set_attribute(el, bn.WEIGHT, lsgo.weight)
    for like in lsgo.likelihoods
        add_ref_el(el, like)
    end
    add_ref_el(el, get_loadings_prior(lsgo.prior))
    lsgo.el = el
    return el
end

function get_parameter(lsg::LoadingsScaleGibbsOperator)
    for likelihood in lsg.likelihoods
        T = typeof(likelihood)
        if T <: IntegratedFactorsXMLElement || T <: LatentFactorModelXMLElement
            return get_loadings_scale(likelihood)
        end
    end

    error("Could not find parameter")
end


text:
## Array.jl ##

for (T, suff) in ((Float64, ""), (Float32, "f"))

    # 1-arg functions
    onearg_funcs = (:ceil,:floor,:sqrt,:rsqrt,:rec,
              :exp,:exp2,:expm1,:log,:log1p,:log2,:log10,
              :sin,:sinpi,:cos,:cospi,:tan,:tanpi,:asin,:acos,:atan,
              :sinh,:cosh,:tanh,:asinh,:acosh,:atanh)
    # Combine this with the list of renamed 1-arg funcs
    onearg_funcs = (
        ((x,x) for x in onearg_funcs)...,
        (:trunc,:int),
        (:round,:nint),
        (:exponent,:logb),
        (:abs,:fabs)
    )

    for (f,fa) in onearg_funcs
        f! = Symbol("$(f)!")
        @eval begin
            # Allocating variant
            function ($f)(X::Array{$T})
                out = Array{$T}(undef, size(X))
                ($f!)(out, X)
            end

            # In-place mutating variant
            function ($f!)(out::Array{$T}, X::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ref{Cint}),out,X,length(X))
                out
            end
        end
    end

    # 2 arg functions
    twoarg_funcs = (
        (:copysign,:copysign),
        (:rem,:fmod),
        (:div_float,:div),
        (:atan,:atan2)
    )

    for (f, fa) in twoarg_funcs
        f! = Symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = Array{$T}(undef, size(X))
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ref{Cint}),out,X,Y,length(X))
                out
            end
        end
    end

    # for some bizarre reason, vvpow/vvpowf reverse the order of arguments.
    for (f, fa) in ((:pow,:pow),)
        f! = Symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = Array{$T}(undef, size(X))
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ref{Cint}),out,Y,X,length(X))
                out
            end
        end
    end

    # two-arg return
    for (f, fa) in ((:sincos,:sincos),)
        f! = Symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T})
                out1 = Array{$T}(undef, size(X))
                out2 = Array{$T}(undef, size(X))
                ($f!)(out1, out2, X)
            end
            function ($f!)(out1::Array{$T}, out2::Array{$T}, X::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Cvoid,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ref{Cint}),out1,out2,X,length(X))
                out1, out2
            end
        end
    end

    # complex return
    for (f,fa) in ((:cis,:cosisin),)
        f! = Symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T})
                out = Array{Complex{$T}}(undef, size(X))
                ($f!)(out, X)
            end
            function ($f!)(out::Array{Complex{$T}}, X::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Cvoid,
                      (Ptr{Complex{$T}},Ptr{$T},Ref{Cint}),out,X,length(X))
                out
            end
        end
    end

    for (f, fa) in onearg_funcs
        f! = Symbol("$(f)!")
        @eval begin
            # Broadcasting override such that f.(X) turns into f(X)
            Base.copy(bc::Base.Broadcast.Broadcasted{Style, Axes, typeof($f), Tuple{Array{$T, N}}}) where {Style, Axes, N} = ($f)(bc.args...)
            Base.copyto!(dest::Array{$T, N}, bc::Base.Broadcast.Broadcasted{Style, Axes, typeof($f), Tuple{Array{$T, N}}}) where {Style, Axes, N} = ($f!)(dest, bc.args...)
        end
    end
    for (f, fa) in (twoarg_funcs...,(:pow,:pow))
        f! = Symbol("$(f)!")
        @eval begin
            # Broadcasting override such that f.(X) turns into f(X)
            Base.copy(bc::Base.Broadcast.Broadcasted{Style, Axes, typeof($f), Tuple{Array{$T, N},Array{$T,N}}}) where {Style, Axes, N} = ($f)(bc.args...)
            Base.copyto!(dest::Array{$T, N}, bc::Base.Broadcast.Broadcasted{Style, Axes, typeof($f), Tuple{Array{$T,N},Array{$T,N}}}) where {Style, Axes, N} = ($f!)(dest, bc.args...)
        end
    end
end

# Functions over single vectors that return scalars/tuples
for (T, suff) in ((Float32, ""), (Float64, "D"))

    for (f, fa) in ((:maximum, :maxv), (:minimum, :minv), (:mean, :meanv),
                    (:meansqr, :measqv), (:meanmag,  :meamgv), (:sum, :sve))
        @eval begin
            function ($f)(X::Vector{$T})
                val = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa, suff), libacc)),  Cvoid,
                      (Ptr{$T}, Int64,  Ref{$T}, UInt64),
                      X, 1, val, length(X))
                return val[]
            end
        end
    end

    for (f, fa) in ((:findmax, :maxvi), (:findmin, :minvi))
        @eval begin
            function ($f)(X::Vector{$T})
                index = Ref{Int}(0)
                val = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa, suff), libacc)),  Cvoid,
                      (Ptr{$T}, Int64,  Ref{$T}, Ref{Int}, UInt64),
                      X, 1, val, index, length(X))
                return (val[], index[]+1)
            end
        end
    end
end

# Element-wise operations over two vectors
for (T, suff) in ((Float32, ""), (Float64, "D"))

    for (f, name) in ((:vadd, "addition"),  (:vsub, "subtraction"),
                      (:vdiv, "division"), (:vmul, "multiplication"))
        f! = Symbol("$(f)!")

        @eval begin
            @doc """
            `$($f!)(result::Vector{$($T)}, X::Vector{$($T)}, Y::Vector{$($T)})`

            Implements element-wise **$($name)** over two **Vector{$($T)}** and overwrites
            the result vector with computed value. *Returns:* **Vector{$($T)}** `result`
            """ ->
            function ($f!)(result::Vector{$T}, X::Vector{$T}, Y::Vector{$T})
                ccall(($(string("vDSP_", f, suff), libacc)),  Cvoid,
                      (Ptr{$T}, Int64, Ptr{$T},  Int64, Ptr{$T}, Int64,  UInt64),
                      Y, 1, X, 1, result, 1, length(result))
                return result
            end
        end

        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)}, Y::Vector{$($T)})`

            Implements element-wise **$($name)** over two **Vector{$($T)}**. Allocates
            memory to store result. *Returns:* **Vector{$($T)}**
            """ ->
            function ($f)(X::Vector{$T}, Y::Vector{$T})
                result = similar(X)
                ($f!)(result, X, Y)
                return result
            end
        end
    end
end


text:
__precompile__(false)

module MeshManager

using MathManager
using LogManager
using GraphicsManager

using ModernGL

export Transform
export MeshBuffer
export MeshData
export update
export upload
export linkData
export setAttributes
export setDrawArray
export createBuffer
export createBuffers
export getMeshBufferRefID
export createArrayObjects

const GPU =  GraphicsManager
GPU.addListener(@__MODULE__)

"""
object which holds a transformation matrix and a model reference
"""
mutable struct Transform
  model::Mat4x4f
  mvp::Mat4x4f
end

"""
array which holds raw data of mesh
"""
mutable struct MeshBuffer
  elements::UInt32
  size::UInt32
  length::UInt32
  count::UInt32
  refID::GLuint
  typ::GLenum
  usage::GLenum
  loaded::Bool

  data::Union{Nothing,AbstractArray}

  function MeshBuffer(id::Symbol, typ=GL_ARRAY_BUFFER,usage=GL_STATIC_DRAW,data=nothing; refID=0, elems=1)
    this=new(elems,0,0,0,refID,typ,usage,false,data)
    setData(this, data)
    set(id,this)
    this
  end
end

"""
data which holds various arrays of mesh
"""
mutable struct MeshData
  id::Symbol
  vao::GLuint
  arrays::Dict{Symbol,MeshBuffer}
  draw::Union{Nothing,MeshBuffer}

  function MeshData()
    id = Symbol(length(LIST_DATA))
    this = new(id,0,Dict(),nothing)
    set(id,this)
    this
  end
end

LIST_BUFFER = Dict{Symbol, MeshBuffer}()
LIST_DATA = Dict{Symbol, MeshData}()

get(list::AbstractDict, id::Symbol) = haskey(list, id) ? list[id] : nothing

function get(typ::DataType, id::Symbol)
  list = nothing
  if typ == MeshBuffer list = LIST_BUFFER
  elseif typ == MeshData list = LIST_DATA
  end
  list != nothing ? get(list, id) : nothing
end

set(id::Symbol, buffer::MeshBuffer) = LIST_BUFFER[id] = buffer
set(id::Symbol, data::MeshData) = LIST_DATA[id] = data

function cleanUp()
  global LIST_BUFFER = typeof(LIST_BUFFER)()
  global LIST_DATA = typeof(LIST_DATA)()
end

function getMeshBufferRefID(buffers::Array{MeshBuffer,1})
  [buffer.refID for buffer in buffers]
end

gltypes=Dict(Float32=>GL_FLOAT,Float64=>GL_DOUBLE,UInt32=>GL_UNSIGNED_INT,Int32=>GL_INT)

"""
set attributes for shader
"""
function setAttributes(this::MeshBuffer, program, attrb; flexible=false, typ=0, bindbuffer=true)
  if this.count == 0 return end

  STRIDE = GLsizei(length(attrb)<=1 ? 0 : reduce(+, (x->sizeof(x[2])*x[3]).(attrb)))
  OFFSET =  C_NULL

  if bindbuffer glBindBuffer(typ<=0 ? this.typ : typ, this.refID) end
  #glBufferData(this.typ, sizeof(this.data), this.data,  this.usage)

  index=-1
  for (name,typ,elems,inst) in attrb
    index+=1
    atr = glGetAttribLocation(program, name)
    if flexible && atr <= -1 atr=index end
    if  atr > -1
      #typ = nothing; old = this.data
      #while true typ = eltype(old); if typ == old break; else old=typ; end; end
      glEnableVertexAttribArray(atr)
      glVertexAttribPointer(atr, elems, gltypes[typ], GL_FALSE, STRIDE, OFFSET)
      if inst>0 glVertexAttribDivisor(atr, inst) end
      OFFSET +=  sizeof(typ)*elems
    else warn("Could not load Attribute \"$name\"")
    end
  end

  if bindbuffer glBindBuffer(this.typ, 0) end
end

"""
creates vaos
"""
createArrayObjects(count=1) = GPU.create(:VERTEXARRAY, :DEFAULT, count)

"""
creates gpu buffers
"""
function createBuffers(this::MeshData)
  if this.vao == 0 this.vao = GPU.create(:VERTEXARRAY, this.id) end
  for (s,a) in this.arrays if a.refID == 0 && a.count > 0 a.refID = GPU.create(:BUFFER, this.id) end end
end

function createBuffers(id::Symbol, data::AbstractArray, count=1; size=0, typ=GL_ARRAY_BUFFER, usage=GL_STATIC_DRAW)
  list=Array{MeshBuffer,1}(undef, count)

  buffers=GPU.create(:BUFFER, id, count)

  has_data = length(data) > 0
  if size <= 0 size = sizeof(data) end

  for i=1:count
    refID=buffers[i]

    glBindBuffer(typ, refID)
    glBufferData(typ, size, has_data ? data : C_NULL,  usage)

    list[i] = MeshBuffer(Symbol(id,i), typ, usage, data; refID=refID,elems=1)
  end
  glBindBuffer(typ, 0)

  list
end

function createBuffer(id::Symbol, data::AbstractArray, count=1; typ=GL_ARRAY_BUFFER, usage=GL_STATIC_DRAW)
  if length(data) <= 0 return nothing end

  refID=GPU.create(:BUFFER, id)
  size=sizeof(data)

  this=MeshBuffer(id, typ, usage, nothing; refID=refID,elems=1)

  glBindBuffer(typ, refID)
  glBufferData(typ, count*size, C_NULL, usage)
  for i=1:count glBufferSubData(typ , (i-1)*size, size, data) end
  glBindBuffer(typ, 0)

  this
end

"""
delete gpu buffer
"""
function delete(this::MeshBuffer)
  GPU.delete(:BUFFER, this.refID)
  #GPU.glDeleteBuffer(this.refID)
end

"""
delete gpu buffer
"""
delete(buffers::Array{MeshBuffer,1}) = for buffer in buffers delete(buffer) end

export delete

"""
deletes gpu buffers
"""
function deleteBuffers(this::MeshData)
  if this.vao != 0 GPU.delete(:VERTEXARRAY, this.vao); this.vao=0 end
  for (s,a) in this.arrays if a.refID != 0 GPU.delete(:BUFFER, a.refID) end end
  this.arrays = Dict()
  this.draw = nothing
end

"""
sets attributes for shader
"""
function setAttributes(this::MeshData, program; flexible=false)
  glBindVertexArray(this.vao)
  glCheckError("glBindVertexArray")
  if haskey(this.arrays,:vertices) list=this.arrays[:vertices]; setAttributes(list, program, [("iVertex",Float32,list.elements,0)]) end
  if haskey(this.arrays,:points) list=this.arrays[:points]; setAttributes(list, program, [("iInstancePos",Float32,3,0),("iInstanceFlags",Float32,3,0)]) end
  if haskey(this.arrays,:instances) list=this.arrays[:instances]; setAttributes(list, program, [("iInstancePos",Float32,3,1),("iInstanceFlags",Float32,3,1)]) end
  glBindVertexArray(0)
end

"""
sets attributes for shader
"""
function setAttributes(this::MeshData, buffers::Array{MeshBuffer,1}, program)
  glBindVertexArray(this.vao)
  for buffer in buffers
    setAttributes(buffer, program, [("iInstancePos",Float32,3,0),("iInstanceFlags",Float32,3,0)])
  end
  glBindVertexArray(0)
end

"""
sets attributes for shader
"""
function setAttributes(vaos::Array{GLuint,1}, buffers::Array{MeshBuffer,1}, program)
  len=length(vaos)
  if len < length(buffers) warn("Object (VAO) count cannot be smaller than buffer count"); return end
  for i=1:len
    vao=vaos[i]
    buffer=buffers[i]
    glBindVertexArray(vao)
    setAttributes(buffer, program, [("iInstancePos",Float32,3,0),("iInstanceFlags",Float32,3,0)])
  end
  glBindVertexArray(0)
end

"""
sets array which will be drawn by opengl
"""
function setDrawArray(this::MeshData, key::Symbol)
  this.draw = this.arrays[key]
end

"""
sets data
"""
function setData(this::MeshBuffer, data, elems=0)
  isNotNull = data != nothing
  this.data = data
  this.size = isNotNull ? sizeof(this.data) : 0
  this.length = isNotNull ? length(this.data) : 0
  this.elements = elems
  this.count = isNotNull && elems > 0 ? this.length / this.elements : 0
end

"""
uploads data
"""
function upload(this::MeshBuffer, usage=this.usage)
  if this.loaded || this.count == 0 return end
  glBindBuffer(this.typ, this.refID)
  glBufferData(this.typ, this.size, this.data, usage)
  this.data = nothing
  this.loaded = true
end

"""
reset part data
"""
function resetPart(this::MeshBuffer, id=1, value=0)
  if id < 1 || id > this.length return error("Invalid index for data") end
  glBindBuffer(this.typ, this.refID)
  part=this.data[id]
  glBufferSubData(this.typ , 0, sizeof(part), typeof(part)[value])
  glBindBuffer(this.typ, 0)
end

"""
uploads data
"""
function upload(this::MeshData)
  glBindVertexArray(this.vao)
  for (s,a) in this.arrays upload(a) end
  glBindVertexArray(0)
end

"""
uploads data
"""
function upload(this::MeshData, key::Symbol, data::AbstractArray)
   if !haskey(this.arrays,key) warn("Could not find key $key"); return end
  a=this.arrays[key]
  setData(a, data)
  a.loaded = false
  glBindVertexArray(this.vao)
  upload(a)
  glBindVertexArray(0)
end

"""
links data
"""
function linkData(this::MeshData, args...)
  d=Dict(args)

  deleteBuffers(this)

  for (s,x) in d
    l=isa(x,Tuple) ? length(x) : 0

    data = l>0 ? x[1] : x
    dtyp = eltype(data)
    elems = l>1 ? x[2] : 1
    btyp = l>2 ? x[3] : GL_ARRAY_BUFFER
    usage = l>3 ? x[4] : GL_STATIC_DRAW
    draw = l>4 ? x[5] : false

    if haskey(this.arrays,s)
      a = this.arrays[s]
      a.loaded = false
    else
      this.arrays[s] = a = MeshBuffer(s)
      a.typ = btyp
      a.usage = usage
    end

    setData(a, data, elems)

    if draw this.draw = a
    elseif this.draw == nothing this.draw = a
    end
  end

  createBuffers(this)
  upload(this)
end

end #MeshManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
# Alphabet
# ========
#
# Alphabet of biological symbols.
#
# Subtypes of `Alphabet` represent a domain of biological characters. For
# example, `DNAAlphabet{2}` has a domain of unambiguous nucleotides (i.e. A, C,
# G, and T). These types are used for parameterizing biological sequences and so
# on. A pair of encoder and decoder is associated with an alphabet, which maps
# values between binary and Julia-level representation.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

"""
Alphabet of biological characters.
"""
abstract type Alphabet end

"""
DNA nucleotide alphabet.
"""
struct DNAAlphabet{n} <: Alphabet end

"""
RNA nucleotide alphabet.
"""
struct RNAAlphabet{n} <: Alphabet end

"""
Amino acid alphabet.
"""
struct AminoAcidAlphabet <: Alphabet end

"""
General character alphabet.
"""
struct CharAlphabet <: Alphabet end

"""
Void alphabet (internal use only).
"""
struct VoidAlphabet <: Alphabet end

const NucAlphs = Union{DNAAlphabet,RNAAlphabet}
const TwoBitNucs = Union{DNAAlphabet{2}, RNAAlphabet{2}}
const FourBitNucs = Union{DNAAlphabet{4}, RNAAlphabet{4}}

"""
The number of bits to represent the alphabet.
"""
function bitsof end

for n in (2, 4)
    @eval begin
        bitsof(::Type{DNAAlphabet{$n}}) = $n
        bitsof(::Type{RNAAlphabet{$n}}) = $n
    end
end
bitsof(::Type{AminoAcidAlphabet}) = 8
bitsof(::Type{CharAlphabet}) = 32
bitsof(::Type{VoidAlphabet}) = 0

Base.eltype(::Type{DNAAlphabet}) = DNA
Base.eltype(::Type{RNAAlphabet}) = RNA
Base.eltype(::Type{DNAAlphabet{n}}) where {n} = DNA
Base.eltype(::Type{RNAAlphabet{n}}) where {n} = RNA
Base.eltype(::Type{AminoAcidAlphabet}) = AminoAcid
Base.eltype(::Type{CharAlphabet}) = Char
Base.eltype(::Type{VoidAlphabet}) = Nothing

BioSymbols.alphabet(::Type{DNAAlphabet{2}}) = ACGT
BioSymbols.alphabet(::Type{RNAAlphabet{2}}) = ACGU
BioSymbols.alphabet(::Type{DNAAlphabet{4}}) = alphabet(DNA)
BioSymbols.alphabet(::Type{RNAAlphabet{4}}) = alphabet(RNA)
BioSymbols.alphabet(::Type{AminoAcidAlphabet}) = alphabet(AminoAcid)
# TODO: this alphabet includes invalid Unicode scalar values
BioSymbols.alphabet(::Type{CharAlphabet}) = typemin(Char):typemax(Char)
BioSymbols.alphabet(::Type{VoidAlphabet}) = nothing

# Promotion of Alphabets
# ----------------------

for alph in (DNAAlphabet, RNAAlphabet)
    @eval function Base.promote_rule(::Type{A}, ::Type{B}) where {A<:$alph,B<:$alph}
        return $alph{max(bitsof(A),bitsof(B))}
    end
end

# Encoders & Decoders
# -------------------

"""
Encode biological characters to binary representation.
"""
function encode end

struct EncodeError{A<:Alphabet,T} <: Exception
    val::T
end

EncodeError(::Type{A}, val::T) where {A,T} = EncodeError{A,T}(val)

function Base.showerror(io::IO, err::EncodeError{A}) where {A}
    print(io, "cannot encode ", err.val, " in ", A)
end

"""
Decode biological characters from binary representation.
"""
function decode end

struct DecodeError{A<:Alphabet,T} <: Exception
    val::T
end

DecodeError(::Type{A}, val::T) where {A,T} = DecodeError{A,T}(val)

function Base.showerror(io::IO, err::DecodeError{A}) where {A}
    print(io, "cannot decode ", err.val, " in ", A)
end

# DNA and RNA alphabets
# ---------------------

for A in (DNAAlphabet, RNAAlphabet)
    T = eltype(A)
    @eval begin
        # 2-bit encoding
        @inline function encode(::Type{$(A){2}}, nt::$(T))
            if count_ones(nt) != 1 || !isvalid(nt)
                throw(EncodeError($(A){2}, nt))
            end
            return convert(UInt8, trailing_zeros(nt))
        end
        @inline function decode(::Type{$(A){2}}, x::UInt8)
            if x > 0x03
                throw(DecodeError($(A){2}, x))
            end
            return reinterpret($(T), 0x01 << x)
        end
        @inline decode(::Type{$(A){2}}, x::Unsigned) = decode($(A){2}, UInt8(x))

        # 4-bit encoding
        @inline function encode(::Type{$(A){4}}, nt::$(T))
            if !isvalid(nt)
                throw(EncodeError($(A){4}, nt))
            end
            return reinterpret(UInt8, nt)
        end
        @inline function decode(::Type{$(A){4}}, x::UInt8)
            if !isvalid($(T), x)
                throw(DecodeError($(A){4}, x))
            end
            return reinterpret($(T), x)
        end
        @inline decode(::Type{$(A){4}}, x::Unsigned) = decode($(A){4}, UInt8(x))
    end
end

# AminoAcidAlphabet
# -----------------

@inline function encode(::Type{AminoAcidAlphabet}, aa::AminoAcid)
    if aa > AA_Gap
        throw(EncodeError(AminoAcidAlphabet, aa))
    end
    return reinterpret(UInt8, aa)
end

@inline function decode(::Type{AminoAcidAlphabet}, x::UInt8)
    if x > 0x1b
        throw(DecodeError(AminoAcidAlphabet, x))
    end
    return reinterpret(AminoAcid, x)
end

@inline function decode(::Type{AminoAcidAlphabet}, x::Unsigned)
    return decode(AminoAcidAlphabet, UInt8(x))
end

# CharAlphabet
# ------------

@inline function encode(::Type{CharAlphabet}, char::Char)
    if char > '\U10ffff'
        throw(EncodeError(CharAlphabet, char))
    end
    return reinterpret(UInt32, char)
end

@inline function decode(::Type{CharAlphabet}, x::UInt32)
    c = reinterpret(Char, x)
    if !isvalid(c)
        throw(DecodeError(CharAlphabet, x))
    end
    return c
end

@inline function decode(::Type{CharAlphabet}, x::Unsigned)
    return decode(CharAlphabet, UInt32(x))
end


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

mutable struct AliasType <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    paths::Any # spec type: Union{ Nothing, Vector{AliasPathType} } # spec name: paths

    function AliasType(;name=nothing, paths=nothing)
        o = new()
        validate_property(AliasType, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(AliasType, Symbol("paths"), paths)
        setfield!(o, Symbol("paths"), paths)
        o
    end
end # type AliasType

const _property_map_AliasType = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("paths")=>Symbol("paths"))
const _property_types_AliasType = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("paths")=>"Vector{AliasPathType}")
Base.propertynames(::Type{ AliasType }) = collect(keys(_property_map_AliasType))
Swagger.property_type(::Type{ AliasType }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_AliasType[name]))}
Swagger.field_name(::Type{ AliasType }, property_name::Symbol) =  _property_map_AliasType[property_name]

function check_required(o::AliasType)
    true
end

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