text:
module Variational

import Base: start, done, next

export transposecalls!

# aux methods to hide the separation in head and args
numchds(xp) = length(xp.args) + 1
child(xp, idx) = idx == 1 ? xp.head : xp.args[idx-1]

type DepthFirst
    xp::Expr
end

"""
    depthfirst(xp)

Returns an iterator that visits all nodes in an Expr depth first. head and
args are visited before the Expr they define.
"""
depthfirst(xp) = DepthFirst(xp)

type DepthFirstState
    val
    par
    idx
end

function start(itr::DepthFirst)
    head = DepthFirstState(itr.xp, nothing, -1)
    return DepthFirstState(itr.xp, head, 0)
end

done(itr::DepthFirst, state::DepthFirstState) = (state.par == nothing)

function next(itr::DepthFirst, state::DepthFirstState)

    # if all children processed, move up on level
    if state.idx == numchds(state.val)
        return state.val, state.par
    end

    # move to next child
    state = DepthFirstState(state.val, state.par, state.idx+1)

    # if the next child is a leaf, return it
    chd = child(state.val, state.idx)
    if !isa(chd, Expr)
        return chd, state
    end

    # the next child is an expression; descend into it
    # and find the next valid state. There will always
    # be at least one, pointing to the child itself.
    state = DepthFirstState(chd, state, 0)
    return next(itr, state)
end

"""
    transposecall!(xp, skip=[])

Goes through the syntax tree and replace all function calls
`f(x,p1,p2,...)` with `x(f,p1,p2,...)`.
"""
function transposecalls!(xp, skip=[])
    isa(xp, Expr) || return xp
    for x in depthfirst(xp)
        if isa(x, Expr) && x.head == :call && !(x.args[1] in skip)
            @assert length(x.args) >= 2
            tmp = x.args[1]
            x.args[1] = x.args[2]
            x.args[2] = tmp
        end
    end
    return xp
end

import Base.==

export hilbertspace, @hilbertspace
export @varform, Equation
export LinForm, LinTerm, BilForm, BilTerm

import Base: +, -, *, dot, getindex, ^, print

type HilbertVector
    idx
    space
    opstack
end

type LinForm
  test_space
  terms
end

type LinTerm
  test_id
  test_ops
  coeff
  functional
end

type BilForm
  test_space
  trial_space
  terms
end

type BilTerm
  test_id
  trial_id
  test_ops
  trial_ops
  coeff
  kernel
end

type Equation
    lhs
    rhs
end

"""
    ==(lhs::BilForm, rhs::LinForm)

Build an equation from a left hand and right hand side
"""
==(lhs::BilForm, rhs::LinForm) = Equation(lhs, rhs)

"""
    hilbert_space(type, g1, g2, ...)

Returns generators defining a Hilbert space of field `type`
"""
hilbertspace(vars::Symbol...) = [HilbertVector(i, [vars...], []) for i in 1:length(vars)]

macro hilbertspace(syms...)

    for sym in syms
        @assert isa(sym, Symbol) "@hilbertspace takes a list of Symbols"
    end

    rhs = :(hilbertspace())
    for sym in syms
        push!(rhs.args, QuoteNode(sym))
    end

    vars = gensym()
    xp = quote
        $vars = $rhs
    end
    for (i,s) in enumerate(syms)
        push!(xp.args, :(global $s = $vars[$i]))
    end

    xp
end

"""
    call(u::HilbertVector, f, params...)
    u(f, params...)

Add another operation to the opstack of `u`.
"""
(u::HilbertVector)(f, params...) = HilbertVector(u.idx, u.space, [(f, params...); u.opstack])

"""
    getindex(f, v::HilbertVector)
    f[v]

Return a LinForm corresponding to f[v]
"""
getindex(f, v::HilbertVector) = LinForm(v.space, [LinTerm(v.idx, v.opstack, 1, f)])

"""
    getindex(A, v::HilbertVector, u::HilbertVector)

Create a BilForm corresponding to A[v,u]
"""
function getindex(A, v::HilbertVector, u::HilbertVector)
    terms = [ BilTerm(v.idx, u.idx, v.opstack, u.opstack, 1, A) ]
    BilForm(v.space, u.space, terms)
end

"Add two BilForms together"
function +(a::BilForm, b::BilForm)
    @assert a.test_space == b.test_space
    @assert a.trial_space == b.trial_space
    BilForm(a.test_space, a.trial_space, [a.terms; b.terms])
end

function+(a::LinForm, b::LinForm)
    @assert a.test_space == b.test_space
    LinForm(a.test_space, [a.terms; b.terms])
end

function *(α::Number, a::BilForm)
  b = deepcopy(a)
  for t in b.terms t.coeff *= α end
  return b
end

function *(α::Number, a::LinForm)
    b = deepcopy(a)
    for t in b.terms t.coeff *= α end
    return b
end

-(a::BilForm) = (-1 * a)
-(a::BilForm, b::BilForm) = a + (-b)

-(a::LinForm) = (-1 * a)
-(a::LinForm, b::LinForm) = a + (-b)

function print(io::IO, v::HilbertVector)
    sym = v.space[v.idx]
    ops = v.opstack
    for op in ops
        print(io, op[1], "(")
    end
    print(io, sym)
    for op in reverse(ops)
        for p in op[2:end]
            print(io, ", ", p)
        end
        print(io, ")")
    end
end

function print(io::IO, a::LinForm)
    N = length(a.terms)
    #T = typeof(a.terms[1].coeff)
    for (n,t) in enumerate(a.terms)
      u = HilbertVector(t.test_id, a.test_space, t.test_ops)
      t.coeff != 1 && print(io, t.coeff, "*")
      print(io, t.functional, "[", u, "]")
      n == N || print(io, " + ")
  end
end

function print(io::IO, f::BilForm)
    N = length(f.terms)
    #T = typeof(f.terms[1].coeff)

    for (n,t) in enumerate(f.terms)
        u = HilbertVector(t.test_id, f.test_space, t.test_ops)
        v = HilbertVector(t.trial_id, f.trial_space, t.trial_ops)
        t.coeff != 1 && print(io, t.coeff, "*")
        print(io, t.kernel, "[", u, ", ", v, "]")
        n == N || print(io, " + ")
    end
end

function print(io::IO, eq::Equation)
    print(io, eq.lhs)
    print(io, " == ")
    print(io, eq.rhs)
end

"""
    @varform <form-definition>

The Julia form compiler uses the Julia parser and meta-programming
based traversal of the AST to create a structure containing all
information required for the description of a variational formulation
from an Expr that follows closely widespread mathematical convention.

E.g:

    EFIE = @varform T[k,j] = e[k]
    MFIE = @varform 0.5*I[k,j] + K[k,j] = h[k]
    PMCH = @varform M[k,j] - η*T[k,m] + 1/η*T[l,j] + M[l,m] = e[k] + h[l]
"""
macro varform(x)
    y = transposecalls!(x, [:+, :-, :*, :(==)])
    esc(y)
end

end


text:
using Test, ClusterAnalysis
using CSV, DataFrames
using Random

@testset "ClusterAnalysis" begin

include("kmeans.jl");
include("dbscan.jl");

end

ext: jl


text:
using sb, LinearAlgebra, Statistics, StaticArrays

struct storages
    time :: Vector{Float64}
    energy :: Vector{Float64}
    function storages()
        new(Vector{Float64}(),Vector{Float64}())
    end
end

function store!(storage, t, ops, dotops, sbm)
    append!(storage.time,t)
    append!(storage.energy  ,sb.CEIDenergy(ops,sbm))
end

function my_run(knondiag, thermostat)
    storage = storages()
    # Create a model
    sbm = SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 7.5, α = 0.1, No = 60, β = 5.0))
    # Create containers starting from a localized state
    ops = CEIDOps(SA[1.0 0.0; 0.0 0.0],sbm)
    # run dynamics
    RunCEID!(sbm,ops,2000,0.001,store!, storage, thermostat = thermostat, knondiag = knondiag)
    return storage.time,storage.energy
end

# test non cross K integration
ts, energies = my_run(false,false)
@test mean(energies .- energies[1]) < 0.0001

# test cross K integration
ts, energies = my_run(true,false)
@test mean(energies .- energies[1]) < 0.0001

# test thermostat
ts, energies = my_run(true,true)
@test mean(energies[2000] .- energies[1]) < 0.0

# test zero temperature initiazization

function my_run2(knondiag, thermostat)
    storage = storages()
    # Create a model
    sbm = SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 7.5, α = 0.1, No = 60, β = Inf))
    # Create containers starting from a localized state
    ops = CEIDOps(SA[1.0 0.0; 0.0 0.0],sbm)
    # run dynamics
    RunCEID!(sbm,ops,2000,0.001,store!, storage, thermostat = thermostat, knondiag = knondiag)
    return storage.time,storage.energy
end

ts, energies = my_run2(false,false)
@test mean(energies .- energies[1]) < 0.0001


text:
# Some of the code in this file is taken from Registrator.jl (license: MIT)
# https://github.com/JuliaRegistries/Registrator.jl
#
# Some of the code in this file is taken from BinaryBuilder.jl (license: MIT)
# https://github.com/JuliaPackaging/BinaryBuilder.jl

function parse_github_exception(ex::ErrorException)
    msgs = map(strip, split(ex.msg, '\n'))
    d = Dict()
    for m in msgs
        a, b = split(m, ":"; limit=2)
        d[a] = strip(b)
    end
    return d
end

function is_pr_exists_exception(ex)
    d = parse_github_exception(ex)

    if d["Status Code"] == "422" &&
       match(r"A pull request already exists", d["Errors"]) !== nothing
        return true
    end

    return false
end

function _get_token_from_auth(auth::GitHub.OAuth2)
    return auth.token
end

function whoami_username(; auth)
    io = IOBuffer()
    Downloads.download(
        "https://api.github.com/user",
        io;
        headers = Dict(
            "Authorization" => "token $(_get_token_from_auth(auth))",
        ),
    )
    str = String(take!(io))
    json = JSON3.read(str)
    login = json.login::String
    return login
end

function create_or_update_pull_request(repo::GitHub.Repo, params; auth)
    whoami = whoami_username(; auth)
    try
        @info "Attempting to create a new pull request..."
        new_pr = GitHub.create_pull_request(
            repo;
            params=params,
            auth=auth,
        )
        @info "Created new pull request" new_pr
        return new_pr
    catch ex
        @debug("", exception=(ex, catch_backtrace()))
        if is_pr_exists_exception(ex)
            found_prs, _ = GitHub.pull_requests(
                repo;
                auth=auth,
                params=Dict(
                    "state" => "open",
                    "base" => params["base"],
                    "head" => params["head"],
                ),
            )
            @debug "" found_prs
            if length(found_prs) != 1
                @warn "length(found_prs) != 1" length(found_prs) found_prs [x.number for x in found_prs]
            end
            found_pr = found_prs[1]
            @info "There is already an existing pull request: $(found_pr.number)"
            @info "Attempting to update the existing pull request..."
            existing_pr = GitHub.update_pull_request(
                repo,
                found_pr.number;
                auth=auth,
                params=params,
            )
            @info "Updated existing pull request" existing_pr
            return existing_pr
        else
            rethrow(ex)
        end
    end
end

function create_or_get_fork(fork::AbstractString, upstream::GitHub.Repo; auth)
    try
        return GitHub.repo(fork; auth = auth)
    catch ex
        @debug("", exception=(ex, catch_backtrace()))
        return GitHub.create_fork(upstream; auth = auth)
    end
end

function update_fork_branch(fork::GitHub.Repo, upstream::GitHub.Repo, branch_name::AbstractString; auth)
    mktempdir() do temp_dir
        cd(temp_dir) do
            upstream_clone_url = "https://github.com/$(upstream.full_name).git"
            token = auth.token
            fork_clone_url = "https://x-access-token:$(token)@github.com/$(fork.full_name).git"
            run(`git clone $(fork_clone_url) FORK`)
            cd("FORK") do
                run(`git remote add upstream $(upstream_clone_url)`)
                run(`git fetch --all --prune`)
                run(`git checkout -B $(branch_name) upstream/$(branch_name)`)
                run(`git pull upstream $(branch_name)`)
                run(`git reset --hard upstream/$(branch_name)`)
                run(`git push --force origin $(branch_name)`)
            end # cd("FORK")
        end # cd(temp_dir)
    end # mktempdir()
    return nothing
end


text:
import DataFrames
import Random
import StaticArrays
using DifferentialEquations

seasonality(t, β, ε) = (1 + ε * sin(π * t)^6) .* β

dydt(y, z, w, λ, γ, σ, μ, χ) = λ .* ((1 .- w) .+ (1 .- γ) .* (w .- z)) - σ .* y .- μ .* y .+ χ' * y
dzdt(z, λ, μ, χ) = λ .* (1 .- z) .- μ .* z .+ χ' * z
dwdt(w, ssλ, μ, χ) = ssλ .* (1 .- w) .- μ .* w .+ χ' * w

function mantis!(deriv, state, params, t)
     y = @view state[:,:,1];  z = @view state[:,:,2];  w = @view state[:,:,3]
    dy = @view deriv[:,:,1]; dz = @view deriv[:,:,2]; dw = @view deriv[:,:,3]
    # adjust β to address seasonality
    βi = seasonality(t, params.β, params.ε)
    # force of infection in y,z equations
    # λ = βi .* y
    λ = βi .* y
    # force of infection for strains that share epitopes with each strain
    # [i.e. force of infection for each virus in w equation]
    ssλ = hcat([sum(βi[x] .* view(state, :, params.sharedalleles[x,:], 1), dims=2) for x in 1:size(state, 2)]...)
    # calculate the derivatives
    dy .= dydt(y, z, w, λ, params.γ, params.σ, params.μ, params.χ)
    dz .= dzdt(z, λ, params.μ, params.χ)
    dw .= dwdt(w, ssλ, params.μ, params.χ)
end

function runMANTIS(;strainstructure=[2 2], tmax=1000.0, tstep=missing,
                   beta=4.0, gamma=0.72, epsilon=0.0, sigma=1.0, mu=0.0,
                   chi=[-0.02 0.0; 0.0 -0.02], rseed=1, initvals::Any="random")
    # strainstructure=[2 2];tmax=100.0;tstep=missing;beta=40.0;gamma=0.72;
    # epsilon=0.0;sigma=10.0;mu=0.0;chi=[-0.02 0.0; 0.0 -0.02];rseed=1;initvals="random"

    chi = isempty(size(chi)) ? StaticArrays.SMatrix{1,1}(chi) : StaticArrays.SMatrix{size(chi)...}(chi)

    strains = reshape(collect(Iterators.product([1:x for x in strainstructure]...)), :, 1)

    npops = size(chi, 1)
    nstrains = length(strains)

    if typeof(initvals) == String && initvals == "random"
        # Random.seed!(rseed)
        Random.seed!(trunc(Int, rand() * 1000))
        initvals = rand(npops, nstrains) |> x->x./sum(x, dims=2)
    end

    # initcond = StaticArrays.Size((size(initvals)..., 3))([initvals..., initvals..., min.(3*initvals, 1)...])
    initcond = cat(initvals, initvals, min.(3*initvals, 1), dims=3)

    if length(beta) == 1
        beta_resized = fill(beta, npops, nstrains)
    elseif length(beta) == npops
        beta_resized = hcat(fill(reshape(beta, npops, 1), nstrains)...)
    elseif length(beta) == nstrains
        beta_resized = vcat(fill(reshape(beta, 1, nstrains), npops)...)
    elseif length(beta) == npops * nstrains
        beta_resized = reshape(beta, npops, nstrains)
    else error("Invalid beta value") end

    if length(sigma) == 1
        sigma_resized = fill(sigma, npops, nstrains)
    elseif length(sigma) == npops
        sigma_resized = hcat(fill(reshape(sigma, npops, 1), nstrains)...)
    elseif length(sigma) == nstrains
        sigma_resized = vcat(fill(reshape(sigma, 1, nstrains), npops)...)
    elseif length(sigma) == npops * nstrains
        sigma_resized = reshape(sigma, npops, nstrains)
    else error("Invalid sigma value") end

    if length(mu) == 1
        mu_resized = StaticArrays.@SMatrix fill(mu, npops, nstrains)
    elseif length(mu) == npops
        mu_resized = hcat(fill(reshape(mu, npops, 1), nstrains)...)
    elseif length(mu) == nstrains
        mu_resized = vcat(fill(reshape(mu, 1, nstrains), npops)...)
    elseif length(mu) == npops * nstrains
        mu_resized = reshape(mu, npops, nstrains)
    else error("Invalid mu value") end

    if length(gamma) == 1
        gamma_resized = StaticArrays.@SMatrix fill(gamma, npops, nstrains)
    elseif length(gamma) == npops
        gamma_resized = hcat(fill(reshape(gamma, npops, 1), nstrains)...)
    elseif length(gamma) == nstrains
        gamma_resized = vcat(fill(reshape(gamma, 1, nstrains), npops)...)
    elseif length(gamma) == npops * nstrains
        gamma_resized = reshape(gamma, npops, nstrains)
    else error("Invalid gamma value") end

    params = (β = beta_resized,
              γ = gamma_resized,
              ε = epsilon,
              σ = sigma_resized,
              μ = mu_resized,
              χ = chi,
              sharedalleles = hcat([[any(x .== y) for x in strains] for y in strains]...))

    prob = ODEProblem(mantis!, initcond, (0.0, tmax), params)
    if ismissing(tstep)
        sol = solve(prob, save_everystep=false)
        return(sol)
    else
        sol = solve(prob, AutoTsit5(Rosenbrock23()), reltol=1e-12, abstol=1e-12, saveat=tstep)
        tuplestosave = tuples(sol)#[trunc(Int, tmax / tstep * 0.9):end]
        return(Dict("timeseries" => DataFrames.DataFrame(hcat([[reshape(u, (1,prod(size(u))))..., t] for (u,t) in tuplestosave]...)'),
                    "parameters" => params))
    end
end


text:
#!/usr/bin/env julia
# push!(LOAD_PATH, joinpath(@__DIR__, "..", "src"))
using ArgParse
using Mycelia
import Pkg

function parse_commandline()
    settings = ArgParseSettings()
    settings.description = "description"
    settings.version = Pkg.TOML.parsefile(joinpath(pkgdir(Mycelia), "Project.toml"))["version"]
    settings.add_version = true

    @add_arg_table settings begin
       "construct"
           help = "construct a reference graph"
           action = :command
       "assemble"
           help = "assemble a sequence graph"
           action = :command
       "test"
           help = "run tests"
           action = :command
    end

   settings["construct"].description = ""
   @add_arg_table settings["construct"] begin
        "--k"
            help = "kmer size that will be used. Must be prime and < 64"
            arg_type = Int
            required = true
        "--fastx"
            help = "one or more fasta or fastq file(s). Can be gzipped as long as files end in .gz"
            arg_type = String
            required = true
            nargs = '*'
        "--out"
            help = "outpath for the genome graph that will be constructed"
            arg_type = String
            required = true
    end

#     settings["plot"].description = ""
#     @add_arg_table settings["plot"] begin
#         "histogram"
#             help = "plot a histogram of kmer counts"
#             action = :command
#         "rank-frequency"
#             help = "plot the rank-frequency zipfs law relationship of kmer counts"
#             action = :command
#         "gfa"
#             help = "render a sequence graph as an SVG"
#             action = :command
#     end

#     settings["plot"]["histogram"].description = ""
#     @add_arg_table settings["plot"]["histogram"] begin
#         "--histogram"
#             help = "kmer counts histogram file to plot"
#             arg_type = String
#             required = true
#     end

#     settings["plot"]["rank-frequency"].description = ""
#     @add_arg_table settings["plot"]["rank-frequency"] begin
#         "--kmer-counts"
#             help = "kmer count file to plot"
#             arg_type = String
#             required = true
#     end

#     settings["plot"]["gfa"].description = ""
#     @add_arg_table settings["plot"]["gfa"] begin
#         "--gfa"
#             help = "gfa file to plot"
#             arg_type = String
#             required = true
#     end

   settings["assemble"].description = ""
   @add_arg_table settings["assemble"] begin
       "--kmers"
           help = "sorted file of trusted kmers"
           arg_type = String
           required = true
       "--sequences"
           help = """
                  One or more paths to fasta or fastq files.

                  Sequences will be used for determining edges between kmers
                  and depth of coverage for kmers and resulting segments in the
                  GFA graph
                  """
           arg_type = String
           nargs = '*'
       "--colors"
           help = """
                  Paths to files that represent the color groups for the colored debruijn graphs.

                  Each file contains a list of fastq and/or fasta files. Segments will be colored
                  using the tag CL:Z:# in the output GFA file, where # is an interger representing
                  the order of the color groups as they were presented in the command line.
                  """
           arg_type = String
           nargs = '*'
    end

    return parse_args(settings)
end

function main()
    parsed_args = parse_commandline()
    if parsed_args["%COMMAND%"] == "construct"
        println("command construct chosen!")
    elseif parsed_args["%COMMAND%"] == "assemble"
        Mycelia.assemble(parsed_args["assemble"])
    elseif parsed_args["%COMMAND%"] == "test"
        runtests_file = joinpath(@__DIR__, "..", "test", "runtests.jl")
        run(`julia --color yes $runtests_file`)
    else
        println("command")
    end
#     if parsed_args["%COMMAND%"] == "plot"
#         if parsed_args["plot"]["%COMMAND%"] == "gfa"
#             Eisenia.plot_gfa(parsed_args["plot"]["gfa"])
#         elseif parsed_args["plot"]["%COMMAND%"] == "histogram"
#             Eisenia.plot_histogram(parsed_args["plot"]["histogram"])
#         elseif parsed_args["plot"]["%COMMAND%"] == "rank-frequency"
#             Eisenia.plot_rank_frequency(parsed_args["plot"]["rank-frequency"])
#         end
end

main()

max_stars_repo_name: cjprybol/Mycelia

max_issues_repo_name: cjprybol/Mycelia

max_forks_repo_name: cjprybol/Mycelia


text:
import COBRA

function lpproblem(S::AbstractMatrix, b::AbstractVector, c::AbstractVector,
                   lb::AbstractVector, ub::AbstractVector)
    @assert length(b) == size(S,1)
    @assert length(lb) == length(ub) == length(c) == size(S,2)
    osense = -1 # maximize
    csense = fill('=', length(b))
    mets = "M" .* string.(1 : size(S,1))
    rxns = "R" .* string.(1 : size(S,2))
    COBRA.LPproblem(S, b, c, lb, ub, osense, csense, rxns, mets)
end

"""
Converts a MetNet to COBRA.LPProblem.
The MetNet struct is defined in
https://github.com/anna-pa-m/Metabolic-EP
"""
function lpproblem(metnet)
    lpproblem(metnet.S, metnet.b, metnet.c, metnet.lb, metnet.ub)
end


text:
## AlmostBandedMatrix

immutable AlmostBandedMatrix{T} <: AbstractMatrix{T}
    bands::BandedMatrix{T}
    fill::LowRankMatrix{T}
    function AlmostBandedMatrix(bands::BandedMatrix{T},fill::LowRankMatrix{T})
        if size(bands) ≠ size(fill)
            error("Data and fill must be compatible size")
        end
        new(bands,fill)
    end
end

AlmostBandedMatrix(bands::BandedMatrix,fill::LowRankMatrix) =
    AlmostBandedMatrix{promote_type(eltype(bands),eltype(fill))}(bands,fill)

AlmostBandedMatrix{T}(::Type{T},n::Integer,m::Integer,l::Integer,u::Integer,r::Integer) =
    AlmostBandedMatrix(BandedMatrix(T,n,m,l,u),LowRankMatrix(T,n,m,r))

abzeros{T}(::Type{T},n::Integer,m::Integer,l::Integer,u::Integer,r::Integer) =
    AlmostBandedMatrix(bzeros(T,n,m,l,u),lrzeros(T,n,m,r))

Base.size(A::AlmostBandedMatrix) = size(A.bands)
Base.linearindexing{ABM<:AlmostBandedMatrix}(::Type{ABM}) =
    Base.LinearSlow()

function getindex(B::AlmostBandedMatrix,k::Integer,j::Integer)
    if j > k + bandwidth(B.bands,2)
        B.fill[k,j]
    else
        B.bands[k,j]
    end
end

# can only change the bands, not the fill
function Base.setindex!(B::AlmostBandedMatrix,v,k::Integer,j::Integer)
        B.bands[k,j] = v
end

function pad!(B::AlmostBandedMatrix,n::Integer,m::Integer)
    pad!(B.bands,n,m)
    pad!(B.fill,n,m)
    B
end


text:
using Test

using cfgrib

const dir_tests = abspath(joinpath(dirname(pathof(cfgrib)), "..", "test"))
const dir_testfiles = abspath(joinpath(dir_tests, "sample-data"))

@testset "cfgrib.jl" begin
    @testset "Indexing" begin
        include("indexing.jl")
    end

    @testset "CFMessage" begin
        include("cfmessage.jl")
    end

    @testset "Dataset" begin
        include("dataset.jl")
    end
end


text:
module mt001
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using Test
function test()
    a = VecAttrib([1, 2, 4]);
    # @test show(a) != ""
    c = [(1, 2, 3), (4, 1, 3)]
    cc = [SVector{nvertices(T3)}(c[idx]) for idx in 1:length(c)]
    elements = ShapeColl(T3, length(c))
    vrts = ShapeColl(P1, 4)
    t2v = IncRel(elements, vrts, cc)
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    # println(summary(t2e))
    # @show e2v
    @test t2e[1] == [4, -2, 1]
    # t2v = ShapeColl(T3, 6, "t3s")
    # vrts = ShapeColl(P1, 12, "vrts")
    # ir = IncRel(q4s, vrts, cc)
    # @test summary(ir.left) == "q4s = 6 x Q4"
    # @test summary(ir.right) == "vrts = 12 x P1"
    # @test summary(ir) == "(q4s, vrts): q4s = 6 x Q4, vrts = 12 x P1"
    # @show Q4
    # @show q4s
    true
end
end
using .mt001
mt001.test()

module mt002
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using ShellStructureTopo: orient_surface_mesh
using Test
function test()
    a = VecAttrib([1, 2, 4]);
    # @test show(a) != ""
    c = [(1, 2, 3), (4, 1, 3)]
    cc = [SVector{nvertices(T3)}(c[idx]) for idx in 1:length(c)]
    elements = ShapeColl(T3, length(c))
    vrts = ShapeColl(P1, 4)
    t2v = IncRel(elements, vrts, cc)
    t2v.right.attributes["geom"] = VecAttrib([SVector(0.0, 0.0) for _ in 1:4])
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    orientedt2v, orientable = orient_surface_mesh(t2v)
    @test length(unique(attribute(orientedt2v.left, "surfid"))) == 1
    # @show orientedt2v, orientable
    true
end
end
using .mt002
mt002.test()

module mt003
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using ShellStructureTopo: orient_surface_mesh
using Test
function test()
    c = [(1, 2, 3), (4, 3, 1)]
    cc = [SVector{nvertices(T3)}(c[idx]) for idx in 1:length(c)]
    elements = ShapeColl(T3, length(c))
    vrts = ShapeColl(P1, 4)
    t2v = IncRel(elements, vrts, cc)
    t2v.right.attributes["geom"] = VecAttrib([SVector(0.0, 0.0) for _ in 1:4])
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    orientedt2v, orientable = orient_surface_mesh(t2v)
    @test length(unique(attribute(orientedt2v.left, "surfid"))) == 1
    # @show t2v._v, orientedt2v._v, orientable
    true
end
end
using .mt003
mt003.test()

module mt004
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshSteward: T3block
using MeshSteward: vtkwrite
using ShellStructureTopo: orient_surface_mesh
using Test
function test()
    t2v = T3block(1.0, 2.0, 3, 6)
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    orientedt2v, orientable = orient_surface_mesh(t2v)
    @test length(unique(attribute(orientedt2v.left, "surfid"))) == 1
    # @show t2v._v, orientedt2v._v, orientable
    orientedt2v.right.attributes["geom"] = t2v.right.attributes["geom"]
    vtkwrite("mt004", orientedt2v)
    try rm("mt004.vtu"); catch end
    true
end
end
using .mt004
mt004.test()

module mt005
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshSteward: T3block
using MeshSteward: vtkwrite
using MeshSteward: import_ABAQUS, vtkwrite, export_MESH, import_MESH
using ShellStructureTopo: orient_surface_mesh
using Test
function test()
    connectivities = import_ABAQUS(joinpath("../models", "extrusion.inp"))
    @test length(connectivities) == 1
    t2v = connectivities[1]
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    orientedt2v, orientable = orient_surface_mesh(t2v)
    @test length(unique(attribute(orientedt2v.left, "surfid"))) == 3
    orientedt2v.right.attributes["geom"] = t2v.right.attributes["geom"]
    vtkwrite("mt005", orientedt2v)
    try rm("mt005.vtu"); catch end
    true
end
end
using .mt005
mt005.test()

module mt006
using Random
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshSteward: T3block
using MeshSteward: vtkwrite
using MeshSteward: import_ABAQUS, vtkwrite, export_MESH, import_MESH
using ShellStructureTopo: orient_surface_mesh
using Test
function normal(c, geom)
    StaticArrays.cross(geom[c[2]] - geom[c[1]], geom[c[3]] - geom[c[1]])
end
function test()
    connectivities = import_ABAQUS(joinpath("../models", "extrusion.inp"))
    @test length(connectivities) == 1
    t2v = connectivities[1]
    for k in 1:nshapes(t2v.left)
        t2v._v[k] = shuffle!(Vector(t2v._v[k]))
    end
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    orientedt2v, orientable = orient_surface_mesh(t2v)
    @test length(unique(attribute(orientedt2v.left, "surfid"))) == 3
    orientedt2v.right.attributes["geom"] = t2v.right.attributes["geom"]
    geom = orientedt2v.right.attributes["geom"]
    t2v.right.attributes["normals"] = VecAttrib([normal(t2v[i], geom) for i in 1:nshapes(t2v.left)])
    vtkwrite("mt006-inconsistent", t2v, [(name = "normals", )])
    # try rm("mt006-inconsistent.vtu"); catch end
    orientedt2v.right.attributes["normals"] = VecAttrib([normal(orientedt2v[i], geom) for i in 1:nshapes(orientedt2v.left)])
    vtkwrite("mt006-consistent", orientedt2v, [(name = "normals", )])
    # try rm("mt006-consistent.vtu"); catch end
    true
end
end
using .mt006
mt006.test()

module mt007
using Random
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshSteward: T3block
using MeshSteward: vtkwrite
using MeshSteward: import_ABAQUS, vtkwrite, export_MESH, import_MESH
using ShellStructureTopo: orient_surface_mesh, _classify_mesh_edges
using Test
function normal(c, geom)
    StaticArrays.cross(geom[c[2]] - geom[c[1]], geom[c[3]] - geom[c[1]])
end
function test()
    connectivities = import_ABAQUS(joinpath("../models", "extrusion.inp"))
    @test length(connectivities) == 1
    t2v = connectivities[1]
    for k in 1:nshapes(t2v.left)
        t2v._v[k] = shuffle!(Vector(t2v._v[k]))
    end
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    orientedt2v, orientable = orient_surface_mesh(t2v)
    @test length(unique(attribute(orientedt2v.left, "surfid"))) == 3
    orientedt2v.right.attributes["geom"] = t2v.right.attributes["geom"]
    geom = orientedt2v.right.attributes["geom"]
    orientedt2v.right.attributes["normals"] = VecAttrib([normal(orientedt2v[i], geom) for i in 1:nshapes(orientedt2v.left)])
    vtkwrite("mt007-consistent", orientedt2v, [(name = "normals", )])
    # try rm("mt007-consistent.vtu"); catch end
    t2v = orientedt2v
    e2v = ir_skeleton(t2v)
    e2v.right.attributes["geom"] = t2v.right.attributes["geom"]
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    e2t = ir_transpose(t2e)
    e2t = _classify_mesh_edges(orientedt2v, e2t, 30*pi/180)
    isonte = attribute(e2t.left, "isonte")
    c = [e2v[idx] for idx in 1:nshapes(e2t.left) if isonte[idx]]
    lft = ShapeColl(shapedesc(e2v.left), length(c))
    tee2t = IncRel(lft, e2v.right, c)
    vtkwrite("mt007-te", tee2t)
    true
end
end
using .mt007
mt007.test()

module mt008
using Random
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, T3, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, attribute
using MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets
using MeshSteward: T3block
using MeshSteward: vtkwrite
using MeshSteward: import_ABAQUS, vtkwrite, export_MESH, import_MESH
using ShellStructureTopo: orient_surface_mesh, _normal, make_topo_faces
using Test
function test()
    connectivities = import_ABAQUS(joinpath("../models", "extrusion.inp"))
    @test length(connectivities) == 1
    t2v = connectivities[1]
    for k in 1:nshapes(t2v.left)
        t2v._v[k] = shuffle!(Vector(t2v._v[k]))
    end
    e2v = ir_skeleton(t2v)
    # println(summary(e2v))
    t2e = ir_bbyfacets(t2v, e2v)
    t2v, orientable = orient_surface_mesh(t2v)
    @test orientable == true
    @test length(unique(attribute(t2v.left, "surfid"))) == 3
    geom = t2v.right.attributes["geom"]
    t2v.right.attributes["normals"] = VecAttrib([_normal(t2v[i], geom) for i in 1:nshapes(t2v.left)])
    vtkwrite("mt008-normals", t2v, [(name = "normals", )])
    # try rm("mt008-consistent.vtu"); catch end
    t2v = make_topo_faces(t2v)
    @test length(unique(attribute(t2v.left, "tfid"))) == 4
    vtkwrite("mt008-classified", t2v, [(name = "tfid", )])
    true
end
end
using .mt008
mt008.test()


text:
#include("ConsoleUtil.jl")
include("SortingAlgorithms.jl")

data = [1, 10, 4, 1, 10, 9, 9, 7, 7, 6]

println(data)

quicksort!(data)

println(data)


text:
######################## Kernel Input Transfomations ###########################
############################ Length Scale ######################################
# IDEA: using meta-programming, could write constructors:
# EQ(l::Float) = Lengthscale(EQ(), l)
# for all isotropic kernels?
struct Lengthscale{T, K} <: IsotropicKernel{T}
    k::K
    l::T
    function Lengthscale(k::IsotropicKernel, l)
        all(>(0), l) || throw(DomainError("l = $l is non-positive"))
        length(l) == 1 || throw(DimensionMismatch("lengthscale l has to has length 1"))
        S = promote_type(eltype(k), typeof(l))
        l = convert(S, l)
        new{S, typeof(k)}(k, l)
    end
end
@functor Lengthscale
# IDEA: for optmization with Zygote, could let k.l be vector and write sum(r² ./ k.l)
(k::Lengthscale)(r²::Number) = k.k(r² / k.l[1]^2)

########################### Change of Input Norm  ##############################
# apply a different norm to input radius r of an isotropic kernel
# special case of input transformation functional
# NOTE: needs
struct Normed{T, K, N} <: StationaryKernel{T}
    k::K
    n²::N # norm for r call
    function Normed(k::AbstractKernel{T}, n²) where T
        new{T, typeof(k), typeof(n²)}(k, n²)
    end
    function Normed(k, n²)
        T = typeof(k(n²(0.)))
        new{T, typeof(k), typeof(n²)}(k, n²)
    end
end
@functor Normed
# WARNING: input function must be isotropic
(m::Normed)(τ) = m.k(m.n²(τ))
(m::Normed)(x, y) = m(difference(x, y))

# automatic relevance determination with length scale parameters l
function ARD(k, l::AbstractVector{<:Real})
    f(x) = enorm2(Diagonal(inv.(l)), x)
    Normed(k, f)
end
ARD(k, l::Real) = Lengthscale(k, l)
function Energetic(k, A::AbstractMatOrFac{<:Real})
    f(x) = enorm2(A, x)
    Normed(k, f)
end
############################ periodic kernel ###################################
# derived by David MacKay
# input has to be 1D stationary or isotropic
struct Periodic{T, K<:IsotropicKernel{T}} <: StationaryKernel{T}
    k::K
end
@functor Periodic
# squared euclidean distance of x, y in the space (cos(x), sin(x))
# since τ_new^2 = (cos(x) - cos(y))^2 + (sin(x) - sin(y))^2 = 4*sin((x-y)/2)^2 = 4sin(τ/2)^2
# (p::Periodic)(x, y) = p(difference(x, y))
function (p::Periodic)(τ::Number)
    r² = (2sin(π*τ))^2 # since isotropic kernels are called with r²
    p.k(r²) # without scaling, this is 1-periodic
end
parameters(p::Periodic) = parameters(p.k)
nparameters(p::Periodic) = nparameters(p.k)
Base.similar(p::Periodic, θ::AbstractVector) = similar(p.k, θ)

######################## Kernel Input Transfomations ###########################
######################## general linear input scaling ##########################
struct ScaledInputKernel{T, K, UT} <: AbstractKernel{T}
    k::K
    U::UT
end
@functor ScaledInputKernel
ScaledInputKernel(k, U) = ScaledInputKernel{Float64}(k, U)
ScaledInputKernel{T}(k, U) where T = ScaledInputKernel{T, typeof(k), typeof(U)}(k, U)
(S::ScaledInputKernel)(x, y) = S.k(S.U*x, S.U*y)
(S::ScaledInputKernel)(r) = S.k(S.U*r) # for stationary kernels

# if U has quadratic complexity in d, this reduces from O(n^2d^2) to O(nd^2 + n^2d)
# currently, this acceleration is currently only used if S.U is square and not diagonal
function gramian(S::ScaledInputKernel, x::AbstractVector, y::AbstractVector)
    if size(S.U, 1) < size(S.U, 2)
        Gramian(S, x, y)
    else
        x = (S.U,) .* x
        y = x ≡ y ? x : (S.U,) .* y
        gramian(S.k, x, y)
    end
end
# if the scaling can be carried out efficiently, like for Diagonal, no need to pre-compute scaled data
function gramian(S::ScaledInputKernel{<:Real, <:Any, <:Diagonal},  x::AbstractVector, y::AbstractVector)
    Gramian(S, x, y)
end

############################### input warping ##################################
struct Warped{T, K, U} <: AbstractKernel{T}
    k::K
    u::U
end
@functor Warped
Warped(k, u) = Warped{Float64}(k, u)
Warped{T}(k, u) where T = Warped{T, typeof(k), typeof(u)}(k, u)
function Warped{T}(k, U::AbstractMatOrFac) where T
    u(x) = U*x
    Warped{T, typeof(k), typeof(u)}(k, u)
end
(W::Warped)(x, y) = W.k(W.u(x), W.u(y))

# if U has quadratic complexity in d, this reduces from O(n^2d^2) to O(nd^2 + n^2d)
# currently, this acceleration is currently only used if S.U is square and not diagonal
function gramian(W::Warped, x::AbstractVector, y::AbstractVector)
    x = W.u.(x)
    y = x ≡ y ? x : W.u.(y)
    gramian(W.k, x, y)
end
# if the scaling can be carried out efficiently, like for Diagonal, no need to pre-compute scaled data
function gramian(W::Warped{<:Real, <:Any, <:Diagonal}, x::AbstractVector, y::AbstractVector)
    Gramian(W, x, y)
end

############################ Symmetric Kernel ##################################
# make this useable for multi-dimensional inputs!
# in more dimensions, could have more general axis of symmetry
struct SymmetricKernel{T, K<:AbstractKernel} <: AbstractKernel{T}
    k::K # kernel to be symmetrized
    z::T # center
end
@functor SymmetricKernel
SymmetricKernel(k::AbstractKernel{T}) where T = SymmetricKernel(k, zero(T))

# for 1D axis symmetry
function (k::SymmetricKernel)(x, y)
    x -= k.z; y -= k.z;
    (k.k(x, y) + k.k(-x, y))/2
end

######################## Kernel Output Transformations #########################
############################ Scalar Chain Rule #################################
struct Chained{T, F, K} <: AbstractKernel{T}
    f::F
    k::K # requires k to be a kernel function
end
@functor Chained
Chained(f, k::AbstractKernel{T}) where T = Chained{T}(f, k)
Base.:∘(f, k::AbstractKernel) = Chained(f, k)

(C::Chained)(x, y) = C.f(C.k(x, y))
input_trait(C::Chained) = input_trait(C.k)

############################## rescaled kernel #################################
# IDEA: should be called DiagonalRescaling in analogy to the matrix case
# diagonal rescaling of covariance functions
# generalizes multiplying by constant kernel to multiplying by function
struct VerticalRescaling{T, K<:AbstractKernel{T}, F} <: AbstractKernel{T}
    k::K
    f::F
end
@functor VerticalRescaling

(k::VerticalRescaling)(x, y) = k.f(x) * k.k(x, y) * k.f(y)
# preserve structure if k.k is stationary and x, y are regular grids,
# since that introduces Toeplitz structure
function gramian(k::VerticalRescaling, x::AbstractVector, y::AbstractVector)
    # Diagonal(k.a.(x)) * gramian(k.k, x, y) * Diagonal(k.a.(y)))
    Dx = Diagonal(k.f.(x))
    Dy = Diagonal(k.f.(y))
    K = gramian(k.k, x, y)
    return LazyMatrixProduct(Dx, K, Dy)
end

# normalizes an arbitary kernel so that k(x,x) = 1
normalize(k::AbstractKernel) = VerticalRescaling(k, x->1/√k(x, x))


text:
using NaiveBayes
using RDatasets
using StatsBase

# Example 1
iris = dataset("datasets", "iris")

# observations in columns and variables in rows
X = array(iris[:, 1:4])'
p, n = size(X)
# by default species is a PooledDataArray,
y = [species for species in iris[:, 5]]

# how much data use for training
train_frac = 0.9
k = int(floor(train_frac * n))
idxs = randperm(n)
train = idxs[1:k]
test = idxs[k+1:end]

model = GaussianNB(unique(y), p)
fit(model, X[:, train], y[train])

accuracy = countnz(predict(model, X[:,test]) .== y[test]) / countnz(test)
println("Accuracy: $accuracy")

# Example 2
# 3 classes and 100 random data samples with 5 variables.
n_obs = 100
m = GaussianNB([:a, :b, :c], 5)
X = randn(5, n_obs)
y = sample([:a, :b, :c], n_obs)
fit(m, X, y)
accuracy = sum(predict(m, X) .== y) / n_obs
println("Accuracy: $accuracy")


text:
module Day14

using CombinedParsers
using CombinedParsers.Regexp: whitespace_maybe

export day14

function parsemem(str)
    mem = Sequence("mem[", :ind => Numeric(Int), "] = ", :value => Numeric(Int))
    str |> mem
end

function applytovalue(mask, value)
    result = 0
    for (m,v) in zip(mask, value)
        if m == 'X'
            result = result * 2 + v
        elseif m == '0'
            result = result * 2
        else
            result = result * 2 + 1
        end
    end
    result
end

function day14(::Val{:a}, filename)
    mem = Dict{Int,Int}()
    mask = Char[]
    for line in readlines(filename)
        if line[1:4] == "mask"
            mask = collect(line[8:43])
        else
            data = parsemem(line)
            val = reverse(digits(data[:value], base = 2))
            mem[data[:ind]] = applytovalue(mask, vcat(zeros(Int, 36 - length(val)), val))
        end
    end
    values(mem) |> sum
end

function applytoaddress(mask, addr)
    indices = Int[0]
    for (m,v) in zip(mask, addr)
        new_indices = Int[]
        for ind in indices
            if m == '0'
                push!(new_indices, ind * 2 + v)
            elseif m == '1'
                push!(new_indices, ind * 2 + 1)
            else
                push!(new_indices, ind * 2)
                push!(new_indices, ind * 2 + 1)
            end
        end
        indices = copy(new_indices)
    end
    indices
end

function day14(::Val{:b}, filename)
    mem = Dict{Int,Int}()
    mask = Char[]
    for line in readlines(filename)
        if line[1:4] == "mask"
            mask = collect(line[8:43])
        else
            data = parsemem(line)
            addr = reverse(digits(data[:ind], base = 2))
            indices = applytoaddress(mask, vcat(zeros(Int, 36 - length(addr)), addr))
            for ind in indices
                mem[ind] = data[:value]
            end
        end
    end
    values(mem) |> sum
end

end


text:
"""
 holds algorithm independent information about training procedure. e.g. training cycles, what to record, etc
"""
abstract type trainingInfo end

struct training_and_test{A<:AbstractArray, I<:Integer} <: trainingInfo
    training_data::A
    test_data::A
    batch_size::Union{I,Nothing}
    iters::I
    tr_data_length::I
    tst_data_length::I

    function training_and_test(training_data, test_data, batch_size, iters)
        # automatically compute data length
        new{typeof(training_data), typeof(iters)}(training_data, test_data,
                                        batch_size, iters, size(training_data)[end], size(test_data)[end])
    end
end
"""
    trainingDataOnly(training_data, batch_size, iters)
    trainingDataOnly(training_data, iters)
without batch_size input it takes full batches
"""
trainingDataOnly(tr_data::T, batch_size, iters) where T<:AbstractArray = training_and_test(tr_data, tr_data, batch_size, iters)
trainingDataOnly(tr_data::T, iters) where T<:AbstractArray = trainingDataOnly(tr_data, batch_size, iters)

"""
    get_training_batch(tr::trainingInfo)
get a view of a tr.batch_size'd, randomly chosen subset of tr.training_data
"""
function get_training_batch(tr::trainingInfo)
    function get_the_batch(batch_size::Nothing)
        return @view tr.training_data[:,:]
    end

    function get_the_batch(batch_size::Integer)
        return @view tr.training_data[:,rand(1:tr.tr_data_length, tr.batch_size)]
    end

    get_the_batch(tr.batch_size)
end

"""
    get_training_batch!(existing_data, tr::trainingInfo)

mutating version of get_training_batch. mutates existing_data. recommended over get_training_batch for speed where possible
"""
function get_training_batch!(existing_data, tr::trainingInfo)
    function get_the_batch(batch_size::Nothing)
        # do nothing
    end

    function get_the_batch(batch_size::Integer)
        return @view tr.training_data[:,rand(1:tr.tr_data_length, tr.batch_size)]
    end

    get_the_batch(tr.batch_size)
end

"""
    get_test_data(tr::training_and_test)
get a view of the test data in tr
"""
function get_test_batch(tr::training_and_test)
    return @view tr.test_data[:,:]
end


text:
"""
    stochastic.homodyne_carmichael(H0, C, theta)

Helper function that defines the functions needed to compute homodyne detection
trajectories according to Carmichael with `stochastic.schroedinger_dynamic`.

# Arguments
* `H0`: The deterministic, time-independent system Hamiltonian.
* `C`: Collapse operator (or vector of operators) of the detected output channel(s).
* `theta`: The phase difference between the local oscillator and the signal field.
    Defines the operator of the measured quadrature as
    ``X_\\theta = C e^{-i\\theta} + C^\\dagger e^{i\\theta}``. Needs to be a
    vector of the same length as `C` if `C` is a vector.
* `normalize_expect=true`: Specifiy whether or not to normalize the state vector
    when the expectation value in the nonlinear term is calculated. NOTE:
    should only be set to `false` if the state is guaranteed to be normalized,
    e.g. by setting `normalize_state=true` in `stochastic.schroedinger_dynamic`.

Returns `(fdeterm, fstoch)`, where `fdeterm(t, psi) -> H` and
`fstoch(t, psi) -> Hs` are functions returning the deterministic and stochastic
part of the Hamiltonian required for calling `stochastic.schroedinger_dynamic`.

The deterministic and stochastic parts of the Hamiltonian are constructed as

```math
H_{det} = H_0 + H_{nl},
```

where

```math
H_{nl} = iCe^{-i\\theta} \\langle X_\\theta \\rangle - \\frac{i}{2} C^\\dagger C,
```

and

```math
H_s = iCe^{-i\\theta}.
```
"""
function homodyne_carmichael(H0::AbstractOperator, C::Vector{T}, theta::Vector{R};
            normalize_expect::Bool=true) where {T <: AbstractOperator, R <: Real}
    n = length(C)
    @assert n == length(theta)
    Hs = 1.0im*C .* exp.(-1.0im .* theta)
    X = C .* exp.(-1.0im .* theta) + dagger.(C) .* exp.(1.0im .* theta)
    CdagC = -0.5im .* dagger.(C) .* C

    fstoch(t, psi::StateVector) = Hs
    if normalize_expect
        function H_nl_n(psi::StateVector)
            psi_n = normalize(psi)
            sum(expect(X[i], psi_n)*Hs[i] + CdagC[i] for i=1:n)
        end
        fdeterm_n(t, psi::StateVector) = H0 + H_nl_n(psi)
        return fdeterm_n, fstoch
    else
        H_nl_un(psi::StateVector) =
            sum(expect(X[i], psi)*Hs[i] + CdagC[i] for i=1:n)
        fdeterm_un(t, psi::StateVector) = H0 + H_nl_un(psi)
        return fdeterm_un, fstoch
    end
end
homodyne_carmichael(H0::AbstractOperator, C::AbstractOperator, theta::Real; kwargs...) =
    homodyne_carmichael(H0, [C], [theta]; kwargs...)


text:
"""
LatentND <: AbstractLatent

# Description
A mutable describing a multivariate latent variable.

# Fields
- **`names::Vector{String}`**
- **`val::Vector{Float64}`**
- **`bounds::Vector{Vector{Float64}}`**
- **`prior::Distributions.MultivariateDistribution`**
- **`posterior::Distributions.MultivariateDistribution`**
- **`chain::Vector{Vector{Float64}}`**
- **`expected_information::Matrix{Float64}`**

# Factories
LatentND(names, val, bounds, prior, posterior, chain, expected_information) = new(names, val, bounds, prior, posterior, chain, expected_information)
Creates a new N-dimensional latent variable with custom fields.

# Random Initializers
LatentND()

Randomly generates a value for the N-dimensional latent variable and assigns a default standardized Gaussian prior and posterior
LatentND(dist, bounds)

Randomly generates a value for the N-dimensional latent variable and assigns a custom multivariate distribution to prior and posterior with specific bounds.
"""
mutable struct LatentND <: AbstractLatent
    names::Vector{String}
    val::Vector{Float64}
    bounds::Vector{Vector{Float64}}
    prior::Distributions.MultivariateDistribution
    posterior::Distributions.MultivariateDistribution
    chain::Vector{Vector{Float64}}
    expected_information::Matrix{Float64}

    function LatentND(names, val, bounds, prior, posterior, chain, expected_information)
        new(names, val, bounds, prior, posterior, chain, expected_information)
    end

    # Random Initializers

    function LatentND(
        dist::Distributions.MultivariateDistribution,
        bounds::Vector{Vector{Float64}},
    )
        N = size(bounds, 1)
        val = truncate_rand(dist, bounds)
        new(string.("L_", collect(1:N)), val, bounds, dist, dist, [val], LinearAlgebra.I(N))
    end
end


text:
using FastMarching
using Images
using FileIO

function maze()
  Float64.(channelview(img))
  img = load(joinpath(Pkg.dir("FastMarching"),"examples/images/maze.png"))
end

maze()


text:
addprocs(8)
@everywhere import MultipleShootingPEM
@everywhere ms = MultipleShootingPEM
@everywhere import ParallelTrainingNN
@everywhere nn = ParallelTrainingNN
@everywhere using CSV
@everywhere using JLD2

# Import data
@everywhere df = CSV.read("examples/processed_data/training_set.csv");
@everywhere ti, ui, yi = Vector(df[:time]), Vector(df[:input]), Vector(df[:output]);
@everywhere Ts = (ti[2] - ti[1])/3600;
@everywhere Ni = length(yi)

@everywhere identification_data = nn.IdData(yi, ui, Ts);
@everywhere df = CSV.read("examples/processed_data/test_set.csv");
@everywhere tv, uv, yv = Vector(df[:time]), Vector(df[:input]), Vector(df[:output]);
@everywhere Ts = (tv[2] - tv[1])/3600;
@everywhere validation_data = nn.IdData(yv, uv, Ts);

# Define neural network model
@everywhere mdl = nn.FeedforwardNetwork(2, 1, [10])
@everywhere yterms = [[1]];
@everywhere uterms = [[1]];
@everywhere mdl = nn.learn_normalization(mdl, yterms, uterms, identification_data);
@everywhere dx_aux = Matrix{Float64}(1, 2)
@everywhere function f(y, dx, dθ, x, k, θ)
    x_aux = [x[1], ui[k]]
    nn.evaluate!(mdl, x_aux, θ, y; dx=dx_aux, dΘ=dθ)
    dx[1, 1] = dx_aux[1, 1]
    return
end
@everywhere function g(y, dx, dθ, x, k, θ)
    y .= x
    dx .= 1
    dθ .= 0
    return
end

@everywhere function compute_solution(seed, shoot_len::Int)
    println("shoot_len = " * string(shoot_len))
    srand(seed)
    θ0 = nn.initial_guess(mdl)
    k0 = 1
    k0_list = collect(k0:shoot_len:k0+Ni-1)
    yi_aux = [[element] for element in yi]
    x0_list = yi_aux[k0_list]
    opt = ms.OptimizationProblem(f, g, x0_list, yi_aux, k0_list, θ0)

    res = ms.solve(opt, options=Dict("gtol" => 1e-12,
                                     "xtol" => 1e-12,
                                     "maxiter" => 2000,
                                     "initial_constr_penalty" => 0.01))
    delete!(res, "jac")
    JLD2.@save "solutions/sol"*string(seed)*"_"*string(shoot_len)*".jld2" res
end

function compute_solution(seed, shoot_len::String)
    tic()
    if shoot_len == "ARX"
        dynfit, opt, e = nn.narx(mdl, yterms, uterms, identification_data,
                                 maxIter=2000);
    elseif shoot_len == "OE"
        dynfit, opt, e = nn.noe(mdl, yterms, uterms, identification_data,
                                use_extended=false, maxIter=2000);
    end
    execution_time = toc()
    res = Dict("x" => dynfit, "execution_time" => execution_time)
    JLD2.@save "solutions/sol"*string(seed)*"_"*string(shoot_len)*".jld2" res
end

# Multiple shooting error estimation
@everywhere function compute_solutions(seed)
    println("seed = " * string(seed))
    shoot_len_list = [3, "ARX", "OE"]
    for shoot_len in shoot_len_list
        compute_solution(seed, shoot_len)
    end
end

if ~Base.Filesystem.isdir("solutions")
    Base.Filesystem.mkdir("solutions")
end

#@everywhere compute_solution(shoot_len) = compute_solution(1, shoot_len)
#pmap(compute_solution, [Ni, 50, 20, 10, 5, 3])
pmap(compute_solutions_original_paper, 1:100)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

################################################################################
# Dataset related utilities
################################################################################
function get_data_dir()
  data_dir = joinpath(@__DIR__, "..", "data")
  mkpath(data_dir)
  data_dir
end

function get_mnist_ubyte()
  data_dir  = get_data_dir()
  mnist_dir = joinpath(data_dir, "mnist")
  mkpath(mnist_dir)
  filenames = Dict(:train_data  => "train-images-idx3-ubyte",
                   :train_label => "train-labels-idx1-ubyte",
                   :test_data   => "t10k-images-idx3-ubyte",
                   :test_label  => "t10k-labels-idx1-ubyte")
  filenames = Dict((x[1] => joinpath(mnist_dir, x[2]) for x ∈ pairs(filenames)))
  if !all(isfile, values(filenames))
    cd(mnist_dir) do
      mnist_dir = download("http://data.mxnet.io/mxnet/data/mnist.zip", "mnist.zip")
        try
          run(`unzip -u $mnist_dir`)
        catch
          try
            run(pipe(`7z x $mnist_dir`,stdout = devnull))
          catch
            error("Extraction Failed:No extraction program found in path")
          end
      end
    end
  end
  return filenames
end

function get_cifar10()
  data_dir    = get_data_dir()
  cifar10_dir = joinpath(data_dir, "cifar10")
  mkpath(cifar10_dir)
  filenames = Dict(:train => "cifar/train.rec", :test => "cifar/test.rec")
  filenames = Dict(map((x) -> x[1] => joinpath(cifar10_dir, x[2]), filenames))
  if !all(isfile, values(filenames))
    cd(cifar10_dir) do
      download("http://data.mxnet.io/mxnet/data/cifar10.zip", "cifar10.zip")
        try
          run(`unzip -u cifar10.zip`)
        catch
          try
            run(pipeline(`7z x cifar10.zip`, stdout = devnull))
          catch
            error("Extraction Failed:No extraction program found in path")
          end
      end
    end
  end

  filenames[:mean] = joinpath(cifar10_dir, "cifar/cifar_mean.bin")
  return filenames
end

################################################################################
# Internal Utilities
################################################################################
function _get_libmx_op_names()
  n = Ref{MX_uint}(0)
  names = Ref{char_pp}(0)

  @mxcall(:MXListAllOpNames, (Ref{MX_uint}, Ref{char_pp}), n, names)

  names = unsafe_wrap(Array, names[], n[])
  return [unsafe_string(x) for x in names]
end
function _get_libmx_op_handle(name :: String)
  handle = Ref{MX_handle}(0)
  @mxcall(:NNGetOpHandle, (char_p, Ref{MX_handle}), name, handle)
  return MX_OpHandle(handle[])
end

# We keep a cache and retrieve the address everytime
# we run Julia, instead of pre-compiling with macro,
# because the actual handle might change in different
# runs
const _libmx_op_cache = Dict{String, MX_OpHandle}()
function _get_cached_libmx_op_handle(name :: String)
  if !haskey(_libmx_op_cache, name)
    handle = _get_libmx_op_handle(name)
    _libmx_op_cache[name] = handle
    return handle
  else
    return _libmx_op_cache[name]
  end
end

function _get_libmx_op_description(name::String, handle::MX_OpHandle)
  # get operator information (human readable)
  ref_real_name = Ref{char_p}(0)
  ref_desc = Ref{char_p}(0)
  ref_narg = Ref{MX_uint}(0)

  ref_arg_names = Ref{char_pp}(0)
  ref_arg_types = Ref{char_pp}(0)
  ref_arg_descs = Ref{char_pp}(0)

  ref_key_narg  = Ref{char_p}(0)
  ref_ret_type  = Ref{char_p}(0)

  @mxcall(:MXSymbolGetAtomicSymbolInfo,
         (MX_handle, Ref{char_p}, Ref{char_p}, Ref{MX_uint}, Ref{char_pp},
          Ref{char_pp}, Ref{char_pp}, Ref{char_p}, Ref{char_p}),
          handle, ref_real_name, ref_desc, ref_narg, ref_arg_names,
          ref_arg_types, ref_arg_descs, ref_key_narg, ref_ret_type)

  real_name = unsafe_string(ref_real_name[])
  signature = _format_signature(Int(ref_narg[]), ref_arg_names)
  desc = "    " * name * "(" * signature * ")\n\n"
  if real_name != name
    desc *= name * " is an alias of " * real_name * ".\n\n"
  end

  key_narg = unsafe_string(ref_key_narg[])
  if key_narg != ""
    desc *= "**Note**: " * name * " takes variable number of positional inputs. "
    desc *= "So instead of calling as $name([x, y, z], $key_narg=3), "
    desc *= "one should call via $name(x, y, z), and $key_narg will be "
    desc *= "determined automatically.\n\n"
  end

  desc *= unsafe_string(ref_desc[]) * "\n\n"
  desc *= "# Arguments\n"
  desc *= _format_docstring(Int(ref_narg[]), ref_arg_names, ref_arg_types, ref_arg_descs)
  return desc, key_narg
end

_format_typestring(s::String) = replace(s, r"\bSymbol\b" => "SymbolicNode")

function _format_docstring(narg::Int, arg_names::Ref{char_pp}, arg_types::Ref{char_pp}, arg_descs::Ref{char_pp}, remove_dup::Bool=true)
  param_keys = Set{String}()

  arg_names  = unsafe_wrap(Array, arg_names[], narg)
  arg_types  = unsafe_wrap(Array, arg_types[], narg)
  arg_descs  = unsafe_wrap(Array, arg_descs[], narg)
  docstrings = String[]

  for i = 1:narg
    arg_name = unsafe_string(arg_names[i])
    if arg_name ∈ param_keys && remove_dup
      continue
    end
    push!(param_keys, arg_name)

    arg_type = _format_typestring(unsafe_string(arg_types[i]))
    arg_desc = unsafe_string(arg_descs[i])
    push!(docstrings, "* `$arg_name::$arg_type`: $arg_desc\n")
  end
  return join(docstrings, "\n")
end

function _format_signature(narg::Int, arg_names::Ref{char_pp})
  arg_names  = unsafe_wrap(Array, arg_names[], narg)

  return join([unsafe_string(name) for name in arg_names] , ", ")
end

"""
Extract the line of `Defined in ...`

julia> mx._getdocdefine("sgd_update")
"Defined in src/operator/optimizer_op.cc:L53"
```
"""
function _getdocdefine(name::String)
  op = _get_libmx_op_handle(name)
  str = _get_libmx_op_description(name, op)[1]
  lines = split(str, '\n')
  for m ∈ match.(Ref(r"^Defined in .*$"), lines)
    m != nothing && return m.match
  end
  ""
end

"""
libmxnet operators signature checker.

C/Python have different convernsion of accessing array. Those languages
handle arrays in row-major and zero-indexing which differs from Julia's
colume-major and 1-indexing.

This function scans the docstrings of NDArray's APIs,
filter out the signature which contain `axis`, `axes`, `keepdims` and `shape`
as its function argument.

We invoks this checker in Travis CI build and pop up the warning message
if the functions does not get manually mapped
(imply it's dimension refering may looks weird).

If you found any warning in Travis CI build, please open an issue on GitHub.
"""
function _sig_checker()
  names = filter(n -> ∉(lowercase(n), _op_import_bl), _get_libmx_op_names())
  foreach(names) do name
    op_handle = _get_libmx_op_handle(name)

    desc, key_narg = _get_libmx_op_description(name, op_handle)
    _sig = desc |> s -> split(s, '\n') |> first |> strip
    _m = match(r"(axis|axes|keepdims|shape)", _sig)

    if _m === nothing
      return
    end

    @warn(_sig)

  end
end

"""
Get first position argument from function sig
"""
function _firstarg(sig::Expr)
  if sig.head ∈ (:where, :(::))
    _firstarg(sig.args[1])
  elseif sig.head == :call
    i = if sig.args[2] isa Expr && sig.args[2].head == :parameters
      # there are some keyward arguments locate at args[2]
      3
    elseif sig.args[1] === :broadcast_
      # case of broadcasting, skip the first arg `::typeof(...)`
      3
    else
      2
    end
    _firstarg(sig.args[i])
  end
end

_firstarg(s::Symbol) = s

const _import_map = Dict{Symbol,Union{Missing,Module}}(
  :diag    => LinearAlgebra,
  :dot     => LinearAlgebra,
  :norm    => LinearAlgebra,

  :shuffle => Random,

  :mean    => Statistics,

  :gamma   => missing,
)

function _import_expr(func_name::Symbol)
  mod = get(_import_map, func_name, Base)
  isdefined(mod, func_name) ? :(import $(Symbol(mod)): $func_name) : :()
end

max_stars_repo_name: bigbigzxl/incubator-mxnet

max_issues_repo_name: bigbigzxl/incubator-mxnet

max_forks_repo_name: bigbigzxl/incubator-mxnet


text:
using TexTables
function merge_dicts(ad::Array{Dict{String,Any},1})
    new_dict = Dict()
    for (i,d) in enumerate(ad)
        for (key,value) in d
            new_dict["$(i)_$key"] = value
        end
    end
    return new_dict
end

function expand_data(xQ::Array{Float64},xP::Array{Float64};order::Int64=2)
    for i = 2:order
        xQ = [xQ xQ.^i]
        xP = [xP xP.^i]
    end
    return [xQ xP]
end

function get_data(data::Dict;return_df::Bool=false,return_params::Bool=true)
    tasks = collect(keys(data))
    task_number = []
    task_order = []
    try
        task_number = parse.(Int64,collect(tasks))
    catch
        task_number = []
    end

    xQ = zeros(length(tasks))
    xP = zeros(length(tasks))
    N = zeros(length(tasks))
    tp = zeros(length(tasks))
    outcome = Array{Bool}(undef,length(tasks))
    for (i,t) in enumerate(tasks)
        xQ[i] = data[t]["xQ"]
        xP[i] = data[t]["xP"]
        #  println("task: $t, xQ: $(xQ[i]), xP: $(xP[i])")
        if return_params
            N[i] = data[t]["num_nodes"]
            tp[i] = data[t]["trans_prob"]
        end
        data[t]["outcome"] == "fail" ? outcome[i] = false : outcome[i] = true
    end

    Xp = [xQ[outcome] xP[outcome]]
    Xn = [xQ[.!outcome] xP[.!outcome]]

    (Xp_m,Xp_c) = mean_and_cov(Xp)
    (Xn_m,Xn_c) = mean_and_cov(Xn)

    if !return_df
        if return_params
            return xQ, xP, outcome, Xp, Xn, Xp_m, Xp_c, Xn_m, Xn_c, N, tp, task_number
        else
            return xQ, xP, outcome, Xp, Xn, Xp_m, Xp_c, Xn_m, Xn_c, task_number
        end
    else
        df = DataFrame(xQ=xQ,xP=xP,outcome=outcome,N=N,tp=tp,task_number=task_number)
        return xQ, xP, df
    end
end

function prop_in_quadrant(xQ,xP)
    num_in_q = Array{Int64,1}()
    for i = 1:4
        if i == 1
            n = sum([q>1.0 && p>0.0 for (q,p) in zip(xQ,xP)])
        elseif i == 2
            n = sum([q<1.0 && p>0.0 for (q,p) in zip(xQ,xP)])
        elseif i == 3
            n = sum([q<1.0 && p<0.0 for (q,p) in zip(xQ,xP)])
        elseif i == 4
            n = sum([q>1.0 && p<0.0 for (q,p) in zip(xQ,xP)])
        end
        push!(num_in_q,n)
    end
    return num_in_q./length(xQ)
end

function num_from_prop(props::Array{Float64},num::Int64;max_num::Int64=-1,min_num::Int64=-1)
    num_float = props.*num

    if max_num != -1
        num_float[num_float.>max_num] .= max_num
    end
    if min_num != -1
        num_float[num_float.<min_num] .= min_num
    end

    return Int.(ceil.(num_float))
end

function choose_equal_portion_set(xQ,xP,outcomes,num::Array{Int64}=[];max_num::Int64=16,min_num::Int64=8)
    if length(num) == 1
        #default to equal examples per quadrant
        props = prop_in_quadrant(xQ,xP)
        num_per_quadrant = num_from_prop(props,num[1],max_num=max_num,min_num=min_num)
        println("Experiment Trial Selection:")
        println("target N: $(sum(num_per_quadrant)), p_num: $num_per_quadrant, props: $props")
    elseif length(num) == 4
        num_per_quadrant = num
    else
        error("need 4 quadrants")
    end
    task_set = Array{Int64,1}()
    q_cnt = Array{Int64,1}()
    s_cnt = Array{Int64,1}()
    sf_ratio = Array{Float64,1}()
    for i = 1:4
        if i == 1
            quad_idx = [q>=1.0 && p>=0.0 for (q,p) in zip(xQ,xP)]
        elseif i == 2
            quad_idx = [q<1.0 && p>=0.0 for (q,p) in zip(xQ,xP)]
        elseif i == 3
            quad_idx = [q<1.0 && p<0.0 for (q,p) in zip(xQ,xP)]
        elseif i == 4
            quad_idx = [q>=1.0 && p<0.0 for (q,p) in zip(xQ,xP)]
        end

        q_outcomes = outcomes[quad_idx]
        q_s_ratio = sum(q_outcomes)/length(q_outcomes)
        q_f_ratio = sum(.!q_outcomes)/length(q_outcomes)
        q_sf_ratio = q_s_ratio/q_f_ratio

        q_subset = collect(1:length(xQ))[quad_idx]

        num_pos = Int(round(num_per_quadrant[i]*q_s_ratio))
        num_neg = Int(round(num_per_quadrant[i]*q_f_ratio))

        println(length(q_outcomes))

        pos_subset = q_subset[q_outcomes][rand(1:sum(q_outcomes),num_pos)]
        neg_subset = q_subset[.!q_outcomes][rand(1:sum(.!q_outcomes),num_neg)]
        println("size pos: $(size(pos_subset)), neg: $(size(neg_subset))")

        subset = [pos_subset; neg_subset]

        println("Q1 subset: $(sort(subset))")

        if length(subset) < num_per_quadrant[i]
            subset = subset
        else
            subset = subset[1:num_per_quadrant[i]]
        end
        #  println("q$i size: $(sum(quad_idx))")
        push!(q_cnt,sum(quad_idx))
        push!(s_cnt,length(subset))
        push!(task_set,subset...)
        push!(sf_ratio,q_sf_ratio)
    end

    println("actual N: $(length(task_set)), p_num_actual: $s_cnt")
    println("sf ratios: $sf_ratio")
    println("total q: $q_cnt")
    println("total q ratio: $(q_cnt./sum(q_cnt))")
    println("total Q12: $(sum(s_cnt[1:2])), total Q14: $(sum(s_cnt[[1,4]]))")
    return task_set
end

function choose_equal_portion_set_lhc(xQ,xP,outcomes,num::Array{Int64}=[];max_num::Int64=16,min_num::Int64=8)
    if length(num) == 1
        #default to equal examples per quadrant
        props = prop_in_quadrant(xQ,xP)
        num_per_quadrant = num_from_prop(props,num[1],max_num=max_num,min_num=min_num)
        println("Experiment Trial Selection:")
        println("target N: $(sum(num_per_quadrant)), p_num: $num_per_quadrant, props: $props")
    elseif length(num) == 4
        num_per_quadrant = num
    else
        error("need 4 quadrants")
    end
    task_set = Array{Int64,1}()
    q_cnt = Array{Int64,1}()
    s_cnt = Array{Int64,1}()
    s_ratio = Array{Float64,1}()
    f_ratio = Array{Float64,1}()
    pos_subset_ary = Array{Any,1}()
    neg_subset_ary = Array{Any,1}()
    q_idx = Array{Any,1}()
    for i = 1:4
        if i == 1
            quad_idx = [q>=1.0 && p>=0.0 for (q,p) in zip(xQ,xP)]
        elseif i == 2
            quad_idx = [q<1.0 && p>=0.0 for (q,p) in zip(xQ,xP)]
        elseif i == 3
            quad_idx = [q<1.0 && p<0.0 for (q,p) in zip(xQ,xP)]
        elseif i == 4
            quad_idx = [q>=1.0 && p<0.0 for (q,p) in zip(xQ,xP)]
        end

        q_outcomes = outcomes[quad_idx]
        q_s_ratio = sum(q_outcomes)/length(q_outcomes)
        q_f_ratio = sum(.!q_outcomes)/length(q_outcomes)
        q_sf_ratio = q_s_ratio/q_f_ratio

        q_subset = collect(1:length(xQ))[quad_idx]

        num_pos = Int(round(num_per_quadrant[i]*q_s_ratio))
        num_neg = Int(round(num_per_quadrant[i]*q_f_ratio))

        #  @info "Q$i" num_pos num_neg
        pos_data = [xQ xP][q_subset[q_outcomes],:]
        neg_data = [xQ xP][q_subset[.!q_outcomes],:]
        pos_subset = q_subset[q_outcomes][bin_sample(pos_data,num_pos)]
        neg_subset = q_subset[.!q_outcomes][bin_sample(neg_data,num_neg)]

        #  @info "Actual Size" size(pos_subset) size(neg_subset)

        subset = [pos_subset; neg_subset]

        println("Q$i subset: $pos_subset, $neg_subset")

        #  println("q$i size: $(sum(quad_idx))")
        push!(q_cnt,sum(quad_idx))
        push!(s_cnt,length(subset))
        push!(task_set,subset...)
        push!(s_ratio,q_s_ratio)
        push!(f_ratio,q_f_ratio)
        push!(pos_subset_ary,pos_subset)
        push!(neg_subset_ary,neg_subset)
        push!(q_idx,quad_idx)
    end

    println("actual N: $(length(task_set)), p_num_actual: $s_cnt")
    println("total q: $q_cnt")
    println("total q: $(q_cnt./sum(q_cnt))")
    println("total Q12: $(sum(s_cnt[1:2])), total Q14: $(sum(s_cnt[[1,4]]))")

    function return_row(row_name,data_ary,first_val,col_prefix::String="Quad.")
        row = TableCol("Total",[row_name],[first_val])
        for i = length(data_ary):-1:1
            row = hcat(TableCol("$col_prefix $i",[row_name],[data_ary[i]]),row)
        end
        return row
    end
    #####
    # save table with the statistics of the data set
    rn_row = return_row("Original N",sum.(q_idx),length(outcomes))
    rn_prop = sum.(q_idx)./length.(q_idx)
    rn_prop_row = return_row("Original N Proportion",rn_prop,sum(rn_prop))

    sf = [sum(outcomes[q])/sum(.!outcomes[q]) for q in q_idx]
    r_sf_ratio = return_row("Original S/F Ratio",sf,sum(outcomes)/length(outcomes))

    n_row = return_row("N",s_cnt,sum(s_cnt))

    p_row = return_row("N Proportion",s_cnt./sum(s_cnt),sum(s_cnt./sum(s_cnt)))

    tot_s = sum(length.(pos_subset_ary))
    s_row = return_row("Successes",length.(pos_subset_ary),tot_s)

    tot_f = sum(length.(neg_subset_ary))
    f_row = return_row("Failures",length.(neg_subset_ary),tot_f)

    sf_ratio = return_row("S/F Ratio",length.(pos_subset_ary)./length.(neg_subset_ary),tot_s/tot_f)

    q_table = vcat(rn_row,rn_prop_row,r_sf_ratio,n_row,p_row,s_row,f_row,sf_ratio)
    display(q_table)
    write_tex("q_table.tex",q_table)
    #####

    return task_set
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
#using Pkg
#cd("c:/Users/francis.smart.ctr/GitDir/AbstractLogicJL")
#Pkg.activate(".")

# Global set variables Ω, Υ hold
integer(x::AbstractString) = parse(Int, strip(x))
Base.range(x::AbstractString) = range(integer(match(r"^[0-9]+", x).match),
                               stop = integer(match(r"[0-9]+$", x).match))
ABoccursin(y::Symbol) = any( [ y ∈ keys(Ω[i]) for i in 1:length(Ω) ] )
ABoccursin(x::Hotcomb, y::Symbol) = y ∈ keys(x)

Base.range(Ω::Hotcomb, ℧) = [Symbol(keys(Ω)[i])=> sort(unique(Ω[℧][:,i])) for i in 1:size(Ω)[2]]

#command = "a,b,c,d ∈ 1:5"
#Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]"]); Ω[℧]

function ABparse(commands::Array{String,1}; Ω::Hotcomb = Hotcomb(0), ℧::AbstractArray{Bool,1} = Bool[0])
  println("")

  for command in commands
      print(command)
      Ω, ℧ = ABparse(command,Ω,℧)

      feasibleoutcomes = (length(℧)>0) ? sum(℧) : 0
      filler = repeat("\t", max(1, 3-Integer(round(length(command)/18))))

      (feasibleoutcomes == 0)  && (check = "X")
      (feasibleoutcomes  > 1)  && (check = "✓")
      (feasibleoutcomes == 1)  && (check = "✓✓")

      println(" $filler feasible outcomes $feasibleoutcomes $check")

  end
  (Ω, ℧)
end

function ABparse(command::String,  Ω::Hotcomb, ℧::AbstractArray{Bool,1})
  # A vector of non-standard operators to ignore
  exclusionlist = ["\bin\b"]

  if occursin(r"∈|\bin\b", command)
      (Ω,℧) = ABassign(command)
      return (Ω,℧)
  end

  # Check for the existance of any symbols in Ω
  varcheck = eachmatch(r"[a-zA-Z][0-9a-zA-Z_.]*", command)

  # Checks if any of the variables does not exist in Ω
  for S in [Symbol(s.match) for s in varcheck if !(s.match ∈ exclusionlist)]
      if (occursin("{{", string(S))) && (!ABoccursin(Ω, S))
          throw("In {$command} variable {:$S} not found in Ω")
      end
  end

  if occursin(r"( |\b)([><=|!+\\-\^\\&]{1,4})(\b| )", command)
      return SuperSuperOperatorEval(command,Ω,℧)
  end
  println("Warning! { $command } not interpretted")
end

function ABassign(command::String)
  vars, valsin = strip.(split(command, r"∈|\bin\b"))
  varsVect = split(vars, ",") .|> strip
  vals0 = split(replace(valsin, r"\[|\]" => ""), ",")

  (length(vals0) > 1) && (vals =  [(occursin(r"^[0-9]+$", i) && integer(i)) for i in vals0])
  (length(vals0) == 1 && occursin(r"^[0-9]+:[0-9]+$", vals0[1])) && (vals = range(vals0[1]))

  outset = (; zip([Symbol(i) for i in varsVect], fill(vals, length(vars)))...)

  Ω  = Hotcomb(outset)
  ℧ = fill(true, size(Ω)[1])

  (Ω ,℧)
end

function grab(argument::AbstractString, Ω::Hotcomb, ℧::AbstractArray{Bool,1}; command = "")
  matcher = r"^([a-zA-z][a-zA-z0-9]*)*([0-9]+)*([+\-*/])*([a-zA-z][a-zA-z0-9]*)*([0-9]+)*$"

  m = match(matcher, argument)
  nvar = 5-sum([i === nothing for i in m.captures])
  (nvar==0) && throw("Argument $argument could not be parsed in $command")

  v1, n1, o1, v2, n2 = m.captures

  !(v1 === nothing) && (left  = Ω[℧, Symbol(v1)])
  !(n1 === nothing) && (left  = fill(integer(n1), length(℧)))

  (nvar==1) && return left

  !(v2 === nothing) && (right = Ω[℧, Symbol(v2)])
  !(n2 === nothing) && (right = fill(integer(n2), length(left)))

  (o1 == "+") && return left .+ right
  (o1 == "-") && return left .- right
  (o1 == "/") && return left ./ right
  (o1 == "*") && return left .* right
end

function OperatorSpawn(command, Ω::Hotcomb, ℧::AbstractArray{Bool,1})
    tempcommand = command
    m = eachmatch(r"(\{\{.*?\}\})", tempcommand)
    matches = [replace(x[1], r"\{|\}"=>"") for x in collect(m)] |> unique

    if occursin(r"^[0-9]+,[0-9]+$", matches[end])
        countrange = (x -> x[1]:x[2])(integer.(split(matches[end], ",")))
        tempcommand = replace(tempcommand, "{{$(matches[end])}}"=>"") |> strip
        matches = matches[1:(end-1)]
    elseif occursin(r"^[0-9]+$", matches[end])
        countrange = (x -> x[1]:x[1])(integer(matches[end]))
        tempcommand = replace(tempcommand, "{{$(matches[end])}}"=>"") |> strip
        matches = matches[1:(end-1)]
    else
        countrange = missing
    end

    mykeys = keys(Ω)
    domain = 1:length(mykeys)

    iSet  = [m for m in matches if m[end:end] ∈ ["i", "j"]]
    iSet1 = [m[1:1] for m in iSet]
    iSet2 = [m[1:2] for m in iSet if length(m)>=2]
    iSetend = unique([m[end:end] for m in iSet])

    (length(iSet)>2) && throw("Only one !i, >i, >=i, <=i, <i wildcard allowed with i (or j)")
    (length(iSetend)>1) && throw("Only one type i or j allowed")

    wild  = matches[length.(matches) .== 1]
    wild2 = matches[length.(matches) .!= 1]

    collection = []

    println()
    for i in domain, j in domain
      ((length(wild2)  == 0) || wild2[1][1] ∈ 'i':'j') && (j>1)  && continue
      ("!"   ∈ iSet1)       && (i==j) && continue
      (">"   ∈ iSet1)       && (i>=j) && continue
      ("<"   ∈ iSet1)       && (i<=j) && continue
      (">="  ∈ iSet2)       && (i<j)  && continue
      ("<="  ∈ iSet2)       && (i>j)  && continue
      #println("$i and $j")

       txtcmd = tempcommand

       (length(wild2)==1) && !(wild2[1][1] ∈ 'i':'j') &&  (txtcmd = subout(txtcmd, j, wild2[1], mykeys))
       (length(wild2)==1) &&  (wild2[1][1] ∈ 'i':'j') &&  (txtcmd = subout(txtcmd, i, wild2[1], mykeys))
       (length(wild)  == 1) && (txtcmd = subout(txtcmd, i, wild[1], mykeys))

       if occursin("~~OUTOFBOUNDS~~", txtcmd)
           ("i" ∈ iSet) && push!(collection, fill(false, length(℧)))
           continue
       end

       ℧∇ = SuperOperatorEval(txtcmd, Ω, ℧)[2]

       println(">>> $txtcmd")

       push!(collection, ℧∇)
    end

    collector = hcat(collection...)

    if (countrange === missing)
      ℧Δ = ℧ .& [all(collector[i,:]) for i in 1:size(collector)[1]]
    else
      ℧Δ = ℧ .& [sum(collector[i,:]) ∈ countrange for i in 1:size(collector)[1]]
    end
    (Ω, ℧Δ)
end

function subout(txtcmd, i, arg, mykeys)
  lookup(vect, i) = i ∈ 1:length(vect) ? vect[i] : "~~OUTOFBOUNDS~~"

  (arg[end] ∈ ['i', 'j'])  && return replace(txtcmd, "{{$arg}}"=>lookup(mykeys,i))

  mod = integer(match(r"([0-9]+$)", arg).match)

  occursin("+", arg) && return replace(txtcmd, "{{$arg}}"=>lookup(mykeys,i+mod))
  occursin("-", arg) && return replace(txtcmd, "{{$arg}}"=>lookup(mykeys,i-mod))

  txtcmd
end

function SuperSuperOperatorEval(command, Ω::Hotcomb, ℧::AbstractArray{Bool,1})
    #println("OperatorEval($command)")
    (sum(℧) == 0) && return (Ω, ℧)
    (!occursin(r"([><=|!+\\-\^&]{4})", command)) && return SuperOperatorEval(command, Ω, ℧)

    m = match(r"(^.*?)([><=|!+\\-\^&]{4})(.*?$)", replace(command, " "=>""))
    left, supersuperoperator, right = m.captures

    υ = copy(℧); ℧η = copy(℧)

    ℧left  = SuperSuperOperatorEval(left ,Ω,℧)[2]
    ℧right = SuperSuperOperatorEval(right,Ω,℧)[2]

    (supersuperoperator == "====") && (℧η = υ .& (℧left .& ℧right))
    (supersuperoperator == "||||") && (℧η = υ .& (℧left .| ℧right))
    (supersuperoperator == "^^^^") && (℧η = υ .& ((℧left .& .!℧right) .| (.!℧left .& ℧right)))

    (Ω, ℧η)
end

function SuperOperatorEval(command, Ω::Hotcomb, ℧::AbstractArray{Bool,1})
    #println("OperatorEval($command)")
    (sum(℧) == 0) && return (Ω, ℧)
    (!occursin(r"([><=|!+\\-\^&]{3})", command)) && return OperatorEval(command, Ω, ℧)
    occursin(r"\{\{.*\}\}", command) && return OperatorSpawn(command, Ω, ℧)

    m = match(r"(^.*?)([><=|!+\\-\\^&]{3})(.*?$)",replace(command, " "=>""))
    left, superoperator, right = m.captures

    υ = copy(℧); ℧η = copy(℧)

    ℧left  = SuperOperatorEval(left ,Ω,℧)[2]
    ℧right = SuperOperatorEval(right,Ω,℧)[2]

    if superoperator == "&&&"
        ℧η = υ .& (℧left .& ℧right)

    elseif superoperator == "^^^"
        ℧η = υ .& ((℧left .& .!℧right) .| (.!℧left .& ℧right))

    # this can be dangerous, false equal to false such as with previous exclusions will cause inconsistencies
    elseif superoperator == "==="
        ℧η = υ .& (℧left .== ℧right)

    #warning this generally will not work
    elseif superoperator == "---"
        ℧η[υ] = (℧left .- ℧right)[υ]

    #warning this generally will not work
    elseif superoperator == "+++"
        ℧η = ℧left .+ ℧right

    elseif superoperator == "|||"
        ℧η = υ .& (℧left .| ℧right)

    elseif superoperator ∈ ["|=>","==>"]
        ℧η[℧left] .= ℧[℧left]  .& ℧right[℧left]

    elseif superoperator ∈ ["<=|","<=="]
        ℧η[℧right] = ℧[℧right] .& ℧left[℧right]

    elseif superoperator ∈ ["<=>","<=>"] # ???????????????????????????????????
        ℧η[℧right]     .=  ℧[℧right]   .&   ℧left[℧right]
        ℧η[.!℧right]   .=  ℧[.!℧right] .& .!℧left[.!℧right]
    end
    (Ω, ℧η)
end

function OperatorEval(command, Ω::Hotcomb, ℧::AbstractArray{Bool,1})
    #println("OperatorEval($command)")

    (sum(℧) == 0) && return (Ω, ℧)
    occursin(r"\{\{.*\}\}", command) && return OperatorSpawn(command, Ω, ℧)

    n = 1:sum(℧); ℧Δ = copy(℧); ℧η = copy(℧)

    # convert a = b|c to a |= b,c
    if occursin("|", command) & occursin(r"(\b| )[|]*=+[|]*(\b| )", command)
      command = replace(command, "|"=>",")
      command = replace(command, r",*=+,*"=>"|=")
    end

    m = match(r"^(.*?)(([><=|!^&]{1,2})\b)(.*?)(\{([0-9]+),?([0-9]+)*\})?$",replace(command, " "=>""))
    left, right, operator, nmin, nmax  = m.captures[[1,4,3,6,7]]

    (nmin === nothing)  && (nmax === nothing)  && (nrange = 1:999)
    !(nmin === nothing) && (nmax === nothing)  && (nrange = integer(nmin):integer(nmin))
    !(nmin === nothing) && !(nmax === nothing) && (nrange = integer(nmin):integer(nmax))

    leftarg  = strip.(split(left,  r"[,&|!]"))
    rightarg = strip.(split(right, r"[,&|!]"))

    leftvals  = hcat([grab(L, Ω, ℧, command=command) for L in leftarg]...)
    rightvals = hcat([grab(R, Ω, ℧, command=command) for R in rightarg]...)

    if operator == "!="
        lcheck = [any(leftvals[i,j] .== rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [!all(lcheck[i,:]) for i in n]

    elseif operator == "|"
        lcheck = [any(leftvals[i,:]  .== 1) for i in n]
        rcheck = [any(rightvals[i,:] .== 1) for i in n]
        ℧Δ = [all(lcheck[i,:]) | all(rcheck[i,:]) for i in n]

    elseif operator == "&"
        lcheck = [all(leftvals[i,:]  .== 1) for i in n]
        rcheck = [all(rightvals[i,:] .== 1) for i in n]
        ℧Δ = [all(lcheck[i,:]) & all(rcheck[i,:]) for i in n]

    elseif operator == "!"
        lcheck = [all(leftvals[i,:]  .== 0) for i in n]
        rcheck = [all(rightvals[i,:] .== 0) for i in n]
        ℧Δ = [all(lcheck[i,:]) & all(rcheck[i,:]) for i in n]

    elseif operator == "^="
         lcheck = [sum(leftvals[i,j]  .== rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
         rcheck = [sum(rightvals[i,j] .== leftvals[i,:])  for i in n, j in 1:size(rightvals)[2]]
         ℧Δ = [sum(lcheck[i,:]) + sum(rcheck[i,:]) == 2 for i in n]

    elseif operator  ∈ ["==","="]
        lcheck = [all(leftvals[i,j] .== rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator ∈ ["|=", "=|"]
        lcheck = [any(leftvals[i,j] .== rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [sum(lcheck[i,:]) ∈ nrange for i in n]

    elseif operator == "<="
        lcheck = [all(leftvals[i,j] .<= rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == "<"
        lcheck = [all(leftvals[i,j] .< rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == "<<"
        lcheck = [all(leftvals[i,j] .< rightvals[i,:].-1) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == ">="
        lcheck = [all(leftvals[i,j] .>= rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == ">"
        lcheck = [all(leftvals[i,j] .> rightvals[i,:]) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    elseif operator == ">>"
        lcheck = [all(leftvals[i,j] .> rightvals[i,:].+1) for i in n, j in 1:size(leftvals)[2]]
        ℧Δ = [all(lcheck[i,:]) for i in n]

    end

    ℧η[℧η] = ℧Δ

    (Ω, ℧η)
end

Ω,℧ = ABparse(["a, b, c, d, e, f  ∈ 1:6", "{{i}} != {{!i}}"]); Ω[℧]

#Ω,℧ = ABparse(["a, b, c, d, e, f  ∈ 6"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c ∈ 1:6", "{{i}} != {{!i}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "b != a,c", "c =| 1,2"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "{{i}} != {{!i}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a|b = 1"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c ∈ 1:3", "a | b"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c ∈ 1:3", "a|b|c"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c ∈ 1:3", "a | b | c"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c ∈ 1:3", "a,b,c = 1"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c ∈ 1:3", "a & b & c"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c ∈ 0:2", "a,b,c = 0"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c ∈ 0:2", "a ! b ! c"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]"]); Ω[℧]

# At least 2 less/more >> or <<
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a << b"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a >> b"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a,b ^= 1"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a=1 ^^^ b=1"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a=1 ^^^^ b=1"]); Ω[℧]

# All equal commands
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a,b |=   1"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a=1 |||  b = 1"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a=1 |||| b = 1"]); Ω[℧]

# Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a,b ^=   1"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a|b = 1 ||| c == 1"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{i}} == {{!i}}"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{i}} != {{!i}}"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{i}} > {{>i}}"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{i}} > {{i+1}} {{2,3}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "{{i}} == 1 {{2}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{j}} > {{j+1}}"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{i}} > {{i+1}} {{0}}"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{j}} <= {{j+1}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c, d  ∈  [1,2,3,4]", "{{i}} != {{!i}}"]); Ω[℧]
Ω,℧ = ABparse(["a, b, c, d  ∈  [1,2,3,4]", "{{i}} != {{!i}} {{0}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c, d  ∈  [1,2,3,4]", "{{i}} = 2 {{2,3}}"]); Ω[℧] # ??????????????????????
Ω,℧ = ABparse(["a, b, c, d  ∈  [1,2,3,4]", "{{i}} != {{!i}} {{1,5}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3,4]", "{{j}} > {{j+1}}"]); Ω[℧]

Ω,℧ = ABparse(["a, b  ∈  [1,2,3]", "a|b = 1 {1}"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b  ∈  [1,2,3]", "a|b = 1 {0}"]); Ω[℧,:]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "b != a,c", "c =| 1,2 {1}"]); Ω[℧,:]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "b != a,c", "c == 1|2"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a < b,c",  "c |= 1,2"]);Ω[℧,:]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "b < 3 |=> a = b"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a == b <=| b << 3"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a == 1 <=> b == c"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "a <= b,c", "c |= 1,2"]); Ω[℧,:]

Ω,℧ = ABparse(["a, b, c     ∈  [1,2,3]", "b , a == c-1", "c |= 1,2"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c     ∈  [1,2,3]", "a == c+b"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c, d  ∈  [1,2,3,4]", "a , b |= c, d"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c, d  ∈  [1,2,3,4]", "a != b, c, d", "b != c,d", "c != d"]); Ω[℧,:]
Ω,℧ = ABparse(["a, b, c, d  ∈  [1,2,3,4]", "a != b, c, d", "b != c,d", "c != d", "a == c+1", "d == a*2"]); Ω[℧,:]

Ω,℧ = ABparse(["a, b, c  ∈  [1,2,3]", "b != a,c", "c =| 1,2"]); Ω[℧,:]

Ω,℧ = ABparse(["a, b, c, d, e, f, g  ∈  [1,2,3,4]", "a,b,c,d,e,f,g == 1 +++ a,b,c,d,e,f,g == 1 ==== 2"]); Ω[℧,:]

Ω,℧ = ABparse(["a, b, c     ∈  [1,2,3]", "b , a == c-1"]); Ω[℧,:]


text:
using LinearLeastSquares
using Gadfly

line_mat_x = readdlm("tux_sparse_x.txt")
line_mat_y = readdlm("tux_sparse_y.txt")
line_mat_val = readdlm("tux_sparse_val.txt")
line_vals = readdlm("tux_sparse_lines.txt")

# Form the sparse matrix from the data
# Image is 50 x 50
img_size = 50
# The number of pixels in the image
num_pixels = img_size * img_size

line_mat = spzeros(3300, num_pixels)

num_vals = length(line_mat_val)

for i in 1:num_vals
  x = int(line_mat_x[i])
  y = int(line_mat_y[i])
  line_mat[x + 1, y + 1] = line_mat_val[i]
end

x = Variable(num_pixels)
objective = sum_squares(line_mat * x - line_vals)
optval = minimize!(objective)

rows = zeros(img_size*img_size)
cols = zeros(img_size*img_size)
for i = 1:img_size
  for j = 1:img_size
    rows[(i-1)*img_size + j] = i
    cols[(i-1)*img_size + j] = img_size + 1 - j
  end
end

p = plot(
  x=rows, y=cols, color=reshape(evaluate(x), img_size, img_size), Geom.rectbin,
  Scale.ContinuousColorScale(Scale.lab_gradient(color("black"), color("white")))
)
#draw(PNG("tomography.png", 16cm, 14cm), p)


text:
# Default link for recording the information about routing links.
"""
Default implementation of [`RoutingLink`](@ref)

Simple container for channel indices, cost, and capacity.

Constructors
------------
$(METHODLIST)
"""
struct BasicRoutingLink <: RoutingLink
    "Vector of channels curently assigned to the link."
    channels::Vector{ChannelIndex}

    "Base cost of using this link."
    cost::Float64

    """
    Number of channels that can be mapped to this link without it being
    considered congested.
    """
    capacity::Int64
end
function BasicRoutingLink(; cost = 1.0, capacity = 1)
    return BasicRoutingLink(ChannelIndex[], cost, capacity)
end

# API Methods
@inline channels(a::RoutingLink) = a.channels
@inline cost(a::RoutingLink) = a.cost
@inline capacity(a::RoutingLink) = a.capacity
@inline occupancy(a::RoutingLink) = length(a.channels)

addchannel(a::RoutingLink, c) = push!(channels(a), c)

function remchannel(a::RoutingLink, c)
    return deleteat!(a.channels, findfirst(x -> x == c, a.channels))
end

# Derived Methods
iscongested(a::RoutingLink) = occupancy(a) > capacity(a)
iscongested(links::Vector{<:RoutingLink}) = any(iscongested, links)

################################################################################
# Routing Link Documentation
################################################################################

@doc """
    channels(link::RoutingLink) :: Vector{ChannelIndex}

Return list of channels currently occupying `link`.
""" channels

@doc """
    cost(link::RoutingLink) :: Float64

Return the base cost of a channel using `link` as a routing resource.
""" cost

@doc """
    capacity(link::RoutingLink) :: Real

Return the capacity of `link`.
""" capacity

@doc """
    occupancy(link::RoutingLink)

Return the number of channels currently using `link`.
""" occupancy

@doc """
    addchannel(link::RoutingLink, channel::ChannelIndex)

Record that `channel` is using `link`.
""" addchannel

@doc """
    remchannel(link::RoutingLink, channel::ChannelIndex)

Remove `channel` from the list of channels using `link`.
""" remchannel


text:
### Transformtion
abstract DataTransform

# apply the transform
transform!{T<:AbstractFloat, S<:DataTransform}(t::S, x::DenseArray{T,1}) = transform!(x, t, x)
transform!{T<:AbstractFloat, S<:DataTransform}(t::S, x::DenseArray{T,2}) = transform!(x, t, x)

transform{T<:Real, S<:DataTransform}(t::S, x::DenseArray{T,1}) = transform!(Array(T, size(x)), t, x)
transform{T<:Real, S<:DataTransform}(t::S, x::DenseArray{T,2}) = transform!(Array(T, size(x)), t, x)

# reconstruct the original data from transformed values
reconstruct!{T<:AbstractFloat, S<:DataTransform}(t::S, x::DenseArray{T,1}) = reconstruct!(x, t, x)
reconstruct!{T<:AbstractFloat, S<:DataTransform}(t::S, x::DenseArray{T,2}) = reconstruct!(x, t, x)

reconstruct{T<:Real, S<:DataTransform}(t::S, y::DenseArray{T,1}) = reconstruct!(Array(T, size(y)), t, y)
reconstruct{T<:Real, S<:DataTransform}(t::S, y::DenseArray{T,2}) = reconstruct!(Array(T, size(y)), t, y)

# Z-score transformation
immutable ZScoreTransform{T<:Real} <: DataTransform
    dim::Int
    mean::AbstractVector{T}
    scale::AbstractVector{T}

    ZScoreTransform() = new(0, Float64[], Float64[])

    function ZScoreTransform{T}(d::Int, m::AbstractVector{T}, s::AbstractVector{T})
        lenm = length(m)
        lens = length(s)
        lenm == d || lenm == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lens == d || lens == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new(d, m, s)
    end
end

indim(t::ZScoreTransform) = t.dim
outdim(t::ZScoreTransform) = t.dim

# fit a z-score transform
function fit{T<:Real}(::Type{ZScoreTransform}, X::DenseArray{T,2}; center::Bool=true, scale::Bool=true)
    d, n = size(X)
    n >= 2 || error("X must contain at least two columns.")

    m = vec(mean(X, 2))
    s = std(X, 2, mean=m)

    return ZScoreTransform{T}(d, (center ? map(T, vec(m)) : Array(T, 0)),
                              (scale ? map(T, vec(s)) : Array(T, 0)))
end

function transform!{YT<:Real,XT<:Real}(y::DenseArray{YT,1}, t::ZScoreTransform, x::DenseArray{XT,1})
    d = t.dim
    length(x) == length(y) == d || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(x, y)
                copy!(y, x)
            end
        else
            for i = 1:d
                @inbounds y[i] = x[i] / s[i]
            end
        end
    else
        if isempty(s)
            for i = 1:d
                @inbounds y[i] = x[i] - m[i]
            end
        else
            for i = 1:d
                @inbounds y[i] = (x[i] - m[i]) / s[i]
            end
        end
    end
    return y
end

function transform!{YT<:Real,XT<:Real}(y::DenseArray{YT,2}, t::ZScoreTransform, x::DenseArray{XT,2})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(x,2)
    size(y,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(x, y)
                copy!(y, x)
            end
        else
            for j = 1:n
                for i = 1:d
                    @inbounds y[i, j] = x[i, j] / s[i]
                end
            end
        end
    else
        if isempty(s)
            for j = 1:n
                for i = 1:d
                    @inbounds y[i,j] = x[i,j] - m[i]
                end
            end
        else
            for j = 1:n
                for i = 1:d
                    @inbounds y[i, j] = (x[i, j] - m[i]) / s[i]
                end
            end
        end
    end
    return y
end

function reconstruct!{YT<:Real,XT<:Real}(x::DenseArray{XT,1}, t::ZScoreTransform, y::DenseArray{YT,1})
    d = t.dim
    length(x) == length(y) == d || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(y, x)
                copy!(x, y)
            end
        else
            for i = 1:d
                @inbounds x[i] = y[i] * s[i]
            end
        end
    else
        if isempty(s)
            for i = 1:d
                @inbounds x[i] = y[i] + m[i]
            end
        else
            for i = 1:d
                @inbounds x[i] = y[i] * s[i] + m[i]
            end
        end
    end
    return x
end

function reconstruct!{YT<:Real,XT<:Real}(x::DenseArray{XT,2}, t::ZScoreTransform, y::DenseArray{YT,2})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(y,2)
    size(x,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(y, x)
                copy!(x, y)
            end
        else
            for j = 1:n
                for i = 1:d
                    @inbounds x[i, j] = y[i, j] * s[i]
                end
            end
        end
    else
        if isempty(s)
            for j = 1:n
                for i = 1:d
                    @inbounds x[i, j] = y[i, j] + m[i]
                end
            end
        else
            for j = 1:n
                for i = 1:d
                    @inbounds x[i, j] = y[i, j] * s[i] + m[i]
                end
            end
        end
    end
    return x
end

# UnitRangeTransform normalization

immutable UnitRangeTransform  <: DataTransform
    dim::Int
    unit::Bool
    min::Vector{Float64}
    scale::Vector{Float64}

    UnitRangeTransform() = new(0, false, Float64[], Float64[])

    function UnitRangeTransform(d::Int, unit::Bool, min::Vector{Float64}, max::Vector{Float64})
        lenmin = length(min)
        lenmax = length(max)
        lenmin == d || lenmin == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lenmax == d || lenmax == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new(d, unit, min, max)
    end
end

indim(t::UnitRangeTransform) = t.dim
outdim(t::UnitRangeTransform) = t.dim

function fit{T<:Real}(::Type{UnitRangeTransform}, X::DenseArray{T,2}; unit::Bool=true)
    d, n = size(X)

    tmin = Array(Float64, d)
    tmax = Array(Float64, d)
    copy!(tmin, X[:, 1])
    copy!(tmax, X[:, 1])
    for j = 2:n
        @inbounds for i = 1:d
            if X[i, j] < tmin[i]
                tmin[i] = X[i, j]
            elseif X[i, j] > tmax[i]
                tmax[i] = X[i, j]
            end
        end
    end
    for i = 1:d
        @inbounds tmax[i] = 1.0 / (tmax[i] - tmin[i])
    end
    return UnitRangeTransform(d, unit, tmin, tmax)
end

function transform!{YT<:Real,XT<:Real}(y::DenseArray{YT,1}, t::UnitRangeTransform, x::DenseArray{XT,1})
    d = t.dim
    length(x) == length(y) == d || throw(DimensionMismatch("Inconsistent dimensions."))

    tmin = t.min
    tscale = t.scale

    if t.unit
        for i = 1:d
            @inbounds y[i] = (x[i] - tmin[i]) * tscale[i]
        end
    else
        for i = 1:d
            @inbounds y[i] = x[i] * tscale[i]
        end
    end
    return y
end

function transform!{YT<:Real,XT<:Real}(y::DenseArray{YT,2}, t::UnitRangeTransform, x::DenseArray{XT,2})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(x,2)
    size(y,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    tmin = t.min
    tscale = t.scale

    if t.unit
        for j = 1:n
            for i = 1:d
                @inbounds y[i, j] = (x[i,j] - tmin[i]) * tscale[i]
            end
        end
    else
        for j = 1:n
            for i = 1:d
                @inbounds y[i, j] = x[i, j] * tscale[i]
            end
        end
    end
    return y
end

function reconstruct!{YT<:Real,XT<:Real}(x::DenseArray{XT,1}, t::UnitRangeTransform, y::DenseArray{YT,1})
    d = t.dim
    length(x) == length(y) == d || throw(DimensionMismatch("Inconsistent dimensions."))

    tmin = t.min
    tscale = t.scale

    if t.unit
        for i = 1:d
            @inbounds x[i] = y[i] / tscale[i] +  tmin[i]
        end
    else
        for i = 1:d
            @inbounds x[i] = y[i] / tscale[i]
        end
    end
    return x
end

function reconstruct!{YT<:Real,XT<:Real}(x::DenseArray{XT,2}, t::UnitRangeTransform, y::DenseArray{YT,2})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(y,2)
    size(x,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    tmin = t.min
    tscale = t.scale

    if t.unit
        for j = 1:n
            for i = 1:d
                @inbounds x[i,j] = y[i,j] / tscale[i] + tmin[i]
            end
        end
    else
        for j = 1:n
            for i = 1:d
                @inbounds x[i,j] = y[i,j] / tscale[i]
            end
        end
    end
    return x
end