text:
const CViewer = Ptr{Cvoid}
struct PetscViewer
    ptr::Ref{CViewer}
    comm::MPI.Comm

    PetscViewer(comm::MPI.Comm) = new(Ref{CViewer}(), comm)
end

# allows us to pass PetscViewer objects directly into CViewer ccall signatures
Base.cconvert(::Type{CViewer}, viewer::PetscViewer) = viewer.ptr[]

"""
    PetscViewerASCIIGetStdout(comm::MPI.Comm = MPI.COMM_WORLD)

Wrapper for `PetscViewerASCIIGetStdout`
"""
function PetscViewerASCIIGetStdout(comm::MPI.Comm = MPI.COMM_WORLD)
    viewer = PetscViewer(comm)
    error = ccall((:PetscViewerASCIIGetStdout, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CViewer}), comm, viewer.ptr)
    @assert iszero(error)
    return viewer
end

const PetscViewerStdWorld = PetscViewerASCIIGetStdout

"""
    PetscViewerCreate(comm::MPI.Comm = MPI.COMM_WORLD)

Wrapper for `PetscViewerCreate`
"""
function PetscViewerCreate(comm::MPI.Comm = MPI.COMM_WORLD)
    viewer = PetscViewer(comm)
    error = ccall((:PetscViewerCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CViewer}), comm, viewer.ptr)
    @assert iszero(error)
    return viewer
end

"""
    PetscViewerPushFormat(viewer::PetscViewer, format::PetscViewerFormat)

Wrapper for `PetscViewerPushFormat`
"""
function PetscViewerPushFormat(viewer::PetscViewer, format::PetscViewerFormat)
    error = ccall((:PetscViewerPushFormat, libpetsc), PetscErrorCode, (CViewer, PetscViewerFormat), viewer, format)
    @assert iszero(error)
end

"""
    PetscViewerPopFormat(viewer::PetscViewer)

Wrapper for `PetscViewerPopFormat`
"""
function PetscViewerPopFormat(viewer::PetscViewer)
    error = ccall((:PetscViewerPopFormat, libpetsc), PetscErrorCode, (CViewer,), viewer)
    @assert iszero(error)
end

"""
    PetscViewerASCIIOpen(comm::MPI.Comm, filename)

Wrapper for `PetscViewerASCIIOpen`
"""
function PetscViewerASCIIOpen(comm::MPI.Comm, filename)
    viewer = PetscViewer(comm)
    error = ccall((:PetscViewerASCIIOpen, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Cstring, Ptr{CViewer}), comm, filename, viewer.ptr)
    @assert iszero(error)
    return viewer
end

"""
PetscViewerFileSetName(viewer::PetscViewer, filename)

Wrapper for `PetscViewerFileSetName`
"""
function PetscViewerFileSetName(viewer::PetscViewer, filename::String)
    error = ccall((:PetscViewerFileSetName, libpetsc), PetscErrorCode, (CViewer, Cstring), viewer, filename)
    @assert iszero(error)
end

"""
PetscViewerFileSetMode(viewer::PetscViewer, mode::PetscFileMode = FILE_MODE_WRITE)

Wrapper for `PetscViewerFileSetMode`
"""
function PetscViewerFileSetMode(viewer::PetscViewer, mode::PetscFileMode = FILE_MODE_WRITE)
    error = ccall((:PetscViewerFileSetMode, libpetsc), PetscErrorCode, (CViewer, PetscFileMode), viewer, mode)
    @assert iszero(error)
end

"""
    PetscViewerHDF5Open(comm::MPI.Comm, filename::String, type::PetscFileMode)

Wrapper for `PetscViewerHDF5Open`
"""
function PetscViewerHDF5Open(comm::MPI.Comm, filename::String, type::PetscFileMode)
    viewer = PetscViewer(comm)
    error = ccall((:PetscViewerHDF5Open, libpetsc), PetscErrorCode,
        (MPI.MPI_Comm, Cstring, PetscFileMode, Ptr{CViewer}),
        comm, filename, type, viewer.ptr)
    @assert iszero(error)
    return viewer
end

"""
    PetscViewerSetType(viewer::PetscViewer, type::String)

Wrapper for `PetscViewerSetType`. Values for `type` alors available here:
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Viewer/PetscViewerType.html#PetscViewerType
"""
function PetscViewerSetType(viewer::PetscViewer, type::String)
    error = ccall((:PetscViewerSetType, libpetsc), PetscErrorCode, (CViewer, Cstring), viewer, type)
    @assert iszero(error)
end

"""
    PetscViewerView(v::PetscViewer, viewer::PetscViewer = PetscViewerStdWorld())

Wrapper to `PetscViewerView`
"""
function PetscViewerView(v::PetscViewer, viewer::PetscViewer = PetscViewerStdWorld())
    error = ccall((:PetscViewerView, libpetsc), PetscErrorCode, (CViewer, CViewer), v, viewer)
    @assert iszero(error)
end

"""
    PetscViewerDestroy(viewer::PetscViewer)

Wrapper for `PetscViewerDestroy`

Warning : from what I understand, all viewers must not be destroyed explicitely using `PetscViewerDestroy`.

"""
function PetscViewerDestroy(viewer::PetscViewer)
    error = ccall((:PetscViewerDestroy, libpetsc), PetscErrorCode, (Ptr{CViewer}, ), viewer.ptr)
    @assert iszero(error)
end

max_stars_repo_name: bmxam/JuliaPetsc

max_issues_repo_name: bmxam/JuliaPetsc

max_forks_repo_name: bmxam/JuliaPetsc


text:
module PlanningVsCompensatoryCaching
using SpecialFunctions, PGFPlotsX, LinearAlgebra, Distributions, DataStructures, YAML, HCubature, StaticArrays
push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usetikzlibrary{patterns}")

include("plotting.jl")
include("utils.jl")

end # module

ext: jl


text:
@testset "equal_to_prune_two" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:4] <= 5, Int)
    @constraint(m, sum(x) == 10)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]
    @test constraint.is_equal
    @test CS.count_unfixed(com, constraint) == 4
    x_var = com.search_space[constraint.indices]
    # fix two and test that it is still feasible
    @test CS.fix!(com, x_var[1], 3)
    @test CS.fix!(com, x_var[2], -3)

    # here 3 and 4 must be equal to 5 and we can prune with
    # prune_is_equal_two_var!
    @test CS.prune_is_equal_two_var!(com, constraint, constraint.fct)
    @test CS.values(x_var[3]) == [5]
    @test CS.values(x_var[4]) == [5]
end

@testset "equal_to_prune_two in all different" begin
    # test if 3 and 4 are also in all different constraint
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:4] <= 5, Int)
    @constraint(m, sum(x) == 10)
    @constraint(m, x[3:4] in CS.AllDifferentSet())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]
    x_var = com.search_space[constraint.indices]
    # fix two and test that it is still feasible
    @test CS.fix!(com, x_var[1], 3)
    @test CS.fix!(com, x_var[2], -3)
    @test CS.count_unfixed(com, constraint) == 2
    # not feasible as 3 and 4 would need to be 5 but there is an all different constraint
    @test !CS.prune_is_equal_two_var!(com, constraint, constraint.fct)
end

@testset "equal_to_prune_two not in same all different" begin
    # test if 3 and 4 are also in all different constraint
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:4] <= 5, Int)
    @constraint(m, sum(x) == 10)
    @constraint(m, x[1:3] in CS.AllDifferentSet())
    @constraint(m, [x[2], x[4]] in CS.AllDifferentSet())
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]
    x_var = com.search_space[constraint.indices]
    # fix two and test that it is still feasible
    @test CS.fix!(com, x_var[1], 3)
    @test CS.fix!(com, x_var[2], -3)
    @test CS.count_unfixed(com, constraint) == 2
    # feasible as 3 and 4 can be 5
    @test CS.prune_is_equal_two_var!(com, constraint, constraint.fct)
    @test CS.values(x_var[3]) == [5]
    @test CS.values(x_var[4]) == [5]
end

@testset "equal_to" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, y[1:3], CS.Integers([-3, 1, 2, 3]))
    @constraint(m, sum(y) + 1 == 5)
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]

    # doesn't check the length
    # 1+2+1 + constant (1) == 5
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 1])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 2])

    constr_indices = constraint.indices
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        -3,
    )
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        1,
    )
    # not actually feasible but will not be tested fully here
    CS.fix!(com, com.search_space[constr_indices[1]], 3)
    # will be tested when setting the next
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        1,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1
    # now setting it to 1 should be feasible
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        1,
    )

    # feasible but remove -3 and 3
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == [1, 2]
    end
    @test CS.fix!(com, com.search_space[constr_indices[1]], 2)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[2:3]
        @test CS.values(com.search_space[ind]) == [1]
    end
    CS.values(com.search_space[constr_indices[1]]) == [2]
end

@testset "eqsum is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, sum(x) == 10)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]

    variables = com.search_space
    @test CS.fix!(com, variables[1], 1; check_feasibility = false)
    @test CS.fix!(com, variables[2], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, sum(x) == 10)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]

    variables = com.search_space
    @test CS.fix!(com, variables[1], 5; check_feasibility = false)
    @test CS.fix!(com, variables[2], 5; check_feasibility = false)
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end

@testset "constraint without variables" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, x[2] - x[1] + (-x[2]) + x[1] == 10)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, x[2] - x[1] + (-x[2]) + x[1] == 0)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
end


text:
@testset "Meataxe" begin

  @testset "cleanvect" begin

    F, a = Nemo.FiniteField(3, 1, "a")
    M=MatrixSpace(F,2,3)([1,1,0,0,1,0])
    v=MatrixSpace(F,1,3)([2,2,0])
    @test iszero(Hecke.cleanvect(M,v))
    v[1,3]=1
    @test !iszero(Hecke.cleanvect(M,v))

  end

  @testset "closure and spinning" begin

   F, a = Nemo.FiniteField(3, 1, "a")
   G=[MatrixSpace(F,4,4)([1,2,0,0,1,1,0,0,0,0,1,2,0,0,1,1])]
   M=FqGModule(G)
   v=MatrixSpace(F,1,4)([1,0,0,0])
   X=Hecke.closure(v,M.G)
   Y=Hecke.spinning(v,M.G)
   @test rows(X)==2
   @test rows(Y)==2
   rref!(Y)
   @test X==Y
   w=MatrixSpace(F,1,4)([0,0,0,1])
   v=vcat(v,w)
   @test rows(Hecke.closure(v,M.G))==4

  end

  @testset "meataxe" begin

    F, a = Nemo.FiniteField(3, 1, "a")
    G=[MatrixSpace(F,4,4)([1,2,0,0,1,1,0,0,0,0,1,2,0,0,1,1])]
    M=FqGModule(G)
    bool,B=meataxe(M)
    @test !bool
    @test rows(B)==2
    @test rows(Hecke.closure(B, M.G))==2

    N=Hecke.actsub(B,G)
    bool,B=meataxe(N)
    @test bool

    G=[MatrixSpace(F,3,3)([1,0,0,0,0,1,0,1,0]), MatrixSpace(F,3,3)([0,0,1,1,0,0,0,1,0])]
    M=FqGModule(G)
    bool,B=meataxe(M)
    @test !bool
    @test rows(Hecke.closure(B, M.G))==rows(B)

  end

  @testset "composition factors and series" begin

    F, a = Nemo.FiniteField(3, 1, "a")
    G=[MatrixSpace(F,4,4)([1,2,0,0,1,1,0,0,0,0,1,2,0,0,1,1])]
    M=FqGModule(G)
    lf=composition_factors(M)
    cs=composition_series(M)
    @test length(lf)==1
    @test length(cs)==2
    x=Hecke.actsub(cs[1],M.G)
    @test Hecke.isisomorphic(lf[1][1],x)
    x,_=Hecke.actquo(cs[1],M.G)
    @test Hecke.isisomorphic(lf[1][1],x)

    M=FqGModule([matrix(F,2,2,[0,1,2,0])])
    N=FqGModule([matrix(F,2,2,[0,2,1,0])])
    M.isirreducible=true
    @test Hecke.isisomorphic(M,N)

    M1=matrix(F,2,2,[1,0,1,1])
    M2=matrix(F,2,2,[1,1,0,1])
    M=FqGModule([M1,M2])
    M.isirreducible=true

    N1=matrix(F,2,2,[2,2,1,0])
    N2=matrix(F,2,2,[1,1,0,1])
    N=FqGModule([N1,N2])

    @test Hecke.isisomorphic(M,N)

  end

  @testset "Submodules" begin

    F, a = Nemo.FiniteField(3, 1, "a")
    A=MatrixSpace(F,3,3)(1)
    M=FqGModule([A])
    ls=minimal_submodules(M)
    @test length(ls)==13

    F, a = Nemo.FiniteField(2, 1, "a")
    A=MatrixSpace(F,6,6)(1)
    A[5,6]=1
    M=FqGModule([A])
    ls=minimal_submodules(M)
    @test length(ls)==31
    ls=submodules(M,4)
    @test length(ls)==171

  end

end


text:
"""
    WB.jl

This is a Julia port of a MATLAB implementation of batch and incremental
WB-Index (WB) Cluster Validity Index.

Authors:
MATLAB implementation: Leonardo Enzo Brito da Silva
Julia port: Sasha Petrenko <sap625@mst.edu>

REFERENCES
[1] L. E. Brito da Silva, N. M. Melton, and D. C. Wunsch II, "Incremental
Cluster Validity Indices for Hard Partitions: Extensions  and  Comparative
Study," ArXiv  e-prints, Feb 2019, arXiv:1902.06711v1 [cs.LG].
[2] Q. Zhao, M. Xu, and P. Franti, "Sum-of-Squares Based Cluster Validity
Index and Significance Analysis," in Adaptive and Natural Computing Algorithms,
M. Kolehmainen, P. Toivanen, and B. Beliczynski, Eds. Berlin, Heidelberg:
Springer Berlin Heidelberg, 2009, pp. 313–322.
[3] Q. Zhao and P. Franti, "WB-index: A sum-of-squares based index for
cluster validity," Data Knowledge Engineering, vol. 92, pp. 77–89, 2014.
[4] M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, and J. Bailey,
"Online Cluster Validity Indices for Streaming Data," ArXiv e-prints, 2018,
arXiv:1801.02937v1 [stat.ML].
[5] M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, J. Bailey, "Online
cluster validity indices for performance monitoring of streaming data clustering,"
Int. J. Intell. Syst., pp. 1–23, 2018.
"""

"""
    WB

The stateful information of the WB-Index (WB) Cluster Validity Index.

# References
1. L. E. Brito da Silva, N. M. Melton, and D. C. Wunsch II, "Incremental Cluster Validity Indices for Hard Partitions: Extensions  and  Comparative Study," ArXiv  e-prints, Feb 2019, arXiv:1902.06711v1 [cs.LG].
2. Q. Zhao, M. Xu, and P. Franti, "Sum-of-Squares Based Cluster Validity Index and Significance Analysis," in Adaptive and Natural Computing Algorithms, M. Kolehmainen, P. Toivanen, and B. Beliczynski, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009, pp. 313–322.
3. Q. Zhao and P. Franti, "WB-index: A sum-of-squares based index for cluster validity," Data Knowledge Engineering, vol. 92, pp. 77–89, 2014.
4. M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, and J. Bailey, "Online Cluster Validity Indices for Streaming Data," ArXiv e-prints, 2018, arXiv:1801.02937v1 [stat.ML].
5. M. Moshtaghi, J. C. Bezdek, S. M. Erfani, C. Leckie, J. Bailey, "Online cluster validity indices for performance monitoring of streaming data clustering," Int. J. Intell. Syst., pp. 1–23, 2018.
"""
mutable struct WB <: CVI
    label_map::LabelMap
    dim::Integer
    n_samples::Integer
    mu::RealVector          # dim
    n::IntegerVector        # dim
    v::RealMatrix           # dim x n_clusters
    CP::RealVector          # dim
    SEP::RealVector         # dim
    G::RealMatrix           # dim x n_clusters
    BGSS::Float
    WGSS::Float
    n_clusters::Integer
    criterion_value::Float
end # WB <: CVI

"""
    WB()

Default constructor for the WB-Index (WB) Cluster Validity Index.
"""
function WB()
    WB(
        LabelMap(),                     # label_map
        0,                              # dim
        0,                              # n_samples
        Array{Float, 1}(undef, 0),      # mu
        Array{Integer, 1}(undef, 0),    # n
        Array{Float, 2}(undef, 0, 0),   # v
        Array{Float, 1}(undef, 0),      # CP
        Array{Float, 1}(undef, 0),      # SEP
        Array{Float, 2}(undef, 0, 0),   # G
        0.0,                            # BGSS
        0.0,                            # WGSS
        0,                              # n_clusters
        0.0                             # criterion_value
    )
end # WB()

"""
    setup!(cvi::WB, sample::Vector{T}) where {T<:Real}
"""
function setup!(cvi::WB, sample::Vector{T}) where {T<:Real}
    # Get the feature dimension
    cvi.dim = length(sample)
    # Initialize the augmenting 2-D arrays with the correct feature dimension
    # NOTE: R is emptied and calculated in evaluate!, so it is not defined here
    cvi.v = Array{T, 2}(undef, cvi.dim, 0)
    cvi.G = Array{T, 2}(undef, cvi.dim, 0)
end # setup!(cvi::WB, sample::Vector{T}) where {T<:Real}

function param_inc!(cvi::WB, sample::RealVector, label::Integer)
    # Get the internal label
    i_label = get_internal_label!(cvi.label_map, label)

    n_samples_new = cvi.n_samples + 1
    if isempty(cvi.mu)
        mu_new = sample
        setup!(cvi, sample)
    else
        mu_new = (1 - 1/n_samples_new) .* cvi.mu + (1/n_samples_new) .* sample
    end

    if i_label > cvi.n_clusters
        n_new = 1
        v_new = sample
        CP_new = 0.0
        G_new = zeros(cvi.dim)
        # Update 1-D parameters with a push
        cvi.n_clusters += 1
        push!(cvi.CP, CP_new)
        push!(cvi.n, n_new)
        # Update 2-D parameters with appending and reassignment
        cvi.v = [cvi.v v_new]
        cvi.G = [cvi.G G_new]
    else
        n_new = cvi.n[i_label] + 1
        v_new = (1 - 1/n_new) .* cvi.v[:, i_label] + (1/n_new) .* sample
        delta_v = cvi.v[:, i_label] - v_new
        diff_x_v = sample .- v_new
        CP_new = cvi.CP[i_label] + transpose(diff_x_v)*diff_x_v + cvi.n[i_label]*transpose(delta_v)*delta_v + 2*transpose(delta_v)*cvi.G[:, i_label]
        G_new = cvi.G[:, i_label] + diff_x_v + cvi.n[i_label].*delta_v
        # Update parameters
        cvi.n[i_label] = n_new
        cvi.v[:, i_label] = v_new
        cvi.CP[i_label] = CP_new
        cvi.G[:, i_label] = G_new
    end
    cvi.n_samples = n_samples_new
    cvi.mu = mu_new
    cvi.SEP = [cvi.n[ix] * sum((cvi.v[:, ix] - cvi.mu).^2) for ix=1:cvi.n_clusters]
end # param_inc!(cvi::WB, sample::RealVector, label::Integer)

function param_batch!(cvi::WB, data::RealMatrix, labels::IntegerVector)
    cvi.dim, cvi.n_samples = size(data)
    # Take the average across all samples, but cast to 1-D vector
    cvi.mu = mean(data, dims=2)[:]
    # u = findfirst.(isequal.(unique(labels)), [labels])
    u = unique(labels)
    cvi.n_clusters = length(u)
    cvi.n = zeros(Integer, cvi.n_clusters)
    cvi.v = zeros(cvi.dim, cvi.n_clusters)
    cvi.CP = zeros(cvi.n_clusters)
    cvi.SEP = zeros(cvi.n_clusters)
    for ix = 1:cvi.n_clusters
        subset = data[:, findall(x->x==u[ix], labels)]
        cvi.n[ix] = size(subset, 2)
        cvi.v[1:cvi.dim, ix] = mean(subset, dims=2)
        diff_x_v = subset - cvi.v[:, ix] * ones(1, cvi.n[ix])
        cvi.CP[ix] = sum(diff_x_v.^2)
        cvi.SEP[ix] = cvi.n[ix] * sum((cvi.v[:, ix] - cvi.mu).^2);
    end
end # param_batch!(cvi::WB, data::RealMatrix, labels::IntegerVector)

function evaluate!(cvi::WB)
    # Within group sum of scatters
    cvi.WGSS = sum(cvi.CP)
    if cvi.n_clusters > 1
        # Between groups sum of scatters
        cvi.BGSS = sum(cvi.SEP)
        # WB index value
        cvi.criterion_value = (cvi.WGSS / cvi.BGSS) * cvi.n_clusters
    else
        cvi.BGSS = 0.0
        cvi.criterion_value = 0.0;
    end
end # evaluate(cvi::WB)


text:
module AdaptiveFilters
using OnlineStats, DSP

export adaptive_filter, focused_adaptive_filter, OMAP, MSPI, OMAS, ADAM, ExponentialWeight, EqualWeight

"""
    yh = adaptive_filter(y, alg=MSPI; order=4, lr=0.1)

Filters `y` with an adaptive AR (only poles) filter with specified order.
Returns `yh` which is the predicted output from an adaptive line enhancer (ALE). If your noise is wideband and signal narrowband, `yh` is your desired filtered signal. If the noise is narrowband and the signal is wideband, then `y-yh` is your desired filtered signal.

The first `order` samples of `yh` will be copies of `y`. The signals will thus have the same length.

#Arguments:
- `alg`: Stochastic approximation algorithm or weight function. Examples: `OMAP, MSPI, OMAS, ADAM, ExponentialWeight, EqualWeight`
- `y`: Input signal
- `order`: Filter order
- `lr`: Learning rate or weight depending on `alg`
"""
function adaptive_filter(y, alg::Type{<:OnlineStats.Algorithm}=MSPI; order=4, lr=0.1)
    T = length(y)
    model = StatLearn(order, alg(), rate=LearningRate(lr))
    yh = similar(y)
    yh[1:order] = y[1:order]
    for t in 1:T-order
        x = @view(y[t:t+order-1])
        o  = fit!(model, (x, y[t+order]))
        yh[t+order] = OnlineStats.predict(o, x)
    end
    yh
end

function adaptive_filter(y, alg::Type{<:OnlineStats.Weight}; order=6, lr=0.2)
    T = length(y)
    model = LinReg(weight=alg(lr))
    yh = similar(y)
    yh[1:order] = y[1:order]
    for t in 1:T-order
        x = @view(y[t:t+order-1])
        o  = fit!(model, (x, y[t+order]))
        if t > order
            value(model)
            yh[t+order] = OnlineStats.predict(model, x)
        else
            yh[t+order] = y[t]
        end
    end
    yh
end

"""
    focused_adaptive_filter(y, band, fs, args...; kwargs...)

An adaptive filter that focuses its attention to a specific frequency range.

#Arguments:
- `y`: Input signal
- `band`: Frequency band typle
- `fs`: Sample rate, e.g., 44100
- `args`: Passed through to `adaptive_filter`
- `kwargs`: Passed through to `adaptive_filter`
"""
function focused_adaptive_filter(y,band,fs,args...; kwargs...)
    dm = Butterworth(3)
    rt = Bandpass(band...,fs=fs)
    f = digitalfilter(rt, dm)
    yh = filt(f, reverse(y))
    yh = filt(f, reverse(yh))
    yh = adaptive_filter(yh, args...; kwargs...)
    yh
end

end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    struct NetworkError
        url::String
    end

Simple structure that is thrown internally if there was an error retrieving data from the network.

__Values:__

| Name | Description
| :--- | :----------
| url  | The URL that caused the network error.
"""
struct NetworkError
    url::String
end

# =================================================================================================
"""
    struct FileError
        file::String
    end

Simple structure that is thrown internally if there was an error accessing a file.

__Values:__

| Name | Description
| :--- | :----------
| file | The full path of the file that caused the error.
"""
struct FileError
    file::String
end


text:
using Plots
gr()

v = rand(10)
w = rand(10)

print(v + w)

using LinearAlgebra
# soma dos vetores v e w
print(dot(v, w))


text:
# HeaderSearch
function PrintStats(x::HeaderSearch)
    @check_ptrs x
    return clang_HeaderSearch_PrintStats(x)
end


text:
"""
Package: Forecast

    pacf(x::{AbstractVector,DataFrame};
         type = "step-real",
         lag = Integer(ceil(10*log10(length(x)))),
         alpha = (0.95,0.99))

Compute the partial auto-correlation for an univariate series.

There are two versions; the "step" version estimates the auto-regressive parameters of an increasing model, the "real" version estimates the actual partial auto-correlation by eliminating the linear information provided by the lags. When using the default type "stepwise-real" both versions are calculated.

The distribution used to estimate the confidence intervals is an approximation of a Fisher Transformation via a Normal Distribution. There is a plot recipe for the returned object.

# Arguments
- `x`: Vector or uni-dimensional DataFrame of data.
- `type` = Valid values are "stepwise", "real" and "stepwise-real".
- `lag`: Maximum number of lags.
- `alpha`: A tuple with two thresholds (t1,t2) with t1 <= t2 to plot confidence intervals. The default values are 0.95 and 0.99.

# Returns
A CCF object

# Examples
```julia-repl
julia> x = rand(100);
res = pacf(x);
plot(res)
```
"""
function pacf(df::DataFrame;
              type::String = "stepwise-real",
              lag::Integer = Integer(ceil(10*log10(length(ta)))),
              alpha::Tuple{Real,Real} = (0.95,0.99))

    x = Array(df[:,eltype.(eachcol(df)) .<: Real])

    pacf(x; type, lag, alpha)

end

function pacf(x::AbstractVector{<:Real};
              type::String = "stepwise-real",
              lag::Integer = Integer(ceil(10*log10(length(x)))),
              alpha::Tuple{Real,Real} = (0.95,0.99))

    N = length(x)
    @assert type in  ["stepwise","real","stepwise-real"] "The options for `type` are `stepwise`, `real` and `both`"
    @assert 1 <= lag <= N-4
    @assert length(alpha) == 2
    @assert 0.0 < alpha[1] < alpha[2] < 1.0

    # Confidence Intervals
    function fse(v::Real)::AbstractFloat
        1/sqrt(v-3)
    end

    a1 = alpha[1]
    a2 = alpha[2]
    z1 = quantile(Normal(), a1 + (1-a1)/2)
    ci1 = z1*fse(N-lag)
    z2 = quantile(Normal(), a2 + (1-a2)/2)
    ci2 = z2*fse(N-lag)
    ci = (ci1,ci2)

    call = "pacf(x"*
        "; type=\""*type*
        "\", lag="*string(lag)*
        ", alpha="*string(alpha)*")"

    if (type::String == "stepwise")
        pac = pacf_stepwise(x; lag, alpha)
        return CCF(pac, N, "pacf_stepwise", lag, alpha, ci, true, call)
    end

    if (type::String == "real")
        pac = pacf_real(x; ag, alpha)
        return CCF(pac, N, "pacf_real", lag, alpha, ci, true, call)
    end

    if type::String == "stepwise-real"
        pacs = pacf_stepwise(x; lag, alpha)
        pacr = pacf_real(x; lag, alpha)
        return CCF(hcat(pacs, pacr),
                   N, "pacf_stepwise-real", lag, alpha, ci, true, call)
    end

end

function pacf_lag(x::AbstractVector{<:Real};
                   lag::Integer = Integer(ceil(10*log10(length(x)))),
                   alpha::Tuple{Real,Real} = (0.95,0.99))

    R = vcat(cov(x,x),acf(x; lag, type="cov").ccf)
    p = length(R)-1

    A = zeros(p,p)
    for i in p:-1:1
        A[i,i:end] = R[1:p-i+1]
        A[i:end,i] = R[1:p-i+1]
    end
    b = R[2:end]

    lambda = 10^-36
    (A'*A+lambda*I(p))\(A'*b)

end

function pacf_stepwise(x::AbstractVector{<:Real};
                   lag::Integer = Integer(ceil(10*log10(length(x)))),
                   alpha::Tuple{Real,Real} = (0.95,0.99))

    #TODO optimize with Levinson-Durbin algorithm
    pac = []
    for i in 1:lag
        push!(pac, pacf_lag(x; lag = i)[end])
    end
    pac

end

function pacf_real(x::AbstractVector{<:Real};
             lag::Integer = Integer(ceil(10*log10(length(x)))),
             alpha::Tuple{Real,Real} = (0.95,0.99))

    # prevent singularities
    lambda = 10^-36

    N = length(x)
    A = zeros(N-lag,lag+1)
    for i in 1:(N-lag)
        A[i,:] = x[i:i+lag]
    end
    bt1 = A[:,1]
    A = A[:,2:end]
    A = hcat(A,repeat([1.0],N-lag))

    pacf_res = []
    for i in 2:lag+1
        A1i = drop(A, c=i-1)
        lsq = (A1i'*A1i+lambda*I(lag))\(A1i'*bt1)
        r1 = round.(bt1.-A1i*lsq, digits=10)

        bt0 = A[:,i-1]
        lsq = (A1i'*A1i+lambda*I(lag))\(A1i'*bt0)
        r0 = round.(bt0.-A1i*lsq, digits=10)
        push!(pacf_res,cor(r0,r1))

    end

    pacf_res
end


text:
using Knet
libknet8handle = Libdl.dlopen(Knet.libknet8)

SIZE = 100000
ITER = 100000
x32 = KnetArray(rand(Float32,SIZE))
y32 = KnetArray(rand(Float32,SIZE))
z32 = similar(x32)
x64 = KnetArray(rand(Float64,SIZE))
y64 = KnetArray(rand(Float64,SIZE))
z64 = similar(x64)

function cuda11test(fname, jname=fname, o...)
    println(fname)
    fcpu = eval(parse(jname))
    f32 = Libdl.dlsym(libknet8handle, fname*"_32_11")
    @time cuda11rep(f32,x32,y32,z32)
    isapprox(Array(z32),fcpu(Array(x32),Array(y32))) || warn("$fname 32")
    f64 = Libdl.dlsym(libknet8handle, fname*"_64_11")
    @time cuda11rep(f64,x64,y64,z64)
    isapprox(Array(z64),fcpu(Array(x64),Array(y64))) || warn("$fname 64")
end

function cuda11rep{T}(f,x::KnetArray{T},y::KnetArray{T},z::KnetArray{T})
    n = Cint(length(z))
    for i=1:ITER
        ccall(f,Void,(Cint,Ptr{T},Ptr{T},Ptr{T}),n,x,y,z)
    end
    Knet.@cuda(cudart,cudaDeviceSynchronize,())
    Knet.@cuda(cudart,cudaGetLastError,())
end

for f in Knet.cuda11
    isa(f,Tuple) || (f=(f,))
    cuda11test(f...)
end


text:
# Problem 3

using CairoMakie

const G = 6.67 * 1e-8
const M = 1.989 * 1e33
const R = 8.3 * 1e7
const μe = 2.0

function Tc(ρc, K, γ, C)
	(μe/R) * (C * G * M^(2/3) * ρc^(1/3) - (K * ρc^(γ - 1)/μe^γ))
end

ρc = 10 .^LinRange(6, 9, 50)
C = 10
ks = [1.24 * 1e15, 1e13]
γs = [4/3, 5/3]
colors = [:orange, :cyan]
fig = Figure(resolution=(600, 400), font = "Times New Roman", fontsize = 12)
ax = fig[1,1] = Axis(fig, xlabel = "rho_c", ylabel = "T_c")
lins = [lines!(ρc, Tc.(ρc, kval, γval, C), color = color, linewidth = 2) for (kval, γval, color) in zip(ks, γs, colors)]
legs = ["gamma = $γval, K = $kval" for (kval, γval) in zip(ks, γs)]
leg = fig[1,1] = Legend(fig, lins, legs)
fig


text:
__precompile__()
module ConstantQ

using Compat

export
    GeometricFrequency,   # geometrically spaced frequency
    KernelProperty,       # kernel property
    SpectralKernelMatrix, # frequency-domain kernel matrix
    TemporalKernelMatrix, # time-domain kernel matrix
    property,
    nbins_per_octave,     # number of frequency bins per octave
    nfreqs,               # number of frequency bins
    freqs,                # generate array of frequencies
    q,                    # Q-factor
    speckernel,           # construct frequency-domain kernel
    tempkernel,           # construct time-domain kernel
    cqt                   # A fast constant-Q transform

include("cqt.jl")

end # module


text:
#=

MCMC utility functions for Tapestree

Ignacio Quintero Mächler

February 6 2017

t(-_-t)

=#

"""
    uniupt(p::Float64, tn::Float64)

Uniform parameter window move.
"""
uniupt(p::Float64, tn::Float64) = abs(p + (rand()-0.5) * tn)

"""
    addupt(p::Float64, tn::Float64)

Gaussian parameter window move.
"""
addupt(p::Float64, tn::Float64) = p + randn() * tn

"""
    addupt_lims(p::Float64, tn::Float64, xmin::Float64, xmax::Float64)

Gaussian parameter window move within a region of interest using rejection.
"""
function addupt_lims(p::Float64, tn::Float64, xmin::Float64, xmax::Float64)

  s = p + randn() * tn
  if !(xmin < s < xmax)
    s = p
  end

  return s
end

"""
    addupt(p::Float64, tn::Float64)

Gaussian parameter window move to vector.
"""
function addupt!(p::Vector{Float64}, tn::Vector{Float64}, j::Int64, i::Int64)

  @inbounds p[j] += randn() * tn[i]

  return nothing
end

"""
    duoupd(xd1::Float64,
           xd2::Float64,
           td1::Float64,
           td2::Float64,
           σ²ϕ::Float64)

Duo of Gaussians parameter update.
"""
function duoupd(xd1::Float64,
                xd2::Float64,
                td1::Float64,
                td2::Float64,
                σ²ϕ::Float64)
  invt = 1.0/(td1 + td2)
  return randn()*sqrt(td1 * td2 * invt * σ²ϕ) +
    (td2 * invt * xd1 + td1 * invt * xd2)
end

"""
    trioupd(xpr::Float64,
            xd1::Float64,
            xd2::Float64,
            tpr::Float64,
            td1::Float64,
            td2::Float64,
            σ²ϕ::Float64)

Trio of Gaussians parameter update.
"""
function trioupd(xpr::Float64,
                 xd1::Float64,
                 xd2::Float64,
                 tpr::Float64,
                 td1::Float64,
                 td2::Float64,
                 σ²ϕ::Float64)

    t = 1.0/(1.0/tpr + 1.0/td1 + 1.0/td2)
    return randn()*sqrt(t*σ²ϕ) + (xpr/tpr + xd1/td1 + xd2/td2)*t
end

"""
    absaddupt(p::Float64, tn::Float64)

Non-negative Gaussian parameter window move.
"""
absaddupt(p::Float64, tn::Float64) = abs(p + randn() * tn)

"""
    mulupt(p::Float64, tn::Float64)

Multiplicative parameter window move.
"""
mulupt(p::Float64, tn::Float64) = p * exp((rand() - 0.5) * tn)

"""
    makescalef(obj_ar::Float64)

Make scaling function given the objective acceptance rates.
"""
function makescalef(obj_ar::Float64)
  nar::Float64 = 1.0 - obj_ar

  function f(window::Float64, rate::Float64)
    if rate > obj_ar
      window *= (1.0 + (rate - obj_ar) / nar)::Float64
    else
      window /= (2.0 - rate / obj_ar)::Float64
    end

    return window::Float64
  end

  return f
end

"""
    globalscalef(λ::Float64, grate::Float64, stepsize::Float64, obj_ar::Float64)

Estimate global scaling factor.
"""
function globalscalef(λ       ::Float64,
                      grate   ::Float64,
                      stepsize::Float64,
                      obj_ar  ::Float64)
  return (exp(log(λ) + stepsize * (grate - obj_ar)))::Float64
end

"""
    adaptiveupd!(Σ::Array{Float64,2}, psam::Array{Float64,1}, pmean::Array{Float64,1}, stepsize::Float64)

Adaptive update for parameter mean and Σ in place.
"""
function adaptiveupd!(Σ       ::Array{Float64,2},
                      psam    ::Array{Float64,1},
                      pmean   ::Array{Float64,1},
                      stepsize::Float64)

  @inbounds begin
    for i in Base.OneTo(length(pmean))
      psam[i]  -= pmean[i]::Float64
      pmean[i] += (stepsize * psam[i])::Float64
    end

    BLAS.axpy!(stepsize,
               BLAS.gemm!('N', 'T', 1.0, psam, psam, -1.0, copy(Σ)),
               Σ)
  end
end

"""
    makestepsize(C::Float64, η::Float64)

Make function for the stepsize for the adaptive update.
"""
function makestepsize(η::Float64)

  β::Float64 = rand(range((1.0/(1.0 + η)),1))

  function f(t::Float64, C::Float64)
    return (C/(t^β))::Float64
  end

  return f
end

"""
    makemvnproposal(Σ::Array{Float64,2})

Make the multivariate update given the covariance matrix.
"""
function makemvnproposal(Σ::Array{Float64,2})

  spde = *(eigvecs(Σ),sqrt.(diagm(eigvals(Σ))))
  ln   = size(spde,1)

  function f(pvec::Array{Float64,1})
    (pvec .+
     BLAS.gemv('N', spde, randn(ln))::Array{Float64,1}
     )::Array{Float64,1}
  end

  return f
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 0b8acfe4-5073-11eb-1630-33963891fa65
begin
	using Revise

	using MAT

	using Clustering
	using ISODATA
	using MLDataUtils
	using MultivariateStats
	using Plots
	using StatsBase
end

# ╔═╡ 9d50e9dc-5112-11eb-261e-fb54e7238993
gradient = cgrad([:black, :purple, :red, :yellow])

# ╔═╡ 24a82cb0-5073-11eb-198e-a9de05e0c7a3
begin
	dfile = matopen("data/SalinasA.mat")
	Xraw = read(dfile)["salinasA"]
	close(dfile)
	Xraw
end

# ╔═╡ 103c6df2-5077-11eb-1ee3-137d3269b6eb
heatmap(Xraw[:, :, 4])

# ╔═╡ a688fe58-5073-11eb-2968-bd2d796d1f2f
begin
	tfile = matopen("data/SalinasA_gt.mat")
	Yraw = read(tfile)["salinasA_gt"]
	close(tfile)
	Yraw
end

# ╔═╡ b58883cc-5075-11eb-0c36-87c4d7f584f3
heatmap(Yraw, c=gradient)

# ╔═╡ 02200e30-5074-11eb-18c1-bbbd91537f47
X = reshape(deepcopy(Xraw), 7138, 224)'

# ╔═╡ 97b50c02-5074-11eb-1974-47e0028a67bb
begin
	Y = reshape(deepcopy(Yraw), 7138)
	for i in 10:14
		Y[findall(x -> x == Float64(i), Y)] .= i - 8.0
	end
	Y = round.(Int, Y)
end

# ╔═╡ efe37d6e-50f6-11eb-1e0d-c3b097dac6d6
heatmap(reshape(Y, 83, 86), c=gradient)

# ╔═╡ 59fa6a82-510d-11eb-06a9-c729dc0b2920
struct ExperimentResult
	iter::Float64
	centers::Float64
	adj_idx::Float64
	rand_idx::Float64
	mirkin_idx::Float64
	hubert_idx::Float64
end

# ╔═╡ 32c3f282-5110-11eb-050b-e3dc9548a51e
function avg_result(results::Vector{ExperimentResult})
	n = size(results, 1)

	iters = 0.0
	centers = 0.0
	adjs = 0.0
	rands = 0.0
	mirkins = 0.0
	huberts = 0.0
	for i in 1:n
		iters += results[i].iter
		centers += results[i].centers
		adjs += results[i].adj_idx
		rands += results[i].rand_idx
		mirkins += results[i].mirkin_idx
		huberts += results[i].hubert_idx
	end

	avg_iters = iters/n
	avg_centers = centers/n
	avg_adjs = adjs/n
	avg_rands = rands/n
	avg_mirkins = mirkins/n
	avg_huberts = huberts/n
	return ExperimentResult(
		avg_iters, avg_centers, avg_adjs, avg_rands, avg_mirkins, avg_huberts)
end

# ╔═╡ 1db73a8c-510d-11eb-089b-7b92b54526e5
# kmeans on salinasA
function experiment4(
	data,
	k,
	num_trials,
	iseeds;
	normalize=true
)
	results = Vector{ExperimentResult}(undef, num_trials)

	best_score = 0.0
	best_centers = Array{Float64}
	best_idx = 0

	for i in 1:num_trials
		train, test = splitobs(shuffleobs(data), at=0.8)
		train_x, train_y = train[1:end-1, :], round.(Int, train[end, :])
		test_x, test_y = test[1:end-1, :], round.(Int, test[end, :])

		if normalize
			train_z = fit(ZScoreTransform, train_x, dims=2)
			train_x = StatsBase.transform(train_z, train_x)
			test_x = StatsBase.transform(train_z, test_x)
		end

		# run kmeans
		res = ISODATA.kmeans(train_x, 6, iseeds=iseeds)

		# fit centers to test data
		fit_y = clusterfit(test_x, res.centers)

		# cluster analysis
		adj_idx, rand_idx, mirkin_idx, hubert_idx = randindex(test_y, fit_y)
		num_centers = size(res.centers, 2)

		results[i] = ExperimentResult(
			res.iter, num_centers, adj_idx, rand_idx, mirkin_idx, hubert_idx)

		if adj_idx > best_score
			best_score = adj_idx
			best_centers = res.centers
			best_idx = i
		end
	end

	return results[best_idx], best_centers, avg_result(results)
end

# ╔═╡ 415ae730-50f5-11eb-03e3-6365247ed76b
begin
	inds = findall(x -> x != 0.0, Y)

	M = fit(PCA, X; maxoutdim=3)
	Xout = transform(M, X)

	full = transpose(cat(Xout', Y, dims=2))
	filtered = full[:, inds]
end

# ╔═╡ f65644f8-5114-11eb-2fd9-09e75a642aca
transpose(cat(X', Y, dims=2))

# ╔═╡ f9a6065e-510d-11eb-2a36-e72f2d507a66
rkmeans_res, rkmeans_centers, rkmeans_avg = experiment4(
	filtered, 6, 100, randomcenters)

# ╔═╡ 226e3d78-510f-11eb-2d8d-155aabe1918d
begin
	rkmeans_pred = clusterfit(filtered[1:end-1, :], rkmeans_centers)
	rkmeans_Y = deepcopy(Y)
	rkmeans_Y[inds] .= rkmeans_pred
	heatmap(reshape(rkmeans_Y, 83, 86), c=gradient)
end

# ╔═╡ afa88d40-510f-11eb-327f-ab765d0a3c8b
kkmeans_res, kkmeans_centers, kkmeans_avg = experiment4(
	filtered, 6, 100, kmppcenters)

# ╔═╡ baae817c-510f-11eb-372b-05a0b920d188
begin
	kkmeans_pred = clusterfit(filtered[1:end-1, :], kkmeans_centers)
	kkmeans_Y = deepcopy(Y)
	kkmeans_Y[inds] .= kkmeans_pred
	heatmap(reshape(kkmeans_Y, 83, 86), c=gradient)
end

# ╔═╡ e8d00128-510f-11eb-359e-9964c959ce74
ikmeans_res, ikmeans_centers, ikmeans_avg = experiment4(
	filtered, 6, 100, isodata_s, normalize=true)

# ╔═╡ f17da848-510f-11eb-04b6-7d129dec43fb
begin
	ikmeans_pred = clusterfit(filtered[1:end-1, :], ikmeans_centers)
	ikmeans_Y = deepcopy(Y)
	ikmeans_Y[inds] .= ikmeans_pred
	heatmap(reshape(ikmeans_Y, 83, 86), c=gradient)
end

# ╔═╡ b159bdf4-5110-11eb-00a0-1bece66bc22b
# isodata on salinasA
function experiment5(
	data,
	k,
	num_trials,
	iseeds;
	normalize=true,
	_θn_bias::Float64=0.20,
	_θe_bias::Float64=0.20,
	_θc_bias::Float64=0.20
)
	results = Vector{ExperimentResult}(undef, num_trials)

	best_score = 0.0
	best_centers = Array{Float64}
	best_idx = 0

	for i in 1:num_trials
		train, test = splitobs(shuffleobs(data), at=0.8)
		train_x, train_y = train[1:end-1, :], round.(Int, train[end, :])
		test_x, test_y = test[1:end-1, :], round.(Int, test[end, :])

		if normalize
			train_z = fit(ZScoreTransform, train_x, dims=2)
			train_x = StatsBase.transform(train_z, train_x)
			test_x = StatsBase.transform(train_z, test_x)
		end

		# calculate isodata params
		summaries = summarize(train)
		θn, θe, θc = iso_param(
			train_x,
			summaries,
			_θn_bias=_θn_bias,
			_θe_bias=_θe_bias,
			_θc_bias=_θc_bias)

		# run kmeans
		res = base_conv(train_x, k, θn=θn, θe=θe, θc=θc, iseeds=iseeds)

		# fit centers to test data
		fit_y = clusterfit(test_x, res.centers)

		# cluster analysis
		adj_idx, rand_idx, mirkin_idx, hubert_idx = randindex(test_y, fit_y)
		num_centers = size(res.centers, 2)

		results[i] = ExperimentResult(
			res.iter, num_centers, adj_idx, rand_idx, mirkin_idx, hubert_idx)

		if adj_idx > best_score
			best_score = adj_idx
			best_centers = res.centers
			best_idx = i
		end
	end

	return results[best_idx], best_centers, avg_result(results)
end

# ╔═╡ 5b33ed70-5111-11eb-0795-015a8472911b
risodata_res, risodata_centers, risodata_avg = experiment5(
	filtered, 6, 1000, randomcenters, normalize=true,
	_θn_bias=0.0,
	_θe_bias=0.0,
	_θc_bias=0.0)

# ╔═╡ 02b383c0-5116-11eb-2748-d371c5aa42b6
begin
	risodata_pred = clusterfit(filtered[1:end-1, :], risodata_centers)
	risodata_Y = deepcopy(Y)
	risodata_Y[inds] .= risodata_pred
	heatmap(reshape(risodata_Y, 83, 86), c=gradient)
end

# ╔═╡ 9bc5aa2e-5111-11eb-06d7-b3cb072d95cd
kisodata_res, kisodata_centers, kisodata_avg = experiment5(
	filtered, 6, 1000, kmppcenters, normalize=true,
	_θn_bias=0.0,
	_θe_bias=0.0,
	_θc_bias=0.0)

# ╔═╡ aae23512-5115-11eb-0b01-4764155fb2d1
begin
	kisodata_pred = clusterfit(filtered[1:end-1, :], kisodata_centers)
	kisodata_Y = deepcopy(Y)
	kisodata_Y[inds] .= kisodata_pred
	heatmap(reshape(kisodata_Y, 83, 86), c=gradient)
end

# ╔═╡ a66f7ffc-5111-11eb-12a4-91874657e647
iisodata_res, iisodata_centers, iisodata_avg = experiment5(
	filtered, 6, 1000, isodata_s, normalize=true,
	_θn_bias=0.0,
	_θe_bias=0.0,
	_θc_bias=0.0)

# ╔═╡ d1d323d6-5113-11eb-3835-c13bc51b8e13
begin
	iisodata_pred = clusterfit(filtered[1:end-1, :], iisodata_centers)
	iisodata_Y = deepcopy(Y)
	iisodata_Y[inds] .= iisodata_pred
	heatmap(reshape(iisodata_Y, 83, 86), c=gradient)
end

# ╔═╡ Cell order:
# ╠═0b8acfe4-5073-11eb-1630-33963891fa65
# ╠═9d50e9dc-5112-11eb-261e-fb54e7238993
# ╠═24a82cb0-5073-11eb-198e-a9de05e0c7a3
# ╠═103c6df2-5077-11eb-1ee3-137d3269b6eb
# ╠═a688fe58-5073-11eb-2968-bd2d796d1f2f
# ╠═b58883cc-5075-11eb-0c36-87c4d7f584f3
# ╠═02200e30-5074-11eb-18c1-bbbd91537f47
# ╠═97b50c02-5074-11eb-1974-47e0028a67bb
# ╠═efe37d6e-50f6-11eb-1e0d-c3b097dac6d6
# ╠═59fa6a82-510d-11eb-06a9-c729dc0b2920
# ╠═32c3f282-5110-11eb-050b-e3dc9548a51e
# ╠═1db73a8c-510d-11eb-089b-7b92b54526e5
# ╠═415ae730-50f5-11eb-03e3-6365247ed76b
# ╠═f65644f8-5114-11eb-2fd9-09e75a642aca
# ╠═f9a6065e-510d-11eb-2a36-e72f2d507a66
# ╠═226e3d78-510f-11eb-2d8d-155aabe1918d
# ╠═afa88d40-510f-11eb-327f-ab765d0a3c8b
# ╠═baae817c-510f-11eb-372b-05a0b920d188
# ╠═e8d00128-510f-11eb-359e-9964c959ce74
# ╠═f17da848-510f-11eb-04b6-7d129dec43fb
# ╠═b159bdf4-5110-11eb-00a0-1bece66bc22b
# ╠═5b33ed70-5111-11eb-0795-015a8472911b
# ╠═02b383c0-5116-11eb-2748-d371c5aa42b6
# ╠═9bc5aa2e-5111-11eb-06d7-b3cb072d95cd
# ╠═aae23512-5115-11eb-0b01-4764155fb2d1
# ╠═a66f7ffc-5111-11eb-12a4-91874657e647
# ╠═d1d323d6-5113-11eb-3835-c13bc51b8e13


text:
using Documenter, ForwardDiff

makedocs(modules=[ForwardDiff],
         doctest = false,
         format = :html,
         sitename = "ForwardDiff",
         pages = ["Introduction" => "index.md",
                  "User Documentation" => [
                    "Limitations of ForwardDiff" => "user/limitations.md",
                    "Differentiation API" => "user/api.md",
                    "Advanced Usage Guide" => "user/advanced.md",
                    "Upgrading from Older Versions" => "user/upgrade.md"],
                  "Developer Documentation" => [
                    "How ForwardDiff Works" => "dev/how_it_works.md",
                    "How to Contribute" => "dev/contributing.md"]])

deploydocs(repo = "github.com/JuliaDiff/ForwardDiff.jl.git",
           osname = "linux",
           julia = "0.7",
           target = "build",
           deps = nothing,
           make = nothing)


text:
using Jeeves
using DataFrames:DataFrame
using Test
using RDatasets

rtol = 0.0001

@testset "OLS perfect fit" begin
    X = rand(1_000, 5)
    β = [ 1; 2; 3; 4; 5]
    y = X*β

    my_model = Jeeves.OLSModel(y, X, vcov = Jeeves.vcovIID())
    fitted_model = fit(my_model)
    model_coefs = coef(fitted_model)
    @test isapprox(model_coefs[1], β[1], rtol = rtol)
    @test isapprox(model_coefs[2], β[2], rtol = rtol)
    @test isapprox(model_coefs[3], β[3], rtol = rtol)
    @test isapprox(model_coefs[4], β[4], rtol = rtol)
    @test isapprox(model_coefs[5], β[5], rtol = rtol)
end

@testset "OLS perfect fit, default vcov" begin
    X = rand(1_000, 5)
    β = [ 1; 2; 3; 4; 5]
    y = X*β

    my_model = Jeeves.OLSModel(y, X)
    fitted_model = fit(my_model)
    inference = Jeeves.inference(fitted_model, vcovIID())
    SEs = inference[1]
    model_coefs = coef(fitted_model)
    @test isapprox(model_coefs[1], β[1], rtol = rtol)
    @test isapprox(model_coefs[2], β[2], rtol = rtol)
    @test isapprox(model_coefs[3], β[3], rtol = rtol)
    @test isapprox(model_coefs[4], β[4], rtol = rtol)
    @test isapprox(model_coefs[5], β[5], rtol = rtol)

    @test isless(SEs[1], 1e-10)
    @test isless(SEs[2], 1e-10)
    @test isless(SEs[3], 1e-10)
    @test isless(SEs[4], 1e-10)
    @test isless(SEs[5], 1e-10)
end

rtol = 0.0001
@testset "OLS perfect fit, DataFrame version" begin
    X = rand(1_000, 5)
    β = [ 1; 2; 3; 4; 5]
    y = X*β

    regression_df =  DataFrame(hcat(y, X), ["y", "x1", "x2", "x3", "x4", "x5"])
    my_model = Jeeves.OLSModel(regression_df[!, "y"], regression_df[!, ["x1", "x2", "x3", "x4", "x5"]], vcov = Jeeves.vcovIID())

    fitted_model = fit(my_model)
    model_coefs = coef(fitted_model)
    @test isapprox(model_coefs[1], β[1], rtol = rtol)
    @test isapprox(model_coefs[2], β[2], rtol = rtol)
    @test isapprox(model_coefs[3], β[3], rtol = rtol)
    @test isapprox(model_coefs[4], β[4], rtol = rtol)
    @test isapprox(model_coefs[5], β[5], rtol = rtol)
end

@testset "TSLS and Anderson Rubin works" begin
    N = 10_000
    confounder = randn(N)
    X_exog = randn(N, 3)
    X_beta = [1; 2; 3]
    β_endog = [1; -1]
    Z = randn(N, 2)
    X_endog = randn(N, 2) .+ confounder .+ Z
    Y = X_exog * X_beta + X_endog*β_endog + confounder + randn(N, 1)
    IV_model = Jeeves.TSLSModel(Y[:], X_endog, X_exog, Z)
    IV_fit = fit(IV_model)
    separate_inference = Jeeves.inference(IV_fit, vcovIID())
    model_coefs = coef(IV_fit)
    @test isapprox(model_coefs[1], 1, rtol = 1e-1)
    @test isapprox(model_coefs[2], -1, rtol = 1e-1)
    @test isapprox(model_coefs[3], 1, rtol = 1e-1)
    @test isapprox(model_coefs[4], 2, rtol = 1e-1)
    @test isapprox(model_coefs[5], 3, rtol = 1e-1)

    # Testing AndersonRubinCI returns results
    β_grid = Iterators.product(-1.:0.5:0., -1.:0.5:0.)
    AR_test = Jeeves.AndersonRubinCI(β_grid, IV_fit)
    @test size(AR_test, 1) == 9
    @test size(AR_test[1], 1) == 5
end

@testset "JIVE works" begin
    N = 10_000
    confounder = randn(N)
    X_exog = randn(N, 3)
    X_beta = [1; 2; 3]
    β_endog = [1; -1]
    Z = randn(N, 2)
    X_endog = randn(N, 2) .+ confounder .+ Z
    Y = X_exog * X_beta + X_endog*β_endog + confounder + randn(N, 1)
    IV_model = Jeeves.TSLSModel(Y[:], X_endog, X_exog, Z)
    JIVE_fit = Jeeves.jive(IV_model)

    JIVE_fit.modelfit.se_β

    separate_inference = Jeeves.inference(JIVE_fit, vcovIID())
    @test isequal(JIVE_fit.modelfit.se_β, separate_inference[1])
    model_coefs = coef(JIVE_fit)
    @test isapprox(model_coefs[1], 1, rtol = 1e-1)
    @test isapprox(model_coefs[2], -1, rtol = 1e-1)
    @test isapprox(model_coefs[3], 1, rtol = 1e-1)
    @test isapprox(model_coefs[4], 2, rtol = 1e-1)
    @test isapprox(model_coefs[5], 3, rtol = 1e-1)
end

# Kinda hard to test these functions.
@testset "tidying" begin
    X = rand(1_000, 5)
    β = [ 1; 2; 3; 4; 5]
    y = X*β

    my_model = broadcast(x -> Jeeves.OLSModel(y, X), 1:5)
    fitted_models = fit.(my_model)
    single_tidy = tidy(fitted_models[1])
    many_tidy = tidy(fitted_models, collect(1:5))

    tex_table = Jeeves.TableCol("ed", fitted_models[1], drop = ["x_1"])
    @test length(tex_table.data) == 6
    @test typeof(single_tidy) == DataFrame
    @test typeof(many_tidy) == DataFrame
end

@testset "F-test matches t-stat squared" begin
    X = hcat(repeat([1.0], 1_000), rand(1_000, 1))
    β = [1, 5]
    y =  X*β + randn(1_000)
    model = Jeeves.OLSModel(y, X)
    ols_fit = fit(model)
    tidy_fit = tidy(ols_fit)
    F_stat, _, _ = Jeeves.Ftest(ols_fit)
    t_stat_sq = tidy_fit[2, "t_stat"]^2

    @test isapprox(F_stat, t_stat_sq)
end

# Literally just testing if we get an answer at this stage lol...
@testset "Clustered SEs Run" begin
    N = 4000
    cluster_var_1 = repeat(["Kentucky", "Mass", "Illinois", "Florida"], Int(N/4))
    cluster_var_2 = repeat(["t_1", "t_2", "t_3", "t_4", "t_5"],Int(N/5))
    cluster_var_3 = repeat(1:50, Int(N/50))

    using Random
    cluster_matrix = Matrix(hcat(shuffle(cluster_var_1),
                                 shuffle(cluster_var_2)))

    X = randn(N, 5)
    ϵ = randn(N, 1)
    β = [1, 2, 3, 4, 5]
    y = X * β    + ϵ

    model = Jeeves.OLSModel(y[:], X)
    model_fit = fit(model)
    SEs = Jeeves.inference(model_fit, Jeeves.vcovCluster(cluster_matrix))
    SEs_single_cluster = Jeeves.inference(model_fit, Jeeves.vcovCluster(cluster_matrix[:, 1:1]))
    @test typeof(SEs[1]) == Vector{Float64}
    @test typeof(SEs_single_cluster[1]) == Vector{Float64}
end

@testset "Hetero SEs run" begin
    N = 1000
    X = randn(N, 5)
    ϵ = randn(N, 1)
    β = [1, 2, 3, 4, 5]
    y = X * β    + ϵ

    model = Jeeves.OLSModel(y[:], X)
    model_fit = fit(model)
    SEs = Jeeves.inference(model_fit, Jeeves.vcovHetero())
    @test typeof(SEs[1]) == Vector{Float64}
end

@testset "Clustered SEs Match Sandwich" begin
    PublicSchools = dropmissing(dataset("sandwich", "PublicSchools"))
    PublicSchools[!, "fake_cluster"] .= repeat(1:5, 10)
    PublicSchools[!, "const"] .= 1.0
    PS_model = Jeeves.OLSModel(PublicSchools[!, "Expenditure"],
                               PublicSchools[!, ["const", "Income"]],
                               vcov = Jeeves.vcovCluster(
                                   Matrix(PublicSchools[!, ["State",
                                                            "fake_cluster"]])))
    PS_fit = fit(PS_model)
    PS_coef = coef(PS_fit)
    PS_SE = PS_fit.modelfit.se_β
    # Quick coef checks vs R package
    @test isapprox(-151.265090, PS_coef[1], rtol = rtol)
    @test isapprox(0.068939, PS_coef[2], rtol = rtol)

    # SEs
    # R's sandwich defaults to HC1 169.594852, 0.023020
    @test isapprox(169.594852, PS_SE[1], rtol = rtol)
    @test isapprox(0.023020, PS_SE[2], rtol = rtol)
end

@testset "Boot SEs Roughly Match IID when DGP IID" begin
    N = 4000
    cluster_var = repeat(1:50, Int(N/50))

    using Random

    X = randn(N, 5)
    ϵ = randn(N, 1)
    β = [1, 2, 3, 4, 5]
    y = X * β    + ϵ

    model = Jeeves.OLSModel(y[:], X)
    model_fit = fit(model)
    Jeeves.inference(model_fit, Jeeves.vcovBoot(1, cluster_var[:,:]))
    SEs_boot = Jeeves.inference(model_fit, Jeeves.vcovBoot(1000, cluster_var[:,:]))
    SEs = Jeeves.inference(model_fit, Jeeves.vcovIID())
    SE_comp = hcat(SEs[1], SEs_boot[1])

    @test typeof(SEs_boot[1]) == Vector{Float64}
    # SEs match
    for i in 1:5
        @test isapprox(SE_comp[i, 1], SE_comp[i, 2], atol = 0.05)
    end
    # Pvals match
    for i in 1:5
        @test isapprox(SEs[2][i], SEs_boot[2][i], atol = 0.05)
    end
end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
struct SimpleNesterovSequence{R} end

function Base.iterate(::SimpleNesterovSequence{R}, k=R(0)) where R
    val = iszero(k) ? R(0) : (k - 1) / (k + 2)
    return val, k + 1
end

struct NesterovSequence{R} end

function Base.iterate(::NesterovSequence{R}, thetas=(R(1), R(1))) where R
    theta_prev, theta = thetas
    return (theta_prev - 1) / theta, (theta, (1 + sqrt(1 + 4 * theta^2)) / 2)
end

struct NesterovExtrapolation{S} end

NesterovExtrapolation(S) = NesterovExtrapolation{S}()
NesterovExtrapolation() = NesterovExtrapolation{SimpleNesterovSequence}()

acceleration_style(::Type{<:NesterovExtrapolation}) = NesterovStyle()

initialize(::NesterovExtrapolation{S}, x) where S = Iterators.Stateful(S{real(eltype(x))}())


text:
@testsuite "indexing scalar" AT->begin
    AT <: AbstractGPUArray && @allowscalar @testset "errors and warnings" begin
        x = AT([0])

        allowscalar(true, false)
        x[1] = 1
        @test x[1] == 1

        @disallowscalar begin
            @test_throws ErrorException x[1]
            @test_throws ErrorException x[1] = 1
        end

        x[1] = 2
        @test x[1] == 2

        allowscalar(false)
        @test_throws ErrorException x[1]
        @test_throws ErrorException x[1] = 1

        @allowscalar begin
            x[1] = 3
            @test x[1] == 3
        end

        allowscalar() do
            x[1] = 4
            @test x[1] == 4
        end

        @test_throws ErrorException x[1]
        @test_throws ErrorException x[1] = 1

        allowscalar(true, false)
        x[1]

        allowscalar(true, true)
        @test_logs (:warn, r"Performing scalar operations on GPU arrays: .*") x[1]
        @test_logs x[1]

        # NOTE: this inner testset _needs_ to be wrapped with allowscalar
        #       to make sure its original value is restored.
    end

    @allowscalar @testset "getindex with $T" for T in supported_eltypes()
        x = rand(T, 32)
        src = AT(x)
        for (i, xi) in enumerate(x)
            @test src[i] == xi
        end
        @test Array(src[1:3]) == x[1:3]
        @test Array(src[3:end]) == x[3:end]
    end

    @allowscalar @testset "setindex! with $T" for T in supported_eltypes()
        x = fill(zero(T), 7)
        src = AT(x)
        for i = 1:7
            src[i] = i
        end
        @test Array(src) == T[1:7;]
        src[1:3] = T[77, 22, 11]
        @test Array(src[1:3]) == T[77, 22, 11]
        src[1] = T(0)
    end

    @allowscalar @testset "issue #42 with $T" for T in supported_eltypes()
        Ac = rand(Float32, 2, 2)
        A = AT(Ac)
        @test A[1] == Ac[1]
        @test A[end] == Ac[end]
        @test A[1, 1] == Ac[1, 1]
    end

    @allowscalar @testset "get/setindex!" begin
        # literal calls to get/setindex! have different return types
        @test compare(x->getindex(x,1), AT, zeros(Int, 2))
        @test compare(x->setindex!(x,1,1), AT, zeros(Int, 2))

        # issue #319
        @test compare(x->setindex!(x,1,1,1), AT, zeros(Float64, 2, 2))
    end
end

@testsuite "indexing multidimensional" AT->begin
    @testset "sliced setindex" for T in supported_eltypes()
        x = AT(zeros(T, (10, 10, 10, 10)))
        y = AT(rand(T, (5, 5, 10, 10)))
        x[2:6, 2:6, :, :] = y
        @test Array(x[2:6, 2:6, :, :]) == Array(y)
    end

    @testset "sliced setindex, CPU source" for T in supported_eltypes()
        x = AT(zeros(T, (2,3,4)))
        y = AT(rand(T, (2,3)))
        x[:, :, 2] = y
        @test Array(x[:, :, 2]) == Array(y)
    end

    @allowscalar @testset "empty array" begin
        @testset "1D" begin
            Ac = zeros(Float32, 10)
            A = AT(Ac)
            @test typeof(A[[]]) == typeof(AT(Ac[[]]))
            @test size(A[[]]) == size(Ac[[]])
        end

        @testset "2D with other index $other" for other in (Colon(), 1:5, 5)
            Ac = zeros(Float32, 10, 10)
            A = AT(Ac)

            @test typeof(A[[], other]) == typeof(AT(Ac[[], other]))
            @test size(A[[], other]) == size(Ac[[], other])

            @test typeof(A[other, []]) == typeof(AT(Ac[other, []]))
            @test size(A[other, []]) == size(Ac[other, []])
        end

        @test compare(AT, rand(Float32, 2,2)) do a
            a[:, 2:end-2] = AT{Float32}(undef,2,0)
        end
    end

    @testset "GPU source" begin
        a = rand(3)
        i = rand(1:3, 2)
        @test compare(getindex, AT, a, i)
        @test compare(getindex, AT, a, i')
        @test compare(getindex, AT, a, view(i, 2:2))
    end

    @testset "CPU source" begin
        # JuliaGPU/CUDA.jl#345
        a = rand(3,4)
        i = rand(1:3,2,2)
        @test compare(a->a[i,:], AT, a)
        @test compare(a->a[i',:], AT, a)
        @test compare(a->a[view(i,1,:),:], AT, a)
    end

    @testset "JuliaGPU/CUDA.jl#461: sliced setindex" begin
        @test compare((X,Y)->(X[1,:] = Y), AT, zeros(2,2), ones(2))
    end
end


text:
relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    # rootdet
    ((2, 2, MatNegRtdet()),),
    ((2, 2, MatNegRtdetEFExp()),),
    ((2, 2, MatNegRtdetEFPow()),),
    # tr negentropyconj (domain not positive)
    ((1, 2, MatNegEntropyConj()),),
    ((1, 2, MatNegEntropyConjEigOrd()),),
    # tr negsqrtconj
    ((1, 2, MatNegSqrtConj()),),
    ((1, 2, MatNegSqrtConjDirect()),),
    # tr negpower01conj
    ((1, 2, MatNegPower01Conj(0.3)),),
    ((1, 2, MatNegPower01ConjEigOrd(0.3)), nothing, relaxed_tols),
    # tr power12conj (domain not positive)
    ((1, 2, MatPower12Conj(1.7)),),
    ((1, 2, MatPower12ConjEigOrd(1.7)),),
    ]
insts["fast"] = [
    ((1, 10, MatNegRtdetEFExp()),),
    ((1, 15, MatNegRtdet()),),
    ((2, 3, MatNegRtdetEFPow()),),
    ((2, 3, MatNegSqrtEigOrd()),),
    ((2, 6, MatNegSqrt()),),
    ((2, 5, MatNegSqrtConjDirect()),),
    ((2, 6, MatNegSqrtConj()),),
    ((2, 6, MatNegEntropy()),),
    ((2, 7, MatLogdetCone()),),
    ((3, 2, MatNegEntropyEigOrd()),),
    ((3, 2, MatNegEntropyConj()),),
    ((3, 2, MatNegEntropyConjEigOrd()),),
    ((3, 2, MatNegPower01Conj(0.3)),),
    ((3, 2, MatNegPower01ConjEigOrd(0.3)),),
    ((3, 2, MatPower12EigOrd(1.5)),),
    ((3, 4, MatPower12(1.5)),),
    ((3, 4, MatNegRtdet()),),
    ((7, 2, MatNegLog()),),
    ((7, 2, MatPower12(1.5)),),
    ((7, 2, MatNegPower01(0.7)),),
    ]
insts["various"] = [
    ((2, 5, MatNegSqrtEigOrd()),),
    ((2, 5, MatNegLogEigOrd()),),
    ((2, 4, MatPower12EigOrd(1.5)), nothing, relaxed_tols),
    ((2, 4, MatNegEntropyEigOrd()),),
    ((2, 10, MatNegRtdetEFExp()),),
    ((2, 8, MatNegRtdetEFPow()),),
    ((2, 14, MatNegRtdet()),),
    ((2, 12, MatNegSqrtConj()),),
    ((2, 10, MatLogdetCone()),),
    ((2, 8, MatNegEntropy()),),
    ((2, 6, MatNegPower01(0.7)),),
    ((3, 5, MatNegLogDirect()),),
    ((3, 6, MatNegEntropy()),),
    ((3, 6, MatNegLog()),),
    ((3, 5, MatPower12(1.5)),),
    ((3, 5, MatNegEntropyConj()),),
    ((3, 4, MatNegPower01Conj(0.3)),),
    ((5, 2, MatNegEntropyConjEigOrd()), nothing, relaxed_tols),
    ((5, 2, MatPower12ConjEigOrd(1.5)),),
    ((5, 3, MatNegSqrtConjDirect()),),
    ((5, 3, MatNegLogDirect()),),
    ((6, 3, MatNegRtdet()),),
    ((6, 3, MatNegSqrt()),),
    ]
return (CentralPolyMatJuMP, insts)


text:
using Term
import Term: chars

prt(pan) = begin
    print(" " * hLine(pan.measure.w; style = "red"))
    print(vLine(pan.measure.h; style = "red") * pan)
    println(pan.measure, "   ", length(pan.segments))
end

print("\n"^5)

# check that unreasonable widhts are ignored
@time prt(TextBox("nofit"^25; width = 1000))

@time prt(TextBox("truncate"^25; width = 100, fit = :truncate))

@time prt(TextBox("truncate"^25; width = 100))

@time prt(TextBox("truncate"^8; fit = :fit))

@time prt(TextBox("[red]truncate[/red]"^8; fit = :fit))

@time prt(TextBox("[red]tru\nncate[/red]test"^1; fit = :fit))

@time prt(TextBox("[red]truncate[/red]test"^8; fit = :fit, justify = :left))


text:
module TestRelease

using Test
using Pkg
using IonBase
using UUIDs
using IonBase.ReleaseCmd

include(joinpath(pkgdir(IonBase), "test", "utils.jl"))

with_test_ion() do
    @testset "release" begin
        @test ReleaseCmd.is_version_number("0.1.0")
        @test ReleaseCmd.is_version_number("v0.1.0")

        basic_project = joinpath(test_dir, "Basic")
        project = ReleaseCmd.Project(basic_project, quiet=true)
        @test project.pkg.version == v"0.1.0"
        @test project.path == basic_project

        ReleaseCmd.update_version!(project, "0.2.0")
        @test project.pkg.version == v"0.2.0"
        @test Pkg.Types.read_project(joinpath(basic_project, "Project.toml")).version == v"0.2.0"

        ReleaseCmd.update_version!(project, "patch")
        @test project.pkg.version == v"0.2.1"
        @test Pkg.Types.read_project(joinpath(basic_project, "Project.toml")).version == v"0.2.1"

        ReleaseCmd.update_version!(project, "minor")
        @test project.pkg.version == v"0.3.0"
        @test Pkg.Types.read_project(joinpath(basic_project, "Project.toml")).version == v"0.3.0"

        ReleaseCmd.update_version!(project, "major")
        @test project.pkg.version == v"1.0.0"
        @test Pkg.Types.read_project(joinpath(basic_project, "Project.toml")).version == v"1.0.0"

        project = ReleaseCmd.Project(pkgdir(IonBase), quiet=true)
        matches = ReleaseCmd.query_project_registry(project)
        @test length(matches) == 1
        @test first(matches).uuid == UUID("23338594-aafe-5451-b93e-139f81909106")
    end
end

end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

#TODO: robustness, isopen checks, reconnection

FLAGS = [:_echo, :_info, :reset!, :actions, :observe, :act!, :terminated]
FMAP = Bijection(Dict(UInt8(i) => f for (i, f) in enumerate(FLAGS))) # maps flags <-> bytes

"""
    ASTClient

Solver-side client. Interacts with server via TCP.
"""
Base.@kwdef mutable struct ASTClient
    ip::IPAddr                      = ip"::1" # address(es) to connect to
    port::Int64                     = 1812    # port to connect to
    conn::Union{TCPSocket, Nothing} = nothing
    tunnel::Bool                    = false
    verbose::Bool                   = false
    caching::Bool                   = false   # accelerates episodic problems
    cache::Vector                   = []
end

"""
    RemoteASTMDP{S<:State, A<:Action} <: AbstractASTMDP{S, A}

MDP object that virtually references remote MDP. Autogenerated by `ASTClient`.
"""
Base.@kwdef mutable struct RemoteASTMDP{S<:State, A<:Action} <: AbstractASTMDP{S, A}
    client::ASTClient
    episodic::Bool   = false
    num_steps::Int64 = 0
end

"""
Exception type for server-side error. Allows exceptions to propogate to client to avoid
hanging connection.
"""
Base.@kwdef struct ASTServerError <: Exception
    msg::String = "ASTServer encountered an error."
end
Base.showerror(io::IO, e::ASTServerError) = print(io, e.msg)

"""
Perform function call on server MDP. Call sends request to server with provided function
and arguments and receives return value. Blocks until complete.
"""
function call(client::ASTClient, f::Function, args...; kwargs...)
    flag = Symbol(f)
    request = Dict(:f => FMAP(flag), :a => args, :k => NamedTuple(kwargs))
    bson(client.conn, request)
    client.verbose && @info "Sending request to server:" request
    response = BSON.load(client.conn)
    client.verbose && @info "Received response from server:" response
    haskey(response, :e) && response[:e] && throw(ASTServerError())
    return response[:r]
end

function reset!(mdp::RemoteASTMDP)
    mdp.client.caching && empty!(mdp.client.cache)
    sync(mdp.client)
    call(mdp.client, reset!)
end

function actions(mdp::RemoteASTMDP)
    action_type(mdp) == SeedAction ? UInt32 : Dirac(call(mdp.client, actions))
end

function observe(mdp::RemoteASTMDP)
    call(mdp.client, observe)
end

#TODO: check for off-by-one errors
function act!(mdp::RemoteASTMDP, action)
    if mdp.client.caching
        cache = mdp.client.cache
        done = length(cache) == mdp.num_steps
        push!(cache, action)
        return done ? call(mdp.client, act!, cache; batch=true) : 0.0
    else
        return call(mdp.client, act!, action)
    end
end

function terminated(mdp::RemoteASTMDP)
    mdp.client.caching ? length(mdp.client.cache) >= mdp.num_steps : call(mdp.client, terminated)
end

"""
    ping(client::ASTClient)

Request empty echo from ASTServer, returning approximate round-trip time in seconds.
"""
ping(client::ASTClient) = @elapsed call(client, _echo, nothing)

"""
    generate_mdp(client::ASTClient)

Query remote MDP and generate `RemoteASTMDP` with matching attributes to interface with
appropriate solvers.
"""
function generate_mdp(client::ASTClient)
    attr = call(client, _info)
    S = getproperty(Interface, attr[:S])
    A = getproperty(Interface, attr[:A])
    return RemoteASTMDP{S, A}(client, attr[:episodic], attr[:num_steps])
end

"""
Read from socket until communication is synchronized. Typically unnecessary, but mitigates
issue in which IO between client and server can become mismatched.
"""
function sync(client::ASTClient)
    seed = rand(RandomDevice(), UInt32)
    r = call(client, _echo, seed)
    while r != seed
        r = BSON.load(client.conn)[:r]
    end
    return
end

"""
    connect!(client::ASTClient; remote::String="", remote_port::Int64=1812)

Connect client to server, optionally through SSH tunnel. Optional keyword argument `remote`
should be of the form `user@machine`.
"""
function connect!(client::ASTClient; remote::String="", remote_port::Int64=1812)
    disconnect!(client)
    !isempty(remote) && open_tunnel(client, remote, remote_port)
    client.conn = connect(client.ip, client.port)
    Sockets.nagle(client.conn, false)
    t_ms = Int64(1000 * round(ping(client); digits=3))
    @info "ASTServer responded in $t_ms milliseconds."
    return
end

"""
    disconnect!(client::ASTClient)

Disconnect client from server.
"""
function disconnect!(client::ASTClient)
    if client.conn !== nothing
        close(client.conn)
        client.conn = nothing
    end
    client.tunnel && close_tunnel()
    client.tunnel = false
    return
end


text:
using Plots, Statistics, StatsBase, LaTeXStrings, JLD, ProgressBars
#mutable struct -> instances of the composite types can be modified.
mutable struct MDSystem{FT<:AbstractFloat}
    L::FT
    step::FT
    Eᵤ::FT
    Eₖ::FT
    T::FT
    P::FT
    N::Integer
    F::Matrix{FT}
    r::Matrix{FT}
    v::Matrix{FT}
    function MDSystem(N::Integer, T₀::FT, step::FT, L::FT) where {FT<:AbstractFloat}
        r = rand(FT, 2, N) * L
        v = rand(FT, 2, N) .- 0.5
        v *= sqrt((T₀ ÷ mean(v .^ 2)))
        F = Matrix{FT}(undef, 2, N)
        return new{FT}(L, step, FT(NaN), FT(NaN), FT(NaN), FT(NaN), N, F, r, v)
    end
end

function Initialize(N::Integer,L::FT, T₀::FT, step::FT) where {FT<:AbstractFloat}
    md = MDSystem(N,  T₀, step, L)
    alignleft(md)
    for i in 1:md.N
        Force_cal(md, i)
    end
    VelVerlet(md)
    Boundary_condition(md)
    Potential_cal(md)
    Kinetic_cal(md)
    Temprature_cal(md)
    Pressure_cal(md)
    return md
end

function Sim_process(md :: MDSystem)
    VelVerlet(md)
    Boundary_condition(md)
    Potential_cal(md)
    Kinetic_cal(md)
    Temprature_cal(md)
    Pressure_cal(md)
end

function alignleft(md::MDSystem)
    xcount = ceil(Integer, (2 * md.N)^(1 / 2) / 2)
    ycount = 2 * xcount

    xs = collect(range(0.1 * md.L / 2, 0.9 * md.L / 2, length = xcount))
    md.r[1, :] .= repeat(xs, outer = ycount)[1:md.N]

    ys = collect(range(0.1 * md.L, 0.9 * md.L, length = ycount))
    md.r[2, :] .= repeat(ys, inner = xcount, outer = 1)[1:md.N]
end

function Force_cal(md :: MDSystem, i :: Integer)
    md.F[:, i] *= 0
    for j in setdiff(1:md.N, i)
        rp = md.r[:, j]
        Δx = md.r[1, i] - rp[1]
        Δy = md.r[2, i] - rp[2]
        if Δx < -md.L / 2
            rp[1] -= md.L
        elseif Δx > md.L / 2
            rp[1] += md.L
        end
        if Δy < -md.L / 2
            rp[2] -= md.L
        elseif Δy > md.L / 2
            rp[2] += md.L
        end
        Δr = sqrt(sum((md.r[:, i] - rp) .^ 2))
        md.F[:, i] += 48 * (((Δr)^-14) - 0.5 * ((Δr)^-8)) * (md.r[:, i] - rp)
    end
end

function VelVerlet(md :: MDSystem)
    for i in 1:md.N
        a₁ = md.F[:, i]
        md.r[:, i] += md.step * md.v[:, i] + 1/2 * md.step  * md.step * a₁
        Force_cal(md, i)
        a₂ = md.F[:, i]
        md.v[:, i] += md.step * (a₁ + a₂) / 2
    end
end

function Boundary_condition(md :: MDSystem)
    md.r = (md.r .+ md.L) .% md.L
end

function Potential_cal(md :: MDSystem)
    ΣΔr = 0.0
    for i in 1:(md.N - 1)
        for j in (i+1):md.N
            rp = md.r[:, j]
            Δx = md.r[1, i] - rp[1]
            Δy = md.r[2, i] - rp[2]
            if Δx < -md.L / 2
                rp[1] -= md.L
            elseif Δx > md.L / 2
                rp[1] += md.L
            end
            if Δy < -md.L / 2
                rp[2] -= md.L
            elseif Δy > md.L / 2
                rp[2] += md.L
            end
            Δr = sqrt(sum((md.r[:, i] - rp) .^ 2))
            ΣΔr += ((Δr) ^ (-12)) - ((Δr) ^ (-6))
        end
    end
    md.Eᵤ = 4 * ΣΔr
end

function Kinetic_cal(md :: MDSystem)
    md.Eₖ = sum(md.v .^2) / 2 #kinetic Energy per molecule
end

function Temprature_cal(md :: MDSystem)
    md.T = md.Eₖ / md.N
end

function Pressure_cal(md :: MDSystem)
    Σrterm = 0.0
    Σvᵢ² = sum(md.v .^2)
    for i in 1:(md.N-1)
        for j in (i+1):md.N
            rp = md.r[:, j]
            Δx = md.r[1, i] - rp[1]
            Δy = md.r[2, i] - rp[2]
            if Δx < -md.L / 2
                rp[1] -= md.L
            elseif Δx > md.L / 2
                rp[1] += md.L
            end
            if Δy < -md.L / 2
                rp[2] -= md.L
            elseif Δy > md.L / 2
                rp[2] += md.L
            end
            Δr = sqrt(sum(md.r[:, i] - rp) .^ 2)
            Σrterm += ((Δr) ^ (-12)) - (1/2) * ((Δr) ^ (-6))
        end
    end
    md.P = (Σvᵢ² + 48 * Σrterm)/(2 * md.L * md.L)
end

N = 100; T₀ = 1.0; h = 0.001; L = 50.0; step_num = 20000
md = Initialize(N, L, T₀ , h)

Poses = zeros(step_num, 2, N)
Velocities = zeros(step_num, 2, N)
Potentials = zeros(step_num)
Kinetics = zeros(step_num)
Temps = zeros(step_num)
pressures = zeros(step_num)

for i in ProgressBar(1:step_num)
    Sim_process(md)
    Poses[i, :, :] = md.r
    Velocities[i, :, :] = md.v
    Potentials[i] = md.Eᵤ
    Kinetics[i] = md.Eₖ
    Temps[i] = md.T
    pressures[i] = md.P
end
Poses
save("../../MD/SingleData.jld",
    "Poses", Poses, "Velocities", Velocities, "Potentials", Potentials,
    "Kinetics", Kinetics, "Temps", Temps, "Pressures", pressures)

plot(collect(0:h:(h*step_num))[1:20000], [count(<=(md.L / 2), Poses[i, 1, :]) / md.N for i ∈ 1:step_num][1:20000],
    ylabel=L"N_l_e_f_t/N", xlabel=L"t/\tau", dpi=400, title=L"Molecules'\ Densities\ at\ the\ left\ half", legend=false)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\MD\\Figs\\LeftHalfDensity.png")

plot(collect(0:h:(h*step_num))[1:2000],Kinetics[1:2000], label = L"E_k")
plot!(collect(0:h:(h*step_num))[1:2000],Potentials[1:2000], label = L"E_U")
plot!(collect(0:h:(h*step_num))[1:2000],Potentials[1:2000]+Kinetics[1:2000], label = L"E_t_o_t_a_l"
    , title = L"Enery\ of\ the\ System",dpi = 400)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\MD\\Figs\\Energies.png")

plot(collect(0:h:(h*step_num))[1:20000],120*Temps[1:20000], label = L"Temprature", title=L"Temprature\ of\ the\ System")
plot!([0.0, 20.0], 120 * [mean(Temps[10000:end]), mean(Temps[10000:end])],
    label = L"\langle T_{eq} \rangle = %$(round(120 * mean(Temps[10000:end]), digits = 3))")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\MD\\Figs\\Temprature.png")


text:
# This file is a part of TextClassification.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export EntModelConfigSpace, VectorModelConfigSpace

@with_kw struct EntModelConfig{W_<:LocalWeighting}
    local_weighting::W_ = TfWeighting()
    mindocs::Int = 1.0
    smooth::Float64 = 0.0
    keeptop::Float64 = 1.0
    weights::Union{Symbol,Vector{Float64}} = :balance
end

function create(c::EntModelConfig, train_X::AbstractVector{BOW}, train_y::CategoricalArray)
    model = VectorModel(EntropyWeighting(), c.local_weighting, train_X, train_y, mindocs=c.mindocs, smooth=c.smooth, weights=c.weights)
    c.keeptop < 1.0 ? prune_select_top(model, c.keeptop) : model
end

@with_kw struct EntModelConfigSpace <: AbstractSolutionSpace
    local_weighting = [TfWeighting(), TpWeighting(), FreqWeighting(), BinaryLocalWeighting()]
    mindocs = 1:5
    smooth = 0:3 # by default=1 in favor of mindocs
    keeptop = [0.5, 1.0]
    weights = [:balance, :none]
    scale_mindocs = (lower=1, s=1.3, upper=11)
    scale_smooth = (lower=0.0, s=1.3, upper=11)
    scale_keeptop = (lower=0.01, s=1.1, upper=1.0)
end

Base.eltype(::EntModelConfigSpace) = EntModelConfig

function Base.rand(space::EntModelConfigSpace)
    EntModelConfig(
        rand(space.local_weighting),
        rand(space.mindocs),
        rand(space.smooth),
        rand(space.keeptop),
        rand(space.weights)
    )
end

function combine(a::EntModelConfig, b::EntModelConfig)
    EntModelConfig(
        a.local_weighting,
        div(a.mindocs + b.mindocs, 2),
        (a.smooth + b.smooth) / 2,
        (a.keeptop + b.keeptop) / 2,
        b.weights
    )
end

function mutate(space::EntModelConfigSpace, c::EntModelConfig, iter)
    mindocs = space.scale_mindocs === nothing ? c.mindocs : SearchModels.scale(c.mindocs; space.scale_mindocs...)
    smooth = space.scale_smooth === nothing ? c.smooth : SearchModels.scale(c.smooth; space.scale_smooth...)
    keeptop = space.scale_keeptop === nothing ? c.keeptop : SearchModels.scale(c.keeptop; space.scale_keeptop...)
    lw = SearchModels.change(c.local_weighting, space.local_weighting, p1=0.3)

    EntModelConfig(
        lw,
        mindocs,
        smooth,
        keeptop,
        c.weights
    )
end

@with_kw struct VectorModelConfig{L_<:LocalWeighting, G_<:GlobalWeighting}
    global_weighting::G_ = IdfWeighting()
    local_weighting::L_ = TfWeighting()
    mindocs = 1
    keeptop::Float64 = 1.0
end

function create(c::VectorModelConfig, train_X::AbstractVector{BOW}, train_y::CategoricalArray)
    model = VectorModel(c.global_weighting, c.local_weighting, train_X)
    c.keeptop < 1.0 ? prune_select_top(model, c.keeptop) : model
end

@with_kw struct VectorModelConfigSpace <: AbstractSolutionSpace
    global_weighting = [IdfWeighting(), BinaryGlobalWeighting()]
    local_weighting = [TfWeighting(), TpWeighting(), FreqWeighting(), BinaryLocalWeighting()]
    mindocs = [1, 4, 7]
    keeptop = [0.5, 1.0]
    scale_mindocs = (lower=1, s=1.3, upper=30)
    scale_keeptop = (lower=0.01, s=1.3, upper=1.0)
end

Base.eltype(::VectorModelConfigSpace) = VectorModelConfig

function Base.rand(space::VectorModelConfigSpace)
    VectorModelConfig(
        rand(space.global_weighting),
        rand(space.local_weighting),
        rand(space.mindocs),
        rand(space.keeptop)
    )
end

function combine(a::VectorModelConfig, b::VectorModelConfig)
    VectorModelConfig(
        b.global_weighting,
        a.local_weighting,
        div(a.mindocs + b.mindocs, 2),
        (a.keeptop + b.keeptop) / 2
    )
end

function mutate(space::VectorModelConfigSpace, c::VectorModelConfig, iter)
    mindocs = SearchModels.scale(c.mindocs; space.scale_mindocs...)
    keeptop = SearchModels.scale(c.keeptop; space.scale_keeptop...)
    lw = SearchModels.change(c.local_weighting, space.local_weighting, p1=0.3)
    gw = SearchModels.change(c.global_weighting, space.global_weighting, p1=0.3)

    VectorModelConfig(gw, lw, mindocs, keeptop)
end


text:
# auxiliary.jl

# getindex
for TM in (:TaylorModel1, :RTaylorModel1, :TaylorModelN)
    @eval begin
        copy(f::$TM) = $TM(copy(f.pol), f.rem, f.x0, f.I)

        getindex(a::$TM, n::Integer) = a.pol[n]
        getindex(a::$TM, u::UnitRange) = view(a.pol, u .+ 1 )
        getindex(a::$TM, c::Colon) = view(a.pol, c)
    end
end

for TM in tupleTMs
    @eval begin
        setindex!(a::$TM{T,S}, x::T, n::Integer) where {T<:Number, S} = a.pol[n] = x
        setindex!(a::$TM{T,S}, x::T, u::UnitRange) where {T<:Number, S} = a.pol[u] .= x
        function setindex!(a::$TM{T,S}, x::Array{T,1}, u::UnitRange) where {T<:Number, S}
            @assert length(u) == length(x)
            for ind in eachindex(x)
                a.pol[u[ind]] = x[ind]
            end
        end
        setindex!(a::$TM{T,S}, x::T, c::Colon) where {T<:Number, S} = a[c] .= x
        setindex!(a::$TM{T,S}, x::Array{T,1}, c::Colon) where {T<:Number, S} = a[c] .= x
    end
end


text:
__precompile__()

module AstrobeeRobot

using RigidBodyDynamics
using StaticArrays

export Astrobee,
  packagepath,
  urdfpath

packagepath() = joinpath(@__DIR__, "..", "deps")
urdfpath() = joinpath(packagepath(), "Astrobee", "astrobee.urdf")

mutable struct Astrobee{T}
  mechanism::Mechanism{T}
  basejoint::Joint
  body::RigidBody

  function Astrobee{T}(;floating=true) where T
    mechanism = RigidBodyDynamics.parse_urdf(T,urdfpath())

    #create floating joint
    body = findbody(mechanism, "body")
    basejoint = joint_to_parent(body,mechanism)
    if floating
     println("instantiated a floating joint")
     floatingjoint = Joint(basejoint.name, frame_before(basejoint), frame_after(basejoint), QuaternionFloating{T}())
     replace_joint!(mechanism, basejoint, floatingjoint)
     basejoint=floatingjoint
   end

  new{T}(mechanism, basejoint, body)
  end
end
Astrobee(::Type{T} = Float64; kwargs...) where {T} = Astrobee{T}(; kwargs...)

function __init__()
  if !isfile(urdfpath())
    error("Could not find $(urdfpath()). Please run `import Pkg; Pkg.build(\"Astrobee\")`.")
  end
end

end # module


text:
using Printf

struct Shell
  shell_id::Int64
  atom_id::Int64
  atomic_number::Int64

  exponents::SVector{MAX_CONTRACTION,Float64}
  coefficients::SVector{2*MAX_CONTRACTION,Float64}

  atom_center::SVector{3,Float64}

  class::Cstring
  am::Int64
  nbas::Int64
  nprim::Int64
  pos::Int64
  sp::Int64
end
export Shell

Shell(shell_id, atom_id, atomic_number, exponents,
  coefficients, atom_center, am, nprim,
  pos, unnormalize) = Shell(shell_id, atom_id, atomic_number,
  create_static_vector_small(exponents),
  calculate_coefficients(coefficients,exponents,am,nprim,
  am_to_nbas_cart(am),unnormalize),
  atom_center, shell_class(am), am, am_to_nbas_cart(am), nprim, pos,
  Int64(am_to_nbas_cart(am) == 4))

@inline function am_to_nbas_cart(am::Int64)
  return (am == -1) ? 4 : Int64(am*(am+1)/2)
end
export am_to_nbas_cart

@inline function create_static_vector_small(
  input::Vector{T}) where {T<:Number}

  output_dynamic::Vector{T} = Vector{T}(undef,MAX_CONTRACTION)
  @views output_dynamic[1:size(input)[1]] = input[:]
  return SVector{MAX_CONTRACTION,T}(output_dynamic)
end

@inline function create_static_vector_large(
  input::Vector{T}) where {T<:Number}

  output_dynamic::Vector{T} = Vector{T}(undef,2*MAX_CONTRACTION)
  @views output_dynamic[1:size(input)[1]] = input[:]
  return SVector{2*MAX_CONTRACTION,T}(output_dynamic)
end

@inline function calculate_coefficients(
  coef::Vector{T}, exp::Vector{T}, am::Int64, nprim::Int64, nbas::Int64,
  unnormalize::Bool) where {T<:Number}

  #== unnormalize basis functions ==#
  #println("UNNORMALIZE")
  coef_factor::Float64 = 0.0
  if unnormalize
    for iprim::Int64 in 1:nprim
      two_exp::Float64 = 2*exp[iprim]
      coef_factor = (pi/two_exp)^1.5
      coef_factor *= am == 2 ? 0.5/two_exp : 1.0
      coef_factor *= am == 3 ? 0.75/(two_exp^2) : 1.0
      coef_factor *= am == 4 ? 1.875/(two_exp^3) : 1.0
      coef[iprim] /= sqrt(coef_factor)
    end

    if nbas == 4  #account for sp shells
      for iprim::Int64 in (nprim+1):(2*nprim)
        two_exp::Float64 = 2*exp[iprim-nprim]
        unnorm = (pi/two_exp)^1.5
        unnorm *=  0.5/two_exp
        coef[iprim] /= sqrt(coef_factor)
      end
    end
  end

  #for icoef in coef
  #  @printf("%5.15f\n",icoef)
  #end
  #println("")

  #println("RENORMALIZE")
  #== renormalize basis functions ==#
  coef_factor = 0.0
  for iprim::Int64 in 1:nprim, jprim::Int64 in 1:iprim
    exp_j = exp[iprim] + exp[jprim]
    fac_j = exp_j^1.5
    dummy_coef_factor = coef[iprim] * coef[jprim]/fac_j
    dummy_coef_factor *= am == 2 ? 0.5/exp_j : 1.0
    dummy_coef_factor *= am == 3 ? 0.75/(exp_j^2.0) : 1.0
    dummy_coef_factor *= am == 4 ? 1.875/(exp_j^3.0) : 1.0
    if iprim != jprim
      dummy_coef_factor *= 2.0
    end
    coef_factor += dummy_coef_factor
  end

  if coef_factor > 1E-10 coef_factor = 1.0/sqrt(coef_factor*(pi^1.5)) end

  for icoef in coef
    icoef *= coef_factor
  #  @printf("%5.15f\n",icoef)
  end
  #println("")

  return create_static_vector_large(coef)
end

function shell_class(am)
  sh_class = ""
  if am == -1 sh_class = "L"
  elseif am == 1 sh_class = "s"
  elseif am == 2 sh_class = "p"
  elseif am == 3 sh_class = "d"
  elseif am == 4 sh_class = "f"
  end

  return pointer(sh_class)
end

#=
struct ShPair
  sh_a::Shell
  sh_b::Shell

  class::Cstring
  am2::Int64
  nbas2::Int64
  nprim2::Int64
end
export ShPair

ShPair(sh_a,sh_b) = ShPair(sh_a, sh_b, shpair_class(sh_a, sh_b),
  sh_a.am+sh_b.am, sh_a.nbas*sh_b.nbas, sh_a.nprim*sh_b.nprim)

function shpair_class(sh_a, sh_b)
  class_a = unsafe_string(sh_a.class)
  class_b = unsafe_string(sh_b.class)

  shpair_class = ""
  if class_a == "s" && class_b == "s" shpair_class = "ss"
  elseif class_a == "p" && class_b == "s" shpair_class = "ps"
  elseif class_a == "s" && class_b == "p" shpair_class = "sp"
  elseif class_a == "d" && class_b == "s" shpair_class = "ds"
  elseif class_a == "s" && class_b == "d" shpair_class = "sd"
  elseif class_a == "f" && class_b == "s" shpair_class = "fs"
  elseif class_a == "s" && class_b == "f" shpair_class = "sf"
  elseif class_a == "p" && class_b == "p" shpair_class = "pp"
  elseif class_a == "d" && class_b == "p" shpair_class = "dp"
  elseif class_a == "p" && class_b == "d" shpair_class = "pd"
  elseif class_a == "f" && class_b == "p" shpair_class = "fp"
  elseif class_a == "p" && class_b == "f" shpair_class = "pf"
  elseif class_a == "d" && class_b == "d" shpair_class = "dd"
  elseif class_a == "f" && class_b == "d" shpair_class = "fd"
  elseif class_a == "d" && class_b == "f" shpair_class = "df"
  elseif class_a == "f" && class_b == "f" shpair_class = "ff"
  end

  return pointer(shpair_class)
end

struct ShQuartet
  bra::ShPair
  ket::ShPair
end
export ShQuartet
=#
struct Basis
  shells::Vector{Shell}
  basis_cxx::JERI.BasisSet
  shpdata_cxx::StdVector{JERI.ShellPair}

  model::String
  norb::Int64
  nels::Int64
end
export Basis

function max_ang_mom(basis_set::Basis)
  max_am = 0
  for shell in basis_set.shells
    max_am = shell.am > max_am ? shell.am : max_am
  end
  return max_am
end
export max_ang_mom

function Base.getindex(basis_set::Basis, index)
  return basis_set.shells[index]
end

function Base.length(basis_set::Basis)
  return length(basis_set.shells)
end

function Base.iterate(basis_set::Basis)
  return iterate(basis_set.shells)
end

function Base.iterate(basis_set::Basis, state)
  return iterate(basis_set.shells, state)
end


text:
#!/usr/local/bin/julia

using ArgParse

function main(args)
	# Argument parsing
	configs = ArgParseSettings()

	add_arg_group!(configs, "Options:")
	@add_arg_table! configs begin
		"--input", "-i"
			help	 = "Input filename"
			action	 = :store_arg
			arg_type = String
		"--output", "-o"
			help	 = "Output filename"
			action	 = :store_arg
			arg_type = String
		end
	# Store arguments to variables
	parsed_args = parse_args(args, configs)
	infile      = parsed_args["input"]
	outfile     = parsed_args["output"]

	# Read file and compute problem via function
	input = chomp(open(f->read(f, String), infile))
	@show output =
	open(f->write(f, "$output"), outfile, "w+")

end

main(ARGS)


text:
#=
# MNIST

We begin by importing the required packages.
We load MNIST via the MLDatasets.jl package.
=#

import Makie
import CairoMakie
import MLDatasets
import Flux
import RestrictedBoltzmannMachines as RBMs
using Statistics: mean, std, var
using Random: bitrand
using ValueHistories: MVHistory
using RestrictedBoltzmannMachines: visible, BinaryRBM, sample_from_inputs, minibatch_count
using RestrictedBoltzmannMachines: initialize!, log_pseudolikelihood, pcd!
nothing #hide

# Useful function to plot grids of MNIST digits.

"""
    imggrid(A)

Given a four dimensional tensor `A` of size `(width, height, ncols, nrows)`
containing `width x height` images in a grid of `nrows x ncols`, this returns
a matrix of size `(width * ncols, height * nrows)`, that can be plotted in a heatmap
to display all images.
"""
imggrid(A::AbstractArray{<:Any,4}) =
    reshape(permutedims(A, (1,3,2,4)), size(A,1)*size(A,3), size(A,2)*size(A,4))

#=
Load the MNIST dataset.
We will train an RBM with binary (0,1) visible and hidden units.
Therefore we binarize the data.
In addition, we consider only one kind of digit so that training is faster.
=#

Float = Float32
train_x, train_y = MLDatasets.MNIST.traindata()
train_x = Array{Float}(train_x[:, :, train_y .== 0] .≥ 0.5)
nothing #hide

# Let's visualize some random digits.

nrows, ncols = 10, 15
fig = Makie.Figure(resolution=(40ncols, 40nrows))
ax = Makie.Axis(fig[1,1], yreversed=true)
idx = rand(1:size(train_x,3), nrows * ncols) # random indices of digits
digits = reshape(train_x[:,:,idx], 28, 28, ncols, nrows)
Makie.image!(ax, imggrid(digits), colorrange=(1,0))
Makie.hidedecorations!(ax)
Makie.hidespines!(ax)
fig

# Initialize an RBM with 400 hidden units.

rbm = BinaryRBM(Float, (28,28), 400)
initialize!(rbm, train_x) # match single-site statistics
nothing #hide

#=
Initially, the RBM assigns a poor pseudolikelihood to the data.
=#

println("log(PL) = ", mean(@time log_pseudolikelihood(rbm, train_x)))

#=
(Incidentally, note how long it takes to evaluate the pseudolikelihood on the full
dataset.)

Now we train the RBM on the data.
This returns a [MVHistory](https://github.com/JuliaML/ValueHistories.jl)
collecting some info during training.
=#

batchsize = 256
batchcount = minibatch_count(train_x; batchsize)
epochs = 500
history = MVHistory()
@time pcd!(
    rbm, train_x; epochs, batchsize,
    callback = function(; rbm, epoch, batch_idx, _...)
        if batch_idx == batchcount && epoch % 5 == 0
            push!(history, :lpl, mean(log_pseudolikelihood(rbm, train_x)))
        end
    end
)
nothing #hide

# After training, the pseudolikelihood score of the data improves significantly.
# Plot of log-pseudolikelihood of trian data during learning.

fig = Makie.Figure(resolution=(500,300))
ax = Makie.Axis(fig[1,1], xlabel = "train time", ylabel="pseudolikelihood")
Makie.lines!(ax, get(history, :lpl)...)
fig

# Sample digits from the RBM starting from a random condition.

nsteps = 3000
fantasy_F = zeros(nrows*ncols, nsteps)
fantasy_x = bitrand(28,28,nrows*ncols)
fantasy_F[:,1] .= RBMs.free_energy(rbm, fantasy_x)
@time for t in 2:nsteps
    fantasy_x .= RBMs.sample_v_from_v(rbm, fantasy_x)
    fantasy_F[:,t] .= RBMs.free_energy(rbm, fantasy_x)
end
nothing #hide

# Check equilibration of sampling

fig = Makie.Figure(resolution=(400,300))
ax = Makie.Axis(fig[1,1], xlabel="sampling time", ylabel="free energy")
fantasy_F_μ = vec(mean(fantasy_F; dims=1))
fantasy_F_σ = vec(std(fantasy_F; dims=1))
Makie.band!(ax, 1:nsteps, fantasy_F_μ - fantasy_F_σ/2, fantasy_F_μ + fantasy_F_σ/2)
Makie.lines!(ax, 1:nsteps, fantasy_F_μ)
fig

# Plot the sampled digits.

fig = Makie.Figure(resolution=(40ncols, 40nrows))
ax = Makie.Axis(fig[1,1], yreversed=true)
Makie.image!(ax, imggrid(reshape(fantasy_x, 28, 28, ncols, nrows)), colorrange=(1,0))
Makie.hidedecorations!(ax)
Makie.hidespines!(ax)
fig