text:
using Statistics
using LinearAlgebra
using Distributions
using Random

function log_bayesian_posterior(s_param, θ::Array{Float64,1}, forward::Function,
    y::Array{Float64,1},  Σ_η::Array{Float64,2},
    μ0::Array{Float64,1}, Σ0::Array{Float64,2})

    Gu = forward(s_param, θ)
    Φ = - 0.5*(y - Gu)'/Σ_η*(y - Gu) - 0.5*(θ - μ0)'/Σ0*(θ - μ0)
    return Φ

end

function d_log_bayesian_posterior(s_param, θ::Array{Float64,1}, forward::Function,
    y::Array{Float64,1},  Σ_η::Array{Float64,2},
    μ0::Array{Float64,1}, Σ0::Array{Float64,2})

    Gu, dGu = forward(s_param, θ)
    Φ = - 0.5*(y - Gu)'/Σ_η*(y - Gu) - 0.5*(θ - μ0)'/Σ0*(θ - μ0)
    dΦ = dGu'* (Σ_η\(y - Gu)) - Σ0\(θ - μ0)
    return Φ, dΦ

end

function log_likelihood(s_param, θ::Array{Float64,1}, forward::Function,
    y::Array{Float64,1},  Σ_η::Array{Float64,2})

    Gu = forward(s_param, θ)
    Φ = - 0.5*(y - Gu)'/Σ_η*(y - Gu)
    return Φ

end

"""
When the density function is Φ/Z,
The f_density function return log(Φ) instead of Φ
"""
function RWMCMC_Run(log_bayesian_posterior::Function, θ0::Array{FT,1}, τ::FT, n_ite::IT; seed::IT=11) where {FT<:AbstractFloat, IT<:Int}

    N_θ = length(θ0)
    θs = zeros(Float64, n_ite, N_θ)
    fs = zeros(Float64, n_ite)

    θs[1, :] .= θ0
    fs[1] = log_bayesian_posterior(θ0)

    Random.seed!(seed)
    for i = 2:n_ite
        θ_p = θs[i-1, :]
        θ = θ_p + τ * rand(Normal(0, 1), N_θ)

        fs[i] = log_bayesian_posterior(θ)
        α = min(1.0, exp(fs[i] - fs[i-1]))
        if α > rand(Uniform(0, 1))
            # accept
            θs[i, :] = θ
            fs[i] = fs[i]
        else
            # reject
            θs[i, :] = θ_p
            fs[i] = fs[i-1]
        end

    end

    return θs
end

"""
When the density function is Φ/Z,
The f_density function return log(Φ) instead of Φ
"""
function PCN_Run(log_likelihood::Function, θ0::Array{FT,1}, θθ0_cov::Array{FT,2}, β::FT, n_ite::IT; seed::IT=11) where {FT<:AbstractFloat, IT<:Int}

    N_θ = length(θ0)
    θs = zeros(Float64, n_ite, N_θ)
    fs = zeros(Float64, n_ite)

    θs[1, :] .= θ0
    fs[1] = log_likelihood(θ0)

    Random.seed!(seed)
    for i = 2:n_ite
        θ_p = θs[i-1, :]
        θ = sqrt(1 - β^2)*θ_p + β * rand(MvNormal(zeros(size(θ0)), θθ0_cov))

        fs[i] = log_likelihood(θ)
        α = min(1.0, exp(fs[i] - fs[i-1]))
        if α > rand(Uniform(0, 1))
            # accept
            θs[i, :] = θ
            fs[i] = fs[i]
            @info "accept i = ", i
        else
            # reject
            θs[i, :] = θ_p
            fs[i] = fs[i-1]
        end

    end

    return θs
end

# Stein variational gradient descent
function SVGD_Kernel(θ; h = -1)
    J, nd = size(θ)

    XY = θ*θ';
    x2= sum(θ.^2, dims=2);
    X2e = repeat(x2, 1, J);
    pairwise_dists = X2e + X2e' - 2*XY

    if h < 0
        h = median(pairwise_dists)
        h = sqrt(0.5 * h / log(J+1))
    end
    # compute the rbf kernel
    Kxy = exp.( -pairwise_dists / h^2 / 2)

    dxkxy = -Kxy * θ
    sumkxy = sum(Kxy, dims=2)
    for i = 1:nd
        dxkxy[:, i] = dxkxy[:,i] + θ[:,i].*sumkxy
    end

    dxkxy = dxkxy / (h^2)
    return Kxy, dxkxy
end

# lnprob(θ) = ∇ log p(θ)
# θs0 is a J by n_dim matrix
function SVGD_Run(θs0, lnprob, n_ite::IT = 1000, stepsize::FT = 1e-3, alpha = 0.9) where {FT<:AbstractFloat, IT<:Int}
        # Check input

    θs = copy(θs0)
    lnpgrad = similar(θs0)
    J, nd = size(θs0)
    # adagrad with momentum
    fudge_factor = 1e-6
    historical_grad = 0

    for iter = 1:n_ite
        for i = 1:J
            _, lnpgrad[i, :] = lnprob(θs[i, :])
        end
        # calculating the kernel matrix
        kxy, dxkxy = SVGD_Kernel(θs; h = -1)
        grad_theta = (kxy * lnpgrad + dxkxy) / J

        # adagrad
        if iter == 1
            historical_grad = historical_grad .+ grad_theta .^ 2
        else
            historical_grad .= alpha * historical_grad .+ (1 - alpha) * (grad_theta .^ 2)
        end

        adj_grad = grad_theta ./ (fudge_factor .+ sqrt.(historical_grad))
        θs = θs + stepsize * adj_grad
    end

    return θs
end

"""
Metropolis Adjusted Langevin Algorithm
When the density function is Φ/Z,
The f_density function return log(Φ) instead of Φ
"""
function MALA_Run(d_log_bayesian_posterior::Function, θ0::Array{FT,1}, τ::FT, n_ite::IT; seed::IT=11) where {FT<:AbstractFloat, IT<:Int}

    N_θ = length(θ0)
    θs = zeros(Float64, n_ite, N_θ)
    fs = zeros(Float64, n_ite)
    dfs = zeros(Float64, n_ite, N_θ)

    θs[1, :] .= θ0
    fs[1], dfs[1, :] = d_log_bayesian_posterior(θ0)

    Random.seed!(seed)
    for i = 2:n_ite
        θ_p = θs[i-1, :]
        θ = θ_p + τ *dfs[i-1, :] + sqrt(2τ)*rand(Normal(0,1), N_θ)

        fs[i], dfs[i, :] = d_log_bayesian_posterior(θ)
        α = min(1.0, exp( (fs[i]   - norm(θ_p - θ - τ *dfs[i, :])^2/(4τ)  ) -
                          (fs[i-1] - norm(θ - θ_p - τ *dfs[i-1, :])^2/(4τ))
                        )
                )
        # @info α , θ,  fs[i], dfs[i,:], θ_p, fs[i-1], dfs[i-1,:],
        if α > rand(Uniform(0, 1))
            # accept
            θs[i, :] = θ
            fs[i] = fs[i]
            dfs[i, :] = dfs[i, :]
        else
            # reject
            θs[i, :] = θ_p
            fs[i] = fs[i-1]
            dfs[i, :] = dfs[i-1, :]
        end

    end

    return θs
end

# θ_s: N_ens/2 by N_θ matrix, the set
# θ_c: N_ens/2 by N_θ matrix, the complementary set
function emcee_Propose(θ_s::Array{FT,2}, θ_c::Array{FT,2}; a::FT = 2.0) where {FT<:AbstractFloat, IT<:Int}
    Ns, N_θ = size(θ_s)

    zz = ((a - 1.0) * rand(Uniform(0, 1), Ns) .+ 1).^2.0 / a
    factors = (N_θ - 1.0) * log.(zz)

    rint = rand(1:Ns, Ns,)
    return θ_c[rint, :] - (θ_c[rint, :] - θ_s) .* zz, factors

end

"""
When the density function is Φ/Z,
The log_bayesian_posterior function return log(Φ) instead of Φ
θ0 : initial ensemble of size N_ens by N_θ
"""
function emcee_Run(log_bayesian_posterior::Function, θ0::Array{FT,2}, n_ite::IT; random_split::Bool = true, a::FT = 2.0, seed::IT=11) where {FT<:AbstractFloat, IT<:Int}
    Random.seed!(seed)

    N_ens, N_θ = size(θ0)
    @assert(N_ens >= 2N_θ)
    @assert(N_ens % 2 == 0)

    θs = zeros(Float64, n_ite, N_ens, N_θ)
    fs = zeros(Float64, n_ite, N_ens)

    θs[1, :, :] .= θ0
    for k = 1:N_ens
        fs[1, k] = log_bayesian_posterior(θ0[k, :])
    end

    nsplit = 2
    N_s = div(N_ens, 2)

    all_inds = Array(1:N_ens)
    inds = all_inds .% nsplit # 2 group
    log_probs = zeros(Float64, N_s)
    for i_t = 2:n_ite

        if random_split
            shuffle!(inds)
        end
        for split = 0:1
            # boolean array of the current set
            s_inds = (inds .== split)
            # boolean array for the complementary set
            c_inds = (inds .!= split)

            s, c = θs[i_t - 1, s_inds, :], θs[i_t - 1, c_inds, :]

            q, factors = emcee_Propose(s, c; a = a)
            for i = 1:N_s
                log_probs[i] = log_bayesian_posterior(q[i, :])
            end

            # Loop over the walkers and update them accordingly.
            for i = 1:N_s
                j = all_inds[s_inds][i]

                # j is the index
                # @info i, factors[i] , log_probs[i] , fs[i_t - 1, j]

                α = min(1.0, exp(factors[i] + log_probs[i] - fs[i_t - 1, j]))
                if α > rand(Uniform(0, 1))
                    # accept
                    θs[i_t, j, :] = q[i, :]
                    fs[i_t, j] = log_probs[i]
                else
                    # reject
                    θs[i_t, j, :] = θs[i_t - 1, j, :]
                    fs[i_t, j] = fs[i_t - 1, j]
                end
            end
        end
    end

    return θs
end


text:
using FromFile
import LineSearches
import Optim
@from "Core.jl" import CONST_TYPE, Node, Options, Dataset
@from "Utils.jl" import getTime
@from "EquationUtils.jl" import getConstants, setConstants, countConstants
@from "LossFunctions.jl" import scoreFunc
@from "PopMember.jl" import PopMember

# Proxy function for optimization
function optFunc(x::Vector{CONST_TYPE}, dataset::Dataset{T}, baseline::T,
                 tree::Node, options::Options; allow_diff=false)::T where {T<:Real}
    setConstants(tree, x)
    return scoreFunc(dataset, baseline, tree, options; allow_diff=allow_diff)
end

# Use Nelder-Mead to optimize the constants in an equation
function optimizeConstants(dataset::Dataset{T},
                           baseline::T, member::PopMember,
                           options::Options)::PopMember where {T<:Real}

    nconst = countConstants(member.tree)
    if nconst == 0
        return member
    end
    x0 = getConstants(member.tree)
    f(x::Vector{CONST_TYPE})::T = optFunc(x, dataset, baseline, member.tree, options)
    differentiable_f(x::Vector{CONST_TYPE})::T = optFunc(x, dataset, baseline, member.tree, options; allow_diff=true)
    use_differentiable = false
    if nconst == 1
        algorithm = Optim.Newton()
    else
        if options.optimizer_algorithm == "NelderMead"
            algorithm = Optim.NelderMead()
        elseif options.optimizer_algorithm == "BFGS"
            use_differentiable = true
            algorithm = Optim.BFGS(linesearch=LineSearches.BackTracking())#order=3))
        else
            error("Optimization function not implemented.")
        end
    end
    result = if !use_differentiable
        Optim.optimize(f, x0, algorithm, Optim.Options(iterations=options.optimizer_iterations))
    else
         Optim.optimize(differentiable_f, x0, algorithm, Optim.Options(iterations=options.optimizer_iterations))
    end
    # Try other initial conditions:
    for i=1:options.optimizer_nrestarts
        new_start = x0 .* (convert(CONST_TYPE, 1) .+ convert(CONST_TYPE, 1//2)*randn(CONST_TYPE, size(x0, 1)))
        tmpresult = if !use_differentiable
             Optim.optimize(f, new_start, algorithm, Optim.Options(iterations=options.optimizer_iterations))
        else
             Optim.optimize(differentiable_f, new_start, algorithm, Optim.Options(iterations=options.optimizer_iterations))
        end

        if tmpresult.minimum < result.minimum
            result = tmpresult
        end
    end

    if Optim.converged(result)
        setConstants(member.tree, result.minimizer)
        member.score = convert(T, result.minimum)
        member.birth = getTime()
    else
        setConstants(member.tree, x0)
    end
    return member
end

ext: jl


text:
using SeisIO, Noise, JLD2, PlotlyJS

include("../../src/utils.jl")

# input parameters
InputDict = Dict( "corrname"           => "outputData/BPnetworkxcorr_neq.jld2",
                  "refname"            => "refXCorrs/reference_xcorr_neq.jld2",
                  "method"             => "MWCS",
                  "freqmin"            => 0.1,
                  "freqmax"            => 9.9,
                  "fs"                 => 20.0,
                  "maxtimelag"         => 100.0,
                  "mintimelag"         => -100.0,
                  "dtt_lag"            => "static",
                  "dtt_v"              => 1.0,
                  "dtt_minlag"         => -100.0,
                  "dtt_width"          => 10.0,
                  "window_length"      => 10.0,
                  "window_step"        => 5.0,
                  "dtt_sides"          => "Both",
                  "smoothing_half_win" => 5)

# read daily xcorrs and reference xcorr
xcorrs = jldopen(InputDict["corrname"])
reference = jldopen(InputDict["refname"])

tstamplist = xcorrs["info/timestamplist"][:]
stationpairs = reference["info/reference_xcorrnames"]

#iterate over time steps
for tstamp in tstamplist[1:3]
    for stnpair in stationpairs[37:end]
        ref = reference[stnpair]
        cur = xcorrs["$tstamp/$stnpair"]
        println("$tstamp/$stnpair")
        stack!(cur, allstack=true)

        dist = xcorrs["$tstamp/$stnpair"].misc["dist"]

        if InputDict["method"] == "MWCS"
            time_axis, dt, err, coh = mwcs(ref,
                                           cur.corr,
                                           InputDict["freqmin"],
                                           InputDict["freqmax"],
                                           InputDict["fs"],
                                           InputDict["mintimelag"],
                                           InputDict["window_length"],
                                           InputDict["window_step"],
                                           InputDict["smoothing_half_win"])

            m, em, a, ea, m0, em0 = mwcs_dvv(time_axis,
                                             dt,
                                             err,
                                             coh,
                                             InputDict["dtt_lag"],
                                             dist/1000,
                                             InputDict["dtt_v"],
                                             InputDict["dtt_minlag"],
                                             InputDict["dtt_width"],
                                             InputDict["dtt_sides"])

        elseif InputDict["method"] == "stretching"
            # set up time axis
            time_axis = InputDict["mintimelag"]:1/InputDict["fs"]:InputDict["maxtimelag"]

            window_length_samples = convert(Int,InputDict["window_length"] * InputDict["fs"])
            window_step_samples = convert(Int,InputDict["window_step"] * InputDict["fs"])
            minind = 1:window_step_samples:length(ref) - window_length_samples

            N = length(minind)
            dv_list = zeros(N)
            for ii=1:N
                window = collect(minind[ii]:minind[ii]+window_length_samples-1)
                dv, cc, cdp, eps, err, C = stretching(ref[:, 1],
                                                      cur.corr[:, 1],
                                                      time_axis,
                                                      window,
                                                      InputDict["freqmin"],
                                                      InputDict["freqmax"],
                                                      dvmin=-0.03,
                                                      dvmax=0.03)
                dv_list[ii] = dv
            end
            # average dv/V over all windows
            dvV=sum(dv_list)/N
        end
        # save dv/v to a GeoDataFrame (?)
    end
end

close(xcorrs)
close(ref)

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FlowCutterPACE17_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FlowCutterPACE17")
JLLWrappers.@generate_main_file("FlowCutterPACE17", UUID("008204e2-cd5c-5c6d-9360-d31f32b5f6c2"))
end  # module FlowCutterPACE17_jll


text:
using MetaLearning
using Test

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


text:
module HierarchicalGeometry

# Write your package code here.
using StaticArrays
using LinearAlgebra
using Polyhedra
using LazySets
# using QHull

"""
    cross_product_operator(x)

Computes `Ωₓ` such that `Ωₓy = x ⨯ y`
"""
cross_product_operator(x) = SMatrix{3,3}(
    [0.0    -x[3]   x[2];
     x[3]   0.0     -x[1];
     -x[2]  x[1]    0.0]
)

abstract type OverapproxModel end

export
    RegularPolyhedronOverapprox

"""
    RegularPolyhedronOverapprox{D,N}

Used to overrapproximate convex sets with a pre-determined set of support
vectors.
"""
struct RegularPolyhedronOverapprox{D,N} <: OverapproxModel
    support_vectors::NTuple{N,SVector{D,Float64}}
end
get_support_vecs(model::RegularPolyhedronOverapprox) = [v for v in model.support_vectors]

"""
    RegularPolyhedronOverapprox(dim::Int,N::Int,epsilon=0.1)

Construct a regular polyhedron overapproximation model by specifying the number
of dimensions and the number of support vectors to be arranged radially about
the axis formed by the unit vector along each dimension.
epsilon if the distance between support vector v and an existing support vector
is less than epsilon, the new vector will not be added.
"""
function RegularPolyhedronOverapprox(dim::Int,N::Int,epsilon=0.1)
    vecs = Vector{SVector{dim,Float64}}()
    # v - the initial support vector for a given dimension
    v = zeros(dim)
    v[end] = 1.0
    for i in 1:dim
        d = zeros(dim)
        d[i] = 1.0
        A = cross_product_operator(d)
        for j in 1:N
            v = normalize(exp(A*j*2*pi/N)*v)
            add = true
            for vp in vecs
                if norm(v-vp) < epsilon
                    add = false
                    break
                end
            end
            if add
                @show v
                push!(vecs,v)
            end
        end
        v = d
    end
    RegularPolyhedronOverapprox(tuple(vecs...))
end

export
    equatorial_overrapprox_model

"""
    equatorial_overapprox_model(lat_angles=[-π/4,0.0,π/4],lon_angles=collect(0:π/4:2π),epsilon=0.1)

Returns a RegularPolyhedronOverapprox model generated with one face for each
combination of pitch and yaw angles specified by lat_angles and lon_angles,
respectively. There are also two faces at the two poles
"""
function equatorial_overapprox_model(lat_angles=[-π/4,0.0,π/4],lon_angles=collect(0:π/4:2π),epsilon=0.1)
    vecs = Vector{SVector{3,Float64}}()
    push!(vecs,[0.0,0.0,1.0])
    push!(vecs,[0.0,0.0,-1.0])
    for phi in lat_angles
        for theta in lon_angles
            v = normalize([
                cos(phi)*cos(theta),
                cos(phi)*sin(theta),
                sin(phi)
            ])
            add = true
            for vp in vecs
                if norm(v-vp) < epsilon
                    add = false
                    break
                end
            end
            if add
                push!(vecs,v)
            end
        end
    end
    RegularPolyhedronOverapprox(tuple(vecs...))
end

function LazySets.overapproximate(lazy_set,model::RegularPolyhedronOverapprox{N},epsilon::Float64=0.1) where {D,N}
    halfspaces = map(v->LazySets.HalfSpace(v,ρ(v,lazy_set)),get_support_vecs(model))
    sort!(halfspaces; by=h->ρ(h.a, lazy_set))
    poly = HPolyhedron()
    while !isempty(halfspaces) #&& !isbounded(poly)
        poly = intersection(poly,halfspaces[1])
        deleteat!(halfspaces,1)
    end
    @assert isbounded(poly)
    hpoly = convert(HPolytope,poly)
    # vpoly = tovrep(hpoly)
    # continue adding new halfspace constraints until the distance "clipped off"
    # by the new halfspace constraint is less than epsilon
    # while !isempty(halfspaces)
    #     for v in vertices_list(hpoly)
    #
    #     end
    # end
    hpoly
end

# """
#     sort_facet_vtxs
#
# Returns a counterclockwise ordering of 'vtx_ids' around the face defined by
# the halfspace associated with h_idx. It is assumed that the vertices form a
# convex polygon
# """
# function sort_facet_vtxs(poly,h_idx,vtx_idxs)
#     if !isempty(vtx_idxs)
#         vtxs = map(v->get(poly,v), vtx_idxs)
#         h = get(poly,h_idx) # halfspace
#         ctr = sum(vtxs)
#         # axes
#         a1 = normalize(vtxs[1]-ctr)
#         a2 = normalize(cross(h.a,a1))
#         # angles = map(v->atan(dot(a2,v),dot(a1,v)))
#         return sort(vtx_idxs;by=v->atan(dot(a2,get(poly,v)),dot(a1,get(poly,v))))
#     end
#     return vtx_idxs
# end

# """
#     polytope_vtxs_and_edges
#
# Identify the edges of each facet of a polytope
# """
# function polytope_vtxs_and_edges(hpolytope::AbstractPolyhedron)
#     poly = polyhedron(hpolytope)
#     vtxs = collect(points(poly))
#     edges = Vector{Tuple{Int,Int}}()
#     for idx in eachindex(halfspaces(poly))
#         vtx_idxs = incidentpointindices(poly,idx)
#         sorted_vtx_ids = sort_facet_vtxs(poly,idx,vtx_idxs)
#         for i in 1:length(sorted_vtx_ids)-1
#             v1 = sorted_vtx_ids[i]
#             v2 = sorted_vtx_ids[i+1]
#             push!(edges, (v1.value, v2.value))
#         end
#         push!(edges,(sorted_vtx_ids[end].value,sorted_vtx_ids[1].value))
#     end
#     return vtxs, edges
# end
#
# function write_vertices(filename,vtxs::Vector{Vector{Float64}})
#     open(filename, "w") do io
#         for v in vtxs
#             for c in v
#                 print(io,c," ")
#             end
#             print(io,"\n")
#         end
#     end
# end
#
# function write_edges(filename,edge_list::Vector{Tuple{Int,Int}})
#     open(filename, "w") do io
#         for e in edge_list
#             for v in e
#                 print(io,v," ")
#             end
#             print(io,"\n")
#         end
#     end
# end
#
# function write_vtxs_and_edges(filename,poly)
#     vtx_file = join([filename,".nodes"])
#     edge_file = join([filename,".edges"])
#     vtxs, edge_list = polytope_vtxs_and_edges(poly)
#     write_vertices(vtx_file,vtxs)
#     write_edges(edge_file,edge_list)
# end

end


text:
# imag(conj(x) * y) avoiding computing the real part if possible
@inline _imagconjtimes(x, y) = imag(conj(x) * y)
@inline function _imagconjtimes(x::Complex, y::Complex)
    return muladd(-imag(x), real(y), real(x) * imag(y))
end
@inline _imagconjtimes(x::Real, y::Complex) = x * imag(y)
@inline _imagconjtimes(x::Complex, y::Real) = -imag(x) * y
@inline _imagconjtimes(x::Real, y::Real) = ZeroTangent()


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct StringToInt64Map_MapEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StringToInt64Map_MapEntry(; kwargs...)
        obj = new(meta(StringToInt64Map_MapEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StringToInt64Map_MapEntry (mapentry)
const __meta_StringToInt64Map_MapEntry = Ref{ProtoMeta}()
function meta(::Type{StringToInt64Map_MapEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StringToInt64Map_MapEntry)
            __meta_StringToInt64Map_MapEntry[] = target = ProtoMeta(StringToInt64Map_MapEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Int64]
            meta(target, StringToInt64Map_MapEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StringToInt64Map_MapEntry[]
    end
end
function Base.getproperty(obj::StringToInt64Map_MapEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct StringToInt64Map <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StringToInt64Map(; kwargs...)
        obj = new(meta(StringToInt64Map), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StringToInt64Map
const __meta_StringToInt64Map = Ref{ProtoMeta}()
function meta(::Type{StringToInt64Map})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StringToInt64Map)
            __meta_StringToInt64Map[] = target = ProtoMeta(StringToInt64Map)
            allflds = Pair{Symbol,Union{Type,String}}[:map => Base.Dict{AbstractString,Int64}]
            meta(target, StringToInt64Map, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StringToInt64Map[]
    end
end
function Base.getproperty(obj::StringToInt64Map, name::Symbol)
    if name === :map
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Int64}
    else
        getfield(obj, name)
    end
end

mutable struct Int64ToStringMap_MapEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Int64ToStringMap_MapEntry(; kwargs...)
        obj = new(meta(Int64ToStringMap_MapEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Int64ToStringMap_MapEntry (mapentry)
const __meta_Int64ToStringMap_MapEntry = Ref{ProtoMeta}()
function meta(::Type{Int64ToStringMap_MapEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Int64ToStringMap_MapEntry)
            __meta_Int64ToStringMap_MapEntry[] = target = ProtoMeta(Int64ToStringMap_MapEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => Int64, :value => AbstractString]
            meta(target, Int64ToStringMap_MapEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Int64ToStringMap_MapEntry[]
    end
end
function Base.getproperty(obj::Int64ToStringMap_MapEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct Int64ToStringMap <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Int64ToStringMap(; kwargs...)
        obj = new(meta(Int64ToStringMap), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Int64ToStringMap
const __meta_Int64ToStringMap = Ref{ProtoMeta}()
function meta(::Type{Int64ToStringMap})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Int64ToStringMap)
            __meta_Int64ToStringMap[] = target = ProtoMeta(Int64ToStringMap)
            allflds = Pair{Symbol,Union{Type,String}}[:map => Base.Dict{Int64,AbstractString}]
            meta(target, Int64ToStringMap, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Int64ToStringMap[]
    end
end
function Base.getproperty(obj::Int64ToStringMap, name::Symbol)
    if name === :map
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{Int64,AbstractString}
    else
        getfield(obj, name)
    end
end

mutable struct StringToDoubleMap_MapEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StringToDoubleMap_MapEntry(; kwargs...)
        obj = new(meta(StringToDoubleMap_MapEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StringToDoubleMap_MapEntry (mapentry)
const __meta_StringToDoubleMap_MapEntry = Ref{ProtoMeta}()
function meta(::Type{StringToDoubleMap_MapEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StringToDoubleMap_MapEntry)
            __meta_StringToDoubleMap_MapEntry[] = target = ProtoMeta(StringToDoubleMap_MapEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => Float64]
            meta(target, StringToDoubleMap_MapEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StringToDoubleMap_MapEntry[]
    end
end
function Base.getproperty(obj::StringToDoubleMap_MapEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct StringToDoubleMap <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StringToDoubleMap(; kwargs...)
        obj = new(meta(StringToDoubleMap), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StringToDoubleMap
const __meta_StringToDoubleMap = Ref{ProtoMeta}()
function meta(::Type{StringToDoubleMap})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StringToDoubleMap)
            __meta_StringToDoubleMap[] = target = ProtoMeta(StringToDoubleMap)
            allflds = Pair{Symbol,Union{Type,String}}[:map => Base.Dict{AbstractString,Float64}]
            meta(target, StringToDoubleMap, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StringToDoubleMap[]
    end
end
function Base.getproperty(obj::StringToDoubleMap, name::Symbol)
    if name === :map
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,Float64}
    else
        getfield(obj, name)
    end
end

mutable struct Int64ToDoubleMap_MapEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Int64ToDoubleMap_MapEntry(; kwargs...)
        obj = new(meta(Int64ToDoubleMap_MapEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Int64ToDoubleMap_MapEntry (mapentry)
const __meta_Int64ToDoubleMap_MapEntry = Ref{ProtoMeta}()
function meta(::Type{Int64ToDoubleMap_MapEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Int64ToDoubleMap_MapEntry)
            __meta_Int64ToDoubleMap_MapEntry[] = target = ProtoMeta(Int64ToDoubleMap_MapEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => Int64, :value => Float64]
            meta(target, Int64ToDoubleMap_MapEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Int64ToDoubleMap_MapEntry[]
    end
end
function Base.getproperty(obj::Int64ToDoubleMap_MapEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct Int64ToDoubleMap <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Int64ToDoubleMap(; kwargs...)
        obj = new(meta(Int64ToDoubleMap), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Int64ToDoubleMap
const __meta_Int64ToDoubleMap = Ref{ProtoMeta}()
function meta(::Type{Int64ToDoubleMap})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Int64ToDoubleMap)
            __meta_Int64ToDoubleMap[] = target = ProtoMeta(Int64ToDoubleMap)
            allflds = Pair{Symbol,Union{Type,String}}[:map => Base.Dict{Int64,Float64}]
            meta(target, Int64ToDoubleMap, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Int64ToDoubleMap[]
    end
end
function Base.getproperty(obj::Int64ToDoubleMap, name::Symbol)
    if name === :map
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{Int64,Float64}
    else
        getfield(obj, name)
    end
end

mutable struct StringVector <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function StringVector(; kwargs...)
        obj = new(meta(StringVector), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct StringVector
const __meta_StringVector = Ref{ProtoMeta}()
function meta(::Type{StringVector})
    ProtoBuf.metalock() do
        if !isassigned(__meta_StringVector)
            __meta_StringVector[] = target = ProtoMeta(StringVector)
            allflds = Pair{Symbol,Union{Type,String}}[:vector => Base.Vector{AbstractString}]
            meta(target, StringVector, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_StringVector[]
    end
end
function Base.getproperty(obj::StringVector, name::Symbol)
    if name === :vector
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    else
        getfield(obj, name)
    end
end

mutable struct Int64Vector <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Int64Vector(; kwargs...)
        obj = new(meta(Int64Vector), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Int64Vector
const __meta_Int64Vector = Ref{ProtoMeta}()
function meta(::Type{Int64Vector})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Int64Vector)
            __meta_Int64Vector[] = target = ProtoMeta(Int64Vector)
            pack = Symbol[:vector]
            allflds = Pair{Symbol,Union{Type,String}}[:vector => Base.Vector{Int64}]
            meta(target, Int64Vector, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Int64Vector[]
    end
end
function Base.getproperty(obj::Int64Vector, name::Symbol)
    if name === :vector
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    else
        getfield(obj, name)
    end
end

mutable struct FloatVector <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function FloatVector(; kwargs...)
        obj = new(meta(FloatVector), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct FloatVector
const __meta_FloatVector = Ref{ProtoMeta}()
function meta(::Type{FloatVector})
    ProtoBuf.metalock() do
        if !isassigned(__meta_FloatVector)
            __meta_FloatVector[] = target = ProtoMeta(FloatVector)
            pack = Symbol[:vector]
            allflds = Pair{Symbol,Union{Type,String}}[:vector => Base.Vector{Float32}]
            meta(target, FloatVector, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_FloatVector[]
    end
end
function Base.getproperty(obj::FloatVector, name::Symbol)
    if name === :vector
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float32}
    else
        getfield(obj, name)
    end
end

mutable struct DoubleVector <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function DoubleVector(; kwargs...)
        obj = new(meta(DoubleVector), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct DoubleVector
const __meta_DoubleVector = Ref{ProtoMeta}()
function meta(::Type{DoubleVector})
    ProtoBuf.metalock() do
        if !isassigned(__meta_DoubleVector)
            __meta_DoubleVector[] = target = ProtoMeta(DoubleVector)
            pack = Symbol[:vector]
            allflds = Pair{Symbol,Union{Type,String}}[:vector => Base.Vector{Float64}]
            meta(target, DoubleVector, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_DoubleVector[]
    end
end
function Base.getproperty(obj::DoubleVector, name::Symbol)
    if name === :vector
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    else
        getfield(obj, name)
    end
end

mutable struct Int64Range <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Int64Range(; kwargs...)
        obj = new(meta(Int64Range), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Int64Range
const __meta_Int64Range = Ref{ProtoMeta}()
function meta(::Type{Int64Range})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Int64Range)
            __meta_Int64Range[] = target = ProtoMeta(Int64Range)
            allflds = Pair{Symbol,Union{Type,String}}[:minValue => Int64, :maxValue => Int64]
            meta(target, Int64Range, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Int64Range[]
    end
end
function Base.getproperty(obj::Int64Range, name::Symbol)
    if name === :minValue
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :maxValue
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct Int64Set <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Int64Set(; kwargs...)
        obj = new(meta(Int64Set), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Int64Set
const __meta_Int64Set = Ref{ProtoMeta}()
function meta(::Type{Int64Set})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Int64Set)
            __meta_Int64Set[] = target = ProtoMeta(Int64Set)
            pack = Symbol[:values]
            allflds = Pair{Symbol,Union{Type,String}}[:values => Base.Vector{Int64}]
            meta(target, Int64Set, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Int64Set[]
    end
end
function Base.getproperty(obj::Int64Set, name::Symbol)
    if name === :values
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    else
        getfield(obj, name)
    end
end

mutable struct DoubleRange <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function DoubleRange(; kwargs...)
        obj = new(meta(DoubleRange), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct DoubleRange
const __meta_DoubleRange = Ref{ProtoMeta}()
function meta(::Type{DoubleRange})
    ProtoBuf.metalock() do
        if !isassigned(__meta_DoubleRange)
            __meta_DoubleRange[] = target = ProtoMeta(DoubleRange)
            allflds = Pair{Symbol,Union{Type,String}}[:minValue => Float64, :maxValue => Float64]
            meta(target, DoubleRange, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_DoubleRange[]
    end
end
function Base.getproperty(obj::DoubleRange, name::Symbol)
    if name === :minValue
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :maxValue
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

export StringToInt64Map_MapEntry, StringToInt64Map, Int64ToStringMap_MapEntry, Int64ToStringMap, StringToDoubleMap_MapEntry, StringToDoubleMap, Int64ToDoubleMap_MapEntry, Int64ToDoubleMap, StringVector, Int64Vector, FloatVector, DoubleVector, Int64Range, Int64Set, DoubleRange
# mapentries: "StringToInt64Map_MapEntry" => ("AbstractString", "Int64"), "Int64ToDoubleMap_MapEntry" => ("Int64", "Float64"), "Int64ToStringMap_MapEntry" => ("Int64", "AbstractString"), "StringToDoubleMap_MapEntry" => ("AbstractString", "Float64")


text:
#*******************************************************/
#* Copyright(c) 2018 by Artelys                        */
#* This source code is subject to the terms of the     */
#* MIT Expat License (see LICENSE.md)                  */
#*******************************************************/

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use Knitro to solve the following
# simple mathematical program with equilibrium/complementarity
# constraints(MPEC/MPCC).
#
#  min  (x0 - 5)^2 +(2 x1 + 1)^2
#  s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2
#        x2 complements(3 x0 - x1 - 3)
#        x3 complements(-x0 + 0.5 x1 + 4)
#        x4 complements(-x0 - x1 + 7)
#        x0, x1, x2, x3, x4 >= 0
#
# The complementarity constraints must be converted so that one
# nonnegative variable complements another nonnegative variable.
#
#  min  (x0 - 5)^2 +(2 x1 + 1)^2
#  s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2  (c0)
#        3 x0 - x1 - 3 - x5 = 0                 (c1)
#        -x0 + 0.5 x1 + 4 - x6 = 0              (c2)
#        -x0 - x1 + 7 - x7 = 0                  (c3)
#        x2 complements x5
#        x3 complements x6
#        x4 complements x7
#        x0, x1, x2, x3, x4, x5, x6, x7 >= 0
#
# The solution is(1, 0, 3.5, 0, 0, 0, 3, 6), with objective value 17.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using KNITRO, Test

function example_mpec1(; verbose=true)
    # Create a new Knitro solver instance.
    kc = KNITRO.KN_new()

    # Illustrate how to override default options by reading from
    # the knitro.opt file.
    options = joinpath(dirname(@__FILE__), "..", "examples", "knitro.opt")
    KNITRO.KN_load_param_file(kc, options)

    # Initialize Knitro with the problem definition.

    # Add the variables and set their bounds and initial values.
    # Note: unset bounds assumed to be infinite.
    KNITRO.KN_add_vars(kc, 8)
    KNITRO.KN_set_var_lobnds(kc, zeros(Float64, 8))
    KNITRO.KN_set_var_primal_init_values(kc, zeros(Float64, 8))

    # Add the constraints and set their bounds.
    KNITRO.KN_add_cons(kc, 4)
    KNITRO.KN_set_con_eqbnds(kc, Float64[2, 3, -4, -7])

    # Add coefficients for all linear constraints at once.

    # c0
    lconIndexCons = Int32[0, 0, 0, 0, 0]
    lconIndexVars = Int32[0, 1, 2, 3, 4]
    lconCoefs = [-1.5, 2.0, 1.0, -0.5, 1.0]

    # c1
    lconIndexCons = [lconIndexCons; Int32[1, 1, 1]]
    lconIndexVars = [lconIndexVars; Int32[0, 1, 5]]
    lconCoefs = [lconCoefs; [3.0, -1.0, -1.0]]

    # c2
    lconIndexCons = [lconIndexCons; Int32[2, 2, 2]]
    lconIndexVars = [lconIndexVars; Int32[0, 1, 6]]
    lconCoefs = [lconCoefs; [-1.0, 0.5, -1.0]]

    # c3
    lconIndexCons = [lconIndexCons; Int32[3, 3, 3]]
    lconIndexVars = [lconIndexVars; Int32[0, 1, 7]]
    lconCoefs = [lconCoefs; [-1.0, -1.0, -1.0]]

    KNITRO.KN_add_con_linear_struct(kc, lconIndexCons, lconIndexVars, lconCoefs)

    # Note that the objective(x0 - 5)^2 +(2 x1 + 1)^2 when
    # expanded becomes:
    #    x0^2 + 4 x1^2 - 10 x0 + 4 x1 + 26

    # Add quadratic coefficients for the objective
    qobjIndexVars1 = Int32[0, 1]
    qobjIndexVars2 = Int32[0, 1]
    qobjCoefs = [1.0, 4.0]
    KNITRO.KN_add_obj_quadratic_struct(kc, qobjIndexVars1, qobjIndexVars2, qobjCoefs)

    # Add linear coefficients for the objective
    lobjIndexVars = Int32[0, 1]
    lobjCoefs = [-10.0, 4.0]
    KNITRO.KN_add_obj_linear_struct(kc, lobjIndexVars, lobjCoefs)

    # Add constant to the objective
    KNITRO.KN_add_obj_constant(kc, 26.0)

    # Set minimize or maximize(if not set, assumed minimize)
    KNITRO.KN_set_obj_goal(kc, KNITRO.KN_OBJGOAL_MINIMIZE)

    # Now add the complementarity constraints
    ccTypes = [KNITRO.KN_CCTYPE_VARVAR, KNITRO.KN_CCTYPE_VARVAR, KNITRO.KN_CCTYPE_VARVAR]
    indexComps1 = Int32[2, 3, 4]
    indexComps2 = Int32[5, 6, 7]
    KNITRO.KN_set_compcons(kc, ccTypes, indexComps1, indexComps2)

    kn_outlev = verbose ? KNITRO.KN_OUTLEV_ALL : KNITRO.KN_OUTLEV_NONE
    KNITRO.KN_set_param(kc, KNITRO.KN_PARAM_OUTLEV, kn_outlev)

    # Solve the problem.
    #
    # Return status codes are defined in "knitro.h" and described
    # in the Knitro manual.
    nStatus = KNITRO.KN_solve(kc)

    # An example of obtaining solution information.
    nStatus, objSol, x, lambda_ = KNITRO.KN_get_solution(kc)

    if verbose
        println("Knitro converged with final status = ", nStatus)
        println("  optimal objective value  = ", objSol)
        println("  optimal primal values x0=", x[1])
        println("                        x1=", x[2])
        println("                        x2=", (x[3], x[6]))
        println("                        x3=", (x[4], x[7]))
        println("                        x4=", (x[5], x[8]))
        println("  feasibility violation    = ", KNITRO.KN_get_abs_feas_error(kc))
        println("  KKT optimality violation = ", KNITRO.KN_get_abs_opt_error(kc))
    end

    # Delete the Knitro solver instance.
    KNITRO.KN_free(kc)

    @testset "Example MPEC 1" begin
        @test nStatus == 0
        @test objSol ≈ 17.
        @test x ≈ [1., 0., 3.5, 0., 0., 0., 3., 6.]
    end
end

example_mpec1(; verbose=isdefined(Main, :KN_VERBOSE) ? KN_VERBOSE : true)


text:
export k_antibodies_extinction, empirical_log_ext, k_antibodies_extinction_penalized

# since ab_profile is measured in transitions, actually need ne to be 1/2 θ_ts
# θ_ts = 2 n/λ * (β *mu) = 2 n mu /2 = n mu
# fitness: σf = (2 * ne) .* f = θ_ts f  => f = β0/(β mu) =>
# This implies f = should be 2 β_0 /λ mu *
function k_antibodies_extinction(par::PopRates, k, ab_profile; samples = 10^5)
    f = mean_fit(par)
    θ_ref = pop_θ(par)
    λ = par.λ
    rf = (1 - (λ - f)/(λ + f)) * par.κ / θ_ref # rf = σf / θ_ref
    k_antibodies_extinction( θ_ref, k, rf, ab_profile; samples = samples)
end

function k_antibodies_extinction_penalized(par::PopRates, k, ab_profile; samples = 10^5)
    f = mean_fit(par)
    θ_ref = pop_θ(par)
    λ = par.λ
    rf = (1 - (λ - f)/(λ + f)) * par.κ / θ_ref # rf = σf / θ_ref
    Δ_beta = (1/f - 1 /(λ+f)) * (par.ν * λ) / 2 #
    sampler_array = map(x->WrightSampler(x...),
        hcat([θ_ref.* ab_profile for  _ in 1:k]...))
    fit_array = map(x->x[3],
        hcat([θ_ref.* ab_profile for  _ in 1:k]...)) #return an array of fitnesses
    mean( exp(- θ_ref * rf * penalized_escape_freq(sample!.(sampler_array), fit_array ; Δ_beta = Δ_beta)) for ii in 1:samples)
end

function penalized_escape_freq(freq_array, fit_array; Δ_beta = 0) # penalize according to fitness loss
    prod(
        prod(freq_array[ii] .+ (1 .- freq_array[ii]) .* exp.(-Δ_beta .* fit_array[ii]) ) -
        prod(freq_array[ii])
        for ii in eachindex(freq_array))
end

function prob_of_extinction(σf, sampler_array, samples)
    mean( exp(- σf * escape_freq(sample!.(sampler_array))) for ii in 1:samples)
end

#  θ_ref is the base diversity parameter measured in the reference temporal units.
# rf is the ratio of the bare fitness (growth rate) to the mutation rate

function k_antibodies_extinction(θ_ref,k,rf,ab_profile; samples = 10^5)
    sampler_array = map(x->WrightSampler(x...),
        hcat([θ_ref.* ab_profile for  _ in 1:k]...))
    prob_of_extinction(θ_ref .* rf, sampler_array, samples)
end

function escape_freq(freq_array)
    prod(1 .- prod(freq_array, dims = 1))
end

## function of simulations using each individuals exact growth rate

function empirical_log_ext(pop::PopState, par::PopRates;
    sites = [[1,2],[3]] # Example encoding: must escape from and(or(1,2),3)
    )
    escape_fn(genome::BitArray) =  Bool(prod(
        1 - prod(1-genome[ii] for ii in jj) # 1 if one of the sites is escaped
        for jj in sites)) # 1 if all the sites have atleast 1 escape
    survivors = filter(escape_fn, pop.individuals)
    birth_rates = (par.λ .+ par.f.(survivors)) ./ 2
    death_rates = (par.λ .- par.f.(survivors)) ./ 2
    log_ext = sum(log.(death_rates) .- log.(birth_rates))
end

# Only applicable to a single site, after simulation
# function df_log_ext(df, par::PopRates, fit;
#     sites = 1 # must escape from and(or(1,2),3)
#     )
#     out = log(1 - 2 fit / par.λ) * df.pop_size[] .* df.freq[ind]
# end

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
using Test
using OMEinsum
using OMEinsum: subindex, dynamic_indexer, DynamicEinCode, StaticEinCode, getixs, getiy, labeltype

@testset "EinCode" begin
    code = EinCode(((1,2), (2,3)), (1,3))
    @test code isa EinCode
    @test OMEinsum.getixs(code) == [[1,2], [2,3]]
    @test OMEinsum.getiy(code) == [1,3]

    code1 = ein"ab,bc->ac"
    code2 = EinCode((('a', 'b'), ('b', 'c')), ('a', 'c'))
    @test code2 isa DynamicEinCode
    @test code1 isa StaticEinCode
    @test DynamicEinCode(code1) == code2
    @test StaticEinCode(code2) === code1
    @test collect(collect.(getixs(code1))) == getixs(code2)
    @test collect(getiy(code1)) == getiy(code2)
    @test labeltype(code1) == labeltype(code2)

    code1 = ein",->"
    code2 = EinCode(((), ()), ())
    @test collect(collect.(getixs(code1))) == getixs(code2)
    @test collect(getiy(code1)) == getiy(code2)
    @test labeltype(code1) == labeltype(code2)

    code1 = ein"->"
    code2 = EinCode(((),), ())
    @test collect(collect.(getixs(code1))) == getixs(code2)
    @test collect(getiy(code1)) == getiy(code2)
    @test labeltype(code1) == labeltype(code2)

    @test_throws ErrorException EinCode((), ())
end

@testset "indexer" begin
    si = EinIndexer{()}(())
    @test subindex(si, (1,2,3)) == 1
    si = EinIndexer{(3,2)}((7,6))
    @test OMEinsum.getlocs(si) == (3,2)
    a = randn(7,6)
    @test a[subindex(si, (4,5,2))] == a[2,5]
end

@testset "dynamic indexer" begin
    si = dynamic_indexer((), ())
    @test subindex(si, (1,2,3)) == 1
    si = dynamic_indexer((3,2), (7,6))
    @test OMEinsum.getlocs(si) == (3,2)
    a = randn(7,6)
    @test a[subindex(si, (4,5,2))] == a[2,5]
end

@testset "EinArray" begin
    locs_xs = (EinIndexer{(1,2)}((8,8)), EinIndexer{(2,3)}((8,8)))
    ixs = ((1,2), (2,3))
    iy = (1,3)
    x1 = randn(8, 8)
    x2 = randn(8, 8)
    arr = einarray(Val(ixs), Val(iy), (x1, x2), OMEinsum.get_size_dict!(ixs, (x1, x2), Dict{Int,Int}()))
    @test size(arr) == (8,8,8)
    # inner first, then outer
    @test arr[CartesianIndex((3,4,2))] == x1[4,3]*x2[3,2]
    @test arr[3,4,2] == x1[4,3]*x2[3,2]
end


text:
@testset "Test storage and update_storage!" begin

    Nx = 2
    Ne = 8

    X    =  [0.267333   1.43021;
              0.364979   0.607224;
             -1.23693    0.249277;
             -2.0526     0.915629;
             -0.182465   0.415874;
              0.412907   1.01672;
              1.41332   -0.918205;
              0.766647  -1.00445]';

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]

    for b in Blist

      B = MultiBasis(b, Nx)

      idx = [0 0; 0 1; 1 0; 1 1; 1 2]
      truncidx = idx[1:2:end,:]
      Nψ = 5

      coeff =  [0.6285037650645056;
       -0.4744029092496623;
        1.1405280011620331;
       -0.7217760771894809;
        0.11855056306742319]
      f = ExpandedFunction(B, idx, coeff)
      R = IntegratedFunction(f)

      H = HermiteMapComponent(R; α = 1e-2)

      S = Storage(H.I.f, X);

      @test norm(S.ψoff     -  evaluate_offdiagbasis(f, X)) < 1e-8
      @test norm(S.ψoffψd   -  evaluate_offdiagbasis(f, X) .*
                               evaluate_diagbasis(f, X)) < 1e-8
      @test norm(S.ψoffψd0  -  evaluate_offdiagbasis(f, X) .*
                               repeated_evaluate_basis(f, zeros(Ne))) < 1e-8
      @test norm(S.ψoffdψxd -  evaluate_offdiagbasis(f, X) .*
                               repeated_grad_xk_basis(f, X[end,:])) < 1e-8
      @test norm(S.ψnorm - 1/sqrt(Ne)*norm.(eachcol(evaluate_basis(f, X)))) < 1e-8

       #test cache dimension
      @test size(S.cache_dcψxdt) == (Ne, Nψ)
      @test size(S.cache_dψxd) == (Ne,)
      @test size(S.cache_integral) == (Ne + Ne * Nψ,)
      @test size(S.cache_g) == (Ne,)

      # Test add new components via update_storage!
      addedidx = [2 1; 2 2; 2 3; 3 2]

      S = update_storage(S, X, addedidx)

      addednψ = size(addedidx,1)
      fnew = ExpandedFunction(f.B, vcat(f.idx, addedidx), vcat(f.coeff, zeros(addednψ)))

      @test norm(S.ψoff     -  evaluate_offdiagbasis(fnew, X)) < 1e-8
      @test norm(S.ψoffψd   -  evaluate_offdiagbasis(fnew, X) .*
                               evaluate_diagbasis(fnew, X)) < 1e-8
      @test norm(S.ψoffψd0  -  evaluate_offdiagbasis(fnew, X) .*
                               repeated_evaluate_basis(fnew, zeros(Ne))) < 1e-8
      @test norm(S.ψoffdψxd -  evaluate_offdiagbasis(fnew, X) .*
                               repeated_grad_xk_basis(fnew, X[end,:])) < 1e-8
      @test norm(S.ψnorm - 1/sqrt(Ne)*norm.(eachcol(evaluate_basis(fnew, X)))) < 1e-8

      newNψ = addednψ + Nψ
      #test cache dimension
     @test size(S.cache_dcψxdt) == (Ne, newNψ)
     @test size(S.cache_dψxd) == (Ne,)
     @test size(S.cache_integral) == (Ne + Ne * newNψ,)
     @test size(S.cache_g) == (Ne,)
   end
end


text:
using ExprOptimization
using Test

include("pcfg_test.jl")
include("ppt_test.jl")

include("mc_test.jl")
include("gp_test.jl")
include("ge_test.jl")
include("ce_test.jl")
include("pipe_test.jl")


text:
using LatticeQCD
using Test

@testset "LatticeQCD.jl" begin
    eps = 1e-1

    @testset "quenched HMC" begin

        @testset "SU(2)" begin
            @time plaq = run_LQCD("test02.jl")
            plaq_comparison = 0.5575491312570713
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end

        @testset "SU(3)" begin
            @time plaq = run_LQCD("test01.jl")
            plaq_comparison = 0.6190393357419764
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end

        @testset "SU(4)" begin
            @time plaq = run_LQCD("test03.jl")
            plaq_comparison = 0.4966683811089479
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end
    end

    @testset "Heatbath" begin
        @testset "SU(2)" begin
            @time plaq = run_LQCD("test02-hb.jl")
            plaq_comparison = 0.5287735727118359
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end

        @testset "SU(3)" begin
            @time plaq = run_LQCD("test01-hb.jl")
            plaq_comparison = 0.5821680570717788
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end

        @testset "SU(4)" begin
            @time plaq = run_LQCD("test03-hb.jl")
            plaq_comparison = 0.5467724338528576
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end
    end

    @testset "HMC" begin

        @testset "WilsonClover SU(3)" begin
            @time plaq = run_LQCD("test_wilsonclover.jl")
            plaq_comparison = 0.00624053999484795
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end

        @testset "Staggered SU(3) with 4 tastes" begin
            @time plaq = run_LQCD("test_staggered.jl")
            plaq_comparison = 0.00624053999484795
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end

        @testset "Staggered SU(3) with 2 tastes" begin
            @time plaq = run_LQCD("test_Nf2.jl")
            plaq_comparison = 0.5630198767336069
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end

        @testset "Staggered SU(3) with 3 tastes" begin
            @time plaq = run_LQCD("test_Nf3.jl")
            plaq_comparison = 0.565176584402352
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end
    end

    @testset "SLMC" begin

        @testset "Staggered SU(3)" begin
            @time plaq = run_LQCD("test06_slmc_ks.jl")
            plaq_comparison = 0.5279623287035208
            @test abs(plaq - plaq_comparison)/plaq_comparison < eps
        end
    end

end


text:
module Curves

import ..CryptoGroups: <|, tobits, a, b, cofactor, order, name

include("elliptic_curves.jl")
include("ecpoint.jl")

end


text:
module AutoTM

# Setup environmental variables for nGraph
function __init__()
    setup_affinities()
    setup_codegen()
    setup_passes()

    # Setup PMEM if that option is enabled in nGraph.jl
    settings = nGraph.settings()
    pmm_enabled = get(settings, "PMDK", false)
    pmm_enabled && setup_pmem()

    # just construct a CPU backend to make sure the CPU library is properly loaded
    #
    # Otherwise, some code that lives in the CPU portion of the nGraph code base will not
    # be loaded/links and we error :(
    backend = nGraph.Backend("CPU")

    # DEPRECATRED: This is needed for a newer version of ngraph - but I reverted back
    # to the ASPLOS version of ngraph because the newer version was too unstable for
    # development.
    # # If pmm is enabled, configure the backend to use our custom allocator.
    # # pmm_enabled && nGraph.Lib.set_pmm_allocator(nGraph.getpointer(backend))
end

# stdlibs
using Dates, Random, Serialization, Statistics

# deps
using nGraph
using nGraph: TensorDescriptor, NodeDescriptor, inputs, outputs, description

using JuMP
using JSON
using LightGraphs
using IterTools
using ProgressMeter
using DataStructures
using Flux

# Determine if we're using Gurobi as the ILP solver.
settings = JSON.parsefile(joinpath(dirname(@__DIR__), "deps", "build.json"))
if settings["GUROBI"] == true
    using Gurobi
    getoptimizer() = with_optimizer(Gurobi.Optimizer; TimeLimit = 3600, MIPGap = 0.01)
else
    using Cbc
    getoptimizer() = with_optimizer(Cbc.Optimizer; seconds = 3600, ratioGap = 1.01)
end

# docstringing it up!
using DocStringExtensions

Backend(args...) = nGraph.Backend(args...)

#####
##### Core Functionality
#####

include("setup.jl")
include("utils/utils.jl")
include("profiler/profiler.jl")
include("optimizer/optimizer.jl")
include("verifier/verifier.jl")

# Bring in some utils stuff
using .Utils: Actualizer, actualize, @closure
using .Verifier: verify

#####
##### Predefined models
#####

include("zoo/zoo.jl")

#####
##### Plotting
#####

#include("visualizer/visualizer.jl")

#####
##### Experiments
#####

include("experiments/experiments.jl")

end # module

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
#=
This file implements StatsModels constructors for TableCols.  That is,
it allows one to pass a StatsModels object directly to the TableCol
constructor (or, since the Table constructor falls back to TableCols on
a vararg, the Table constructor).
=#

########################################################################
#################### General Regression API Implementation #############
########################################################################

RegCol{M} = TableCol{3,M}

function RegCol(header::Printable)
    return TableCol(TableIndex(1, header), TableDict{3, FormattedNumber}())
end

"""
```
getnext(t::TableCol{3,M}, group::Int, level::Int) where M
```
In a TableCol `t` of row depth 3, computes the next index on the third level
given that the first dimension of the index is `group` and the second is
`level`.
"""
function getnext(t::RegCol, group::Int, level::Int)
    max_idx = 0
    for ridx in keys(t.data)
        ridx.idx[1] != group && continue
        ridx.idx[2] != level && continue
        max_idx = max(max_idx, ridx.idx[3])
    end
    return max_idx + 1
end

const block_labs = Dict(:setcoef! =>"Coefficient",
                        :setmeta! =>"Metadata",
                        :setstats! =>"Statistics")

for (block, fname) in zip([:1,:2,:3], [:setcoef!, :setmeta!, :setstats!])

    """
    ```
    $fname(t::$RegCol, key, val[, se]; level=1, name="")
    $fname(t::$RegCol, key=>val; level=1, name="")
    $fname(t::$RegCol, kv::AbstractDict)
    ```
    Inserts into `t` a key/value pair (possibly with a standard error) within
    the block.  Like the `TableCol` constructor, the pairs can be passed as
    either individual key/value[/se] tuples or pairs, as several vectors of
    key/value[/se] pairs, or as an associative.

    To add additional sub-blocks, use the `level` keyword argument.  Integers
    less than 0 will appears in blocks above the standard block, and integers
    greater than 1 will appear below it.

    To name the block or sub-block, pass a nonempty string as the `name` keyword
    argument.
    """
    @eval function ($fname)(t::RegCol, key::Printable, val; level=1, name="")
        next_idx = getnext(t, $block, level)
        index    = TableIndex(($block, level, next_idx), ("", name, key))
        t[index] = val
    end

    @eval function ($fname)(t::RegCol, key::Printable, val, se; level=1, name="")
        next_idx = getnext(t, $block, level)
        index    = TableIndex(($block, level, next_idx), ("", name, key))
        t[index] = val, se
    end

    @eval function ($fname)(t::RegCol, p::Pair; level=1, name="")
        next_idx = getnext(t, $block, level)
        key      = p.first
        val      = p.second
        index    = TableIndex(($block, level, next_idx), ("", name, key))
        t[index] = val
    end

    # Handle associatives
    @eval function ($fname)(t::RegCol, args...)
        for kv in zip(args)
            ($fname)(t, kv...)
        end
    end

    @eval function ($fname)(t::RegCol, ps::AbstractDict; level=1, name="")
        ($fname)(t, ps...)
    end
end

########################################################################
#################### Linear Model Interface ############################
########################################################################

function TableCol(header, m::RegressionModel;
                  stats=(:N=>Int∘nobs, "\$R^2\$"=>r2),
                  meta=(), stderror::Function=StatsBase.stderror, kwargs...)

    # Compute p-values
    pval(m) = ccdf.(FDist(1, dof_residual(m)),
                    abs2.(coef(m)./stderror(m)))

    # Initialize the column
    col  = RegCol(header)

    # Add the coefficients
    for (name, val, se, p) in zip(coefnames(m), coef(m), stderror(m), pval(m))
        setcoef!(col, name, val, se)
        0.05 <  p <= .1  && star!(col[name], 1)
        0.01 <  p <= .05 && star!(col[name], 2)
                p <= .01 && star!(col[name], 3)
    end

    # Add in the fit statistics
    setstats!(col, OrderedDict(p.first=>p.second(m) for p in tuplefy(stats)))

    # Add in the metadata
    setmeta!(col, OrderedDict(p.first=>p.second(m) for p in tuplefy(meta)))

    return col
end

########################################################################
#################### regtable Interface ###############################
########################################################################

TableAble = Union{RegressionModel, TexTable, Pair, Tuple}

function regtable(args::Vararg{TableAble}; num=1, kwargs...)
    cols = TexTable[]
    for arg in args
        new_tab = regtable(arg; num=num, kwargs...)
        n, m = size(new_tab)
        num += m
        push!(cols, new_tab)
    end
    return hcat(cols...)
end

function regtable(t::RegressionModel; num=1, kwargs...)
    return TableCol("($num)", t; kwargs...)
end

function regtable(t::TexTable; kwargs...)
    return t
end

function regtable(p::Pair; kwargs...)
    return join_table(p.first=>regtable(p.second; kwargs...))
end

function regtable(p::Tuple; kwargs...)
    return regtable(p...; kwargs...)
end


text:
#### Pattern matching

@inline alwaystrue(x) = true

# matches one term
# syntax:  ~x
struct Slot{P}
    name::Symbol
    predicate::P
end

Slot(s) = Slot(s, alwaystrue)

Base.isequal(s1::Slot, s2::Slot) = s1.name == s2.name

Base.show(io::IO, s::Slot) = (print(io, "~"); print(io, s.name))

# matches zero or more terms
# syntax: ~~x
struct Segment{F}
    name::Symbol
    predicate::F
end

ismatch(s::Segment, t) = s.predicate(t)

Segment(s) = Segment(s, alwaystrue)

Base.show(io::IO, s::Segment) = (print(io, "~~"); print(io, s.name))

makesegment(s::Symbol, keys) = (push!(keys, s); Segment(s))

function makesegment(s::Expr, keys)
    if !(s.head == :(::))
        error("Syntax for specifying a segment is ~~x::\$predicate, where predicate is a boolean function")
    end

    name = s.args[1]

    push!(keys, name)
    :(Segment($(QuoteNode(name)), $(esc(s.args[2]))))
end

makeslot(s::Symbol, keys) = (push!(keys, s); Slot(s))

function makeslot(s::Expr, keys)
    if !(s.head == :(::))
        error("Syntax for specifying a slot is ~x::\$predicate, where predicate is a boolean function")
    end

    name = s.args[1]

    push!(keys, name)
    :(Slot($(QuoteNode(name)), $(esc(s.args[2]))))
end

function makepattern(expr, keys)
    if expr isa Expr
        if expr.head === :call
            if expr.args[1] === :(~)
                if expr.args[2] isa Expr && expr.args[2].args[1] == :(~)
                    # matches ~~x::predicate
                    makesegment(expr.args[2].args[2], keys)
                else
                    # matches ~x::predicate
                    makeslot(expr.args[2], keys)
                end
            else
                :(term($(map(x->makepattern(x, keys), expr.args)...); type=Any))
            end
        else
            error("Unsupported Expr of type $(expr.head) found in pattern")
        end
    else
        # treat as a literal
        return esc(expr)
    end
end

function makeconsequent(expr)
    if expr isa Expr
        if expr.head === :call
            if expr.args[1] === :(~)
                if expr.args[2] isa Symbol
                    return :(getindex(__MATCHES__, $(QuoteNode(expr.args[2]))))
                elseif expr.args[2] isa Expr && expr.args[2].args[1] == :(~)
                    @assert expr.args[2].args[2] isa Symbol
                    return :(getindex(__MATCHES__, $(QuoteNode(expr.args[2].args[2]))))
                end
            else
                return Expr(:call, map(makeconsequent, expr.args)...)
            end
        else
            return Expr(expr.head, map(makeconsequent, expr.args)...)
        end
    else
        # treat as a literal
        return esc(expr)
    end
end

### Matching procedures
# A matcher is a function which takes 3 arguments
# 1. Expression
# 2. Dictionary
# 3. Callback: takes arguments Dictionary × Number of elements matched
#
function matcher(val::Any)
    function literal_matcher(data, bindings, next)
        !isempty(data) && isequal(car(data), val) ? next(bindings, 1) : nothing
    end
end

function matcher(slot::Slot)
    function slot_matcher(data, bindings, next)
        isempty(data) && return
        val = bindings[slot.name]
        if val !== nothing
            if isequal(val, car(data))
                return next(bindings, 1)
            end
        else
            if slot.predicate(car(data))
                next(assoc(bindings, slot.name, car(data)), 1)
            end
        end
    end
end

# returns n == offset, 0 if failed
function trymatchexpr(data, value, n)
    if isempty(value)
        return n
    elseif islist(value) && islist(data)
        if isempty(data)
            # didn't fully match
            return nothing
        end

        while isequal(car(value), car(data))
            n += 1
            value = cdr(value)
            data = cdr(data)

            if isempty(value)
                return n
            elseif isempty(data)
                return nothing
            end
        end

        return isempty(value) ? n : nothing
    elseif isequal(value, data)
        return n + 1
    end
end

function matcher(segment::Segment)
    function segment_matcher(data, bindings, success)
        val = bindings[segment.name]

        if val !== nothing
            n = trymatchexpr(data, val, 0)
            if n !== nothing
                success(bindings, n)
            end
        else
            res = nothing

            for i=length(data):-1:0
                subexpr = take_n(data, i)

                if segment.predicate(subexpr)
                    res = success(assoc(bindings, segment.name, subexpr), i)
                    if res !== nothing
                        break
                    end
                end
            end

            return res
        end
    end
end

function matcher(term::Term)
    matchers = (matcher(operation(term)), map(matcher, arguments(term))...,)
    function term_matcher(data, bindings, success)

        isempty(data) && return nothing
        !(car(data) isa Term) && return nothing

        function loop(term, bindings′, matchers′) # Get it to compile faster
            if isempty(matchers′)
                if  isempty(term)
                    return success(bindings′, 1)
                end
                return nothing
            end
            res = car(matchers′)(term, bindings′,
                                 (b, n) -> loop(drop_n(term, n), b, cdr(matchers′)))
        end

        loop(car(data), bindings, matchers) # Try to eat exactly one term
    end
end


text:
using LaTeXStrings
include("workspace.jl")
pyplot()

N = 27
a = 4

#Calculate the plotting data
xvals = 0:N-1
yvals = [mod(a^x,N) for x ∈ xvals]

plot(xvals, yvals, linewidth=1, ls = :dot, marker = :x,
xlabel = "\$x\$", ylabel = "\$$a^x\$ mod \$$N\$",
title = "Periodicity of modular exponentiation", legend = false, grid = false, mc = :red, lc = :red)
xticks!(0:5:xvals[end]+1)
gui()

r = 0
i = 2
while i < length(xvals)
    if yvals[1] == yvals[i]
        return r = i-1
        break
    end
    if i == length(xvals)
        println("No period was found")
        break
    end
    i += 1
end
plot!([0:1:r,r:-1:0],ones(length(0:1:r),2),
arrow = :arrow, lc = :black)
annotate!((r/2,1.5, Plots.text("\$r = $r\$", 8,:black,:center)))

##
function controlbits(a::Int, n::Int)
    # a is the starting random integer, n is the number of qubits in the circuit
    digits(a,base=2,pad=n)
end

function phaseadd!(C::Circuit, controlbits::AbstractArray{Int})
    n = length(controlbits)
    for m in 1:n
        for i in m:n
            if controlbits[n-i+1] == 1
                k = i-m+1
                λ = 2π/(2^(k))
                U!(C,λ,n-m+1)
            end
        end
    end

   # for qubit in 1:div(n,2)-1
   #     Swap!(C, qubit, n-qubit)
   # end
end

function phaseadd_dagger!(C::Circuit, controlbits::AbstractArray{Int})
    n = length(controlbits)
    for m in n:-1:1
        for i in n:-1:m
            if controlbits[n-i+1] == 1
                λ = 2π/2^(i-m+1)
                U!(C,-λ,n-m+1)
            end
        end
    end
end

#for m in 1:n
#    for i in m:n
#        println(n-i+1," ",i-m+1)
#    end
#end
#
#for m in 5:-1:1
#    for i in 5:-1:m
#        println(5-i+1," ",i-m+1)
#    end
#end

function qft_dagger!(C::Circuit, n::Int)
    #Don't forget the swaps!
    #for qubit in 1:div(n,2)-1
    #    Swap!(C, qubit, n-qubit)
    #end
    for j in 1:n-1
        for m in 1:j-1
            CU!(C,-π/(2^(j-m)), m, j)
        end
        H!(C,j)
    end
end

function qft!(C::Circuit, n::Int)
    #Don't forget the swaps!
    for j in n-1:-1:1
        H!(C,j)
        for m in j-1:-1:1
            CU!(C,π/(2^(j-m)), m, j)
        end
    end
    #for qubit in 1:div(n,2)-1
    #    Swap!(C, qubit, n-qubit)
    #end
end

n = 3
c1 = Circuit(n)
a = 3

#X!(c1,1)
#X!(c1,2)
controls = controlbits(a,n-1)
qft!(c1,n)
phaseadd!(c1,controls)
#phaseadd_dagger!(c1,controls)
qft_dagger!(c1,n)
#qft_dagger!(c1,n)
#qft!(c1,n)
dat = shots(c1, 10_000)
graph!(dat)


text:
include("generateSyntheticData.jl")

function compute_objective_value(X, Y, U, mu, lambda)
    """
    This function computes the objective value of the optimization problem
    for the solution given by (X, Y). The objective value is given by:
        ||U - X - Y||_F^2 + lambda * ||X||_F^2 + mu * ||Y||_F^2

    :param X: An n-by-n matrix.
    :param Y: An n-by-n matrix.
    :param U: An n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
               (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).

    :return: The objective value (Float64).
    """
    val = norm(U-X-Y)^2 + lambda * norm(X)^2 + mu * norm(Y)^2

    return val

end;

function is_feasible(X, Y, k_sparse, k_rank; epsilon = 1e-10)
    """
    This function verifies that the solution (X, Y) is feasible under the
    constraints rank(X) <= k_rank and ||Y||_0 <= k_sparse. Returns true if
    these constraints are satisfied, false otherwise.

    :param X: An n-by-n matrix.
    :param Y: An n-by-n matrix.
    :paeam k_sparse: The maximum number of non-zeros elements in Y (Int64).
    :param k_rank: The maximum rank of the X (Int64).

    :return: True if the constraints are satisfied, false otherwise (Bool).
    """
    # Verify the sparsity constraint
    num_nonzero = count(abs.(Y) .> epsilon)
    if num_nonzero > k_sparse
        return false
    end

    # Verify the rank constraint
    rank_X = rank(X)
    if rank_X > k_rank
        return false
    end

    return true

end;

function project_matrix(A, k_rank; exact_svd=true)
    """
    This function projects the matrix A onto its first k_rank principal
    components.

    :param A: An arbitrary n-by-n matrix.
    :param k_rank: The desired rank of the projection (Int64).
    :param exact_svd: If true, compute an exact truncated SVD. If false,
                      approximate the SVD using a randomized algorithm (Bool).

    :return: The n-by-n low rank projection of A.
    """
    if exact_svd

        factors = tsvd(A, k_rank)
        U = factors[1]
        S = factors[2]
        V = factors[3]

        projected_matrix = U * Diagonal(S) * V'

    else

        projected_matrix = pqrfact(A, rank=k_rank)
        projected_matrix = Matrix(projected_matrix)

    end

    return projected_matrix

end;

function construct_binary_matrix(A, k_sparse; zero_indices=[], one_indices=[])
    """
    This function constructs a binary matrix S of the same shape as A where
    S_ij = 1 if and only if A_ij is one of the k_sparse largest entries of A in
    absolute value.

    :param A: An arbitrary n-by-n matrix.
    :param k_sparse: The desired sparsity of the output matrix (Int64).
    :param zero_indices: List of 2-tuples of Int64 where each entry consists of
                         an index (i, j) such that S_ij is constrained to take
                         value 0.
    :param one_indices: List of 2-tuples of Int64 where each entry consists of
                        an index (i, j) such that S_ij is constrained to take
                        value 1.

    :return: An n-by-n binary matrix.
    """
    n = size(A)[1]

    num_zeros = size(zero_indices)[1]
    num_ones = size(one_indices)[1]

    # Verify that zero_indices and one_indices have a valid length
    @assert num_ones <= k_sparse
    @assert num_zeros <= n * n - k_sparse

    reverse_sorted_entries = sortperm(vec(broadcast(abs, A)), rev=true)

    S = zeros(n, n)
    for (i, j) in one_indices
        S[i, j] = 1
    end

    num_selected_indices = num_ones
    reverse_sorted_index = 1

    while num_selected_indices < k_sparse
        index = reverse_sorted_entries[reverse_sorted_index]
        i = Int64((index - 1) % n + 1)
        j = Int64(floor((index - 1) / n)) + 1

        if S[i, j] == 0 && !((i, j) in zero_indices)
            S[i, j] = 1
            num_selected_indices = num_selected_indices + 1
        end
        reverse_sorted_index = reverse_sorted_index + 1
    end

    return S

end;

function solve_sparse_problem(U_tilde, mu, k_sparse;
    zero_indices=[], one_indices=[])
    """
    This function exactly solves the problem of approximating (in squared
    frobenius norm) an arbitrary matrix U_tilde by a sparse matrix and under a
    regularization penalty.

    :param U_tilde: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
               (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).
    :param zero_indices: List of 2-tuples of Int64 where each entry consists of
                         an index (i, j) such that Y_ij is constrained to take
                         value 0.
    :param one_indices: List of 2-tuples of Int64 where each entry consists of
                        an index (i, j) such that S_ij is constrained to take
                        value 1 where S is the binary matrix denoting the
                        sparsity pattern of Y.

    :return: The n-by-n sparse matrix that solves the problem.
    """
    S = construct_binary_matrix(U_tilde, k_sparse,
        zero_indices = zero_indices, one_indices = one_indices)

    Y = (S .* U_tilde) ./ (1 + mu)

    return Y

end;

function solve_rank_problem(U_tilde, lambda, k_rank; exact_svd=true)
    """
    This function solves the problem of approximating (in squared
    frobenius norm) an arbitrary matrix U_tilde by a low rank matrix and under
    a regularization penalty.

    :param U_tilde: An arbitrary n-by-n matrix.
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).
    :param exact_svd: If true, compute an exact truncated SVD. If false,
                      approximate the SVD using a randomized algorithm (Bool).

    :return: The n-by-n low rank matrix that solves the problem.
    """
    return project_matrix(U_tilde ./ (1 + lambda), k_rank, exact_svd=exact_svd)

end;

function iterate_X_Y(U, mu, lambda, k_sparse, k_rank, X_init, Y_init;
    zero_indices = [], one_indices = [], min_improvement=0.001, exact_svd=true)
    """
    This function computes a feasible solution to the problem
        min ||U - X - Y||_F^2 + lambda * ||X||_F^2 + mu * ||Y||_F^2
        subject to rank(X) <= k_rank, ||Y||_0 <= k_sparse

    by iteratively solving the sparse subproblem and the rank subproblem
    starting from an initial feasible solution (X_init, Y_init).

    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
                   (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).
    :param X_init: An n-by-n matrix with rank at most k_rank.
    :param Y_init: An n-by-n matrix with sparsity at most k_sparse.
    :param zero_indices: List of 2-tuples of Int64 where each entry consists of
                         an index (i, j) such that Y_ij is constrained to take
                         value 0.
    :param one_indices: List of 2-tuples of Int64 where each entry consists of
                        an index (i, j) such that S_ij is constrained to take
                        value 1 where S is the binary matrix denoting the
                        sparsity pattern of Y.
    :param min_improvement: The minimal fractional decrease in the objective
                            value required for the procedure to continue
                            iterating.
    :param exact_svd: If true, compute exact SVDs. If false, compute
                      approximate SVDs (Bool).

    :return: This function returns three values. The first value is a tuple of
             two n-by-n arrays that correspond to the feasible solution found
             by this method (the first element in the tuple is the matrix X and
             the second element is the matrix Y). The second value is the
             objective value (Float64) achieved by the returned feasible
             solution. The third value is the number of iterations performed by
             this method (Int64).
    """
    old_objective = compute_objective_value(X_init, Y_init, U, mu, lambda)

    Y = solve_sparse_problem(U - X_init, mu, k_sparse,
                             zero_indices=zero_indices, one_indices=one_indices)
    X = solve_rank_problem(U - Y, lambda, k_rank, exact_svd=exact_svd)

    new_objective = compute_objective_value(X, Y, U, mu, lambda)
    step_count = 1

    # If the last step improved the solution more than min_provement,
    # take another step
    while((old_objective - new_objective) / old_objective > min_improvement)

        Y = solve_sparse_problem(U - X, mu, k_sparse,
                        zero_indices=zero_indices, one_indices=one_indices)
        X = solve_rank_problem(U - Y, lambda, k_rank, exact_svd=exact_svd)

        old_objective = new_objective
        new_objective = compute_objective_value(X, Y, U, mu, lambda)
        step_count = step_count + 1

    end

    return (X, Y), new_objective, step_count

end;

# TODO: edit function description to describe zero_indices and one_indices params
function SLR_AM(U, mu, lambda, k_sparse, k_rank; zero_indices=[],
    one_indices=[], random_restarts=100, min_improvement=0.001, exact_svd=true)
    """
    This function computes a feasible solution to the problem
        min ||U - X - Y||_F^2 + lambda * ||X||_F^2 + mu * ||Y||_F^2
        subject to rank(X) <= k_rank, ||Y||_0 <= k_sparse

    by randomly initializing a feasible solution (X, Y) and iteratively solving
    the rank subproblem and the sparsity subproblem until convergence. This is
    repeated random_restarts times and the best solution is returned.

    :param U: An arbitrary n-by-n matrix.
    :param mu: The regularization parameter for the sparse matrix penalty
                   (Float64).
    :param lambda: The regularization parameter for the low rank matrix penalty
                   (Float64).
    :param k_sparse: The maximum sparsity of the sparse matrix (Int64).
    :param k_rank: The maximum rank of the low rank matrix (Int64).
    :param zero_indices: List of 2-tuples of Int64 where each entry consists of
                         an index (i, j) such that Y_ij is constrained to take
                         value 0.
    :param one_indices: List of 2-tuples of Int64 where each entry consists of
                        an index (i, j) such that S_ij is constrained to take
                        value 1 where S is the binary matrix denoting the
                        sparsity pattern of Y.
    :param random_restarts: The number of times to repeat the optimizization
                            process from a random initial feasible solution
                            (Int64).
    :param min_improvement: The minimal fractional decrease in the objective
                            value required for the minimization procedure to
                            continue iterating.
    :param exact_svd: If true, compute exact SVDs. If false, compute
                      approximate SVDs (Bool).

    :return: This function returns two values. The first value is a tuple of two
             n-by-n arrays that correspond to the best feasible solution found
             by this method (the first element in the tuple is the matrix X and
             the second element is the matrix Y). The second value is the
             objective value (Float64) achieved by the returned feasible
             solution.
    """
    # Initialize X and Y
    n = size(U)[1]

    X_init = zeros(n, n)
    Y_init = zeros(n, n)

    best_sol, best_obj, _ = iterate_X_Y(U, mu, lambda, k_sparse, k_rank,
                                        X_init, Y_init,
                                        zero_indices=zero_indices,
                                        one_indices=one_indices,
                                        min_improvement=min_improvement,
                                        exact_svd=exact_svd)

    # Repeat the search process random_restarts times
    for i=1:random_restarts
        X_random = rand(Float64, (n, n))
        X_init = project_matrix(X_random, k_rank)
        Y_random = rand(Float64, (n, n))
        Y_init = construct_binary_matrix(Y_random, k_sparse) .* Y_random

        new_sol, new_obj, _ = iterate_X_Y(U, mu, lambda, k_sparse, k_rank,
                                          X_init, Y_init,
                                          zero_indices=zero_indices,
                                          one_indices=one_indices,
                                          min_improvement=min_improvement,
                                          exact_svd=exact_svd)

        # Store the newly found solution if it is better than the best found
        # solution to this point
        if new_obj < best_obj
            best_obj = new_obj
            best_sol = new_sol
        end
    end

    # Verify that the solution to be returned is in fact feasible
    @assert is_feasible(best_sol[1], best_sol[2], k_sparse, k_rank)

    return best_sol, best_obj

end;

function validate_AM_params(n, param_frac, sigma; signal_to_noise=10,
                            exact_svd=true, train_size=20,
                            candidate_lambdas=[10, 1, 0.1, 0.01],
                            candidate_mus=[10, 1, 0.1, 0.01],
                            fixed_rank_sparse=false,
                            fixed_rank=1, fixed_sparse=0,
                            symmetric_flag=false)
    """
    This function performs cross validation to select the regularization
    parameters lambda and mu for the optimization problem given by
        min ||U - X - Y||_F^2 + lambda * ||X||_F^2 + mu * ||Y||_F^2
        subject to rank(X) <= k_rank, ||Y||_0 <= k_sparse

    :param n: The row and column dimension of the data (Int64).
    :param param_frac: Parameter that controls the rank and sparsity of the
                       sampled low rank and sparse matrices when
                       fixed_rank_sparse=false (Float64).
    :param sigma: Parameter that controls the absolute magnitude of the noise
                  (Float64).
    :param signal_to_noise: Parameter that controls the signal to noise ratio of
                            the sampled data (Float64).
    :param exact_svd: If true, compute exact SVDs. If false, compute
                      approximate SVDs (Bool).
    :param train_size: The number of samples to draw per parameter when
                       performing cross validation (Int64).
    :param candidate_lambdas: List of values of lambda (Float64) to be
                              considered during cross validation.
    :param candidate_mus: List of values of mu (Float64) to be considered during
                          cross validation.
    :param fixed_rank_sparse: If true, use fixed_rank and fixed_sparse as the
                              rank and sparsity respectively of the low rank
                              matrix and the sparse matrix when sampling (Bool).
    :param fixed_rank: Rank of the low rank matrices to be sampled when
                       fixed_rank_sparse=true.
    :param fixed_sparse: Sparsity of the sparse matrices to be sampled when
                         fixed_rank_sparse=true
    :param symmetric_flag: If true, sampling symmetric data (Bool).

    :return: This function returns 2 values.
             1) The best performing lambda value (Float64).
             2) The best performing mu value (Float64).
    """

    param_scores = Dict()
    for lambda_mult in candidate_lambdas, mu_mult in candidate_mus

        param_scores[(lambda_mult, mu_mult)] = 0
        this_lambda = lambda_mult / n^0.5
        this_mu = mu_mult / n^0.5

        for trial = 1:train_size

            if symmetric_flag

                D, L_0, S_0, k_sparse, k_rank = generate_synthetic_data_symmetric(n,
                                                    param_frac, sigma,
                                                    signal_to_noise=signal_to_noise,
                                                    fixed_rank_sparse=fixed_rank_sparse,
                                                    fixed_rank=fixed_rank,
                                                    fixed_sparse=fixed_sparse)

            else

                D, L_0, S_0, k_sparse, k_rank = generate_synthetic_data(n,
                                                    param_frac, sigma,
                                                    signal_to_noise=signal_to_noise,
                                                    fixed_rank_sparse=fixed_rank_sparse,
                                                    fixed_rank=fixed_rank,
                                                    fixed_sparse=fixed_sparse)
            end

            sol, _ = SLR_AM(D, this_mu, this_lambda, k_sparse, k_rank,
                            exact_svd=exact_svd)

            D_error = norm(D - sol[1] - sol[2])^2 / norm(D)^2
            L_error = norm(L_0 - sol[1])^2 / norm(L_0)^2
            S_error = norm(S_0 - sol[2])^2 / norm(S_0)^2
            score = D_error + L_error + S_error

            param_scores[(lambda_mult, mu_mult)] += score

        end

        param_scores[(lambda_mult, mu_mult)] /= train_size

    end

    best_score = 1000
    best_params = ()
    for (param, score) in param_scores
        if score < best_score
            best_score = score
            best_params = param
        end
    end

    return best_params[1], best_params[2]

end;


text:
"""
Original publication:
V. Barnabé-Lortie, C. Bellinger, and N. Japkowicz.
Active learning for one-class classification. In 2015
IEEE 14th International Conference on Machine
Learning and Applications (ICMLA), pages 390–395,
Dec 2015.
"""
struct HighConfidencePQs <: ModelBasedPQs
    occ::SVDD.OCClassifier
end

function qs_score(qs::HighConfidencePQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    return SVDD.predict(qs.occ, x)
end

function qs_score(qs::HighConfidencePQs,
                  x::Array{T, 2},
                  labels::Dict{Symbol, Array{Int, 1}},
                  subspaces::Vector{Vector{Int}}) where T <: Real
    size(x,1) == size(qs.occ.data,1) || throw(DimensionMismatch("Number of dimensions in x must match the number of dimensions of the classifier."))
    return map(idx -> SVDD.predict(qs.occ, x[subspaces[idx], :], idx), eachindex(subspaces))
end


text:
using SurrogateModelOptim
using Test

include("infill.jl")
include("interface.jl")
include("smopt.jl")
include("surrogate.jl")
include("utils.jl")


text:
function cardinal_rotations()
    six_dirs = [let o=zeros(3);o[i] = sign*1.0;o end for sign in [-1,1] for i in 1:3]
    cardinal_rots = []
    for dir in six_dirs
        for dir2 in six_dirs
            if dot(dir, dir2) == 0.0
                dir3 = cross(dir, dir2)
                rot = vcat(dir', dir2', dir3')
                push!(cardinal_rots, RotMatrix3{Float64}(rot))
            end
        end
    end
    cardinal_rots
end

function rotations_90()
    [
        let
            rot = zeros(3)
            rot[i] = sign * pi/2
            RotXYZ(rot...)
        end
        for i in 1:3 for sign in [-1, 1]
    ]
end

export cardinal_rotations, rotations_90