text:
# # Approximate Inference using ESS

# Loading the necessary packages and setting seed.

using AbstractGPs, Plots, Random
Random.seed!(1234);

# Loading toy regression
# [dataset](https://github.com/GPflow/docs/blob/master/doc/source/notebooks/basics/data/regression_1D.csv)
# taken from GPFlow examples.

x = [0.8658165855998895, 0.6661700880180962, 0.8049218148148531, 0.7714303440386239,
    0.14790478354654835, 0.8666105548197428, 0.007044577166530286, 0.026331737288148638,
    0.17188596617099916, 0.8897812990554013, 0.24323574561119998, 0.028590102134105955];
y = [1.5255314337144372, 3.6434202968230003, 3.010885733911661, 3.774442382979625,
    3.3687639483798324, 1.5506452040608503, 3.790447985799683, 3.8689707574953,
    3.4933565751758713, 1.4284538820635841, 3.8715350915692364, 3.7045949061144983];
scatter(x, y, xlabel="x", ylabel="y")

# Split the observations into train and test set.

(x_train, y_train) = (x[begin:8], y[begin:8]);
(x_test, y_test) = (x[9:end], y[9:end]);

# Instantiating the kernel.

k = Matern52Kernel()

# Instantiating a Gaussian Process with the given kernel `k`.

f = GP(k)

# Instantiating a `FiniteGP`, a finite dimentional projection at the inputs of the dataset
# observed under Gaussian Noise with $\sigma = 0.001$ .

fx = f(x_train, 0.001)

# Data's log-likelihood w.r.t prior `GP`.

logpdf(fx, y_train)

# Calculating the exact posterior over `f` given `y`. The GP's kernel currently has some
# arbitrary fixed parameters.

p_fx = posterior(fx, y_train)

# Data's log-likelihood under the posterior `GP`. We see that it drastically increases.

logpdf(p_fx(x_test), y_test)

# Plot the posterior `p_fx` along with the observations.

plt = scatter(x_train, y_train; label = "Train data")
scatter!(plt, x_test, y_test; label = "Test data")
plot!(plt, p_fx, 0:0.001:1; label="Posterior")

# # Elliptical Slice Sampler
# Previously we computed the log likelihood of the untuned kernel parameters
# of the GP, $-1.285$. We now also perform approximate inference over said
# kernel parameters using the
# [Elliptical Slice Sampling](http://proceedings.mlr.press/v9/murray10a/murray10a.pdf)
# provided by
# [EllipticalSliceSampling.jl](https://github.com/TuringLang/EllipticalSliceSampling.jl/).
# We start of by loading necessary packages.

using EllipticalSliceSampling, Distributions

# We define a function which returns log-probability of the data under the
# GP / log-likelihood of the parameters of the GP.

function logp(params; x=x_train, y=y_train)
    kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(exp(params[1]))
        ),
        exp(params[2])
    )
    f = GP(kernel)
    fx = f(x, 0.1)
    return logpdf(fx, y)
end

# We define a Gaussian prior over the joint distribution on kernel parameters
# space. Since we have only two parameters, we define a multi-variate Gaussian
# of dimension two.

prior = MvNormal(2, 1)

# Sanity check for the defined `logp` function and `prior` distribution.

logp(rand(prior))

# Generate 2,000 samples using `ESS_mcmc` provided by `EllipticalSliceSampling.jl`.

samples = sample(ESSModel(prior, logp), ESS(), 2_000; progress=false)
samples_mat = reduce(hcat, samples)';

# Mean of samples of both the parameters.

mean_params = mean(samples_mat, dims=1)

# Plot a histogram of the samples for the two parameters. The vertical line in each
# graph indicates the mean of the samples.

plt = histogram(samples_mat; layout=2, labels="Param")
vline!(plt, mean_params; layout=2, label="Mean")

# Average log-marginal-probability of data with posterior kernel parameter samples
# sampled using ESS. We can observe that there is significant improvement over
# exact posterior with default kernel parameters.

mean(logp(param; x=x_test, y=y_test) for param in samples)

# Plot sampled functions from posterior with tuned parameters

plt = scatter(x_train, y_train; label="Train data")
scatter!(plt, x_train, y_train; label="Test data")
for params in @view(samples[(end-100):end,:])
    opt_kernel = ScaledKernel(
        transform(
            Matern52Kernel(),
            ScaleTransform(exp(params[1]))
        ),
        exp(params[2])
    )
    f = GP(opt_kernel)
    p_fx = posterior(f(x, 0.1), y)
    sampleplot!(plt, p_fx(collect(0:0.02:1)), 1)
end
plt


text:
using Documenter, QuadraticTools

makedocs(sitename="QuadraticTools Docs")

ext: jl


text:
__precompile__()
module JAMSDWriter

import MathProgBase
const MPB = MathProgBase
#importall MathProgBase.SolverInterface

using Compat
using Compat.LinearAlgebra, Compat.SparseArrays

if VERSION >= v"0.7"
    # quickfix for Nullable
    using Nullables
end

export JAMSDSolver, getsolvername, getsolveresult, getsolveresultnum, getsolvemessage, getsolveexitcode, LinearQuadraticModel

# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("JAMSDWriter not installed properly, run Pkg.build(\"JAMSDWriter\"), restart Julia and try again")
end

include(depsjl_path)

const CONFIG = Dict(
:debug => false,
:export_gms => false,
:solver_log => false
)

const solverdata_dir = joinpath(dirname(@__DIR__), ".solverdata")

include("jamsd_linearity.jl")
include("jamsd_params.jl")
include("jamsd_convert.jl")
include("jamsd_fun.jl")

const solver_stat = [
   :Optimal,
   :Iteration,
   :Resource,
   :Solver,
   :EvalError,
   :Capability,
   :License,
   :User,
   :SetupErr,
   :SolverErr,
   :InternalErr,
   :Skipped,
   :SystemErr
]

const model_stat = [
   :OptimalGlobal,
   :OptimalLocal,
   :Unbounded,
   :InfeasibleGlobal,
   :InfeasibleLocal,
   :InfeasibleIntermed,
   :Feasible,
   :Integer,
   :NonIntegerIntermed,
   :IntegerInfeasible,
   :LicenseError,
   :ErrorUnknown,
   :ErrorNoSolution,
   :NoSolutionReturned,
   :SolvedUnique,
   :Solved,
   :SolvedSingular,
   :UnboundedNoSolution,
   :InfeasibleNoSolution
]

mutable struct JAMSDSolver <: MPB.AbstractMathProgSolver
    solver_name::String
    options::Dict{String,Any}
    emp::Nullable{Function}
end

# TODO(xhub) write a better struct/enum here
@enum MODEL_TYPE qcp=7 nlp=2 miqcp=6 minlp=5 emp=9

"change the debug state"
setdebug(b::Bool) = global CONFIG[:debug] = b
"export the problem to a GAMS Model file (.gms)"
setexport(b::Bool) = global CONFIG[:export_gms] = b
"printout the log from the solver"
setsolverlog(b::Bool) = global CONFIG[:solver_log] = b

"""
Create a JAMSDSolver Solver for MPB. The optional arguments are:

# Optional Arguments
- `solver_name::String=""`: solver used for this problem
- `options::Dict{String,Any}=Dict{String,Any}()`: the JAMSD options

"""
function JAMSDSolver(solver_name::String="",
                     options::Dict{String,Any}=Dict{String,Any}())
    JAMSDSolver(solver_name, options, Nullable{Function}())
end

getsolvername(s::JAMSDSolver) = basename(s.solver_name)

mutable struct JAMSDMathProgModel <: MPB.AbstractMathProgModel
    options::Dict{String, Any}

    solver_name::String

    x_l::Vector{Float64}
    x_u::Vector{Float64}
    g_l::Vector{Float64}
    g_u::Vector{Float64}

    nvar::Int
    ncon::Int

    obj
    constrs::Vector{Any}

    lin_constrs::Vector{Dict{Int, Float64}}
    lin_obj::Dict{Int, Float64}

    r_codes::Vector{Int}
    j_counts::Vector{Int}

    vartypes::Vector{Symbol}
    varlinearities_con::Vector{Symbol}
    varlinearities_obj::Vector{Symbol}
    conlinearities::Vector{Symbol}
    objlinearity::Symbol

    v_index_map::Dict{Int, Int}
    v_index_map_rev::Dict{Int, Int}

    # was:
    # c_index_map::Dict{Int, Int}
    # c_index_map_rev::Dict{Int, Int}

    nonquad_idx::Dict{Int, Int}
    quad_idx::Dict{Int, Int}

    sense::Symbol

    x_0::Vector{Float64}

    objval::Float64
    solution::Vector{Float64}

    status::Symbol
    solve_exitcode::Int
    solve_result_num::Int
    solve_result::String
    model_result_num::Int
    model_result::String
    solve_message::String
    solve_time::Float64

    model_type::MODEL_TYPE

    quad_equs::Vector{Any}
#    quad_obj::Tuple{Vector{Int}, Vector{Int}, Vector{Float64}}
    quad_obj::Tuple
    offset::Int
    emp::Nullable{Function}

    d::Nullable{MPB.AbstractNLPEvaluator}

    jamsd_ctx::Ptr{context}
    jamsd_ctx_dest::Ptr{context}
    jamsd_options::Ptr{jamsd_options}
    gams_dir::String

    function JAMSDMathProgModel(solver_name::String,
                                options::Dict{String,Any},
                                model_type::MODEL_TYPE,
                                emp)
        o = new(options,
            solver_name,
            zeros(0),
            zeros(0),
            zeros(0),
            zeros(0),
            0,
            0,
            :(0),
            [],
            Dict{Int, Float64}[],
            Dict{Int, Float64}(),
            Int[],
            Int[],
            Symbol[],
            Symbol[],
            Symbol[],
            Symbol[],
            :Lin,
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            :Min,
            zeros(0),
            NaN,
            zeros(0),
            :NotSolved,
            -1,
            -1,
            "?",
            -1,
            "?",
            "",
            NaN,
            model_type,
            Vector{Any}(),
            (),
            0,
            emp,
            Nullable{MPB.AbstractNLPEvaluator}(),
            Ptr{context}(C_NULL),
            Ptr{context}(C_NULL),
            Ptr{jamsd_options}(C_NULL),
            "")
        @compat finalizer(jamsd_cleanup, o)
        o
    end
end

struct JAMSDLinearQuadraticModel <: MPB.AbstractLinearQuadraticModel
    inner::JAMSDMathProgModel
end

struct JAMSDNonlinearModel <: MPB.AbstractNonlinearModel
    inner::JAMSDMathProgModel
end

struct JAMSDConicModel <: MPB.AbstractNonlinearModel
    inner::JAMSDMathProgModel
end

include("jamsd_write.jl")

MPB.NonlinearModel(s::JAMSDSolver) = JAMSDNonlinearModel(
    JAMSDMathProgModel(s.solver_name, s.options, nlp, s.emp)
)

MPB.LinearQuadraticModel(s::JAMSDSolver) = JAMSDLinearQuadraticModel(
    JAMSDMathProgModel(s.solver_name, s.options, qcp, s.emp)
)

function MPB.ConicModel(s::JAMSDSolver)
    error("ConicModel is not yet supported")
end

function MPB.loadproblem!(outer::JAMSDNonlinearModel, nvar::Integer, ncon::Integer,
                      x_l, x_u, g_l, g_u, sense::Symbol,
                      d::MPB.AbstractNLPEvaluator)
    m = outer.inner

    m.nvar, m.ncon = nvar, ncon
    loadcommon!(m, x_l, x_u, g_l, g_u, sense)

    m.d = d

    MPB.initialize(m.d.value, [:ExprGraph])

    # Process constraints
    m.constrs = map(1:m.ncon) do i
        c = MPB.constr_expr(m.d.value, i)

        # Remove relations and bounds from constraint expressions
        if length(c.args) == 3
            if VERSION < v"0.5-"
                expected_head = :comparison
                expr_index = 1
                rel_index = 2
            else
                expected_head = :call
                expr_index = 2
                rel_index = 1
            end

            @assert c.head == expected_head
            # Single relation constraint: expr rel bound
            rel = c.args[rel_index]
            m.r_codes[i] = relation_to_jamsd[rel]
            if rel == [:<=, :(==)]
                m.g_u[i] = c.args[3]
            end
            if rel in [:>=, :(==)]
                m.g_l[i] = c.args[3]
            end
            c = c.args[expr_index]
        else
            # Double relation constraint: bound <= expr <= bound
            @assert c.head == :comparison
            m.r_codes[i] = relation_to_jamsd[:multiple]
            m.g_u[i] = c.args[5]
            m.g_l[i] = c.args[1]
            c = c.args[3]
        end

        # Convert non-linear expression to non-linear, linear and constant
        c, constant, m.conlinearities[i] = process_expression!(
            c, m.lin_constrs[i], m.varlinearities_con)

        # Update bounds on constraint
        m.g_l[i] -= constant
        m.g_u[i] -= constant

        # Update jacobian counts using the linear constraint variables
        for j in keys(m.lin_constrs[i])
            m.j_counts[j] += 1
        end
        c
    end

    # Process objective
    m.obj = MPB.obj_expr(m.d.value)

    if length(m.obj.args) < 2
        m.obj = 0
    else
        # Convert non-linear expression to non-linear, linear and constant
        m.obj, constant, m.objlinearity = process_expression!(
            m.obj, m.lin_obj, m.varlinearities_obj)

        # Add constant back into non-linear expression
        if constant != 0
            m.obj = add_constant(m.obj, constant)
        end
    end
    m
end

function MPB.loadproblem!(outer::JAMSDLinearQuadraticModel, A::AbstractMatrix,
                      x_l, x_u, c, g_l, g_u, sense)
    m = outer.inner
    m.ncon, m.nvar = size(A)

    loadcommon!(m, x_l, x_u, g_l, g_u, sense)

    # Load A into the linear constraints
    @assert (m.ncon, m.nvar) == size(A)
    load_A!(m, A)
    m.constrs = zeros(m.ncon)  # Dummy constraint expression trees

    # Load c
    for (index, val) in enumerate(c)
        m.lin_obj[index] = val
    end
    # TODO(xhub) see if we can get rid of that
    m.obj = 0  # Dummy objective expression tree

    # Process variables bounds
    for j = 1:m.ncon
        lower = m.g_l[j]
        upper = m.g_u[j]
        if lower == -Inf
            if upper == Inf
                error("Neither lower nor upper bound on constraint $j")
            else # <=
                m.r_codes[j] = 2
            end
        else
            if lower == upper  # ==
                m.r_codes[j] = 0
            elseif upper == Inf # >=
                m.r_codes[j] = 1
            else # lb <= expr <= ub
                m.r_codes[j] = -1
            end
        end
    end
    m
end

function load_A!(m::JAMSDMathProgModel, A::SparseMatrixCSC{Float64})
    for var = 1:A.n, k = A.colptr[var] : (A.colptr[var + 1] - 1)
        m.lin_constrs[A.rowval[k]][var] = A.nzval[k]
        m.j_counts[var] += 1
    end
end

function load_A!(m::JAMSDMathProgModel, A::Matrix{Float64})
    for con = 1:m.ncon, var = 1:m.nvar
        val = A[con, var]
        if val != 0
            m.lin_constrs[A.rowval[k]][var] = A.nzval[k]
            m.j_counts[var] += 1
        end
    end
end

function loadcommon!(m::JAMSDMathProgModel, x_l, x_u, g_l, g_u, sense)
    m.x_l, m.x_u = x_l, x_u
    m.g_l, m.g_u = g_l, g_u
    setsense!(m, sense)

    m.lin_constrs = [Dict{Int, Float64}() for _ in 1:m.ncon]
    m.j_counts = zeros(Int, m.nvar)

    m.r_codes = Vector{Int}(undef, m.ncon)

    m.varlinearities_con = fill(:Lin, m.nvar)
    m.varlinearities_obj = fill(:Lin, m.nvar)
    m.conlinearities = fill(:Lin, m.ncon)
    m.objlinearity = :Lin

    m.vartypes = fill(:Cont, m.nvar)
    m.x_0 = zeros(m.nvar)
end

getvartype(m::JAMSDMathProgModel) = copy(m.vartypes)
function setvartype!(m::JAMSDMathProgModel, cat::Vector{Symbol})
    @assert all(x-> (x in [:Cont,:Bin,:Int,:external]), cat)
    m.vartypes = copy(cat)
end

getsense(m::JAMSDMathProgModel) = m.sense
function setsense!(m::JAMSDMathProgModel, sense::Symbol)
    @assert sense == :Min || sense == :Max
    m.sense = sense
end

setwarmstart!(m::JAMSDMathProgModel, v::Vector{Float64}) = m.x_0 = v

function MPB.addquadconstr!(m::JAMSDLinearQuadraticModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
   # we have to do a little translation of the sense here ...
   push!(m.inner.quad_equs, tuple(linearidx,linearval,quadrowidx,quadcolidx,quadval, quad_relation_sense[sense], rhs))
end

function MPB.setquadobj!(m::JAMSDLinearQuadraticModel, rowidx, colidx, quadval)
   m.inner.quad_obj = tuple(rowidx, colidx, quadval)
end

function optimize!(m::JAMSDMathProgModel)
    m.status = :NotSolved
    m.solve_exitcode = -1
    m.solve_result_num = -1
    m.solve_result = "?"
    m.model_result_num = -1
    m.model_result = "?"
    m.solve_message = ""

    # There is no non-linear binary type, only non-linear discrete, so make
    # sure binary vars have bounds in [0, 1]
    for i in 1:m.nvar
        if m.vartypes[i] == :Bin
            if m.x_l[i] < 0
                m.x_l[i] = 0
            end
            if m.x_u[i] > 1
                m.x_u[i] = 1
            end
        end
    end

    m.jamsd_options = jamsd_options_set(m.options)

    if m.emp.hasvalue
       return m.emp.value()
    end

    make_var_index!(m)
    make_con_index!(m)

    # Run solver and save exitcode
    t = time()
    m.jamsd_ctx = create_jamsd_ctx(m)
    jamsd_set_modeltype(m)
    # Solve via gams for now
    m.jamsd_ctx_dest, m.gams_dir = jamsd_setup_gams()

    m.solve_exitcode = jamsd_solve(m.jamsd_ctx, m.jamsd_ctx_dest, m.solver_name)

#    ccall((:print_model, libjamsd), Cint, (Ptr{context},), m.jamsd_ctx)
    m.solve_time = time() - t

    if m.solve_exitcode == 0
        report_results(m)
    else
        println("JAMSD: solver failed with status $(m.solve_exitcode)")
        m.status = :Error
        m.solution = fill(NaN, m.nvar)
        m.solve_result = "failure"
        m.solve_result_num = 999
    end
end

function getconstrduals(m::JAMSDMathProgModel)

    ctx = m.jamsd_ctx

    x = fill(NaN, numconstr(m))

   if has_objective(m)
        offset = m.offset
    else
        offset = m.offset - 1
    end

    for idx in 1:m.ncon
       eidx = m.nonquad_idx[idx] + offset
       x[idx] = ctx_getmultiplierval(ctx, eidx)
    end

    return x
end

function getquadconstrduals(quadm::JAMSDLinearQuadraticModel)

    m = quadm.inner
    ctx = m.jamsd_ctx

    x = fill(NaN, numconstr(m))

    if has_objective(m)
        offset = m.offset
    else
        offset = m.offset - 1
    end

    for (idx, equ) in enumerate(m.quad_equs)
       eidx = m.quad_idx[idx] + offset
       x[idx] = ctx_getmultiplierval(ctx, eidx)
    end

    return x
end

function getreducedcosts(m::JAMSDMathProgModel)

    ctx = m.jamsd_ctx

    x = fill(NaN, numvar(m))

    for idx in 1:numvar(m)
       x[idx] = ctx_getvarmult(ctx, idx-1)
    end

    return x
end

MPB.getreducedcosts(nlpm::JAMSDNonlinearModel) = getreducedcosts(nlpm.inner)
MPB.getreducedcosts(quadm::JAMSDLinearQuadraticModel) = getreducedcosts(quadm.inner)
MPB.getconstrduals(nlpm::JAMSDNonlinearModel) = getconstrduals(nlpm.inner)
MPB.getconstrduals(quadm::JAMSDLinearQuadraticModel) = getconstrduals(quadm.inner)

function process_expression!(nonlin_expr::Expr, lin_expr::Dict{Int, Float64},
                             varlinearities::Vector{Symbol})
    # Get list of all variables in the expression
    extract_variables!(lin_expr, nonlin_expr)
    # Extract linear and constant terms from non-linear expression
    tree = LinearityExpr(nonlin_expr)
    tree = pull_up_constants(tree)
    _, tree, constant = prune_linear_terms!(tree, lin_expr)
    # Make sure all terms remaining in the tree are .nl-compatible
    nonlin_expr = convert_formula(tree)

    # Track which variables appear nonlinearly
    nonlin_vars = Dict{Int, Float64}()
    extract_variables!(nonlin_vars, nonlin_expr)
    for j in keys(nonlin_vars)
        varlinearities[j] = :Nonlin
    end

    # Remove variables at coeff 0 that aren't also in the nonlinear tree
    for (j, coeff) in lin_expr
        if coeff == 0 && !(j in keys(nonlin_vars))
            delete!(lin_expr, j)
        end
    end

    # Mark constraint as nonlinear if anything is left in the tree
    linearity = nonlin_expr != 0 ? :Nonlin : :Lin

    return nonlin_expr, constant, linearity
end
function process_expression!(nonlin_expr::Real, lin_expr, varlinearities)
    # Special case where body of constraint is constant
    # Return empty nonlinear and linear parts, and use the body as the constant
    0, nonlin_expr, :Lin
end

status(m::JAMSDMathProgModel) = m.status
getsolution(m::JAMSDMathProgModel) = copy(m.solution)
getobjval(m::JAMSDMathProgModel) = m.objval
numvar(m::JAMSDMathProgModel) = m.nvar
numconstr(m::JAMSDMathProgModel) = m.ncon + length(m.quad_equs)
getsolvetime(m::JAMSDMathProgModel) = m.solve_time

# Access to solve results
get_solve_result(m::JAMSDMathProgModel) = m.solve_result
get_solve_result_num(m::JAMSDMathProgModel) = m.solve_result_num
get_model_result(m::JAMSDMathProgModel) = m.model_result
get_model_result_num(m::JAMSDMathProgModel) = m.model_result_num
get_solve_message(m::JAMSDMathProgModel) = m.solve_message
get_solve_exitcode(m::JAMSDMathProgModel) = m.solve_exitcode

# We need to track linear coeffs of all variables present in the expression tree
extract_variables!(lin_constr::Dict{Int, Float64}, c) = c
extract_variables!(lin_constr::Dict{Int, Float64}, c::LinearityExpr) =
    extract_variables!(lin_constr, c.c)
function extract_variables!(lin_constr::Dict{Int, Float64}, c::Expr)
    if c.head == :ref
        CONFIG[:debug] && println("DEBUG: extract_variables :: variable case: $c")
        if c.args[1] == :x
            @assert isa(c.args[2], Int)
            lin_constr[c.args[2]] = 0
        else
            error("Unrecognized reference expression $c")
        end
    else
        map(arg -> extract_variables!(lin_constr, arg), c.args)
    end
end

add_constant(c, constant::Real) = c + constant
add_constant(c::Expr, constant::Real) = Expr(:call, :+, c, constant)

function make_var_index!(m::JAMSDMathProgModel)
    m.v_index_map = Dict(zip(1:m.nvar, 0:(m.nvar-1)))
    m.v_index_map_rev = Dict(zip(0:(m.nvar-1), 1:m.nvar))
end

#function make_var_index!(m::JAMSDMathProgModel)
#    nonlin_cont = Int[]
#    nonlin_int = Int[]
#    lin_cont = Int[]
#    lin_int = Int[]
#    lin_bin = Int[]
#
#    # TODO(xhub) we do that multiple times in the EMP context
#    for i in 1:length(m.vartypes)
#        if m.varlinearities_obj[i] == :Nonlin ||
#           m.varlinearities_con[i] == :Nonlin
#            if m.vartypes[i] == :Cont || m.vartypes[i] == :external
#                push!(nonlin_cont, i)
#            else
#                push!(nonlin_int, i)
#            end
#        else
#           if m.vartypes[i] == :Cont || m.vartypes[i] == :external
#                push!(lin_cont, i)
#            elseif m.vartypes[i] == :Int
#                push!(lin_int, i)
#            else
#                push!(lin_bin, i)
#            end
#        end
#    end
#
#    # Index variables in required order
#    for var_list in (nonlin_cont, nonlin_int, lin_cont, lin_bin, lin_int)
#        add_to_index_maps!(m.v_index_map, m.v_index_map_rev, var_list, 0)
#    end
#    CONFIG[:debug] && println("DEBUG: $(m.v_index_map)")
#end

function make_con_index!(m::JAMSDMathProgModel)
    nonlin_cons = Int[]
    lin_cons = Int[]

    for i in 1:m.ncon
        if m.conlinearities[i] == :Nonlin
            push!(nonlin_cons, i)
        else
            push!(lin_cons, i)
        end
    end
    for con_list in (nonlin_cons, lin_cons)
        add_to_index_maps!(m.nonquad_idx, con_list, 1)
    end

    if length(m.quad_idx) == 0
       m.quad_idx = Dict(enumerate(Vector{Int}(Compat.range(1+m.ncon, length=length(m.quad_equs)))))
    end

    CONFIG[:debug] && println("DEBUG: make_con_index: nonquad_idx: $(m.nonquad_idx)\n quad_idx: $(m.quad_idx)")
end

function add_to_index_maps!(forward_map::Dict{Int, Int},
                            backward_map::Dict{Int, Int},
                            inds::Array{Int},
                            offset::Int)
    for i in inds
        index = length(forward_map) + offset
        forward_map[i] = index
        backward_map[index] = i
    end
end

function add_to_index_maps!(forward_map::Dict{Int, Int},
                            inds::Array{Int},
                            offset::Int)
    for i in inds
        index = length(forward_map) + offset
        forward_map[i] = index
    end
end

function report_results_common(m::JAMSDMathProgModel)
    x = fill(NaN, m.nvar)
    m.objval = NaN

    for index in 0:(m.nvar - 1)
        i = m.v_index_map_rev[index]
        x[i] = ctx_getvarval(m.jamsd_ctx, index)
    end

    m.solution = x

    ###########################################################################
    # Convert solve_result
    #
    # GAMS return two information:
    #  - the solve status
    #  - the model status
    #
    # - :Optimal
    # - :Infeasible
    # - :Unbounded
    # - :UserLimit (iteration limit or timeout)
    # - :Error (and maybe others)
    ###########################################################################

    tmpCint = Ref{Cint}(0)
    res = ccall((:ctx_getsolvestat, libjamsd), Cint, (Ptr{context}, Ref{Cint}), m.jamsd_ctx_dest, tmpCint)
    res != 0 && error("return code $res from JAMSD")
    m.solve_result_num = tmpCint.x
    m.solve_result = unsafe_string(ccall((:ctx_getsolvestattxt, libjamsd), Cstring, (Ptr{context}, Cint), m.jamsd_ctx_dest, m.solve_result_num))

    res = ccall((:ctx_getmodelstat, libjamsd), Cint, (Ptr{context}, Ref{Cint}), m.jamsd_ctx_dest, tmpCint)
    res != 0 && error("return code $res from JAMSD")

    m.model_result_num = tmpCint.x

    m.model_result = unsafe_string(ccall((:ctx_getmodelstattxt, libjamsd), Cstring, (Ptr{context}, Cint), m.jamsd_ctx_dest, m.model_result_num))

    # GAMS already uses an 1-indices
    solver_code = solver_stat[m.solve_result_num]
    model_code = model_stat[m.model_result_num]

    CONFIG[:debug] && println("solver stat $(m.solve_result) ($(m.solve_result_num)); model stat $(m.model_result) ($(m.model_result_num))")

    if solver_code == :Optimal
       if model_code == :OptimalGlobal || model_code == :OptimalLocal || model_code == :Integer
          m.status = :Optimal
       elseif model_code == :Unbounded || model_code == :UnboundedNoSolution
          m.status = :Unbounded
       elseif model_code == :InfeasibleGlobal || model_code == :InfeasibleLocal || model_code == :InfeasibleIntermed || model_code == :InfeasibleNoSolution
          m.status = :Infeasible
       elseif model_code == :Feasible
          # TODO investigate that. Baron is weird
          m.status = :Optimal
       elseif model_code == :NoSolutionReturned
          gams_solver = ctx_get_solvername(m.jamsd_ctx_dest)
          m.status = :Optimal
#          if gams_solver == "jams" || gams_solver == "JAMS"
             # This is fine, we have a kludge in the code
#          else
#             println("JAMSD: Solve successed, but no solution was returned by solver $(gams_solver)!")
#             m.status = :Error
#          end
       else
          println("JAMSD: unhandle case: solver stat $(m.solve_result); model stat $(m.model_result)")
          m.status = :Error
       end
    elseif solver_code == :Iteration || solver_code == :Resource
       m.status == :UserLimit
    elseif solver_code == :License || model_code == :LicenseError
       println("JAMSD: License error. Check that you have a valid license")
       m.status == :Error
    elseif solver_code == :Capability
       if m.solver_name == ""
          sname = default
       else
          sname = m.solver_name
       end

       println("JAMSD: solver $(name) cannot solve the specific problem")
       m.status = :Error

    else
       println("JAMSD: solver stat is $(m.solve_result) and model stat is $(m.model_result)")
       m.status = :Error
    end

    CONFIG[:debug] && println("status is $(m.status)")
 end

function report_results(m::JAMSDMathProgModel)
    # TODO(Xhub) fix this hack
    res = ccall((:model_eval_eqns, libjamsd), Cint, (Ptr{context}, Ptr{context}), m.jamsd_ctx, m.jamsd_ctx_dest)
    res != 0 && error("JAMSD: error code $res")

    # Next, read for the variable values
    report_results_common(m)

    # TODO(xhub) this should not be necessary
    if m.solve_exitcode == 0
        if m.objlinearity == :Nonlin
            # Try to use NLPEvaluator if we can.
            # Can fail due to unsupported functions so fallback to eval
            try
                m.objval = eval_f(m.d.value, m.solution)
            catch
                CONFIG[:debug] && println("Error: could not evaluate the objective function")
            end
        end

        # Calculate objective value from nonlinear and linear parts
        obj_nonlin = eval(substitute_vars!(deepcopy(m.obj), m.solution))
        obj_lin = evaluate_linear(m.lin_obj, m.solution)
        if (length(m.quad_obj) == 3)
           ridx, cidx, vals = m.quad_obj
           obj_quad = evaluate_quad(ridx, cidx, vals, m.solution)
        else
           obj_quad = 0.
        end
        m.objval = obj_nonlin + obj_lin + obj_quad
    end
end

substitute_vars!(c, x::Array{Float64}) = c
function substitute_vars!(c::Expr, x::Array{Float64})
    if c.head == :ref
        if c.args[1] == :x
            index = c.args[2]
            @assert isa(index, Int)
            c = x[index]
        else
            error("Unrecognized reference expression $c")
        end
    else
        if c.head == :call
            # Convert .nl unary minus (:neg) back to :-
            if c.args[1] == :neg
                c.args[1] = :-
            # Convert .nl :sum back to :+
            elseif c.args[1] == :sum
                c.args[1] = :+
            end
        end
        map!(arg -> substitute_vars!(arg, x), c.args, c.args)
    end
    c
end

function evaluate_linear(linear_coeffs::Dict{Int, Float64}, x::Array{Float64})
    total = 0.0
    for (i, coeff) in linear_coeffs
        total += coeff * x[i]
    end
    total
end

function evaluate_quad(rowidx, colidx, qvals, x::Array{Float64})
   n = length(x)
   mat = sparse(rowidx, colidx, qvals, n, n)
   # This is soooo ugly --xhub
   Q = (mat + mat') - Diagonal(diag(mat))
   total = .5*x'*Q*x
   return total
end

# Wrapper functions
for f in [:getvartype,:getsense,:optimize!,:status,:getsolution,:getobjval,:numvar,:numconstr,:getsolvetime]
    @eval MPB.$f(m::JAMSDNonlinearModel) = $f(m.inner)
    @eval MPB.$f(m::JAMSDLinearQuadraticModel) = $f(m.inner)
end

for f in [:get_solve_result,:get_solve_result_num,:get_solve_message,:get_solve_exitcode]
    @eval $f(m::JAMSDNonlinearModel) = $f(m.inner)
    @eval $f(m::JAMSDLinearQuadraticModel) = $f(m.inner)
end
for f in [:setvartype!,:setsense!,:setwarmstart!]
    @eval MPB.$f(m::JAMSDNonlinearModel, x) = $f(m.inner, x)
    @eval MPB.$f(m::JAMSDLinearQuadraticModel, x) = $f(m.inner, x)
end

# Utility method for deleting any leftover debug files
function clean_solverdata()
    for file in readdir(solverdata_dir)
        ext = splitext(file)[2]
        (ext == ".nl" || ext == ".sol") && rm(joinpath(solverdata_dir, file))
    end
end

include("jamsd_mathprgm.jl")
include("jamsd_ovf.jl")
include("jamsd_solve.jl")

function jamsd_cleanup(o::JAMSDMathProgModel)
    ctx_dealloc(o.jamsd_ctx)
    ctx_dealloc(o.jamsd_ctx_dest)
    jamsd_options_dealloc(o.jamsd_options)
    if (!isempty(o.gams_dir))
        try
            rm(o.gams_dir, recursive=true, force=true)
        catch
            iswin && run(`cmd /C RMDIR /s /q $(o.gams_dir)`)
        end
    end
end

function jamsd_options_set(opt::Dict{String,Any})
    jopt = jamsd_options_alloc()

    for (k,v) in opt
        jamsd_option_set(jopt, k, v)
    end

    return jopt
end

function jamsd_setemp!(m::JAMSDSolver, emp)
   m.emp = emp
end

end


text:
@testset "compute_absorbing_boundary_condition_matrix" begin
    mmesh = Mesh(10,10,0.1)
    bdedge = bcedge(mmesh)
    N = size(bdedge, 1)
    ρ = ones(4N)
    cs = ones(4N)
    cp = ones(4N)
    A = compute_absorbing_boundary_condition_matrix(ρ, cs, cp, bdedge, mmesh)
    @test false
end


text:
 	const centuries = ("I", "II", "III", "IV", "V", "VI", "VII", "VIII", "IX", "X", "XI", "XII", "XIII", "XIV", "XV", "XVI","XVII", "XVIII", "XIX", "XX", "XXI")

 	const year_top =2100

 	const ap = ("AM","PM")

 	unix_time()::Int= rand(0: Dates.datetime2unix(Dates.DateTime(Dates.now())))
 	iso8601()::DateTime = Dates.DateTime(rand(1970:year_top),rand(01:12),1,rand(01:23),rand(00:59),rand(00:59)) + Dates.Day(rand(01:31))
 	date_time(pattern::String="Y-m-d HH:MM:SS")::String = Dates.format(iso8601(), pattern)
 	date(pattern::String="Y-m-d")::String = Dates.format(iso8601(), pattern)
 	time(pattern::String="H:M:S")::String = Dates.format(Dates.DateTime(1,1,1,rand(01:23),rand(00:59),rand(00:59)), pattern)

function date_time_between(start_date::String="-30y", finish_date::String="now", pattern::String="Y-m-d HH:MM:SS")::String

 	init_date::DateTime = Dates.DateTime(Dates.now())
    end_date::DateTime = Dates.DateTime(Dates.now())

    if start_date != "now"
    	init_date = init_date + Dates.Year(parse(Int,start_date[1:end-1]))
    end

    if finish_date != "now"
    	end_date = end_date + Dates.Year(parse(Int,start_date[1:end-1]))
    end

	 return Dates.format(rand(collect(init_date:Day(1):end_date)), pattern)
end

function random_datetime(bf_now::Bool,af_now::Bool,time_start::DateTime,now_time::DateTime,time_finish::DateTime)::DateTime
    init_date::DateTime = now_time
    end_date::DateTime = now_time

    if bf_now
    	init_date = time_start
    end

    if af_now
    	end_date = time_finish
    end

    return rand(collect(time_start:Day(1):now_time))
end

function date_time_this_century(;before_now::Bool=true, after_now::Bool=false, pattern::String="Y-m-d HH:MM:SS")::String
	r::DateTime = Dates.DateTime(Dates.now())
	this_century_start::DateTime = Dates.DateTime((Int(Dates.year(r)) - (Dates.year(r) % 100)), 1, 1)
	this_century_finish::DateTime = (Dates.DateTime(Dates.year(this_century_start) + 100, 1, 1))-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_century_start,r,this_century_finish), pattern)
end

function  date_time_this_decade(;before_now::Bool=true, after_now::Bool=false, pattern::String="Y-m-d HH:MM:SS")::String
	r::DateTime = Dates.DateTime(Dates.now())
	this_decade_start::DateTime = Dates.DateTime((Int(Dates.year(r)) - (Dates.year(r) % 10)), 1, 1)
	this_decade_finish::DateTime = (Dates.DateTime(Dates.year(this_decade_start) + 10, 1, 1))-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_decade_start,r,this_decade_finish), pattern)
end

function date_time_this_year(;before_now::Bool=true, after_now::Bool=false, pattern::String="Y-m-d HH:MM:SS")::String
	r::DateTime = Dates.DateTime(Dates.now())
	this_year_start::DateTime = Dates.DateTime(Dates.year(r), 1, 1)
	this_year_finish::DateTime = this_year_start+Dates.Year(1)-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_year_start,r,this_year_finish), pattern)
end

function date_time_this_month(;before_now::Bool=true, after_now::Bool=false, pattern::String="Y-m-d HH:MM:SS")::String
	r::DateTime = Dates.DateTime(Dates.now())
	this_month_start::DateTime = r-Dates.Day(r)+Dates.Day(1)
	this_month_finish::DateTime = this_month_start+Dates.Month(1)-Dates.Day(1)

	return Dates.format(random_datetime(before_now,after_now,this_month_start,r,this_month_finish), pattern)
end

am_pm()::String = rand(ap)
day_of_month()::Int = rand(1:31)
day_of_week()::String = Dates.dayname(Date(2015,07,rand(12:18)))
months()::Int =rand(1:12)
month_name()::String = Dates.monthname(Date(2015,rand(1:12),15))
year()::Int =rand(1920:year_top)
century()::String = rand(centuries)
function timezone()::String
	executor(data["faker"]["address"]["time_zone"])
end


text:
import Pkg;
Pkg.activate(".")
Pkg.add("Test")
Pkg.add("Zygote")
Pkg.add("ForwardDiff")
Pkg.add("Flux")
Pkg.add("Rotations")
# Pkg.add("JLD2")
# Pkg.add("Documenter")
# Pkg.add("Example")
Pkg.add("DSP")
# Pkg.add("Makie")
Pkg.add("GLMakie")
# Pkg.add("CairoMakie")
# Pkg.add("FFTW")

max_stars_repo_name: paulxshen/eno

max_issues_repo_name: paulxshen/eno

max_forks_repo_name: paulxshen/eno


text:
using Distances
# kernelweights.jl: kernel weights for nonparametrics
function kernelweights(x, xeval, bandwidth, prewhiten=true, kernel="gaussian", neighbors=1)
    if prewhiten
        v = cov(x)
        if isposdef(v)
            P = inv(chol(v))
        else
            P = sqrt(inv(diag(v)))
        end
        x *= P
        xeval *= P
    end
    neval = size(xeval,1)
    n = size(x,1)
    weights = zeros(n,neval)
    # Gaussian product kernel, add others if desired
    if kernel=="gaussian"
        for i=1:neval
            z = (x.-xeval[[i],:])/bandwidth
            weights[:,i] = exp.(-0.5*sum(z.*z,2))
        end
    end

    if kernel=="knngaussian"
        distances = pairwise(Euclidean(),x', xeval') # get all distances
        @inbounds  for i = 1:neval
            di = view(distances,:,[i])
            ind = sortperm(di) # indices of k nearest neighbors
            selected = vec(ind[1:neighbors,:])
            z = (x[selected,:].-xeval[[i],:])/bandwidth
            weights[selected,i] = exp.(-0.5*sum(z.*z,2))
        end
    end
    weights ./= sum(weights,1)
    return weights
end


text:
using OrderedCollections
"""
A simple struct to contain a vector of energies, and utilities to print them in a nice format.
"""
struct Energies{T <: Number}
    # energies["TermName"]
    # parametrization on T acts as a nice check that all terms return correct type
    energies::OrderedDict{String, T}
end

function Base.show(io::IO, energies::Energies)
    print(io, "Energies(total = $(energies.total))")
end
function Base.show(io::IO, ::MIME"text/plain", energies::Energies)
    println(io, "Energy breakdown (in Ha):")
    for (name, value) in energies.energies
        @printf io "    %-20s%-10.7f\n" string(name) value
    end
    @printf io "\n    %-20s%-15.12f" "total" energies.total
end
Base.getindex(energies::Energies, i) = energies.energies[i]
Base.values(energies::Energies)      = values(energies.energies)
Base.keys(energies::Energies)        = keys(energies.energies)
Base.pairs(energies::Energies)       = pairs(energies.energies)
Base.iterate(energies::Energies)     = iterate(energies.energies)
Base.iterate(energies::Energies, state) = iterate(energies.energies, state)
Base.haskey(energies::Energies, key) = haskey(energies.energies, key)

function Energies(term_types::Vector, energies::Vector{T}) where {T}
    # nameof is there to get rid of parametric types
    Energies{T}(OrderedDict([string(nameof(typeof(term))) => energies[i]
                             for (i, term) in enumerate(term_types)]...))
end

function Base.propertynames(energies::Energies, private::Bool=false)
    ret = keys(energies)
    append!(ret, "total")
    private && append!(ret, "energies")
end
function Base.getproperty(energies::Energies, x::Symbol)
    x == :total && return sum(values(energies))
    x == :energies && return getfield(energies, x)
    energies.energies[string(x)]
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestAdjointTranspose

using Test, LinearAlgebra, SparseArrays

@testset "Adjoint and Transpose inner constructor basics" begin
    intvec, intmat = [1, 2], [1 2; 3 4]
    # Adjoint/Transpose eltype must match the type of the Adjoint/Transpose of the input eltype
    @test_throws ErrorException Adjoint{Float64,Vector{Int}}(intvec)
    @test_throws ErrorException Adjoint{Float64,Matrix{Int}}(intmat)
    @test_throws ErrorException Transpose{Float64,Vector{Int}}(intvec)
    @test_throws ErrorException Transpose{Float64,Matrix{Int}}(intmat)
    # Adjoint/Transpose wrapped array type must match the input array type
    @test_throws MethodError Adjoint{Int,Vector{Float64}}(intvec)
    @test_throws MethodError Adjoint{Int,Matrix{Float64}}(intmat)
    @test_throws MethodError Transpose{Int,Vector{Float64}}(intvec)
    @test_throws MethodError Transpose{Int,Matrix{Float64}}(intmat)
    # Adjoint/Transpose inner constructor basic functionality, concrete scalar eltype
    @test (Adjoint{Int,Vector{Int}}(intvec)::Adjoint{Int,Vector{Int}}).parent === intvec
    @test (Adjoint{Int,Matrix{Int}}(intmat)::Adjoint{Int,Matrix{Int}}).parent === intmat
    @test (Transpose{Int,Vector{Int}}(intvec)::Transpose{Int,Vector{Int}}).parent === intvec
    @test (Transpose{Int,Matrix{Int}}(intmat)::Transpose{Int,Matrix{Int}}).parent === intmat
    # Adjoint/Transpose inner constructor basic functionality, abstract scalar eltype
    anyvec, anymat = Any[1, 2], Any[1 2; 3 4]
    @test (Adjoint{Any,Vector{Any}}(anyvec)::Adjoint{Any,Vector{Any}}).parent === anyvec
    @test (Adjoint{Any,Matrix{Any}}(anymat)::Adjoint{Any,Matrix{Any}}).parent === anymat
    @test (Transpose{Any,Vector{Any}}(anyvec)::Transpose{Any,Vector{Any}}).parent === anyvec
    @test (Transpose{Any,Matrix{Any}}(anymat)::Transpose{Any,Matrix{Any}}).parent === anymat
    # Adjoint/Transpose inner constructor basic functionality, concrete array eltype
    intvecvec = [[1, 2], [3, 4]]
    intmatmat = [[[1 2]] [[3 4]] [[5 6]]; [[7 8]] [[9 10]] [[11 12]]]
    @test (X = Adjoint{Adjoint{Int,Vector{Int}},Vector{Vector{Int}}}(intvecvec);
            isa(X, Adjoint{Adjoint{Int,Vector{Int}},Vector{Vector{Int}}}) && X.parent === intvecvec)
    @test (X = Adjoint{Adjoint{Int,Matrix{Int}},Matrix{Matrix{Int}}}(intmatmat);
            isa(X, Adjoint{Adjoint{Int,Matrix{Int}},Matrix{Matrix{Int}}}) && X.parent === intmatmat)
    @test (X = Transpose{Transpose{Int,Vector{Int}},Vector{Vector{Int}}}(intvecvec);
            isa(X, Transpose{Transpose{Int,Vector{Int}},Vector{Vector{Int}}}) && X.parent === intvecvec)
    @test (X = Transpose{Transpose{Int,Matrix{Int}},Matrix{Matrix{Int}}}(intmatmat);
            isa(X, Transpose{Transpose{Int,Matrix{Int}},Matrix{Matrix{Int}}}) && X.parent === intmatmat)
end

@testset "Adjoint and Transpose outer constructor basics" begin
    intvec, intmat = [1, 2], [1 2; 3 4]
    # the wrapped array's eltype strictly determines the Adjoint/Transpose eltype
    # so Adjoint{T}/Transpose{T} constructors are somewhat unnecessary and error-prone
    # so ascertain that such calls throw whether or not T and the input eltype are compatible
    @test_throws MethodError Adjoint{Int}(intvec)
    @test_throws MethodError Adjoint{Int}(intmat)
    @test_throws MethodError Adjoint{Float64}(intvec)
    @test_throws MethodError Adjoint{Float64}(intmat)
    @test_throws MethodError Transpose{Int}(intvec)
    @test_throws MethodError Transpose{Int}(intmat)
    @test_throws MethodError Transpose{Float64}(intvec)
    @test_throws MethodError Transpose{Float64}(intmat)
    # Adjoint/Transpose outer constructor basic functionality, concrete scalar eltype
    @test (Adjoint(intvec)::Adjoint{Int,Vector{Int}}).parent === intvec
    @test (Adjoint(intmat)::Adjoint{Int,Matrix{Int}}).parent === intmat
    @test (Transpose(intvec)::Transpose{Int,Vector{Int}}).parent === intvec
    @test (Transpose(intmat)::Transpose{Int,Matrix{Int}}).parent === intmat
    # the tests for the inner constructors exercise abstract scalar and concrete array eltype, forgoing here
end

@testset "Adjoint and Transpose add additional layers to already-wrapped objects" begin
    intvec, intmat = [1, 2], [1 2; 3 4]
    @test (A = Adjoint(Adjoint(intvec))::Adjoint{Int,Adjoint{Int,Vector{Int}}}; A.parent.parent === intvec)
    @test (A = Adjoint(Adjoint(intmat))::Adjoint{Int,Adjoint{Int,Matrix{Int}}}; A.parent.parent === intmat)
    @test (A = Transpose(Transpose(intvec))::Transpose{Int,Transpose{Int,Vector{Int}}}; A.parent.parent === intvec)
    @test (A = Transpose(Transpose(intmat))::Transpose{Int,Transpose{Int,Matrix{Int}}}; A.parent.parent === intmat)
end

@testset "Adjoint and Transpose basic AbstractArray functionality" begin
    # vectors and matrices with real scalar eltype, and their adjoints/transposes
    intvec, intmat = [1, 2], [1 2 3; 4 5 6]
    tintvec, tintmat = [1 2], [1 4; 2 5; 3 6]
    @testset "length methods" begin
        @test length(Adjoint(intvec)) == length(intvec)
        @test length(Adjoint(intmat)) == length(intmat)
        @test length(Transpose(intvec)) == length(intvec)
        @test length(Transpose(intmat)) == length(intmat)
    end
    @testset "size methods" begin
        @test size(Adjoint(intvec)) == (1, length(intvec))
        @test size(Adjoint(intmat)) == reverse(size(intmat))
        @test size(Transpose(intvec)) == (1, length(intvec))
        @test size(Transpose(intmat)) == reverse(size(intmat))
    end
    @testset "indices methods" begin
        @test axes(Adjoint(intvec)) == (Base.OneTo(1), Base.OneTo(length(intvec)))
        @test axes(Adjoint(intmat)) == reverse(axes(intmat))
        @test axes(Transpose(intvec)) == (Base.OneTo(1), Base.OneTo(length(intvec)))
        @test axes(Transpose(intmat)) == reverse(axes(intmat))
    end
    @testset "IndexStyle methods" begin
        @test IndexStyle(Adjoint(intvec)) == IndexLinear()
        @test IndexStyle(Adjoint(intmat)) == IndexCartesian()
        @test IndexStyle(Transpose(intvec)) == IndexLinear()
        @test IndexStyle(Transpose(intmat)) == IndexCartesian()
    end
    # vectors and matrices with complex scalar eltype, and their adjoints/transposes
    complexintvec, complexintmat = [1im, 2im], [1im 2im 3im; 4im 5im 6im]
    tcomplexintvec, tcomplexintmat = [1im 2im], [1im 4im; 2im 5im; 3im 6im]
    acomplexintvec, acomplexintmat = conj.(tcomplexintvec), conj.(tcomplexintmat)
    # vectors and matrices with real-vector and real-matrix eltype, and their adjoints/transposes
    intvecvec = [[1, 2], [3, 4]]
    tintvecvec = [[[1 2]] [[3 4]]]
    intmatmat = [[[1 2]] [[3  4]] [[ 5  6]];
                 [[7 8]] [[9 10]] [[11 12]]]
    tintmatmat = [[hcat([1, 2])] [hcat([7, 8])];
                  [hcat([3, 4])] [hcat([9, 10])];
                  [hcat([5, 6])] [hcat([11, 12])]]
    # vectors and matrices with complex-vector and complex-matrix eltype, and their adjoints/transposes
    complexintvecvec, complexintmatmat = im .* (intvecvec, intmatmat)
    tcomplexintvecvec, tcomplexintmatmat = im .* (tintvecvec, tintmatmat)
    acomplexintvecvec, acomplexintmatmat = conj.(tcomplexintvecvec), conj.(tcomplexintmatmat)
    @testset "getindex methods, elementary" begin
        # implicitly test elementary definitions, for arrays with concrete real scalar eltype
        @test Adjoint(intvec) == tintvec
        @test Adjoint(intmat) == tintmat
        @test Transpose(intvec) == tintvec
        @test Transpose(intmat) == tintmat
        # implicitly test elementary definitions, for arrays with concrete complex scalar eltype
        @test Adjoint(complexintvec) == acomplexintvec
        @test Adjoint(complexintmat) == acomplexintmat
        @test Transpose(complexintvec) == tcomplexintvec
        @test Transpose(complexintmat) == tcomplexintmat
        # implicitly test elementary definitions, for arrays with concrete real-array eltype
        @test Adjoint(intvecvec) == tintvecvec
        @test Adjoint(intmatmat) == tintmatmat
        @test Transpose(intvecvec) == tintvecvec
        @test Transpose(intmatmat) == tintmatmat
        # implicitly test elementary definitions, for arrays with concrete complex-array type
        @test Adjoint(complexintvecvec) == acomplexintvecvec
        @test Adjoint(complexintmatmat) == acomplexintmatmat
        @test Transpose(complexintvecvec) == tcomplexintvecvec
        @test Transpose(complexintmatmat) == tcomplexintmatmat
    end
    @testset "getindex(::AdjOrTransVec, ::Colon, ::AbstractArray{Int}) methods that preserve wrapper type" begin
        # for arrays with concrete scalar eltype
        @test Adjoint(intvec)[:, [1, 2]] == Adjoint(intvec)
        @test Transpose(intvec)[:, [1, 2]] == Transpose(intvec)
        @test Adjoint(complexintvec)[:, [1, 2]] == Adjoint(complexintvec)
        @test Transpose(complexintvec)[:, [1, 2]] == Transpose(complexintvec)
        # for arrays with concrete array eltype
        @test Adjoint(intvecvec)[:, [1, 2]] == Adjoint(intvecvec)
        @test Transpose(intvecvec)[:, [1, 2]] == Transpose(intvecvec)
        @test Adjoint(complexintvecvec)[:, [1, 2]] == Adjoint(complexintvecvec)
        @test Transpose(complexintvecvec)[:, [1, 2]] == Transpose(complexintvecvec)
    end
    @testset "getindex(::AdjOrTransVec, ::Colon, ::Colon) methods that preserve wrapper type" begin
        # for arrays with concrete scalar eltype
        @test Adjoint(intvec)[:, :] == Adjoint(intvec)
        @test Transpose(intvec)[:, :] == Transpose(intvec)
        @test Adjoint(complexintvec)[:, :] == Adjoint(complexintvec)
        @test Transpose(complexintvec)[:, :] == Transpose(complexintvec)
        # for arrays with concrete array elype
        @test Adjoint(intvecvec)[:, :] == Adjoint(intvecvec)
        @test Transpose(intvecvec)[:, :] == Transpose(intvecvec)
        @test Adjoint(complexintvecvec)[:, :] == Adjoint(complexintvecvec)
        @test Transpose(complexintvecvec)[:, :] == Transpose(complexintvecvec)
    end
    @testset "getindex(::AdjOrTransVec, ::Colon, ::Int) should preserve wrapper type on result entries" begin
        # for arrays with concrete scalar eltype
        @test Adjoint(intvec)[:, 2] == intvec[2:2]
        @test Transpose(intvec)[:, 2] == intvec[2:2]
        @test Adjoint(complexintvec)[:, 2] == conj.(complexintvec[2:2])
        @test Transpose(complexintvec)[:, 2] == complexintvec[2:2]
        # for arrays with concrete array eltype
        @test Adjoint(intvecvec)[:, 2] == Adjoint.(intvecvec[2:2])
        @test Transpose(intvecvec)[:, 2] == Transpose.(intvecvec[2:2])
        @test Adjoint(complexintvecvec)[:, 2] == Adjoint.(complexintvecvec[2:2])
        @test Transpose(complexintvecvec)[:, 2] == Transpose.(complexintvecvec[2:2])
    end
    @testset "setindex! methods" begin
        # for vectors with real scalar eltype
        @test (wv = Adjoint(copy(intvec));
                wv === setindex!(wv, 3, 2) &&
                 wv == setindex!(copy(tintvec), 3, 1, 2)    )
        @test (wv = Transpose(copy(intvec));
                wv === setindex!(wv, 4, 2) &&
                 wv == setindex!(copy(tintvec), 4, 1, 2)    )
        # for matrices with real scalar eltype
        @test (wA = Adjoint(copy(intmat));
                wA === setindex!(wA, 7, 3, 1) &&
                 wA == setindex!(copy(tintmat), 7, 3, 1)    )
        @test (wA = Transpose(copy(intmat));
                wA === setindex!(wA, 7, 3, 1) &&
                 wA == setindex!(copy(tintmat), 7, 3, 1)    )
        # for vectors with complex scalar eltype
        @test (wz = Adjoint(copy(complexintvec));
                wz === setindex!(wz, 3im, 2) &&
                 wz == setindex!(copy(acomplexintvec), 3im, 1, 2)   )
        @test (wz = Transpose(copy(complexintvec));
                wz === setindex!(wz, 4im, 2) &&
                 wz == setindex!(copy(tcomplexintvec), 4im, 1, 2)   )
        # for  matrices with complex scalar eltype
        @test (wZ = Adjoint(copy(complexintmat));
                wZ === setindex!(wZ, 7im, 3, 1) &&
                 wZ == setindex!(copy(acomplexintmat), 7im, 3, 1)   )
        @test (wZ = Transpose(copy(complexintmat));
                wZ === setindex!(wZ, 7im, 3, 1) &&
                 wZ == setindex!(copy(tcomplexintmat), 7im, 3, 1)   )
        # for vectors with concrete real-vector eltype
        @test (wv = Adjoint(copy(intvecvec));
                wv === setindex!(wv, Adjoint([5, 6]), 2) &&
                 wv == setindex!(copy(tintvecvec), [5 6], 2))
        @test (wv = Transpose(copy(intvecvec));
                wv === setindex!(wv, Transpose([5, 6]), 2) &&
                 wv == setindex!(copy(tintvecvec), [5 6], 2))
        # for matrices with concrete real-matrix eltype
        @test (wA = Adjoint(copy(intmatmat));
                wA === setindex!(wA, Adjoint([13 14]), 3, 1) &&
                 wA == setindex!(copy(tintmatmat), hcat([13, 14]), 3, 1))
        @test (wA = Transpose(copy(intmatmat));
                wA === setindex!(wA, Transpose([13 14]), 3, 1) &&
                 wA == setindex!(copy(tintmatmat), hcat([13, 14]), 3, 1))
        # for vectors with concrete complex-vector eltype
        @test (wz = Adjoint(copy(complexintvecvec));
                wz === setindex!(wz, Adjoint([5im, 6im]), 2) &&
                 wz == setindex!(copy(acomplexintvecvec), [-5im -6im], 2))
        @test (wz = Transpose(copy(complexintvecvec));
                wz === setindex!(wz, Transpose([5im, 6im]), 2) &&
                 wz == setindex!(copy(tcomplexintvecvec), [5im 6im], 2))
        # for matrices with concrete complex-matrix eltype
        @test (wZ = Adjoint(copy(complexintmatmat));
                wZ === setindex!(wZ, Adjoint([13im 14im]), 3, 1) &&
                 wZ == setindex!(copy(acomplexintmatmat), hcat([-13im, -14im]), 3, 1))
        @test (wZ = Transpose(copy(complexintmatmat));
                wZ === setindex!(wZ, Transpose([13im 14im]), 3, 1) &&
                 wZ == setindex!(copy(tcomplexintmatmat), hcat([13im, 14im]), 3, 1))
    end
end

@testset "Adjoint and Transpose convert methods that convert underlying storage" begin
    intvec, intmat = [1, 2], [1 2 3; 4 5 6]
    @test convert(Adjoint{Float64,Vector{Float64}}, Adjoint(intvec))::Adjoint{Float64,Vector{Float64}} == Adjoint(intvec)
    @test convert(Adjoint{Float64,Matrix{Float64}}, Adjoint(intmat))::Adjoint{Float64,Matrix{Float64}} == Adjoint(intmat)
    @test convert(Transpose{Float64,Vector{Float64}}, Transpose(intvec))::Transpose{Float64,Vector{Float64}} == Transpose(intvec)
    @test convert(Transpose{Float64,Matrix{Float64}}, Transpose(intmat))::Transpose{Float64,Matrix{Float64}} == Transpose(intmat)
end

@testset "Adjoint and Transpose similar methods" begin
    intvec, intmat = [1, 2], [1 2 3; 4 5 6]
    # similar with no additional specifications, vector (rewrapping) semantics
    @test size(similar(Adjoint(intvec))::Adjoint{Int,Vector{Int}}) == size(Adjoint(intvec))
    @test size(similar(Transpose(intvec))::Transpose{Int,Vector{Int}}) == size(Transpose(intvec))
    # similar with no additional specifications, matrix (no-rewrapping) semantics
    @test size(similar(Adjoint(intmat))::Matrix{Int}) == size(Adjoint(intmat))
    @test size(similar(Transpose(intmat))::Matrix{Int}) == size(Transpose(intmat))
    # similar with element type specification, vector (rewrapping) semantics
    @test size(similar(Adjoint(intvec), Float64)::Adjoint{Float64,Vector{Float64}}) == size(Adjoint(intvec))
    @test size(similar(Transpose(intvec), Float64)::Transpose{Float64,Vector{Float64}}) == size(Transpose(intvec))
    # similar with element type specification, matrix (no-rewrapping) semantics
    @test size(similar(Adjoint(intmat), Float64)::Matrix{Float64}) == size(Adjoint(intmat))
    @test size(similar(Transpose(intmat), Float64)::Matrix{Float64}) == size(Transpose(intmat))
    # similar with element type and arbitrary dims specifications
    shape = (2, 2, 2)
    @test size(similar(Adjoint(intvec), Float64, shape)::Array{Float64,3}) == shape
    @test size(similar(Adjoint(intmat), Float64, shape)::Array{Float64,3}) == shape
    @test size(similar(Transpose(intvec), Float64, shape)::Array{Float64,3}) == shape
    @test size(similar(Transpose(intmat), Float64, shape)::Array{Float64,3}) == shape
end

@testset "Adjoint and Transpose parent methods" begin
    intvec, intmat = [1, 2], [1 2 3; 4 5 6]
    @test parent(Adjoint(intvec)) === intvec
    @test parent(Adjoint(intmat)) === intmat
    @test parent(Transpose(intvec)) === intvec
    @test parent(Transpose(intmat)) === intmat
end

@testset "Adjoint and Transpose vector vec methods" begin
    intvec = [1, 2]
    @test vec(Adjoint(intvec)) == intvec
    @test vec(Transpose(intvec)) === intvec
    cvec = [1 + 1im]
    @test vec(cvec')[1] == cvec[1]'
end

@testset "horizontal concatenation of Adjoint/Transpose-wrapped vectors and Numbers" begin
    # horizontal concatenation of Adjoint/Transpose-wrapped vectors and Numbers
    # should preserve the Adjoint/Transpose-wrapper to preserve semantics downstream
    vec, tvec, avec = [1im, 2im], [1im 2im], [-1im -2im]
    vecvec = [[1im, 2im], [3im, 4im]]
    tvecvec = [[[1im 2im]] [[3im 4im]]]
    avecvec = [[[-1im -2im]] [[-3im -4im]]]
    # for arrays with concrete scalar eltype
    @test hcat(Adjoint(vec), Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == hcat(avec, avec)
    @test hcat(Adjoint(vec), 1, Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == hcat(avec, 1, avec)
    @test hcat(Transpose(vec), Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == hcat(tvec, tvec)
    @test hcat(Transpose(vec), 1, Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == hcat(tvec, 1, tvec)
    # for arrays with concrete array eltype
    @test hcat(Adjoint(vecvec), Adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == hcat(avecvec, avecvec)
    @test hcat(Transpose(vecvec), Transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == hcat(tvecvec, tvecvec)
end

@testset "map/broadcast over Adjoint/Transpose-wrapped vectors and Numbers" begin
    # map and broadcast over Adjoint/Transpose-wrapped vectors and Numbers
    # should preserve the Adjoint/Transpose-wrapper to preserve semantics downstream
    vec, tvec, avec = [1im, 2im], [1im 2im], [-1im -2im]
    vecvec = [[1im, 2im], [3im, 4im]]
    tvecvec = [[[1im 2im]] [[3im 4im]]]
    avecvec = [[[-1im -2im]] [[-3im -4im]]]
    # unary map over wrapped vectors with concrete scalar eltype
    @test map(-, Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == -avec
    @test map(-, Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == -tvec
    # unary map over wrapped vectors with concrete array eltype
    @test map(-, Adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -avecvec
    @test map(-, Transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -tvecvec
    # binary map over wrapped vectors with concrete scalar eltype
    @test map(+, Adjoint(vec), Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec
    @test map(+, Transpose(vec), Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec
    # binary map over wrapped vectors with concrete array eltype
    @test map(+, Adjoint(vecvec), Adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == avecvec + avecvec
    @test map(+, Transpose(vecvec), Transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == tvecvec + tvecvec
    # unary broadcast over wrapped vectors with concrete scalar eltype
    @test broadcast(-, Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == -avec
    @test broadcast(-, Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == -tvec
    # unary broadcast over wrapped vectors with concrete array eltype
    @test broadcast(-, Adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -avecvec
    @test broadcast(-, Transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == -tvecvec
    # binary broadcast over wrapped vectors with concrete scalar eltype
    @test broadcast(+, Adjoint(vec), Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec
    @test broadcast(+, Transpose(vec), Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec
    # binary broadcast over wrapped vectors with concrete array eltype
    @test broadcast(+, Adjoint(vecvec), Adjoint(vecvec))::Adjoint{Adjoint{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == avecvec + avecvec
    @test broadcast(+, Transpose(vecvec), Transpose(vecvec))::Transpose{Transpose{Complex{Int},Vector{Complex{Int}}},Vector{Vector{Complex{Int}}}} == tvecvec + tvecvec
    # trinary broadcast over wrapped vectors with concrete scalar eltype and numbers
    @test broadcast(+, Adjoint(vec), 1, Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec .+ 1
    @test broadcast(+, Transpose(vec), 1, Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec .+ 1
    @test broadcast(+, Adjoint(vec), 1im, Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec .+ 1im
    @test broadcast(+, Transpose(vec), 1im, Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec .+ 1im
    # ascertain inference friendliness, ref. https://github.com/JuliaLang/julia/pull/25083#issuecomment-353031641
    sparsevec = SparseVector([1.0, 2.0, 3.0])
    @test map(-, Adjoint(sparsevec), Adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}}
    @test map(-, Transpose(sparsevec), Transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}}
    @test broadcast(-, Adjoint(sparsevec), Adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}}
    @test broadcast(-, Transpose(sparsevec), Transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}}
    @test broadcast(+, Adjoint(sparsevec), 1.0, Adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}}
    @test broadcast(+, Transpose(sparsevec), 1.0, Transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}}
end

@testset "Adjoint/Transpose-wrapped vector multiplication" begin
    realvec, realmat = [1, 2, 3], [1 2 3; 4 5 6; 7 8 9]
    complexvec, complexmat = [1im, 2, -3im], [1im 2 3; 4 5 -6im; 7im 8 9]
    # Adjoint/Transpose-vector * vector
    @test Adjoint(realvec) * realvec == dot(realvec, realvec)
    @test Transpose(realvec) * realvec == dot(realvec, realvec)
    @test Adjoint(complexvec) * complexvec == dot(complexvec, complexvec)
    @test Transpose(complexvec) * complexvec == dot(conj(complexvec), complexvec)
    # vector * Adjoint/Transpose-vector
    @test realvec * Adjoint(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3)))
    @test realvec * Transpose(realvec) == broadcast(*, realvec, reshape(realvec, (1, 3)))
    @test complexvec * Adjoint(complexvec) == broadcast(*, complexvec, reshape(conj(complexvec), (1, 3)))
    @test complexvec * Transpose(complexvec) == broadcast(*, complexvec, reshape(complexvec, (1, 3)))
    # Adjoint/Transpose-vector * matrix
    @test (Adjoint(realvec) * realmat)::Adjoint{Int,Vector{Int}} ==
        reshape(copy(Adjoint(realmat)) * realvec, (1, 3))
    @test (Transpose(realvec) * realmat)::Transpose{Int,Vector{Int}} ==
        reshape(copy(Transpose(realmat)) * realvec, (1, 3))
    @test (Adjoint(complexvec) * complexmat)::Adjoint{Complex{Int},Vector{Complex{Int}}} ==
        reshape(conj(copy(Adjoint(complexmat)) * complexvec), (1, 3))
    @test (Transpose(complexvec) * complexmat)::Transpose{Complex{Int},Vector{Complex{Int}}} ==
        reshape(copy(Transpose(complexmat)) * complexvec, (1, 3))
    # Adjoint/Transpose-vector * Adjoint/Transpose-matrix
    @test (Adjoint(realvec) * Adjoint(realmat))::Adjoint{Int,Vector{Int}} ==
        reshape(realmat * realvec, (1, 3))
    @test (Transpose(realvec) * Transpose(realmat))::Transpose{Int,Vector{Int}} ==
        reshape(realmat * realvec, (1, 3))
    @test (Adjoint(complexvec) * Adjoint(complexmat))::Adjoint{Complex{Int},Vector{Complex{Int}}} ==
        reshape(conj(complexmat * complexvec), (1, 3))
    @test (Transpose(complexvec) * Transpose(complexmat))::Transpose{Complex{Int},Vector{Complex{Int}}} ==
        reshape(complexmat * complexvec, (1, 3))
end

@testset "Adjoint/Transpose-wrapped vector pseudoinversion" begin
    realvec, complexvec = [1, 2, 3, 4], [1im, 2, 3im, 4]
    rowrealvec, rowcomplexvec = reshape(realvec, (1, 4)), reshape(complexvec, (1, 4))
    # pinv(Adjoint/Transpose-vector) should match matrix equivalents
    # TODO tighten type asserts once pinv yields Transpose/Adjoint
    @test pinv(Adjoint(realvec))::Vector{Float64} ≈ pinv(rowrealvec)
    @test pinv(Transpose(realvec))::Vector{Float64} ≈ pinv(rowrealvec)
    @test pinv(Adjoint(complexvec))::Vector{Complex{Float64}} ≈ pinv(conj(rowcomplexvec))
    @test pinv(Transpose(complexvec))::Vector{Complex{Float64}} ≈ pinv(rowcomplexvec)
end

@testset "Adjoint/Transpose-wrapped vector left-division" begin
    realvec, complexvec = [1., 2., 3., 4.,], [1.0im, 2., 3.0im, 4.]
    rowrealvec, rowcomplexvec = reshape(realvec, (1, 4)), reshape(complexvec, (1, 4))
    # \(Adjoint/Transpose-vector, Adjoint/Transpose-vector) should mat matrix equivalents
    @test Adjoint(realvec)\Adjoint(realvec) ≈ rowrealvec\rowrealvec
    @test Transpose(realvec)\Transpose(realvec) ≈ rowrealvec\rowrealvec
    @test Adjoint(complexvec)\Adjoint(complexvec) ≈ conj(rowcomplexvec)\conj(rowcomplexvec)
    @test Transpose(complexvec)\Transpose(complexvec) ≈ rowcomplexvec\rowcomplexvec
end

@testset "Adjoint/Transpose-wrapped vector right-division" begin
    realvec, realmat = [1, 2, 3], [1 0 0; 0 2 0; 0 0 3]
    complexvec, complexmat = [1im, 2, -3im], [2im 0 0; 0 3 0; 0 0 -5im]
    rowrealvec, rowcomplexvec = reshape(realvec, (1, 3)), reshape(complexvec, (1, 3))
    # /(Adjoint/Transpose-vector, matrix)
    @test (Adjoint(realvec) / realmat)::Adjoint ≈ rowrealvec / realmat
    @test (Adjoint(complexvec) / complexmat)::Adjoint ≈ conj(rowcomplexvec) / complexmat
    @test (Transpose(realvec) / realmat)::Transpose ≈ rowrealvec / realmat
    @test (Transpose(complexvec) / complexmat)::Transpose ≈ rowcomplexvec / complexmat
    # /(Adjoint/Transpose-vector, Adjoint matrix)
    @test (Adjoint(realvec) / Adjoint(realmat))::Adjoint ≈ rowrealvec / copy(Adjoint(realmat))
    @test (Adjoint(complexvec) / Adjoint(complexmat))::Adjoint ≈ conj(rowcomplexvec) / copy(Adjoint(complexmat))
    @test (Transpose(realvec) / Adjoint(realmat))::Transpose ≈ rowrealvec / copy(Adjoint(realmat))
    @test (Transpose(complexvec) / Adjoint(complexmat))::Transpose ≈ rowcomplexvec / copy(Adjoint(complexmat))
    # /(Adjoint/Transpose-vector, Transpose matrix)
    @test (Adjoint(realvec) / Transpose(realmat))::Adjoint ≈ rowrealvec / copy(Transpose(realmat))
    @test (Adjoint(complexvec) / Transpose(complexmat))::Adjoint ≈ conj(rowcomplexvec) / copy(Transpose(complexmat))
    @test (Transpose(realvec) / Transpose(realmat))::Transpose ≈ rowrealvec / copy(Transpose(realmat))
    @test (Transpose(complexvec) / Transpose(complexmat))::Transpose ≈ rowcomplexvec / copy(Transpose(complexmat))
end

@testset "norm and opnorm of Adjoint/Transpose-wrapped vectors" begin
    # definitions are in base/linalg/generic.jl
    realvec, complexvec = [3, -4], [3im, -4im]
    # one norm result should be sum(abs.(realvec)) == 7
    # two norm result should be sqrt(sum(abs.(realvec))) == 5
    # inf norm result should be maximum(abs.(realvec)) == 4
    for v in (realvec, complexvec)
        @test norm(Adjoint(v)) ≈ 5
        @test norm(Adjoint(v), 1) ≈ 7
        @test norm(Adjoint(v), Inf) ≈ 4
        @test norm(Transpose(v)) ≈ 5
        @test norm(Transpose(v), 1) ≈ 7
        @test norm(Transpose(v), Inf) ≈ 4
    end
    # one opnorm result should be maximum(abs.(realvec)) == 4
    # two opnorm result should be sqrt(sum(abs.(realvec))) == 5
    # inf opnorm result should be sum(abs.(realvec)) == 7
    for v in (realvec, complexvec)
        @test opnorm(Adjoint(v)) ≈ 5
        @test opnorm(Adjoint(v), 1) ≈ 4
        @test opnorm(Adjoint(v), Inf) ≈ 7
        @test opnorm(Transpose(v)) ≈ 5
        @test opnorm(Transpose(v), 1) ≈ 4
        @test opnorm(Transpose(v), Inf) ≈ 7
    end
end

@testset "adjoint and transpose of Numbers" begin
    @test adjoint(1) == 1
    @test adjoint(1.0) == 1.0
    @test adjoint(1im) == -1im
    @test adjoint(1.0im) == -1.0im
    @test transpose(1) == 1
    @test transpose(1.0) == 1.0
    @test transpose(1im) == 1im
    @test transpose(1.0im) == 1.0im
end

@testset "adjoint!(a, b) return a" begin
    a = fill(1.0+im, 5)
    b = fill(1.0+im, 1, 5)
    @test adjoint!(a, b) === a
    @test adjoint!(b, a) === b
end

@testset "aliasing with adjoint and transpose" begin
    A = collect(reshape(1:25, 5, 5)) .+ rand.().*im
    B = copy(A)
    B .= B'
    @test B == A'
    B = copy(A)
    B .= transpose(B)
    @test B == transpose(A)
    B = copy(A)
    B .= B .* B'
    @test B == A .* A'
end

@testset "test show methods for $t of Factorizations" for t in (Adjoint, Transpose)
    A = randn(4, 4)
    F = lu(A)
    Fop = t(F)
    @test "LinearAlgebra."*sprint(show, Fop) ==
                  "$t of "*sprint(show, parent(Fop))
    @test "LinearAlgebra."*sprint((io, t) -> show(io, MIME"text/plain"(), t), Fop) ==
                  "$t of "*sprint((io, t) -> show(io, MIME"text/plain"(), t), parent(Fop))
end

const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl"))
using .Main.OffsetArrays

@testset "offset axes" begin
    s = Base.Slice(-3:3)'
    @test axes(s) === (Base.OneTo(1), Base.IdentityUnitRange(-3:3))
    @test collect(LinearIndices(s)) == reshape(1:7, 1, 7)
    @test collect(CartesianIndices(s)) == reshape([CartesianIndex(1,i) for i = -3:3], 1, 7)
    @test s[1] == -3
    @test s[7] ==  3
    @test s[4] ==  0
    @test_throws BoundsError s[0]
    @test_throws BoundsError s[8]
    @test s[1,-3] == -3
    @test s[1, 3] ==  3
    @test s[1, 0] ==  0
    @test_throws BoundsError s[1,-4]
    @test_throws BoundsError s[1, 4]
end

end # module TestAdjointTranspose

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
function nsieve(limit::Int)
    primes = trues(limit)
    primes[1] = false
    p = 2
    while !isnothing(p) && p <= isqrt(limit)
        primes[p*p:p:limit] .= false
        p = findnext(primes, p + 1)
    end
    println("Primes up to ", lpad(limit, 8), lpad(count(primes), 9))
end

function main()
    n = isempty(ARGS) ? 4 : parse(Int, ARGS[1])
    for i in n:-1:n-2
        nsieve(10000 << i)
    end
end

function real_main()
    isinteractive() || main()
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end


text:
export zoo_jlc, zoo_jlc_file,
    JlcFormat, JlcVtreeFormat

using DirectedAcyclicGraphs: label_nodes

struct JlcFormat <: FileFormat end

const JlcVtreeFormat = Tuple{JlcFormat,VtreeFormat}
Tuple{JlcFormat,VtreeFormat}() = (JlcFormat(),VtreeFormat())

##############################################
# Read JLC (Juice Logic Circuit)
##############################################

zoo_jlc_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/jlcs/$name"

"""
    zoo_jlc(name)

Loads JLC file with given name from model zoo. See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
zoo_jlc(name) =
    read(zoo_jlc_file(name), LogicCircuit, JlcFormat())

const jlc_grammar = raw"""
    start: header (_NL node)+ _NL?

    header : "jlc" _WS INT

    node : "F" _WS INT -> false_node
         | "T" _WS INT -> true_node
         | "L" _WS INT _WS INT _WS SIGNED_INT -> literal_node
         | "O" _WS INT _WS INT _WS INT _WS child_nodes -> or_node
         | "A" _WS INT _WS INT _WS INT _WS child_nodes -> and_node

    child_nodes : INT (_WS INT)*

    %import common.INT
    %import common.SIGNED_INT
    %import common.WS_INLINE -> _WS
    %import common.NEWLINE -> _NL
    """ * dimacs_comments

const jlc_parser = Lark(jlc_grammar)

abstract type JlcParse <: JuiceTransformer end

@inline_rule header(t::JlcParse, x) =
    Base.parse(Int,x)

@rule start(t::JlcParse, x) = begin
    @assert num_nodes(x[end]) == x[1]
    x[end]
end

@rule child_nodes(t::JlcParse, x) =
    map(id -> t.nodes[id], x)

#  parse unstructured
struct PlainJlcParse <: JlcParse
    nodes::Dict{String,PlainLogicCircuit}
    PlainJlcParse() = new(Dict{String,PlainLogicCircuit}())
end

@rule literal_node(t::PlainJlcParse, x) =
    t.nodes[x[1]] = PlainLiteralNode(Base.parse(Lit,x[3]))

@rule false_node(t::PlainJlcParse, x) =
    t.nodes[x[1]] = PlainConstantNode(false)

@rule true_node(t::PlainJlcParse, x) =
    t.nodes[x[1]] = PlainConstantNode(true)

@rule or_node(t::PlainJlcParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    t.nodes[x[1]] = Plain⋁Node(x[4])
end

@rule and_node(t::PlainJlcParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    t.nodes[x[1]] = Plain⋀Node(x[4])
end

function Base.parse(::Type{PlainLogicCircuit}, str, ::JlcFormat)
    ast = Lerche.parse(jlc_parser, str)
    Lerche.transform(PlainJlcParse(), ast)
end

Base.read(io::IO, ::Type{PlainLogicCircuit}, ::JlcFormat) =
    parse(PlainLogicCircuit, read(io, String), JlcFormat())

#  parse structured
struct StructJlcParse <: JlcParse
    id2vtree::Dict{String,<:Vtree}
    nodes::Dict{String,StructLogicCircuit}
    StructJlcParse(id2vtree) =
        new(id2vtree,Dict{String,StructLogicCircuit}())
end

@rule literal_node(t::StructJlcParse, x) = begin
    lit = Base.parse(Lit,x[3])
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = PlainStructLiteralNode(lit, vtree)
end

@rule false_node(t::StructJlcParse, x) =
    t.nodes[x[1]] = PlainStructConstantNode(false)

@rule true_node(t::StructJlcParse, x) =
    t.nodes[x[1]] = PlainStructConstantNode(true)

@rule or_node(t::StructJlcParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = PlainStruct⋁Node(x[4], vtree)
end

@rule and_node(t::StructJlcParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3]) == 2
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = PlainStruct⋀Node(x[4]..., vtree)
end

function Base.parse(::Type{PlainStructLogicCircuit}, str::AbstractString, ::JlcFormat, id2vtree)
    ast = Lerche.parse(jlc_parser, str)
    Lerche.transform(StructJlcParse(id2vtree), ast)
end

function Base.parse(::Type{PlainStructLogicCircuit}, strings, format::JlcVtreeFormat)
    id2vtree = parse(Dict{String,Vtree}, strings[2], format[2])
    parse(PlainStructLogicCircuit, strings[1], format[1], id2vtree)
end

Base.read(io::IO, ::Type{PlainStructLogicCircuit}, ::JlcFormat, id2vtree) =
    parse(PlainStructLogicCircuit, read(io, String), JlcFormat(), id2vtree)

function Base.read(ios::Tuple{IO,IO}, ::Type{PlainStructLogicCircuit}, ::JlcVtreeFormat)
    circuit_str = read(ios[1], String)
    vtree_str = read(ios[2], String)
    parse(PlainStructLogicCircuit, (circuit_str,vtree_str), JlcVtreeFormat())
end

# Note: no parse as SDD since we are not sure whether
# the original file satisfied the extra SDD properties

##############################################
# Write JLCs
##############################################

const JLC_FORMAT = """c this file was saved by LogicCircuits.jl
c ids of jlc nodes start at 0
c jlc nodes appear bottom-up, children before parents
c
c file syntax:
c jlc count-of-jlc-nodes
c F id-of-false-jlc-node
c T id-of-true-jlc-node
c L id-of-literal-jlc-node id-of-vtree literal
c O id-of-or-jlc-node id-of-vtree number-of-children {child-id}+
c A id-of-and-jlc-node id-of-vtree number-of-children {child-id}+
c"""

function Base.write(io::IO, circuit::LogicCircuit, ::JlcFormat, vtreeid::Function = (x -> 0))

    labeling = label_nodes(circuit)
    map!(x -> x-1, values(labeling)) # vtree nodes are 0-based indexed

    println(io, JLC_FORMAT)
    println(io, "jlc $(num_nodes(circuit))")
    foreach(circuit) do n
        if isliteralgate(n)
            println(io, "L $(labeling[n]) $(vtreeid(n)) $(literal(n))")
        elseif isconstantgate(n)
            sign = isfalse(n) ? "F" : "T"
            println(io, "$sign $(labeling[n])")
        else
            t = is⋀gate(n) ? "A" : "O"
            print(io, "$t $(labeling[n]) $(vtreeid(n)) $(num_children(n))")
            for child in children(n)
                print(io, " $(labeling[child])")
            end
            println(io)
        end
    end
    nothing
end

function Base.write(ios::Tuple{IO,IO}, circuit::StructLogicCircuit, format::JlcVtreeFormat)
    vtree2id = write(ios[2], vtree(circuit), format[2])
    write(ios[1], circuit, format[1], n -> vtree2id[vtree(n)])
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

# # Cylinder in torsion using distributed coupling

# Distributed coupling can be used to distribute point moments as point forces
# to nodes which dont have rotational degrees of freedom.

# The model is a 3d cylinder, shown in picture.

# ![](example_dcoupling_cylindrical_beam/example_dcoupling_cylindrical_beam.png)

# Distributed coupling is used to distribute a torque load to the end of the
# cylinder. The other end is fixed. Node sets for coupling nodes, reference
# node and fixed face were made in the ABAQUS input file. Reference node is
# located in the centrum of the coupling nodes.

using JuliaFEM
using JuliaFEM: add_elements!, Problem
using JuliaFEM.Preprocess
using JuliaFEM.Abaqus: create_surface_elements
using FEMBase
using FEMCoupling
using FEMCoupling: add_reference_node!, add_coupling_nodes!

# reading mesh from ABAQUS input file

datadir = Pkg.dir("FEMCoupling", "examples", "example_dcoupling_cylindrical_beam")
mesh = abaqus_read_mesh(joinpath(datadir, "example_dcoupling_cylindrical_beam.inp"))
println("Number of nodes in a model: ", length(mesh.nodes))

# # Elements
# Creating elements for the body. Mesh and element types are defined in
# the ABAQUS input file. The cylinder body is named "Body1" in the input file.

cylinder_body = create_elements(mesh,"Body1")

# Updating values for the elements.

update!(cylinder_body, "youngs modulus", 210e3)
update!(cylinder_body, "poissons ratio", 0.3)

# Creating an elasticity problem and adding the elements to it.

cylinder_problem = Problem(Elasticity,"cylinder_problem",3)
add_elements!(cylinder_problem, cylinder_body)

# # Boundary conditions
# Creating Poi1-type elements as boundary condition elements to nodes of the
# node set Fixed_face_set.

bc_elements = [Element(Poi1, [j]) for j in mesh.node_sets[:Fixed_face_set]]

# Updating geometry for the bc elements

update!(bc_elements, "geometry", mesh.nodes)

# Fixing all displacements for the bc elements.

for i=1:3
    update!(bc_elements, "displacement $i", 0.0)
end

# Creating a bc problem and adding the bc elements to it.

bc = Problem(Dirichlet, "fixed", 3, "displacement")
add_elements!(bc, bc_elements)

# # Distributed coupling
# Creating Poi1 elements to nodes in coupling nodes set.

coupling_nodes = [Element(Poi1, [j]) for j in mesh.node_sets[:Coupling_nodes_set]]

# Updating geometry for the coupling nodes.

update!(coupling_nodes, "geometry", mesh.nodes)

# Creating Poi1 element for the reference node.

reference_node_id = collect(mesh.node_sets[:ref_node_set])
reference_node = Element(Poi1, reference_node_id)

# Updating geometry and applying a point moment for the reference node.

update!(reference_node, "geometry", mesh.nodes)
update!(reference_node, "point moment 3", 1500.0)

# Creating a coupling problem and adding coupling nodes and reference nodes to
# it.

coupling = Problem(Coupling, "cylind", 3, "displacement")
add_coupling_nodes!(coupling, coupling_nodes)
add_reference_node!(coupling, reference_node)

# # Analysis
# Creating a step and running the analysis. The cylinder_problem contains
# information about the body, its elements and their values. The bc contains
# information about boundary conditions and coupling contains information
# about distributed coupling.

step = Analysis(Nonlinear)
add_problems!(step, [cylinder_problem, bc, coupling])
run!(step)

# # Results

# Comparing calculated results with ABAQUS results. The node set
# circlenodes_set contains nodes which are on the outer face radius.
# These circle nodes should have the maximum displacement magnitude (norm(u)).

node_on_circle = first(mesh.node_sets[:circlenodes_set])

# declaring displacements at time 0.0 to variable u

time=0.0
u = cylinder_problem("displacement", time)[node_on_circle]
u_mag = norm(u)

# Making a testset.

using FEMBase.Test
@testset "displacement magnitude" begin
u_mag_expected=6.306e-4
@test isapprox(u_mag, u_mag_expected, rtol=1e-3)
end

# Printing node ids
println("reference node id = $(reference_node_id[1])")
println("node on circle id = $node_on_circle")


text:
using WinReg
using Base.Test

@test querykey(WinReg.HKEY_LOCAL_MACHINE,"System\\CurrentControlSet\\Control\\Session Manager\\Environment","OS") == "Windows_NT"


text:
using Gadfly

plot(layer(x=1:10, y=1:10, Geom.point, Theme(default_color=color("orange"))),
     layer(x=collect(1:10) .+ 0.1, y=collect(1:10) .+ 0.1, Geom.point,
           Theme(default_color=color("blue")), order=1))


text:
# default reset
reset!(π::AbstractPolicy) = π

mutable struct RandomPolicy <: AbstractPolicy end
action(policy::RandomPolicy, r, s′, A′) = rand(A′)

# include("online_gae.jl")
# include("actor_critic.jl")


text:
# StrICU.jl - Wrapper for ICU (International Components for Unicode) library

# Some content of the documentation strings was derived from the ICU header files.
# (Those portions copyright (C) 1996-2015, International Business Machines Corporation and others)

"""
    StrICU (International Components for Unicode) Wrapper
"""
module StrICU

using ModuleInterfaceTools

@api extend! StrBase

@static if !V6_COMPAT
    const is_windows = Sys.iswindows
    finalizer(o, f::Function) = Base.finalizer(f, o)
end

import Base: parse, get, close

export ICU
const ICU = StrICU

const cvt_utf8 = utf8
const cvt_utf16 = utf16
export cvt_utf8, cvt_utf16

const ByteStr   = Union{ASCIIStr, UTF8Str, String}

const WordStringCSE = Union{UCS2CSE, _UCS2CSE, UTF16CSE}
const WordStrings = Str{<:WordStringCSE}

export set_locale!

include("../deps/deps.jl")
include("../deps/versions.jl")

@static if is_windows()
    # make sure versions match
    v1 = last(matchall(r"\d{2}", iculib))
    v2 = last(matchall(r"\d{2}", iculibi18n))
    v1 == v2 ||
        error("ICU library version mismatch $v1 != $v2 -- please correct $(realpath("../deps/deps.jl"))")
end

function __init__()
    set_locale!("")
end

global version
global suffix

dliculib = Libdl.dlopen(iculib)

for (suf,ver) in [("",0);
                         [("_$i",i) for i in versions];
                         [("_$(string(i)[1])_$(string(i)[2])",i) for i in versions]]
    if Libdl.dlsym_e(dliculib, "u_strToUpper"*suf) != C_NULL
        @eval const version = $ver
        @eval const suffix  = $suf
        break
    end
end

_libicu(s, lib, p) = ( Symbol(string(p, s, suffix)), lib )

const UBool      = Int8
const UChar      = UInt16
const UErrorCode = Int32
const U_PARSE_CONTEXT_LEN = 16
const U_PARSE_TUPLE = ntuple((i)->0x20, U_PARSE_CONTEXT_LEN)

struct UParseError
    line::UErrorCode                                 # The line on which the error occured
    offset::UErrorCode                               # The character offset to the error
    preContext::NTuple{U_PARSE_CONTEXT_LEN, UInt8}   # Textual context before the error
    postContext::NTuple{U_PARSE_CONTEXT_LEN, UInt8}  # The error itself and/or textual context after the error
end
UParseError() = UParseError(0, 0, U_PARSE_TUPLE, U_PARSE_TUPLE)

FAILURE(x::Integer) = x > 0
SUCCESS(x::Integer) = x <= 0
U_BUFFER_OVERFLOW_ERROR = 15

const locale   = ASCIIStr[""]

include("utext.jl")
include("ustring.jl")
include("ubrk.jl")
include("ucnv.jl")
include("ucol.jl")
include("ucsdet.jl")
include("udat.jl")
include("ucal.jl")
include("ucasemap.jl")

@api freeze

end # module StrICU


text:
using ITensors, Printf

function main()
  i = Index(100,"i");

  A = randomITensor(i);
  B = randomITensor(i);

  # Dry run (JIT)
  time_dry = @elapsed begin
    R = A*B
  end
  reset!(timer)
  # End dry run

  time = @elapsed begin
    for n =1:10
    R = A*B
  end
  end

  @printf "time = %.12f\n" time
  printTimes(timer)

end; main()


text:
import Gadfly
Gadfly.set_default_plot_size(20Gadfly.cm, 12Gadfly.cm)

# make energy conservation and velocity distribution plots

pe = readdlm("pe.dat")
vir = readdlm("v.dat")

function mdplot(values::Array, title::String, name::String, N::Integer)
    p = Gadfly.plot(
        x=[10:10:N],
        y=values,
        Gadfly.Geom.line,
        Gadfly.Guide.xlabel("Simulation Step"),
        Gadfly.Guide.ylabel(name),
        Gadfly.Guide.title("$title vs. Simulation Step for Liquid Argon")
    )

    # save as svg and png
    Gadfly.draw(Gadfly.SVG("$name.svg", 20Gadfly.cm, 12Gadfly.cm), p)
    Gadfly.draw(Gadfly.PNG("$name.png", 12Gadfly.cm, 8Gadfly.cm), p)

    return p
end

# m = 1, k = 1 in our natural units
function maxboltzbuilder(T::Real)
    return v -> √(2/π*T^3) * v.^2 .* exp(-v.^2 / (2T))
end

function mdplotspeeds(speeds::Array, name::String, T::Real)
    xs = [0:0.01:1]
    p = Gadfly.plot(
        Gadfly.layer(
            x = speeds,
            Gadfly.Geom.histogram(bincount=100, density=true),
            color = ["Measured speed histogram for T = $T"]
        ),
###         Gadfly.layer(
###             x = xs,
###             y = maxboltzbuilder(T)(xs),
###             Gadfly.Geom.line,
###             color = ["Predicted by Maxwell-Boltzmann"]
###         ),
###         color = ["Legend"],
        Gadfly.Guide.title("Temperature Distribution for T=$T vs. Boltzmann Distribution"),
        Gadfly.Guide.xlabel("Velocity"),
        Gadfly.Guide.ylabel("Probability"),
        Gadfly.Guide.xticks(ticks=[0:0.1:1])
    )

    # save as svg and png
    Gadfly.draw(Gadfly.SVG("$name.svg", 20Gadfly.cm, 12Gadfly.cm), p)
    Gadfly.draw(Gadfly.PNG("$name.png", 12Gadfly.cm, 8Gadfly.cm), p)

    return p
end

pe_plots = Gadfly.Plot[]
push!(pe_plots, mdplot(pe, "Potential Energies (T = 1.069)", "temp-1.069", 20000))

e_plots = Gadfly.Plot[]
push!(e_plots, mdplot(e_1069, "Total Energies (T = 1.069)", "temp-1.069", 20000))


text:
using SafeTestsets

@safetestset "Updating Prior" begin
    using Test, Distributions, AdaptiveDesignOptimization

    α = 5
    β = 5
    prior = [Beta(α, β)]
    loglike(θ, d, data) = logpdf(Bernoulli(θ), data)
    model = Model(;prior, loglike)

    parm_list = (
        θ = range(0, 1, length=200),
    )

    design_list = (d = [1,2],)

    data_list = (choice=[true, false],)

    randomizer = Optimizer(;design_list, parm_list, data_list, model,
        approach=Randomize)

    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β)) atol = 5e-3

    update!(randomizer, false)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β+1)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β+1)) atol = 5e-3

    update!(randomizer, false)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β+2)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β+2)) atol = 5e-3

    update!(randomizer, true)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α+1,β+2)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α+1,β+2)) atol = 5e-3
end

@safetestset "Updating Default Prior" begin
    using Test, Distributions, AdaptiveDesignOptimization

    α = 1
    β = 1
    loglike(θ, d, data) = logpdf(Bernoulli(θ), data)
    model = Model(; loglike)

    parm_list = (
        θ = range(0, 1, length=200),
    )

    design_list = (d = [1,2],)

    data_list = (choice=[true, false],)

    randomizer = Optimizer(;design_list, parm_list, data_list, model,
        approach=Randomize)

    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β)) atol = 5e-3

    update!(randomizer, false)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β+1)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β+1)) atol = 5e-3

    update!(randomizer, false)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β+2)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β+2)) atol = 5e-3

    update!(randomizer, true)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α+1,β+2)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α+1,β+2)) atol = 5e-3
end

@safetestset "Updating Prior Dynamic Model" begin
    using Test, Distributions, AdaptiveDesignOptimization

    struct State1 end

    α = 5
    β = 5
    prior = [Beta(α, β)]
    loglike(θ, d, data, args...) = logpdf(Bernoulli(θ), data)
    model = Model(;prior, loglike)

    parm_list = (
        θ = range(0, 1, length=200),
    )

    design_list = (d = [1,2],)

    data_list = (choice=[true, false],)

    update_state!(args...) = nothing

    randomizer = Optimizer(;design_list, parm_list, data_list, model,
    update_state!, state_type=State1, model_type=Dynamic)

    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β)) atol = 5e-3

    update!(randomizer, false)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β+1)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β+1)) atol = 5e-3

    update!(randomizer, false)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α,β+2)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α,β+2)) atol = 5e-3

    update!(randomizer, true)
    @test mean_post(randomizer)[1] ≈ mean(Beta(α+1,β+2)) atol = 5e-3
    @test std_post(randomizer)[1] ≈ std(Beta(α+1,β+2)) atol = 5e-3
end

@safetestset "Dimension Check" begin
    using Test, Distributions, AdaptiveDesignOptimization

    loglike(μ, σ, data, args...) = logpdf(Normal(μ, σ), data)
    model = Model(; loglike)

    parm_list = (
        μ = range(0, 10, length=10),
        σ = range(1, 5, length=10),
    )

    design_list = (
        d = range(1, 5, length=10),
        z = range(1, 5, length=10),
        v = range(1, 5, length=10)
    )

    data_list = (
        choice=[true, false],
        v = 1:3
    )

    optimizer = Optimizer(;design_list, parm_list, data_list, model)

    @test optimizer.parm_grid |> size == (10,10)
    @test optimizer.design_grid |> size == (10,10,10)
    @test optimizer.data_grid |> size == (2,3)
    @test optimizer.log_like |> size == (100,1000,6)
    @test optimizer.log_post |> size == (100,)
    @test optimizer.priors |> size == (10,10)
    @test optimizer.mutual_info |> size == (1000,)
    @test optimizer.marg_log_like |> size == (1,1000,6)
    @test optimizer.marg_entropy |> size == (1000,)
    @test optimizer.cond_entropy |> size == (1000,)
end


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

import Test
Test.@testset "Package GERDAMetadata" begin

include("test_filekey.jl")

end # testset


max_stars_repo_name: Physolia/-THE-Magic-Button

max_issues_repo_name: Physolia/-The-System-

max_forks_repo_name: Physolia/-The-System-


text:
@testset "Maximization for quantum search" begin

    n = 5
    g = complete_graph(n)
    qss_penalty = ceil(Int, log(nv(g)))

    @testset "Continuous mode" begin

        cqss = QWSearch(CTQW(g), [1], qss_penalty)
        cqss0 = QWSearch(CTQW(g), [1])

        @test_throws AssertionError maximize_quantum_search(cqss, -1.)
        @test_throws AssertionError maximize_quantum_search(cqss, 1., -1.)
        @test runtime(maximize_quantum_search(cqss)) > 0
        @test runtime(maximize_quantum_search(cqss, 1., 1.)) > 0
        @test_warn r".*" maximize_quantum_search(cqss0)
        @test_nowarn maximize_quantum_search(cqss)

        @test isapprox(runtime(maximize_quantum_search(cqss, 1., 0.2)), 1., atol=0.001)

    end

    @testset "Discrete mode" begin
        dqss = QWSearch(Szegedy(g), [1], qss_penalty)
        dqss0 = QWSearch(Szegedy(g), [1])

        @test_throws AssertionError maximize_quantum_search(dqss0, -1)
        @test_throws AssertionError maximize_quantum_search(dqss0, 1, :unknown_mode)
        @test runtime(maximize_quantum_search(dqss)) > 0
        @test runtime(maximize_quantum_search(dqss, nv(g), :firstmaxprob)) > 0
        @test runtime(maximize_quantum_search(dqss, nv(g), :firstmaxeff)) > 0
        @test runtime(maximize_quantum_search(dqss, nv(g), :maxtimeeff)) > 0
        @test runtime(maximize_quantum_search(dqss, nv(g), :maxeff)) > 0
        @test runtime(maximize_quantum_search(dqss, nv(g), :maxtimeprob)) > 0

        @test_warn r".*" maximize_quantum_search(dqss0)
        @test_nowarn maximize_quantum_search(dqss)
    end
end


text:
"""
    k_i = get_k_i(pb::Problem)

    Compute the degree of each constraint of `pb`.
"""
function get_ki(problem::Problem)
    println("\n=== get_ki(problem::Problem)")
    println("Compute the degree of each constraint of `pb`.")
    return Dict{String, Int}()
end

"""
    check_di_ki!(d_i, k_i)

    Enforce the condition d_i - k_i ≥ 0 on each constraint (if needed).
"""
function check_di_ki!(d_i, k_i)
    println("\n=== check_di_ki!(d_i, k_i)")
    println("Enforce the condition d_i - k_i ≥ 0 on each constraint (if needed).")
end

"""
    SparsityPattern

    Type for storing and working on sparsitty patterns.
"""
type SparsityPattern end

"""
    sparsity_pattern = compute_sparsitypattern(problem, d_i, k_i)

    Compute the sparsity_pattern corresponding to the given partial orders.
"""
function compute_sparsitypattern(problem::Problem, d_i, k_i)
    println("\n=== compute_sparsitypattern(problem::Problem, d_i, k_i)")
    println("Compute the sparsity_pattern corresponding to the given partial orders.")
    sparsity_pattern = SparsityPattern()
    return sparsity_pattern::SparsityPattern
end

"""
    compute_chordalextension!(sparsity_pattern)

    Compute the chordal extension on the provided sparsity_pattern.
"""
function compute_chordalextension!(sparsity_pattern::SparsityPattern)
    println("\n=== compute_chordalextension!(sparsity_pattern::SparsityPattern)")
    println("Compute the chordal extension on the provided sparsity_pattern.")
end

"""
    maxcliques = compute_maxcliques(sparsity_pattern)

    Compute a `Array{Set{Variable}}` describing the maximum cliques on the provided sparsity_pattern.
"""
function compute_maxcliques(sparsity_pattern::SparsityPattern)
    println("\n=== compute_maxcliques(sparsity_pattern::SparsityPattern)")
    println("Compute a `Array{Set{Variable}}` describing the maximum cliques on the provided sparsity_pattern.")
    maxcliques = Array{Set{Variable}}()
    return maxcliques::Array{Set{Variable}}
end

"""
    varsbycstr = compute_varsbycstr(problem)

    Compute a `Dict{String, Set{Variable}}` providing the set of variables involved in each constraint.
"""
function compute_varsbycstr(problem::Problem)
    println("\n=== compute_varsbycstr(problem::Problem)")
    println("Compute a `Dict{String, Set{Variable}}` providing the set of variables involved in each constraint.")
    varsbycstr = Dict{String, Set{Variable}}()
    return varsbycstr::Dict{String, Set{Variable}}
end

"""
    cliquevarsbycstr = compute_varsbycstr(sparsity_pattern, max_cliques, varsbycstr)

    Compute a `Dict{String, Set{Variable}}` providing the set of variables involved in the SDP localizing matrix corresponding to each constraint.
"""
function compute_varsbycstr(sparsity_pattern, max_cliques, varsbycstr)
    println("\n=== compute_varsbycstr(sparsity_pattern, max_cliques, varsbycstr)")
    println("Compute a `Dict{String, Set{Variable}}` providing the set of variables involved in the SDP localizing matrix corresponding to each constraint.")
    cliquevarsbycstr = Dict{String, Set{Variable}}()
    return cliquevarsbycstr::Dict{String, Set{Variable}}
end

"""
    orderbyclique = compute_cliqueorders(sparsity_pattern, d_i, varsbycstr, max_cliques)

    Compute a `Array{Int}` providing the relaxation order corresponding to each clique.
"""
function compute_cliqueorders(sparsity_pattern, d_i, varsbycstr, max_cliques)
    println("\n=== compute_cliqueorders(sparsity_pattern, d_i, varsbycstr, max_cliques)")
    println("Compute a `Array{Int}` providing the relaxation order corresponding to each clique.")
    orderbyclique = Dict{Int, Int}()
    return orderbyclique::Dict{Int, Int}
end

"""
    B_i_dict = compute_Bibycstr(problem, max_cliques, cliquevarsbycstr, orderbyclique)

    Compute a the decomposition of each SDP constraint on the moment basis, yielding several matrices B_i,α,β for each constraint i.
"""
function compute_Bibycstr(problem, max_cliques, cliquevarsbycstr, orderbyclique)
    println("\n=== compute_Bibycstr(problem, max_cliques, cliquevarsbycstr, orderbyclique)")
    println("Compute a the decomposition of each SDP constraint on the moment basis, yielding several matrices B_i,α,β for each constraint i.")
    return
end

"""
    SDP_SOS = build_SDP_SOS(problem, d_i, max_cliques, B_i, cliquevarsbycstr, orderbyclique)

    Build the primal SDP corresponding to the dual SOS hierarchy for the provided problem.
"""
function build_SDP_SOS(problem, d_i, max_cliques, B_i, cliquevarsbycstr, orderbyclique)
    println("\n=== build_SDP_SOS(problem, d_i, max_cliques, B_i, cliquevarsbycstr, orderbyclique)")
    println("Build the primal SDP corresponding to the dual SOS hierarchy for the provided problem.")
    return
end

"""
    m = make_JuMPproblem(SDP_SOS)

    Convert the SDP_SOS problem into a JuMP problem
"""
function make_JuMPproblem(SDP_SOS)
    println("\n=== make_JuMPproblem(SDP_SOS)")
    println("Convert the SDP_SOS problem into a JuMP problem")
    return
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
struct Translation{Dim,T} <: Transform
  vec::Point{Dim,T}
end
Translation(vals::T...) where {T} = Translation(Point{length(vals),T}(vals))
Translation(vec::AbstractVector) where {Dim,T} = Translation(Point{length(vec),eltype(vec)}(vals))
(t::Translation)(p::AbstractVector) = t.vec + p

Base.:∘(t1::Translation, t2::Translation) = Translation(t1.vec .+ t2.vec)

(≈)(x::Translation, y::Translation) = x.vec ≈ y.vec

inv(tr::Translation) = Translation(-tr.vec)

const Translated{O,Dim,T} = Transformed{O,Translation{Dim,T}}
Translated(obj::O, transf::Translation{Dim,T}) where {O,Dim,T} = Translated{O,Dim,T}(obj, transf)