text:
import StackTraces
using FactCheck

@noinline child() = StackTraces.stacktrace()
@noinline parent() = child()
@noinline grandparent() = parent()

@noinline bad_function() = nonexistent_var
@noinline function good_function()
    try
        bad_function()
    catch
        return StackTraces.catch_stacktrace()
    end
end

format_stack = [
    StackTraces.StackFrame(:frame1, "path/file.1", 10, "path/file.inline", 0, false, 0),
    StackTraces.StackFrame(:frame2, "path/file.2", 20, Symbol(""), -1, false, 0)
]

facts() do
    context("basic") do
        stack = grandparent()
        @fact [:child, :parent, :grandparent] --> [f.func for f in stack[1:3]]
        for (line, frame) in zip(5:7, stack[1:3])
            @fact [Symbol(@__FILE__), line] -->
                anyof([frame.file, frame.line], [frame.inlined_file, frame.inlined_line])
        end
        @fact [false, false, false] --> [f.from_c for f in stack[1:3]]
    end

    context("from_c") do
        default, with_c, without_c = StackTraces.stacktrace(), StackTraces.stacktrace(true), StackTraces.stacktrace(false)
        @fact default --> without_c
        @fact length(with_c) --> greater_than(length(without_c))
        @fact filter(frame -> frame.from_c, with_c) --> not(isempty)
        @fact filter(frame -> frame.from_c, without_c) --> isempty
    end

    context("remove_frames!") do
        stack = StackTraces.remove_frames!(grandparent(), :parent)
        @fact stack[1] --> StackTraces.StackFrame(:grandparent, @__FILE__, 7, Symbol(""), -1, false, 0)

        stack = StackTraces.remove_frames!(grandparent(), [:child, :something_nonexistent])
        @fact stack[1:2] --> [
            StackTraces.StackFrame(:parent, @__FILE__, 6, Symbol(""), -1, false, 0),
            StackTraces.StackFrame(:grandparent, @__FILE__, 7, Symbol(""), -1, false, 0)
        ]
    end

    context("try...catch") do
        stack = good_function()
        @fact stack[1:2] --> [
            StackTraces.StackFrame(:bad_function, @__FILE__, 9, Symbol(""), -1, false, 0),
            StackTraces.StackFrame(:good_function, @__FILE__, 12, Symbol(""), -1, false, 0)
        ]
    end

    context("unkown") do
        @fact StackTraces.lookup(C_NULL) --> StackTraces.UNKNOWN
    end

    context("formatting") do
        context("frame") do
            @fact StackTraces.format_stackframe(format_stack[1]) -->
                "[inlined code from file.1:10] frame1 at file.inline:0"
            @fact StackTraces.format_stackframe(format_stack[1]; full_path=true) -->
                "[inlined code from path/file.1:10] frame1 at path/file.inline:0"

            @fact StackTraces.format_stackframe(format_stack[2]) --> "frame2 at file.2:20"
            @fact StackTraces.format_stackframe(format_stack[2]; full_path=true) -->
                "frame2 at path/file.2:20"
        end

        context("stack") do
            @fact StackTraces.format_stacktrace(format_stack, ", ") -->
                "[inlined code from file.1:10] frame1 at file.inline:0, frame2 at file.2:20"
            @fact StackTraces.format_stacktrace(format_stack, ", "; full_path=true) -->
                string(
                    "[inlined code from path/file.1:10] ",
                    "frame1 at path/file.inline:0, frame2 at path/file.2:20"
                )

            @fact StackTraces.format_stacktrace(format_stack, ", ", "Stack: ") -->
                string(
                    "Stack: [inlined code from file.1:10] ",
                    "frame1 at file.inline:0, frame2 at file.2:20"
                )
            @fact StackTraces.format_stacktrace(format_stack, ", ", "Stack: "; full_path=true) -->
                string(
                    "Stack: [inlined code from path/file.1:10] ",
                    "frame1 at path/file.inline:0, frame2 at path/file.2:20"
                )

            @fact StackTraces.format_stacktrace(format_stack, ", ", "{", "}") -->
                string(
                    "{[inlined code from file.1:10] ",
                    "frame1 at file.inline:0, frame2 at file.2:20}"
                )
            @fact StackTraces.format_stacktrace(format_stack, ", ", "{", "}", full_path=true) -->
                string(
                    "{[inlined code from path/file.1:10] ",
                    "frame1 at path/file.inline:0, frame2 at path/file.2:20}"
                )
        end

        context("empty") do
            @fact StackTraces.format_stacktrace(StackTraces.StackTrace(), ", ") --> ""
            @fact StackTraces.format_stacktrace(StackTraces.StackTrace(), ", ", "Stack: ") --> ""
            @fact StackTraces.format_stacktrace(StackTraces.StackTrace(), ", ", "{", "}") --> ""
        end
    end

    context("output") do
        io = IOBuffer()
        StackTraces.show_stacktrace(io, format_stack)
        @fact takebuf_string(io) -->
            """
            StackTrace with 2 StackFrames:
              [inlined code from file.1:10] frame1 at file.inline:0
              frame2 at file.2:20
            """
        StackTraces.show_stacktrace()   # Improves code coverage.
    end
end

FactCheck.exitstatus()


text:
"""
     text3D(str, anchor::Point3D;
        halign=:left,
        valign=:baseline,
        about=Point3D(0., 0., 0.),
        rotation::Rotation=RotXYZ(0, 0, 0))

Draw text at point `pt`, lying in the plane of the x axis.
Angles in `rotation` rotate the text about the `about` point,
defaulting to `Point3D(0, 0, 0)`.

Uses Luxor's `fontface()` and `fontsize()` settings.

Specify rotations using functions from Rotations.jl, such as:

- `RotX(a)`
- `RotZ(a)`
- `RotXZ(a1, a2)`
- `RotXYZ(a1, a2, a3)`
"""
function text3D(str, anchor::Point3D;
        halign=:left,
        valign=:baseline,
        about=Point3D(0., 0., 0.),
        rotation::Rotation=RotXYZ(0, 0, 0))

    textoutlines(str, O, :path,
        halign=halign,
        valign=valign,
        startnewpath=true)

    o = getpathflat()
    newpath() # otherwise the path would be drawn twice
     for e in o
        if e.element_type == 0
            (x, y) = e.points
            newpt = rotateby(Point3D(anchor.x + x, anchor.y -y, anchor.z), about, rotation)
            pin(newpt,  gfunction = (p3, p2) -> move(p2))
        elseif e.element_type == 1
            (x, y) = e.points
            newpt = rotateby(Point3D(anchor.x + x, anchor.y -y, anchor.z), about, rotation)
            pin(newpt,  gfunction = (p3, p2) -> line(p2))
        elseif e.element_type == 3
            closepath()
        else
            error("unknown path element " * repr(e.element_type) * repr(e.points))
        end
    end
    fillpath()
end

ext: jl


text:
using ITensors
using OptimKit
using Random
using Zygote

nsites = 4 # Number of sites
nlayers = 2 # Layers of gates in the ansatz
gradtol = 1e-4 # Tolerance for stopping gradient descent

# The Hamiltonian we are minimizing
function ising_hamiltonian(nsites; h)
  ℋ = OpSum()
  for j in 1:(nsites - 1)
    ℋ += -1, "Z", j, "Z", j + 1
  end
  for j in 1:nsites
    ℋ += h, "X", j
  end
  return ℋ
end

# A layer of the circuit we want to optimize
function layer(nsites, θ⃗)
  RY_layer = [("Ry", (n,), (θ=θ⃗[n],)) for n in 1:nsites]
  CX_layer = [("CX", (n, n + 1)) for n in 1:2:(nsites - 1)]
  return [RY_layer; CX_layer]
end

# The variational circuit we want to optimize
function variational_circuit(nsites, nlayers, θ⃗)
  range = 1:nsites
  circuit = layer(nsites, θ⃗[range])
  for n in 1:(nlayers - 1)
    circuit = [circuit; layer(nsites, θ⃗[range .+ n * nsites])]
  end
  return circuit
end

s = siteinds("Qubit", nsites)

h = 1.3
ℋ = ising_hamiltonian(nsites; h=h)
H = MPO(ℋ, s)
ψ0 = MPS(s, "0")

#
# The loss function, a function of the gate parameters
# and implicitly depending on the Hamiltonian and state:
#
# loss(θ⃗) = ⟨0|U(θ⃗)† H U(θ⃗)|0⟩ = ⟨θ⃗|H|θ⃗⟩
#
function loss(θ⃗)
  nsites = length(ψ0)
  s = siteinds(ψ0)
  𝒰θ⃗ = variational_circuit(nsites, nlayers, θ⃗)
  Uθ⃗ = ops(𝒰θ⃗, s)
  ψθ⃗ = apply(Uθ⃗, ψ0; cutoff=1e-8)
  return inner(ψθ⃗, H, ψθ⃗; cutoff=1e-8)
end

Random.seed!(1234)
θ⃗₀ = 2π * rand(nsites * nlayers)

@show loss(θ⃗₀)

println("\nOptimize circuit with gradient optimization")

loss_∇loss(x) = (loss(x), convert(Vector, loss'(x)))
algorithm = LBFGS(; gradtol=1e-3, verbosity=2)
θ⃗ₒₚₜ, lossₒₚₜ, ∇lossₒₚₜ, numfg, normgradhistory = optimize(loss_∇loss, θ⃗₀, algorithm)

@show loss(θ⃗ₒₚₜ)

println("\nRun DMRG as a comparison")

sweeps = Sweeps(5)
setmaxdim!(sweeps, 10)
e_dmrg, ψ_dmrg = dmrg(H, ψ0, sweeps)

println("\nCompare variational circuit energy to DMRG energy")
@show loss(θ⃗ₒₚₜ), e_dmrg

nothing


text:
    # import Base: show

"""
The structure used for detailed session requests.
"""
mutable struct SessionDetailsRequest
    id::String
    description::String
    initialPoseType::String
    shouldInitialize::Bool
    SessionDetailsRequest(id::String, description::String, initialPoseType::String="Pose2", shouldInitialize::Bool=false) = new(id, description, initialPoseType, shouldInitialize)
end

"""
The structure used for detailed session responses.
"""
mutable struct SessionDetailsResponse
    sessionId::String
    userId::Union{Nothing, String}
    robotId::Union{Nothing, String}
    environmentIds::Union{Nothing, Vector{String}}
    initialPoseType::String
    nodeCount::Int
    solveCount::Int
    lastSolvedTimestamp::Union{Nothing, String}
    isSolverEnabled::Bool # If 1 then the ad-hoc solver will pick up on it, otherwise will ignore this session.
    shouldInitialize::Bool
    id::String
    description::String
    # lastSolvedResult::Union{Nothing, String}
    # solveTimes::Vector{Any}
    createdTimestamp::String
    lastUpdatedTimestamp::String
    links::Dict{String, Any}
end

function show(io::IO, c::SessionDetailsResponse)
    println(io, "GraffSDK Session:")
    println(io, " - ID: $(c.id)")
    println(io, " - Description: $(c.description)")
    println(io, " - User ID: $(c.userId != nothing ? c.userId : "")")
    println(io, " - Robot ID: $(c.robotId != nothing ? c.robotId : "")")
    println(io, " - Environments: $(c.environmentIds != nothing ? c.environmentIds : "")")
    println(io, " - Initial Pose Type: $(c.initialPoseType)")
    println(io, " - Node Count: $(c.nodeCount)")
    println(io, " - Solver Enabled: $(c.isSolverEnabled)")
    println(io, " - Solve Count: $(c.solveCount)")
    println(io, " - Created: $(c.createdTimestamp)")
    println(io, " - Last Solved: $(c.lastSolvedTimestamp != nothing ? c.lastSolvedTimestamp : "")")
    # println(io, " - Last Solved Result: $(c.lastSolvedResult != nothing ? c.lastSolvedResult : "")")
end

"""
A list of session response summaries.
"""
mutable struct SessionsResponse
    sessions::Vector{SessionDetailsResponse}
    links::Dict{String, String}
end

function show(io::IO, c::SessionsResponse)
    println(io, "GraffSDK Sessions (count = $(length(c.sessions))):")
    for s in c.sessions
        print(io, " - $s")
    end
end

"""
The structure used to briefly describe a node in a response.
"""
mutable struct NodeResponse
    id::Int
    label::String
    mapEst::Union{Nothing, Vector{Float64}}
    createdTimestamp::String
    lastUpdatedTimestamp::String
    links::Dict{String, Any}
end

function show(io::IO, n::NodeResponse)
    println(io, "GraffSDK Node - ID: $(n.id), label: $(n.label), mapEst: $(n.mapEst != nothing ? n.mapEst : "nothing")")
end

"""
The structure describing a complete node in a response.
"""
mutable struct NodeDetailsResponse
    dataEntries::Union{Nothing, Vector{BigDataEntryResponse}}
    type::String
    packed::Dict{String, Any}
    smalldata::String
    ready::Int
    backendset::Int
    id::Int
    label::String
    mapEst::Union{Nothing, Vector{Float64}}
    createdTimestamp::String
    lastUpdatedTimestamp::String
    links::Union{Nothing, Dict{String, Any}}
    tags::Vector{String}
end

function show(io::IO, n::NodeDetailsResponse)
    println(io, "GraffSDK Node:")
    println(io, " - ID: $(n.id)")
    println(io, " - Label: $(n.label)")
    println(io, " - Graph Labels: $(n.tags)")
end

"""
The structure describing a high-level add-odometry request.
"""
mutable struct AddOdometryRequest
    timestamp::String
    deltaMeasurement::Vector{Float64}
    pOdo::Array{Float64, 2}
    N::Int64
    AddOdometryRequest(deltaMeasurement::Vector{Float64}, pOdo::Array{Float64, 2}, N::Int64) = new(string(now(UTC)), deltaMeasurement, pOdo, N)
    AddOdometryRequest(deltaMeasurement::Vector{Float64}, pOdo::Array{Float64, 2}) = new(string(now(UTC)), deltaMeasurement, pOdo, 100)
end

"""
Result of an AddOdometryRequest - provides the variable and factor details for the created node.
"""
struct AddOdometryResponse
    variable::NodeResponse
    factor::NodeResponse
end

"""
The parameters structure for CreateVariable request.
"""
mutable struct VariableRequest
    label::String
    variableType::String
    N::Int64
    labels::Vector{String}
    VariableRequest(label::String, variableType::String, N::Int64, labels::Vector{String}) = new(label, variableType, N, labels)
    VariableRequest(label::String, variableType::String, labels::Vector{String}) = new(label, variableType, 100, labels)
    VariableRequest(label::String, variableType::String) = new(label, variableType, 100, String[])
end

"""
Result of a CreateVariableRequest.
"""
struct VariableResponse
end

"""
Parameters for a general distribution request - the distribution type and the accompanying parameters.
"""
mutable struct DistributionRequest
    distType::String
    params::Vector{Float64}
end

"""
A 2D bearing+range request body.
"""
mutable struct BearingRangeRequest
    pose2Id::String
    point2Id::String
    bearing::DistributionRequest
    range::DistributionRequest
end

"""
Parameter for a CreateFactor request - the factor type and packed factor details.
"""
mutable struct FactorBody
    factorType::String
    packedFactorType::String
    encoding::String
    body::String
end

"""
The body of a CreateFactor request - the variables to be linked, the body of the factor, and whether it should be autoinitialized and is ready for solving.
"""
mutable struct FactorRequest
    variables::Vector{String}
    body::FactorBody
    autoinit::Bool
    ready::Bool
    FactorRequest(variables::Vector{String}, factorType::String, packedFactor; autoinit::Bool = false, ready::Bool = false ) = begin
        # try
            #TODO: Simplify this
            factBody = FactorBody(factorType, string(typeof(packedFactor)), "application/json", JSON.json(packedFactor))
            return FactorRequest(variables, factBody; autoinit=autoinit, ready=ready)
        # catch ex
        #     error("Unable to serialize the factor body - $ex")
        # end
    end
    FactorRequest(variables::Vector{String}, body::FactorBody; autoinit::Bool = false, ready::Bool = false ) = new(variables, body, autoinit, ready)
end


text:
abstract type AbstractRegression end

struct LocalLinear <: AbstractRegression

    k::Int64
    ivar::Vector{Int64}
    ivar_neighboor::Vector{Int64}
    res::Array{Float64,2}
    beta::Array{Float64,2}
    Xm::Array{Float64,2}
    Xr::Array{Float64,2}
    Cxx::Array{Float64,2}
    Cxx2::Array{Float64,2}
    Cxy::Array{Float64,2}
    pred::Array{Float64,2}
    X0r::Array{Float64,2}

    function LocalLinear(k, ivar, ivar_neighboor)

        #X = zeros(Float64,(nvn, k))
        #Y = zeros(Float64,(nv, k))
        #w = zeros(Float64, k)

        nv = length(ivar)
        nvn = length(ivar_neighboor)

        res = zeros(Float64, (nv, k))
        beta = zeros(Float64, (nv, nvn + 1))
        Xm = zeros(Float64, (nvn, 1))
        Xr = ones(Float64, (nvn + 1, k))
        Cxx = zeros(Float64, (nvn + 1, nvn + 1))
        Cxx2 = zeros(Float64, (nvn + 1, nvn + 1))
        Cxy = zeros(Float64, (nv, nvn + 1))
        pred = zeros(Float64, (nv, k))
        X0r = ones(Float64, (nvn + 1, 1))

        new(k, ivar, ivar_neighboor, res, beta, Xm, Xr, Cxx, Cxx2, Cxy, pred, X0r)

    end

end

function compute(ll::LocalLinear, x, xf_tmp, xf_mean, ip, X, Y, w)

    ivar = ll.ivar
    ivar_neighboor = ll.ivar_neighboor

    ll.Xm .= sum(X .* w', dims = 2)
    ll.Xr[2:end, :] .= X .- ll.Xm
    mul!(ll.Cxx, (ll.Xr .* w'), ll.Xr')
    mul!(ll.Cxx2, (ll.Xr .* w' .^ 2), ll.Xr')
    mul!(ll.Cxy, (Y .* w'), ll.Xr')

    ll.Cxx .= pinv(ll.Cxx, rtol = 0.01)
    ll.Cxx2 .= ll.Cxx2 * ll.Cxx
    # regression on principal components
    mul!(ll.beta, ll.Cxy, ll.Cxx)
    ll.X0r[2:end, :] .= x[ivar_neighboor, ip] .- ll.Xm
    # weighted mean
    xf_mean[ivar, ip] = ll.beta * ll.X0r
    mul!(ll.pred, ll.beta, ll.Xr)
    Y .-= ll.pred
    xf_tmp[ivar, :] .= xf_mean[ivar, ip] .+ Y
    # weigthed covariance
    cov_xf = (Y * (w .* Y')) ./ (1 .- tr(ll.Cxx2))
    cov_xf .= Symmetric(cov_xf .* (1 .+ tr(ll.Cxx2 * ll.X0r * ll.X0r' * ll.Cxx)))

    return cov_xf

end


text:
function solve_bvp(problem::SpectralBVPProblem{TF, TBCL, TBCR}) where {TF, TBCL<:PeriodicBC, TBCR<:PeriodicBC}
    rhshat = fft(problem.rhs);
    uhat = similar(rhshat);
    @. uhat[2:end] = rhshat[2:end]/problem.λ[2:end];
    u = ifft(uhat);
    return u
end

function solve_bvp(problem::SpectralBVPProblem{TF, TBCL, TBCR}) where {TF, TBCL<:DirichletBC, TBCR<:DirichletBC}
    rhshat = fft(problem.rhs);
    uhat = similar(rhshat);
    @. uhat[2:end] = rhshat[2:end]/problem.λ[2:end];
    u_ = ifft(uhat)
    u = u_[problem.N÷2 + 2:end];
    return u
end

function solve_bvp(problem::SpectralBVPProblem{TF, TBCL, TBCR}) where {TF, TBCL<:NeumannBC, TBCR<:NeumannBC}
    rhshat = fft(problem.rhs);
    uhat = similar(rhshat);
    @. uhat[2:end] = rhshat[2:end]/problem.λ[2:end];
    u_ = ifft(uhat);
    u = [u_[problem.N÷2 + 1:end]; u_[1]];
    return u
end


text:
# FileDataset

function FileDataset(dir, pattern = "*")
    return rglob(pattern, string(dir))
end

pathparent(p::String) = splitdir(p)[1]
pathname(p::String) = splitdir(p)[2]

# File utilities

"""
    rglob(filepattern, dir = pwd(), depth = 4)

Recursive glob up to 6 layers deep.
"""
function rglob(filepattern = "*", dir = pwd(), depth = 4)
    patterns = [
        "$filepattern",
        "*/$filepattern",
        "*/*/$filepattern",
        "*/*/*/$filepattern",
        "*/*/*/*/$filepattern",
        "*/*/*/*/*/$filepattern",
    ]
    return vcat([glob(pattern, dir) for pattern in patterns[1:depth]]...)
end

"""
    loadfile(file)

Load a file from disk into the appropriate format.
"""
function loadfile(file::String)
    if isimagefile(file)
        # faster image loading
        return FileIO.load(file, view = true)
    elseif endswith(file, ".csv")
        return DataFrame(CSV.File(file))
    else
        return FileIO.load(file)
    end
end

loadfile(file::AbstractPath) = loadfile(string(file))

#TableDataset

struct TableDataset{T}
    table::T #Should implement Tables.jl interface
    TableDataset{T}(table::T) where {T} =
        Tables.istable(table) ? new{T}(table) :
        error("Object doesn't implement Tables.jl interface")
end

TableDataset(table::T) where {T} = TableDataset{T}(table)
TableDataset(path::AbstractPath) = TableDataset(DataFrame(CSV.File(path)))

function LearnBase.getobs(dataset::FastAI.Datasets.TableDataset, idx)
    if Tables.rowaccess(dataset.table)
        row, _ = Iterators.peel(Iterators.drop(Tables.rows(dataset.table), idx - 1))
        return row
    elseif Tables.columnaccess(dataset.table)
        colnames = Tables.columnnames(dataset.table)
        rowvals = [Tables.getcolumn(dataset.table, i)[idx] for i = 1:length(colnames)]
        return (; zip(colnames, rowvals)...)
    else
        error(
            "The Tables.jl implementation used should have either rowaccess or columnaccess.",
        )
    end
end

function LearnBase.nobs(dataset::TableDataset)
    if Tables.columnaccess(dataset.table)
        return length(Tables.getcolumn(dataset.table, 1))
    elseif Tables.rowaccess(dataset.table)
        return length(Tables.rows(dataset.table)) # length might not be defined, but has to be for this to work.
    else
        error(
            "The Tables.jl implementation used should have either rowaccess or columnaccess.",
        )
    end
end

LearnBase.getobs(dataset::TableDataset{<:DataFrame}, idx) = dataset.table[idx, :]
LearnBase.nobs(dataset::TableDataset{<:DataFrame}) = nrow(dataset.table)

LearnBase.getobs(dataset::TableDataset{<:CSV.File}, idx) = dataset.table[idx]
LearnBase.nobs(dataset::TableDataset{<:CSV.File}) = length(dataset.table)

# ## Tests

@testset "TableDataset" begin
    @testset "TableDataset from rowaccess table" begin
        Tables.columnaccess(::Type{<:Tables.MatrixTable}) = false
        Tables.rowaccess(::Type{<:Tables.MatrixTable}) = true

        testtable = Tables.table([1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"])
        td = TableDataset(testtable)

        @test all(getobs(td, 1) .== [1, 4.0, "7"])
        @test nobs(td) == 3
    end

    @testset "TableDataset from columnaccess table" begin
        Tables.columnaccess(::Type{<:Tables.MatrixTable}) = true
        Tables.rowaccess(::Type{<:Tables.MatrixTable}) = false

        testtable = Tables.table([1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"])
        td = TableDataset(testtable)

        @test [data for data in getobs(td, 2)] == [2, 5.0, "8"]
        @test nobs(td) == 3

        @test getobs(td, 1) isa NamedTuple
    end

    @testset "TableDataset from DataFrames" begin
        testtable = DataFrame(
            col1 = [1, 2, 3, 4, 5],
            col2 = ["a", "b", "c", "d", "e"],
            col3 = [10, 20, 30, 40, 50],
            col4 = ["A", "B", "C", "D", "E"],
            col5 = [100.0, 200.0, 300.0, 400.0, 500.0],
            split = ["train", "train", "train", "valid", "valid"],
        )
        td = TableDataset(testtable)
        @test td isa TableDataset{<:DataFrame}

        @test [data for data in getobs(td, 1)] == [1, "a", 10, "A", 100.0, "train"]
        @test nobs(td) == 5
    end

    @testset "TableDataset from CSV" begin
        open("test.csv", "w") do io
            write(io, "col1,col2,col3,col4,col5, split\n1,a,10,A,100.,train")
        end
        testtable = CSV.File("test.csv")
        td = TableDataset(testtable)
        @test td isa TableDataset{<:CSV.File}
        @test [data for data in getobs(td, 1)] == [1, "a", 10, "A", 100.0, "train"]
        @test nobs(td) == 1
        rm("test.csv")
    end
end


text:
#
#  Use git:                    git status ....;   git add <filenames>;   git commit -m "..";   git push;   git rm  <filenames>
#  Use Jupyter notebooks:      using IJulia;   notebook()
#  Activation:                 ];   pkg> up;   pkg> activate
#  Working with JAC:           using Revise;   using JAC;   include("../src/jac.jl");   pkg> test
#
#  Copy to desktop             scp -r JAC.jl/ fritzsch@10.140.119.236:~/fri/.
"""
`module JAC`
    ... Jena Atomic Calculator (JAC) provides tools for performing atomic (structure) calculations at various degrees of complexity
        and sophistication. It has been designed to not only calculate atomic level structures and properties [such as g-factors or
        hyperfine and isotope-shift parameters] but also transition amplitudes between bound-state levels [for the anapole moment, dipole
        operator, electron electric-dipole moment, parity non-conservation, etc.] and, in particular, (atomic) transition probabilities,
        Auger rates, photoionization cross sections, radiative and dielectronic recombination rates as well as cross sections for many
        other (elementary) processes. JAC also facilitates interactive computations, the simulation of atomic cascades, the time-evolution
        of statistical tensors, a few semi-empirical estimates of atomic properties as well as the simplification of symbolic expressions
        from Racah's algebra. -- In addition, the JAC module supports the display of level energies, electron and photon spectra,
        radial orbitals and other atomic data.

**`Perform (atomic) computations of different complexity:`**
    JAC will eventually support **ten kinds** of computations which can be summarized as follows:

   + Atomic computations, based on explicitly specified electron configurations.
   + Restricted active-space computations (RAS).
   + Interactive computations.
   + Atomic cascade computations (partly implemented).
   + Atomic representations (Green and close-coupling functions, complex rotation; partly implemented).
   + Atomic responses (partly implemented).
   + Atomic descriptors for machine learning algorithms (not yet implemented).
   + Time-evolution of statistical tensors in (intense) light pusles (not yet implemented).
   + Semi-empirical estimates of cross sections, etc. (partly implemented).
   + Symbolic evaluation of expressions from Racah's algebra, etc.

**`Further details and information`**

        + Kinds of atomic implementation                                       [cf. ? Details.kindsOfComputation]
        + Atomic amplitudes (partly) implemented in JAC                        [cf. ? Details.amplitudes]
        + Atomic level properties (partly) implemented in JAC                  [cf. ? Details.properties]
        + Atomic processes (partly) implemented in JAC                         [cf. ? Details.processes]
        + Interactive use of JAC procedures                                    [cf. ? Details.interactive]
        + Design principles and limitations of the JAC program                 [cf. ? Details.design]
        + Data types, structs and name conventions of the JAC module           [cf. ? Details.datatypes]
        + Atomic cascade computations and approximations                       [cf. ? Details.decayCascades]
        + Use of (em) light pulses in the time evolution of statist. tensors   [cf. ? Details.pulses]
        + Why Julia ?                                                          [cf. ? Details.whyJulia]

"""
module JAC

using  Dates, Printf,  LinearAlgebra, IJulia, SpecialFunctions, FortranFiles, QuadGK, GSL, JLD2, SymEngine,
       HypergeometricFunctions  ## , Interact, GaussQuadrature

export AbstractConfigurationRestriction, AbstractEeInteraction, AbstractCImethod, AbstractPotential, AbstractQedModel, AbstractStartOrbitals,
       AbstractProcessSettings, AbstractPropertySettings,
       add, analyze, AlphaX, AlphaVariation, AnapoleMoment,
       AngularJ64, AngularM64, AngularJ, AngularMomentum,
       AsfSettings, Atomic, AtomicState, AtomicStructure, Auger, AugerInPlasma, AutoIonization,
       Basics, Basis, BreitInteraction, Bsplines,
       CartesianVector, CiSettings, CiExpansion, CloseCoupling, compute, convertUnits, Cascade, Compton, Configuration, ConfigurationR,
       Continuum, CsfR, Coulex, CoulombExcitation, Coulion, CoulombIonization, CoulombBreit, CoulombInteraction, ClebschGordan,
       CorePolarization,
       diagonalize, Defaults, DecayYield, Details, Dielectronic, Dierec, Djpq, DoubleAutoIonization, DoubleAuger,
       Eimex, ElectronCapture, ElecCapture, estimate, ElectricDipoleMoment, Einstein, EinsteinX, EmMultipole, evaluate, ExpStokes,
       E1, M1, E2, M2, E3, M3, E4, M4,
       FormFactor, FormF, FullCIeigen,
       generate, GreenSettings, GreenChannel, GreenExpansion, getDefaults, Green, Gui,
       Hfs, HighHarmonic, HFS, HydrogenicIon,
       interpolate, integrate, Integral, ImpactExcAuto, ImpactExcitation, ImpactExcitationAutoion, ImpactIonization,
       InternalConv, InternalConversion, Isotope, IsotopeShift,
       Kronecker,
       LandeZeeman, LandeJ,  LandeF, Level, LevelSymmetry, LSjj, LSjjSettings, LevelSelection, LineSelection,
       ManyElectron, MeanFieldSettings, MeanFieldBasis, minus, Model, modify, MultiPhotonDE, MultiPhotonDeExcitation, MultiPhotonDoubleIon,
       MultiPI, MultiPDI, MultiPhotonIonization, MultipoleMoment, MultipolePolarizibility, Multiplet,
       NoAmplitude, NoProcess, Nuclear, NoneQed, NoProperty,
       OneElectronSettings, OneElectronSpectrum, Orbital,
       PathwaySelection, perform, provide, PairA1P, PairAnnihilation1Photon, PairAnnihilation2Photon,
       PairProduction, ParityNonConservation, PeriodicTable, Parity, plus,
       Photo, PhotoExc, PhotoExcAuto, PhotoExcFluor, PhotoEmission, PhotoExcitation, PhotoExcitationAutoion, PhotoExcitationFluores,
       PhotoIonAuto, PhotoIonFluor, PhotoDouble, PhotoIonization, PhotoDoubleIonization, PhotoIonizationFluores, PhotoIonizationAutoion,
       PhotoRecombination, PlasmaShift, Plasma, Polarity, Pulse,
       QedPetersburg, QedSydney,
       RacahAlgebra, RacahExpression, Radial, RadialIntegrals, Radiative, RadiativeAuger, RAuger, RasSettings, RasStep,
       RasExpansion, RayleighCompton, recast, Rec, REDA, READI, Representation, ReducedDensityMatrix, RadiativeOpacity,
       RestrictNoElectronsTo, RestrictParity, RestrictToShellDoubles, RequestMinimumOccupation, RequestMaximumOccupation,
       SchiffMoment, Semiempirical, setDefaults, Shell, SolidAngle, Spectroscopy, SpinAngular, SphericalTensor,
       StartFromHydrogenic, StartFromPrevious, StrongField, Subshell,
       tabulate, TestFrames, Triangle, tools,
       UseCoulomb, UseBabushkin, UseGauge,
       WeightedCartesian, W3j, W6j, W9j,
       Yields, Ylm,
       Zeeman

# Basic data and data structures
include("module-Basics.jl");            using ..Basics
include("module-Radial.jl");            using ..Radial
include("module-Math.jl");              using ..Math
include("module-Defaults.jl");          using ..Defaults
include("module-ManyElectron.jl");      using ..ManyElectron
include("module-Nuclear.jl");           using ..Nuclear
include("module-BiOrthogonal.jl");      using ..BiOrthogonal

# Specialized functions/methods to manipulate these data
include("module-AngularMomentum.jl")
## include("module-AngularCoefficients-Ratip2013.jl")  ## keep for internal test purposes only
include("module-SpinAngular.jl");       using ..SpinAngular
include("module-Bsplines.jl")
include("module-Continuum.jl")
include("module-Details.jl")
include("module-RadialIntegrals.jl")
include("module-HydrogenicIon.jl")
include("module-InteractionStrength.jl")
include("module-InteractionStrengthQED.jl")
include("module-PeriodicTable.jl")
include("module-TableStrings.jl")
include("module-Tools.jl")
include("module-AtomicState.jl");       using ..AtomicState
include("module-LSjj.jl")

# Functions/methods for atomic amplitudes
include("module-MultipoleMoment.jl")
include("module-ParityNonConservation.jl")

include("module-PhotoEmission.jl")

# Functions/methods for atomic properties
include("module-Einstein.jl")
include("module-Hfs.jl")
include("module-IsotopeShift.jl")
include("module-LandeZeeman.jl")
include("module-FormFactor.jl")
include("module-ReducedDensityMatrix.jl")
include("module-PlasmaShift.jl")
include("module-AlphaVariation.jl")
include("module-DecayYield.jl")
include("module-RadiativeOpacity.jl")
##x include("module-CloseCoupling.jl")
include("module-MultipolePolarizibility.jl")

# Functions/methods for atomic processes
include("module-PhotoExcitation.jl")
include("module-PhotoIonization.jl")
include("module-PhotoRecombination.jl")
include("module-AutoIonization.jl")
include("module-ElectronCapture.jl")
include("module-Dielectronic.jl")
include("module-PhotoExcitationFluores.jl")
include("module-PhotoExcitationAutoion.jl")
include("module-RayleighCompton.jl")
include("module-MultiPhotonDeExcitation.jl")
include("module-DoubleAutoIonization.jl")
include("module-PhotoDoubleIonization.jl")
include("module-CoulombExcitation.jl")
include("module-PhotoIonizationFluores.jl")
include("module-PhotoIonizationAutoion.jl")
include("module-CoulombIonization.jl")
include("module-ImpactExcitation.jl")
include("module-ImpactExcitationAutoion.jl")
include("module-RadiativeAuger.jl")
include("module-MultiPhotonIonization.jl")
include("module-MultiPhotonDoubleIon.jl")
include("module-InternalConversion.jl")
#= Further processes, not yet included into the code
include("module-REDA.jl")
include("module-READI.jl")
include("module-PairProduction.jl")
include("module-PairAnnihilation1Photon.jl")
include("module-PairAnnihilation2Photon.jl")  =#

# Functions/methods for atomic responses and time evolutions
include("module-Pulse.jl")
# include("module-Statistical.jl")

# Functions/methods for the computation of atomic responses
## include("module-HighHarmonic.jl")
include("module-StrongField.jl")

# Functions/methods for semi-empirical estimations
# include("module-ImpactIonization.jl")
include("module-Semiempirical.jl")

# Functions/methods for atomic and cascade computations
include("module-Atomic.jl");            using ..Atomic
include("module-Cascade.jl");           using ..Cascade

# Functions/methods for symbolic computations
include("module-RacahAlgebra.jl");      using ..RacahAlgebra
include("module-SphericalTensor.jl");   using ..SphericalTensor

# Basic functions/methods to manipulate these data
include("module-BasicsAZ.jl")
include("module-ManyElectronAZ.jl")

# Specialized macros
##x include("macro-racahsum.jl")

# All test functions/methods stay with the JAC root module
include("module-TestFrames.jl");        using ..TestFrames

function __init__()
    # The following variables need to be initialized at runtime to enable precompilation
    global JAC_SUMMARY_IOSTREAM = stdout
    global JAC_TEST_IOSTREAM    = stdout
end

println("\nWelcome to JAC:  A community approach to the computation of atomic structures, cascades and time evolutions " *
        "[(C) Copyright by Stephan Fritzsche, Jena (2018-2021)].")

end


text:
using StructuresKit

#362S162-54, Mcrd, no springs, no moment gradient, Lm>>Lcrd

CorZ=0
t=0.0566

#out to out dimensions
ho=3.625
b=1.625
d=0.5
θc=90.0

E=29500.0
μ=0.30
G=E/(2*(1+μ))
f1=1
f2=-1
M1=1.0
M2=1.0
Lm=100.0
kϕ=0.0

Sf=0.9085/(3.625/2)

CurvatureSign = -1.0

Mcrd = AISIS10016.app23331(CorZ, t, ho, b, d, θc, E, μ, G, f1, f2, M1, M2, CurvatureSign, Lm, kϕ, Sf)

#compare to CFSEI TechNote G100-08 (or G100-09?), Table 3
#https://cfsei.memberclicks.net/technical-notes
@test Mcrd ≈ 89.2*Sf atol=0.01*Mcrd


text:
"""
Automated generator for constraint append function
"""

function genetate_equalityconstraints_function(nh=100)
	open("append_equalityconstraints.jl", "w+") do io
		write(io, "\"\"\"\n")
		write(io, "    append_hs!(model::JuMP.Model, x, nh::Int)\n\n")
	    write(io, "Append equalty constraints (auto-generated)\n\n")
	    write(io, "# Arguments\n")
	    write(io, "    `model::JuMP.Model`: optimization model\n")
	    write(io, "    `x::JuMP.variables`: optimization model variables\n")
	    write(io, "    `nh::Int`: number of equality constraints\n")
	    write(io, "\"\"\"\n")
	    write(io, "function append_hs!(model::JuMP.Model, x, nh::Int)\n")
	    write(io, "    for ih = 1:nh \n")
		write(io, "        if ih == 1 \n")
		write(io, "            @NLconstraint(model, h1(x...) == 0)\n")
	    for j = 2:nh
		    write(io, "        elseif ih == $j \n")
		    write(io, "            @NLconstraint(model, h$j(x...) == 0)\n")
	    end
	     write(io, "        end\n")
	    write(io, "    end\n")
	    write(io, "end\n")
	end
end

function genetate_inequalityconstraints_function(ng=100)
	open("append_inequalityconstraints.jl", "w+") do io
		write(io, "\"\"\"\n")
		write(io, "    append_gs!(model::JuMP.Model, x, ng::Int)\n\n")
	    write(io, "Append inequality constraints (auto-generated)\n\n")
	    	    write(io, "# Arguments\n")
	    write(io, "    `model::JuMP.Model`: optimization model\n")
	    write(io, "    `x::JuMP.variables`: optimization model variables\n")
	    write(io, "    `ng::Int`: number of inequality constraints\n")
	    write(io, "\"\"\"\n")
	    write(io, "function append_gs!(model::JuMP.Model, x, ng::Int)\n")
	    write(io, "    for ig = 1:ng \n")
		write(io, "        if ig == 1 \n")
		write(io, "            @NLconstraint(model, -Inf <= g1(x...) <= 0)\n")
	    for j = 2:ng
		    write(io, "        elseif ig == $j \n")
		    write(io, "            @NLconstraint(model, -Inf <= g$j(x...) <= 0)\n")
	    end
	     write(io, "        end\n")
	    write(io, "    end\n")
	    write(io, "end\n")
	end
end

genetate_equalityconstraints_function()
genetate_inequalityconstraints_function()

max_stars_repo_name: Yuricst/NLPSaUT

max_issues_repo_name: Yuricst/NLPSaUT

max_forks_repo_name: Yuricst/NLPSaUT


text:
#=
    There are six general cases for processing a/b (when b excludes 0)
    and six special cases (where the denominator includes 0 as a bound)
    and two exceptional cases (where the denominator is 0 or straddles 0)
                  b <= 0       b >= 0
       a <=  0    LteLte       LteGte
       a <0< a    ZerLte       ZerGte
       a >=  0    GteLte       GteGte
=#

function (/){G1<:Grasp,G2<:Grasp,R<:Real}(a::Rvl{G1,R}, b::Rvl{G2,R})
    z = zero(R)
    if (b.lo == z) | (b.hi == z) | ((b.lo < z) & (b.hi > z))
      divisorContainsZero(a,b)
   elseif a.hi <= z
       if     b.hi <= z
           divLteLte(a,b)
       else # b.lo.fp >= z
           divLteGte(a,b)
       end
   elseif a.lo >= z
       if     b.hi <= z
           divGteLte(a,b)
       else # b.lo >= 0
           divGteGte(a,b)
       end
   else   # a straddles 0
       if     b.hi <= z
           divZerLte(a,b)
       else # b.lo >= 0
           divZerGte(a,b)
       end
   end
end

for (fn,loa,lob,hia,hib) in [ (:divLteLte, :(a.hi), :(b.lo), :(a.lo), :(b.hi)),
                              (:divLteGte, :(a.lo), :(b.lo), :(a.hi), :(b.hi)),
                              (:divGteLte, :(a.hi), :(b.hi), :(a.lo), :(b.lo)),
                              (:divGteGte, :(a.lo), :(b.hi), :(a.hi), :(b.lo)),
                              (:divZerLte, :(a.hi), :(b.hi), :(a.lo), :(b.hi)),
                              (:divZerGte, :(a.lo), :(b.lo), :(a.hi), :(b.lo)),
                            ]
  @eval begin
    function ($fn){G1<:Grasp,G2<:Grasp,R<:Real}(a::Rvl{G1,R}, b::Rvl{G2,R})
        aLoIsOpen, aHiIsOpen = boundaries(G1)
        bLoIsOpen, bHiIsOpen = boundaries(G2)
        abGrasp = boundaries( (aLoIsOpen|bLoIsOpen), (aHiIsOpen|bHiIsOpen) )

        lo = (/)(($loa), ($lob), RoundDown)
        hi = (/)(($hia), ($hib), RoundUp)

        Rvl{abGrasp,R}(lo, hi)
    end
  end
end

(/){G<:Grasp,R<:Real}(a::Rvl{G,R}, b::R) = (/)(a, Rvl{G,R}(b))
(/){G<:Grasp,R<:Real}(a::R, b::Rvl{G,R}) = (/)(Rvl{G,R}(a), b)

(/){G<:Grasp,I<:WorkInt}(a::Rvl{G,Float64}, b::I) = (/)(a, Rvl{G,Float64}(convert(Float64,b)))
(/){G<:Grasp,I<:WorkInt}(a::I, b::Rvl{G,Float64}) = (/)(Rvl{G,Float64}(convert(Float64,a)), b)
(/){G<:Grasp,I<:WorkInt}(a::Rvl{G,Float32}, b::I) = (/)(a, Rvl{G,Float32}(convert(Float32,b)))
(/){G<:Grasp,I<:WorkInt}(a::I, b::Rvl{G,Float32}) = (/)(Rvl{G,Float32}(convert(Float32,a)), b)

function divisorContainsZero{G1<:Grasp,G2<:Grasp,R<:Real}(a::Rvl{G1,R}, b::Rvl{G2,R})
    aLoIsOpen, aHiIsOpen = boundaries(G1)
    bLoIsOpen, bHiIsOpen = boundaries(G2)
    abGrasp = boundaries( (aLoIsOpen|bLoIsOpen), (aHiIsOpen|bHiIsOpen) )

    if (b.lo < zero(R)) & (b.hi > zero(R))
       return divisorStraddlesZero(a,b)
    end
    if b.lo == zero(R)
       lo = one(R) / b.hi
       hi = (R)(Inf)
    else b.hi == zero(R)
       lo = (R)(-Inf)
       hi = one(R) / b.lo
    end
    Rvl{abGrasp,R}(lo, hi)
end

function divisorStraddlesZero{G1<:Grasp,G2<:Grasp,R<:Real}(a::Rvl{G1,R}, b::Rvl{G2,R})
    ErrorException("Divisor straddles zero: $(a) / $(b)")
end


text:
include("lib_for_practic_10.jl")
#r = Robot()
#r = BorderRobot{CoordRobot}(r)
#show!(r.robot)

function movements_mark!(r,side::HorizonSide)
    while (isborder(r,side) ==  false) && (ismarker(r) == false)
        putmarker!(r)
        move!(r,side)
    end
end

function kvadro!(r)

    count = movements!(r,Ost)

    for side in [Sud, West, Nord, Ost, Sud]
        movements_mark!(r.robot,side)
    end

    movements!(r, West, count)

end

kvadro!(r)
show(r.robot)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
# Autogenerated wrapper script for mlpack_jll for x86_64-apple-darwin
export libmlpack, libmlpack_julia_adaboost, libmlpack_julia_approx_kfn, libmlpack_julia_bayesian_linear_regression, libmlpack_julia_cf, libmlpack_julia_dbscan, libmlpack_julia_decision_stump, libmlpack_julia_decision_tree, libmlpack_julia_det, libmlpack_julia_emst, libmlpack_julia_fastmks, libmlpack_julia_gmm_generate, libmlpack_julia_gmm_probability, libmlpack_julia_gmm_train, libmlpack_julia_hmm_generate, libmlpack_julia_hmm_loglik, libmlpack_julia_hmm_train, libmlpack_julia_hmm_viterbi, libmlpack_julia_hoeffding_tree, libmlpack_julia_image_converter, libmlpack_julia_kde, libmlpack_julia_kernel_pca, libmlpack_julia_kfn, libmlpack_julia_kmeans, libmlpack_julia_knn, libmlpack_julia_krann, libmlpack_julia_lars, libmlpack_julia_linear_regression, libmlpack_julia_linear_svm, libmlpack_julia_lmnn, libmlpack_julia_local_coordinate_coding, libmlpack_julia_logistic_regression, libmlpack_julia_lsh, libmlpack_julia_mean_shift, libmlpack_julia_nbc, libmlpack_julia_nca, libmlpack_julia_nmf, libmlpack_julia_pca, libmlpack_julia_perceptron, libmlpack_julia_preprocess_binarize, libmlpack_julia_preprocess_describe, libmlpack_julia_preprocess_one_hot_encoding, libmlpack_julia_preprocess_scale, libmlpack_julia_preprocess_split, libmlpack_julia_radical, libmlpack_julia_random_forest, libmlpack_julia_softmax_regression, libmlpack_julia_sparse_coding, libmlpack_julia_util

using boost_jll
using armadillo_jll
using OpenBLAS_jll
JLLWrappers.@generate_wrapper_header("mlpack")
JLLWrappers.@declare_library_product(libmlpack, "@rpath/libmlpack.3.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_adaboost, "@rpath/libmlpack_julia_adaboost.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_approx_kfn, "@rpath/libmlpack_julia_approx_kfn.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_bayesian_linear_regression, "@rpath/libmlpack_julia_bayesian_linear_regression.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_cf, "@rpath/libmlpack_julia_cf.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_dbscan, "@rpath/libmlpack_julia_dbscan.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_decision_stump, "@rpath/libmlpack_julia_decision_stump.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_decision_tree, "@rpath/libmlpack_julia_decision_tree.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_det, "@rpath/libmlpack_julia_det.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_emst, "@rpath/libmlpack_julia_emst.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_fastmks, "@rpath/libmlpack_julia_fastmks.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_generate, "@rpath/libmlpack_julia_gmm_generate.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_probability, "@rpath/libmlpack_julia_gmm_probability.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_train, "@rpath/libmlpack_julia_gmm_train.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_generate, "@rpath/libmlpack_julia_hmm_generate.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_loglik, "@rpath/libmlpack_julia_hmm_loglik.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_train, "@rpath/libmlpack_julia_hmm_train.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_viterbi, "@rpath/libmlpack_julia_hmm_viterbi.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_hoeffding_tree, "@rpath/libmlpack_julia_hoeffding_tree.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_image_converter, "@rpath/libmlpack_julia_image_converter.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_kde, "@rpath/libmlpack_julia_kde.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_kernel_pca, "@rpath/libmlpack_julia_kernel_pca.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_kfn, "@rpath/libmlpack_julia_kfn.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_kmeans, "@rpath/libmlpack_julia_kmeans.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_knn, "@rpath/libmlpack_julia_knn.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_krann, "@rpath/libmlpack_julia_krann.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_lars, "@rpath/libmlpack_julia_lars.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_linear_regression, "@rpath/libmlpack_julia_linear_regression.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_linear_svm, "@rpath/libmlpack_julia_linear_svm.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_lmnn, "@rpath/libmlpack_julia_lmnn.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_local_coordinate_coding, "@rpath/libmlpack_julia_local_coordinate_coding.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_logistic_regression, "@rpath/libmlpack_julia_logistic_regression.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_lsh, "@rpath/libmlpack_julia_lsh.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_mean_shift, "@rpath/libmlpack_julia_mean_shift.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_nbc, "@rpath/libmlpack_julia_nbc.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_nca, "@rpath/libmlpack_julia_nca.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_nmf, "@rpath/libmlpack_julia_nmf.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_pca, "@rpath/libmlpack_julia_pca.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_perceptron, "@rpath/libmlpack_julia_perceptron.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_binarize, "@rpath/libmlpack_julia_preprocess_binarize.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_describe, "@rpath/libmlpack_julia_preprocess_describe.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_one_hot_encoding, "@rpath/libmlpack_julia_preprocess_one_hot_encoding.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_scale, "@rpath/libmlpack_julia_preprocess_scale.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_split, "@rpath/libmlpack_julia_preprocess_split.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_radical, "@rpath/libmlpack_julia_radical.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_random_forest, "@rpath/libmlpack_julia_random_forest.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_softmax_regression, "@rpath/libmlpack_julia_softmax_regression.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_sparse_coding, "@rpath/libmlpack_julia_sparse_coding.dylib")
JLLWrappers.@declare_library_product(libmlpack_julia_util, "@rpath/libmlpack_julia_util.dylib")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, armadillo_jll, OpenBLAS_jll)
    JLLWrappers.@init_library_product(
        libmlpack,
        "lib/libmlpack.3.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Interact

import Interact.statedict

statedict(s::Union(Slider, Progress)) =
    {:value=>s.value,
     :min=>first(s.range),
     :step=>step(s.range),
     :max=>last(s.range)}

# when we say value to javascript, it really means value label
statedict(d::Options) =
    {:value_name=>d.value_label,
     :value_names=>collect(keys(d.options))}


text:
using Pkg
Pkg.instantiate()
Pkg.activate(".")
using Psychometrics
using Distributions
using LinearAlgebra
using Dates
using Random
using SharedArrays
using JLD2

I = 1020
I_to_calibrate = 20
N = 2_000
test_length = 35
field_items = 5
oper_items  = test_length - field_items
iter_mcmc_latent = 2_000
iter_mcmc_item = 4_000
#after how many responses update item parameter estimates
batch_size = 5
N_T = 200

@load "test/simulation/batch5    NT200 Dgain/true values/true_items_1.jld2" items
@load "test/simulation/batch5    NT200 Dgain/true values/true_examinees_1.jld2" examinees
a_s = zeros(1020,100)
b_s = zeros(1020,100)
for rep = 1:100
    @load string("test/simulation/batch5    NT200 Dgain/rep_", rep, "_items_est.jld2") items_est
    a = map(i -> i.parameters.a, items_est)
    a_s[:, rep] = copy(a)
    b = map(i -> i.parameters.b, items_est)
    b_s[:, rep] = copy(b)
end

a_true = map(i -> i.parameters.a, items)
a_true = hcat([a_true for rep=1:100]...)
b_true = map(i -> i.parameters.b, items)
b_true = hcat([b_true for rep=1:100]...)
diff_a = (a_s[1001:1020,:] - a_true[1001:1020,:]).^2
diff_b = (b_s[1001:1020,:] - b_true[1001:1020,:]).^2

println(sqrt(sum(diff_a)/20/100))
println(sqrt(sum(diff_b)/20/100))

#examinees

@load "test/simulation/batch5    NT200 Dgain/true values/true_items_1.jld2" items
@load "test/simulation/batch5    NT200 Dgain/true values/true_examinees_1.jld2" examinees
theta_s = zeros(2000,100)
for rep = 1:100
    @load string("test/simulation/batch5    NT200 Dgain/examinees_est_rep_", rep, ".jld2") examinees_est
    t = map(e -> e.latent.val, examinees_est)
    theta_s[:, rep] = copy(t)
end

t_true = map(e -> e.latent.val, examinees)
t_true = hcat([t_true for rep=1:100]...)
diff_t = (theta_s - t_true).^2

println(sqrt(sum(diff_t[1:800])/800/100))

@load "test/simulation/batch5    NT200 Dgain/examinees_est_theta_rep_1.jld2" examinees_est_theta

rmses=Matrix{Float64}(undef, test_length, N)
for n=1:N
    rmses[:,n] = sqrt.(((examinees_est_theta[n].-examinees[n].latent.val)).^2)
end

plot(1:test_length, rmses[:,1:20], legend=false)

rmses_k=Matrix{Float64}(undef, 35, size(collect(-3.0:1:3.0),1))
thetas = map(e -> e.latent.val, examinees)[1:700]
thetas_k = collect(-3.0:1:3.0)
for k=2:size(thetas_k, 1)
    rmses_k[:,k] = mean(rmses[:, findall((thetas .< thetas_k[k]) .& (thetas .>= thetas_k[max(1,k-1)]))], dims=2)[:,1]
end


text:
using Polyhedra
using Base.Test
using CDDLib
using GeometryTypes

ext1 = SimpleVRepresentation([0 1 2 3;0 2 1 3; 1 0 2 3; 1 2 0 3; 2 0 1 3; 2 1 0 3;
                              0 1 3 2;0 3 1 2; 1 0 3 2; 1 3 0 2; 3 0 1 2; 3 1 0 2;
                              0 3 2 1;0 2 3 1; 3 0 2 1; 3 2 0 1; 2 0 3 1; 2 3 0 1;
                              3 1 2 0;3 2 1 0; 1 3 2 0; 1 2 3 0; 2 3 1 0; 2 1 3 0])
poly1 = CDDPolyhedron{4,Rational{BigInt}}(ext1)

poly2 = project(poly1, [1 1 1; -1 1 1; 0 -2 1; 0 0 -3])


text:
add_axis(x) = reshape(x, size(x)..., 1)

function get_subdir_path(sub_dir::String)
    # Get the directory this module is located in
    abs_path_module = realpath(@__DIR__)
    module_dir, _ = splitdir(abs_path_module)

    joinpath(module_dir, sub_dir)
end

function get_file_path(sub_dir::String, file_name::String, default_ext::String)
    #=
    Get the absolute path of a resource file, which may be relative to
    the Duckietown.jl module directory, or an absolute path.

    This function is necessary because the simulator may be imported by
    other packages, and we need to be able to load resources no matter
    what the current working directory is.
    =#

    @assert '.' ∉ default_ext
    @assert '/' ∉ default_ext

    # If this is already a real path
    ispath(file_name) && (return file_name)

    subdir_path = get_subdir_path(sub_dir)
    file_path = joinpath(subdir_path, file_name)

    if '.' ∉ file_name
        file_path *= "." * default_ext
    end

    return file_path
end


text:
using Documenter, ACME

makedocs(
    modules = [ACME],
    sitename = "ACME.jl",
    pages = Any[
        "Home" => "index.md",
        "gettingstarted.md",
        "ug.md",
        "elements.md",
    ],
)

if v"1.0" ≤ VERSION < v"1.1" # deploy from 1.0
    deploydocs(
        repo = "github.com/HSU-ANT/ACME.jl.git",
        devbranch = "main",
    )
end


text:
part1(input) = join(next_password(collect(readchomp(input))))
part2(input) = join(next_password(next_password(collect(readchomp(input)))))

function next_password(password)
    next_string!(password)
    while !is_valid(password)
        next_string!(password)
    end
    return password
end

function is_valid(password)
    'i' in password && return false
    'l' in password && return false
    'o' in password && return false
    has_triple = any(password[i + 1] == password[i] + 1 && password[i + 2] == password[i] + 2 for i = 1:6)
    has_triple || return false
    i = 1
    pairs = 0
    while i <= 7
        if password[i] == password[i + 1]
            pairs += 1
            i += 2
        else
            i += 1
        end
    end
    return pairs >= 2
end

function next_string!(password)
    for i = 8:-1:1
        password[i] += 1
        if password[i] > 'z'
            password[i] = 'a'
        else
            break
        end
    end
end


text:
using GraphNeuralNetworks
using Statistics

"""
    GcnHP

Hyperparameters for the gcn architecture.

"""

@kwdef struct GcnHP
  depth_common :: Int = 6
  depth_phead :: Int = 2
  depth_vhead :: Int = 2
  hidden_size :: Int = 64
end

"""
    Gin <: TwoHeadGraphNeuralNetwork

A simple two-headed GNN architecture with only dense layers.
"""
mutable struct Gcn <: TwoHeadGraphNeuralNetwork
  gspec
  hyper
  common
  vhead
  phead
end

function Gcn(gspec::AbstractGameSpec, hyper::GcnHP)
  indim, n_nodes  = GI.state_dim(gspec)
  GCN_layers(depth) = [GNNChain(GCNConv(((i==1) ? indim : hyper.hidden_size) => hyper.hidden_size, relu, add_self_loops=true),
                                BatchNorm(hyper.hidden_size, relu)) for i in 1:depth]
  Dense_layers(size, depth) = [Chain(Dense(size, size, relu), BatchNorm(size)) for _ in 1:depth]
  common = GNNChain(GCN_layers(hyper.depth_common)..., BatchNorm(hyper.hidden_size))
  vhead = GNNChain(GlobalPool(mean),
                   Dense_layers(hyper.hidden_size, hyper.depth_vhead)...,
                   Dense(hyper.hidden_size, 1, identity))
  phead = GNNChain(Parallel(vcat, GNNChain(GlobalPool(mean), x -> repeat(x, 1, n_nodes)), identity),
                  Dense_layers(hyper.hidden_size*2, hyper.depth_phead)...,
                  Dense(hyper.hidden_size*2, 1),
                  x->reshape(x, n_nodes,:),
                  softmax)
  return Gcn(gspec, hyper, common, vhead, phead)
end

Network.HyperParams(::Type{Gcn}) = GcnHP

function Base.copy(nn::Gcn)
  return Gcn(
    nn.gspec,
    nn.hyper,
    deepcopy(nn.common),
    deepcopy(nn.vhead),
    deepcopy(nn.phead)
  )
end


text:
"""
    Options(;
        radius = .1,
        bounds,
        n_iters,
        init_parms,
        parallel = false,
        adapt_radius! = adapt!,
        kwargs...
    )

An object that holds configuration options for the MCMC sampler.

# Fields

- `parallel=false`: runs model on threads if true. A speed up is observed if the evaluation
time of the function is 1 ms or greater.
- `radius = .10`: the initial radius length for each chain
- `bounds`: a vector of tuples representing (lowerbound, upperbound) for each dimension in
the parameter space
- `x_range`: the range of allowable values for each parameter
- `n_iters`: number of iterations to perform
- `p_eval`: the function that evalues the model and pattern functions
- `adapt_radius!=adapt!`: a function in the form of `func(chain, options; kwargs...)` that adapts
the radius.
- `init_parms`: a vector of starting points, such as [[.3,.4],[.3,.5]] in the 2 dimensional case.
"""
@concrete struct Options
    parallel
    radius
    bounds
    x_range
    n_iters
    p_eval
    adapt_radius!
    init_parms
    n_dims
end

function Options(;
        radius = .10,
        bounds,
        n_iters,
        init_parms,
        parallel = false,
        adapt_radius! = adapt!,
        kwargs...
    )

    p_eval = parallel ? t_eval_patterns : eval_patterns
    _adapt_radius! = (x,y) -> adapt_radius!(x, y; kwargs...)
    x_range = map(x -> x[2] - x[1], bounds)
    n_dims = length(bounds)

    return Options(
        parallel,
        radius,
        bounds,
        x_range,
        n_iters,
        p_eval,
        _adapt_radius!,
        init_parms,
        n_dims
    )
end

Broadcast.broadcastable(x::Options) = Ref(x)

"""
    Chain(;parms, pattern, radius)

An MCMC chain object.

# Arguments

- `parms`: a vector of parameters, i.e the current state of the chain
- `pattern`: the target pattern of the chain
- `radius`: the radius for the jump distribution

# Fields

- `parms`: a vector of parameters, i.e the current state of the chain
- `n_dims`: the number of dimensions
- `pattern`: the target pattern of the chain
- `radius`: the radius for the jump distribution
- `acceptance`: a Boolean vector indicating whether a proposal was accepted
- `init_parms`: a vector of starting points, such as [[.3,.4],[.3,.5]] in the 2 dimensional case.
"""
@concrete mutable struct Chain
    parms
    n_dims
    pattern
    radius
    acceptance
end

function Chain(parms, pattern, radius)
    return Chain(
        parms,
        length(parms),
        pattern,
        radius,
        [true]
    )
end

"""
    Results(chain_id, chain, iter)

Stores results for parameter space partitioning.

# Arguments

- `chain_id`: integer index for chain
- `chain`: chain object
- `iter`: current iteration

An MCMC chain object.

# Fields

- `iter`: current iteration
- `chain_id`: integer index for chain
- `parms`: parameter vector
- `pattern`: target pattern of chain
- `acceptance`: a Boolean vector indicating whether a proposal was accepted
"""
@concrete struct Results
    iter
    chain_id
    parms
    pattern
    acceptance
end

function Results(chain_id, chain, iter)
    acceptance = isempty(chain.acceptance) ? false : chain.acceptance[end]
    return Results(
        iter,
        chain_id,
        chain.parms,
        chain.pattern,
        acceptance
    )
end

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


text:
module TrafficAssignment

# package code goes here
using LightGraphs, Optim, BinDeps
using Distributed, Printf, LinearAlgebra, SparseArrays

TNTP_SHA = "ca933f74c03af688dfffa33f4b00826cba672479"

include("load_network.jl")
include("frank_wolfe.jl")

export
        load_ta_network, download_tntp, read_ta_network,
        ta_frank_wolfe,
        TA_Data

end # module


text:
# -----------------------------------------------------------------------
# Test set for critical functions
# -----------------------------------------------------------------------
using Test
include("GA.jl")

# Test set for stepper movements
function coordsteps_test()
    @testset "Coordinates and steps" begin
        @testset "One cylinder" begin
            p1 = OneCylType(0.0,0.0)
            p2 = OneCylType(10.0,0.0)
            p3 = OneCylType(20.0,0.0)
            p2_cor = steps2coord(StepperMotionType([41,16],[0]),xorig=[-8.8],d=2.4, Rh=80.0, Re=10.0)
            p3_cor = steps2coord(StepperMotionType([82,-10],[0]),xorig=[-8.8],d=2.4, Rh=80.0, Re=10.0)
            @test coord2steps(p1,p2,xorig=[-8.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [41,16] atol=1
            @test coord2steps(p1,p3,xorig=[-8.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [82,-10] atol=1
            @test coord2steps(p2,p3,xorig=[-8.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [41,-25] atol=1
            @test coord2steps(p1,p3,xorig=[-8.8],d=2.4, Rh=80.0, Re=10.0)[1].strategy == [1,2]
            @test coord2steps(p3,p2,xorig=[-8.8],d=2.4, Rh=80.0, Re=10.0)[1].strategy == [1,2]
            @test convert(Array, p2_cor) ≈ convert(Array, p2) atol=0.2
            @test convert(Array, p3_cor) ≈ convert(Array, p3) atol=0.2
        end
        @testset "Two cylinders" begin
            p1 = TwoCylType(0.0,0.0,2.4,0.0)
            p2 = TwoCylType(10.0,0.0,2.4,0.0)
            p3 = TwoCylType(20.0,0.0,2.4,0.0)
            p2_cor = steps2coord(StepperMotionType([41,16,41,47],[0]),xorig=[-8.8,-18.8],d=2.4, Rh=80.0, Re=10.0)
            p3_cor = steps2coord(StepperMotionType([82,-10,82,53],[0]),xorig=[-8.8,-18.8],d=2.4, Rh=80.0, Re=10.0)
            @test coord2steps(p1,p2,xorig=[-8.8,-18.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [41,16,41,47] atol=1
            @test coord2steps(p1,p3,xorig=[-8.8,-18.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [82,-10,82,53] atol=1
            @test coord2steps(p2,p3,xorig=[-8.8,-18.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [41,-25,41,6] atol=1
            @test coord2steps(p1,p3,xorig=[-8.8,-18.8],d=2.4, Rh=80.0, Re=10.0)[1].strategy == [3,4,1,2]
            @test coord2steps(p3,p2,xorig=[-8.8,-18.8],d=2.4, Rh=80.0, Re=10.0)[1].strategy == [1,2,3,4]
            @test convert(Array, p2_cor) ≈ convert(Array, p2) atol=0.2
            @test convert(Array, p3_cor) ≈ convert(Array, p3) atol=0.2
        end
        @testset "Three cyliders" begin
            p1 = ThreeCylType(0.0,0.0,2.4,0.0,2.4,0.0)
            p2 = ThreeCylType(10.0,0.0,2.4,0.0,2.4,0.0)
            p3 = ThreeCylType(20.0,0.0,2.4,0.0,2.4,0.0)
            p2_cor = steps2coord(StepperMotionType([41,16,41,47,42,80],[0]),xorig=[-8.8,-18.8,-28.8],d=2.4, Rh=80.0, Re=10.0)
            p3_cor = steps2coord(StepperMotionType([82,-10,82,53,83,117],[0]),xorig=[-8.8,-18.8,-28.8],d=2.4,Rh=80.0, Re=10.0)
            @test coord2steps(p1,p2,xorig=[-8.8,-18.8,-28.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [41,16,41,47,42,80] atol=1
            @test coord2steps(p1,p3,xorig=[-8.8,-18.8,-28.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [82,-10,82,53,83,117] atol=1
            @test coord2steps(p2,p3,xorig=[-8.8,-18.8,-28.8],d=2.4, Rh=80.0, Re=10.0)[1].step ≈ [41,-25,41,6,41,37] atol=1
            @test coord2steps(p1,p3,xorig=[-8.8,-18.8,-28.8],d=2.4, Rh=80.0, Re=10.0)[1].strategy ≈ [5,6,3,4,1,2]
            @test coord2steps(p3,p2,xorig=[-8.8,-18.8,-28.8],d=2.4, Rh=80.0, Re=10.0)[1].strategy ≈ [1,2,3,4,5,6]
            @test convert(Array, p2_cor) ≈ convert(Array, p2) atol=0.2
            @test convert(Array, p3_cor) ≈ convert(Array, p3) atol=0.2
        end
    end
end

function bitsconversion_tests()
    @testset "Bits and coordinates transformations" begin

        c11 = OneCylType(1.0, 2.0)
        c11_bit = cylpos2bits(c11, size=[6,5])
        @test bits2cylpos(c11_bit) == c11
        c12 = OneCylType(1.0, -2.0)
        c12_bit = cylpos2bits(c12, size=[6,5])
        @test bits2cylpos(c12_bit) == c12

        c21 = TwoCylType(2.0,1.0,3.0,4.0)
        c21_bit = cylpos2bits(c21, size=[6,5,4,4])
        @test bits2cylpos(c21_bit) == c21
        c22 = TwoCylType(2.0,-1.0,3.0,-4.0)
        c22_bit = cylpos2bits(c22, size=[6,5,4,4])
        @test bits2cylpos(c22_bit) == c22

        c31 = ThreeCylType(10.0,2.0,2.0,3.0,2.0,1.0)
        c31_bit = cylpos2bits(c31, size=[6,5,4,4,4,4])
        @test bits2cylpos(c31_bit) == c31
        c32 = ThreeCylType(10.0,-2.0,2.0,-3.0,2.0,-1.0)
        c32_bit = cylpos2bits(c32, size=[6,5,4,4,4,4])
        @test bits2cylpos(c32_bit) == c32
        c33 = ThreeCylType(0,0,0,0,0,0)
        c33_bit = cylpos2bits(c33, size=[6,5,4,4,4,4])
        @test bits2cylpos(c33_bit) == c33

    end
end

# TEST ALL
function testall()
    coordsteps_test()
    bitsconversion_tests()
end


text:
import Base: size, length, show, iterate, lastindex, isempty, convert, ndims, float, eltype, copy, ==
if VERSION >= v"0.7-"
    using Dates
else
    using Base.Dates
end

################################################################################
# TYPE DEFINITION ##############################################################
################################################################################

abstract type AbstractTS end

"""
Time series type aimed at efficiency and simplicity.

Motivated by the `xts` package in R and the `pandas` package in Python.
"""
mutable struct TS{V<:Real,T<:TimeType}
    values::Matrix{V}
    index::Vector{T}
    fields::Vector{Symbol}
    function TS{V,T}(values::AbstractArray{V}, index::AbstractVector{T}, fields::Vector{Symbol}) where {V<:Real,T<:TimeType}
        @assert size(values,1)==length(index) "Length of index not equal to number of value rows."
        @assert size(values,2)==length(fields) "Length of fields not equal to number of columns in values."
        order = sortperm(index)
        return new(values[order,:], index[order], SANITIZE_NAMES ? namefix.(fields) : fields)
    end
end

TS(v::AbstractArray{V}, t::AbstractVector{T}, f::Union{Symbol,String,Char}) where {V,T} = TS{V,T}(v, t, [Symbol(f)])
TS(v::AbstractArray{V}, t::AbstractVector{T}, f::Union{AbstractVector{Symbol},AbstractVector{String},AbstractVector{Char}}) where {V,T} = TS{V,T}(v, t, Symbol.(f))
TS(v::AbstractArray{V}, t::AbstractVector{T}) where {V,T} = TS{V,T}(v, t, autocol(1:size(v,2)))
TS(v::AbstractArray{V}, t::T, f) where {V,T} = TS{V,T}(v, [t], f)
TS(v::V, t::AbstractVector{T}, f) where {V,T} = TS{V,T}([v], t, f)
TS(v::V, t::T, f) where {V,T} = TS{V,T}([v][:,:], [t], f)
TS(v::V, t::T) where {V,T} = TS{V,T}([v], [t], [:A])
TS(v::AbstractArray{V}) where {V} = TS{V,Date}(v, autoidx(size(v,1)), autocol(1:size(v,2)))
TS() = TS{Float64,Date}(Matrix{Float64}(UndefInitializer(),0,0), Date[], Symbol[])
TS(X::TS) = TS(X.values, X.index, X.fields)

# Conversions ------------------------------------------------------------------
convert(::Type{TS{Float64}}, x::TS{Bool}) = TS{Float64}(map(Float64, x.values), x.index, x.fields)
convert(::Type{TS{Int}}, x::TS{Bool}) = TS{Int}(map(Int, x.values), x.index, x.fields)
convert(::Type{TS{Bool}}, x::TS{V}) where {V<:Real} = TS{Bool}(map(V, x.values), x.index, x.fields)
convert(x::TS{Bool}) = convert(TS{Int}, x::TS{Bool})
const ts = TS

################################################################################
# BASIC UTILITIES ##############################################################
################################################################################
size(x::TS) = size(x.values)
size(x::TS, dim::Int) = size(x.values, dim)
ndims(::Type{TS{V,T}}) where {V,T} = 2
length(x::TS) = prod(size(x))::Int
lastindex(x::TS) = lastindex(x.values)
lastindex(x::TS, d) = lastindex(x.values, d)
iterate(x::TS) = size(x,1) == 0 ? nothing : (x.index[1], x.values[1,:]), 2
iterate(x::TS, i::Int) = i == lastindex(x, 1) + 1 ? nothing : ((x.index[i], x.values[i,:]), i+1)
isempty(x::TS) = (isempty(x.index) && isempty(x.values))
first(x::TS) = x[1]
last(x::TS) = x[end]
ndims(::TS) = 2
float(x::TS) = ts(float(x.values), x.index, x.fields)
eltype(x::TS) = eltype(x.values)
copy(x::TS) = TS(x.values, x.index, x.fields)

################################################################################
# RENAME #######################################################################
################################################################################
# in place
function rename!(ts::TS, args::Pair{Symbol, Symbol}...)
    d = Dict{Symbol, Symbol}(args...)
    flag = false
    for (i, field) in enumerate(ts.fields)
        if field in keys(d)
            ts.fields[i] = d[field]
            flag = true
        end
    end
    flag
end

function rename!(f::Base.Callable, ts::TS, colnametyp::Type{Symbol} = Symbol)
    for (i, field) in enumerate(ts.fields)
        ts.fields[i] = f(field)
    end
    true
end

function rename!(f::Base.Callable, ts::TS, colnametyp::Type{String})
    f = Symbol ∘ f ∘ string
    rename!(f, ts)
end

# not in place
function rename(ts::TS, args::Pair{Symbol, Symbol}...)
    ts2 = copy(ts)
    rename!(ts2, args...)
    ts2
end

function rename(f::Base.Callable, ts::TS, colnametyp::Type{Symbol} = Symbol)
    ts2 = copy(ts)
    rename!(f, ts2, colnametyp)
    ts2
end

function rename(f::Base.Callable, ts::TS, colnametyp::Type{String})
    ts2 = copy(ts)
    rename!(f, ts2, colnametyp)
    ts2
end


text:
"""
This struct is passed to `mysolve`.

fields:
* alg
* f!
* u0
* tspan
* dt
* p
"""
struct MyProbInplace
    alg!
    f!
    u0
    tspan
    dt
    p
end

"""
	mysolve_inplace(prob::MyProbInplace)
"""
function mysolve_inplace(prob::MyProbInplace)
    m = Matrix{eltype(prob.u0)}(undef, length(prob.u0), Int((prob.tspan[2] - prob.tspan[1]) / prob.dt))
    prob.alg!(m, prob.f!, prob.u0, prob.tspan, prob.dt, prob.p)
end

"""
	mysolve(alg, f, u0, tspan, dt; p=nothing)

Allocating solve.
The difference between `mysolve` and `mysolve_inplace` is
`mysolve` expects `f` to return an array, not mutate one.
"""
function mysolve(alg, f, u0, tspan, dt; p=nothing)
    m = Matrix{eltype(u0)}(undef, length(u0), Int((tspan[2] - tspan[1]) / dt))
    alg(m, f, u0, tspan, dt, p)
end

"""
	mysolve_inplace(alg!, f!, u0, tspan, dt; p=nothing)

Non-allocating solve.
The difference between `mysolve` and `mysolve_inplace` is
`mysolve` expects `f` to return an array, not mutate one.
"""
function mysolve_inplace(alg!, f!, u0, tspan, dt; p=nothing)
    m = Matrix{eltype(u0)}(undef, length(u0), Int((tspan[2] - tspan[1]) / dt))
    alg!(m, f!, u0, tspan, dt, p)
end

"""
	eulers_inplace_f!(m, f!, u0, tspan, dt, p)

Non-allocating explicit first order method. Fixed dt.
"""
function eulers_inplace_f!(m, f!, u0, tspan, dt, p)
    du = similar(u0)
    m[:, 1] = 	u0
    for i in 2:size(m, 2)
        t = (i - 1) * dt + tspan[1]
        f!(du, m[:, i - 1],	 p, t)
        m[:, i] = m[:, i - 1] + dt * du
    end
    m
end

"""
	midpoint!(m, f, u0, tspan, dt, p)

Explicit second order method. Fixed dt.
Allocates for intermediate `k` values.
This means `f` does not return `nothing`.
"""
function midpoint!(m, f, u0, tspan, dt, p)
# du = similar(u0)
    m[:, 1] = u0
    for i in 2:size(m, 2)
        t = (i - 1) * dt + tspan[1]
        k1 = f(m[:, i - 1], p, t)
        k2 = f(m[:, i - 1] + (dt / 2) * k1, p, t + dt / 2)
        m[:, i] = m[:, i - 1] + dt * k2
    end
    m
end

"""
	midpoint_inplace_f!(m, f, u0, tspan, dt, p)

Explicit second order method. Fixed dt.

TODO: have du allocated beforehand, so that this function
has zero allocations.
"""
function midpoint_inplace_f!(m, f!, u0, tspan, dt, p)
    du = similar(u0)
    m[:, 1] = u0
    for i in 2:size(m, 2)
        t = (i - 1) * dt + tspan[1]
        f!(du, m[:, i - 1], p, t)
        f!(du, m[:, i - 1] + (dt / 2) * du, p, t + dt / 2)
        m[:, i] = m[:, i - 1] + dt * du
    end
    m
end

"""
	rk4!(m, f, u0, tspan, dt, p)

Explicit fourth-order Runge-Kutta method. Fixed dt.
`m` is mutated inplace, but `f` is not an in-place function.
"""
function rk4!(m, f, u0, tspan, dt, p)
    m[:, 1] = u0
    for i in 2:size(m, 2)
        u = m[:, i - 1]
        t = (i - 1) * dt + tspan[1]

        k1 = f(u, p, t)
        k2 = f(u + (dt / 2) * k1, p, t + dt / 2)
        k3 = f(u + (dt / 2) * k2, p, t + dt / 2)
        k4 = f(u + dt * k3, p, t + dt)

        m[:, i] = u + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6
    end
    m
end

"""
	rk4_inplace_f!(m, f!, u0, tspan, dt, p)

Explicit fourth-order Runge-Kutta method. Fixed dt.
`m` is mutated inplace and `f` is an in-place function.
"""
function rk4_inplace_f!(m, f!, u0, tspan, dt, p)

	du = similar(u0)
	m[:, 1] = u0
	# k = Matrix{eltype(u0)}(undef, length(u0), 3)

	for i in 2:size(m, 2) # what is the better way to write this, prealloc ks and overwrite?
		u = m[:, i-1]
		t = (i-1)*dt + tspan[1]
		f!(du, u, p, t)
		k1 = copy(du)
		f!(du, u + (dt/2)*k1, p, t + dt/2)
		k2 = copy(du)
		f!(du, u + (dt/2)*k2, p, t + dt/2)
		k3 = copy(du)
		f!(du, u + dt*k3, p, t + dt)

		m[:, i] = u + dt*(k1 + 2*k2 + 2*k3 + du)/6
	end
	m
end

"""
	rk4_inplace_f_prealloc_k!(m, f!, u0, tspan, dt, p)

Explicit fourth-order Runge-Kutta method. Fixed dt.
`m` is mutated inplace and `f` is an in-place function.

The difference with the other `rk4` methods is that
the `k` values are allocated before the loop.

I want to see which of these is fastest.
"""
function rk4_inplace_f_prealloc_k!(m, f!, u0, tspan, dt, p)

    du = similar(u0)
    m[:, 1] = u0
    k = Matrix{eltype(u0)}(undef, length(u0), 3)

    for i in 2:size(m, 2)
        u = m[:, i - 1]
        t = (i - 1) * dt + tspan[1]
        f!(du, u, p, t)
        k[:, 1] = copy(du)
        f!(du, u + (dt / 2) * k[:, 1], p, t + dt / 2)
        k[:, 2] = copy(du)
        f!(du, u + (dt / 2) * k[:, 2], p, t + dt / 2)
        k[:, 3] = copy(du)
        f!(du, u + dt * k[:, 3], p, t + dt)

        m[:, i] = u + dt * (k[:, 1] + 2 * k[:, 2] + 2 * k[:, 3] + du) / 6
    end
    m
end

"Dummy struct for the Dormand-Prince tableau"
struct DP5Cache
    a
    b
    c
end

"""
	dp5_cache()

allocates the DP5 tableau.
"""
function dp5_cache()
    dp_a_1 = [0]
    dp_a_2 = [1 / 5]
    dp_a_3 = [3 / 40, 9 / 40]
    dp_a_4 = [44 / 45, -56 / 15, 32 / 9]
    dp_a_5 = [19372 / 6561, -25360 / 2187, 64448 / 6561, -212 / 729]
    dp_a_6 = [9017 / 3168, -355 / 33, 46732 / 5247, 49 / 176, -5103 / 18656]
    dp_a_7 = [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84]
    a = [dp_a_1, dp_a_2, dp_a_3, dp_a_4, dp_a_5, dp_a_6, dp_a_7]
    b = [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84, 0]
    c = [0, 1 / 5, 3 / 10, 4 / 5, 8 / 9, 1, 1]
    DP5Cache(a, b, c)
end

"""
	dp5(f, u0, tspan, dt, p)

Uses `dp5!` to solve f and store in `m`.
"""
function dp5(f, u0, tspan, dt, p)
    m = Matrix{eltype(u0)}(undef, length(u0), Int((tspan[2] - tspan[1]) / dt))
    dp5!(m, f, u0, tspan, dt, p)
end

"""
	dp5!(m, f, u0, tspan, dt, p)

Explicit Dormand-Prince method. Fixed dt.
Non-allocating for the solves.
Allocates the intermediate `f` calls though.
"""
function dp5!(m, f, u0, tspan, dt, p)
    cache = dp5_cache() # allocates, inefficient
    a, b, c = cache.a, cache.b, cache.c
    S = length(c) # number of stages
    N = length(u0)
    m[:, 1] = u0
    K = Matrix{eltype(u0)}(undef, N, S)
    for i in 2:size(m, 2)
        u = m[:, i - 1]
        t = (i - 1) * dt + tspan[1]
        K[:, 1] = f(u, p, t)
        for j in 2:S # calculate K vals
            tmp = sum([a[j][k] * K[:, k] for k in 1:j - 1])
            K[:, j] = f(u + dt * tmp, p, t + dt * c[j])
        end
        m[:, i] = m[:, i - 1] + dt * sum([b[k] .* K[:, k] for k in 1:S])
    end
    m
end

"""
	dp5_inplace!(m, f!, u0, tspan, dt, p)

Explicit Dormand-Prince method. Fixed dt.

"""
function dp5_inplace!(m, f!, u0, tspan, dt, p; cache=dp5_cache())
    # cache = dp5_cache() # allocates, inefficient
    a, b, c = cache.a, cache.b, cache.c
    S = length(c) # number of stages
    N = length(u0)
    m[:, 1] = u0
    K = Matrix{eltype(u0)}(undef, N, S)
    for i in 2:size(m, 2)
        u = m[:, i - 1]
        t = (i - 1) * dt + tspan[1]
        # K[:, 1] = f(u, p, t)
        f!(K[:, 1], u, p, t)
        for j in 2:S # calculate K vals
            tmp = sum([a[j][k] * K[:, k] for k in 1:j - 1])
            f!(K[:, j], u + dt * tmp, p, t + dt * c[j])
        end
        m[:, i] = m[:, i - 1] + dt * sum([b[k] .* K[:, k] for k in 1:S])
    end
    m
end


text:
"""
    guesswork_lower_bound(
        p::AbstractVector{T},
        ρBs::AbstractVector{<:AbstractMatrix};
        solver,
        c = T[1:length(p)..., 10_000],
        verbose::Bool = false,
    )

See [`guesswork`](@ref) for the meaning of the arguments. Computes a lower bound
to the optimal expected number of guesses by solving a relaxed version of the
primal SDP. For `J` states, only needs `J^2` PSD variables subject to two linear
constraints.
"""
function guesswork_lower_bound(
    p::AbstractVector{T},
    ρBs::AbstractVector{<:AbstractMatrix};
    solver,
    c = T[1:length(p)..., 10_000],
    verbose::Bool = false,
) where {T<:Number}
    length(p) == length(ρBs) ||
    throw(ArgumentError("Length of prior and vector of side information must match J"))
    J = length(p)
    dB = size(ρBs[1], 1)
    all(ρB -> size(ρB) == (dB, dB), ρBs) ||
    throw(ArgumentError("All side-information states must be square matrices of the same size"))

    ℰs = [ComplexVariable(dB, dB) for j = 1:J, k = 1:J]

    objective = real(sum(c[j] * tr(ℰs[j, k] * p[k] * ρBs[k]) for j = 1:J for k = 1:J))

    constraints = Convex.Constraint[ℰ ⪰ 0 for ℰ in ℰs] |> vec

    for j = 1:J
        push!(constraints, sum(ℰs[j, :]) == I(dB))
        push!(constraints, sum(ℰs[:, j]) == I(dB))
    end

    problem = compatible_problem(Convex.minimize, objective, constraints, T)

    solve!(problem, solver; verbose = verbose)

    ℰs = evaluate.(ℰs)

    data = (J = J, K = J, c = c, p = p, ρBs = ρBs, dB = dB)

    return (optval = problem.optval, status = problem.status, ℰs = ℰs, data...)
end


text:
using Test
using IsConstFoldable

# Simple functions with constant arguments should definitely fold away!
@testset "simple functions const args should fold" begin
    @test @is_const_foldable 5+5
    @test @is_const_foldable zero(Int)
end

# An operation that has non-const arguments cannot fold away.
@testset "nonconst args shouldn't fold" begin
    # Global args aren't const
    global x,y = 3,2
    @test !@is_const_foldable x+y
end

@testset "called from inside functions" begin
    # Can be called from inside functions (typically via an @assert)
    foo(x) = @is_const_foldable typemax(x)
    @test foo(x)

    # Operations on runtime argument will not const-fold.
    foo(x) = @is_const_foldable x+1
    @test !foo(x)
end

@testset "empty functions" begin
    function empty_f() end
    @test @is_const_foldable empty_f()
end

# Whether g() will const-fold depends entirely on the function's body.
@generated function g(x)
    if x <: Int
        # Will const-fold
        return :(typemax(x))
    else
        # Won't const-fold (because relies on runtime values)
        return :(x+x)
    end
end

@testset "@generated functions" begin
    @assert g(0) == typemax(Int)
    @assert g(1.0) == 2.0

    # Use a runtime-argument to prevent const-folding of complex expressions.
    f(x) = @is_const_foldable g(x)
    # The foldable case
    @test f(0)
    # The non-foldable case
    @test !f(1.0)
end


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

module Grisu

export print_shortest
export DIGITS, grisu

const SHORTEST = 1
const FIXED = 2
const PRECISION = 3

const DIGITS = Vector{UInt8}(uninitialized, 309+17)

include(joinpath("grisu", "float.jl"))
include(joinpath("grisu", "fastshortest.jl"))
include(joinpath("grisu", "fastprecision.jl"))
include(joinpath("grisu", "fastfixed.jl"))
include(joinpath("grisu", "bignums.jl"))
include(joinpath("grisu", "bignum.jl"))

const BIGNUMS = [Bignums.Bignum(),Bignums.Bignum(),Bignums.Bignum(),Bignums.Bignum()]

function grisu(v::AbstractFloat,mode,requested_digits,buffer=DIGITS,bignums=BIGNUMS)
    if signbit(v)
        neg = true
        v = -v
    else
        neg = false
    end
    if mode == PRECISION && requested_digits == 0
        buffer[1] = 0x00
        len = 0
        return 0, 0, neg
    end
    if v == 0.0
        buffer[1] = 0x30
        buffer[2] = 0x00
        len = point = 1
        return len, point, neg
    end
    if mode == SHORTEST
        status,len,point = fastshortest(v,buffer)
    elseif mode == FIXED
        status,len,point = fastfixedtoa(v,0,requested_digits,buffer)
    elseif mode == PRECISION
        status,len,point = fastprecision(v,requested_digits,buffer)
    end
    status && return len-1, point, neg
    status, len, point = bignumdtoa(v,mode,requested_digits,buffer,bignums)
    return len-1, point, neg
end

nanstr(x::AbstractFloat) = "NaN"
nanstr(x::Float32) = "NaN32"
nanstr(x::Float16) = "NaN16"
infstr(x::AbstractFloat) = "Inf"
infstr(x::Float32) = "Inf32"
infstr(x::Float16) = "Inf16"

function _show(io::IO, x::AbstractFloat, mode, n::Int, typed, compact)
    isnan(x) && return write(io, typed ? nanstr(x) : "NaN")
    if isinf(x)
        signbit(x) && write(io,'-')
        write(io, typed ? infstr(x) : "Inf")
        return
    end
    typed && isa(x,Float16) && write(io, "Float16(")
    (len,pt,neg),buffer = grisu(x,mode,n),DIGITS
    pdigits = pointer(buffer)
    if mode == PRECISION
        while len > 1 && buffer[len] == 0x30
            len -= 1
        end
    end
    neg && write(io,'-')
    exp_form = pt <= -4 || pt > 6
    exp_form = exp_form || (pt >= len && abs(mod(x + 0.05, 10^(pt - len)) - 0.05) > 0.05) # see issue #6608
    if exp_form # .00001 to 100000.
        # => #.#######e###
        unsafe_write(io, pdigits, 1)
        write(io, '.')
        if len > 1
            unsafe_write(io, pdigits+1, len-1)
        else
            write(io, '0')
        end
        write(io, (typed && isa(x,Float32)) ? 'f' : 'e')
        write(io, dec(pt-1))
        typed && isa(x,Float16) && write(io, ")")
        return
    elseif pt <= 0
        # => 0.00########
        write(io, "0.")
        while pt < 0
            write(io, '0')
            pt += 1
        end
        unsafe_write(io, pdigits, len)
    elseif pt >= len
        # => ########00.0
        unsafe_write(io, pdigits, len)
        while pt > len
            write(io, '0')
            len += 1
        end
        write(io, ".0")
    else # => ####.####
        unsafe_write(io, pdigits, pt)
        write(io, '.')
        unsafe_write(io, pdigits+pt, len-pt)
    end
    typed && !compact && isa(x,Float32) && write(io, "f0")
    typed && isa(x,Float16) && write(io, ")")
    nothing
end

function Base.show(io::IO, x::Union{Float64,Float32})
    if get(io, :compact, false)
        _show(io, x, PRECISION, 6, x isa Float64, true)
    else
        _show(io, x, SHORTEST, 0, get(io, :typeinfo, Any) !== typeof(x), false)
    end
end

function Base.show(io::IO, x::Float16)
    hastypeinfo = Float16 === get(io, :typeinfo, Any)
    # if hastypeinfo, the printing would be more compact using `SHORTEST`
    # while still retaining all the information
    # BUT: we want to print all digits in `show`, not in display, so we rely
    # on the :compact property to make the decision
    # (cf. https://github.com/JuliaLang/julia/pull/24651#issuecomment-345535687)
    if get(io, :compact, false) && !hastypeinfo
        _show(io, x, PRECISION, 5, false, true)
    else
        _show(io, x, SHORTEST, 0, !hastypeinfo, false)
    end
end

Base.print(io::IO, x::Float32) = _show(io, x, SHORTEST, 0, false, false)
Base.print(io::IO, x::Float16) = _show(io, x, SHORTEST, 0, false, false)

# normal:
#   0 < pt < len        ####.####           len+1
#   pt <= 0             0.000########       len-pt+1
#   len <= pt (dot)     ########000.        pt+1
#   len <= pt (no dot)  ########000         pt
# exponential:
#   pt <= 0             ########e-###       len+k+2
#   0 < pt              ########e###        len+k+1

function _print_shortest(io::IO, x::AbstractFloat, dot::Bool, mode, n::Int)
    isnan(x) && return write(io, "NaN")
    x < 0 && write(io,'-')
    isinf(x) && return write(io, "Inf")
    (len,pt,neg),buffer = grisu(x,mode,n),DIGITS
    pdigits = pointer(buffer)
    e = pt-len
    k = -9<=e<=9 ? 1 : 2
    if -pt > k+1 || e+dot > k+1
        # => ########e###
        unsafe_write(io, pdigits+0, len)
        write(io, 'e')
        write(io, dec(e))
        return
    elseif pt <= 0
        # => 0.000########
        write(io, "0.")
        while pt < 0
            write(io, '0')
            pt += 1
        end
        unsafe_write(io, pdigits+0, len)
    elseif e >= dot
        # => ########000.
        unsafe_write(io, pdigits+0, len)
        while e > 0
            write(io, '0')
            e -= 1
        end
        if dot
            write(io, '.')
        end
    else # => ####.####
        unsafe_write(io, pdigits+0, pt)
        write(io, '.')
        unsafe_write(io, pdigits+pt, len-pt)
    end
    nothing
end

"""
    print_shortest(io::IO, x)

Print the shortest possible representation, with the minimum number of consecutive non-zero
digits, of number `x`, ensuring that it would parse to the exact same number.
"""
print_shortest(io::IO, x::AbstractFloat, dot::Bool) = _print_shortest(io, x, dot, SHORTEST, 0)
print_shortest(io::IO, x::Union{AbstractFloat,Integer}) = print_shortest(io, float(x), false)

end # module

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
"""This module implements error losses, their gradients and equivalents."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

function error_loss(θ, z, regressionModel::RegressionModel)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    error("'error_loss' not defined for $(typedof(s))")
end

function error_loss(θ, z, linearRegression::LinearRegression)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    return (θ'*z[1:end-1] - z[end])^2
end

function error_loss(θ, z, logisticRegression::LogisticRegression)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------
    Name of constraint
    """
    return log(1+exp(-z[end]*(θ'*z[1:end-1])))
end

function error_prime(θ, z, k, regressionModel::RegressionModel)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    error("'error_loss' not defined for $(typedof(s))")
end

function error_prime(θ, z, k, regressionModel::LinearRegression)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    return 2*z[k]*(θ'*z[1:end-1] - z[end])
end

function error_prime(θ, z, k, regressionModel::LogisticRegression)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    return -(z[end]*z[k])/(1+exp(z[end]*(θ'*z[1:end-1])))
end

function equivalent(x, i, Z, regressionModel::RegressionModel)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    error("'error_loss' not defined for $(typedof(s))")
end

function equivalent(x, i, Z, regressionModel::LinearRegression)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    lt = size(Z)[2]-1

    λ = x[1]
    μ = x[2]
    θ = x[3:3+lt-1]
    s = x[3+lt:end]
    normZ = Z[i,1:end-1]'*Z[i,1:end-1]
    gthemu = (error_loss(θ,Z[i,:],regressionModel)-μ)^2
    step = zeros(length(x))
    step[1] = -λ^2(error_loss(θ,Z[i,:],regressionModel)-μ)/(4*normZ*error_loss(θ,Z[i,:],regressionModel)+gthemu-1)
    step[2] = -λ^2/(4*normZ*error_loss(θ,Z[i,:],regressionModel)+gthemu-1)
    for k in 1:lt-1
        step[3-1+k] = 2*λ^2*Z[i,k]*(θ'*Z[i,1:end-1] - Z[i,end])/(4*normZ*error_loss(θ,Z[i,:],regressionModel)+gthemu-1)
    end
    step[3+lt-1+i] = -λ^3/(4*normZ*error_loss(θ,Z[i,:],regressionModel)+gthemu-1)

    return step
end

function equivalent(x, i, Z, regressionModel::LogisticRegression)
    """
    Aggregator of projection

    Parameters
    ----------
    i::Int

    robustModel::RobustModel

    Returns
    -------

    """
    lt = size(Z)[2]-1

    λ = x[1]
    μ = x[2]
    θ = x[3:3+lt-1]
    s = x[3+lt:end]
    cik = (Z[i,end]^2)/(1+exp(-Z[i,end]*θ'Z[i,1:end-1]))^2*Z[i,1:end-1]'*Z[i,1:end-1]
    gthemu = (error_loss(θ,Z[i,:],regressionModel)-μ)^2
    step = zeros(length(x))
    step[1] = -λ^2(error_loss(θ,Z[i,:],regressionModel)-μ)/(cik+gthemu-1)
    step[2] = -λ^2/(cik+gthemu-1)
    for k in 1:lt-1
        step[3-1+k] = -λ^2/(cik+gthemu-1)*(Z[i,end]*Z[i,k])/(exp(error_loss(θ,Z[i,:],regressionModel)))
    end
    step[3+lt-1+i] = -λ^3/(cik+gthemu-1)

    return step
end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
# This file is part of RLinearAlgebra.jl

"""
    LSStopMaxIterations <: LinSysStopCriterion

A structure that specifies a maximum iteration stopping criterion. That is, once a method
    achieves a certain number of iterations, it stops.

# Fields
- `max_iter::Int64`, the maximum number of iterations.
"""
struct LSStopMaxIterations <: LinSysStopCriterion
    max_iter::Int64
end

# Common interface for stopping criteria
function check_stop_criterion(
    log::LinSysSolverLog,
    stop::LSStopMaxIterations
)
    return log.iterations == stop.max_iter ? true : false
end