text:
""" AST and pretty printer for Graphviz's DOT language.

References:

- DOT grammar: http://www.graphviz.org/doc/info/lang.html
- DOT language guide: http://www.graphviz.org/pdf/dotguide.pdf
"""
module Graphviz
export Expression, Statement, Attributes, Graph, Digraph, Subgraph,
  Node, NodeID, Edge, pprint, run_graphviz, to_graphviz

using DataStructures: OrderedDict
import LightGraphs, MetaGraphs
using LightGraphs: edges, vertices, src, dst
using MetaGraphs: get_prop, props
using Parameters

# AST
#####

abstract type Expression end
abstract type Statement <: Expression end

""" AST type for Graphviz's "HTML-like" node labels.

For now, the HTML content is just a string.
"""
struct Html
  content::String
end
Base.print(io::IO, html::Html) = print(io, html.content)

const AttributeValue = Union{String,Html}
const Attributes = OrderedDict{Symbol,AttributeValue}

as_attributes(attrs::Attributes) = attrs
as_attributes(d::OrderedDict) = Attributes(Symbol(k) => d[k] for k in keys(d))
as_attributes(d::AbstractDict) =
  Attributes(Symbol(k) => d[k] for k in sort!(collect(keys(d))))

@with_kw_noshow struct Graph <: Expression
  name::String
  directed::Bool
  stmts::Vector{Statement}=Statement[]
  graph_attrs::Attributes=Attributes()
  node_attrs::Attributes=Attributes()
  edge_attrs::Attributes=Attributes()
end

Graph(name::String, stmts::Vector{Statement}; kw...) =
  Graph(; name=name, directed=false, stmts=stmts, kw...)
Graph(name::String, stmts::Vararg{Statement}; kw...) =
  Graph(; name=name, directed=false, stmts=collect(stmts), kw...)
Digraph(name::String, stmts::Vector{Statement}; kw...) =
  Graph(; name=name, directed=true, stmts=stmts, kw...)
Digraph(name::String, stmts::Vararg{Statement}; kw...) =
  Graph(; name=name, directed=true, stmts=collect(stmts), kw...)

@with_kw_noshow struct Subgraph <: Statement
  name::String="" # Subgraphs can be anonymous
  stmts::Vector{Statement}=Statement[]
  graph_attrs::Attributes=Attributes()
  node_attrs::Attributes=Attributes()
  edge_attrs::Attributes=Attributes()
end

Subgraph(stmts::Vector{Statement}; kw...) = Subgraph(; stmts=stmts, kw...)
Subgraph(stmts::Vararg{Statement}; kw...) = Subgraph(; stmts=collect(stmts), kw...)
Subgraph(name::String, stmts::Vector{Statement}; kw...) =
  Subgraph(; name=name, stmts=stmts, kw...)
Subgraph(name::String, stmts::Vararg{Statement}; kw...) =
  Subgraph(; name=name, stmts=collect(stmts), kw...)

struct Node <: Statement
  name::String
  attrs::Attributes
end
Node(name::String, attrs::AbstractDict) = Node(name, as_attributes(attrs))
Node(name::String; attrs...) = Node(name, attrs)

struct NodeID <: Expression
  name::String
  port::String
  anchor::String
  NodeID(name::String, port::String="", anchor::String="") = new(name, port, anchor)
end

struct Edge <: Statement
  path::Vector{NodeID}
  attrs::Attributes
end
Edge(path::Vector{NodeID}, attrs::AbstractDict) = Edge(path, as_attributes(attrs))
Edge(path::Vector{NodeID}; attrs...) = Edge(path, attrs)
Edge(path::Vararg{NodeID}; attrs...) = Edge(collect(path), attrs)
Edge(path::Vector{String}, attrs) = Edge(map(NodeID, path), attrs)
Edge(path::Vector{String}; attrs...) = Edge(map(NodeID, path), attrs)
Edge(path::Vararg{String}; attrs...) = Edge(map(NodeID, collect(path)), attrs)

# Bindings
##########

""" Run a Graphviz program.

Assumes that Graphviz is installed on the local system and invokes Graphviz
through its command-line interface.

For bindings to the Graphviz C API, see the the GraphViz.jl package
(https://github.com/Keno/GraphViz.jl). GraphViz.jl is unmaintained at the time
of this writing.
"""
function run_graphviz(graph::Graph; prog::String="dot", format::String="json0")
  gv = open(`$prog -T$format`, "r+")
  pprint(gv.in, graph)
  close(gv.in)
  result = read(gv.out, String)
  if !success(gv)
    error("Graphviz $prog failed with exit code $(gv.exitcode) and signal $(gv.termsignal)")
  end
  result
end

function Base.show(io::IO, ::MIME"image/svg+xml", graph::Graph)
  println(io, run_graphviz(graph, format="svg"))
end

# MetaGraphs
############

""" Convert an attributed graph (MetaGraph) to a Graphviz graph.

This method is usually more convenient than direct AST manipulation for creating
Graphviz graphs. It supports graphs that are directed or undirected, simple or
multi-edged. For more advanced features, like nested subgraphs, you must use
the Graphviz AST directly.
"""
function to_graphviz(g::MetaGraphs.AbstractMetaGraph; multigraph::Bool=false)::Graph
  gv_name(v::Int) = "n$v"
  gv_path(e::LightGraphs.Edge) = [gv_name(src(e)), gv_name(dst(e))]

  # Add Graphviz node for each vertex.
  stmts = Graphviz.Statement[]
  for v in vertices(g)
    push!(stmts, Node(gv_name(v), props(g, v)))
  end

  # Add Graphviz edge for each (multi)edge.
  if multigraph
    for e in edges(g)
      path = gv_path(e)
      append!(stmts, [ Edge(path, attrs) for attrs in get_prop(g, e, :edges) ])
    end
  else
    for e in edges(g)
      push!(stmts, Edge(gv_path(e), props(g, e)))
    end
  end

  attrs = props(g)
  Graph(
    name = get(attrs, :name, "G"),
    directed = LightGraphs.is_directed(g),
    stmts = stmts,
    graph_attrs = as_attributes(get(attrs, :graph, Dict())),
    node_attrs = as_attributes(get(attrs, :node, Dict())),
    edge_attrs = as_attributes(get(attrs, :edge, Dict())),
  )
end

# Pretty-print
##############

""" Pretty-print the Graphviz expression.
"""
pprint(expr::Expression) = pprint(stdout, expr)
pprint(io::IO, expr::Expression) = pprint(io, expr, 0)

function pprint(io::IO, graph::Graph, n::Int)
  indent(io, n)
  print(io, graph.directed ? "digraph " : "graph ")
  print(io, graph.name)
  println(io, " {")
  pprint_attrs(io, graph.graph_attrs, n+2; pre="graph", post=";\n")
  pprint_attrs(io, graph.node_attrs, n+2; pre="node", post=";\n")
  pprint_attrs(io, graph.edge_attrs, n+2; pre="edge", post=";\n")
  for stmt in graph.stmts
    pprint(io, stmt, n+2, directed=graph.directed)
    println(io)
  end
  indent(io, n)
  println(io, "}")
end

function pprint(io::IO, subgraph::Subgraph, n::Int; directed::Bool=false)
  indent(io, n)
  if isempty(subgraph.name)
    println(io, "{")
  else
    print(io, "subgraph ")
    print(io, subgraph.name)
    println(io, " {")
  end
  pprint_attrs(io, subgraph.graph_attrs, n+2; pre="graph", post=";\n")
  pprint_attrs(io, subgraph.node_attrs, n+2; pre="node", post=";\n")
  pprint_attrs(io, subgraph.edge_attrs, n+2; pre="edge", post=";\n")
  for stmt in subgraph.stmts
    pprint(io, stmt, n+2, directed=directed)
    println(io)
  end
  indent(io, n)
  print(io, "}")
end

function pprint(io::IO, node::Node, n::Int; directed::Bool=false)
  indent(io, n)
  print(io, node.name)
  pprint_attrs(io, node.attrs)
  print(io, ";")
end

function pprint(io::IO, node::NodeID, n::Int)
  print(io, node.name)
  if !isempty(node.port)
    print(io, ":")
    print(io, node.port)
  end
  if !isempty(node.anchor)
    print(io, ":")
    print(io, node.anchor)
  end
end

function pprint(io::IO, edge::Edge, n::Int; directed::Bool=false)
  indent(io, n)
  for (i, node) in enumerate(edge.path)
    if i > 1
      print(io, directed ? " -> " : " -- ")
    end
    pprint(io, node, n)
  end
  pprint_attrs(io, edge.attrs)
  print(io, ";")
end

function pprint_attrs(io::IO, attrs::Attributes, n::Int=0;
                      pre::String="", post::String="")
  if !isempty(attrs)
    indent(io, n)
    print(io, pre)
    print(io, " [")
    for (i, (key, value)) in enumerate(attrs)
      if (i > 1) print(io, ",") end
      print(io, key)
      print(io, "=")
      print(io, isa(value, Html) ? "<" : "\"")
      print(io, value)
      print(io, isa(value, Html) ? ">" : "\"")
    end
    print(io, "]")
    print(io, post)
  end
end

indent(io::IO, n::Int) = print(io, " "^n)

end


text:
using SearchLight, SearchLight.QueryBuilder

function initialize_searchlight()
  try
    SearchLight.Configuration.load() |> SearchLight.Database.connect!
    SearchLight.load_resources()
  catch ex
    @error ex
  end

  nothing
end

@async initialize_searchlight()

ext: jl


text:
function _fill_index_compare!(x, r)
    @simd for i in r
        x[i] = i
    end
end
function _compare(dsl, dsr, ::Val{T}; onleft, onright, cols_left, cols_right, check = true, mapformats = false, on_mapformats = [true, true], stable = false, alg = HeapSort, accelerate = false, method = :sort, threads = true, eq = isequal, obs_id_name = :obs_id, multiple_match = false, multiple_match_name = :multiple, drop_obs_id = true, makeunique = false) where T
    names_left = names(dsl)[cols_left]
    names_right = names(dsr)[cols_right]
    if !(mapformats isa AbstractVector)
        mapformats = repeat([mapformats], 2)
    else
        length(mapformats) !== 2 && throw(ArgumentError("`mapformats` must be a Bool or a vector of Bool with size two"))
    end

    if onleft == nothing
        n_dsl = nrow(dsl)
        n_dsr = nrow(dsr)
        total_length = max(n_dsl, n_dsr)
        obs_id_left = _missings(T, total_length)
        obs_id_right = _missings(T, total_length)
        _fill_index_compare!(obs_id_left, 1:n_dsl)
        _fill_index_compare!(obs_id_right, 1:n_dsr)
        res = Dataset(x1=obs_id_left, x2=obs_id_right, copycols = false)
        rename!(res, :x1=>Symbol(obs_id_name, "_left"), :x2=>Symbol(obs_id_name, "_right"))
    else
        res = outerjoin(dsl[!, onleft], dsr[!, onright], on = 1:length(onleft) .=> 1:length(onright), check = check, mapformats = on_mapformats, stable = stable, alg = alg, accelerate = accelerate, method = method, threads = threads, obs_id = true, obs_id_name = obs_id_name, multiple_match = multiple_match, multiple_match_name = multiple_match_name, makeunique = makeunique)
        total_length = nrow(res)
        obs_cols = index(res)[[Symbol(obs_id_name, "_left"), Symbol(obs_id_name, "_right")]]
        obs_id_left = _columns(res)[obs_cols[1]]
        obs_id_right = _columns(res)[obs_cols[2]]
    end
    _info_cols = ncol(res)
    for j in 1:length(cols_left)
            fl = identity
            if mapformats[1]
                fl = getformat(dsl, cols_left[j])
            end
            fr = identity
            if mapformats[2]
                fr = getformat(dsr, cols_right[j])
            end
            _left_type = Core.Compiler.return_type(fl, Tuple{eltype(_columns(dsl)[cols_left[j]])})
            _right_type = Core.Compiler.return_type(fr, Tuple{eltype(_columns(dsr)[cols_right[j]])})
            _res = allocatecol(Core.Compiler.return_type(eq, Tuple{_left_type, _right_type}), total_length)

            _compare_barrier_function!(_res, _columns(dsl)[cols_left[j]], _columns(dsr)[cols_right[j]], fl, fr, eq, obs_id_left, obs_id_right, threads)

            push!(_columns(res), _res)
            push!(index(res),  Symbol(names_left[j]* "=>" * names_right[j]))
    end
    if drop_obs_id
        select!(res, Not([Symbol(obs_id_name, "_left"), Symbol(obs_id_name, "_right")]))
    end
    res
end

function _compare_barrier_function!(_res, xl, xr, fl, fr, eq_fun, obs_id_left, obs_id_right, threads)
    @_threadsfor threads for i in 1:length(_res)
        if ismissing(obs_id_left[i]) || ismissing(obs_id_right[i])
            _res[i] = missing
        else
            _res[i] = eq_fun(fl(xl[obs_id_left[i]]), fr(xr[obs_id_right[i]]))
        end
    end
    _res
end


text:
export FreeLSProcessTomo,
       LSProcessTomo

function build_liou_process_predictor(prep::Vector, obs::Vector)
    exps = [ vec(p)*vec(o)' for o in obs, p in prep ]
    return reduce(vcat, map(m->vec(m)', vec(exps)) )
end

function build_choi_process_predictor(prep::Vector, obs::Vector)
    exps = [ QuantumInfo.choi_liou_involution(vec(o)*vec(p)') for o in obs, p in prep ]
    return reduce(vcat, map(m->vec(m)', vec(exps)) )
end

function isstate(m::Matrix)
    LinearAlgebra.ishermitian(m) && isapprox(LinearAlgebra.tr(m),1.0)
end

"""
`FreeLSProcessTomo(states::Vector,obs::Vector)`

Free (unconstrained) least-squares process tomography algorithm. It is
constructed from a vector of states that are used to probe a process, and
a vector of observables that are used to measure the probe states.
"""
struct FreeLSProcessTomo
    statecount::Int
    obscount::Int
    inputdim::Int
    outputdim::Int
    pred::Matrix{ComplexF64}
    function FreeLSProcessTomo(states::Vector,obs::Vector)
        @assert all([LinearAlgebra.ishermitian(o) for o in obs]) "Observables must be Hermitian"
        @assert all([isstate(o) for o in states])  "States must be valid density matrices"

        pred = build_liou_process_predictor(states, obs)
        outputdim = size(pred,1)
        inputdim = size(pred,2)
        new(length(states),length(obs),inputdim,outputdim,pred)
    end
end

function predict(method::FreeLSProcessTomo, process::Matrix)
    reshape(real(method.pred*vec(process)), (method.statecount,method.obscount))
end

function fit(method::FreeLSProcessTomo,
             means::Vector)
    if length(means) != method.outputdim
        error("The number of sample means does not match the required number of experiments")
    end
    d = round(Int, size(method.pred, 2) |> sqrt)
    reg = method.pred \ means
    return (reshape(reg,d,d),
            LinearAlgebra.norm(method.pred*reg-means,2)^2,
            :Optimal)
end

function fit(method::FreeLSProcessTomo,
             means::Vector,
             vars::Vector{Float64})
    if length(means) != method.outputdim
        error("The number of sample means does not match the required number of experiments")
    end
    if length(vars) != method.outputdim
        error("The number of sample variances does not match the required number of experiments")
    end
    if any(vars.<0)
        error("Samples variances must be positive for generalized least squares process tomography.")
    end
    d = round(Int, size(method.pred,2) |> sqrt)
    # weighted least squares
    reg = (LinearAlgebra.Diagonal(1 ./ sqrt.(vars)) * method.pred) \ (LinearAlgebra.Diagonal(1 ./ sqrt.(vars)) * means)
    return (reshape(reg,d,d),
            LinearAlgebra.norm((method.pred*reg-means) ./ sqrt.(vars),2)^2,
            :Optimal)
end

function fit(method::FreeLSProcessTomo,
             means::Matrix)
    if size(means) != (method.statecount, method.obscount)
        error("Mean matrix must be $((method.statecount, method.obscount)), but is $(size(means))")
    end
    fit(method, vec(means))
end

function fit(method::FreeLSProcessTomo,
             means::Matrix,
             vars::Matrix{Float64})
    if size(means) != (method.statecount, method.obscount)
        error("Mean matrix must be $((method.statecount, method.obscount)), but is $(size(means))")
    elseif size(vars) != (method.statecount, method.obscount)
        error("Pointwise variance matrix must be $((method.statecount, method.obscount)), but is $(size(means))")
    end
    fit(method, vec(means), vec(vars))
end

"""
`LSProcessTomo(states::Vector,obs::Vector)`

Least-squares process tomography algorithm where the result is
contrained to be completely positive (CP) and trace perserving
(TP). It is constructed from a vector of states that are used to probe
a process, and a vector of observables that are used to measure the
probe states.

"""
struct LSProcessTomo
    statecount::Int
    obscount::Int
    inputdim::Int
    outputdim::Int
    pred::Matrix{ComplexF64}
    function LSProcessTomo(states::Vector,obs::Vector)
        @assert all([LinearAlgebra.ishermitian(o) for o in obs]) "Observables must be Hermitian"
        @assert all([isstate(o) for o in states]) "States must be valid density matrices"

        pred = build_choi_process_predictor(states, obs)
        outputdim = size(pred,1)
        inputdim = size(pred,2)
        new(length(states),length(obs),inputdim,outputdim,pred)
    end
end

function predict(method::LSProcessTomo, process::Matrix)
    reshape(real(method.pred*vec(QuantumInfo.choi_liou_involution(process))), (method.statecount,method.obscount))
end

begin
    global fit
    tomo_dim = 0
    ptrb = SparseArrays.spzeros(0,0)

    function trb_sop(da,db)
        sop = zeros(da^2,(da*db)^2)
        for i=0:da-1
            for j=0:da-1
                for k=0:db-1
                    sop += vec(QuantumInfo.ketbra(i,j,da))*vec(kron(QuantumInfo.ketbra(i,j,da),QuantumInfo.ketbra(k,k,db)))'
                end
            end
        end
        return SparseArrays.sparse(sop)
    end

    function fit(method::LSProcessTomo,
                 means::Vector{Float64},
                 vars = ones(length(means));
                 #solver = MosekSolver(LOG=0))
                 solver = SCS.SCSSolver(verbose=0, max_iters=10_000, eps = 1e-8))

        if length(means) != length(vars) || size(method.pred,1) != length(means)
            error("Size of observations and predictons do not match.")
        end
        d4 = size(method.pred,2)
        d2 = Int(sqrt(d4))
        d  = Int(sqrt(d2))

        # We assume that the predictions are always real-valued
        # and we need to do the complex->real translation manually since
        # Convex.jl does not support complex numbers yet
        rpred = real(method.pred)
        ipred = -imag(method.pred)
        ivars = 1 ./ sqrt.(vars)

        if d != tomo_dim
            ptrb = trb_sop(d,d)
        end

        ρr = Convex.Variable(d2,d2)
        ρi = Convex.Variable(d2,d2)

        # least squares problem
#       NOTE: Notation "^2" is not working here.
#        problem = Convex.minimize( Convex.norm( (means - (rpred*vec(ρr) + ipred*vec(ρi))) .* ivars, 2 )^2 )
        problem = Convex.minimize( Convex.square(Convex.norm( (means - (rpred*vec(ρr) + ipred*vec(ρi))) .* ivars, 2 ) ))

        # CP constraint
        problem.constraints += Convex.isposdef([ρr ρi; -ρi ρr])
        # TP constraint
        problem.constraints += reshape(ptrb*vec(ρr),d,d) == QuantumInfo.eye(d)
        problem.constraints += reshape(ptrb*vec(ρi),d,d) == zeros(d,d)

        Convex.solve!(problem, solver)

        return QuantumInfo.choi_liou_involution((ρr.value + 1im*ρi.value)), problem.optval, problem.status
    end
end

function fit(method::LSProcessTomo,
             means::Matrix,
             vars::Matrix{Float64}=ones(size(means)))
    if size(means) != (method.statecount, method.obscount)
        error("Matrix of sample means must be $((method.statecount, method.obscount)), but is $(size(means))")
    end
    if size(vars) != (method.statecount, method.obscount)
        error("Matrix of sample variances must be $((method.statecount, method.obscount)), but is $(size(means))")
    end
    fit(method, vec(means), vec(vars))
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -5.992552298143351 <= q <= 4.750344613921084)

                     add_NL_constraint(m, :(sigmoid(0.47930083070023954 + -0.2811056889538768*sigmoid(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])) + -0.8918376604733185*sigmoid(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])) + -0.8083246377220759*sigmoid(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))) + sigmoid(-0.740550943643274 + -0.007926005977438866*sigmoid(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])) + 0.328386397895855*sigmoid(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])) + 0.6266373732109729*sigmoid(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))) + sigmoid(-0.3322291763851517 + -0.3879865017411772*sigmoid(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])) + 0.13951112377587105*sigmoid(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])) + -0.16669374723112762*sigmoid(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
    setup_problem(problem, parameters::Parameters)

Set up a fixed-dimensional optimization problem.

There are several `setup_problem()` method that accept different type of
`problem` argument:
    * `OptimizationProblem`
    * function (`:NumDimensions` has to be specified in `parameters`)
    * `FunctionBasedProblemFamily` (`:NumDimensions` has to be specified in `parameters`)
"""
function setup_problem(problem::OptimizationProblem, parameters::Parameters)
    return problem, parameters
end

function setup_problem(family::FunctionBasedProblemFamily, parameters::Parameters)
    # If an anydim problem was given the dimension param must have been specified.
    if parameters[:NumDimensions] == :NotSpecified
        throw(ArgumentError("You MUST specify NumDimensions= when a problem family is given"))
    end
    problem = instantiate(family, parameters[:NumDimensions])

    return problem, parameters
end

function setup_problem(func, parameters::Parameters)
    ss = check_and_create_search_space(parameters)

    # Now create an optimization problem with the given information. We currently reuse the type
    # from our pre-defined problems so some of the data for the constructor is dummy.
    problem = FunctionBasedProblem(func, "<unknown>", parameters[:FitnessScheme], ss,
                                   parameters[:TargetFitness])

    # validate fitness: create a random solution from the search space and ensure that fitness(problem) returns fitness_type(problem).
    ind = rand_individual(search_space(problem))
    res = fitness(ind, problem)
    fitnessT = fitness_type(problem)
    if !isa(res, fitnessT)
      throw(ArgumentError("The supplied fitness function does NOT return the expected fitness type $(fitnessT)"*
                          "when called with a potential solution "*
                          "(when called with $(ind) it returned $(res) of type $(typeof(res)) so we cannot optimize it!"))
    end

    return problem, parameters
end

"""
    bboptimize(problem[; parameters::Associative, kwargs...])

Solve given optimization `problem`.

See `setup_problem()` for the types of `problem` supported.
In addition, the `problem` could be `OptController` containing the
results of the previous optimization runs.

The optimization method parameters could be specified via `kwargs` or `parameters` arguments.

Returns `OptimizationResults` instance.

See also `bbsetup()` and [`BlackBoxOptim.OptRunController`](@ref) for a full list of supported parameters.
"""
function bboptimize(optctrl::OptController; kwargs...)
    if length(kwargs) > 0
        update_parameters!(optctrl, kwargs2dict(kwargs))
    end
    run!(optctrl)
end

function bboptimize(functionOrProblem, parameters::Parameters = EMPTY_PARAMS; kwargs...)
    optctrl = bbsetup(functionOrProblem, parameters; kwargs...)
    run!(optctrl)
end

"""
    bbsetup(problem[; parameters::Associative, kwargs...])

Set up optimization method for a given problem.

See `setup_problem()` for the types of `problem` supported.
The optimization method parameters could be specified via `kwargs` or `parameters` arguments.

Returns the initialized `OptController` instance. To actually run the method
call `bboptimize()` or `run!()`.

See also [BlackBoxOptim.OptRunController](@ref) for a full list of supported parameters.
"""
function bbsetup(functionOrProblem, parameters::Parameters = EMPTY_PARAMS; kwargs...)
    parameters = chain(convert(ParamsDict, parameters), kwargs2dict(kwargs))
    problem, params = setup_problem(functionOrProblem, chain(DefaultParameters, parameters))
    check_valid!(params)

    optimizer_func = chain(SingleObjectiveMethods, MultiObjectiveMethods)[params[:Method]]
    optimizer = optimizer_func(problem, params)

    # Now set up a controller for this problem. This will handle
    # application of optimizer, checking for termination conditions
    # as well as collecting the statistics about the number of function evals,
    # keep an archive and top list of candidates.
    ctrl = OptController(optimizer, problem, params)

    return ctrl
end


text:
########################### Parallel FFT tests ################################
# using Distributed
# TODO: make a functional test here
# TODO: define a simple way of establishing Scattering transform on n workers
# addprocs(2)
# @everywhere using Revise
# @everywhere using ScatteringTransform
# m1=2; n1=50; p1=55
# X1 = randn(n1,p1)
# layer1 = stParallel(m1, size(X1), typeBecomes=eltype(X1))
# n, q, dataSizes, outputSizes, resultingSize =
#     ScatteringTransform.calculateSizes(layer1, (-1,-1), size(X1))

# plans = createFFTPlans(layer1, dataSizes, verbose = true, T=eltype(X1))
# fetch(plans[1,1])
# layers = layer1
# padBy = getPadBy(layers.shears[1])
# thing = createRemoteFFTPlan(dataSizes[1], padBy, Float64,false); size(thing)
# ptrThing = remotecall(createRemoteFFTPlan, 2, dataSizes[1][(end-2):(end-1)], padBy, Float64, false)
# tmp = fetch(ptrThing); typeof(tmp)
# FFTs = Array{Future,2}(undef, nworkers(), layers.m+1)

# TODO: integrate these test in some way
# plot(heatmap(abs.((cwt(reshape(f,(1,size(f)...)), layers.shears[1],nScales=50))[1,:,:])'),heatmap(abs.(cwt(f))[5:end,:]))
# cwt(f)
# c = CFW(WT.Morlet())
# ω = [0:ceil(Int, 200/2); -floor(Int,200/2)+1:-1]*2π
# daughters = zeros(200,54)
# for a1 in 0:53
#   daughters[:,a1+1] = WT.Daughter(c, 2.0^(a1/c.scalingFactor), ω)
# end
# heatmap((daughters./[norm(daughters[:,i]) for i=1:54]')')
# norm(daughters[:,53])
# localDaughters = computeWavelets(f, layers.shears[1],nScales=50)
# layers.shears[1]
# daughters[:,5:end]'
# abs.(localDaughters[:,2:end])'
# heatmap(daughters[:,5:end]')
# heatmap(abs.(localDaughters[:,2:end])')
# (daughters[:,5:end]-localDaughters)
# heatmap(abs.(localDaughters'))
# heatmap(abs.(computeWavelets(f, layers.shears[1],nScales=49))')

#TODO fix the 1D tests
# function testConstruction1D(lay::stParallel, m::Int, shearType, averagingLength, averagingType, nShears::Vector{Int}, subsampling::Array{Float64,1}, n::Int)
#   @test lay.m == m
#   @test lay.subsampling==subsampling
#   sizes = ScatteringTransform.sizes(bsplineType(), bspline,lay.subsampling,n)
#   for i=1:length(lay.shears)
#     @test typeof(lay.shears[i]) == shearType[i]
#     @test lay.shears[i].averagingLength == averagingLength[i]
#     @test lay.shears[i].averagingType == averagingType[i]
#     @test numScales(lay.shears[i], sizes[i]) == nShears[i]
#   end
# end
# n = 10214
# f = randn(n)
# lay = stParallel(1, size(f)[end])
# m=3
# lay = stParallel(m, f, subsampling=[8,4,2,1], nScales=[16,8,8,8], CWTType=WT.dog2, averagingLength=[16,4,4,2],averagingType=[:Mother for i=1:(m+1)],boundary=[WT.DEFAULT_BOUNDARY for i=1:(m+1)])
# @testset "construction tests" begin
#   testConstruction1D(lay, 3, [ScatteringTransform.CFWA{Wavelets.WT.PerBoundary} for i=1:4], [16, 4, 4, 2], [:Mother for i=1:4], [165,62,46,40], [8,4,2,1]*1.0, n)
# end
# function testTransform(f::Array{T}, lay::stParallel) where T<:Number
#   n=size(f)[end]
#   @time output = st(f,lay)
#   sizes = ScatteringTransform.sizes(bsplineType(), bspline,lay.subsampling,n)
#   actualSizes = fill(0, size(sizes))
#   actualSizes[1] = size(output.data[1])[end-1]
#   for i=1:length(actualSizes)-2
#     @test size(output.data[i+1])[end-1] == size(output.output[i])[end-1]
#     actualSizes[i+1] = size(output.data[i+1])[end-1]
#   end
#   actualSizes[end] = size(output.output[end])[end-1]
#   @test actualSizes== sizes
#   # check the number of paths grows correctly
#   effectiveSize = length(size(f))==1 ? (1,size(f)[1]) : size(f)
#   numTransformed = ScatteringTransform.getQ(lay, effectiveSize)
#   @test numTransformed[2:end]==([size(out)[end] for out in output.output])[2:end]
#   @test eltype(output.data) == Array{Complex{eltype(f)},1+length(effectiveSize)}
#   @test eltype(output.output) == Array{eltype(f),1+length(effectiveSize)}
# end
# @testset "transform tests" begin
#   f=1.0 .*[1:50; 50:-1:1]
#   lay = stParallel(3,f)
#   # check that various constructions are well-defined
#   testTransform(f,lay)
#   f = randn(1020)
#   m=3
#   lay = stParallel(m, f, subsampling=[8,4,2,1], nScales=[16,8,8,8], CWTType=WT.dog2, averagingLength=[16,4,4,2],averagingType=[:Mother for i=1:(m+1)],boundary=[WT.DEFAULT_BOUNDARY for i=1:(m+1)])
#   testTransform(f, lay)
#   stParallel(3, length(f), [2, 2, 2, 2], [2, 2, 2, 2], WT.dog2)
#   stParallel(3, length(f), 8, [2, 2, 2, 2], WT.dog2)
#   stParallel(3, length(f), [2, 2, 2, 2], 2, WT.dog2)
#   stParallel(3, length(f), 8, 2, WT.dog2)

#   # test on input with multiple samples
#   testTransform(randn(2,3,1020),lay)
#   testTransform(randn(2,1020),lay)
# end
# Scattered(lay, randn(3,1020), "full")

# # Scattered construction tests
# Scattered{Float64,2}(2,1,[im*randn(10,10) for i=1:3],[randn(10,10) for i=1:3])
# layers = stParallel(3, length(f), 8, 4)

# results = Scattered(layers, f)

# f = randn(10,30,500) #there are 10 examples of length 500
# stParallel(3,f[1,1,:])
# #TODO: make a version that can handle input of arbitrary number of initial dimensions
# tmp = stParallel(1,f[1,1,:])
# ScatteringTransform.cwt(f, tmp.shears[1])
# averageNoiseMat = cwt(f,tmp.shears[1],J1=1)
# using FFTW
# A = randn(10,10,30)
# eachindex
# CI = CartesianIndex(size(A))
# for i in Tuple(CI)
#   print(i)
# end
# axes(A,3)
# for x in axes(A)[1:end-1]
#   print(x)
# end
# for x in eachindex.(axes(A)[1:end-1])
#   print("$(size(A[x]))")
# end
# A = randn(10,10,30)
# using LinearAlgebra
# for x in eachindex(view(A,axes(A)[1:end-1]..., 1))
#   print(x)
#   print("$(size(A[x,:]))")
#   println("     ")
# end
# CartesianIndex.(axes(A, 1), axes(A, 2))
# view(A,100,100, :)
# for i in eachindex(view(A,:,1))
#   @show i
# end
# randn(100,100,100).*
# using LinearAlgebra
#
# dw = zeros(10000,10000)
# v = rand(10000)
# h = rand(10000)
# @time for i = 1:10
#   dw += h*v'
# end
# @time for i = 1:10
#   BLAS.gemm!('N', 'T', 1.0, h, v, 1.0, dw)
# end

# TODO write tests for thinST
#      no outputSubsample
#      fixed size outputSubsample
#      fraction outputSubsample

# TODO: finish test folder transformations
# testing folder transformations
# using Distributed
# addprocs(8)
# @everywhere using Revise
# @everywhere using HDF5, JLD, ScatteringTransform
# using ScatteringTransform
# mkpath("tmpData")
# data = randn(20, 3, 5, 100)
# for i=1:20
#   mkpath("tmpData/tmp$(i)")
#   for j=1:3
#     h5write("tmpData/tmp$(i)/tmp$(j).h5", "data",data[i,j,:,:])
#   end
# end
# # direct transformation
# layers = stParallel(2,data)
# @time thinOutput = thinSt(data,layers,outputSubsample=(-1,3))
# # auxillary functions
# targetDir = "tmpOutput"
# function defaultLoadFunction(filename)
#   return (h5read(filename,"data"), true)
# end
# transformFolder("tmpData", targetDir, layers; separate=true, loadThis = defaultLoadFunction, postSubsample=(-1,3))
# load("tmpOutput/settings.jld","layers")
# load("/VastExpanse/data/SSAM2_motionCompensatedData/newData/lipschitz2Layer3Samples/settings.jld  ")
# stParallel
# @testset "folder transformations"
# end
# for i=1:20
#   mkpath("tmpData/tmp$(i)")
#   for j=1:3
#     h5write("tmpData/tmp$(i)/tmp$(j).h5", "data",data[i,j,:,:])
#   end
# end
# "α_"[1]
# layers
# using JLD
# module ScatteringTransform
# include("src/modifiedTransforms.jl")
# export stParallel
# struct stParallel{T}
#   m::Int # the number of layers, not counting the zeroth layer
#   n::Int # the length of a single entry
#   shears::Array{T} # the array of the transforms; the final of these is used only for averaging, so it has length m+1
#   subsampling::Array{Float64,1} # for each layer, the rate of subsampling. There is one of these for layer zero as well, since the output is subsampled, so it should have length m+1
#   stParallel{T}(m::Int, n::Int, shears::Array{T}, subsampling::Array{Float64,1}) where {T} = new(m,n, shears, subsampling)
# end
# end
# @everywhere
# ¥
# "¥_"[2]


text:
using RDatasets: dataset
using SVM
using JLD

iris = dataset("datasets", "iris")
X = convert(Array, iris[:, 1:4])'
p, n = size(X)
y = [species == "setosa" ? 1.0 : -1.0 for species in iris[:Species]]
model = svm(X, y)

save("model.jld", "model", model)


text:
@testset "nlsolve" begin
    f(x) = x^2 + 3x + 1
    j(x) = 2x + 3

    r, x, s = nlsolve(f, j, -2.5)
    @test s == :success

    r, x, s = nlsolve(f, j, -1.0)
    @test s == :success
end


text:
#==============================================================================================#
#                                      TypeTree.jl module                                      #
#==============================================================================================#

module TypeTree

# Includes:
include("interface.jl")
include("functions.jl")

end # module


text:
"""
    readwdir(fun::Function, gl::GitLink; tout = 60.0))

Allow to access the working directory of the GitLink.
The function `fun(working_dir)` will be executed and it should access
the files into `working_dir`.
Any changes to those files won't be commited (please do not mess up the `.git` folder).
It is recommended that `fun` not to be an expensive function.
This method will sleep till (or timeout `tout`) the GitLink lock is free (which must by must of the time, but...).
Returns `true` if the action was successful.
"""
function readwdir(fun::Function, gl::GitLink;
        tout = 60.0
    )

    ok_flag = false
    wdir = repo_dir(gl)
    mkpath(wdir)

    # TODO: control lock vtime
    lock(gl; tout) do
        fun(wdir)
        ok_flag = true
    end
    return ok_flag

end


text:
#
# Dimension: translates sequences of Symbols, Strings, or Ints to
# ordinal values, so the Dimension keys can be used for lookup in
# Arrays in Mimi. A Dimension can be declared one of several ways:
# dim = Dimension(:foo, :bar, :baz)      # varargs
# dim = Dimension([:foo, :bar, :baz])    # Vector
# dim = Dimension(2010:2100)             # Range
# dim = Dimension(4)                     # Same as 1:4
#
# Similarly, values can be referenced several ways:
# value  = Dimension[:foo]               # single lookup
# values = Dimension[(:foo, :bar)]       # Tuple of values
# values = Dimension[[2010, 2020, 2030]] # Vector of values
# values = Dimension[2010, 2020, 2030]   # Varargs
# values = Dimension[2010:10:2030]       # AbstractRange
#
using DataStructures

key_type(dim::Dimension{T}) where {T <: DimensionKeyTypes} = T

#
# Iteration and basic dictionary methods are delegated to the internal dict
#
Base.length(dim::Dimension) = length(dim.dict)
Base.iterate(dim::Dimension, state...) = iterate(dim.dict, state...)

Base.keys(dim::Dimension)   = keys(dim.dict)
Base.values(dim::Dimension) = values(dim.dict)

Base.get(dim::Dimension, key::Union{Number, Symbol, String}, default::Any) = get(dim.dict, key, default)

Base.getindex(dim::Dimension, key::Colon) = collect(values(dim.dict))

Base.getindex(dim::AbstractDimension, key::Union{Number, Symbol, String}) = getindex(dim.dict, key)

# Support dim[[:foo, :bar, :baz]], dim[(:foo, :bar, :baz)], and dim[2010:2020]
Base.getindex(dim::AbstractDimension, keys::Union{Vector{T} where T, Tuple, AbstractRange}) = [getindex(dim.dict, key) for key in keys]

Base.getindex(dim::AbstractDimension, keys...) = getindex(dim, keys)

Base.length(dim::RangeDimension) = length(dim.range)
Base.iterate(dim::RangeDimension, state...)  = iterate(dim.range, state...)

Base.keys(dim::RangeDimension)   = collect(dim.range)
Base.values(dim::RangeDimension) = collect(1:length(dim.range))

# Get first/last value from OrderedDict of keys
Base.firstindex(dim::AbstractDimension) = dim.dict.keys[1]
Base.lastindex(dim::AbstractDimension)  = dim.dict.keys[length(dim)]

#
# Compute the index of a "key" (e.g., a year) in the range.
#
function Base.get(dim::RangeDimension, key::Int, default::Any=0)
    r = dim.range
    i = key - r.start
    return i == 0 ? 1 : (i % r.step != 0 ? default : 1 + div(i, r.step))
end

# Support dim[[2010, 2020, 2030]], dim[(:foo, :bar, :baz)], and dim[2010:2050]
Base.getindex(dim::RangeDimension, keys::Union{Vector{Int}, Tuple, AbstractRange}) = [get(dim, key, 0) for key in keys]

# Symbols are added to the dim_dict in @defcomp (with value of nothing), but are set later using set_dimension!
has_dim(obj::AbstractCompositeComponentDef, name::Symbol) = (haskey(obj.dim_dict, name) && obj.dim_dict[name] !== nothing)

isuniform(obj::AbstractCompositeComponentDef) = obj.is_uniform

set_uniform!(obj::AbstractCompositeComponentDef, value::Bool) = (obj.is_uniform = value)

dimension(obj::AbstractCompositeComponentDef, name::Symbol) = obj.dim_dict[name]

dim_names(obj::AbstractCompositeComponentDef, dims::Vector{Symbol}) = [dimension(obj, dim) for dim in dims]

dim_count_dict(obj::AbstractCompositeComponentDef) = Dict([name => length(value) for (name, value) in dim_dict(obj)])

dim_counts(obj::AbstractCompositeComponentDef, dims::Vector{Symbol}) = [length(dim) for dim in dim_names(obj, dims)]
dim_count(obj::AbstractCompositeComponentDef, name::Symbol) = length(dimension(obj, name))

dim_keys(obj::AbstractCompositeComponentDef, name::Symbol)   = collect(keys(dimension(obj, name)))
dim_values(obj::AbstractCompositeComponentDef, name::Symbol) = collect(values(dimension(obj, name)))

"""
    set_dimension!(ccd::CompositeComponentDef, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange})

Set the values of `ccd` dimension `name` to integers 1 through `count`, if `keys` is
an integer; or to the values in the vector or range if `keys` is either of those types.
"""
function set_dimension!(ccd::AbstractCompositeComponentDef, name::Symbol, keys::Union{Int, Vector, Tuple, AbstractRange})
    redefined = has_dim(ccd, name)
    # if redefined
    #     @warn "Redefining dimension :$name"
    # end

    dim = Dimension(keys)

    if name == :time
        _set_run_period!(ccd, keys[1], keys[end])
        propagate_time!(ccd, dim)
        set_uniform!(ccd, isuniform(keys))
    end

    return set_dimension!(ccd, name, dim)
end

function set_dimension!(obj::AbstractComponentDef, name::Symbol, dim::Dimension)
    dirty!(obj)
    obj.dim_dict[name] = dim

    if name == :time
        for subcomp in compdefs(obj)
            set_dimension!(subcomp, :time, dim)
        end
    end
    return dim
end

function add_dimension!(comp::AbstractComponentDef, name)
    # generally, we add dimension name with nothing instead of a Dimension instance,
    # but in the case of an Int name, we create the "anonymous" dimension on the fly.
    dim = (name isa Int) ? Dimension(name) : nothing
    comp.dim_dict[Symbol(name)] = dim                         # TBD: test this
end

# Note that this operates on the registered comp, not one added to a composite
add_dimension!(comp_id::ComponentId, name) = add_dimension!(compdef(comp_id), name)

function dim_names(ccd::AbstractCompositeComponentDef)
    dims = OrderedSet{Symbol}()             # use a set to eliminate duplicates
    for cd in compdefs(ccd)
        union!(dims, keys(dim_dict(cd)))    # TBD: test this
    end

    return collect(dims)
end

dim_names(comp_def::AbstractComponentDef, datum_name::Symbol) = dim_names(datumdef(comp_def, datum_name))

dim_count(def::AbstractDatumDef) = length(dim_names(def))


text:
import Base

"""
Struct for WilsonFermion
"""
struct WilsonFermion_4D_wing{NC,NDW} <: WilsonFermion_4D{NC} #AbstractFermionfields_4D{NC}

    f::Array{ComplexF64,6}
    NC::Int64
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    NG::Int64
    NDW::Int64
    Dirac_operator::String
    #BoundaryCondition::Vector{Int8}

    function WilsonFermion_4D_wing(NC::T,NX::T,NY::T,NZ::T,NT::T) where T<: Integer
        NG = 4
        NDW = 1
        #@assert NDW == 1 "only NDW = 1 is supported. Now NDW = $NDW"
        f = zeros(ComplexF64,NC,NX+2NDW,NY+2NDW,NZ+2NDW,NT+2NDW,NG)
        Dirac_operator = "Wilson"
        return new{NC,NDW}(f,NC,NX,NY,NZ,NT,NG,NDW,Dirac_operator)
    end

    function WilsonFermion_4D_wing{NC}(NX::T,NY::T,NZ::T,NT::T) where {T<: Integer,NC}
        NG = 4
        NDW = 1
        #@assert NDW == 1 "only NDW = 1 is supported. Now NDW = $NDW"
        f = zeros(ComplexF64,NC,NX+2NDW,NY+2NDW,NZ+2NDW,NT+2NDW,NG)
        Dirac_operator = "Wilson"
        return new{NC,NDW}(f,NC,NX,NY,NZ,NT,NG,NDW,Dirac_operator)
    end

end

function Base.length(x::WilsonFermion_4D_wing{NC,NDW}) where {NC,NDW}
    return NC*x.NX*x.NY*x.NZ*x.NT*x.NG
end

function Base.setindex!(x::WilsonFermion_4D_wing{NC,NDW},v,i1,i2,i3,i4,i5,i6)  where {NC,NDW}
    @inbounds x.f[i1,i2 + NDW,i3 + NDW,i4 + NDW,i5 + NDW,i6] = v
end

function Base.getindex(x::WilsonFermion_4D_wing{NC,NDW},i1,i2,i3,i4,i5,i6) where {NC,NDW}
    @inbounds return x.f[i1,i2 .+ NDW,i3 .+ NDW,i4 .+ NDW,i5 .+ NDW,i6]
end

function Base.getindex(x::WilsonFermion_4D_wing{NC,NDW},i1::N,i2::N,i3::N,i4::N,i5::N,i6::N) where {NC,NDW,N<: Integer}
    @inbounds return  x.f[i1,i2 + NDW,i3 + NDW,i4 + NDW,i5 + NDW,i6]
end

function setindex_global!()
end

#=
function Base.getindex(F::Shifted_fermionfields_4D{NC,WilsonFermion_4D_wing{NC,NDW}},i1::N,i2::N,i3::N,i4::N,i5::N,i6::N)  where {NC,NDW,N <: Integer}

    @inbounds begin
    si2 = i2 + NDW  + F.shift[1]
    si3 = i3 + NDW  + F.shift[2]
    si4 = i4 + NDW  + F.shift[3]
    si5 = i5 + NDW  + F.shift[4]
    end
    #v = F.parent.f
    #return  v[i1,si2,si3,si4,si5,i6]

   return  @inbounds F.parent.f[i1,si2,si3,si4,si5,i6]
end
=#

#=

function Base.getindex(x::WilsonFermion_4D_wing{NC,NDW},i1,i2,i3,i4,i5,i6) where {NC,NDW}
    #=
    i2new = i2 .+ x.NDW
    i3new = i3 .+ x.NDW
    i4new = i4 .+ x.NDW
    i5new = i5 .+ x.NDW
    @inbounds return x.f[i1,i2new,i3new,i4new,i5new,i6]
    =#
    #return @inbounds Base.getindex(x.f,i1,i2 .+ NDW,i3 .+ NDW,i4 .+ NDW,i5 .+ NDW,i6)
    @inbounds return x.f[i1,i2 .+ NDW,i3 .+ NDW,i4 .+ NDW,i5 .+ NDW,i6]
end

=#

function Base.similar(x::T) where T <: WilsonFermion_4D_wing
    return WilsonFermion_4D_wing(x.NC,x.NX,x.NY,x.NZ,x.NT)
end

#=
function Base.getindex(x::T,i1,i2,i3,i4,i5,i6) where T <: WilsonFermion_4D_wing{NC} where NC
    i2new = i2 + x.NDW
    i3new = i3 + x.NDW
    i4new = i4 + x.NDW
    i5new = i5 + x.NDW
    f1 = x.f
    #println((i1,i2new,i3new,i4new,i5new,i6))
    @inbounds v = f1[i1,i2new,i3new,i4new,i5new,i6]
    #x.f[i1,i2new,i3new,i4new,i5new,i6]
    return v
    @inbounds return x.f[i1,i2new,i3new,i4new,i5new,i6]
    #@inbounds return x.f[i1,i2 .+ x.NDW,i3 .+ x.NDW,i4 .+ x.NDW,i5 .+ x.NDW,i6]
end
=#

function set_wing_fermion!(a::WilsonFermion_4D_wing{NC,NDW},boundarycondition) where {NC,NDW}
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX

    #!  X-direction
    for ialpha=1:4
        for it=1:NT
            for iz = 1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,0,iy,iz,it,ialpha] = boundarycondition[1]*a[k,NX,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,NX+1,iy,iz,it,ialpha] = boundarycondition[1]*a[k,1,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    #Y-direction
    for ialpha = 1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,0,iz,it,ialpha] = boundarycondition[2]*a[k,ix,NY,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,NY+1,iz,it,ialpha] = boundarycondition[2]*a[k,ix,1,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        # Z-direction
        for it=1:NT
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,0,it,ialpha] = boundarycondition[3]*a[k,ix,iy,NZ,it,ialpha]
                        a[k,ix,iy,NZ+1,it,ialpha] = boundarycondition[3]*a[k,ix,iy,1,it,ialpha]

                    end
                end
            end
        end

        #T-direction
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,iz,0,ialpha] = boundarycondition[4]*a[k,ix,iy,iz,NT,ialpha]
                        a[k,ix,iy,iz,NT+1,ialpha] =boundarycondition[4]*a[k,ix,iy,iz,1,ialpha]
                    end
                end
            end
        end

    end

end

function set_wing_fermion!(a::WilsonFermion_4D_wing{NC,NDW},boundarycondition,iseven::Bool) where {NC,NDW}
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX

    #!  X-direction
    @inbounds for ialpha=1:4
        for it=1:NT
            for iz = 1:NZ
                for iy=1:NY
                    ix = NX
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        @simd for k=1:NC
                            a[k,0,iy,iz,it,ialpha] = boundarycondition[1]*a[k,NX,iy,iz,it,ialpha]
                        end
                    end
                end
            end
        end
    end

    @inbounds for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    ix = 1
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        @simd for k=1:NC
                            a[k,NX+1,iy,iz,it,ialpha] = boundarycondition[1]*a[k,1,iy,iz,it,ialpha]
                        end
                    end
                end
            end
        end
    end

    #Y-direction
    @inbounds for ialpha = 1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    iy = NY
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven

                        @simd for k=1:NC
                            a[k,ix,0,iz,it,ialpha] = boundarycondition[2]*a[k,ix,NY,iz,it,ialpha]
                        end
                    end
                end

            end
        end
    end

    @inbounds for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    iy = 1
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven

                        @simd for k=1:NC
                            a[k,ix,NY+1,iz,it,ialpha] = boundarycondition[2]*a[k,ix,1,iz,it,ialpha]
                        end
                    end
                end
            end
        end
    end

    @inbounds for ialpha=1:4
        # Z-direction
        for it=1:NT
            for iy=1:NY
                for ix=1:NX
                    iz = NX
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        @simd for k=1:NC
                            a[k,ix,iy,0,it,ialpha] = boundarycondition[3]*a[k,ix,iy,NZ,it,ialpha]
                        end
                    end

                    iz = 1
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        @simd for k=1:NC
                            a[k,ix,iy,NZ+1,it,ialpha] = boundarycondition[3]*a[k,ix,iy,1,it,ialpha]
                        end
                    end
                end
            end
        end

        #T-direction
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    it = NT
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        @simd for k=1:NC
                            a[k,ix,iy,iz,0,ialpha] = boundarycondition[4]*a[k,ix,iy,iz,NT,ialpha]
                        end
                    end

                    it = 1
                    evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                    if evenodd == iseven
                        @simd for k=1:NC
                            a[k,ix,iy,iz,NT+1,ialpha] =boundarycondition[4]*a[k,ix,iy,iz,1,ialpha]
                        end
                    end

                end
            end
        end

    end

end

function add_fermion!(c::WilsonFermion_4D_wing{NC,NDW},α::Number,a::T1,β::Number,b::T2,iseven) where {NC,NDW,T1 <: Abstractfermion,T2 <: Abstractfermion}#c += alpha*a + beta*b
    n1,n2,n3,n4,n5,n6 = size(c.f)

    @inbounds for i6=1:n6
        for i5=1:n5
            it = i5 -NDW
            for i4=1:n4
                iz = i4 -NDW
                for i3=1:n3
                    iy = i3 - NDW
                    for i2=1:n2
                        ix = i2 - NDW
                        evenodd = ifelse((ix + iy + iz + it) % 2 == 0,true,false)
                        if evenodd == iseven
                            @simd for i1=1:NC
                                #println(a.f[i1,i2,i3,i4,i5,i6],"\t",b.f[i1,i2,i3,i4,i5,i6] )
                                c.f[i1,i2,i3,i4,i5,i6] += α*a.f[i1,i2,i3,i4,i5,i6] + β*b.f[i1,i2,i3,i4,i5,i6]
                            end
                        end
                    end
                end
            end
        end
    end
    return
end

function add_fermion!(c::WilsonFermion_4D_wing{NC,NDW},α::Number,a::T1,iseven::Bool) where {NC,NDW,T1 <: Abstractfermion,T2 <: Abstractfermion}#c += alpha*a + beta*b
    n1,n2,n3,n4,n5,n6 = size(c.f)

    @inbounds for i6=1:n6
        for i5=1:n5
            it = i5 -NDW
            for i4=1:n4
                iz = i4 -NDW
                for i3=1:n3
                    iy = i3 - NDW
                    for i2=1:n2
                        ix = i2 - NDW
                        evenodd = ifelse((ix + iy + iz + it) % 2 == 0,true,false)
                        if evenodd == iseven
                            @simd for i1=1:NC
                                #println(a.f[i1,i2,i3,i4,i5,i6],"\t",b.f[i1,i2,i3,i4,i5,i6] )
                                c.f[i1,i2,i3,i4,i5,i6] += α*a.f[i1,i2,i3,i4,i5,i6]
                            end
                        end
                    end
                end
            end
        end
    end
    return
end

function WWx!(xout::T,U::Array{G,1},x::T,A)  where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields} #(1 - K^2 Teo Toe) xe
    iseven = true
    isodd = false
    temp = A._temporary_fermi[7]#temps[4]
    temp2 = A._temporary_fermi[6]#temps[4]

    #println("Wx")
    #@time Wx!(xout,U,x,A)
    clear_fermion!(xout,iseven)

    #Tx!(temp,U,x,A)
    Toex!(temp,U,x,A,iseven) #Toe
    #Tx!(temp2,U,temp,A)
    Toex!(temp2,U,temp,A,isodd) #Teo

    #set_wing_fermion!(temp,A.boundarycondition)
    #add_fermion!(xout,1,x,-1,temp2)
    add_fermion!(xout,1,x,-1,temp2,iseven)
    set_wing_fermion!(xout,A.boundarycondition,iseven)

    #Wx!(xout,U,x,A)
    return

    #Tx!(temp2,U,x,A) #Toe

    #Toex!(temp2,U,x,A,iseven) #Toe

    #Toex!(temp,U,temp2,A,isodd) #Teo
    #Tx!(temp,U,temp2,A) #Toe

    Tx!(temp,U,x,A) #Toe

    add_fermion!(xout,1,x,-1,temp)
    #add_fermion!(xout,1,x,-1,temp,iseven)

    iseven = true
    set_wing_fermion!(xout,A.boundarycondition)

    return
end

function WWdagx!(xout::T,U::Array{G,1},x::T,A)  where  {T <: WilsonFermion_4D_wing,G <: AbstractGaugefields} #(1 - K^2 Teo Toe) xe
    iseven = true
    isodd = false
    temp = A._temporary_fermi[7]#temps[4]
    temp2 = A._temporary_fermi[6]#temps[4]

    clear_fermion!(xout)

    #Tx!(temp,U,x,A)
    Tdagoex!(temp,U,x,A,iseven) #Toe
    #Tx!(temp2,U,temp,A)
    Tdagoex!(temp2,U,temp,A,isodd) #Teo

    #set_wing_fermion!(temp,A.boundarycondition)
    #add_fermion!(xout,1,x,-1,temp2)
    add_fermion!(xout,1,x,-1,temp2,iseven)
    set_wing_fermion!(xout,A.boundarycondition)

    return
end

function clear_fermion!(a::WilsonFermion_4D_wing{NC,NDW} ,iseven) where {NC,NDW}
    n1,n2,n3,n4,n5,n6 = size(a.f)
    @inbounds for i6=1:n6
        for i5=1:n5
            it = i5-NDW
            for i4=1:n4
                iz = i4-NDW
                for i3=1:n3
                    iy = i3 - NDW
                    for i2=1:n2
                        ix = i2 - NDW
                        evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                        if evenodd == iseven
                            @simd for i1=1:NC
                                a.f[i1,i2,i3,i4,i5,i6]= 0
                            end
                        end
                    end
                end
            end
        end
    end
end

function LinearAlgebra.mul!(x::WilsonFermion_4D_wing{NC,NDW},A::TA) where {TA <: AbstractMatrix, NC,NDW}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            x[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            x[ic,ix,iy,iz,it,2] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            x[ic,ix,iy,iz,it,3] = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            x[ic,ix,iy,iz,it,4] = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(x::WilsonFermion_4D_wing{NC,NDW},A::TA,iseven::Bool) where {TA <: AbstractMatrix, NC,NDW}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                        evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                        if evenodd == iseven
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            x[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            x[ic,ix,iy,iz,it,2] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            x[ic,ix,iy,iz,it,3] = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            x[ic,ix,iy,iz,it,4] = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4
                        end
                    end
                end
            end
        end
    end

end

#function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC,NDW},A::TA,x::WilsonFermion_4D_wing{NC}) where {TA <: AbstractMatrix, NC,NDW}
function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC,NDW},A::TA,x::Abstractfermion) where {TA <: AbstractMatrix, NC,NDW}

    NX = xout.NX
    NY = xout.NY
    NZ = xout.NZ
    NT = xout.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            xout[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            xout[ic,ix,iy,iz,it,2] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            xout[ic,ix,iy,iz,it,3] = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            xout[ic,ix,iy,iz,it,4] = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC,NDW},A::TA,x::WilsonFermion_4D_wing{NC},iseven::Bool) where {TA <: AbstractMatrix, NC,NDW}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                        evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                        if evenodd == iseven
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            xout[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[1,2]*e2+A[1,3]*e3+A[1,4]*e4
                            xout[ic,ix,iy,iz,it,2] = A[2,1]*e1+A[2,2]*e2+A[2,3]*e3+A[2,4]*e4
                            xout[ic,ix,iy,iz,it,3] = A[3,1]*e1+A[3,2]*e2+A[3,3]*e3+A[3,4]*e4
                            xout[ic,ix,iy,iz,it,4] = A[4,1]*e1+A[4,2]*e2+A[4,3]*e3+A[4,4]*e4
                        end

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC,NDW},x::WilsonFermion_4D_wing{NC},A::TA) where {TA <: AbstractMatrix, NC,NDW}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            xout[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[2,1]*e2+A[3,1]*e3+A[4,1]*e4
                            xout[ic,ix,iy,iz,it,2] = A[1,2]*e1+A[2,2]*e2+A[3,2]*e3+A[4,2]*e4
                            xout[ic,ix,iy,iz,it,3] = A[1,3]*e1+A[2,3]*e2+A[3,3]*e3+A[4,3]*e4
                            xout[ic,ix,iy,iz,it,4] = A[1,4]*e1+A[2,4]*e2+A[3,4]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC,NDW},x::AbstractFermionfields_4D{NC},A::TA) where {TA <: AbstractMatrix, NC,NDW}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            xout[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[2,1]*e2+A[3,1]*e3+A[4,1]*e4
                            xout[ic,ix,iy,iz,it,2] = A[1,2]*e1+A[2,2]*e2+A[3,2]*e3+A[4,2]*e4
                            xout[ic,ix,iy,iz,it,3] = A[1,3]*e1+A[2,3]*e2+A[3,3]*e3+A[4,3]*e4
                            xout[ic,ix,iy,iz,it,4] = A[1,4]*e1+A[2,4]*e2+A[3,4]*e3+A[4,4]*e4

                    end
                end
            end
        end
    end

end

function LinearAlgebra.mul!(xout::WilsonFermion_4D_wing{NC,NDW},x::WilsonFermion_4D_wing{NC},A::TA,iseven::Bool) where {TA <: AbstractMatrix, NC,NDW}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT

    #n6 = size(x.f)[6]
    #f = zeros(ComplexF64,4)
    #e = zeros(ComplexF64,4)

    for ic=1:NC
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                        evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                        if evenodd == iseven
                            e1 = x[ic,ix,iy,iz,it,1]
                            e2 = x[ic,ix,iy,iz,it,2]
                            e3 = x[ic,ix,iy,iz,it,3]
                            e4 = x[ic,ix,iy,iz,it,4]

                            xout[ic,ix,iy,iz,it,1] = A[1,1]*e1+A[2,1]*e2+A[3,1]*e3+A[4,1]*e4
                            xout[ic,ix,iy,iz,it,2] = A[1,2]*e1+A[2,2]*e2+A[3,2]*e3+A[4,2]*e4
                            xout[ic,ix,iy,iz,it,3] = A[1,3]*e1+A[2,3]*e2+A[3,3]*e3+A[4,3]*e4
                            xout[ic,ix,iy,iz,it,4] = A[1,4]*e1+A[2,4]*e2+A[3,4]*e3+A[4,4]*e4
                        end
                    end
                end
            end
        end
    end

end

function LinearAlgebra.dot(a::WilsonFermion_4D_wing{NC,NDW},b::WilsonFermion_4D_wing{NC,NDW},iseven::Bool) where {NC,NDW}
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX
    NG = a.NG

    c = 0.0im
    @inbounds for α=1:NG
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                        if evenodd == iseven
                            @simd for ic=1:NC
                                c+= conj(a[ic,ix,iy,iz,it,α])*b[ic,ix,iy,iz,it,α]
                            end
                        end
                    end
                end
            end
        end
    end
    return c
end

#Overwrite Y with X*a + Y*b, where a and b are scalars. Return Y.
function LinearAlgebra.axpby!(a::Number, X::T, b::Number, Y::WilsonFermion_4D_wing{NC,NDW},iseven::Bool) where {NC,NDW,T <: AbstractFermionfields_4D}
    n1,n2,n3,n4,n5,n6 = size(Y.f)
    #println("axpby")

    @inbounds for i6=1:n6
        for i5=1:n5
            it = i5+NDW
            for i4=1:n4
                iz = i4+NDW
                for i3=1:n3
                    iy = i3+NDW
                    for i2=1:n2
                        ix = i2+NDW
                        evenodd = ifelse((ix+iy+iz+it) % 2 == 0,true,false)
                        if evenodd == iseven
                            @simd for i1=1:NC
                                Y.f[i1,i2,i3,i4,i5,i6] = a*X.f[i1,i2,i3,i4,i5,i6] + b*Y.f[i1,i2,i3,i4,i5,i6]
                            end
                        end
                    end
                end
            end
        end
    end
end

#=
function set_wing_fermion!(a::WilsonFermion_4D_wing{NC},boundarycondition) where NC
    NT = a.NT
    NZ = a.NZ
    NY = a.NY
    NX = a.NX

    #!  X-direction
    for ialpha=1:4
        for it=1:NT
            for iz = 1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,0,iy,iz,it,ialpha] = boundarycondition[1]*a[k,NX,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for k=1:NC
                        a[k,NX+1,iy,iz,it,ialpha] =boundarycondition[1]*a[k,1,iy,iz,it,ialpha]
                    end
                end
            end
        end
    end

    #Y-direction
    for ialpha = 1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,0,iz,it,ialpha] =boundarycondition[2]*a[k,ix,NY,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        for it=1:NT
            for iz=1:NZ
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,NY+1,iz,it,ialpha] = boundarycondition[2]*a[k,ix,1,iz,it,ialpha]
                    end
                end
            end
        end
    end

    for ialpha=1:4
        # Z-direction
        for it=1:NT
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,0,it,ialpha] = boundarycondition[3]*a[k,ix,iy,NZ,it,ialpha]
                        a[k,ix,iy,NZ+1,it,ialpha] = boundarycondition[3]*a[k,ix,iy,1,it,ialpha]

                    end
                end
            end
        end

        #T-direction
        for iz=1:NZ
            for iy=1:NY
                for ix=1:NX
                    @simd for k=1:NC
                        a[k,ix,iy,iz,0,ialpha] = boundarycondition[4]*a[k,ix,iy,iz,NT,ialpha]
                        a[k,ix,iy,iz,NT+1,ialpha] =boundarycondition[4]*a[k,ix,iy,iz,1,ialpha]
                    end
                end
            end
        end

    end

end

=#

"""
c--------------------------------------------------------------------------c
c     y = gamma_5 * x
c     here
c                  ( -1       )
c        GAMMA5 =  (   -1     )
c                  (     +1   )
c                  (       +1 )
c--------------------------------------------------------------------------c
    """
    function mul_γ5x!(y::WilsonFermion_4D_wing{NC},x::WilsonFermion_4D_wing{NC}) where NC
        n1,n2,n3,n4,n5,n6 = size(x.f)
        @inbounds for i6=1:n6
            for i5=1:n5
                #it = i5+NDW
                for i4=1:n4
                    #iz = i4+NDW
                    for i3=1:n3
                        #iy = i3+NDW
                        for i2=1:n2
                            #ix = i2+NDW
                            @simd for ic=1:NC
                                y.f[i1,i2,i3,i4,i5,i6]=x.f[i1,i2,i3,i4,i5,i6]*ifelse(i6 <= 2,-1,1)
                            end
                        end
                    end
                end
            end
        end

        #=
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        for ig=1:4
            for ic=1:NC
                for it=1:NT
                    for iz=1:NZ
                        for iy=1:NY
                            for ix=1:NX
                                @simd for ic=1:NC
                                    y[ic,ix,iy,iz,it,ig] =x[ic,ix,iy,iz,it,ig]*ifelse(ig <= 2,-1,1)
                                end
                            end
                        end
                    end
                end
            end
        end
        =#

    end

    function apply_γ5!(x::WilsonFermion_4D_wing{NC}) where NC
        n1,n2,n3,n4,n5,n6 = size(x.f)
        #println("axpby")

        @inbounds for i6=1:n6
            for i5=1:n5
                #it = i5+NDW
                for i4=1:n4
                    #iz = i4+NDW
                    for i3=1:n3
                        #iy = i3+NDW
                        for i2=1:n2
                            #ix = i2+NDW
                            @simd for ic=1:NC
                                x.f[i1,i2,i3,i4,i5,i6]=x.f[i1,i2,i3,i4,i5,i6]*ifelse(i6 <= 2,-1,1)
                            end
                        end
                    end
                end
            end
        end

    end

    function mul_1plusγ5x!(y::WilsonFermion_4D_wing{NC},x::WilsonFermion_4D_wing{NC})  where NC#(1+gamma_5)/2
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        #NC = x.NC
        for ic=1:NC
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        @simd for ix=1:NX
                            y[ic,ix,iy,iz,it,1] = 0#-1*x[ic,ix,iy,iz,it,1]
                            y[ic,ix,iy,iz,it,2] = 0#-1*x[ic,ix,iy,iz,it,2]
                            y[ic,ix,iy,iz,it,3] = x[ic,ix,iy,iz,it,3]
                            y[ic,ix,iy,iz,it,4] = x[ic,ix,iy,iz,it,4]
                        end
                    end
                end
            end
        end
    end

    function mul_1plusγ5x_add!(y::WilsonFermion_4D_wing{NC},x::WilsonFermion_4D_wing{NC},factor) where NC#x = x +(1+gamma_5)/2
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        #NC = x.NC
        for ic=1:NC
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        @simd for ix=1:NX
                            y[ic,ix,iy,iz,it,3] += factor*x[ic,ix,iy,iz,it,3]
                            y[ic,ix,iy,iz,it,4] += factor*x[ic,ix,iy,iz,it,4]
                        end
                    end
                end
            end
        end
    end

    function mul_1minusγ5x!(y::WilsonFermion_4D_wing{NC},x::WilsonFermion_4D_wing{NC}) where NC#(1-gamma_5)/2
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        #NC = x.NC
        for ic=1:NC
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        @simd for ix=1:NX
                            y[ic,ix,iy,iz,it,1] = x[ic,ix,iy,iz,it,1]
                            y[ic,ix,iy,iz,it,2] = x[ic,ix,iy,iz,it,2]
                            y[ic,ix,iy,iz,it,3] = 0#x[ic,ix,iy,iz,it,3]
                            y[ic,ix,iy,iz,it,4] = 0#x[ic,ix,iy,iz,it,4]
                        end
                    end
                end
            end
        end
    end

    function mul_1minusγ5x_add!(y::WilsonFermion_4D_wing{NC},x::WilsonFermion_4D_wing{NC},factor) where NC#+(1-gamma_5)/2
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        #NC = x.NC
        for ic=1:NC
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        @simd for ix=1:NX
                            y[ic,ix,iy,iz,it,1] += factor*x[ic,ix,iy,iz,it,1]
                            y[ic,ix,iy,iz,it,2] += factor*x[ic,ix,iy,iz,it,2]
                        end
                    end
                end
            end
        end
    end

    """
               (       -i )              (       -1 )
     GAMMA1 =  (     -i   )     GAMMA2 = (     +1   )
               (   +i     )              (   +1     )
               ( +i       )              ( -1       )
    """

    function mul_1minusγ1x!(y::WilsonFermion_4D_wing{NC},x) where NC#(1-gamma_5)/2
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for it=1:NT
            for iz=1:NZ
                for iy=1:NY
                    @simd for ix=1:NX
                        for ic=1:NC
                            v1 = x[ic,ix,iy,iz,it,1] + im*x[ic,ix,iy,iz,it,4]
                            v2 = x[ic,ix,iy,iz,it,2] + im*x[ic,ix,iy,iz,it,3]
                            v3 = x[ic,ix,iy,iz,it,3] - im*x[ic,ix,iy,iz,it,2]
                            v4 = x[ic,ix,iy,iz,it,4] - im*x[ic,ix,iy,iz,it,1]
                            y[ic,ix,iy,iz,it,1] = v1
                            y[ic,ix,iy,iz,it,2] = v2
                            y[ic,ix,iy,iz,it,3] = v3
                            y[ic,ix,iy,iz,it,4] = v4
                        end
                    end
                end
            end

        end
    end

    function mul_1plusγ1x!(y::WilsonFermion_4D_wing{NC},x) where NC#(1-gamma_5)/2
        NX = y.NX
        NY = y.NY
        NZ = y.NZ
        NT = y.NT
        #NC = x.NC
        @inbounds for ic=1:NC
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        @simd for ix=1:NX
                            v1 = x[ic,ix,iy,iz,it,1] - im*x[ic,ix,iy,iz,it,4]
                            v2 = x[ic,ix,iy,iz,it,2] - im*x[ic,ix,iy,iz,it,3]
                            v3 = x[ic,ix,iy,iz,it,3] + im*x[ic,ix,iy,iz,it,2]
                            v4 = x[ic,ix,iy,iz,it,4] + im*x[ic,ix,iy,iz,it,1]
                            y[ic,ix,iy,iz,it,1] = v1
                            y[ic,ix,iy,iz,it,2] = v2
                            y[ic,ix,iy,iz,it,3] = v3
                            y[ic,ix,iy,iz,it,4] = v4
                        end
                    end
                end
            end
        end
    end


text:
include("lib_for_robot.jl")
function task16(r)
    x,y = tostart(r)
    side=Ost
    while !isborder(r,Nord)||!isborder(r,Ost)
        movetoside(r,side)
        if !isborder(r,Nord)
            move!(r,Nord)
        end
        side=reverse(side)
    end
    tostart(r)
    todefault(r,Ost,x)
    todefault(r,Nord,y)
end


text:
module Callbacks

import .. Pages

import ..HTTP.WebSockets: WebSocket

export callbacks, connections, conditions

export Callback

mutable struct Callback
    callback::Function
    name::String

    function Callback(callback,name)
        cb = new(callback,name)
        callbacks[name] = cb
        cb
    end
end
const callbacks = Dict{String,Callback}() # name => args -> f(args...)
const connections = Dict{String,Dict{String,WebSocket}}()
const conditions = Dict{String,Condition}()

# Callback to notify the Server that a new page is loaded and its WebSocket is ready.
Callback("connected") do client, route, id, data
    if haskey(connections,route)
        connections[route][id] = client
    else
        connections[route] = Dict(id => client)
    end
    println("New connection established (ID: $(id)).")
    notify(conditions["connected"])
end

# Callback used to cleanup when the browser navigates away from the page.
Callback("unloaded") do client, route, id, data
    delete!(connections[route],id)
end

# Callback used to message a specified connected browser.
Callback("message") do client, route, id, data
    message(route,data["args"])
end

# Callback used to broadcast to all connected browsers.
Callback("broadcast") do client, route, id, data
    Pages.broadcast(route,data["args"])
end

# Callback used for blocking Julia control flow until notified by the WebSocket.
Callback("notify") do client, route, id, data
    name = data["args"]
    if haskey(conditions,name)
        notify(conditions[name])
    else
        @warn "Condition '$name' was not found."
    end
end

end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

using SlimOptim, LinearAlgebra, JOLI, TestImages, PyPlot
import TestImages: Gray

img = Float32.(Gray.(testimage("lighthouse.png")[1:2:end, 129:2:end-128]))

n = 256
k = 4
# Sparse in wavelet domain
W = joDWT(n, n; DDT=Float32, RDT=Float32)
# Or with curvelet ifi nstalled
# W = joCurvelet2D(128, 128; DDT=Float32, RDT=Float32)
A = vcat([joRomberg(n, n; DDT=Float32, RDT=Float32) for i=1:k]...)

# Make noisy data
imgn= img .+ .01f0*randn(Float32, size(img))
b = A*vec(imgn)

# setup bregamn
opt = bregman_options(maxIter=200, verbose=2, quantile=.5, alpha=1, antichatter=true)
opt2 = bregman_options(maxIter=200, verbose=2, quantile=.5, alpha=1, antichatter=true, spg=true)

sol = bregman(A, W, zeros(Float32, n*n), b, opt)
sol2 = bregman(A, W, zeros(Float32, n*n), b, opt2)

figure()
subplot(121)
plot(sol.ϕ_trace, label="std");plot(sol2.ϕ_trace, label="spg");
legend()
subplot(122)
loglog(sol.r_trace, label="std");loglog(sol2.r_trace, label="spg");
legend()

figure();
subplot(221)
imshow(img, cmap="Greys", vmin=-.5, vmax=.5)
title("True")
subplot(222)
imshow(reshape(b, n, k*n), cmap="Greys", vmin=-.5, vmax=.5)
title("Measurment")
subplot(223)
imshow(reshape(sol2.x, n, n), cmap="Greys", vmin=-.5, vmax=.5)
title("recovered")
subplot(224)
imshow(img - reshape(sol2.x, n, n), cmap="Greys", vmin=-.05, vmax=.05)
title("Difference x10")
tight_layout()


text:
function _info(model::Optimizer, key::MOI.VariableIndex)
    if haskey(model.variable_info, key)
        return model.variable_info[key]
    end
    throw(MOI.InvalidIndex(key))
    return
end

function _make_var(model::Optimizer, variable::Variable)
    # Initialize `VariableInfo` with a dummy `VariableIndex` and a column,
    # because we need `add_item` to tell us what the `VariableIndex` is.
    index = CleverDicts.add_item(
        model.variable_info,
        VariableInfo(MOI.VariableIndex(0), variable),
    )
    _info(model, index).index = index
    return index
end

function _make_var(
    model::Optimizer,
    variable::Variable,
    set::MOI.AbstractScalarSet,
)
    index = _make_var(model, variable)
    S = typeof(set)
    return index, MOI.ConstraintIndex{MOI.VariableIndex, S}(index.value)
end

function _make_vars(model::Optimizer, variables::Vector{<:Variable})
    # Barely used, because add_constrained_variables may have variable sets (except for AbstractVectorSet).
    # Only implemented in the unconstrained case.
    indices = Vector{MOI.VariableIndex}(undef, length(variables))
    for i in 1:length(variables)
        indices[i] = CleverDicts.add_item(
            model.variable_info,
            VariableInfo(MOI.VariableIndex(0), variables[i]),
        )
        _info(model, indices[i]).index = indices[i]
    end
    return indices
end

function _sanitise_bounds(lb::Real, ub::Real, T)
    if lb === nothing
        lb = typemin(T)
    end
    if ub === nothing
        ub = typemax(T)
    end
    return lb, ub
end

function _make_floatvar(
    model::Optimizer,
    set::MOI.AbstractScalarSet;
    lb::Union{Nothing, Float64}=nothing,
    ub::Union{Nothing, Float64}=nothing,
)
    v = if lb === nothing && ub === nothing
        FloatVar(
            (Store,),
            model.inner,
        )
    else
        lb_, ub_ = _sanitise_bounds(lb, ub, Float64)
        FloatVar(
            (Store, jdouble, jdouble),
            model.inner,
            lb_,
            ub_,
        )
    end

    vindex, cindex = _make_var(model, v, set)
    _info(model, vindex).type = CONTINUOUS
    _info(model, vindex).lb = lb
    _info(model, vindex).ub = ub
    return vindex, cindex
end

function _make_intvar(
    model::Optimizer,
    set::MOI.AbstractScalarSet;
    lb::Int32=typemin(Int32),
    ub::Int32=typemax(Int32),
)
    v = if lb === nothing && ub === nothing
        FloatVar(
            (Store,),
            model.inner,
        )
    else
        lb_, ub_ = _sanitise_bounds(lb, ub, Int32)
        IntVar(
            (Store, jint, jint),
            model.inner,
            lb_,
            ub_,
        )
    end

    vindex, cindex = _make_var(model, v, set)
    _info(model, vindex).type = INTEGER
    _info(model, vindex).lb = lb
    _info(model, vindex).ub = ub
    return vindex, cindex
end

function _make_boolvar(model::Optimizer, set::MOI.AbstractScalarSet)
    vindex, cindex = _make_var(model, BoolVar((Store,), model.inner), set)
    _info(model, vindex).type = BINARY
    return vindex, cindex
end

function supports_add_constrained_variable(
    ::Optimizer,
    ::Type{F},
) where {
    F <: Union{
        MOI.EqualTo{Float64},
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.Interval{Float64},
        MOI.EqualTo{Int},
        MOI.LessThan{Int},
        MOI.GreaterThan{Int},
        MOI.Interval{Int},
        MOI.ZeroOne,
        MOI.Integer,
    },
}
    return true
end

function MOI.add_variable(model::Optimizer)
    v = FloatVar((Store,), model.inner)
    vindex = _make_var(model, v)
    _info(model, vindex).type = CONTINUOUS
    _info(model, vindex).lb = lb
    _info(model, vindex).ub = ub
    return vindex
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.GreaterThan{T},
) where {T <: Real}
    return _make_floatvar(model, set, lb=set.lower)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.LessThan{T},
) where {T <: Real}
    return _make_floatvar(model, set, lb=set.upper)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.EqualTo{T},
) where {T <: Real}
    return _make_floatvar(model, set, lb=set.value, ub=set.value)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.Interval{T},
) where {T <: Real}
    return _make_floatvar(model, set, lb=set.lower, ub=set.upper)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.GreaterThan{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.lower)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.LessThan{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.upper)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.EqualTo{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.value, ub=set.value)
end

function MOI.add_constrained_variable(
    model::Optimizer,
    set::MOI.Interval{T},
) where {T <: Integer}
    return _make_intvar(model, set, lb=set.lower, ub=set.upper)
end

function MOI.add_constrained_variable(model::Optimizer, set::MOI.ZeroOne)
    return _make_boolvar(model, set)
end

function MOI.add_constrained_variable(model::Optimizer, set::MOI.Integer)
    return _make_intvar(model, set)
end

function MOI.is_valid(model::Optimizer, v::MOI.VariableIndex)
    return haskey(model.variable_info, v)
end

# No deletion allowed.


text:
function check_in_interval(num,interval)
    (num >= interval[1]) & (num <= interval[2])
end

function compute_posterior(vec_num,alpha)
    quantile(vec_num,[alpha/2,1-alpha/2])
end

function distanceFunction(ytilde, y)
    return mean((ytilde[3,:]-y[3,:]) .^ 2)
end

function identity_mapping(y)
    y
end

function compute_norm(y,yhat;p=2)
    norm(y[2:3,2:end] - yhat[2:3,2:end], p)/size(y)[1]
end

function compute_l2_norm(y,yhat)
    sol1=hcat(y[:]...)'
    sol2=hcat(yhat[:]...)'
    l2_error = sqrt(sum(sum((sol1.-sol2).^2,dims=2)*0.01,dims=1))[1]
end

function compute_l2_norm_decay(y,yhat)
    sol1=hcat(y[:]...)'
    sol2=hcat(yhat[:]...)'
    T=size(sol1)[1]
    l2_error = sqrt(sum(sum((sol1.-sol2).^2,dims=2).*exp.(-0.001*(1:1:T))*0.01,dims=1))[1]
end

function compute_full_norm(y,yhat;p=Inf)
    norm(y - yhat, p)/size(y)[1]
end

function compute_at_steps(y,yhat;index=3:3:18,p=Inf)
    norm((y .- yhat)[2:3,index], p)/length(index)
end

function random_walk(p_prev, sd)
    # return proposed p, log(q(p_n+1|p_n))

    p_propose = map( (x,y) -> x + rand((Normal(0,y))) ,p_prev,sd)

    log_pdf = map( (x,y,z) -> logpdf((Normal(x,z)), y ), p_prev, p_propose, sd )
    return p_propose, log_pdf
end

function proposal_Normal_density(p1, p2, sd)
    return map( (x,y,z) -> pdf(Normal(x,z),y), p1, p2, sd )
end

function proposal_LN(log_p_prev,sd)
    log_p_prev .+ rand.(Normal.(0,sd))
end

function proposalRatio_LN(p_prev,p_cand,sd)
    0.5
end

function proposal_Gamma(log_p_prev,k)
    #TODO
end

function proposalRatio_Gamma(p_prev,p_cand,k)
    # TODO
end

function proposal_Normal(p_prev,sd)
    p_prev .+ rand.(Normal.(0,sd))
end

function proposalRatio_Normal(p_prev,p_cand,sd)
    1
end

function systematic_resample(weights)
    N = length(weights)
    partitions = rand()/N .+ (0:1:N-1)./N
    idx = zeros(Int,N)
    Q = cumsum(weights)
    i=1;j=1;
    while i<=N
        if partitions[i] < Q[j]
            idx[i]=j
            i +=1
        else
            j+=1
        end
    end
    return idx
end

function CPUtoc_modified(verbose)
    t = CPUTime.CPUtoq()
    if verbose
        println("elpased CPU time: ",t," seconds")
    end
    return t
end

function show_calibration(simulation)
    @printf("Calibration: ABC %f, ABC MCMC %f, ABC SMC %f\n", mean(simulation[2],dims=1)[1], mean(simulation[2],dims=1)[2],mean(simulation[2],dims=1)[3])
    @printf("CPU_time: ABC %f, ABC MCMC %f, ABC SMC %f\n", mean(simulation[4],dims=1)[1], mean(simulation[4],dims=1)[2],mean(simulation[4],dims=1)[3])
    @printf("CPU_time: ABC %f, ABC MCMC %f, ABC SMC %f\n", median(simulation[4],dims=1)[1], median(simulation[4],dims=1)[2],median(simulation[4],dims=1)[3])
    @printf("ESS: ABC %f, ABC MCMC %f, ABC SMC %f\n", median(simulation[5],dims=1)[1], mean(simulation[5],dims=1)[2],mean(simulation[5],dims=1)[3])
    @printf("ESS: ABC %f, ABC MCMC %f, ABC SMC %f\n", mean(simulation[5],dims=1)[1], median(simulation[5],dims=1)[2],median(simulation[5],dims=1)[3])
    @printf("ESS Time: ABC %f, ABC MCMC %f, ABC SMC %f\n", mean(simulation[6],dims=1)[1], mean(simulation[6],dims=1)[2],mean(simulation[6],dims=1)[3])
    @printf("ESS Time: ABC %f, ABC MCMC %f, ABC SMC %f\n", median(simulation[6],dims=1)[1], median(simulation[6],dims=1)[2],median(simulation[6],dims=1)[3])
end


text:
using ProximalAlgorithms

const ForwardBackwardSolver = Union{
    ProximalAlgorithms.ForwardBackward,
    ProximalAlgorithms.ZeroFPR,
    ProximalAlgorithms.PANOC,
}

const default_solver = ProximalAlgorithms.PANOC


text:
==(x::Rule, y::Rule) = x.p == y.p && x.q == y.q

function unique(v::Vector{Rule})
    out = Vector{Rule}(0)
    for i = 1:length(v)
        if !in(v[i], out)
            push!(out, v[i])
        end
    end
    return out
end

"""
This function converts a data matrix `X` into a
vector of transactions, where each row in the orginal
matrix is now an element of the vector. This is used
for pre-processing data prior to applying the a-priori
algorithm for association rule mining.
"""
function make_transactions(X::Array{Any,2})
    n, p = size(X)
    X = map(string, X)

    out = Array{Array{String, 1}, 1}(n)
    for i = 1:n
        out[i] = convert(Array{String, 1}, X[i, :])
    end
    return out
end

function make_transactions(X::DataFrame)
    n, p = size(X)
    Xstr = Array{String}(n, p)
    for j = 1:p
        Xstr[:, j] = map(string, X[:, j])
    end

    out = Array{Array{String, 1}, 1}(n)
    for i = 1:n
        out[i] = convert(Array{String, 1}, Xstr[i, :])
    end
    return out
end

function wrap_line(v, maxchar)
    if isempty(v)
        res = "{}"
    else
        res = replace(string(v), "[", "{")
        res = replace(res, "]", "}")
    end
    res = replace(res, r"\"|String|Any", "")
    n = length(res)
    out = ""
    str = 1
    stp = maxchar
    if n > maxchar
        keep_splitting = true

        while keep_splitting
            if stp < n
                # this regex finds the last comma
                chk = match(r",[^,]*$", res[str:stp])

                # In cases where we have single elements longer than
                # maxchar, `offset` will be of type Void.
                if isa(chk, RegexMatch)
                    offset = chk.offset
                elseif isa(chk, Void)
                    offset = maxchar
                end

                stp = str + offset - 1
            else
                stp = n
                keep_splitting = false
            end
            if keep_splitting
                out = string(out, res[str:stp], "\n")
            else
                n_char = maxchar - (stp - str) - 1
                padding = rpad(" ", n_char, " ")
                out = string(out, res[str:stp], padding)
            end
            str = stp + 1
            stp = str + maxchar - 1
            stp = stp > n ? n : stp
        end
    elseif n ≤ maxchar
        out = res
    end
    return out
end

# v1 = ["this here", "is", "my", "vector", "and", "this", "vector", "contains", "quite", "a", "number", "of", "variable-length", "terms", "such", "as", "these."]
# v2 = ["a", "number", "of", "variable-length", "terms", "such", "as", "these."]
# w1 = wrap_line(v1, 25)
# w2 = wrap_line(v2, 25)
# println(w1)

# Given two vectors of split lines (p and q), where the vectors
# are of the same length, this function prints out the lines and
# aligns them appropriately.
function print_equal_length(p_lines, q_lines, colwidth, i = 1, j = 1)
    np = length(p_lines)

    while i ≤ np
        if i ≠ np
            print(rpad(p_lines[i], colwidth, " "))
            println("    ", q_lines[j])
        else
            print(rpad(p_lines[i], colwidth, " "))
            print(" => ", rpad(q_lines[j], colwidth+1, " "))
        end
        i += 1
        j += 1
    end
end

function print_wrapped_lines(p_line, q_line, colwidth)
    p_lines = split(p_line, '\n', keep = true)
    q_lines = split(q_line, '\n', keep = true)
    np = length(p_lines)
    nq = length(q_lines)

    if np == nq
        print_equal_length(p_lines, q_lines, colwidth)
    elseif np > nq
        i = 1
        ndiff = np - nq

        while ndiff > 0
            println(rpad(p_lines[i], colwidth, " "))
            ndiff -= 1
            i += 1
        end
        print_equal_length(p_lines, q_lines, colwidth, i, 1)
    elseif np < nq
        j = 1
        ndiff = nq - np
        padding = lpad("", colwidth+4, " ")

        while ndiff > 0
            println(padding, rpad(q_lines[j], colwidth, " "))
            ndiff -= 1
            j += 1
        end
        print_equal_length(p_lines, q_lines, colwidth, 1, j)
    end
end

# print_wrapped_lines(wrap_line(v1, 15), wrap_line(v2, 15), 20)

function prettyprint_rulestats(r::Rule, plen, qlen, colwidth)
    p = wrap_line(r.p, colwidth)
    q = wrap_line(r.q, colwidth)
    fmt = "%0.4f"

    # print antecedent and consequent
    print_wrapped_lines(p, q, colwidth)

    # print the rule's stats
    println("$(sprintf1(fmt, r.supp)) $(sprintf1(fmt, r.conf)) $(sprintf1(fmt, r.lift))", )
end

# prettyprint_rulestats(rules[1], 10, 10, 30)

function display(rs::Vector{Rule}, colwidth = 30)
    if isempty(rs)
        return nothing
    end

    n = length(rs)

    maxlen_p = maximum(map(x -> length(wrap_line(x.p, colwidth)), rs))
    maxlen_q = maximum(map(x -> length(wrap_line(x.q, colwidth)), rs))

    # When column width greater than the max string length,
    # we use colwidth as the width of our columns.
    maxlen_p = maxlen_p < colwidth ? maxlen_p : colwidth
    maxlen_q = maxlen_q < colwidth ? maxlen_q : colwidth

    colwidth = max(maxlen_p, maxlen_q)

    println(rpad("lhs", colwidth, " "), "    ", rpad("rhs", colwidth, " "), " supp.  conf.  lift")
    for i = 1:n
        prettyprint_rulestats(rs[i], maxlen_p, maxlen_q, colwidth)
    end
end

# display(rules)

"""
    dataset(dataset_name::AbstractString)

This is a convenience function whose only purpose is to load example datasets provided with this package.

### Arguments
* `dataset_name`: A string specifying the dataset to be loaded (e.g., "zaki_data", "adult", etc).

The returned value will be a dataframe.
"""
function dataset(dataset_name::AbstractString)
    filename = joinpath(dirname(@__FILE__), "..", "data", "$dataset_name.csv")

    if !isfile(filename)
        error(@sprintf "Unable to locate file %s \n" filename)
    else
        return readtable(filename)
    end
end