text:
module Poppler

import Base.open ## Extends open

const libpoppler = "libpoppler-glib.so"

typealias FileURI ASCIIString

type PopplerDocument
    ptr::Ptr{Cint}
end

type PopplerPage
    ptr::Ptr{Cint}
end

## PopplerDocument — Information about a document
open(uri::FileURI) = open(uri, "") ## WORKS
function open(uri::FileURI, password::ASCIIString) ## WORKS
    uri = ismatch(r"^file://", uri) ? uri : "file://$uri"
    document = ccall((:poppler_document_new_from_file, libpoppler),
                     PopplerDocument, ## FIXME: type document
                     (Cstring, Cint, Cint),
                     uri, 0, 0) ## FIXME: Error
    document.ptr == C_NULL ? error("Could not open file") : document
end

function pdfversion(document::PopplerDocument) ## WORKS
    version = bytestring(ccall((:poppler_document_get_pdf_version_string, libpoppler),
                               Ptr{UInt8},
                               (PopplerDocument,),
                               document))
end

function title(document) ## WORKS
    title = ccall((:poppler_document_get_title, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)
    title != C_NULL ? bytestring(title) : ""
end

function author(document) ## WORKS
    author = ccall((:poppler_document_get_author, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)
    author != C_NULL ? bytestring(author) : ""
end

function subject(document) ## WORKS
    subject = ccall((:poppler_document_get_subject, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)
    subject != C_NULL ? bytestring(subject) : ""
end

function keywords(document) ## DON'T KNOW IF WORKS
    keywords = ccall((:poppler_document_get_keywords, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)
    keywords != C_NULL ? bytestring(keywords) : ""
end

function creator(document) ## WORKS
    creator = ccall((:poppler_document_get_creator, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)
    creator != C_NULL ? bytestring(creator) : ""
end

function producer(document) ## WORKS
    producer = ccall((:poppler_document_get_producer, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)
    producer != C_NULL ? bytestring(producer) : ""
end

## WORKS
creationdate(document) = Dates.unix2datetime(Int(ccall((:poppler_document_get_creation_date, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)))

## WORKS
modificationdate(document) = Dates.unix2datetime(Int(ccall((:poppler_document_get_modification_date, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)))

function metadata(document) ## WORKS
    metadata = ccall((:poppler_document_get_metadata, libpoppler), Ptr{UInt8}, (PopplerDocument,), document)
    metadata != C_NULL ? bytestring(metadata) : ""
end

function islinearized(document) ## TO BE TESTED
    Bool(ccall((:poppler_document_is_linearized, libpoppler), Ptr{Bool}, (PopplerDocument,), document))
end

function hasattachments(document) ## TO BE TESTED
    Bool(ccall((:poppler_document_has_attachments, libpoppler), Ptr{Bool}, (PopplerDocument,), document))
end

## WORKS
pages(document) = Int(ccall((:poppler_document_get_n_pages, libpoppler), Ptr{Int64}, (PopplerDocument,), document))

"""
   page(document, index)

Returns a page object. Indexing in julia mode (starts from 1)
"""
function page(document::PopplerDocument, idx::Int) ## WORKS
    ccall((:poppler_document_get_page, libpoppler),
          PopplerPage,
          (PopplerDocument, Cint),
          document, idx - 1)
end

pagecontent(page::PopplerPage) = bytestring(ccall((:poppler_page_get_text, libpoppler), Ptr{UInt8}, (PopplerPage,), page))
pagecontent(document::PopplerDocument, pagen::Int) = pagecontent(page(document, pagen))

"""
   text(document)

Returns the text of all the pages of the document
"""
function text(document) ## FIXME: overload document object/text
    content = []
    for pagen in 1:pages(document)
        push!(content, pagecontent(document, pagen))
    end
    content
end

export open, pdfversion, title, author, subject, keywords, creator, producer, creationdate, modificationdate, metadata, islinearized, hasattachments, pages, page, pagecontent, text

end # module


text:
"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, D},
    ::Type{S},
) where {T <: PSY.ThermalGen, D <: AbstractThermalFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    reactivepower_variables!(psi_container, devices)
    commitment_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, D)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactivepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    commitment_constraints!(psi_container, devices, model, S, get_feedforward(model))
    ramp_constraints!(psi_container, devices, model, S, get_feedforward(model))
    time_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S, get_feedforward(model))

    return
end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    commitment_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, D)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    commitment_constraints!(psi_container, devices, model, S, get_feedforward(model))
    ramp_constraints!(psi_container, devices, model, S, get_feedforward(model))
    time_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S, get_feedforward(model))

    return
end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, ThermalBasicUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    reactivepower_variables!(psi_container, devices)
    commitment_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactivepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    commitment_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, model.formulation, S, get_feedforward(model))

    return
end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, ThermalBasicUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    commitment_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    commitment_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, model.formulation, S, get_feedforward(model))

    return
end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, ThermalRampLimited},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    reactivepower_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactivepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    ramp_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, model.formulation, S, get_feedforward(model))

    return
end

"""
This function creates the model for a full themal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, ThermalRampLimited},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)

    #Initial Conditions
    initial_conditions!(psi_container, devices, model.formulation)

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    ramp_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, model.formulation, S, get_feedforward(model))

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)
    reactivepower_variables!(psi_container, devices)

    #Initial Conditions

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    reactivepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S, get_feedforward(model))

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    #Variables
    activepower_variables!(psi_container, devices)

    #Initial Conditions

    #Constraints
    activepower_constraints!(psi_container, devices, model, S, get_feedforward(model))
    feedforward!(psi_container, T, get_feedforward(model))

    #Cost Function
    cost_function(psi_container, devices, D, S, get_feedforward(model))

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{T, FixedOutput},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    if !validate_available_devices(T, devices)
        return
    end

    nodal_expression!(psi_container, devices, S)

    return
end

ext: jl


text:
# dummy module for testing assignment in parcel
module A
    module B
        module C
        struct CT
            x::Int
        end
        end
        module D
        end
    end
    f(a) = 1
end


text:
using Test
using IncludeGuards

@testset "Working tests" begin

	@includeonce "TestModuleWorking.jl"

	@test true

	@debug "runtests: GUARDS = " IncludeGuards.GUARDS
	println()

	@includeonce "TestModuleWorking.jl"

end

max_stars_repo_name: tgflynn/IncludeGuards

max_issues_repo_name: tgflynn/IncludeGuards

max_forks_repo_name: tgflynn/IncludeGuards


text:
using CompScienceMeshes, BEAST

width, height, h = 1.0, 0.5, 0.05
G1 = meshrectangle(width, height, h)
G2 = CompScienceMeshes.rotate(G1, 0.5π * x̂)
G3 = CompScienceMeshes.rotate(G1, 1.0π * x̂)

G12 = weld(G1,-G2)
G23 = weld(G2,-G3)
G31 = weld(G3,-G1)

X12 = raviartthomas(G12)
X23 = raviartthomas(G23)
X31 = raviartthomas(G31)

Y12 = buffachristiansen(G12)
Y23 = buffachristiansen(G23)
Y31 = buffachristiansen(G31)

κ = 1.0
SL = Maxwell3D.singlelayer(wavenumber=κ)
N = NCross()

X = X12 × X23 × X31
Y = Y12 × Y23 × Y31

E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
e = (n × E) × n

@hilbertspace j
@hilbertspace k
mtefie = @discretise(
    SL[k,j] == e[k],
    j ∈ X, k ∈ X)

# u = gmres(mtefie)
#
# offset = 1; u12 = u[offset:offset+numfunctions(X12)-1]
# offset += numfunctions(X12); u23 = u[offset:offset+numfunctions(X23)-1]
# offset += numfunctions(X23); u31 = u[offset:offset+numfunctions(X31)-1]
#
# fcr12, _ = facecurrents(u12, X12)
# fcr23, _ = facecurrents(u23, X23)
# fcr31, _ = facecurrents(u31, X31)

import PlotlyJS
using LinearAlgebra
# p1 = PlotlyJS.Plot(patch(G12, norm.(fcr12)))
# p2 = PlotlyJS.Plot(patch(G23, norm.(fcr23)))
# p3 = PlotlyJS.Plot(patch(G31, norm.(fcr31)))

# PlotlyJS.plot([p1, p2, p3])

function blkdiagm(blks...)
    m = sum(size(b,1) for b in blks)
    n = sum(size(b,2) for b in blks)
    A = zeros(eltype(first(blks)), m, n)
    o1 = 1
    o2 = 1
    for b in blks
        m1 = size(b,1)
        n1 = size(b,2)
        A[o1:o1+m1-1,o2:o2+n1-1] .= b
        o1 += m1
        o2 += n1
    end
    A
end

Sxx = BEAST.sysmatrix(mtefie)
N1 = assemble(N, X12, Y12)
N2 = assemble(N, X23, Y23)
N3 = assemble(N, X31, Y31)
Nxy = blkdiagm(N1,N2,N3)
Syy = assemble(SL, Y, Y)
iNxy = inv(Nxy)

cond(Sxx)
ex = assemble(e,X)
Q = transpose(iNxy) * Syy * iNxy * Sxx;
R = transpose(iNxy) * Syy * iNxy * ex;

u1, ch1 = solve(BEAST.GMRESSolver(Sxx),ex)
u2, ch2 = solve(BEAST.GMRESSolver(Q),R)

@show ch1.iters
@show ch2.iters
# gmres()


text:
"""
    DensityMatrix(state::AbstractArray{T, 3}; nlevel=2)
    DensityMatrix(state::AbstractMatrix{T}; nlevel=2)

Create a `DensityMatrix` with a state represented by array.
"""
YaoBase.DensityMatrix{D}(state::AbstractMatrix{T}) where {T,D} = DensityMatrix{D,T,typeof(state)}(state)
YaoBase.DensityMatrix(state::AbstractMatrix{T}; nlevel=2) where T = DensityMatrix{nlevel}(state)

"""
    state(ρ::DensityMatrix)

Return the raw state of density matrix `ρ`.
"""
state(ρ::DensityMatrix) = ρ.state

YaoBase.nqubits(ρ::DensityMatrix) = nqudits(ρ)
YaoBase.nqudits(ρ::DensityMatrix{D}) where {D} = logdi(size(state(ρ), 1), D)
YaoBase.nactive(ρ::DensityMatrix) = nqudits(ρ)

"""
    density_matrix(reg, qubits)

Get the reduced density matrix on given `locs`. See also [`focus!`](@ref).
"""
function YaoBase.density_matrix(reg::ArrayReg, qubits)
    freg = focus!(copy(reg), qubits)
    return density_matrix(freg)
end
YaoBase.density_matrix(reg::ArrayReg{D}) where D = DensityMatrix{D}(reg.state * reg.state')
YaoBase.tracedist(dm1::DensityMatrix{D}, dm2::DensityMatrix{D}) where {D} = trnorm(dm1.state .- dm2.state)

# TODO: use batch_broadcast in the future
"""
    probs(ρ)

Returns the probability distribution from a density matrix `ρ`.
"""
YaoBase.probs(m::DensityMatrix) = diag(m.state)

function YaoBase.purify(r::DensityMatrix{D}; num_env::Int = nactive(r)) where {D}
    Ne = D ^ num_env
    Ns = size(r.state, 1)
    R, U = eigen!(r.state)
    state = view(U, :, Ns-Ne+1:Ns) .* sqrt.(abs.(view(R, Ns-Ne+1:Ns)'))
    return ArrayReg{D}(state)
end

# obtaining matrix from Yao.DensityMatrix
LinearAlgebra.Matrix(d::DensityMatrix) = d.state

von_neumann_entropy(dm::DensityMatrix) = von_neumann_entropy(Matrix(dm))
function von_neumann_entropy(dm::AbstractMatrix)
    p = max.(eigvals(dm), eps(real(eltype(dm))))
    return von_neumann_entropy(p)
end
von_neumann_entropy(v::AbstractVector) = -sum(x->x*log(x), v)


text:
module CumulantsFeatures
  using Distributions
  using SymmetricTensors
  using Cumulants
  using Distributed
  using StatsBase
  using LinearAlgebra
  using SharedArrays
  using CumulantsUpdates
  #using DistributedArrays
  if VERSION >= v"1.3"
    using CompilerSupportLibraries_jll
  end

  import SymmetricTensors: issymetric, getblock, pyramidindices

  include("optimizationalgorithms.jl")
  include("detectors.jl")
  include("symten2mat.jl")

  export cumfsel, rxdetect, hosvdc4detect, cum2mat
end


text:
#!/usr/bin/env julia

export pe4

function pe4(st::Integer=100,ed::Integer=999)
    ispalindrome(n)=digits(n)==reverse(digits(n))
    maximum(i*j for i in st:ed, j in st:ed if ispalindrome(i*j))
end

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe4()|>print
end

max_stars_repo_path: 4. Largest palindrome product/4.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 4. Largest palindrome product/4.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 4. Largest palindrome product/4.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
"""
```
computeDistPars( attrition::Attrition )
```
This function pre-computes some reoccurring parameters of the distribution of
the time to attrition of the attrition scheme `attrition`.

This function returns `nothing`.
"""
function computeDistPars( attrition::Attrition )

    attrition.lambdas = - log.( 1 .- attrition.rates ) / attrition.period
    gammas = attrition.curvePoints[2:end] -
        attrition.curvePoints[1:(end - 1)]
    gammas .*= - attrition.lambdas[1:(end - 1 )]
    gammas = exp.( gammas )
    attrition.gammas = vcat( 1.0, cumprod( gammas ) )
    return

end  # computeDistPars( attrition )


text:
#= definitions.jl
Philip Barrett, pobarrett@gmail.com
30may2016, Chicago

Tests the polygon definitions
=#

### NEED TO RE-CHECK ALL THESE TESTS

Z = [ -1 -1; 1.0 -1; 1 1; -1 1 ]::Matrix{Float64}
G = [ 0 -1; 1.0 0; 0 1; -1 0 ]
m = [ 1.0, 1, 1, 1 ]
a = Polygon( Z, G, m ) ;
b = Polygon( pts=Z ) ;
c = Polygon( dirs=G, dists=m ) ;

badZ = Z[ [1, 3, 4, 2], : ]
    # Mis-ordered Z
d = Polygon(pts=badZ)

intZ = vcat( Z, [0.0 0.5 ] )
    # Interior point added
e = Polygon(pts=intZ)
    # Works well enough

dupeZ = Z[ [1, 1, 2, 3, 4, 1, 2 ], : ]
    # Duplicate points
f = Polygon(pts=dupeZ)
    # Also good

nearZ = vcat( Z, Z[ [3,3,3,3,3],:] +
          [0 -1e-11 ; -1e-11 0 ; -1e-11 -1e-11 ; 0 -1e-11 ; -1e-11 0 ] )
g = Polygon(pts=nearZ)
    # Works nicely

for( compare in [b c d e f g] )
    @test a.pts == compare.pts
    @test a.dirs == compare.dirs
    @test a.dists == compare.dists
end


text:
using CompScienceMeshes
using BEAST

Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
X = raviartthomas(Γ)
Y = buffachristiansen(Γ, ibscaled=true)

N = BEAST.NCross()

κ = 1.0
E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
e = (n × E) × n

b = assemble(e, X)
G = assemble(N, X, Y)

x = G \ b
fcr, geo = facecurrents(x, Y)

using Plotly
using LinearAlgebra
Plotly.plot(patch(geo.mesh, norm.(fcr), (0.0, 1.0)))


text: push!(LOAD_PATH, "./jl/")

max_stars_repo_name: wilsonkl/sfmflex-release

max_issues_repo_name: wilsonkl/sfmflex-release

max_forks_repo_name: wilsonkl/sfmflex-release


text:
# SimLynx/src/rendezvous.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

"""
    Message <: AbstractMessage

Represents a message being send to a process.

# Fields
- `time::Float64` - the time the message was sent
- `caller::Process` - the process sending the message
- `name::Symbol` - the type of message being send
- `args::Array{Any, 1}` - the arguments for the message
"""
struct Message <: AbstractMessage
    time::Float64
    caller::Process
    name::Symbol
    args::Array{Any, 1}
end

"""
Check for interprocess communications rendezvous. Returns the message and
acceptor and a Boolean indicating whether a rendezvous will occur. The first
two returned values are nothing if the Boolean is false.
"""
function rendezvous(process::Process)
    for (i, message) in enumerate(process.queue)
        for acceptor in process.acceptors
            if acceptor.name == message.name
                deleteat!(process.queue, i)
                return (message, acceptor, true)
            end
        end
    end
    return (nothing, nothing, false)
end

"""
    accept <caller> <signature>

Accept a message from `caller` with the given `signature`.
"""
macro accept(caller::Symbol, sig)
    @capture(sig, f_Symbol()) ||
        throw(ArgumentError("second argument must be a signature with no arguments, " *
                            "given $sig"))
    acceptors = [Acceptor(f, nothing)]
    quote
        current_process().acceptors = $acceptors
        while true
            message, acceptor, accepted = rendezvous(current_process())
            if accepted
                message.caller.response = nothing
                @schedule immediate message.caller
                break
            else
                yieldto(control_task())
            end
        end
    end
end

"""
    accept <caller> <signature> begin
        <body>
    end

Accept a message from `caller` with the given `signature`.

The `caller` will be bound to the process instance that sent the message.

The `signature` is <name>(<arg>...) where `name` is the name associated with
the message and the `arg`s are bound to the corresponding values in the
message.
"""
macro accept(caller::Symbol, sig, body)
    @capture(sig, f_Symbol(xs__)) ||
        throw(ArgumentError("second argument must be a signature, " *
                            "given $sig"))
    args = Expr(:tuple, caller, xs...)
    proc = Expr(:->, args, body)
    acceptors = [Acceptor(f, eval(proc))]
    quote
        current_process().acceptors = $acceptors
        while true
            message, acceptor, accepted = rendezvous(current_process())
            if accepted
                if !isnothing(acceptor.proc)
                    message.caller.response = acceptor.proc(message.args...)
                else
                    message.caller.response = nothing
                end
                @schedule immediate message.caller
                break
            else
                yieldto(control_task())
            end
        end
    end
end

macro send(callee, sig)
    if !(isa(callee, Symbol) ||
         isa(callee, Expr))
        throw(ArgumentError("first argument must be a symbol or expression, " *
                            "given $callee"))
    end
    if !(isa(sig, Expr) &&
             sig.head === :call &&
             length(sig.args) > 0 &&
             isa(sig.args[1], Symbol))
        throw(ArgumentError("second argument must be a signature, " *
                            "given $sig"))
    end
    name = QuoteNode(sig.args[1])
    quote
        if isnothing(current_process())
            Error("call to send must be within a process")
        end
        let callee_process = $(esc(callee)),
            message = Message(current_time(),
                              current_process(),
                              $name,
                              [current_process(), $(esc.(sig.args[2:end])...)])
            push!(callee_process.queue, message)
            process_state!(current_process(), :delayed)
            @schedule immediate callee_process
            yieldto(control_task())
            current_process().response
        end
    end
end


text:
module MomentDistances

export distance, NamedSum, Weighted, ElementwiseMean, AbsoluteRelative, summarize

using ArgCheck: @argcheck
using DocStringExtensions: FUNCTIONNAME, SIGNATURES, TYPEDEF
using LinearAlgebra: LinearAlgebra
using Base.Multimedia: MIME, @MIME_str
using Statistics: mean
using UnPack: @unpack

####
#### core functionality
####

"""
`$(FUNCTIONNAME)(metric, x, y)`

Calculate the distance (a real number) between `x` and `y` using `metric`.
"""
function distance end

struct NamedSum{T <: NamedTuple}
    named_metrics::T
end

"""
$(SIGNATURES)

Helper function for named sums.
"""
function _named_distance_sum(named_metrics::NamedTuple{K}, x, y) where K
    if @generated
        mapreduce(k -> :(distance(named_metrics.$(k), x.$(k), y.$(k))),
                  (a, b) -> :($(a) + $(b)), K)
    else
        mapreduce((k, v) -> distance(v, getproperty(x, k), getproperty(y, k)),
                  +, pairs(named_metrics))
    end
end

distance(metric::NamedSum, x, y) = _named_distance_sum(metric.named_metrics, x, y)

struct Weighted{M,T <: Real}
    metric::M
    weight::T
    function Weighted(metric::M, weight::T) where {M, T <: Real}
        @argcheck weight > 0
        new{M,T}(metric, weight)
    end
end

function distance(metric::Weighted, x, y)
    metric.weight * distance(metric.metric, x, y)
end

struct ElementwiseMean{M}
    elementwise_metric::M
end

function distance(metric::ElementwiseMean, x, y)
    @argcheck axes(x) == axes(y) DimensionMismatch
    @unpack elementwise_metric = metric
    mean(((x, y),) -> distance(elementwise_metric, x, y), zip(x, y))
end

Base.@kwdef struct AbsoluteRelative{T <: Union{Nothing,Real},F}
    relative_adjustment::T = nothing
    norm::F = LinearAlgebra.norm
end

function distance(metric::AbsoluteRelative, x, y)
    @unpack relative_adjustment, norm = metric
    Δ = norm(x - y)
    if relative_adjustment ≡ nothing
        Δ
    else
        A = max(norm(x), norm(y))
        Δ / (relative_adjustment == Inf ? A : max(1, relative_adjustment * A))
    end
end

####
#### summaries
####
#### FIXME the whole code below is probably a weird abuse of the multimedia I/O system,
#### revamp when summarizing for other outputs. `summary` should be fully internal, an
#### implementation detail constructing a summary, `show` methods should be defined for a
#### wrapper type?

"""
    summary([options], [mime], metric, x, y)

Return summary of how a metric was calculated, according to the specific MIME type.
"""
summary(metric, x, y) = summary(MIME("text/plain"), metric, x, y)

###
### plain text summaries
###

function summary(mime::MIME"text/plain", metric, x, y)
    summary(TextSummaryOptions(), mime, metric, x, y)
end

"Options for printing text summaries."
Base.@kwdef struct TextSummaryOptions
    sigdigits::Int = 3
end

function _dotted_repr(options, x)
    lines = split(repr(x), '\n')
    line1 = first(lines)
    length(lines) > 1 ? line1 * "…" : line1
end

function _dotted_repr(options::TextSummaryOptions, x::Real)
    repr(round(x; sigdigits = options.sigdigits))
end

function _indent(str, n = 2)
    join(map(line -> ' '^n * line, split(str, '\n')), '\n')
end

function summary(options, ::MIME"text/plain", metric, x, y)
    # this is the fallback, and also the indended summary for scalar metrics like
    # AbsoluteRelative
    x_, y_, d_ = _dotted_repr.(Ref(options), (x, y, distance(metric, x, y)))
    "‹" * x_ * " ↔ " * y_ * ": " * d_ * "›"
end

function summary(options, mime::MIME"text/plain", metric::ElementwiseMean, x, y)
    @argcheck axes(x) == axes(y)
    header = "elementwise mean distance: " * _dotted_repr(options, distance(metric, x, y))
    digits_by_axis = ntuple(ndims(x)) do i
        max(length(string(firstindex(x, i))), length(string(lastindex(x, i))))
    end
    body = mapreduce(*, CartesianIndices(x), x, y) do i, x, y
        padded_index = mapreduce((d, i) -> lpad(string(i), d, ' '), (a, b) -> a * "," * b,
                                 digits_by_axis, Tuple(i))
        "\n" * _indent("[" * padded_index * "]  " *
                       summary(options, mime, metric.elementwise_metric, x, y))
    end
    header * body
end

function summary(options, mime::MIME"text/plain", metric::Weighted, x, y)
    "weighted: " * _dotted_repr(options, distance(metric, x, y)) * '\n' *
        _indent(summary(options, mime, metric.metric, x, y))
end

function summary(options, mime::MIME"text/plain", metric::NamedSum, x, y)
    str = "total: " * _dotted_repr(options, distance(metric, x, y))
    for (key, metric) in pairs(metric.named_metrics)
        str *= "\n" * _indent("from $(key):\n" *
                              _indent(summary(options, mime, metric,
                                              getproperty(x, key), getproperty(y, key))))
    end
    str
end

function summarize(io, mime::MIME"text/plain", metric, x, y)
    print(io, summary(mime, metric, x, y))
end

summarize(io, metric, x, y) = summarize(io, MIME("text/plain"), metric, x, y)

summarize(metric, x, y) = summarize(stdout, metric, x, y)

end # module


text:
export get_hessian

function get_hessian(params0::AbstractArray{T}) where T
    N = length(params0)
    params = Dual.(params0, zero(T))
    hes = zeros(T, N, N)
    for i=1:N
        @inbounds i !== 1 && (params[i-1] = Dual(params0[i-1], zero(T)))
        @inbounds params[i] = Dual(params0[i], one(T))
        res = NiLang.AD.gradient(Val(1), embedding_loss, (Dual(0.0, 0.0), params))[2]
        hes[:,i] .= vec(ForwardDiff.partials.(res, 1))
    end
    hes
end


text:
function jt07a()
#
#  test polynomial suggested by Jenkins and Traub
#
   a = 10.0^(-10);
   z = [.001; .01; .1; .1+a*im; .1-a*im; 1; -10]
   p = poly(z)
   z = [z ones(7)]
   y = [-10  1; 1  1; 0.01  1; 0.001  1; 0.1  3];
   p, PolyZeros(z), PolyZeros(y)
end


text:
using LinearAlgebra
using DataFrames
using Test
println("LATHE.PREPROCESS TESTS")
println("......................")
# DECOMPOSITION
using Lathe.preprocess: SVD
@testset "svd test for $T" for T in (Float32, Float64)
  A = randn(T, 100, 100)
  x = SVD(A)
  tol = T == Float32 ? 1e-5 : 1e-10
  @test (norm(A - (x.U * diagm(0 => x.S) * x.V'), 2) / norm(A, 2)) <= tol
end
using Lathe.preprocess: OneHotEncoder, FloatEncoder, OrdinalEncoder
# ENCODING
df = DataFrame(:A => [5, 10, 15, 20], :B => ["test", "test2", "test2", "test"])
@testset "OneHotEncoder Tests" begin
encoder = OneHotEncoder()
ant = size(df)[2]
data = encoder.predict(df, :A)
len = Set(df[!, :A])
@test size(data)[2] - ant == length(len)

end
@testset "OrdinalEncoder Tests" begin
encoder = OrdinalEncoder(df[!, :B])
encoded = encoder.predict(df[!, :B])
inverse = [z => x for (x, z) in encoder.lookup]
encoded = [inverse[i] for i in encoded]
for (i, c) in zip(encoded, df[!, :B])
  @test i == c
end
end
@testset "FloatEncoder Tests" begin
  encoder = FloatEncoder()
  FloatEncoder.predict(df[!, :B])
  @test typeof(df[!, :B][1]) == Float64
end


text:
# This is a test file to make codecov happy, technically all of the
# tests here are already done / integrated within other tests.

@testset "strings" begin
    st = "blah"

    @test F.str(st) == "blah"

    sst = SubString("blahblah", 1:4)
    @test sst == "blah"
    @test F.str(sst) == "blahblah"

    sst = SubString("blah✅💕and etcσ⭒ but ∃⫙∀ done", 1:27)
    @test F.to(sst) == 27

    s = "aabccabcdefabcg"
    for m ∈ eachmatch(r"abc", s)
        @test s[F.matchrange(m)] == "abc"
    end
end

@testset "ocblock" begin
    st = "This is a block  and done"
    τ = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    ocb = F.OCBlock(:COMMENT, (τ[1]=>τ[2]))
    @test F.otok(ocb) == τ[1]
    @test F.ctok(ocb) == τ[2]
end

@testset "isexactly" begin
    steps, b, λ = F.isexactly("


text:
using OnlinePCA
using DelimitedFiles: writedlm
writecsv(filename::AbstractString, data) = writedlm(filename, data, ',')

# Setting
chunksize = parse(Int64, ARGS[1])
outdir = ARGS[2]

# Randomized SVD
out = OnlinePCA.tenxpca(
	tenxfile="../Data/Brain/1M_neurons_filtered_gene_bc_matrices_h5.h5",
	rowmeanlist="../Data/Brain/Feature_LogMeans.csv",
	scale="log",
	dim=20,
	niter=3,
	chunksize=chunksize,
	group="mm10")

# Output
writecsv(joinpath(outdir, "Eigen_vectors.csv"), out[1])
writecsv(joinpath(outdir, "Eigen_values.csv"), out[2])
writecsv(joinpath(outdir, "Loadings.csv"), out[3])

max_stars_repo_name: rikenbit/onlinePCA-experiments

max_issues_repo_name: rikenbit/onlinePCA-experiments

max_forks_repo_name: rikenbit/onlinePCA-experiments


text:
module SIMeasurementsTest
using SIMeasurements
using FactCheck: @fact, facts, context, exitstatus, roughly, exactly,
                 @fact_throws, greater_than, not

contains(x) = y -> x ∈ y
all_integers(x::Array, ε::AbstractFloat=1e-8) = all(abs(x - round(Integer, x)) .< ε)
all_integers(ε::AbstractFloat=1e-8) = y -> all_integers(y, ε)
is_subtype(x::Type) = y -> (y <: x)

include("units.jl")
include("quantity.jl")

exitstatus()
end


text:
function correct_grid_size_regular(FT)
    grid = RegularCartesianGrid(FT; size=(4, 6, 8), length=(2π, 4π, 9π))

    # Checking ≈ as the grid could be storing Float32 values.
    return (grid.Nx ≈ 4  && grid.Ny ≈ 6  && grid.Nz ≈ 8 &&
            grid.Lx ≈ 2π && grid.Ly ≈ 4π && grid.Lz ≈ 9π)
end

function correct_halo_size_regular(FT)
    grid = RegularCartesianGrid(FT; size=(4, 6, 8), length=(2π, 4π, 9π), halo=(1, 2, 3))
    return (grid.Hx == 1  && grid.Hy == 2  && grid.Hz == 3)
end

function faces_start_at_zero_regular(FT)
    grid = RegularCartesianGrid(FT; size=(10, 10, 10), length=(2π, 2π, 2π))
    return grid.xF[1] == 0 && grid.yF[1] == 0 && grid.zF[end] == 0
end

function end_faces_match_grid_length_regular(FT)
    grid = RegularCartesianGrid(FT; size=(12, 13, 14), length=(π, π^2, π^3))
    return (grid.xF[end] - grid.xF[1] ≈ π   &&
            grid.yF[end] - grid.yF[1] ≈ π^2 &&
            grid.zF[end] - grid.zF[1] ≈ π^3)
end

function ranges_have_correct_length_regular(FT)
    Nx, Ny, Nz = 8, 9, 10
    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1, 1, 1))
    return (length(grid.xC) == Nx && length(grid.xF) == Nx+1 &&
            length(grid.yC) == Ny && length(grid.yF) == Ny+1 &&
            length(grid.zC) == Nz && length(grid.zF) == Nz+1)
end

# See: https://github.com/climate-machine/Oceananigans.jl/issues/480
function no_roundoff_error_in_ranges_regular(FT)
    Nx, Ny, Nz = 1, 1, 64
    grid = RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(1, 1, π/2))
    return length(grid.zC) == Nz && length(grid.zF) == Nz+1
end

function grid_properties_are_same_type_regular(FT)
    grid = RegularCartesianGrid(FT; size=(10, 10, 10), length=(1, 1//7, 2π))
    return all(isa.([grid.Lx, grid.Ly, grid.Lz, grid.Δx, grid.Δy, grid.Δz], FT)) &&
           all(eltype.([grid.xF, grid.yF, grid.zF, grid.xC, grid.yC, grid.zC]) .== FT)
end

function correct_constant_grid_spacings(FT)
    grid = VerticallyStretchedCartesianGrid(FT; size=(16, 16, 16), x=(0,1), y=(0,1), z=(0,16), zF=collect(0:16))
    return all(grid.ΔzF .== 1) && all(grid.ΔzC .== 1)
end

function correct_quadratic_grid_spacings(FT)
    Nx = Ny = Nz = 16
    grid = VerticallyStretchedCartesianGrid(FT; size=(Nx, Ny, Nz), x=(0,1), y=(0,1), z=(0,Nz^2), zF=collect(0:Nz).^2)

     zF(k) = (k-1)^2
     zC(k) = ((k-1)^2 + k^2) / 2
    ΔzF(k) = k^2 - (k-1)^2
    ΔzC(k) = 2k

     zF_is_correct = all(isapprox.(grid.zF,   zF.(1:Nz+1)))
     zC_is_correct = all(isapprox.(grid.zC,   zC.(1:Nz)))
    ΔzF_is_correct = all(isapprox.(grid.ΔzF, ΔzF.(1:Nz)))
    ΔzC_is_correct = all(isapprox.(grid.ΔzC, ΔzC.(1:Nz-1)))

    return zF_is_correct && zC_is_correct && ΔzF_is_correct && ΔzC_is_correct
end

function correct_tanh_grid_spacings(FT)
    Nx = Ny = Nz = 16

    S = 3  # Stretching factor
    zF(k) = tanh(S * (2*(k-1)/Nz - 1)) / tanh(S)

    grid = VerticallyStretchedCartesianGrid(FT; size=(Nx, Ny, Nz), x=(0,1), y=(0,1), z=(-1,1), zF=zF)

     zC(k) = (zF(k) + zF(k+1)) / 2
    ΔzF(k) = zF(k+1) - zF(k)
    ΔzC(k) = zC(k+1) - zC(k)

    zF_is_correct = all(isapprox.(grid.zF,   zF.(1:Nz+1)))
    zC_is_correct = all(isapprox.(grid.zC,   zC.(1:Nz)))
   ΔzF_is_correct = all(isapprox.(grid.ΔzF, ΔzF.(1:Nz)))
   ΔzC_is_correct = all(isapprox.(grid.ΔzC, ΔzC.(1:Nz-1)))

   return zF_is_correct && zC_is_correct && ΔzF_is_correct && ΔzC_is_correct
end

function grid_properties_are_same_type_stretched(FT)
    Nx, Ny, Nz = 16, 16, 16
    grid = VerticallyStretchedCartesianGrid(FT; size=(16, 16, 16), x=(0,1), y=(0,1), z=(0.0,16.0), zF=collect(0:16))
    return all(isa.([grid.Lx, grid.Ly, grid.Lz, grid.Δx, grid.Δy], FT)) &&
           all(eltype.([grid.ΔzF, grid.ΔzC, grid.xF, grid.yF, grid.zF, grid.xC, grid.yC, grid.zC]) .== FT)
end

@testset "Grids" begin
    @info "Testing grids..."

    @testset "Regular Cartesian grid" begin
        @info "  Testing regular Cartesian grid..."

        @testset "Grid initialization" begin
            @info "    Testing grid initialization..."

            for FT in float_types
                @test correct_grid_size_regular(FT)
                @test correct_halo_size_regular(FT)
                @test faces_start_at_zero_regular(FT)
                @test end_faces_match_grid_length_regular(FT)
                @test ranges_have_correct_length_regular(FT)
                @test no_roundoff_error_in_ranges_regular(FT)
                @test grid_properties_are_same_type_regular(FT)
            end
        end

        @testset "Grid dimensions" begin
            @info "    Testing grid constructor errors..."

            for FT in float_types
                @test isbitstype(typeof(RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 1, 1))))

                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32,), length=(1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 64), length=(1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32, 16), length=(1, 1, 1))

                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32.0), length=(1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(20.1, 32, 32), length=(1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, nothing, 32), length=(1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, "32", 32), length=(1, 1, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, nothing, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, "1", 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, 1, 1), halo=(1, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(32, 32, 32), length=(1, 1, 1), halo=(1.0, 1, 1))

                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=2)
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), y=[1, 2])
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), z=(-π, π))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=1, y=2, z=3)
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(0, 2), z=4)
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=("0", "1"))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(-1//2, 1), y=(1//7, 5//7), z=(1, 2, 3))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(1, 0), y=(1//7, 5//7), z=(1, 2))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), x=(0, 1), y=(1, 5), z=(π, -π))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 2, 3), x=(0, 1))
                @test_throws ArgumentError RegularCartesianGrid(FT; size=(16, 16, 16), length=(1, 2, 3), x=(0, 1), y=(1, 5), z=(-π, π))
            end
        end
    end

    @testset "Vertically stretched Cartesian grid" begin
        @info "  Testing vertically stretched Cartesian grid..."

        @testset "Grid initialization" begin
            @info "    Testing grid initialization..."

            for FT in float_types
                @test correct_constant_grid_spacings(FT)
                @test correct_quadratic_grid_spacings(FT)
                @test correct_tanh_grid_spacings(FT)
                @test grid_properties_are_same_type_stretched(FT)
            end
        end
    end
end


text:
export lmnn

using mlpack._Internal.io

import mlpack_jll
const lmnnLibrary = mlpack_jll.libmlpack_julia_lmnn

# Call the C binding of the mlpack lmnn binding.
function lmnn_mlpackMain()
  success = ccall((:lmnn, lmnnLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module lmnn_internal
  import ..lmnnLibrary

end # module

"""
    lmnn(input; [batch_size, center, distance, k, labels, linear_scan, max_iterations, normalize, optimizer, passes, print_accuracy, range, rank, regularization, seed, step_size, tolerance, verbose])

This program implements Large Margin Nearest Neighbors, a distance learning
technique.  The method seeks to improve k-nearest-neighbor classification on a
dataset.  The method employes the strategy of reducing distance between similar
labeled data points (a.k.a target neighbors) and increasing distance between
differently labeled points (a.k.a impostors) using standard optimization
techniques over the gradient of the distance between data points.

To work, this algorithm needs labeled data.  It can be given as the last row of
the input dataset (specified with `input`), or alternatively as a separate
matrix (specified with `labels`).  Additionally, a starting point for
optimization (specified with `distance`can be given, having (r x d)
dimensionality.  Here r should satisfy 1 <= r <= d, Consequently a Low-Rank
matrix will be optimized. Alternatively, Low-Rank distance can be learned by
specifying the `rank`parameter (A Low-Rank matrix with uniformly distributed
values will be used as initial learning point).

The program also requires number of targets neighbors to work with ( specified
with `k`), A regularization parameter can also be passed, It acts as a trade of
between the pulling and pushing terms (specified with `regularization`), In
addition, this implementation of LMNN includes a parameter to decide the
interval after which impostors must be re-calculated (specified with `range`).

Output can either be the learned distance matrix (specified with `output`), or
the transformed dataset  (specified with `transformed_data`), or both.
Additionally mean-centered dataset (specified with `centered_data`) can be
accessed given mean-centering (specified with `center`) is performed on the
dataset. Accuracy on initial dataset and final transformed dataset can be
printed by specifying the `print_accuracy`parameter.

This implementation of LMNN uses AdaGrad, BigBatch_SGD, stochastic gradient
descent, mini-batch stochastic gradient descent, or the L_BFGS optimizer.

AdaGrad, specified by the value 'adagrad' for the parameter `optimizer`, uses
maximum of past squared gradients. It primarily on six parameters: the step size
(specified with `step_size`), the batch size (specified with `batch_size`), the
maximum number of passes (specified with `passes`). Inaddition, a normalized
starting point can be used by specifying the `normalize` parameter.

BigBatch_SGD, specified by the value 'bbsgd' for the parameter `optimizer`,
depends primarily on four parameters: the step size (specified with
`step_size`), the batch size (specified with `batch_size`), the maximum number
of passes (specified with `passes`).  In addition, a normalized starting point
can be used by specifying the `normalize` parameter.

Stochastic gradient descent, specified by the value 'sgd' for the parameter
`optimizer`, depends primarily on three parameters: the step size (specified
with `step_size`), the batch size (specified with `batch_size`), and the maximum
number of passes (specified with `passes`).  In addition, a normalized starting
point can be used by specifying the `normalize` parameter. Furthermore,
mean-centering can be performed on the dataset by specifying the
`center`parameter.

The L-BFGS optimizer, specified by the value 'lbfgs' for the parameter
`optimizer`, uses a back-tracking line search algorithm to minimize a function.
The following parameters are used by L-BFGS: `max_iterations`, `tolerance`(the
optimization is terminated when the gradient norm is below this value).  For
more details on the L-BFGS optimizer, consult either the mlpack L-BFGS
documentation (in lbfgs.hpp) or the vast set of published literature on L-BFGS.
In addition, a normalized starting point can be used by specifying the
`normalize` parameter.

By default, the AMSGrad optimizer is used.

Example - Let's say we want to learn distance on iris dataset with number of
targets as 3 using BigBatch_SGD optimizer. A simple call for the same will look
like:

```julia
julia> using CSV
julia> iris = CSV.read("iris.csv")
julia> iris_labels = CSV.read("iris_labels.csv"; type=Int)
julia> _, output, _ = mlpack_lmnn(iris; k=3, labels=iris_labels,
            optimizer="bbsgd")
```

An another program call making use of range & regularization parameter with
dataset having labels as last column can be made as:

```julia
julia> using CSV
julia> letter_recognition = CSV.read("letter_recognition.csv")
julia> _, output, _ = mlpack_lmnn(letter_recognition; k=5, range=10,
            regularization=0.4)
```

# Arguments

 - `input::Array{Float64, 2}`: Input dataset to run LMNN on.
 - `batch_size::Int`: Batch size for mini-batch SGD.  Default value `50`.

 - `center::Bool`: Perform mean-centering on the dataset. It is useful
      when the centroid of the data is far from the origin.  Default value
      `false`.

 - `distance::Array{Float64, 2}`: Initial distance matrix to be used as
      starting point
 - `k::Int`: Number of target neighbors to use for each datapoint.
      Default value `1`.

 - `labels::Array{Int, 1}`: Labels for input dataset.
 - `linear_scan::Bool`: Don't shuffle the order in which data points are
      visited for SGD or mini-batch SGD.  Default value `false`.

 - `max_iterations::Int`: Maximum number of iterations for L-BFGS (0
      indicates no limit).  Default value `100000`.

 - `normalize::Bool`: Use a normalized starting point for optimization.
      Itis useful for when points are far apart, or when SGD is returning NaN.
      Default value `false`.

 - `optimizer::String`: Optimizer to use; 'amsgrad', 'bbsgd', 'sgd', or
      'lbfgs'.  Default value `amsgrad`.

 - `passes::Int`: Maximum number of full passes over dataset for AMSGrad,
      BB_SGD and SGD.  Default value `50`.

 - `print_accuracy::Bool`: Print accuracies on initial and transformed
      dataset  Default value `false`.

 - `range::Int`: Number of iterations after which impostors needs to be
      recalculated  Default value `1`.

 - `rank::Int`: Rank of distance matrix to be optimized.   Default value
      `0`.

 - `regularization::Float64`: Regularization for LMNN objective function
      Default value `0.5`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `step_size::Float64`: Step size for AMSGrad, BB_SGD and SGD (alpha).
      Default value `0.01`.

 - `tolerance::Float64`: Maximum tolerance for termination of AMSGrad,
      BB_SGD, SGD or L-BFGS.  Default value `1e-07`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `centered_data::Array{Float64, 2}`: Output matrix for mean-centered
      dataset.
 - `output::Array{Float64, 2}`: Output matrix for learned distance
      matrix.
 - `transformed_data::Array{Float64, 2}`: Output matrix for transformed
      dataset.

"""
function lmnn(input;
              batch_size::Union{Int, Missing} = missing,
              center::Union{Bool, Missing} = missing,
              distance = missing,
              k::Union{Int, Missing} = missing,
              labels = missing,
              linear_scan::Union{Bool, Missing} = missing,
              max_iterations::Union{Int, Missing} = missing,
              normalize::Union{Bool, Missing} = missing,
              optimizer::Union{String, Missing} = missing,
              passes::Union{Int, Missing} = missing,
              print_accuracy::Union{Bool, Missing} = missing,
              range::Union{Int, Missing} = missing,
              rank::Union{Int, Missing} = missing,
              regularization::Union{Float64, Missing} = missing,
              seed::Union{Int, Missing} = missing,
              step_size::Union{Float64, Missing} = missing,
              tolerance::Union{Float64, Missing} = missing,
              verbose::Union{Bool, Missing} = missing,
              points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, lmnnLibrary), Nothing, ());

  IORestoreSettings("Large Margin Nearest Neighbors (LMNN)")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(batch_size)
    IOSetParam("batch_size", convert(Int, batch_size))
  end
  if !ismissing(center)
    IOSetParam("center", convert(Bool, center))
  end
  if !ismissing(distance)
    IOSetParamMat("distance", distance, points_are_rows)
  end
  if !ismissing(k)
    IOSetParam("k", convert(Int, k))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(linear_scan)
    IOSetParam("linear_scan", convert(Bool, linear_scan))
  end
  if !ismissing(max_iterations)
    IOSetParam("max_iterations", convert(Int, max_iterations))
  end
  if !ismissing(normalize)
    IOSetParam("normalize", convert(Bool, normalize))
  end
  if !ismissing(optimizer)
    IOSetParam("optimizer", convert(String, optimizer))
  end
  if !ismissing(passes)
    IOSetParam("passes", convert(Int, passes))
  end
  if !ismissing(print_accuracy)
    IOSetParam("print_accuracy", convert(Bool, print_accuracy))
  end
  if !ismissing(range)
    IOSetParam("range", convert(Int, range))
  end
  if !ismissing(rank)
    IOSetParam("rank", convert(Int, rank))
  end
  if !ismissing(regularization)
    IOSetParam("regularization", convert(Float64, regularization))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(step_size)
    IOSetParam("step_size", convert(Float64, step_size))
  end
  if !ismissing(tolerance)
    IOSetParam("tolerance", convert(Float64, tolerance))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("centered_data")
  IOSetPassed("output")
  IOSetPassed("transformed_data")
  # Call the program.
  lmnn_mlpackMain()

  return IOGetParamMat("centered_data", points_are_rows),
         IOGetParamMat("output", points_are_rows),
         IOGetParamMat("transformed_data", points_are_rows)
end


text:
using Bootstrap

@testset "recalcuate!" begin
    # SG
    res = Equipercentile(sgtabX, sgtabY)
    res.data.X.raw = collect(1:40)
    recalculate!(res.data.X)
    @test res.data.X.stats.μ == 20.5
    @test res.data.X.table.freq == map(j -> count(i -> i == j, collect(1:40)), res.data.X.table.scale)

    # NEAT
    test = copy(neattabX)
    test.rawX = collect(0:10)
    recalculate!(test)
    @test test.tableX != neattabX.tableX
    @test test.tableV == neattabX.tableV
end

@testset "Base.copy for FreqTab" begin
    res = Linear(sgtabX, sgtabY)
    test = copy(res.data)
    merge(test, (X = freqtab([1,2,3,4,5]), ))
    @test res.data.X != test.X
    test.X.raw = nothing
    @test res.data.X.raw != test.X.raw
end

@testset "Base.copy for SGEquateResult" begin
    res = Linear(sgtabX, sgtabY)
    test = copy(res)
    initialize!(test.data.X)
    @test test.data.X != res.data.X
end

@testset "Base.copy for NEATFreqTab" begin
    test = copy(neattabX)
    test.statsX = 0
    @test neattabX.statsX != test.statsX
end

@testset "Base.copy for NEATEquateResult" begin
    res = BraunHolland(neattabX, neattabY)
    test = copy(res)
    @test true # Not good!
end

@testset "Bootstrap.draw!" begin
    seed!(1234)
    res = Linear(sgtabX, sgtabY)
    o = copy(res.data)
    Bootstrap.draw!(res.data, o)
    @test res.data.X.stats != o.X.stats
    @test res.data.Y.stats != o.Y.stats
end

@testset "Bootstrap.bootstrap for SG" begin
    seed!(1234)
    res = Linear(sgtabX, sgtabY)
    bsres = bootstrap(x -> coef(Linear(x...)), res.data, BasicSampling(10))
    @test size(Bootstrap.estimate_summary(bsres)) == (2, 3)
    seed!(1234)
    res = Mean(sgtabX, sgtabY)
    bsres = bootstrap(x -> coef(Mean(x...)), res.data, BasicSampling(10))
    @test size(Bootstrap.estimate_summary(bsres)) == (2, 3)
    seed!(1234)
    res = Equipercentile(sgtabX, sgtabY)
    bsres = bootstrap(x -> coef(Equipercentile(x...)), res.data, BasicSampling(10))
    @test size(Bootstrap.estimate_summary(bsres)) == (41, 3)
end

@testset "Bootstrap.bootstrap for NEAT" begin
    for f in [Tucker, LevineObservedScore, LevineTrueScore, ChainedLinear, BraunHolland, FrequencyEstimation, ChainedEquipercentile]
        res = f(neattabX, neattabY)
        bsres = bootstrap(x -> coef(f(x...)), res.data, BasicSampling(50))
        @test true
    end
end

# Other bootstrap method

# res = Linear(sgtabX, sgtabY)
# bootstrap(x -> coef(Linear(x...)), res.data, BasicSampling(10)) # OK
# bootstrap(x -> coef(Linear(x...)), res.data, AntitheticSampling(10)) # Only for vectored data.
# bootstrap(x -> coef(Linear(x...)), res.data, BalancedSampling(10)) #
# bootstrap(x -> coef(Linear(x...)), res.data, MaximumEntropySampling(10)) #

# res = ChainedLinear(neattabX, neattabY)
# bsres = bootstrap(x -> coef(ChainedLinear(x...)), res.data, BasicSampling(50))
# bsres = bootstrap(x -> coef(BraunHolland(x...)), res.data, BasicSampling(50))


text:
# How to create your own atom
# This file is a good example of how a user can create his or her own atom
# This antidiag atom is then included in n_queens.jl to solve the n queens problem

#############################################################################
# antidiag.jl
# Returns the kth anti-diagonal (counterdiagonal, secondary diagonal, or minor diagonal) of a matrix expression
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################
import Convex.sign, Convex.monotonicity, Convex.curvature, Convex.evaluate, Convex.conic_form!
using Convex: AbstractExpr, ConstVexity, Nondecreasing, has_conic_form, cache_conic_form!, get_conic_form
export antidiag

### Diagonal
### Represents the kth diagonal of an mxn matrix as a (min(m, n) - k) x 1 vector
struct AntidiagAtom <: AbstractExpr
  head::Symbol
  id_hash::UInt64
  children::Tuple{AbstractExpr}
  size::Tuple{Int, Int}
  k::Int

  function AntidiagAtom(x::AbstractExpr, k::Int=0)
    (num_rows, num_cols) = x.size

    if k >= num_cols || k <= -num_rows
      error("Bounds error in calling diag")
    end

    children = (x, )
    return new(:antidiag, hash((children, k)), children, (minimum(x.size) - k, 1), k)
  end
end

function sign(x::AntidiagAtom)
  return sign(x.children[1])
end

# The monotonicity
function monotonicity(x::AntidiagAtom)
  return (Nondecreasing(),)
end

# If we have h(x) = f o g(x), the chain rule says h''(x) = g'(x)^T f''(g(x))g'(x) + f'(g(x))g''(x);
# this represents the first term
function curvature(x::AntidiagAtom)
  return ConstVexity()
end

function evaluate(x::AntidiagAtom)
  return diag(reverse(evaluate(x.children[1]),dims=1), x.k)
end

antidiag(x::AbstractExpr, k::Int=0) = AntidiagAtom(x, k)

# Finds the "k"-th diagonal of x as a column vector
# If k == 0, it returns the main diagonal and so on
# Let x be of size m x n and d be the diagonal
# Since x is vectorized, the way canonicalization works is:
#
# 1. We calculate the size of the diagonal (sz_diag) and the first index
# of vectorized x that will be part of d
# 2. We create the coefficient matrix for vectorized x, called coeff of size
# sz_diag x mn
# 3. We populate coeff with 1s at the correct indices
# The canonical form will then be:
# coeff * x - d = 0
function conic_form!(x::AntidiagAtom, unique_conic_forms::Convex.UniqueConicForms)
  if !has_conic_form(unique_conic_forms, x)
    (num_rows, num_cols) = x.children[1].size
    k = x.k

    if k >= 0
      start_index = k * num_rows + num_rows
      sz_diag = Base.min(num_rows, num_cols - k)
    else
      start_index = num_rows + k
      sz_diag = Base.min(num_rows + k, num_cols)
    end

    select_diag = spzeros(sz_diag, length(x.children[1]))
    for i in 1:sz_diag
      select_diag[i, start_index] = 1
      start_index += num_rows - 1
    end

    objective = conic_form!(x.children[1], unique_conic_forms)
    new_obj = select_diag * objective
    cache_conic_form!(unique_conic_forms, x, new_obj)
  end
  return get_conic_form(unique_conic_forms, x)
end


text:
struct PowerMap{T<:AbstractOrbit} <: AbstractMatrix{Int}
    cc::Vector{T}
    cache::Matrix{Int}

    function PowerMap(
        ccG::AbstractVector{T},
        exp = exponent(ccG),
    ) where {T<:AbstractOrbit}
        cache = zeros(Int, length(ccG), exp)
        id = one(first(first(ccG)))
        idx = findfirst(cc -> id ∈ cc, ccG)
        cache[:, 1] .= idx
        cache[:, 2] .= 1:length(ccG)
        return new{T}(ccG, cache)
    end
end

Base.axes(pm::PowerMap) = (Base.OneTo(length(pm.cc)), 0:size(pm.cache, 2)-1)
Base.size(pm::PowerMap) = size(pm.cache)

function Base.getindex(p::PowerMap, i::Integer, j::Integer)
    j = mod(j, size(p, 2))
    @boundscheck 1 ≤ i ≤ size(p, 1) || throw(Bounds)

    if iszero(p.cache[i, j+1])
        gj = first(p.cc[i])^j
        p.cache[i, j+1] = findfirst(cc -> gj ∈ cc, p.cc)
    end
    return p.cache[i, j+1]
end


text:
module SandboxTests
import ..Pkg # ensure we are using the correct Pkg

using Test
using UUIDs
using Pkg

using ..Utils
test_test(fn, name; kwargs...) = Pkg.test(name; test_fn=fn, kwargs...)
test_test(fn; kwargs...)       = Pkg.test(;test_fn=fn, kwargs...)

@testset "Basic `test` sandboxing" begin
    # also indirectly checks that test `compat` is obeyed
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "SandboxFallback2")
        Pkg.activate(joinpath(tmp, "SandboxFallback2"))
        test_test("Unregistered") do
            json = get(Pkg.Types.Context().env.manifest, UUID("682c06a0-de6a-54ab-a142-c8b1cf79cde6"), nothing)
            @test json !== nothing
            @test json.version == v"0.20.0"
        end
    end end
    # test dependencies should be preserved, when possible
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "Sandbox_PreserveTestDeps")
        Pkg.activate(joinpath(tmp, "Sandbox_PreserveTestDeps"))
        test_test("Foo") do
            x = get(Pkg.Types.Context().env.manifest, UUID("7876af07-990d-54b4-ab0e-23690620f79a"), nothing)
            @test x !== nothing
            @test x.version == v"0.4.0"
        end
    end end
end

@testset "Basic `build` sandbox" begin
    temp_pkg_dir() do project_path; mktempdir() do tmp
        copy_test_package(tmp, "BasicSandbox")
        Pkg.activate(joinpath(tmp, "BasicSandbox"))
        Pkg.build()
    end end
end

end # module


text:
"""
```
KrusellSmithCT{T} <: AbstractCTModel{T}
```

The `KrusellSmithCT` mutable struct defines the structure of the Krusell Smith (1998) model
in 'Income and Wealth Heterogeneity in the Macroeconomy'
originally written in MATLAB by Ben Moll and SeHyoun Ahn at
https://sehyoun.com/EXAMPLE_PHACT_KS.html

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `endogenous_states_augmented::OrderedDict{Symbol,Int}`: Maps lagged states to
  their columns in the measurement and equilibrium condition equations. These
  are added after `gensys` solves the model.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

* `pseudo_observables::OrderedDict{Symbol,Int}`: Maps each pseudo-observable to
  a row in the model's pseudo-measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"krusell_smith\", cached
  here for filepath computation.

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model units.
  DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` and `Observable` for further details.

* `pseudo_observable_mappings::OrderedDict{Symbol,PseudoObservable}`: A
  dictionary that stores names and transformations to/from model units. See
  `PseudoObservable` for further details.
"""
mutable struct KrusellSmithCT{T} <: AbstractCTModel{T}
    parameters::ParameterVector{T}                          # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                        # model steady-state values
    keys::OrderedDict{Symbol,Int}                           # human-readable names for all the model
                                                            # parameters and steady-states

    endogenous_states::OrderedDict{Symbol,Vector{Int}}      # these fields are used to create matrices in the
    exogenous_shocks::OrderedDict{Symbol,Int}               # measurement and equilibrium condition equations.
    expected_shocks::OrderedDict{Symbol,Vector{Int}}        #
    equilibrium_conditions::OrderedDict{Symbol,Vector{Int}} #
    endogenous_states_augmented::OrderedDict{Symbol,Int}    #
    observables::OrderedDict{Symbol,Int}                    #
    pseudo_observables::OrderedDict{Symbol,Int}             #

    spec::String                                            # Model specification number (eg "m990")
    subspec::String                                         # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                          # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                     # Settings/flags for testing mode
    rng::MersenneTwister                                    # Random number generator
    testing::Bool                                           # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
    pseudo_observable_mappings::OrderedDict{Symbol, PseudoObservable}
end

description(m::KrusellSmithCT) = "Julia implementation of the Krusell Smith (1998) model defined in 'Income and Wealth Heterogeneity in the Macroeconomy' and originally written in MATLAB by Ben Moll and SeHyoun Ahn at https://sehyoun.com/EXAMPLE_PHACT_KS.html: KrusellSmithCT, $(m.subspec)"

"""
`init_model_indices!(m::KrusellSmithCT)`

Arguments:
`m:: KrusellSmithCT`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::KrusellSmithCT)
    # Endogenous states
    endogenous_states = collect([:value_function, :distribution, :log_aggregate_tfp, :K, :r, :w, :output, :C, :investment])

    # Exogenous shocks
    exogenous_shocks = collect([:aggregate_tfp])

    # Expectations shocks
    expected_shocks = collect([:E_V])

    # Equilibrium conditions
    equilibrium_conditions = collect([:value_function, :distribution, :log_aggregate_tfp, :K, :r, :w, :output, :C, :investment])

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = []

    # Observables
    observables = keys(m.observable_mappings)

    # Pseudo-observables
    pseudo_observables = keys(m.pseudo_observable_mappings)

    # Assign dimensions
    m.endogenous_states[:value_function]    = collect(1:200)     # size of state space grid
    m.endogenous_states[:distribution]      = collect(201:399)   # size of grid minus one for ...
    m.endogenous_states[:log_aggregate_tfp] = collect([400])     # number of shocks to distribution
    m.endogenous_states[:K]                 = [401]
    m.endogenous_states[:r]                 = [402]
    m.endogenous_states[:w]                 = [403]
    m.endogenous_states[:output]            = [404]
    m.endogenous_states[:C]                 = [405]
    m.endogenous_states[:investment]        = [406]
    m.expected_shocks[:E_V]                 = m.endogenous_states[:value_function]

    for (i,k) in enumerate(observables);             m.observables[k]             = i end
    for (i,k) in enumerate(exogenous_shocks);        m.exogenous_shocks[k]        = i end
    for (i,k) in enumerate(pseudo_observables);      m.pseudo_observables[k]      = i end
end

function KrusellSmithCT(subspec::String="ss0";
                       custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                       testing = false)

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = KrusellSmithCT{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(),
            OrderedDict{Symbol,Int}(),

            # model indices
            OrderedDict{Symbol,Vector{Int}}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Vector{Int}}(), OrderedDict{Symbol,Vector{Int}}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}(),
            OrderedDict{Symbol,PseudoObservable}())

    # Set parameters
    init_parameters!(m)
    init_model_indices!(m)

    # Set settings
    model_settings!(m)
    default_test_settings!(m)

    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Test to see if needed values are defined, throw error otherwise
    try
        @assert typeof(get_setting(m, :reduce_state_vars)) == Bool
        @assert typeof(get_setting(m, :reduce_v)) == Bool
    catch
        error("Need to specify in settings whether to perform distribution and/or value function reduction")
    end

    if get_setting(m, :reduce_state_vars) || get_setting(m, :reduce_v)
        try
            @assert get_setting(m, :n_jump_vars) > 0
            @assert get_setting(m, :n_state_vars) > 0
            @assert get_setting(m, :n_state_vars_unreduce) >= 0
        catch
            error("Need to enter values in settings for n_jump_vars, n_state_vars, and/or n_state_vars_unreduce in Settings. These indicate the number of jump and aggregate state variables.")
        end
    end

    if get_setting(m, :reduce_state_vars)
        try
            @assert get_setting(m, :krylov_dim) > 0
        catch
            error("Need to specify in settings desired dimension of Krylov subspace, krylov_dim.")
        end
    end

    if get_setting(m, :reduce_v)
        try
            @assert get_setting(m, :n_prior) > 0 && typeof(get_setting(m, :n_prior)) == Int
            @assert get_setting(m, :n_post) > 0 && typeof(get_setting(m, :n_post)) == Int
            @assert in(:knots_dict, keys(m.settings))
            @assert in(:spline_grid, keys(m.settings))
        catch
            error("Need to specify in settings n_prior, n_post, knots_dict and/or spline_grid for value function reduction. See reduction.jl for their uses.")
        end
    end

    # Solve for steady state
    steadystate!(m)
    return m
end

"""
```
init_parameters!(m::KrusellSmithCT)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::KrusellSmithCT)
    # Initialize parameters
    m <= parameter(:γ, 2., fixed=true,
                   description="γ: The coefficient of relative risk aversion.",
                   tex_label="\\gamma")

    m <= parameter(:ρ, 0.01, fixed=true,
                   description="ρ: The rate of time preference.",
                   tex_label="\\rho")

    m <= parameter(:δ, .025, fixed=true,
                   description="δ: Capital depreciation.",
                   tex_label="\\delta")

    m <= parameter(:α, 1/3, fixed=true,
                   description="α: Capital share.",
                   tex_label="\\alpha")

    m <= parameter(:σ_tfp, .007, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), RootInverseGamma(4, 0.5),
                   fixed=false, description="σ_tfp: The standard deviation of the TFP shock.",
                   tex_label="\\sigma_tfp")

    m <= parameter(:ρ_tfp, .95, fixed=true,
                   description="ρ_tfp: The quarterly autocorrelation of TFP shock.",
                   tex_label="\\rho_tfp")

    m <= parameter(:λ1, 1/2, fixed=true,
                   description="λ1: The expected duration of unemployment is 2 quarters.",
                   tex_label="\\lambda_1")

    m <= parameter(:target_ur, .07, fixed=true,
                   description="target_ur: Target unemployment rate.",
                   tex_label="target_ur")

    m <= parameter(:target_er, .93, fixed=true,
                   description="target_er: Target employment rate.",
                   tex_label="target_er")

    m <= parameter(:λ2, (m[:λ1] / (1. * m[:target_er] - 0.)) * (1. - 1. * m[:target_er]), fixed=true,
                   description="λ2: The unemployment rate is 7.",
                   tex_label="\\lambda_2")

    m <= parameter(:μ, .15, fixed=true,
                   description="μ: The UI replacement rate is 15.",
                   tex_label="\\mu")

    m <= parameter(:τ, (m[:μ] / 1.) * (m[:λ2] / m[:λ1]) , fixed=true,
                   description="τ: The labor income tax.",
                   tex_label="\\tau")

    # Steady state parameters
    m <= SteadyStateParameterGrid(:V_ss, Vector{Float64}(undef, 0); description =
                                  "Stacked steady-state value function")
    m <= SteadyStateParameterGrid(:gg_ss, Vector{Float64}(undef, 0); description =
                                  "Stacked steady-state distribution across state space")
    m <= SteadyStateParameter(:log_aggregate_tfp_ss, NaN, description =
                              "Steady-state level of aggregate shocks")
    m <= SteadyStateParameter(:K_ss, NaN, description = "Steady-state capital level")
    m <= SteadyStateParameter(:r_ss, NaN, description = "Steady-state real interest rate")
    m <= SteadyStateParameter(:w_ss, NaN, description = "Steady-state wage")
    m <= SteadyStateParameter(:Y_ss, NaN, description = "Steady-state output level")
    m <= SteadyStateParameter(:C_ss, NaN, description = "Steady-state consumption level")
    m <= SteadyStateParameter(:I_ss, NaN, description = "Steady-state investment level")
    m <= SteadyStateParameterGrid(:If_ss, Matrix{Float64}(undef, 0, 0); description =
                                   "Steady-state forward difference upwind scheme")
    m <= SteadyStateParameterGrid(:Ib_ss, Matrix{Float64}(undef, 0, 0); description =
                                  "Steady-state backward difference upwind scheme")
    m <= SteadyStateParameterGrid(:I0_ss, Matrix{Float64}(undef, 0, 0); description =
                                  "Steady-state no drift upwind scheme")
end

 """
 ```
 steadystate!(m::KrusellSmithCT)
 ```

 Calculates the model's steady-state values. `steadystate!(m)` must be called whenever
 the parameters of `m` are updated.
 """
function steadystate!(m::KrusellSmithCT)

    # Read in parameters
    γ::Float64     = m[:γ].value
    ρ::Float64     = m[:ρ].value
    δ::Float64     = m[:δ].value
    α::Float64     = m[:α].value
    σ_tfp::Float64 = m[:σ_tfp].value
    ρ_tfp::Float64 = m[:ρ_tfp].value
    μ::Float64     = m[:μ].value # UI replacement rate
    τ::Float64     = m[:τ].value # labor income tax

    # Read in grids
    I::Int64             = get_setting(m, :I)
    J::Int64             = get_setting(m, :J)
    z::Vector{Float64}   = get_setting(m, :z)
    a::Vector{Float64}   = get_setting(m, :a)
    da::Float64          = get_setting(m, :da)
    dz::Float64          = get_setting(m, :dz)
    amax::Float64        = get_setting(m, :amax)
    amin::Float64        = get_setting(m, :amin)
    z_avg::Float64       = (m[:λ1].value * z[2] + m[:λ2].value * z[1]) / (m[:λ1].value + m[:λ2].value)
    aa::Array{Float64,2} = get_setting(m, :aa)
    zz::Array{Float64,2} = get_setting(m, :zz)

    # A_switch * v = lambda_z * (v(a, z') - v(a, z)) in Eqn 8 of Why Inequality Matters
    # basically captures the expected change in value function due to jumps in z
    A_switch = [-SparseMatrixCSC{Float64}(LinearAlgebra.I, I, I) * m[:λ1].value SparseMatrixCSC{Float64}(LinearAlgebra.I, I, I) * m[:λ1].value; SparseMatrixCSC{Float64}(LinearAlgebra.I, I, I) * m[:λ2].value -SparseMatrixCSC{Float64}(LinearAlgebra.I, I, I) * m[:λ2].value]

    # Read in initial rates
    r0::Float64   = get_setting(m, :r0)
    rmin::Float64 = get_setting(m, :rmin)
    rmax::Float64 = get_setting(m, :rmax)

    # Read in  approximation parameters
    Ir::Int64          = get_setting(m, :Ir)
    Δ::Float64         = get_setting(m, :Δ_HJB)
    crit_S::Float64    = get_setting(m, :crit_S)
    crit_HJB::Float64  = get_setting(m, :crit_HJB)
    maxit_HJB::Float64 = get_setting(m, :maxit_HJB)

    # Initialize guesses and difference matrices
    dVf::Array{Float64,2}       = zeros(Float64, I, 2) # forward value function difference matrix
    dVb::Array{Float64,2}       = zeros(Float64, I, 2) # backward " "
    dV_upwind::Array{Float64,2} = zeros(Float64, I, 2) # value function difference matrix post upwind
    c::Array{Float64,2}         = similar(dVf) # flow consumption across state space
    u::Array{Float64,2}         = similar(dVf) # flow utility across state space
    g::Array{Float64,2}         = similar(dVf) # distribution vector
    r::Float64                  = r0  # initial interest rate guess
    KD::Float64 = (((α) / (r + δ)) ^ (1 / (1 - α))) * z_avg # capital demand
    w::Float64  = (1 - α) * (KD ^ α) * ((z_avg) ^ (-α))     # wage rate
    v0::Array{Float64,2} = similar(dVf)                     # initial value function guess
    v0[:, 1] = (w * μ * (1 - z[1]) .+ r .* a) .^ (1 - γ)/(1 - γ)/ρ # add guesses
    v0[:, 2] = (w * (1 - τ) * z[2] .+ r .* a) .^ (1 - γ)/(1 - γ)/ρ

    If::Array{Float64,2} = similar(dVf) # 1 if forward difference, 0 if not
    Ib::Array{Float64,2} = similar(dVf) # 1 if backward difference, 0 if not
    I0::Array{Float64,2} = similar(dVf)  # neither forward nor backward difference
    A::SparseMatrixCSC{Float64,Int64} = spzeros(Float64, 2*I, 2*I)

    # Iterate to find steady state interest rate
    for ir = 1:Ir

        # Get capital demand and wages
        KD = (((α) / (r + δ)) ^ (1 / (1 - α))) * z_avg
        w  = (1 - α) * (KD ^ α) * ((z_avg) ^ (-α))
        v  = v0

        # Solve for value function given interest rate
        for n = 1:maxit_HJB
            # Grab guess
            V = v

            # Compute forward difference
            dVf[1:I - 1, :] = (V[2:I, :] - V[1:(I - 1), :])/da
            dVf[I, :] = (w * ((1 - τ) * z .+ μ * (1 .- z)) .+ r .* amax).^(-γ)

            # Compute backward difference
            dVb[2:I, :] = (V[2:I, :] - V[1:(I - 1), :])/da
            dVb[1, :] = (w * ((1 - τ) * z .+ μ * (1 .- z)) .+ r .* amin).^(-γ)

            # Compute consumption and savings with forward difference
            cf  = dVf.^(-1/γ)
            ssf = w*((1 - τ) * zz + μ * (1 .- zz)) + r .* aa - cf

            # Compute consumption and savings with backward difference
            cb = dVb.^(-1/γ)
            ssb = w*((1 - τ) * zz + μ * (1 .- zz)) + r .* aa - cb

            # Compute consumption and derivative of value function for no drift
            c0::Array{Float64,2}  = w*((1 - τ) * zz .+ μ * (1 .- zz)) + r.*aa
            dV0::Array{Float64,2} = c0.^(-γ)

            # Compute upwind differences and matrix
            dV_upwind, If, Ib, I0 = upwind_value_function(dVf, dVb, dV0, ssf, ssb)
            c = dV_upwind .^ (-1 / γ)
            u = c .^ (1 - γ) / (1 - γ)
            savingsSS = w * ((1 - τ) * zz + μ * (1 .- zz)) + r .* aa .- c

            # Create A matrix and solve for hjb
            A = upwind_matrix(A_switch, ssf, ssb, da, I, J)
            V_stacked = solve_hjb(A, ρ, Δ, vec(u), vec(V))::Vector{Float64}
            V = reshape(V_stacked, I, J)

            # Update values and check convergence
            Vchange = V - v
            v = V

            if maximum(abs.(Vchange)) < crit_HJB
                break
            end
        end
        v0 = v

        # Solve for stationary distribution
        gg = solve_kfe(A, [da; dz], [I; J])::Vector{Float64}
        g = reshape(gg, I, J)

        # Compute capital supply, check market-clearing
        KS = Float64(g[:,1]' * a * da + g[:,2]' * a * da)
        S = Float64(KS - KD)

        # Update interest rate
        if S > crit_S
            # Excess Supply (for capital)
            rmax = r
            r = 0.5 * (r + rmin)
        elseif S < -crit_S
            # Excess Demand
            rmin = r
            r = 0.5 * (r + rmax)
        elseif abs(S) < crit_S
            # Set steady state values
            # Must maintain type of declared SteadyState parameter type
            m[:V_ss]  = vec(v)
            m[:gg_ss] = vec(g)[1:2*I - 1]
            m[:log_aggregate_tfp_ss] = 0.
            m[:K_ss]  = KS
            m[:r_ss]  = r
            m[:w_ss]  = w
            m[:Y_ss]  = (m[:K_ss]^α) * (z_avg ^ (1 - α))
            m[:C_ss]  = sum(c[:] .* g[:] * da)
            m[:I_ss]  = δ * m[:K_ss]
            m[:If_ss] = If
            m[:Ib_ss] = Ib
            m[:I0_ss] = I0
            break
        end
    end
    return m
end

function model_settings!(m::KrusellSmithCT)
    default_settings!(m)

    # State space grid
    m <= Setting(:I, 100, "Dimension of wealth grid")
    m <= Setting(:amin, 0., "Minimium wealth")
    m <= Setting(:amax, 100., "Maximum wealth")
    m <= Setting(:a, collect(range(0., stop=100., length=100)), "Wealth grid")
    m <= Setting(:da, 100. / (100 - 1), "Size of partitions in wealth grid")
    m <= Setting(:aa, [get_setting(m, :a) get_setting(m, :a)],
                 "Repetition of wealth grid across income grid dimension")
    m <= Setting(:aaa, reshape(get_setting(m, :aa), 2 * 100, 1), "Stacked repetition of wealth grid")
    m <= Setting(:J, 2, "Dimension of income grid")
    m <= Setting(:z, [0., 1.], "Income grid")
    m <= Setting(:dz, 1., "Size of partitions in income grid")
    m <= Setting(:zz, ones(100, 1) * get_setting(m, :z)', "Repetition of income grid across wealth dimension")
    m <= Setting(:zzz, reshape(get_setting(m, :zz), 2 * 100, 1), "Stacked repetition of income grid")
    endo = m.endogenous_states

    # Number of variables
    m <= Setting(:n_jump_vars, length(endo[:value_function]), "Number of jump variables")
    m <= Setting(:n_state_vars, length(endo[:distribution]) + n_shocks_exogenous(m),
                 "Number of state variables being reduced")
    m <= Setting(:n_state_vars_unreduce, 0, "Number of state variables not being reduced")

    # Approximation of steady state
    m <= Setting(:rmin, .0001, "Lower bound for steady state interest rate during bisection search")
    m <= Setting(:rmax, m[:ρ].value, "Upper bound for steady state interest rate during bisection search")
    m <= Setting(:r0, .005, "Initial guess for steady state interest rate")
    m <= Setting(:maxit_HJB, 100, "Maximum iterations on steady state HJB")
    m <= Setting(:crit_HJB, 1e-6, "Error criterion for convergence of steady state HJB")
    m <= Setting(:Δ_HJB, 1e4, "Update size for implicit scheme for steady state HJB")
    m <= Setting(:Ir, 100, "Maximum iterations for steady state interest rate")
    m <= Setting(:crit_S, 1e-5, "Error criterion for convergence of steady state interest rate")

    # Reduction
    m <= Setting(:reduce_state_vars, true, "Reduce state variables or not")
    m <= Setting(:reduce_v, true, "Reduce value function or not")
    m <= Setting(:krylov_dim, 5, "Krylov reduction dimension")
    m <= Setting(:F, identity, "Function applied during block arnoldi deflation")
    m <= Setting(:n_knots, 12, "Number of knot points")
    m <= Setting(:c_power, 7, "Amount of bending of knot point locations to make them nonuniform")

    # Create knot points for spline reductoin
    knots = collect(range(get_setting(m, :amin),
                          stop=get_setting(m, :amax), length=get_setting(m, :n_knots) .- 1))
    knots = (get_setting(m, :amax) - get_setting(m, :amin))/(2^get_setting(m, :c_power) - 1) *
        ((knots .- get_setting(m, :amin)) / (get_setting(m, :amax) .- get_setting(m, :amin)) .+ 1) .^
        get_setting(m, :c_power) .+ get_setting(m, :amin) .- (get_setting(m, :amax) .-
                                   get_setting(m, :amin))/(2^get_setting(m, :c_power) - 1)

    # Reduction continued
    m <= Setting(:knots_dict, Dict(1 => collect(range(get_setting(m, :amin),
                                   stop=get_setting(m, :amax), length=(12 - 1)))),
                 "Location of knot points for each dimension for value function reduction")
    m <= Setting(:spline_grid, get_setting(m, :a), "Grids for spline basis")
    m <= Setting(:n_prior, 1,
         "Number of dimensions approximated by spline basis that were not used to compute the basis matrix")
    m <= Setting(:n_post, 2, "Number of dimensions that need to be approximated by spline basis")

    # Transition time step
    m <= Setting(:dt, 1. / 2, "Size of time step for transition equation")

    # Track last lag or not?
    m <= Setting(:track_lag, false, "Add first lag when constructing measurement equation")

    # Sampling Method
    m <= Setting(:sampling_method, :SMC, "Set sampling method to SMC")

    # SMC settings
    m <= Setting(:n_particles, 3_000)
    m <= Setting(:n_Φ, 200)
    m <= Setting(:λ, 2.0)
    m <= Setting(:n_smc_blocks, 1)
    m <= Setting(:use_parallel_workers, true)
    # m <= Setting(:use_parallel_workers, false)
    m <= Setting(:step_size_smc, 0.5)
    m <= Setting(:n_mh_steps_smc, 1)
    m <= Setting(:resampler_smc, :polyalgo)
    m <= Setting(:target_accept, 0.25)
    m <= Setting(:mixture_proportion, .9)
    m <= Setting(:tempering_target, 0.95)
    m <= Setting(:resampling_threshold, .5)
    m <= Setting(:use_fixed_schedule, true)
    m <= Setting(:smc_iteration, 0)

    # Forecast
    m <= Setting(:use_population_forecast, true,
                 "Whether to use population forecasts as data")
    m <= Setting(:forecast_zlb_value, 0.13,
        "Value of the zero lower bound in forecast periods, if we choose to enforce it")

    #Simulating states
    m <= Setting(:state_simulation_freq, 2,
                 "How many states you want to simulate between states + 1")

end

"""
Overloading functions so as to properly grab information from CT Hank Models.
"""
n_states(m::KrusellSmithCT) = sum(map(i -> length(collect(m.endogenous_states)[i][2]), 1:length(keys(m.endogenous_states))))
n_shocks_expectational(m::KrusellSmithCT) = sum(map(i -> length(collect(m.expected_shocks)[i][2]), 1:length(keys(m.expected_shocks))))

#=
n_states(m::OneAssetHANK) = sum(map(i -> length(collect(m.endogenous_states)[i][2]), 1:length(keys(m.endogenous_states))))
n_shocks_expectational(m::OneAssetHANK) = sum(map(i -> length(collect(m.expected_shocks)[i][2]), 1:length(keys(m.expected_shocks))))
=#


text: using TreeLas.Dual: wolfe_gap_step

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
using HDF5, FileIO, JLD2, Test, Random, DSGE, DelimitedFiles
path = dirname(@__FILE__)
writing_output = false

m = AnSchorfheide()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, save)

data = h5read("$path/../../reference/smc.h5", "data")

m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, true)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :polyalgo)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:adaptive_tempering_target_smc, false)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:smc_iteration, 0)
m <= Setting(:use_chand_recursion, true)

@everywhere Random.seed!(42)

println("Estimating AnSchorfheide Model... (approx. 2 minutes)")
DSGE.smc2(m, data, verbose = :none) # us.txt gives equiv to periods 95:174 in our current dataset
println("Estimation done!")

test_file = load(rawpath(m, "estimate", "smc_cloud.jld2"))
test_cloud  = test_file["cloud"]
test_w      = test_file["w"]
test_W      = test_file["W"]
test_z = test_file["z"]

if writing_output
    jldopen("$path/../../reference/smc_cloud_fix=true.jld2", true, true, true, IOStream) do file
        write(file, "cloud", test_cloud)
        write(file, "w", test_w)
        write(file, "W", test_W)
        write(file, "z", test_z)
    end
end

saved_file  = load("$path/../../reference/smc_cloud_fix=true.jld2")
saved_cloud = old_to_new_cloud(DSGE.Cloud(saved_file["cloud"]))
saved_w     = saved_file["w"]
saved_W     = saved_file["W"]

####################################################################
cloud_fields = fieldnames(typeof(test_cloud))
@testset "ParticleCloud Fields: AnSchorf" begin
    @test @test_matrix_approx_eq SMC.get_vals(test_cloud) SMC.get_vals(saved_cloud)
    @test @test_matrix_approx_eq SMC.get_loglh(test_cloud) SMC.get_loglh(saved_cloud)
    @test length(test_cloud.particles) == length(saved_cloud.particles)
    @test test_cloud.tempering_schedule == saved_cloud.tempering_schedule
    @test test_cloud.ESS ≈ saved_cloud.ESS
    @test test_cloud.stage_index == saved_cloud.stage_index
    @test test_cloud.n_Φ == saved_cloud.n_Φ
    @test test_cloud.resamples == saved_cloud.resamples
    @test test_cloud.c == saved_cloud.c
    @test test_cloud.accept == saved_cloud.accept
end

test_particle  = test_cloud.particles[1,:]
saved_particle = saved_cloud.particles[1,:]
N = length(test_particle)
@testset "Individual Particle Fields Post-SMC: AnSchorf" begin
    @test test_particle[1:SMC.ind_para_end(N)] ≈ saved_particle[1:SMC.ind_para_end(N)]
    @test test_particle[SMC.ind_loglh(N)]      ≈ saved_particle[SMC.ind_loglh(N)]
    @test test_particle[SMC.ind_logpost(N)]    ≈ saved_particle[SMC.ind_logpost(N)]
    @test test_particle[SMC.ind_logprior(N)]   ≈ saved_particle[SMC.ind_logprior(N)]
    @test test_particle[SMC.ind_old_loglh(N)] == saved_particle[SMC.ind_old_loglh(N)]
    @test test_particle[SMC.ind_accept(N)]    == saved_particle[SMC.ind_accept(N)]
    @test test_particle[SMC.ind_weight(N)]     ≈ saved_particle[SMC.ind_weight(N)]
end

@testset "Weight Matrices: AnSchorf" begin
    @test @test_matrix_approx_eq test_w saved_w
    @test @test_matrix_approx_eq test_W saved_W
end


text:
using MultivariateOrthogonalPolynomials, DifferentialEquations

Z = Zernike(1)
W = Weighted(Z)
xy = axes(W,1)
Δ = Z \ Laplacian(xy) * W