text:
@testset "cme.jl" begin
    @testset "estimate and statistic" begin
        kernel = (ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel()

        function deviation(testprediction, testtarget, prediction, target)
            testlocation = (testprediction, testtarget)
            return mapreduce(+, prediction, 1:length(prediction)) do p, t
                ((t == target) - p) * kernel(testlocation, (prediction, t))
            end
        end

        for nclasses in (2, 10, 100), nsamples in (10, 50, 100)
            # define estimator (sample test locations uniformly)
            dist = Dirichlet(nclasses, 1)
            testpredictions = [rand(dist) for _ in 1:(nsamples ÷ 10)]
            testtargets = rand(1:nclasses, nsamples ÷ 10)
            estimator = UCME(kernel, testpredictions, testtargets)

            # sample predictions and targets
            predictions = [rand(dist) for _ in 1:nsamples]
            targets_consistent = [
                rand(Categorical(prediction)) for prediction in predictions
            ]
            targets_onlyone = ones(Int, length(predictions))

            # compute calibration error estimate and test statistic
            for targets in (targets_consistent, targets_onlyone)
                test = AsymptoticCMETest(estimator, predictions, targets)

                @test test.kernel == kernel
                @test test.nsamples == nsamples
                @test test.ntestsamples == nsamples ÷ 10
                @test test.estimate ≈ estimator(predictions, targets)

                deviations =
                    deviation.(testpredictions', testtargets', predictions, targets)
                mean_deviations = vec(mean(deviations; dims=1))
                @test test.mean_deviations ≈ mean_deviations

                # use of `inv` can lead to slightly different results
                S = cov(deviations; dims=1)
                statistic = nsamples * mean_deviations' * inv(S) * mean_deviations
                @test test.statistic ≈ statistic rtol = 1e-4
            end
        end
    end

    @testset "consistency" begin
        kernel = (ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel()
        αs = 0.05:0.1:0.95
        nsamples = 100

        pvalues_consistent = Vector{Float64}(undef, 100)
        pvalues_onlyone = similar(pvalues_consistent)

        for nclasses in (2, 10)
            # create estimator (sample test locations uniformly)
            dist = Dirichlet(nclasses, 1)
            testpredictions = [rand(dist) for _ in 1:5]
            testtargets = rand(1:nclasses, 5)
            estimator = UCME(kernel, testpredictions, testtargets)

            predictions = [Vector{Float64}(undef, nclasses) for _ in 1:nsamples]
            targets_consistent = Vector{Int}(undef, nsamples)
            targets_onlyone = ones(Int, nsamples)

            for i in eachindex(pvalues_consistent)
                # sample predictions and targets
                for j in 1:nsamples
                    rand!(dist, predictions[j])
                    targets_consistent[j] = rand(Categorical(predictions[j]))
                end

                # define test
                test_consistent = AsymptoticCMETest(
                    estimator, predictions, targets_consistent
                )
                test_onlyone = AsymptoticCMETest(estimator, predictions, targets_onlyone)

                # estimate pvalues
                pvalues_consistent[i] = pvalue(test_consistent)
                pvalues_onlyone[i] = pvalue(test_onlyone)
            end

            # compute empirical test errors
            ecdf_consistent = ecdf(pvalues_consistent)
            @test all(ecdf_consistent(α) < α + 0.15 for α in αs)
            @test all(p < 0.05 for p in pvalues_onlyone)
        end
    end
end


text:
"""
    PoissonLikelihood(l=exp)

Poisson likelihood with rate defined as `l(f)`.

```math
    p(y|f) = \\operatorname{Poisson}(y | θ=l(f))
```

This is to be used if  we assume that the uncertainity associated
with the data follows a Poisson distribution.
"""
struct PoissonLikelihood{L<:AbstractLink} <: AbstractLikelihood
    invlink::L
end

PoissonLikelihood(l=exp) = PoissonLikelihood(link(l))

(l::PoissonLikelihood)(f::Real) = Poisson(l.invlink(f))

(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))

function expected_loglikelihood(
    ::AnalyticExpectation,
    ::PoissonLikelihood{ExpLink},
    q_f::AbstractVector{<:Normal},
    y::AbstractVector{<:Real},
)
    f_μ = mean.(q_f)
    return sum((y .* f_μ) - exp.(f_μ .+ (var.(q_f) / 2)) - loggamma.(y .+ 1))
end

default_expectation_method(::PoissonLikelihood{ExpLink}) = AnalyticExpectation()

ext: jl


text:
function loglikelihood!(
    gc::GaussianCopulaLMMObs{T},
    β::Vector{T},
    τ::T, # inverse of linear regression variance
    Σ::Matrix{T},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where T <: BlasReal
    n, p, q = size(gc.X, 1), size(gc.X, 2), size(gc.Z, 2)
    if needgrad
        fill!(gc.∇β, 0)
        fill!(gc.∇τ, 0)
        fill!(gc.∇Σ, 0)
    end
    if needhess
        fill!(gc.Hβ, 0)
        fill!(gc.Hτ, 0)
        fill!(gc.HΣ, 0)
    end
    # evaluate copula loglikelihood
    sqrtτ = sqrt(τ)
    update_res!(gc, β)
    standardize_res!(gc, sqrtτ)
    rss = abs2(norm(gc.res)) # RSS of standardized residual
    tr = (1//2)dot(gc.ztz, Σ)
    mul!(gc.storage_q1, transpose(gc.Z), gc.res) # storage_q1 = Z' * std residual
    mul!(gc.storage_q2, Σ, gc.storage_q1)        # storage_q2 = Σ * Z' * std residual
    qf = (1//2)dot(gc.storage_q1, gc.storage_q2)
    logl = - (n * log(2π) -  n * log(τ) + rss) / 2 - log(1 + tr) + log(1 + qf)
    # gradient
    if needgrad
        # wrt β
        mul!(gc.∇β, transpose(gc.X), gc.res)
        BLAS.gemv!('N', -inv(1 + qf), gc.xtz, gc.storage_q2, one(T), gc.∇β)
        gc.∇β .*= sqrtτ
        # wrt τ
        gc.∇τ[1] = (n - rss + 2qf / (1 + qf)) / 2τ
        # wrt Σ
        copyto!(gc.∇Σ, gc.ztz)
        BLAS.syrk!('U', 'N', (1//2)inv(1 + qf), gc.storage_q1, (-1//2)inv(1 + tr), gc.∇Σ)
        copytri!(gc.∇Σ, 'U')
    end
    # Hessian: TODO
    if needhess; end;
    # output
    logl
end

function fit!(
    gcm::GaussianCopulaLMMModel,
    solver=Ipopt.IpoptSolver(print_level=0)
    )
    p, q = size(gcm.data[1].X, 2), size(gcm.data[1].Z, 2)
    npar = p + 1 + (q * (q + 1)) >> 1
    optm = MathProgBase.NonlinearModel(solver)
    lb = fill(-Inf, npar)
    ub = fill( Inf, npar)
    MathProgBase.loadproblem!(optm, npar, 0, lb, ub, Float64[], Float64[], :Max, gcm)
    # starting point
    par0 = Vector{Float64}(undef, npar)
    modelpar_to_optimpar!(par0, gcm)
    MathProgBase.setwarmstart!(optm, par0)
    # optimize
    MathProgBase.optimize!(optm)
    optstat = MathProgBase.status(optm)
    optstat == :Optimal || @warn("Optimization unsuccesful; got $optstat")
    # refresh gradient and Hessian
    optimpar_to_modelpar!(gcm, MathProgBase.getsolution(optm))
    loglikelihood!(gcm, true, true)
    gcm
end

"""
    optimpar_to_modelpar!(gcm, par)

Translate optimization variables in `par` to the model parameters in `gcm`.
"""
function optimpar_to_modelpar!(
    gcm::GaussianCopulaLMMModel,
    par::Vector)
    p, q = size(gcm.data[1].X, 2), size(gcm.data[1].Z, 2)
    copyto!(gcm.β, 1, par, 1, p)
    gcm.τ[1] = exp(par[p+1])
    fill!(gcm.ΣL, 0)
    offset = p + 2
    for j in 1:q
        gcm.ΣL[j, j] = exp(par[offset])
        offset += 1
        for i in j+1:q
            gcm.ΣL[i, j] = par[offset]
            offset += 1
        end
    end
    mul!(gcm.Σ, gcm.ΣL, transpose(gcm.ΣL))
    nothing
end

"""
    modelpar_to_optimpar!(gcm, par)

Translate model parameters in `gcm` to optimization variables in `par`.
"""
function modelpar_to_optimpar!(
    par::Vector,
    gcm::GaussianCopulaLMMModel
    )
    p, q = size(gcm.data[1].X, 2), size(gcm.data[1].Z, 2)
    copyto!(par, gcm.β)
    par[p+1] = log(gcm.τ[1])
    Σchol = cholesky(Symmetric(gcm.Σ))
    gcm.ΣL .= Σchol.L
    offset = p + 2
    for j in 1:q
        par[offset] = log(gcm.ΣL[j, j])
        offset += 1
        for i in j+1:q
            par[offset] = gcm.ΣL[i, j]
            offset += 1
        end
    end
    par
end

function MathProgBase.initialize(
    gcm::GaussianCopulaLMMModel,
    requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad])
            error("Unsupported feature $feat")
        end
    end
end

MathProgBase.features_available(gcm::GaussianCopulaLMMModel) = [:Grad]

function MathProgBase.eval_f(
    gcm::GaussianCopulaLMMModel,
    par::Vector)
    optimpar_to_modelpar!(gcm, par)
    loglikelihood!(gcm, false, false)
end

function MathProgBase.eval_grad_f(
    gcm::GaussianCopulaLMMModel,
    grad::Vector,
    par::Vector)
    p, q = size(gcm.data[1].X, 2), size(gcm.data[1].Z, 2)
    optimpar_to_modelpar!(gcm, par)
    loglikelihood!(gcm, true, false)
    # gradient wrt β
    copyto!(grad, gcm.∇β)
    # gradient wrt log(τ)
    grad[p+1] = gcm.∇τ[1] * gcm.τ[1]
    # gradient wrt L
    mul!(gcm.storage_qq, gcm.∇Σ, gcm.ΣL)
    offset = p + 2
    for j in 1:q
        grad[offset] = 2gcm.storage_qq[j, j] * gcm.ΣL[j, j]
        offset += 1
        for i in j+1:q
            grad[offset] = 2gcm.storage_qq[i, j]
            offset += 1
        end
    end
    nothing
end


text:
"""
    Block

    rigid body with fully actuation
"""
struct Block end

function block(;
    representation=:maximal,
    timestep=0.05,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    side=0.5,
    contact=true,
    contact_type=:nonlinear,
    seed=1,
    vis=Visualizer(),
    info=nothing,
    name=:robot,
    control_scaling=Diagonal(ones(3)),
    opts_step=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    opts_grad=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    T=Float64)

    mechanism = get_block(
        timestep=timestep,
        gravity=gravity,
        friction_coefficient=friction_coefficient,
        side=side,
        contact=contact,
        contact_type=contact_type)

    initialize_block!(mechanism)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end
    nu = 3
    no = nx

    aspace = BoxSpace(nu,
        low=(-1.0e8 * ones(nu)),
        high=(1.0e8 * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)

    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z

    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = [I(nu) zeros(nu, nu)]

    build_robot(mechanism,
        vis=vis,
        name=name)

    TYPES = [Block, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    env = Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)

    return env
end

function visualize_block_force!(vis, anim, z_vis, u_vis;
    shift = [0.0; 0.25; 0.0],
    color=orange,
    )

    u_max = maximum([norm(u) for u in u_vis])
    force_vis = Dojo.ArrowVisualizer(vis[:force])
    setobject!(force_vis, Dojo.MeshPhongMaterial(color=color))

    for t = 1:length(z_vis)
        z = z_vis[t]
        u = (t == length(z_vis) ? 0.5 * u_vis[end] ./ u_max : 0.5 * u_vis[t] ./ u_max)
        MeshCat.atframe(anim, t) do
            settransform!(force_vis,
                Dojo.Point(z[1] - u[1] - shift[1], z[2] - u[2] - shift[2], z[3] - u[3] - shift[3]),
                Dojo.Vec(u[1], u[2], u[3]),
                shaft_radius=0.05,
                max_head_radius=0.1)
        end
    end
    MeshCat.setanimation!(vis, anim)

    return vis, anim
end


text:
# This file was generated, do not modify it. # hide
print(
    "Tuned model measurements on test:\n",
    "  brier loss: ", brier_loss(ŷ_tuned, ytest) |> mean, "\n",
    "  auc:        ", auc(ŷ_tuned, ytest),                "\n",
    "  accuracy:   ", accuracy(mode.(ŷ_tuned), ytest)
)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# control_point.jl
@doc """
### controlPoint

Genertates a uniformly distributed control points in the physical space.
Intended only for initializing of Mapping object as of this commit.

**Inputs**

*  `map`  : Object of Mapping type
*  `box`  : Object of BoundingBox type

"""->

function controlPoint(map::AbstractMappingType, box::AbstractBoundingBox)

  # Remember B-spline control points also lie on the interior of the volume
  nctl =sview(map.nctl, :)
  origin =sview(box.origin, :)
  S =sview(box.unitVector, :, 1)*box.lengths[1]
  T =sview(box.unitVector, :, 2)*box.lengths[2]
  U =sview(box.unitVector, :, 3)*box.lengths[3]

  for k = 0:map.nctl[3]-1
    for j = 0:map.nctl[2]-1
      for i = 0:map.nctl[1]-1
        map._cp_xyz[:,i+1,j+1,k+1] = box.origin + (i/(nctl[1]-1))*S +
                                    (j/(nctl[2]-1))*T + (k/(nctl[3]-1))*U
      end
    end
  end

  return nothing
end  # End function controlPoint

@doc """
### writeControlPointsVTS

Writes a `*.vts` file to be viewed in Paraview. Necessary for visualizing the
control points for testing

**Input**

* `map` : Object of AbstractMappingType

"""->

function writeControlPointsVTS(map::AbstractMappingType,
                               fname::AbstractString="control_points")

  vtsfile = vtk_grid(fname, real(map.cp_xyz))
  outfiles = vtk_save(vtsfile)

  return nothing
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

see description in JuMP.jl
=#

struct DOptimalDesignNative{T <: Real} <: ExampleInstanceNative{T}
    q::Int
    p::Int
    n::Int
    n_max::Int
    logdet_obj::Bool
    rootdet_obj::Bool
    geomean_obj::Bool
    use_logdet::Bool # use logdet cone, else use extended formulation
    use_rootdet::Bool # use rootdet cone, else use extended formulation
    use_epinorminf::Bool # use infinity norm cone, else use extended formulation
end

function build(inst::DOptimalDesignNative{T}) where {T <: Real}
    (q, p, n, n_max) = (inst.q, inst.p, inst.n, inst.n_max)
    @assert (p > q) && (n > q) && (n_max <= n)
    @assert +(inst.logdet_obj, inst.geomean_obj, inst.rootdet_obj) == 1
    V = randn(T, q, p)

    # upper bound and nonnegativity of numbers of trials
    if inst.use_epinorminf
        h_norminf = vcat(T(n_max) / 2, fill(-T(n_max) / 2, p))
        G_norminf = vcat(zeros(T, 1, p), -Matrix{T}(I, p, p))
        cones = Cones.Cone{T}[Cones.EpiNormInf{T, T}(p + 1)]
    else
        h_norminf = vcat(zeros(T, p), fill(T(n_max), p))
        G_norminf = vcat(Matrix{T}(-I, p, p), Matrix{T}(I, p, p))
        cones = Cones.Cone{T}[Cones.Nonnegative{T}(p), Cones.Nonnegative{T}(p)]
    end

    # constraint on total number of trials
    A = ones(T, 1, p)
    b = T[n]

    if (inst.logdet_obj && inst.use_logdet) ||
        (inst.rootdet_obj && inst.use_rootdet)
        # maximize the hypograph variable of the cone
        c = vcat(-one(T), zeros(T, p))

        # pad with hypograph variable
        A = hcat(zero(T), A)
        G_norminf = hcat(zeros(T, size(G_norminf, 1)), G_norminf)

        # dimension of vectorized matrix V*diag(np)*V'
        dimvec = Cones.svec_length(q)
        G_detcone = zeros(T, dimvec, p)
        l = 1
        for i in 1:q, j in 1:i
            for k in 1:p
                G_detcone[l, k] = -V[i, k] * V[j, k]
            end
            l += 1
        end
        Cones.scale_svec!(G_detcone, sqrt(T(2)))
        @assert l - 1 == dimvec

        if inst.logdet_obj
            push!(cones, Cones.HypoPerLogdetTri{T, T}(dimvec + 2))
            # pad with hypograph variable and perspective variable
            h_detcone = vcat(zero(T), one(T), zeros(T, dimvec))
            # include perspective variable
            G_detcone = [
                -one(T)    zeros(T, 1, p);
                zeros(T, 1, p + 1);
                zeros(T, dimvec)    G_detcone;
                ]
        else
            push!(cones, Cones.HypoRootdetTri{T, T}(dimvec + 1))
            # pad with hypograph variable
            h_detcone = zeros(T, dimvec + 1)
            # include perspective variable
            G_detcone = [
                -one(T)    zeros(T, 1, p);
                zeros(T, dimvec)    G_detcone;
                ]
        end

        # all conic constraints
        G = vcat(G_norminf, G_detcone)
        h = vcat(h_norminf, h_detcone)
    end

    if inst.geomean_obj
        # auxiliary matrix variable has pq elements stored row-major, auxiliary
        # lower tri variable has svec_length(q) elements, also stored row-major
        pq = p * q
        qq = q ^ 2
        num_trivars = Cones.svec_length(q)
        c = vcat(-one(T), zeros(T, p + pq + num_trivars))

        A_VW = zeros(T, qq, pq)
        A_lowertri = zeros(T, qq, num_trivars)
        # rows index (i, j)
        row_idx = 1
        for i in 1:q
            col_offset = sum(1:(i - 1)) + 1
            A_lowertri[row_idx:(row_idx + i - 1),
                col_offset:(col_offset + i - 1)] = -Matrix{T}(-I, i, i)
            for j in 1:q
                for k in 1:p
                    # columns index (k, j)
                    col_idx = (k - 1) * q + j
                    A_VW[row_idx, col_idx] = V[i, k]
                end
                row_idx += 1
            end
        end
        A = [
            zero(T)    A    zeros(T, 1, pq + num_trivars);
            zeros(T, qq, 1 + p)    A_VW    A_lowertri;
            ]
        b = vcat(b, zeros(T, qq))

        G_geo = zeros(T, q, num_trivars)
        for i in 1:q
            G_geo[i, sum(1:i)] = -1
        end
        push!(cones, Cones.HypoGeoMean{T}(1 + q))
        G_soc_epi = zeros(T, pq + p, p)
        G_soc = zeros(T, pq + p, pq)
        epi_idx = 1
        col_idx = 1
        for i in 1:p
            push!(cones, Cones.EpiNormEucl{T}(q + 1))
            G_soc_epi[epi_idx, i] = -sqrt(T(q))
            G_soc[(epi_idx + 1):(epi_idx + q), col_idx:(col_idx + q - 1)] =
                Matrix{T}(-I, q, q)
            epi_idx += q + 1
            col_idx += q
        end
        zero1 = zeros(T, size(G_norminf, 1), pq + num_trivars)
        zero2 = zeros(T, pq + p, num_trivars)
        G = [
            zeros(T, size(G_norminf, 1))    G_norminf    zero1;
            -one(T)    zeros(T, 1, p + pq + num_trivars); # geomean
            zeros(T, q, 1 + p + pq)    G_geo; # geomean
            zeros(T, pq + p)    G_soc_epi    G_soc    zero2; # epinormeucl
            ]
        h = vcat(h_norminf, zeros(T, p + 1 + q + pq))
    end

    if (inst.rootdet_obj && !inst.use_rootdet) ||
        (inst.logdet_obj && !inst.use_logdet)
        # extended formulations require an upper tri matrix of additional vars
        # we will store this matrix row-major
        num_trivars = Cones.svec_length(q)
        # vectorized dimension of extended psd matrix
        dimvec = q * (2 * q + 1)
        G_psd = zeros(T, dimvec, p + num_trivars)

        # index of diagonal elements in upper triangular matrix
        diag_idx(i::Int) = (i == 1 ? 1 : 1 + sum(q - j for j in 0:(i - 2)))

        # V*diag(np)*V
        l = 1
        for i in 1:q, j in 1:i
            for k in 1:p
                G_psd[l, k] = -V[i, k] * V[j, k]
            end
            l += 1
        end
        # [triangle' diag(triangle)]
        tri_idx = 1
        for i in 1:q
            # triangle'
            # skip zero-valued elements
            l += i - 1
            for j in i:q
                G_psd[l, p + tri_idx] = -1
                l += 1
                tri_idx += 1
            end
            # diag(triangle)
            # skip zero-valued elements
            l += i - 1
            G_psd[l, p + diag_idx(i)] = -1
            l += 1
        end
        Cones.scale_svec!(G_psd, sqrt(T(2)))

        h_psd = zeros(T, dimvec)
        push!(cones, Cones.PosSemidefTri{T, T}(dimvec))
    end

    if inst.rootdet_obj && !inst.use_rootdet
        c = vcat(zeros(T, p + num_trivars), -one(T))
        A = hcat(A, zeros(T, 1, num_trivars + 1))
        h_geo = zeros(T, q + 1)
        G_geo = zeros(T, q, num_trivars)
        for i in 1:q
            G_geo[i, diag_idx(i)] = -1
        end
        push!(cones, Cones.HypoGeoMean{T}(1 + q))
        # all conic constraints
        G = [
            G_norminf   zeros(T, size(G_norminf, 1), num_trivars + 1);
            G_psd    zeros(T, dimvec, 1); # psd
            zeros(T, 1, p)    zeros(T, 1, num_trivars)    -one(T); # hypogeomean
            zeros(T, q, p)    G_geo    zeros(T, q); # hypogeomean
            ]
        h = vcat(h_norminf, h_psd, h_geo)
    end

    if inst.logdet_obj && !inst.use_logdet
        # extended formulation for logdet
        # number of experiments, upper triangular matrix, hypograph variables
        dimx = p + num_trivars + q
        padx = num_trivars + q
        num_hypo = q
        # maximize the sum of hypograph variables of all hypoperlog cones
        c = vcat(zeros(T, p + num_trivars), -ones(T, num_hypo))

        G_log = zeros(T, 3 * q, dimx)
        h_log = zeros(T, 3 * q)
        offset = 1
        for i in 1:q
            # hypograph variable
            G_log[offset, p + num_trivars + i] = -1
            # perspective variable
            h_log[offset + 1] = 1
            # diagonal element in the triangular matrix
            G_log[offset + 2, p + diag_idx(i)] = -1
            push!(cones, Cones.HypoPerLog{T}(3))
            offset += 3
        end

        # pad with triangle matrix variables and q hypoperlog cone hypograph variables
        A = [A zeros(T, 1, padx)]
        # all conic constraints
        G = [
            G_norminf    zeros(T, size(G_norminf, 1), padx);
            G_psd    zeros(T, dimvec, num_hypo);
            G_log;
            ]
        h = vcat(h_norminf, h_psd, h_log)
    end

    model = Models.Model{T}(c, A, b, G, h, cones)
    return model
end


text:
const μ = 3.986e14 * 60^2
const r = 42164e3
const mc = 500.0
const n = sqrt(μ / r^3)

#K₁ = [-28.8287 0.1005 -1449.9754 0.0046;
#      -0.087 -33.2562 0.00462 -1451.5013]
const K₁11 = -28.8287; const K₁12 = 0.1005; const K₁13 = -1449.9754; const K₁14 = 0.0046;
const K₁21 = -0.087; const K₁22 = -33.2562; const K₁23 = 0.00462; const K₁24 = -1451.5013;

#K₂ = [-288.0288 0.1312 -9614.9898 0.0;
#      -0.1312 -288.0 0.0 -9614.9883]
const K₂11 = -288.0288; const K₂12 = 0.1312; const K₂13 = -9614.9898; const K₂14 = 0.0;
const K₂21 = -0.1312; const K₂22 = -288.0; const K₂23 = 0.0; const K₂24 = -9614.9883;

# dynamics in the 'approaching' mode
@taylorize function spacecraft_approaching!(t, x, dx)

    rc = sqrt((r + x[1])^2 + x[2]^2)
    uxy1 = (K₁11 * x[1] + K₁12 * x[2]) + (K₁13 * x[3] + K₁14 * x[4])
    uxy2 = (K₁21 * x[1] + K₁22 * x[2]) + (K₁23 * x[3] + K₁24 * x[4])

    # x' = vx
    dx[1] = x[3]

    # y' = vy
    dx[2] = x[4]

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    dx[3] = (n^2*x[1] + 2*(n*x[4])) + ((μ/(r^2) - μ/(rc^3)*(r + x[1])) + uxy1/mc)

    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    dx[4] = (n^2*x[2] - 2*(n*x[3])) - (μ/(rc^3)*x[2] - uxy2/mc)

    # t' = 1
    dx[5] = 1.0 + 0.0*x[1]

    return dx
end

# dynamics in the 'rendezvous attempt' mode
@taylorize function spacecraft_rendezvous_attempt!(t, x, dx)

    rc = sqrt((r + x[1])^2 + x[2]^2)
    uxy1 = (K₂11 * x[1] + K₂12 * x[2]) + (K₂13 * x[3] + K₂14 * x[4])
    uxy2 = (K₂21 * x[1] + K₂22 * x[2]) + (K₂23 * x[3] + K₂24 * x[4])

    # x' = vx
    dx[1] = x[3]

    # y' = vy
    dx[2] = x[4]

    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    dx[3] = (n^2*x[1] + 2*(n*x[4])) + ((μ/(r^2) - μ/(rc^3)*(r + x[1])) + uxy1/mc)

    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    dx[4] = (n^2*x[2] - 2*(n*x[3])) - (μ/(rc^3)*x[2] - uxy2/mc)

    # t' = 1
    dx[5] = 1.0 + 0.0*x[1]

    return dx
end

# dynamics in the 'aborting' mode
@taylorize function spacecraft_aborting!(t, x, dx)

    rc = sqrt((r + x[1])^2 + x[2]^2)

    # x' = vx
    dx[1] = x[3]

    # y' = vy
    dx[2] = x[4]

    # vx' = n²x + 2n*vy + μ/(r^2) + μ/(rc^3)*(r+x)
    dx[3] = (n^2*x[1] + 2*n*x[4]) + μ/(r^2) + μ/(rc^3)*(r + x[1])

    # vy' = n²y - 2n*vx - μ/(rc^3)y
    dx[4] = (n^2*x[2] - 2*n*x[3]) - μ/(rc^3)*x[2]

    # t' = 1
    dx[5] = 1.0 + 0.0*x[1]

    return dx
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
#= Finding the Nth number in the Fibonacci Sequence using Dynamic Programming
=#

## Function

function fibonacci(n)
    f = Int64[]
    push!(f, 0)
    push!(f, 1)
    for i = 3:n
        temp = f[i-1] + f[i-2]
        push!(f, temp)
    end
    return f[n]
end

## Input

print("Enter the number: ")
n = readline()
n = parse(Int64, n)
println("The $(n)th fibonacci number is ")

## Calling the function

fibonacci(n)

#=
Sample Test Case:
Input:
    Enter the number: 5
Output:
    The 5th fibonacci number is 3

Time Complexity: O ( N )
=#

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
# HashDict

import Base: KeyIterator, ValueIterator, haskey, get, getkey, delete!,
             pop!, empty!, filter!, setindex!, getindex, similar,
             sizehint, length, filter, isempty, start, next, done,
             keys, values, _tablesz, skip_deleted, serialize, deserialize, serialize_type

typealias Unordered Nothing
typealias Ordered   Int

type HashDict{K,V,O<:Union(Ordered,Unordered)} <: Associative{K,V}
    slots::Array{Uint8,1}
    keys::Array{K,1}
    vals::Array{V,1}
    idxs::Array{O,1}
    order::Array{O,1}
    ndel::Int
    count::Int
    deleter::Function

    function HashDict()
        n = 16
        new(zeros(Uint8,n), Array(K,n), Array(V,n), Array(O,n), Array(O,0), 0, 0, identity)
    end
    function HashDict(ks, vs)
        n = length(ks)
        h = HashDict{K,V,O}()
        for i=1:n
            h[ks[i]] = vs[i]
        end
        return h
    end
    function HashDict(kv)
        h = HashDict{K,V,O}()
        for (k,v) in kv
            h[k] = v
        end
        return h
    end
end

HashDict() = HashDict{Any,Any,Unordered}()

HashDict{K,V}(ks::AbstractArray{K}, vs::AbstractArray{V}) = HashDict{K,V,Unordered}(ks,vs)
HashDict(ks, vs) = HashDict{Any,Any,Unordered}(ks, vs)
HashDict{K,V}(kv::AbstractArray{(K,V)}) = HashDict{K,V,Unordered}(kv)

# TODO: these could be more efficient
HashDict{K,V,O}(d::HashDict{K,V,O}) = HashDict{K,V,O}(collect(kv))
HashDict{K,V}(d::Associative{K,V}) = HashDict{K,V,Unordered}(collect(d))

similar{K,V,O}(d::HashDict{K,V,O}) = HashDict{K,V,O}()

function serialize(s, t::HashDict)
    serialize_type(s, typeof(t))
    write(s, int32(length(t)))
    for (k,v) in t
        serialize(s, k)
        serialize(s, v)
    end
end

function deserialize{K,V,O}(s, T::Type{HashDict{K,V,O}})
    n = read(s, Int32)
    t = T(); sizehint(t, n)
    for i = 1:n
        k = deserialize(s)
        v = deserialize(s)
        t[k] = v
    end
    return t
end

hashindex(key, sz) = (int(hash(key)) & (sz-1)) + 1

isslotempty(h::HashDict, i::Int) = h.slots[i] == 0x0
isslotfilled(h::HashDict, i::Int) = h.slots[i] == 0x1
isslotmissing(h::HashDict, i::Int) = h.slots[i] == 0x2

function rehash{K,V}(h::HashDict{K,V,Unordered}, newsz)
    newsz = _tablesz(newsz)

    if h.count == 0
        resize!(h.slots, newsz)
        fill!(h.slots, 0)
        resize!(h.keys, newsz)
        resize!(h.vals, newsz)
        h.ndel = 0
        return h
    end

    olds = h.slots
    oldk = h.keys
    oldv = h.vals
    sz = length(olds)

    slots = zeros(Uint8,newsz)
    keys = Array(K, newsz)
    vals = Array(V, newsz)
    count0 = h.count
    count = 0

    for i = 1:sz
        if olds[i] == 0x1
            k = oldk[i]
            v = oldv[i]
            index = hashindex(k, newsz)
            while slots[index] != 0
                index = (index & (newsz-1)) + 1
            end
            slots[index] = 0x1
            keys[index] = k
            vals[index] = v
            count += 1

            if h.count != count0
                # if items are removed by finalizers, retry
                return rehash(h, newsz)
            end
        end
    end

    h.slots = slots
    h.keys = keys
    h.vals = vals
    h.count = count
    h.ndel = 0

    return h
end

function rehash{K,V}(h::HashDict{K,V,Ordered}, newsz)
    newsz = _tablesz(newsz)

    if h.count == 0
        resize!(h.slots, newsz)
        fill!(h.slots, 0)
        resize!(h.keys, newsz)
        resize!(h.vals, newsz)
        resize!(h.idxs, newsz)
        resize!(h.order, 0)
        h.ndel = 0
        return h
    end

    _compact_order(h)

    olds = h.slots
    oldk = h.keys
    oldv = h.vals
    oldi = h.idxs
    oldo = h.order
    sz = length(olds)

    slots = zeros(Uint8,newsz)
    keys = Array(K, newsz)
    vals = Array(V, newsz)
    idxs = Array(Int, newsz)
    order = Array(Int, h.count)
    count0 = h.count
    count = 0

    for i = 1:sz
        if olds[i] == 0x1
            k = oldk[i]
            v = oldv[i]
            idx = oldi[i]
            index = hashindex(k, newsz)
            while slots[index] != 0
                index = (index & (newsz-1)) + 1
            end
            slots[index] = 0x1
            keys[index] = k
            vals[index] = v
            idxs[index] = idx
            order[idx] = index
            count += 1

            if h.count != count0
                # if items are removed by finalizers, retry
                return rehash(h, newsz)
            end
        end
    end

    h.slots = slots
    h.keys = keys
    h.vals = vals
    h.idxs = idxs
    h.order = order
    h.count = count
    h.ndel = 0

    return h
end

function _compact_order{K,V}(h::HashDict{K,V,Ordered})
    if h.count == length(h.order)
        return
    end

    i = 1
    while h.order[i] > 0;  i += 1; end

    j = i+1
    while h.order[j] == 0; j += 1; end

    for k = j:length(h.order)
        idx = h.order[k]
        if idx > 0
            h.order[i] = idx
            h.idxs[idx] = i
            i += 1
        end
    end

    resize!(h.order, h.count)

    nothing
end

function sizehint(d::HashDict, newsz)
    oldsz = length(d.slots)
    if newsz <= oldsz
        # todo: shrink
        # be careful: rehash() assumes everything fits. it was only designed
        # for growing.
        return d
    end
    # grow at least 25%
    newsz = max(newsz, (oldsz*5)>>2)
    rehash(d, newsz)
end

function empty!{K,V}(h::HashDict{K,V})
    fill!(h.slots, 0x0)
    sz = length(h.slots)
    h.keys = Array(K, sz)
    h.vals = Array(V, sz)
    h.ndel = 0
    h.count = 0
    return h
end

function empty!{K,V}(h::HashDict{K,V,Ordered})
    sz = length(h.slots)
    fill!(h.slots, 0x0)
    h.keys = Array(K, sz)
    h.vals = Array(V, sz)
    h.idxs = Array(Int, sz)
    h.order = Array(Int, 0)
    h.ndel = 0
    h.count = 0
    return h
end

# get the index where a key is stored, or -1 if not present
function ht_keyindex{K,V}(h::HashDict{K,V}, key)
    sz = length(h.keys)
    iter = 0
    maxprobe = max(16, sz>>6)
    index = hashindex(key, sz)
    keys = h.keys

    while true
        if isslotempty(h,index)
            break
        end
        if !isslotmissing(h,index) && isequal(key,keys[index])
            return index
        end

        index = (index & (sz-1)) + 1
        iter+=1
        iter > maxprobe && break
    end

    return -1
end

# get the index where a key is stored, or -pos if not present
# and the key would be inserted at pos
# This version is for use by setindex! and get!
function ht_keyindex2{K,V}(h::HashDict{K,V}, key)
    sz = length(h.keys)

    if h.ndel >= ((3*sz)>>2) || h.count*3 > sz*2
        # > 3/4 deleted or > 2/3 full
        rehash(h, h.count > 64000 ? h.count*2 : h.count*4)
        sz = length(h.keys)  # rehash may resize the table at this point!
    end

    iter = 0
    maxprobe = max(16, sz>>6)
    index = hashindex(key, sz)
    avail = 0
    keys = h.keys

    while true
        if isslotempty(h,index)
            avail < 0 && return avail
            return -index
        end

        if isslotmissing(h,index)
            if avail == 0
                # found an available slot, but need to keep scanning
                # in case "key" already exists in a later collided slot.
                avail = -index
            end
        elseif isequal(key, keys[index])
            return index
        end

        index = (index & (sz-1)) + 1
        iter+=1
        iter > maxprobe && break
    end

    avail < 0 && return avail

    rehash(h, h.count > 64000 ? sz*2 : sz*4)

    return ht_keyindex2(h, key)
end

function _setindex!(h::HashDict, v, key, index)
    h.slots[index] = 0x1
    h.keys[index] = key
    h.vals[index] = v
    h.count += 1
    return h
end

function _setindex!{K,V}(h::HashDict{K,V,Ordered}, v, key, index)
    h.slots[index] = 0x1
    h.keys[index] = key
    h.vals[index] = v
    push!(h.order, index)
    h.idxs[index] = length(h.order)
    h.count += 1
    return h
end

function setindex!{K,V}(h::HashDict{K,V}, v0, key0)
    key = convert(K,key0)
    if !isequal(key,key0)
        error(key0, " is not a valid key for type ", K)
    end
    v   = convert(V,  v0)

    index = ht_keyindex2(h, key)

    if index > 0
        h.vals[index] = v
    else
        _setindex!(h, v, key, -index)
    end

    return h
end

function get!{K,V}(h::HashDict{K,V}, key0, default)
    key = convert(K,key0)
    if !isequal(key,key0)
        error(key0, " is not a valid key for type ", K)
    end

    index = ht_keyindex2(h, key)

    index > 0 && return h.vals[index]

    v = convert(V,  default)
    _setindex!(h, v, key, -index)
    return v
end

function get!{K,V}(h::HashDict{K,V}, key0, default)
    key = convert(K,key0)
    if !isequal(key,key0)
        error(key0, " is not a valid key for type ", K)
    end

    index = ht_keyindex2(h, key)

    index > 0 && return h.vals[index]

    v = convert(V,  default)
    _setindex!(h, v, key, -index)
    return v
end

# TODO: this makes it challenging to have V<:Base.Callable
function get!{K,V,F<:Base.Callable}(h::HashDict{K,V}, key0, default::F)
    key = convert(K,key0)
    if !isequal(key,key0)
        error(key0, " is not a valid key for type ", K)
    end

    index = ht_keyindex2(h, key)

    index > 0 && return h.vals[index]

    v = convert(V,  default())
    _setindex!(h, v, key, -index)
    return v
end

function getindex{K,V}(h::HashDict{K,V}, key)
    index = ht_keyindex(h, key)
    return (index<0) ? throw(KeyError(key)) : h.vals[index]::V
end

function get{K,V}(h::HashDict{K,V}, key, deflt)
    index = ht_keyindex(h, key)
    return (index<0) ? deflt : h.vals[index]::V
end

haskey(h::HashDict, key) = (ht_keyindex(h, key) >= 0)
contains{T<:HashDict}(v::KeyIterator{T}, key) = (ht_keyindex(v.dict, key) >= 0)

function getkey{K,V}(h::HashDict{K,V}, key, deflt)
    index = ht_keyindex(h, key)
    return (index<0) ? deflt : h.keys[index]::K
end

function _pop!(h::HashDict, index)
    val = h.vals[index]
    _delete!(h, index)
    return val
end

function pop!(h::HashDict, key)
    index = ht_keyindex(h, key)
    index > 0 ? _pop!(h, index) : throw(KeyError(key))
end

function pop!(h::HashDict, key, default)
    index = ht_keyindex(h, key)
    index > 0 ? _pop!(h, index) : default
end

function _delete!(h::HashDict, index)
    h.slots[index] = 0x2
    ccall(:jl_arrayunset, Void, (Any, Uint), h.keys, index-1)
    ccall(:jl_arrayunset, Void, (Any, Uint), h.vals, index-1)
    h.ndel += 1
    h.count -= 1
    return h
end

function _delete!{K,V}(h::HashDict{K,V,Ordered}, index)
    h.slots[index] = 0x2
    ccall(:jl_arrayunset, Void, (Any, Uint), h.keys, index-1)
    ccall(:jl_arrayunset, Void, (Any, Uint), h.vals, index-1)
    h.order[h.idxs[index]] = 0
    h.ndel += 1
    h.count -= 1
    return h
end

function delete!(h::HashDict, key)
    index = ht_keyindex(h, key)
    index > 0 && _delete!(h, index)
    return h
end

function skip_deleted{K,V,O}(h::HashDict{K,V,O}, i)
    L = length(h.slots)
    while i<=L && !isslotfilled(h,i)
        i += 1
    end
    return i
end

function skip_deleted{K,V}(h::HashDict{K,V,Ordered}, i)
    L = length(h.order)
    while i<=L && h.order[i] == 0
        i += 1
    end
    return i
end

start(t::HashDict) = skip_deleted(t, 1)
done(t::HashDict, i) = done(t.vals, i)
next(t::HashDict, i) = ((t.keys[i],t.vals[i]), skip_deleted(t,i+1))

done{K,V}(t::HashDict{K,V,Ordered}, i) = done(t.order, i)
next{K,V}(t::HashDict{K,V,Ordered}, i) = ((t.keys[t.order[i]],t.vals[t.order[i]]), skip_deleted(t,i+1))

isempty(t::HashDict) = (t.count == 0)
length(t::HashDict) = t.count

next(v::KeyIterator{HashDict}, i) = (v.dict.keys[i], skip_deleted(v.dict,i+1))
next(v::ValueIterator{HashDict}, i) = (v.dict.vals[i], skip_deleted(v.dict,i+1))

next{K,V}(v::KeyIterator{HashDict{K,V,Ordered}}, i) = (v.dict.keys[v.dict.order[i]], skip_deleted(v.dict,i+1))
next{K,V}(v::ValueIterator{HashDict{K,V,Ordered}}, i) = (v.dict.vals[v.dict.order[i]], skip_deleted(v.dict,i+1))


text:
function activate_github(reponame; tag = nothing, sha = nothing, force = false, add_default_environment = false)
    # Make sure that our .projects environment is kosher.
    projdir = joinpath(pwd(), ".projects")
    mkpath(projdir)
    # For each case of inputs, end up with a concrete URL to download.
    if sha != nothing
        length(sha) >= 6 || throw(ArgumentError("SHA needs to be at least 6 characters."))
        tag == nothing || throw(ArgumentError("You can't give both a tag and an SHA hash."))
        tarprefix = sha
    elseif tag != nothing # Given a version but no SHA.
        tarprefix = tag
    else # Default to master if nothing is given.
        tarprefix = "master"
    end
    # Common objects for all cases.
    repostr = split(reponame, "/")[2]
    target = joinpath(projdir, "$repostr-$tarprefix")
    # Branches of logic.
    if isdir(target) && tarprefix != "master" && force == false # Static prefix that's already downloaded, no force.
        Pkg.activate(target)
    else
        # Make a temporary directory for our use.
        tmpdir = tempdir()
        # Download the tarball to that directory.
        tarurl = "https://github.com/$(reponame)/archive/$(tarprefix).tar.gz"
        tarpath = joinpath(tmpdir, "$repostr-$tarprefix.tar.gz")
        printstyled("Downloading ", bold=true, color=:light_green); println("$reponame-$tarprefix → $projdir")
        try
            run(gen_download_cmd(tarurl, tarpath)) # Download the tarball.
        catch e
            if e isa MethodError
                printstyled("Package installation and activation not supported in this setup. You should add packages manually.", bold = true, color = :red)
                return nothing
            else
                throw(e)
            end
        end
        # Unpack the tarball to that directory.
        run(gen_unpack_cmd(tarpath, tmpdir))
        # Remove the tarball to avoid path conflict with the next steps.
        rm(tarpath)
        # Find the path of the unpacked tarball (could be a full SHA)
        sourcedir = filter(object -> occursin("$repostr", object), readdir(tmpdir))[1] # There will only be one of these.
        # Move to .projects
        mv(joinpath(tmpdir, sourcedir), target, force = true) # Force will overwrite existing dir.
        # Clean.
        isdir(joinpath(tmpdir, sourcedir)) == false || rm(joinpath(tmpdir, sourcedir), recursive = true) # Important for logic.
        # Instantiate and precompile.
        printstyled("Instantiating ", bold=true, color=:light_green); println(target)
        Pkg.activate(target)
        pkg"instantiate"
        pkg"precompile"
    end
    projpath = joinpath(target, "Project.toml")
    packages = parsefile(projpath)["deps"]
    if add_default_environment # seed the default environment with the new packages if true
        printstyled("Adding to the default environment... ", bold=true, color=:light_green);
        pkg"activate "
        @show pkglist = Array{String}(collect(keys(packages)))
        Pkg.add(pkglist)
        Pkg.activate(target) # go back to the activated environment
    end
    tarprefix, target, packages
end

function copy_env(reponame, oldtag, newtag)
    repostr = split(reponame, "/")[2]
    projdir = joinpath(pwd(), ".projects")
    olddir = joinpath(projdir, "$repostr-$oldtag")
    newdir = joinpath(projdir, "$repostr-$newtag")
    cp(olddir, newdir, force = true)
end

# Clone the TOML files from some git repo's _current state_ to the local dir
function activate_github_path(reponame; path = "",
                                tag = "master",
                                force = false,
                                activate = true)
    # download step
    if "Project.toml" ∈ readdir(pwd()) && force == false
        @warn "There's already a Project.toml in the current directory, and force = false."
    else
        # url construction and download
        url_project = (path == "") ? join(["https://raw.githubusercontent.com", reponame, tag, "Project.toml"], "/") : join(["https://raw.githubusercontent.com", reponame, tag, path, "Project.toml"], "/")
        url_manifest = (path == "") ? join(["https://raw.githubusercontent.com", reponame, tag, "Manifest.toml"], "/") : join(["https://raw.githubusercontent.com", reponame, tag, path, "Manifest.toml"], "/")
        try # project
            resp_project = HTTP.get(url_project);
            io = open("Project.toml", "w")
            println(io, String(resp_project.body))
            close(io)
        catch e
            @warn "There was an error retrieving the Project.toml"
            throw(e)
        end
        try # manifest
            resp_manifest = HTTP.get(url_manifest);
            io = open("Manifest.toml", "w")
            println(io, String(resp_manifest.body))
            close(io)
        catch e
            # do nothing, since we aren't too fussed about the Manifest.
        end
    end

    # activation step
    need_activation = (Base.active_project() != joinpath(pwd(), "Project.toml"))
    if activate && need_activation
        pkg"activate ."
        pkg"instantiate"
        pkg"precompile"
    end
end


text:
using BinaryTraitsTestX
using Test

@testset "BinaryTraitsTestX.jl" begin
    # Write your own tests here.
end


text:
using Test, Tables

@testset "utils.jl" begin

    NT = NamedTuple{(), Tuple{}}
    @test Tables.names(NT) === ()
    @test Tables.types(NT) === Tuple{}
    @test isempty(Tables.runlength(Tables.types(NT)))
    @test Tables.columnindex(Tables.names(NT), :i) == 0
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :i) == Union{}

    NT = NamedTuple{(:a, :b, :c), NTuple{3, Int64}}
    @test Tables.names(NT) === (:a, :b, :c)
    @test Tables.types(NT) === Tuple{Int64, Int64, Int64}
    @test Tables.runlength(Tables.types(NT)) == [(Int64, 3)]
    @test Tables.columnindex(Tables.names(NT), :a) == 1
    @test Tables.columnindex(Tables.names(NT), :i) == 0
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :a) == Int64
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :i) == Union{}

    NT = NamedTuple{Tuple(Symbol("a$i") for i = 1:20), Tuple{vcat(fill(Int, 10), fill(String, 10))...}}
    @test Tables.names(NT) === Tuple(Symbol("a$i") for i = 1:20)
    @test Tables.types(NT) === Tuple{vcat(fill(Int, 10), fill(String, 10))...}
    @test Tables.runlength(Tables.types(NT)) == [(Int, 10), (String, 10)]
    @test Tables.columnindex(Tables.names(NT), :a20) == 20
    @test Tables.columnindex(Tables.names(NT), :i) == 0
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :a20) == String
    @test Tables.columntype(Tables.names(NT), Tables.types(NT), :i) == Union{}

    nt = (a=1, b=2, c=3)
    @test getproperty(nt, Int, 1, :a) === 1

    NT = typeof(nt)
    output = [0, 0, 0]
    Tables.eachcolumn(Tables.Schema(Tables.names(NT), Tables.types(NT)), nt, output) do val, col, nm, out
        out[col] = val
    end
    @test output == [1, 2, 3]

    nt = NamedTuple{Tuple(Symbol("a$i") for i = 1:101)}(Tuple(i for i = 1:101))
    NT =typeof(nt)
    @test Tables.runlength(Tables.types(NT)) == [(Int, 101)]
    output = zeros(Int, 101)
    Tables.eachcolumn(Tables.Schema(Tables.names(NT), Tables.types(NT)), nt, output) do val, col, nm, out
        out[col] = val
    end
    @test output == [i for i = 1:101]

    nt = NamedTuple{Tuple(Symbol("a$i") for i = 1:101)}(Tuple(i % 2 == 0 ? i : "$i" for i = 1:101))
    NT = typeof(nt)
    @test Tables.runlength(Tables.types(NT)) == [i % 2 == 0 ? (Int, 1) : (String, 1) for i = 1:101]
    output = Vector{Any}(undef, 101)
    Tables.eachcolumn(Tables.Schema(Tables.names(NT), Tables.types(NT)), nt, output) do val, col, nm, out
        out[col] = val
    end
    @test output == [i % 2 == 0 ? i : "$i" for i = 1:101]

    nt = (a=Ref(0), b=Ref(0))
    Tables.eachcolumn(Tables.Schema((:a, :b), nothing), (a=1, b=2)) do val, col, nm
        nt[nm][] = val
    end
    @test nt.a[] == 1
    @test nt.b[] == 2

    nt = (a=[1,2,3], b=[4,5,6])
    @test collect(Tables.eachcolumn(nt)) == [[1,2,3], [4,5,6]]
end

@testset "namedtuples.jl" begin

    nt = (a=1, b=2, c=3)
    rt = [nt, nt, nt]
    @test Tables.rows(rt) === rt
    @test Tables.schema(rt).names == Tables.names(typeof(nt))
    @test Tables.namedtupleiterator(eltype(rt), rt) === rt

    rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    nt = (a=[1,2,3], b=[4.0, 5.0, 6.0], c=["7", "8", "9"])
    @test Tables.rowcount(nt) == 3
    @test Tables.schema(nt) == Tables.Schema((:a, :b, :c), Tuple{Int, Float64, String})
    @test Tables.columnaccess(typeof(nt))
    @test Tables.columns(nt) === nt
    @test rowtable(nt) == rt
    @test columntable(rt) == nt
    @test rt == (rt |> columntable |> rowtable)
    @test nt == (nt |> rowtable |> columntable)

    @test Tables.buildcolumns(nothing, rt) == nt
    @test Tables.columntable(nothing, nt) == nt

    # append
    nt2 = columntable(nt, rt)
    @test Tables.rowcount(nt2) == 6
    @test Tables.schema(nt2) == Tables.Schema((:a, :b, :c), Tuple{Int, Float64, String})
    @test nt2 == (a = [1, 2, 3, 1, 2, 3], b = [4.0, 5.0, 6.0, 4.0, 5.0, 6.0], c = ["7", "8", "9", "7", "8", "9"])
    rt2 = rowtable(rt, nt)
    @test length(rt2) == 9

    rt = [(a=1, b=4.0, c="7"), (a=2.0, b=missing, c="8"), (a=3, b=6.0, c="9")]
    tt = Tables.buildcolumns(nothing, rt)
    @test isequal(tt, (a = [1.0, 2.0, 3.0], b = Union{Missing, Float64}[4.0, missing, 6.0], c = ["7", "8", "9"]))
    @test tt.a[1] === 1.0
    @test tt.a[2] === 2.0
    @test tt.a[3] === 3.0

    nti = Tables.NamedTupleIterator{Nothing, typeof(rt)}(rt)
    nti2 = collect(nti)
    @test isequal(rt, nti2)
end

@testset "Materializer" begin
    rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    nt = (a=[1,2,3], b=[4.0, 5.0, 6.0], c=["7", "8", "9"])

    @test nt == Tables.materializer(nt)(Tables.columns(nt))
    @test nt == Tables.materializer(nt)(Tables.columns(rt))
    @test nt == Tables.materializer(nt)(rt)
    @test rt == Tables.materializer(rt)(nt)

    function select(table, cols::Symbol...)
        Tables.istable(table) || throw(ArgumentError("select requires a table input"))
        nt = Tables.columntable(table)  # columntable(t) creates a NamedTuple of AbstractVectors
        newcols = NamedTuple{cols}(nt)
        Tables.materializer(table)(newcols)
    end

    @test select(nt, :a, :b, :c) == nt
    @test select(nt, :c, :a) == NamedTuple{(:c, :a)}(nt)
    @test select(rt, :a) == [(a=1,), (a=2,), (a=3,)]
end

import Base: ==
struct GenericRow
    a::Int
    b::Float64
    c::String
end
==(a::GenericRow, b::GenericRow) = a.a == b.a && a.b == b.b && a.c == b.c

struct GenericRowTable
    data::Vector{GenericRow}
end
==(a::GenericRowTable, b::GenericRowTable) = all(a.data .== b.data)

Base.eltype(g::GenericRowTable) = GenericRow
Base.length(g::GenericRowTable) = length(g.data)
Base.size(g::GenericRowTable) = (length(g.data),)
Tables.istable(::Type{GenericRowTable}) = true
Tables.rowaccess(::Type{GenericRowTable}) = true
Tables.rows(x::GenericRowTable) = x
Tables.schema(x::GenericRowTable) = Tables.Schema((:a, :b, :c), Tuple{Int, Float64, String})

function Base.iterate(g::GenericRowTable, st=1)
    st > length(g.data) && return nothing
    return g.data[st], st + 1
end

genericrowtable(x) = GenericRowTable(collect(map(x->GenericRow(x.a, x.b, x.c), Tables.rows(x))))

struct GenericColumn{T} <: AbstractVector{T}
    data::Vector{T}
end
Base.eltype(g::GenericColumn{T}) where {T} = T
Base.length(g::GenericColumn) = length(g.data)
==(a::GenericColumn, b::GenericColumn) = a.data == b.data
Base.getindex(g::GenericColumn, i::Int) = g.data[i]

struct GenericColumnTable
    names::Dict{Symbol, Int}
    data::Vector{GenericColumn}
end

Tables.istable(::Type{GenericColumnTable}) = true
Tables.columnaccess(::Type{GenericColumnTable}) = true
Tables.columns(x::GenericColumnTable) = x
Tables.schema(g::GenericColumnTable) = Tables.Schema(Tuple(keys(getfield(g, 1))), Tuple{(eltype(x) for x in getfield(g, 2))...})
Base.getproperty(g::GenericColumnTable, nm::Symbol) = getfield(g, 2)[getfield(g, 1)[nm]]
Base.propertynames(g::GenericColumnTable) = Tuple(keys(getfield(g, 1)))

function genericcolumntable(x)
    cols = Tables.columns(x)
    sch = Tables.schema(x)
    data = [GenericColumn(getproperty(cols, nm)) for nm in sch.names]
    return GenericColumnTable(Dict(nm=>i for (i, nm) in enumerate(sch.names)), data)
end
==(a::GenericColumnTable, b::GenericColumnTable) = getfield(a, 1) == getfield(b, 1) && getfield(a, 2) == getfield(b, 2)

@testset "Tables.jl interface" begin

    @test !Tables.istable(1)
    @test !Tables.istable(Int)
    @test !Tables.rowaccess(1)
    @test !Tables.rowaccess(Int)
    @test !Tables.columnaccess(1)
    @test !Tables.columnaccess(Int)
    @test Tables.schema(1) === nothing

    sch = Tables.Schema{(:a, :b), Tuple{Int64, Float64}}()
    @test Tables.Schema((:a, :b), Tuple{Int64, Float64}) === sch
    @test Tables.Schema(NamedTuple{(:a, :b), Tuple{Int64, Float64}}) === sch
    @test Tables.Schema((:a, :b), nothing) === Tables.Schema{(:a, :b), nothing}()
    @test Tables.Schema([:a, :b], [Int64, Float64]) === sch
    show(sch)
    @test sch.names == (:a, :b)
    @test sch.types == (Int64, Float64)
    @test_throws ArgumentError sch.foobar

    gr = GenericRowTable([GenericRow(1, 4.0, "7"), GenericRow(2, 5.0, "8"), GenericRow(3, 6.0, "9")])
    gc = GenericColumnTable(Dict(:a=>1, :b=>2, :c=>3), [GenericColumn([1,2,3]), GenericColumn([4.0, 5.0, 6.0]), GenericColumn(["7", "8", "9"])])
    @test gc == (gr |> genericcolumntable)
    @test gr == (gc |> genericrowtable)
    @test gr == (gr |> genericrowtable)
end

@static if :Query in Symbol.(Base.loaded_modules_array())
    rt = (a = Real[1, 2.0, 3], b = Union{Missing, Float64}[4.0, missing, 6.0], c = ["7", "8", "9"])

    dv = Tables.datavaluerows(rt)
    @test eltype(dv) == NamedTuple{(:a, :b, :c),Tuple{Real,DataValue{Float64},String}}
    rt2 = collect(dv)
    @test rt2[1] == (a = 1, b = DataValue{Float64}(4.0), c = "7")

    ei = Tables.DataValueUnwrapper(QueryOperators.EnumerableIterable{eltype(dv), typeof(dv)}(dv))
    nt = ei |> columntable
    @test isequal(rt, nt)
    rt3 = ei |> rowtable
    @test isequal(rt |> rowtable, rt3)

    # rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    mt = Tables.DataValueUnwrapper(ei.x |> y->QueryOperators.map(y, x->(a=x.b, c=x.c), Expr(:block)))
    @inferred (mt |> columntable)
    @inferred (mt |> rowtable)

    # uninferrable case
    mt = Tables.DataValueUnwrapper(ei.x |> y->QueryOperators.map(y, x->(a=x.a, c=x.c), Expr(:block)))
    @test (mt |> columntable) == (a = Real[1, 2.0, 3], c = ["7", "8", "9"])
    @test length(mt |> rowtable) == 3

    rt = (a = Missing[missing, missing], b=[1,2])
    dv = Tables.datavaluerows(rt)
    @test eltype(dv) == NamedTuple{(:a, :b), Tuple{DataValue{Union{}}, Int}}
end


text: include("laplacian.jl")


text:
using Base: @deprecate, @deprecate_binding, depwarn

if !isdefined(Base, :stderr)
    @deprecate stderr(obj::StatisticalModel) stderror(obj)
else
    function (io::typeof(stderr))(obj::StatisticalModel)
        Base.depwarn("stderr(obj::StatisticalModel) is deprecated, use stderror(obj) instead", :stderr)
        io === stderr ? stderror(obj) : throw(MethodErrror(io, (obj,)))
    end
end

@deprecate model_response(obj::StatisticalModel) response(obj)

@deprecate norepeats(a::AbstractArray) allunique(a)

@deprecate(mad!(v::AbstractArray{<:Real}, center;
                constant::Real = BigFloat("1.482602218505601860547076529360423431326703202590312896536266275245674447622701")),
           mad!(v, center=center, constant=constant))

### Deprecated January 2019
@deprecate scattermatm(x::DenseMatrix, mean, dims::Int) scattermat(x, mean=mean, dims=dims)
@deprecate scattermatm(x::DenseMatrix, mean, wv::AbstractWeights, dims::Int) scattermat(x, wv, mean=mean, dims=dims)
@deprecate scattermat(x::DenseMatrix, dims::Int) scattermat(x, dims=dims)
@deprecate scattermat(x::DenseMatrix, wv::AbstractWeights, dims::Int) scattermat(x, wv, dims=dims)
@deprecate scattermat_zm(x::DenseMatrix, dims::Int) scattermat_zm(x, dims=dims)
@deprecate scattermat_zm(x::DenseMatrix, wv::AbstractWeights, dims::Int) scattermat_zm(x::DenseMatrix, wv::AbstractWeights, dims=dims)
@deprecate mean!(R::AbstractArray, A::AbstractArray, w::AbstractWeights, dims::Int) mean!(R, A, w, dims=dims)
@deprecate mean(A::AbstractArray{T}, w::AbstractWeights{W}, dims::Int) where {T<:Number,W<:Real} mean(A, w, dims=dims)

@deprecate wquantile(v::RealVector, w::AbstractWeights{<:Real}, p::RealVector) quantile(v, w, p)
@deprecate wquantile(v::RealVector, w::AbstractWeights{<:Real}, p::Number) quantile(v, w, [p])[1]
@deprecate wquantile(v::RealVector, w::RealVector, p::RealVector) quantile(v, pweights(w), p)
@deprecate wquantile(v::RealVector, w::RealVector, p::Number) quantile(v, pweights(w), [p])[1]
@deprecate wmedian(v::RealVector, w::AbstractWeights{<:Real}) median(v, w)
@deprecate wmedian(v::RealVector, w::RealVector) median(v, weights(w))

@deprecate quantile(v::AbstractArray{<:Real}) quantile(v, [.0, .25, .5, .75, 1.0])

### Deprecated September 2019
@deprecate sum(A::AbstractArray, w::AbstractWeights, dims::Int) sum(A, w, dims=dims)
@deprecate values(wv::AbstractWeights) convert(Vector, wv)

### Deprecated November 2021
@deprecate stdm(x::RealArray, w::AbstractWeights, m::Real; corrected::DepBool=nothing) std(x, w, mean=m, corrected=corrected) false
@deprecate varm(x::RealArray, w::AbstractWeights, m::Real; corrected::DepBool=nothing) var(x, w, mean=m, corrected=corrected) false
@deprecate stdm(x::RealArray, w::AbstractWeights, m::RealArray, dim::Int; corrected::DepBool=nothing) std(x, w, dim, mean=m, corrected=corrected) false
@deprecate varm(x::RealArray, w::AbstractWeights, m::RealArray, dim::Int; corrected::DepBool=nothing) var(x, w, dim, mean=m, corrected=corrected) false
@deprecate varm!(R::AbstractArray, x::RealArray, w::AbstractWeights, m::RealArray, dim::Int; corrected::DepBool=nothing) var!(R, x, w, dim, mean=m, corrected=corrected) false


text:
"""
    ck45()

Returns coefficients rka,rkb,rkc for the 4th order 5-stage low storage Carpenter/Kennedy
Runge Kutta method. Coefficients evolve the residual, solution, and local time, e.g.,

# Example
```julia
res = rk4a[i]*res + dt*rhs # i = RK stage
@. u += rk4b[i]*res
```
"""
function ck45()
    rk4a = [            0.0 ...
    -567301805773.0/1357537059087.0 ...
    -2404267990393.0/2016746695238.0 ...
    -3550918686646.0/2091501179385.0  ...
    -1275806237668.0/842570457699.0];

    rk4b = [ 1432997174477.0/9575080441755.0 ...
    5161836677717.0/13612068292357.0 ...
    1720146321549.0/2090206949498.0  ...
    3134564353537.0/4481467310338.0  ...
    2277821191437.0/14882151754819.0]

    rk4c = [ 0.0  ...
    1432997174477.0/9575080441755.0 ...
    2526269341429.0/6820363962896.0 ...
    2006345519317.0/3224310063776.0 ...
    2802321613138.0/2924317926251.0 ...
    1.0];
    return rk4a, rk4b, rk4c
end


text:
function initialize(::Type{MarkovChain},
                    mcmc::MCMC{T},
                    progress_channel::RemoteChannel;
                    attempts::Int64=1000) where T <: EpidemicModel
  if attempts <= 0
    error("Must have at least 1 initialization attempt")
  end
  bested = 0
  max_lposterior = -Inf
  local markov_chain
  for i in 1:attempts
    events = generate(Events, mcmc)
    rparams = generate(RiskParameters, mcmc)
    @debug "Risk function parameters for attempt $i" θ = convert(Vector{Float64}, rparams)
    lprior = logpriors(rparams, mcmc.risk_priors)
    llikelihood, tr, tx = loglikelihood(rparams,
                                        mcmc.risk_functions,
                                        events,
                                        mcmc.population,
                                        mcmc.starting_states,
                                        early_decision_value = max_lposterior - lprior)
    lposterior = llikelihood + lprior
    if lposterior > max_lposterior
      network = generate(TransmissionNetwork, tr, mcmc, tx)
      markov_chain = MarkovChain(events, network, rparams, lposterior)
      max_lposterior = lposterior
      bested += 1
    end
    put!(progress_channel, true)
  end
  if max_lposterior == -Inf
    error("Failed to initialize Markov Chain")
  end
  @debug "Initialization improved upon $bested times"
  return markov_chain
end

function initialize(::Type{MarkovChain},
                    mcmc::MCMC{T};
                    attempts::Int64=1000) where T <: EpidemicModel
  if attempts <= 0
    error("Must have at least 1 initialization attempt")
  end
  max_lposterior = -Inf
  local markov_chain
  bested = 0
  pmeter = Progress(attempts, "Initialization progress")
  for i in 1:attempts
    next!(pmeter)
    events = generate(Events, mcmc)
    rparams = generate(RiskParameters, mcmc)
    @debug "Risk function parameters for attempt $i" θ = convert(Vector{Float64}, rparams)
    lprior = logpriors(rparams, mcmc.risk_priors)
    llikelihood, tr, tx = loglikelihood(rparams,
                                            mcmc.risk_functions,
                                            events,
                                            mcmc.population,
                                            mcmc.starting_states,
                                            early_decision_value = max_lposterior - lprior)
    lposterior = llikelihood + lprior
    if lposterior > max_lposterior
      network = generate(TransmissionNetwork, tr, mcmc, tx)
      markov_chain = MarkovChain(events, network, rparams, lposterior)
      max_lposterior = lposterior
      bested += 1
    end
  end
  finish!(pmeter)
  if max_lposterior == -Inf
    error("Failed to initialize Markov Chain")
    markov_chain = nothing
  end
  @debug "Initialization improved upon $bested times"
  return markov_chain
end


text:
using Test
using QBase

@testset "./src/evolve.jl" begin

@testset "*(unitary, ket)" begin
    ket = Ket([1.,0], atol=1e-6)
    evo_ket = σx*ket
    @test evo_ket isa Ket{Float64}
    @test evo_ket == [0,1]
    @test evo_ket.atol == 1e-6

    evo_ket2 = σz*σy*Ket([1,0])
    @test evo_ket2 isa Ket{Complex{Int64}}
    @test evo_ket2 == [0,-1im]
    @test evo_ket2.atol == QBase.ATOL
end

@testset "*(bra, unitary)" begin
    bra = Bra([1.,0], atol=1e-6)
    rx = qubit_rotation(π/2)
    evo_bra = bra*rx
    @test evo_bra isa Bra{Complex{Float64}}
    @test evo_bra ≈ [1 -1im]/sqrt(2)
    @test evo_bra.atol == 1e-6

    bra = Bra([1.,0], atol=1e-6)
    evo_bra = bra*rx'
    @test evo_bra isa Bra{Complex{Float64}}
    @test evo_bra ≈[1 1im]/sqrt(2)
end

@testset "*(unitary, state,  unitary)"  begin
    M = σx * State([1. 0;0 0]) * σz
    @test M == [0 0;1 0]
    @test M isa Matrix{Float64}

    M = σy * State([1 0;0 0]) * σy'
    @test M ==  [0 0;0 1]
    @test M isa Matrix{Complex{Int64}}
end

@testset "evolve()" begin
    @testset "qubit evolution" begin
        ρ = evolve(σx,State([1 0;0 0]))

        @test ρ == [0 0;0 1]
        @test ρ isa State{Int64}
    end

    @testset "ket evolution" begin
        ψ = evolve(σz, Ket([1, 1]/sqrt(2)))

        @test ψ isa Ket{Float64}
        @test ψ == [1,-1]/sqrt(2)
    end
end

end


text:
using QuantumOptics
using PyPlot

b = SpinBasis(1//2)

sx = sigmax(b)
sy = sigmay(b)
sz = sigmaz(b)

N = 6

spinchain = b ⊗ b

for ispin in 3:N
    global spinchain
    spinchain = spinchain ⊗ b
end

XX = sx ⊗ sx
YY = sy ⊗ sy
ZZ = sz ⊗ sz

J=1

HHeis = embed(spinchain, [1,2] , J*(XX+YY+ZZ))
for ispin in 2:N-1
    global HHeis
    HHeis += embed(spinchain, [ispin,ispin+1] , J*(XX+YY+ZZ))
end

# make sure H is numerically Hermitian
HHeis = ( HHeis + dagger(HHeis) )/2

s5z = embed(spinchain, [5] , sz)
s6z = embed(spinchain, [6] , sz)

sup = spinup(b)
sdo = spindown(b)

psi0 = sup
for ispin in 2:N
    global psi0
    if ispin == 5
       psi0 = psi0 ⊗ sdo
    else
       psi0 = psi0 ⊗ sup
    end
end

times = 0:.1:1.5

tout, psi_t = timeevolution.schroedinger(times, psi0, HHeis)

m5z = expect(s5z, psi_t)
m6z = expect(s6z, psi_t)

plot(tout, m5z)
plot(tout, m6z)

max_stars_repo_name: BlackBodyRadiation/JupyterTNN

max_issues_repo_name: BlackBodyRadiation/JupyterTNN

max_forks_repo_name: BlackBodyRadiation/JupyterTNN


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Distributed, Random
using ArraysOfArrays, Distributions, PDMats, StatsBase

@testset "mcmc_sample" begin
    mvec = [-0.3, 0.3]
    cmat = [1.0 1.5; 1.5 4.0]
    Σ = @inferred PDMat(cmat)
    mv_dist = MvNormal(mvec, Σ)
    likelihood = @inferred BAT.DistributionDensity(mv_dist)
    bounds = @inferred BAT.HyperRectBounds([-5, -8], [5, 8])
    prior = BAT.ConstDensity(LogDVal(0), bounds)
    nchains = 4
    nsteps = 10^4

    algorithmMW = @inferred(MCMCSampling(mcalg = MetropolisHastings(), trafo = NoDensityTransform(), nchains = nchains, nsteps = nsteps))

    samples = bat_sample(PosteriorDensity(likelihood, prior), algorithmMW).result
    @test BAT.likelihood_pvalue(mv_dist, samples) > 10^-3
    @test (nchains * nsteps - sum(samples.weight)) < 100

    algorithmPW = @inferred MCMCSampling(mcalg = MetropolisHastings(weighting = ARPWeighting()), trafo = NoDensityTransform(), nsteps = 10^5)

    samples, chains = bat_sample(mv_dist, algorithmPW)
    @test BAT.likelihood_pvalue(mv_dist, samples) > 10^-3

    gensamples(rng::AbstractRNG) = bat_sample(rng, PosteriorDensity(mv_dist, prior), algorithmPW).result

    rng = bat_rng()
    @test gensamples(rng) != gensamples(rng)
    @test gensamples(deepcopy(rng)) == gensamples(deepcopy(rng))

    samples = bat_sample(Normal(), MCMCSampling(mcalg = MetropolisHastings(), trafo = NoDensityTransform(), nsteps = 10^4)).result
    @test BAT.likelihood_pvalue(Normal(), samples) > 10^-3
end


text:
cmd = `$(s.graphvizexecutable) "-V"`
gres = run(cmd)

(path,io) = mktemp()
run(gcmd ,stdout=io)

run(pipeline(gcmd,stdout=io));

tmptxt = mktemp()[1]
tmperr = mktemp()[1]
pip = pipeline(cmd, stdout=tmptxt, stderr=tmperr)
rc = run(pip)