text:
module CommentBot

using Sockets
using GitHub
using HTTP
using Distributed
using Pkg
using Logging
using Dates
using JSON
using MbedTLS

import Pkg: TOML
import ..Registrator: post_on_slack_channel, pull_request_contents
import RegistryTools: RegBranch, RegisterParams
import Base: string
using ..Messaging

include("trigger_types.jl")
include("parse_comment.jl")
include("github_utils.jl")
include("verify_projectfile.jl")
include("param_types.jl")
include("approval.jl")

include("../management.jl")

const event_queue = Channel{RequestParams}(1024)
const CONFIG = Dict{String,Any}()
const httpsock = Ref{Sockets.TCPServer}()

function print_entry_log(rp::RequestParams{PullRequestTrigger})
    @info "Creating registration pull request" reponame=rp.reponame prid=rp.trigger_src.prid
end

function print_entry_log(rp::RequestParams{CommitCommentTrigger})
    @info "Creating registration pull request" reponame=rp.reponame sha=get_comment_commit_id(rp.evt)
end

function print_entry_log(rp::RequestParams{IssueTrigger})
    @info "Creating registration pull request" reponame=rp.reponame branch=rp.trigger_src.branch
end

get_trigger_id(rp::RequestParams{PullRequestTrigger}) = rp.trigger_src.prid
get_trigger_id(rp::RequestParams{IssueTrigger}) = get_prid(rp.evt.payload)
get_trigger_id(rp::RequestParams{CommitCommentTrigger}) = get_comment_commit_id(rp.evt)

tag_name(version, subdir) = subdir == "" ? "v$version" : splitdir(subdir)[end] * "-v$version"

function make_pull_request(pp::ProcessedParams, rp::RequestParams, rbrn::RegBranch, target_registry::Dict{String,Any})
    name = rbrn.name
    ver = rbrn.version
    brn = rbrn.branch
    subdir = rp.subdir

    subdir_string = subdir == "" ? "" : ", subdir=$subdir"
    @info("Creating pull request name=$name, ver=$ver, branch=$brn" * subdir_string)
    payload = rp.evt.payload
    creator = get_user_login(payload)
    reviewer = payload["sender"]["login"]
    @debug("Pull request creator=$creator, reviewer=$reviewer")

    trigger_id = get_trigger_id(rp)

    meta = JSON.json(Dict("request_type"=> string(rp),
                          "pkg_repo_name"=> rp.reponame,
                          "trigger_id"=> trigger_id,
                          "tree_sha"=> pp.tree_sha,
                          "version"=> string(ver),
                          "subdir"=> subdir))
    key = CONFIG["enc_key"]
    enc_meta = "<!-- " * bytes2hex(encrypt(MbedTLS.CIPHER_AES_128_CBC, key, meta, key)) * " -->"
    params = Dict("base"=>target_registry["base_branch"],
                  "head"=>brn,
                  "maintainer_can_modify"=>true)
    ref = get_html_url(rp.evt.payload)

    params["title"], params["body"] = pull_request_contents(;
        registration_type=get(rbrn.metadata, "kind", ""),
        package=name,
        repo=string(rp.evt.repository.html_url),
        user="@$creator",
        version=ver,
        commit=pp.sha,
        release_notes=rp.release_notes,
        reviewer="@$reviewer",
        reference=ref,
        meta=enc_meta,
    )

    repo = join(split(target_registry["repo"], "/")[end-1:end], "/")
    pr, msg = create_or_find_pull_request(repo, params, rbrn)
    tag = tag_name(ver, subdir)

    cbody = """
        Registration pull request $msg: [$(repo)/$(pr.number)]($(pr.html_url))

        After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

        This will be done automatically if the [Julia TagBot GitHub Action](https://github.com/marketplace/actions/julia-tagbot) is installed, or can be done manually through the github interface, or via:
        ```
        git tag -a $tag -m "<description of version>" $(pp.sha)
        git push origin $tag
        ```
        """

    if get(rbrn.metadata, "warning", nothing) !== nothing
        cbody *= """
            Also, note the warning: $(rbrn.metadata["warning"])
            This can be safely ignored. However, if you want to fix this you can do so. Call register() again after making the fix. This will update the Pull request.
            """
    end

    @debug(cbody)
    make_comment(rp.evt, cbody)
    return pr
end

string(::RequestParams{PullRequestTrigger}) = "pull_request"
string(::RequestParams{CommitCommentTrigger}) = "commit_comment"
string(::RequestParams{IssueTrigger}) = "issue"

function action(rp::RequestParams{T}, zsock::RequestSocket) where T <: RegisterTrigger
    if rp.target === nothing
        target_registry_name, target_registry = first(CONFIG["targets"])
    else
        filteredtargets = filter(x->(x[1]==rp.target), CONFIG["targets"])
        if length(filteredtargets) == 0
            msg = "Error: target $(rp.target) not found"
            @debug(msg)
            make_comment(rp.evt, msg)
            set_error_status(rp)
            return
        else
            target_registry_name, target_registry = filteredtargets[1]
        end
    end

    pp = ProcessedParams(rp)
    @info("Processing register event", reponame=rp.reponame, target_registry_name)
    try
        if pp.cparams.isvalid
            registry_deps=map(String, get(CONFIG, "registry_deps", String[]))
            regp = RegisterParams(pp.cloneurl,
                                  pp.project,
                                  pp.tree_sha;
                                  subdir=rp.subdir,
                                  registry=target_registry["repo"],
                                  registry_deps=registry_deps,
                                  push=true,
                                  )
            rbrn = sendrecv(zsock, regp; nretry=10)

            if rbrn === nothing || get(rbrn.metadata, "error", nothing) !== nothing
                if rbrn === nothing
                    msg = "ERROR: Registrator backend service unreachable"
                else
                    msg = "Error while trying to register: $(rbrn.metadata["error"])"
                end
                @debug(msg)
                make_comment(rp.evt, msg)
                set_error_status(rp)
            else
                make_pull_request(pp, rp, rbrn, target_registry)
                set_success_status(rp)
            end
        else
            @info("Error while processing event: $(repr(pp.cparams.error))")
            if pp.cparams.report_error
                msg = "Error while trying to register: $(repr(pp.cparams.error))"
                @debug(msg)
                make_comment(rp.evt, msg)
            end
            set_error_status(rp)
        end
    catch ex
        bt = get_backtrace(ex)
        @info("Unexpected error: $bt")
        raise_issue(rp.evt, rp.phrase, bt)
    end
    @info("Done processing register event", reponame=rp.reponame, target_registry_name)
    nothing
end

function comment_handler(event::WebhookEvent, phrase::RegexMatch)
    @debug("Received event for $(event.repository.full_name), phrase: $phrase")
    try
        rp = RequestParams(event, phrase)
        isa(rp.trigger_src, EmptyTrigger) && rp.cparams.error === nothing && return

        if rp.cparams.isvalid && rp.cparams.error === nothing
            print_entry_log(rp)
            put!(event_queue, rp)
            set_pending_status(rp)
        elseif rp.cparams.error !== nothing
            @info("Error while processing event: $(rp.cparams.error)")
            if rp.cparams.report_error
                msg = "Error while trying to register: $(rp.cparams.error)"
                @debug(msg)
                make_comment(event, msg)
            end
            set_error_status(rp)
        end
    catch ex
        bt = get_backtrace(ex)
        @info("Unexpected error: $bt")
        raise_issue(event, phrase, bt)
    end

    return HTTP.Messages.Response(200)
end

make_trigger(cfg=CONFIG) = Regex(cfg["trigger"] * "(.*)", "i")

function github_webhook(http_ip=CONFIG["http_ip"],
                        http_port=get(CONFIG, "http_port", parse(Int, get(ENV, "PORT", "8001"))); kwargs...)
    auth = get_jwt_auth()
    trigger = make_trigger()
    listener = GitHub.CommentListener(comment_handler, trigger; check_collab=false, auth=auth, secret=CONFIG["github"]["secret"])
    httpsock[] = Sockets.listen(IPv4(http_ip), http_port)

    do_action() = GitHub.run(listener, httpsock[], IPv4(http_ip), http_port; kwargs...)
    handle_exception(ex) = (isa(ex, Base.IOError) && (ex.code == -103)) ? :exit : :continue
    keep_running() = isopen(httpsock[])
    @info("GitHub webhook starting...", trigger, http_ip, http_port)
    recover("github_webhook", keep_running, do_action, handle_exception)
end

function request_processor(zsock::RequestSocket)
    do_action() = action(take!(event_queue), zsock)
    handle_exception(ex) = (isa(ex, InvalidStateException) && (ex.state == :closed)) ? :exit : :continue
    keep_running() = isopen(httpsock[])
    recover("request_processor", keep_running, do_action, handle_exception)
end

function main(config::AbstractString=isempty(ARGS) ? "config.toml" : first(ARGS); kwargs...)
    merge!(CONFIG, Pkg.TOML.parsefile(config)["commentbot"])
    if get(CONFIG, "enable_logging", true)
        global_logger(SimpleLogger(stdout, get_log_level(CONFIG["log_level"])))
    end
    zsock = RequestSocket(get(CONFIG, "backend_port", 5555))

    @info("Starting server...")
    t1 = @async request_processor(zsock)
    t2 = @async status_monitor(CONFIG["stop_file"], event_queue, httpsock)
    github_webhook(; kwargs...)
    wait(t1)
    wait(t2)

    # The !stopped! part is grep'ed by restart.sh, do not change
    @warn("Server !stopped!")
end

end    # module


text:
# advection.jl
import PDESolver.evalResidual

@doc """
### AdvectionEquationMod.evalResidual

This function evaluates the Advection equation.

**Inputs**

*  `mesh` : Abstract mesh object
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Advection equation object
*  `opts` : Options dictionary
*  `t`    :

Effectively updates eqn.res -- not eqn.res_vec. To make them consistent, use array3DTo1D on eqn.res and eqn.res_vec

**Outputs**

*  None

"""->
function evalResidual(mesh::AbstractMesh{Tmsh},
sbp::AbstractOperator, eqn::AdvectionData{Tsol, Tres, Tdim},
opts::Dict, t=0.0) where {Tmsh, Tsol, Tres, Tdim}

  myrank = mesh.myrank
  params = eqn.params
  #println("----- entered evalResidual -----")

  eqn.t = t
#  params.time.t_barriers[1] += @elapsed MPI.Barrier(mesh.comm)
  eqn.res = fill!(eqn.res, 0.0)  # Zero eqn.res for next function evaluation

  # start communication right away
  params.time.t_send += @elapsed if opts["parallel_type"] == 1
    setParallelData(eqn.shared_data, opts["parallel_data"])
    startSolutionExchange(mesh, sbp, eqn, opts)
    #  println("send parallel data @time printed above")
  end

  params.time.t_volume += @elapsed evalVolumeIntegrals(mesh, sbp, eqn, opts)
#  println("evalVolumeIntegrals @time printed above")

#  params.time.t_barriers[2] += @elapsed MPI.Barrier(mesh.comm)
  params.time.t_face += @elapsed if mesh.isDG
    evalFaceIntegrals(mesh, sbp, eqn, opts)
  end
#  println("evalFaceIntegrals @time printed above")

#  params.time.t_barriers[3] += @elapsed MPI.Barrier(mesh.comm)
  params.time.t_source += @elapsed evalSRCTerm(mesh, sbp, eqn, opts)
#  println("evalSRCTerm @time printed above")

#  params.time.t_barriers[4] += @elapsed MPI.Barrier(mesh.comm)
  params.time.t_bndry += @elapsed evalBoundaryIntegrals(mesh, sbp, eqn, opts)
#  println("evalBoundaryIntegrals @time printed above")

#  params.time.t_barriers[5] += @elapsed MPI.Barrier(mesh.comm)

  if opts["use_GLS2"]
    applyGLS2(mesh, sbp, eqn, opts, eqn.src_func)
  end
#  println("applyGLS2 @time printed above")

#  params.time.t_barriers[6] += @elapsed MPI.Barrier(mesh.comm)
  # do parallel computation last
  params.time.t_sharedface += @elapsed if mesh.commsize > 1
    evalSharedFaceIntegrals(mesh, sbp, eqn, opts)
  end

#  params.time.t_barriers[7] += @elapsed MPI.Barrier(mesh.comm)

  if opts["use_Minv"]
    applyMassMatrixInverse3D(mesh, sbp, eqn, opts, eqn.res)
  end

  @debug1 flush(params.f)
#  println(params.f, "----- finished evalResidual -----")
  return nothing
end

@doc """
### AdvectionEquationMod.evalVolumeIntegrals

  Evaluates the volume integrals of the weak form.  eqn.res is updated
  with the result.  Both the precompute_volume_flux and non
  precompute_volume_flux versions are contained within this function

  Inputs:

   `mesh` : mesh type
   `sbp`  : Summation-by-parts operator
   `eqn`  : Advection equation object
   `opts` : options dictionary

  Outputs

  None

"""->
function evalVolumeIntegrals(mesh::AbstractMesh{Tmsh},
                   sbp::AbstractOperator,
                   eqn::AdvectionData{Tsol, Tres, Tdim},
                   opts) where {Tmsh, Tsol, Tres, Tdim}

  if opts["use_staggered_grid"]
    calcVolumeIntegralsStaggered(mesh, mesh.mesh2, sbp, mesh.sbp2, eqn, opts)
  elseif opts["precompute_volume_flux"]
    calcAdvectionFlux(mesh, sbp, eqn, opts)
    for i=1:Tdim
      # multiplies flux_parametric by the SBP Q matrix (transposed), stores result in res
      # i is parametric direction
      weakdifferentiate!(sbp, i, sview(eqn.flux_parametric, :, :, :, i), eqn.res, trans=true)
    end
  else  # don't precompute flux

    alphas_xy = zeros(Float64, Tdim)  # advection coefficients in the xy
                                      #directions
    alphas_xy[1] = eqn.params.alpha_x
    alphas_xy[2] = eqn.params.alpha_y
    if Tdim == 3
      alphas_xy[3] = eqn.params.alpha_z
    end

    flux_el = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, Tdim)
    flux_tmp = zeros(Tres, Tdim)

    for i =1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        dxidx_j = ro_sview(mesh.dxidx, :, :, j, i)
        calcAdvectionFlux(eqn.params, eqn.q[1, j, i], alphas_xy, dxidx_j, flux_tmp)

        for k=1:Tdim
          flux_el[1, j, k] = flux_tmp[k]
        end

      end  # end loop j

      # do integration
      res_i = sview(eqn.res, :, :, i)
      for k=1:Tdim
        weakDifferentiateElement!(sbp, k, sview(flux_el, :, :, k), res_i,
                                  SummationByParts.Add(), true)
      end
    end  # end loop i

  end  # end if precompute_face_flux

  return nothing
end

"""
  Interpolate data for a single element from the solution grid to the flux
  grid.

  **Inputs**:

   * params
   * mesh: the solution mesh (needed for the interpolation operator)
   * q_s: the data on the solution mesh, can be either a vector or
          1 x mesh.numNodesPerelement matrix

  **Inputs/Outputs**:

   * q_f: the data on the flux mesh (overwritten)

  Note: q_s *cannot* be params.qL_s
"""
function interpolateElementStaggered(params::ParamType, mesh,
                                     q_s::AbstractArray, q_f::AbstractVector)

  qs_tmp = params.qL_s

  # extract data to vector
  for i=1:mesh.numNodesPerElement
    qs_tmp[i] = q_s[i]
  end

  smallmatvec!(mesh.I_S2F, qs_tmp, q_f)

  return nothing
end

"""
  Volume integrals for staggered grid.

  **Inputs**:

   * mesh_s: solution grid mesh
   * mesh_f: flux grid mesh
   * sbp_s: solution grid SBP
   * sbp_f: flux grid SBP
   * eqn: equation object (eqn.res is updated with the result)
   * opts: options dictionary

"""
function calcVolumeIntegralsStaggered(
             mesh_s::AbstractDGMesh{Tmsh},
             mesh_f::AbstractDGMesh{Tmsh},
             sbp_s::AbstractOperator,
             sbp_f::AbstractOperator,
             eqn::AdvectionData{Tsol, Tres, Tdim},
             opts) where {Tmsh, Tsol, Tres, Tdim}
  q_f = eqn.params.qL_f
  alphas_xy = zeros(Float64, Tdim)  # advection coefficients in the xy
                                    #directions
  alphas_xy[1] = eqn.params.alpha_x
  alphas_xy[2] = eqn.params.alpha_y
  if Tdim == 3
    alphas_xy[3] = eqn.params.alpha_z
  end
  flux_el = zeros(Tres, mesh_f.numNodesPerElement, Tdim)
  flux_tmp = zeros(Tres, Tdim)
  res_f = eqn.params.resL_f
  res_s = eqn.params.resL_s

  for i=1:mesh_f.numEl

    # interpolate to flux grid
    q_s = ro_sview(eqn.q, :, :, i)
    interpolateElementStaggered(eqn.params, mesh_s, q_s, q_f)

    # compute integral
    for j=1:mesh_f.numNodesPerElement
      dxidx_j = ro_sview(mesh_f.dxidx, :, :, j, i)
      calcAdvectionFlux(eqn.params, q_f[j], alphas_xy, dxidx_j, flux_tmp)

      for k=1:Tdim
        flux_el[j, k] = flux_tmp[k]
      end

    end  # end loop j

    # do integration
#    res_i = sview(eqn.res, :, :, i)

    fill!(res_f, 0.0)
    for k=1:Tdim
      weakDifferentiateElement!(sbp_f, k, sview(flux_el, :, k), res_f,
                                SummationByParts.Add(), true)
    end

    # interpolate residual back
    smallmatvec!(mesh_s.I_S2FT, res_f, res_s)

    # update residual
    for j=1:mesh_s.numNodesPerElement
      eqn.res[1, j, i] += res_s[j]
    end

  end  # end loop i

  return nothing
end

"""
  Populates eqn.flux_parametric.  Repeatedly calls the other method of this
  function.

  Inputs:

    mesh
    sbp
    eqn
    opts
"""
function calcAdvectionFlux(mesh::AbstractMesh{Tmsh}, sbp,
   eqn::AdvectionData{Tsol, Tres, Tdim}, opts) where {Tsol, Tres, Tdim, Tmsh}

  flux_parametric = eqn.flux_parametric

  alphas_xy = zeros(Float64, Tdim)      # advection coefficients in the xy directions
#  alphas_param = zeros(Tmsh, Tdim)      # advection coefficients in the parametric directions
  dxidx = mesh.dxidx                    # metric term
  q = eqn.q
  alphas_xy[1] = eqn.params.alpha_x
  alphas_xy[2] = eqn.params.alpha_y
  if Tdim == 3
    alphas_xy[3] = eqn.params.alpha_z
  end

  flux_tmp = zeros(Tres, Tdim)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      dxidx_j = ro_sview(mesh.dxidx, :, :, j, i)
      calcAdvectionFlux(eqn.params, q[1, j, i], alphas_xy, dxidx_j, flux_tmp)

      for k=1:Tdim
        flux_parametric[1, j, i, k] = flux_tmp[k]
      end
    end
  end

  return nothing
end

"""
  Calculates the advection flux in the parametric directions at a node.

  Inputs:

    params: a ParamType object
    q: the solution value at the node
    alphas_xy: the advection velocities in the x-y directions, vector of length
               Tdim
    dxidx: scaled mapping jacobian at the node, Tdim x Tdim matrix

  Inputs/Outputs:

    flux: vector of length Tdim to populate with the flux in the parametric
          directions
"""
function calcAdvectionFlux(
   params::ParamType{Tsol, Tres, Tdim}, q::Tsol,
   alphas_xy::AbstractVector,
   dxidx::AbstractMatrix{Tmsh},
   flux::AbstractVector{Tres}) where {Tsol, Tmsh, Tres, Tdim}

  for k=1:Tdim  # loop over parametric dimensions
    alpha_k = zero(Tmsh)
    for p=1:Tdim  # sum up alpha in the current parametric dimension
      alpha_k += dxidx[k, p]*alphas_xy[p]
    end
    # the first index is the 'equation number'; advection has only one
    # equation so it's always 1
    # for a vector PDE, there would be more
    flux[k] = alpha_k*q
  end

  return nothing
end

@doc """
Evaluate boundary integrals for advection equation, updating eqn.res with
the result.

**Inputs**

*  `mesh` : Abstract mesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Advection equation object
*  `opts` : options dictionary

**Outputs**

*  None

"""->
function evalBoundaryIntegrals(mesh::AbstractMesh{Tmsh},
                   sbp::AbstractOperator, eqn::AdvectionData{Tsol, Tres, Tdim}, opts) where {Tmsh, Tsol, Tres, Tdim}

#  println("----- Entered evalBoundaryIntegrals -----")

  if mesh.isDG && opts["precompute_q_bndry"]
    boundaryinterpolate!(mesh.sbpface, mesh.bndryfaces, eqn.q, eqn.q_bndry)
  end
#  println("    boundaryinterpolate @time printed above")

  for i=1:mesh.numBC
    functor_i = mesh.bndry_funcs[i]
    start_index = mesh.bndry_offsets[i]
    end_index = mesh.bndry_offsets[i+1]
    idx_range_i = start_index:(end_index-1)
    bndry_facenums_i = sview(mesh.bndryfaces, idx_range_i)
    bndryflux_i = sview(eqn.bndryflux, :, :, idx_range_i)

    # call the function that calculates the flux for this boundary condition
    # passing the functor into another function avoid type instability
    if opts["precompute_boundary_flux"]
      calcBoundaryFlux(mesh, sbp, eqn, functor_i, idx_range_i, bndry_facenums_i, bndryflux_i)
   else
     # this does the integration
     calcBoundaryFlux_nopre(mesh, sbp, eqn, functor_i, idx_range_i, bndry_facenums_i, bndryflux_i)
   end
#   println("    calcBoundaryflux @time printed above")
  end

  if opts["precompute_boundary_flux"]
    boundaryintegrate!(mesh.sbpface, mesh.bndryfaces, eqn.bndryflux, eqn.res)
  end

  return nothing
end # end function evalBoundaryIntegrals

@doc """
### AdvectionEquationMod.evalFaceIntegrals

  This function evaluates the interior face integrals for DG methods, using
  the flux function from eqn.flux_func.  The solution variables are interpolated
  to the faces, the flux computed, and then interpolated back to the
  solution points.

  This function also logs some quantities to disk (TODO: move this to
  Utils/logging)

  Inputs:
    mesh:  an AbstractDGMesh
    sbp
    eqn
    opts

"""->
function evalFaceIntegrals(mesh::AbstractDGMesh, sbp::AbstractOperator, eqn::AdvectionData,
                      opts)

#  println("----- Entered evalFaceIntegrals -----")
  if opts["use_staggered_grid"]
    calcFaceIntegralsStaggered_nopre(mesh, mesh.mesh2, sbp, mesh.sbp2, eqn, opts, eqn.flux_func)

  elseif opts["precompute_face_flux"]
    # interpolate solution to faces
    if opts["precompute_q_face"]
      interiorfaceinterpolate!(mesh.sbpface, mesh.interfaces, eqn.q, eqn.q_face)
    end
  #  println("    interiorface interpolate @time printed above")

    myrank = mesh.myrank
    if opts["writeqface"]
      writedlm("qface_$myrank.dat", eqn.q_face)
    end

    # calculate face fluxes
    calcFaceFlux(mesh, sbp, eqn, eqn.flux_func, mesh.interfaces, eqn.flux_face)
  #  println("    calcFaceFlux @time printed above")

    if opts["write_fluxface"]
      writedlm("fluxface_$myrank.dat", eqn.flux_face)
    end

    # integrate and interpolate back to solution points
    if mesh.isDG
      interiorfaceintegrate!(mesh.sbpface, mesh.interfaces, eqn.flux_face, eqn.res)
    else
      error("cannot evalFaceIntegrals for non DG mesh")
    end
  else
    if mesh.isDG
      calcFaceIntegrals_nopre(mesh, sbp, eqn, opts, eqn.flux_func)
    else
      error("cannot evalFaceIntegrals for non DG mesh")
    end
  end  # end if precompute_face_flux

  return nothing
end

@doc """
### AdvectionEquationMod.evalSRCTerm

  This function performs all the actions necessary to update eqn.res
  with the source term.  The source term is stored in eqn.src_func.  It is
  an abstract field, so it cannot be accessed (performantly) direction, so
  it is passed to an inner function.

  Inputs:
    mesh : Abstract mesh type
    sbp  : Summation-by-parts operator
    eqn  : Advection equation object
    opts : options dictonary

  Outputs: none

  Aliasing restrictions: none

"""->
function evalSRCTerm(mesh::AbstractMesh{Tmsh},
sbp::AbstractOperator, eqn::AdvectionData{Tsol, Tres, Tdim},
opts) where {Tmsh, Tsol, Tres, Tdim}

  # placeholder for multiple source term functionality (similar to how
  # boundary conditions are done)
  if opts["use_src_term"]
    applySRCTerm(mesh, sbp, eqn, opts, eqn.src_func)
  end

  return nothing
end  # end function

@doc """
### AdvectionEquationMod.applySRCTerm

  This function updates eqn.res with the source term.

  Inputs:
    mesh
    sbp
    eqn
    opts
    src_func:  the functor that returns the value of the source term at a node
               This functor must have the signature:
               src_func(coords, params, t)
               where coords is a vector of length 2 containing the x and y
               coordinates of the node, alpha_x and alpha_y are the advection
               coefficients, and t is the current time.

  Outputs: none

  Aliasing restrictions: none

"""->
function applySRCTerm(mesh,sbp, eqn, opts, src_func)

  weights = sbp.w
  alpha_x = eqn.params.alpha_x
  alpha_y = eqn.params.alpha_y

  t = eqn.t
  for i=1:mesh.numEl
    jac_i = ro_sview(mesh.jac, :, i)
    res_i = sview(eqn.res, :, :, i)
    for j=1:mesh.numNodesPerElement
      coords_j = ro_sview(mesh.coords, :, j, i)
      src_val = src_func(eqn.params, coords_j, t)
      res_i[j] += weights[j]*src_val/jac_i[j]
    end
  end

  return nothing
end

@doc """
### AdvectionEquationMod.evalSharedFaceIntegrals

  This function does the computation that needs the parallel
  communication to have finished already, namely the face integrals
  for the shared faces

  Inputs:
    mesh
    sbp
    eqn
    opts
"""->
function evalSharedFaceIntegrals(mesh::AbstractDGMesh, sbp, eqn, opts)

  if opts["face_integral_type"] != 1
    throw(ErrorException("unsupported face integral type"))
  end

  if opts["parallel_data"] == PARALLEL_DATA_FACE
#    println(eqn.params.f, "doing face integrals using face data")
    finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data,
                       calcSharedFaceIntegrals)
  elseif opts["parallel_data"] == PARALLEL_DATA_ELEMENT
#    println(eqn.params.f, "doing face integrals using element data")
    finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data,
                       calcSharedFaceIntegrals_element)
  else
    throw(ErrorException("unsupported parallel_data setting"))
  end

  return nothing
end

@doc """
### AdvectionEquationMod.init

  This function initializes the mesh and equation objects with any module
  specific data, such as boundary condition and source term functors.

  Inputs:
    mesh
    sbp
    eqn
    opts

  Outputs: none

  Aliasing restrictions: none
"""->
function init(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
eqn::AbstractAdvectionData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}

  println("Entering Advection Module")
  getBCFunctors(mesh, sbp, eqn, opts)
  getSRCFunctors(mesh, sbp, eqn, opts)
  if mesh.isDG
    getFluxFunctors(mesh, sbp, eqn, opts)
  end

  if opts["use_staggered_grid"]
    mesh2 = mesh.mesh2
    sbp2 = mesh.sbp2
    getBCFunctors(mesh2, sbp2, eqn, opts)
    getSRCFunctors(mesh2, sbp2, eqn, opts)
    if mesh.isDG
      getFluxFunctors(mesh2, sbp2, eqn, opts)
    end
  end

#  initMPIStructures(mesh, opts)
  return nothing
end

@doc """
### AdvectionEquationMod.majorIterationCallback

  This function gets called by the NonlinearSolvers during every major
  iteration.  This provides the opportunity to do output or monitoring.

  Inputs:
    itr: major iteration number
    mesh
    sbp
    eqn
    opts

    Outputs: none

    Aliasing restrictions: none

"""->
function majorIterationCallback(itr::Integer, mesh::AbstractMesh,
                                sbp::AbstractOperator, eqn::AbstractAdvectionData, opts, f::IO)
#=
  if mesh.myrank == 0
    println("Performing major Iteration Callback for iteration ", itr)
  end
=#
  if opts["write_vis"] && (((itr % opts["output_freq"])) == 0 || itr == 1)
    vals = real(eqn.q_vec)  # remove unneded imaginary part
    saveSolutionToMesh(mesh, vals)
#    cd("./SolutionFiles")
    fname = string("solution_", itr)
    writeVisFiles(mesh, fname)
    println(f, "finished writing vis file"); flush(f)

    #=
    # DEBUGGING: write error to file
    q_exact = zeros(eqn.q_vec)
    ex_func = ICDict[opts["IC_name"]]
    ex_func(mesh, sbp, eqn, opts, q_exact)
    q_err = real(eqn.q_vec) - q_exact
    saveSolutionToMesh(mesh, q_err)
    fname = string("error_", itr)
    writeVisFiles(mesh, fname)
    =#

#    cd("../")
  end

  return nothing
end

@doc """
### AdvectionEquationMod.assembleArray

  This function performs an assignment reduction of a 3D array to a vector.
  Note that because this is an assignment reduction, the order in which
  3D array is read matters, because only the last value assigned to a location
  in a vector remains.

  In most cases, what you really want is array3DTo1D().

  Inputs:
    mesh
    sbp
    eqn
    opts
    arr: the 3D array to be reduced into a vector

  Inputs/Outputs:
    res_vec: the vector to reduce the array into

  Keywords:
    zeros_resvec: whether or not to zero res_vec before performing the
                  reduction, default true

   Aliasing restrictions: none

"""->
#TODO: replace this with arrToVecAssign?
function assembleArray(mesh::AbstractMesh{Tmsh},
       sbp::AbstractOperator, eqn::AbstractAdvectionData{Tsol}, opts,
       arr::Abstract3DArray, res_vec::AbstractArray{Tres,1},
       zero_resvec=true) where {Tmsh, Tsol, Tres}
# arr is the array to be assembled into res_vec

#  println("in array3DTo1D")

  if zero_resvec
    fill!(res_vec, 0.0)
  end

  for i=1:mesh.numEl  # loop over elements
    for j=1:mesh.numNodesPerElement
      for k=size(arr, 1)  # loop over dofs on the node

        dofnum_k = mesh.dofs[k, j, i]

        res_vec[dofnum_k] = arr[k,j,i]
      end
    end
  end

  return nothing
end

@doc """
### AdvectionEquationMod.applyMassMatrixInverse3D

  This function applies the 3D inverse mass matrix to an array.
    The array passed in should always be eqn.res

  Inputs:
    mesh: mesh object, needed for numEl and numDofPerNode fields
    sbp: sbp object, needed for numnodes field
    eqn: equation object, needed for Minv3D field
    opts
    arr: the 3D array to have the 3D mass matrix inverse applied to it

"""->
function applyMassMatrixInverse3D(mesh, sbp, eqn, opts, arr)

  for i = 1:mesh.numEl
    for j = 1:sbp.numnodes
      for k = 1:mesh.numDofPerNode
        arr[k, j, i] = eqn.Minv3D[k, j, i] * arr[k, j, i]
      end
    end
  end

  return arr
end

# functions needed to make it compatible with the NonLinearSolvers module
function matVecA0inv(mesh::AbstractMesh{Tmsh},
sbp::AbstractOperator, eqn::AdvectionData{Tsol, Tres, Tdim},
opts, res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}

  return nothing
end

function matVecA0(mesh::AbstractMesh{Tmsh},
sbp::AbstractOperator, eqn::AdvectionData{Tsol, Tres, Tdim}, opts,
res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}

  return nothing
end

ext: jl


text:
# Path function wrapper that uses automatic differentiation to
# compute jacobian matricies
mutable struct ADPathFunction{type, PFT<:Function, SJC, CJC, STJC, TJC} <: PathFunction{type}
    # User defined function
    func!::PFT # SHould be of the form f!(out, x, u, p, t)

    # Number of functions
    nFuncs::Int

    # Number of states, controls, and static params
    nStates::Int
    nControls::Int
    nStatic::Int

    # Name of user defined functions
    funcNames::Vector{String}

    # Allocated output vector
    fOut::Vector{Float64}

    # Jacobian sparsity patterns
    stateSP::SparseMatrixCSC{Bool, Int}
    controlSP::SparseMatrixCSC{Bool, Int}
    staticSP::SparseMatrixCSC{Bool, Int}
    timeSP::SparseMatrixCSC{Bool, Int}

    # ForwardDiff.JacobianConfig objects
    stateJC::SJC
    controlJC::CJC
    staticJC::STJC
    timeJC::TJC
end

# To be fast, need to pass in number of states, controls, and static parameters
# to detect sparsity and allocate forward diff configuration objects.
function ADPathFunction(type::FunctionType, func!::Function, nFuncs::Int, nStates::Int, nControls::Int, nStatic::Int)
    # Temporarally allocate vectors for JacobianConfig creation
    out         = zeros(nFuncs)
    states      = rand(nStates)
    controls    = rand(nControls)
    static      = rand(nStatic)

    # Detect sparsity patterns and generate jacobian configuration objects
    if nStates > 0
        stateSP = jacobian_sparsity((y,x)->func!(y,x,controls,static,[0.0]),out,states)
        stateJC = ForwardDiff.JacobianConfig((y,x)->func!(y,x,controls,static,[0.0]),
            out, states, ForwardDiff.Chunk(states))

    else
        stateSP = sparse(Matrix{Bool}(undef, (0, 0)))
        stateJC = nothing
    end
    if nControls > 0
        controlSP = jacobian_sparsity((y,u)->func!(y,states,u,static,[0.0]),out,controls)
        controlJC = ForwardDiff.JacobianConfig((y,u)->func!(y,states,u,static,[0.0]),
            out, controls, ForwardDiff.Chunk(controls))
    else
        controlSP = sparse(Matrix{Bool}(undef, (0, 0)))
        controlJC = nothing
    end
    if nStatic > 0
        staticSP = jacobian_sparsity((y,p)->func!(y,states,controls,p,[0.0]),out,static)
        staticJC = ForwardDiff.JacobianConfig((y,p)->func!(y,states,controls,p,[0.0]),
            out, static, ForwardDiff.Chunk(static))
    else
        staticSP = sparse(Matrix{Bool}(undef, (0, 0)))
        staticJC = nothing
    end
    timeSP = jacobian_sparsity((y,t)->func!(y,states,controls,static,t[1]),out,rand(1))
    timeJC = ForwardDiff.JacobianConfig((y,t)->func!(y,states,controls,static,t[1]),
        out, [0.0], ForwardDiff.Chunk{1}())

    PFT     = typeof(func!)
    SJC     = typeof(stateJC)
    CJC     = typeof(controlJC)
    STJC    = typeof(staticJC)
    TJC     = typeof(timeJC)
    ADPathFunction{typeof(type),PFT,SJC,CJC,STJC,TJC}(func!,nFuncs,nStates,nControls,nStatic,
        Vector{String}(undef, 0),out,stateSP,controlSP,staticSP,timeSP,stateJC,controlJC,staticJC,timeJC)
end

# Evaluate jacobian functions
function EvaluateJacobian(jacType::State,
                          fp::ADPathFunction{type,PFT,SJC,CJC,STJC,TJC},
                          jac::AbstractMatrix,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat) where {type,PFT,SJC<:ForwardDiff.JacobianConfig,CJC,STJC,TJC}
    # Only evaluate if Jacobian has nonzero entries
    if nnz(fp.stateSP) > 0
        ForwardDiff.jacobian!(jac, (y,x)->fp.func!(y,x,control,static,time),
            fp.fOut, state, fp.stateJC, Val{false}())
    end
    return nothing
end

function EvaluateJacobian(jacType::State,
                          fp::ADPathFunction{type,PFT,SJC,CJC,STJC,TJC},
                          jac::AbstractMatrix,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat) where {type,PFT,SJC<:Nothing,CJC,STJC,TJC}
    return nothing
end

function EvaluateJacobian(jacType::Control,
                          fp::ADPathFunction{type,PFT,SJC,CJC,STJC,TJC},
                          jac::AbstractMatrix,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat) where {type,PFT,SJC,CJC<:ForwardDiff.JacobianConfig,STJC,TJC}
    if nnz(fp.controlSP) > 0
        ForwardDiff.jacobian!(jac, (y,u)->fp.func!(y,state,u,static,time),
            fp.fOut, control, fp.controlJC, Val{false}())
    end
    return nothing
end

function EvaluateJacobian(jacType::Control,
                          fp::ADPathFunction{type,PFT,SJC,CJC,STJC,TJC},
                          jac::AbstractMatrix,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat) where {type,PFT,SJC,CJC<:Nothing,STJC,TJC}
    return nothing
end

function EvaluateJacobian(jacType::Static,
                          fp::ADPathFunction{type,PFT,SJC,CJC,STJC,TJC},
                          jac::AbstractMatrix,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat) where {type,PFT,SJC,CJC,STJC<:ForwardDiff.JacobianConfig,TJC}
    if nnz(fp.stateSP) > 0
        ForwardDiff.jacobian!(jac, (y,p)->fp.func!(y,state,control,p,time),
            fp.fOut, static, fp.staticJC, Val{false}())
    end
    return nothing
end

function EvaluateJacobian(jacType::Static,
                          fp::ADPathFunction{type,PFT,SJC,CJC,STJC,TJC},
                          jac::AbstractMatrix,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat) where {type,PFT,SJC,CJC,STJC<:Nothing,TJC}
    return nothing
end

function EvaluateJacobian(jacType::Time,
                          fp::ADPathFunction{type,PFT,SJC,CJC,STJC,TJC},
                          jac::AbstractMatrix,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat) where {type,PFT,SJC,CJC,STJC,TJC}
    if nnz(fp.timeSP) > 0
        ForwardDiff.jacobian!(jac, (y,t)->fp.func!(y,state,control,static,t[1]),
            fp.fOut, @SVector([time]), fp.timeJC, Val{false}())
    end
    return nothing
end


text:
"""
    blosc_session(f)

Executes f between `blosc_init()` and `blosc_destroy()`. Sutable for blosc1 api calls
"""
function blosc_session(f)
    Lib.blosc_init()
    res = f()
    Lib.blosc_destroy()
    return res
end


text: print("Hello Julia!")

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
using Base: splat
using DataStructures: BinaryHeap

test = """
    #############
    #...........#
    ###B#C#B#D###
      #A#D#C#A#
      #########
    """

# solution for Q1 test
s1 = [(:C, :AB, 1), (:B, :C, 1), (:B, :BC, 1), (:AB, :B, 2), (:A, :B, 1),
      (:D, :CD, 1), (:D, :right, 1), (:CD, :D, 2), (:BC, :D, 1), (:right, :A, 1)]

# solution for Q2 test
s2 = [(:D, :right, 2), (:D, :left, 2), (:C, :right, 1), (:C, :CD, 1),
      (:C, :left, 1), (:B, :C, 3), (:B, :C, 2), (:B, :BC, 1), (:B, :AB, 1),
      (:BC, :B, 4), (:CD, :B, 3), (:right, :B, 2), (:D, :C, 1), (:D, :right, 1),
      (:AB, :D, 4), (:A, :B, 1), (:A, :D, 3), (:A, :AB, 1), (:left, :A, 3),
      (:left, :A, 2), (:AB, :D, 2), (:right, :A, 1), (:right, :D, 1)]

puzzle = """
    #############
    #...........#
    ###D#B#A#C###
      #B#D#A#C#
      #########
    """

extra_lines = """
          #D#C#B#A#
          #D#B#A#C#
        """

function parse_input(input::AbstractString; extra::Bool = false)
    lines = split(input, '\n'; keepempty = false)
    @assert lines[1] == "#############"
    @assert lines[2] == "#...........#"
    @assert lines[5] == "  #########"
    if extra
        lines = [lines[1:3]
                 split(extra_lines, '\n'; keepempty = false)
                 lines[4:5]]
    end
    lines = [split(line[3:11], "") for line in lines[3:end - 1]] |>
        splat(hcat) |> permutedims
    @assert all(==("#"), lines[1:end, 1:2:9])
    lines = lines[1:end, 2:2:8] .|> Symbol |> Matrix{Any}

    return (left  = Any[nothing, nothing],
            A     = lines[:, 1],
            AB    = Any[nothing],
            B     = lines[:, 2],
            BC    = Any[nothing],
            C     = lines[:, 3],
            CD    = Any[nothing],
            D     = lines[:, 4],
            right = Any[nothing, nothing])
end

const ROOMS = [:A, :B, :C, :D]
const HALLS = [:left, :AB, :BC, :CD, :right]
const COSTS = Dict(:A => 1, :B => 10, :C => 100, :D => 1000)
const PATHS = let paths = Dict(),
    places = [:left, :A, :AB, :B, :BC, :C, :CD, :D, :right]
    for from in places, to in places
        (a, b) = string.(places[clamp.(sort(indexin([from, to], places)), 2, 8)])
        path = []
        while a != b && last(a) != first(b)
            a = "$(last(a))$(last(a) + 1)"
            push!(path, Symbol(a))
        end
        paths[from => to] = path
    end
    paths
end

distance(from, to) = 2 * length(PATHS[from => to])
iscomfy(b, t) = all(v -> v === nothing || v == t, b[t])
Base.isempty(b, t) = all(isnothing, b[t])
Base.isempty(b, from, to) = all(h -> isempty(b, h), PATHS[from => to])
top(b, t) = b[t][findfirst(!isnothing, b[t])]

function moves(b)::Vector
    ms = []
    for room in ROOMS
        !iscomfy(b, room) || continue
        # from room to hallway
        for hall in HALLS
            isempty(b, room, hall) && append!(ms, moves(b, room, hall))
        end
        # from room to room - not strictly necessary, but this way we can
        # verify the solver with the official solutions
        t = top(b, room)
        iscomfy(b, t) && isempty(b, room, t) && push!(ms, move(b, room, t))
    end
    for hall in HALLS
        !isempty(b, hall) || continue
        # from hallway to room
        t = top(b, hall)
        iscomfy(b, t) && isempty(b, hall, t) && push!(ms, move(b, hall, t))
    end
    return ms
end

function moves(b, from, to)
    ms = []
    to ∈ [:left, :right] && b[to][1] == nothing &&
        push!(ms, move(b, from, to, 1))
    isempty(b, to) && push!(ms, move(b, from, to))
    return ms
end

function move(b, from, to, j = nothing)
    i = findfirst(!isnothing, b[from])
    @assert i !== nothing "Cannot move from $from to $to:\n$b"
    if j === nothing
        j = findfirst(!isnothing, [b[to]; 1]) - 1
    end
    t = b[from][i]
    ch = (; from => replace(b[from], i, nothing),
          to => replace(b[to], j, t))
    m = merge(b, ch)
    return ((distance(from, to) + i + j) * COSTS[t], m, (from, to, j))
end

function Base.replace(v::AbstractVector{<: T}, i::Integer, x::T) where T
    return [v[1:i - 1]; x; v[i + 1:end]]
end

function isfinal(b)
    return all(hall -> isempty(b, hall), HALLS) &&
        all(room -> iscomfy(b, room), ROOMS)
end

function search(b)
    h = BinaryHeap(Base.By(first), [(0, b, [])])
    closed = Set{UInt64}()
    while true
        (c, b, p) = pop!(h)
        hb = hash(b)
        hb ∈ closed && continue
        isfinal(b) && return (cost = c, path = p)
        push!(closed, hb)
        for (cm, bm, m) in moves(b)
            push!(h, (c + cm, bm, [p; m]))
        end
    end
end

function check(b, path, cost)
    for (step, (from, to, j)) in enumerate(path)
        m = move(b, from, to, j)
        @assert m ∈ moves(b) "#$step $b $from $to[$j]"
        (c, b, _) = m
        cost -= c
    end
    @assert cost == 0 "residual cost"
end

let b1 = parse_input(test),
    b2 = parse_input(test; extra = true)

    check(b1, s1, 12521)
    check(b2, s2, 44169)
    @assert search(b1).cost == 12521
    @assert search(b2).cost == 44169
end

b1 = parse_input(puzzle)
b2 = parse_input(puzzle; extra = true)
println("Q1: ", search(b1).cost)
println("Q2: ", search(b2).cost)

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
abstract type XMLError <: Exception end

struct XMLNoRootError <: XMLError ; end

struct XMLAttributeNotFound <: XMLError ; end

struct XMLParseError{T<:AbstractString} <: XMLError
    msg::T
end

struct XMLWriteError{T<:AbstractString} <: XMLError
    msg::T
end

struct XMLTreeError{T<:AbstractString} <: XMLError
    msg::T
end


text:
import Base.map

# Utility structure for collections of samples.
mutable struct Particles{C}
    calls::Vector{C}
    lws::Vector{Float64}
    lmle::Float64
end

map(fn::Function, ps::Particles) = map(fn, ps.calls)

include("inference/is.jl")
include("inference/pf.jl")
include("inference/mh.jl")
include("inference/vi.jl")
include("inference/hmc.jl")

const hmc = hamiltonian_monte_carlo
const mh = metropolis_hastings
const mh! = metropolis_hastings!
const is = importance_sampling
const advi = automatic_differentiation_variational_inference

# ------------ Documentation (IS) ------------ #

@doc(
"""
Samples from the model prior.
```julia
particles, normalized_weights = importance_sampling(observations::ConstrainedSelection,
                                                    num_samples::Int,
                                                    model::Function,
                                                    args::Tuple)
```
Samples from a programmer-provided proposal function.
```julia
particles, normalized_weights = importance_sampling(observations::ConstrainedSelection,
                                                    num_samples::Int,
                                                    model::Function,
                                                    args::Tuple,
                                                    proposal::Function,
                                                    proposal_args::Tuple)
```

Run importance sampling on the posterior over unconstrained addresses and values. Returns an instance of `Particles` and normalized weights.
""", importance_sampling)

# ------------ Documentation (PF) ------------ #

@doc(
"""
```julia
particles = initialize_filter(observations::ConstrainedHierarchicalSelection,
                              num_particles::Int,
                              fn::Function,
                              args::Tuple)
```

Instantiate a set of particles using a call to `importance_sampling`.
""", initialize_filter)

@doc(
"""
```julia
filter_step!(observations::ConstrainedHierarchicalSelection,
             ps::Particles,
             new_args::Tuple)
```
Perform a single filter step from an instance `ps` of `Particles`, applying the constraints specified by `observations`.

```julia
filter_step!(observations::ConstrainedHierarchicalSelection,
             ps::Particles,
             new_args::Tuple,
             proposal::Function,
             proposal_args::Tuple)
```
Perform a single filter step using a custom proposal function, applying the constraints specified by `observations`.
""", filter_step!)

@doc(
"""
```julia
check_ess_resample!(ps::Particles)
```
Checks the effective sample size using `ess`, then resamples from an existing instance of `Particles` by mutation in place.
""", check_ess_resample!)

# ------------ Documentation (MH) ------------ #

@doc(
"""
```julia
call, accepted, metropolis_hastings(sel::UnconstrainedSelection,
                                    call::HierarchicalCallSite)
```

Perform a Metropolis-Hastings step by proposing new choices using the prior at addressed specified by `sel`. Returns a call site, as well as a Boolean value `accepted` to indicate if the proposal was accepted or rejected.

```julia
call, accepted = metropolis_hastings(sel::UnconstrainedSelection,
                                     call::HierarchicalCallSite,
                                     proposal::Function,
                                     proposal_args::Tuple)
```

Perform a Metropolis-Hastings step by proposing new choices using a custom proposal at addressed specified by `sel`. Returns a call site, as well as a Boolean value `accepted` to indicate if the proposal was accepted or rejected.
""", metropolis_hastings)

# ------------ Documentation (VI) ------------ #

@doc(
"""
```julia
params, elbows, call_sites =  advi(sel::K,
                                   iters::Int,
                                   v_mod::Function,
                                   v_args::Tuple,
                                   mod::Function,
                                   args::Tuple;
                                   opt = ADAM(),
                                   gs_samples = 100) where K <: ConstrainedSelection
```

Given a selection `sel`, perform _automatic-differentiation variational inference_ with a proposal model `v_mod`. The result is a new set of trained parameters `params` for the variational model, the history of ELBO estimates `elbows`, and the call sites `calls` produced by the gradient estimator computation.
""", advi)


text:
#==========================================================
Code for comparing solvers, with time limits

This will run tests of solvers inside a thread so that it can impose a
time limit.

For this reason, julia should either be started with workers "julia -p 2",
or workers should be added before including Laplacians, like
addprocs()
@everywhere using Laplacians
@everywere include(this_file)

Routines for plotting results and producing reports are in
compare_solvers_reports.jl

The timeout features won't work with the matlab solvers.

===========================================================#

"""
    function speedTestLapTL{Tv,Ti}(solvers, dic, a::SparseMatrixCSC{Tv,Ti}, b::Array{Tv,1}; tol::Real=1e-2, maxits=Inf, maxtime=Inf, verbose=false)

Runs many Laplacians solvers.  Puts the build and solve time results into a dictionary dic.  It would be easiest to look at it via DataFrame(dic).  Returns the answer from the last solver.  `solvers` should be an array of `SolverTest`.

Uses threads to time limit all solvers after the first to at most 10x the time.
"""
function speedTestLapTL(solvers, dic, a::SparseMatrixCSC{Tv,Ti}, b::Array{Tv,1}; tol::Real=1e-2, maxits=10000, maxtime=1000, verbose=false, testName="") where {Tv,Ti}

    b = b - mean(b)

    la = lap(a)

    it = Int[1]

    # init the dict, if it is the first time
    initDictCol!(dic, "nv", Int)
    initDictCol!(dic, "ne", Int)
    initDictCol!(dic, "hash_a", UInt64)
    initDictCol!(dic, "testName", String)

    solvecol(name) = "$(name)_solve"
    buildcol(name) = "$(name)_build"
    totcol(name) = "$(name)_tot"
    itscol(name) = "$(name)_its"
    errcol(name) = "$(name)_err"

    dic["names"] = [t.name for t in solvers]

    for solver in solvers
        name = solver.name
        initDictCol!(dic, solvecol(name), Float64)
        initDictCol!(dic, buildcol(name), Float64)
        initDictCol!(dic, totcol(name), Float64)
        initDictCol!(dic, itscol(name), Float64)
        initDictCol!(dic, errcol(name), Float64)
    end

    nv = size(a,1)
    ne = nnz(a)
    hash_a = hash(a)

    push!(dic["nv"],nv)
    push!(dic["ne"],ne)
    push!(dic["hash_a"],hash_a)
    push!(dic["testName"],testName)

    x = []

    i = 1

    solverTest = solvers[i]

    if verbose
        println()
        println(solverTest.name)
    end

    t0 = time()
    ret = testSolver(solverTest.solver, a, b, tol, maxits, verbose)
    tl = 5 + 10*(time()-t0)
    if verbose
      println("time limit: ", tl)
    end

    pushSpeedResult!(dic, solverTest.name, ret)

    for i in 2:length(solvers)
        solverTest = solvers[i]

        if verbose
            println()
            println(solverTest.name)
        end

        #time_limit = 10.0

        ret = testSolverTimed(tl, solverTest.solver, a, b, tol, maxits, verbose)

        pushSpeedResult!(dic, solverTest.name, ret)

    end

    return x

end

function testSolverTimed(timelimit, solver, a, b, tol, maxits, verbose)

  if length(workers()) < 1
    error("Run addprocs() before setting up this code.")
  end
  proc = workers()[1]
  remotecall_fetch(gc, proc)

  c = Channel(1)
  @async put!(c, remotecall_fetch(()->(testSolver(solver, a, b, tol, maxits, verbose)),proc))
  t0 = time()
  while !isready(c) && (time() - t0) < timelimit
    sleep(0.1)
  end

  ret = (Inf,Inf,Inf,Inf)

  if isready(c)
    ret = fetch(c)

  else
    println("Interrupt process at time", time()-t0)
    interrupt(proc)
    sleep(1)
    ret = (Inf,Inf,Inf,Inf)

  end

  return ret

end

import Laplacians.initDictCol!
import Laplacians.testSolver
import Laplacians.pushSpeedResult!

"""
Runs many Laplacians solvers.  Puts the build and solve time results into a dictionary dic.  It would be easiest to look at it via DataFrame(dic).  Returns the answer from the last solver.  `solvers` should be an array of `SolverTest`.

Also compares them against the solvers we have in matlab, with a time limit of 10x the first solver here.
"""
function testVMatlabLap(solvers, dic::Dict, a::SparseMatrixCSC{Tv,Ti}, b::Array{Tv,1};
  tol::Real=1e-8, maxits=1000, maxtime=1000, verbose=false, testName="",
  test_icc=true, test_cmg=true, test_lamg=true, tl_fac=10) where {Tv,Ti}

    b = b - mean(b)

    la = lap(a)

    it = Int[1]

    # init the dict, if it is the first time
    initDictCol!(dic, "nv", Int)
    initDictCol!(dic, "ne", Int)
    initDictCol!(dic, "hash_a", UInt64)
    initDictCol!(dic, "testName", String)

    solvecol(name) = "$(name)_solve"
    buildcol(name) = "$(name)_build"
    totcol(name) = "$(name)_tot"
    itscol(name) = "$(name)_its"
    errcol(name) = "$(name)_err"

    dic["names"] = Array{String}(0)
    for t in solvers
        push!(dic["names"], t.name)
    end

    if test_cmg
      push!(dic["names"], "cmg")
    end
    if test_icc
      push!(dic["names"], "icc")
    end
    if test_lamg
      push!(dic["names"], "lamg")
    end

    for name in dic["names"]
        initDictCol!(dic, solvecol(name), Float64)
        initDictCol!(dic, buildcol(name), Float64)
        initDictCol!(dic, totcol(name), Float64)
        initDictCol!(dic, itscol(name), Float64)
        initDictCol!(dic, errcol(name), Float64)
    end

    nv = size(a,1)
    ne = nnz(a)
    hash_a = hash(a)

    push!(dic["nv"],nv)
    push!(dic["ne"],ne)
    push!(dic["hash_a"],hash_a)
    push!(dic["testName"],testName)

    x = []

    tl = 0

    for i in 1:length(solvers)
        solverTest = solvers[i]

        if verbose
            println()
            println(solverTest.name)
        end

        ret = testSolver(solverTest.solver, a, b, tol, maxits, verbose)

        if i == 1
            x = ret[5]
            tl = round(Int,30 + tl_fac*(ret[1]+ret[2]))
        end

        pushSpeedResult!(dic, solverTest.name, ret)

    end

    if tl == 0
        error("tl is zero")
    end

    if test_cmg
      if verbose
          println("cmg")
      end

      ret = timeLimitCmg(tl, lap(a), b,verbose = true);
      pushSpeedResult!(dic, "cmg", ret)
    end

    if test_icc
      if verbose
          println("icc")
      end
      ret = timeLimitIcc(tl, lap(a), b,verbose = true);
      pushSpeedResult!(dic, "icc", ret)
    end

    if test_lamg
      if verbose
          println("lamg")
      end
      ret = timeLimitLamg(tl, lap(a), b,verbose = true);
      pushSpeedResult!(dic, "lamg", ret)
    end

    return x

end

function testVMatlabLap(solvers::Array, dic::Dict, maker::Function; testName="")
    println(testName)
    a = maker()
    n = size(a,1)
    b = randn(n);
    b = b - mean(b);
    testVMatlabLap(solvers, dic, a, b, testName=testName, verbose=true)
end

function testVMatlabSddm{Tv,Ti}(solvers, dic::Dict, sdd::SparseMatrixCSC{Tv,Ti}, b::Array{Tv,1};
  tol::Real=1e-8, maxits=1000, maxtime=1000, verbose=false, testName="",
  test_icc=true, test_cmg=true, test_lamg=true, tl_fac=10)

    it = Int[1]

    # init the dict, if it is the first time
    initDictCol!(dic, "nv", Int)
    initDictCol!(dic, "ne", Int)
    initDictCol!(dic, "hash_a", UInt64)
    initDictCol!(dic, "testName", String)

    solvecol(name) = "$(name)_solve"
    buildcol(name) = "$(name)_build"
    totcol(name) = "$(name)_tot"
    itscol(name) = "$(name)_its"
    errcol(name) = "$(name)_err"

    dic["names"] = Array{String}(0)
    for t in solvers
        push!(dic["names"], t.name)
    end

    if test_cmg
      push!(dic["names"], "cmg")
    end
    if test_icc
      push!(dic["names"], "icc")
    end
    if test_lamg
      push!(dic["names"], "lamg")
    end

    for name in dic["names"]
        initDictCol!(dic, solvecol(name), Float64)
        initDictCol!(dic, buildcol(name), Float64)
        initDictCol!(dic, totcol(name), Float64)
        initDictCol!(dic, itscol(name), Float64)
        initDictCol!(dic, errcol(name), Float64)
    end

    nv = size(sdd,1)
    ne = nnz(sdd)
    hash_a = hash(sdd)

    push!(dic["nv"],nv)
    push!(dic["ne"],ne)
    push!(dic["hash_a"],hash_a)
    push!(dic["testName"],testName)

    x = []

    tl = 0

    for i in 1:length(solvers)
        solverTest = solvers[i]

        if verbose
            println()
            println(solverTest.name)
        end

        ret = testSolver(solverTest.solver, sdd, b, tol, maxits, verbose)

        if i == 1
            x = ret[5]
            tl = round(Int,30 + tl_fac*(ret[1]+ret[2]))
        end

        pushSpeedResult!(dic, solverTest.name, ret)

    end

    if tl == 0
        error("tl is zero")
    end

    if test_cmg
      if verbose
          println("cmg")
      end

      ret = timeLimitCmg(tl, sdd, b,verbose = true);
      pushSpeedResult!(dic, "cmg", ret)
    end

    if test_icc
      if verbose
          println("icc")
      end
      ret = timeLimitIcc(tl, sdd, b,verbose = true);
      pushSpeedResult!(dic, "icc", ret)
    end

    if test_lamg
      if verbose
          println("lamg")
      end
      ret = timeLimitLamgSddm(tl, sdd, b,verbose = true);
      pushSpeedResult!(dic, "lamg", ret)
    end

    return x

end


text:
module TestSymbolic
    using Calculus2
    using Base.Test

    #
    # Univariate Calculus
    #

    @test isequal(gradexpr(:(2), :x), 0)
    @test isequal(gradexpr(:(x), :x), 1)
    @test isequal(gradexpr(:(x + x), :x), 2)
    @test isequal(gradexpr(:(x - x), :x), 0)
    @test isequal(gradexpr(:(2 * x), :x), 2)
    @test isequal(gradexpr(:(2 / x), :x), :(-2 / x^2))
    @test isequal(gradexpr(:(x / 2), :x), 0.5)
    @test isequal(gradexpr(:(sin(x) / x), :x), :((cos(x) * x - sin(x)) / x^2))
    @test isequal(gradexpr(:(x * 2), :x), 2)
    @test isequal(gradexpr(:(a * x), :x), :a)
    @test isequal(gradexpr(:(x * a), :x), :a)
    @test isequal(gradexpr(:(x ^ 2), :x), :(2 * x))
    @test isequal(gradexpr(:(a * x ^ 2), :x), :(a * (2 * x)))
    @test isequal(gradexpr(:(2 ^ x), :x), :(*(0.6931471805599453, ^(2, x))))
    @test isequal(gradexpr(:(sin(x)), :x), :(cos(x)))
    @test isequal(gradexpr(:(cos(x)), :x), :(-sin(x)))
    @test isequal(gradexpr(:(tan(x)), :x), :(1 + tan(x)^2))
    @test isequal(gradexpr(:(exp(x)), :x), :(exp(x)))
    @test isequal(gradexpr(:(log(x)), :x), :(1 / x))
    @test isequal(gradexpr(:(sin(x) + sin(x)), :x), :(cos(x) + cos(x)))
    @test isequal(gradexpr(:(sin(x) - cos(x)), :x), :(cos(x) + sin(x)))
    @test isequal(gradexpr(:(x * sin(x)), :x), :(sin(x) + x * cos(x)))
    @test isequal(gradexpr(:(x / sin(x)), :x), :((sin(x) - x * cos(x)) / (sin(x)^2)))
    @test isequal(gradexpr(:(sin(sin(x))), :x), :(*(cos(x),cos(sin(x)))))
    @test isequal(gradexpr(:(sin(cos(x) + sin(x))), :x), :(*(+(-sin(x),cos(x)),cos(+(cos(x),sin(x))))))
    @test isequal(gradexpr(:(exp(-x)), :x), :(-exp(-x)))
    @test isequal(gradexpr(:(log(x^2)), :x), :(/(*(2,x),^(x,2))))
    @test isequal(gradexpr(:(x^n), :x), :(*(n, ^(x, -(n, 1)))))
    @test isequal(gradexpr(:(n^x), :x), :(*(^(n, x), log(n))))
    @test isequal(gradexpr(:(n^n), :x), 0)

    #
    # Multivariate Calculus
    #

    @test isequal(gradexpr(:(sin(x) + sin(y)), [:x, :y]), [:(cos(x)), :(cos(y))])
    @test isequal(gradexpr(:(x^2), [:x, :y]), {:(2*x), 0})

    # TODO: Get the generalized power rule right.
    # @test isequal(gradexpr(:(sin(x)^2), :x), :(2 * sin(x) * cos(x)))

    # TODO: Make these work
    # gradexpr(:(sin(x)), :x)(0.0)
    # gradexpr(:(sin(x)), :x)(1.0)
    # gradexpr(:(sin(x)), :x)(pi)

    #
    # SymbolicVariable use
    #

    x = BasicVariable(:x)
    y = BasicVariable(:y)

    @test isequal(@sexpr(x + y), :($x + $y))
    @test isequal(gradexpr(@sexpr(3 * x), x), 3)
    @test isequal(gradexpr(:(sin(sin(x))), :x), :(*(cos(x),cos(sin(x)))))
    @test isequal(gradexpr(@sexpr(sin(sin(x))), x), :(*(cos($x),cos(sin($x)))))

    function testfun(x)
        z = BasicVariable(:z)
        gradexpr(@sexpr(3*x + x^2*z), z)
    end

    @test isequal(testfun(x), :(^($(x),2)))
    @test isequal(testfun(3), 9)
    @test isequal(testfun(@sexpr(x+y)), :(^(+($x,$y),2)))

    #
    # Simplify tests
    #

    @test isequal(simplify(:(x+y)), :(+(x,y)))
    @test isequal(simplify(:(x+3)), :(+(3,x)))
    @test isequal(simplify(:(x+3+4)), :(+(7,x)))
    @test isequal(simplify(:(2+y+x+3)), :(+(5,y,x)))

    @test isequal(simplify(:(x*y)), :(*(x,y)))
    @test isequal(simplify(:(x*3)), :(*(3,x)))
    @test isequal(simplify(:(x*3*4)), :(*(12,x)))
    @test isequal(simplify(:(2*y*x*3)), :(*(6,y,x)))
end


text:
using Test
@time begin
    @time include("test_mesh.jl")
    @time include("test_quadrature.jl")
    @time include("test_basis.jl")
    @time include("test_fe.jl")
    @time include("test_utils.jl")
    @time include("test_ScalarFuncSp.jl")
    @time include("test_FunctionSpace.jl")
    @time include("test_handlers.jl")
    @time include("test_CGExample.jl")
end


text:
import ACE.OrthPolys: TransformedPolys

@doc raw"""
`struct Rn1pBasis <: OneParticleBasis`

One-particle basis of the form $R_n(r_{ij})$, i.e.,
no dependence on species or on $l$.

This does two things: (1) translates the `TransformedPolys` into a valid
one-particle basis; and (2) treat it as having vectorial input, i.e. value is
scalar but gradient is vectorial.

The default symbols are `:rr` for the state and `:n` for the index of the
basis function.
"""
mutable struct Rn1pBasis{T, TT, TJ, VSYM, NSYM} <: OneParticleBasis{T}
   R::TransformedPolys{T, TT, TJ}
end

Rn1pBasis(R::TransformedPolys{T, TT, TJ}; varsym = :rr, nsym = :n
            ) where {T, TT, TJ} =
      Rn1pBasis{T, TT, TJ, varsym, nsym}(R)

# ---------------------- Implementation of Rn1pBasis

Base.length(basis::Rn1pBasis) = length(basis.R)

_varsym(::Rn1pBasis{T, TT, TJ, VSYM, NSYM}) where {T, TT, TJ, VSYM, NSYM} = VSYM
_nsym(::Rn1pBasis{T, TT, TJ, VSYM, NSYM}) where {T, TT, TJ, VSYM, NSYM} = NSYM
_n(b, basis::Rn1pBasis) = getproperty(b, _nsym(basis))
_rr(X, Rn::Rn1pBasis) = getproperty(X, _varsym(Rn))

# -> TODO : figure out how to do this well!!!
# Base.rand(basis::Ylm1pBasis) =
#       AtomState(rand(basis.zlist.list), ACE.Random.rand_vec(basis.J))

function get_spec(basis::Rn1pBasis)
   N = _nsym(basis)
   return [  NamedTuple{(N,)}((n,)) for n = 1:length(basis) ]
end

==(P1::Rn1pBasis, P2::Rn1pBasis) =  ACE._allfieldsequal(P1, P2)

write_dict(basis::Rn1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_Rn1pBasis",
          "R" => write_dict(basis.R),
          "varsym" => _varsym(basis),
          "nsym" => _nsym(basis) )

read_dict(::Val{:ACE_Rn1pBasis}, D::Dict) =
            Rn1pBasis(read_dict(D["R"]),
                      varsym = D["varsym"], nsym = D["nsym"])

fltype(basis::Rn1pBasis{T}) where T = T

gradtype(::Rn1pBasis{T}) where {T} = SVector{3, T}

symbols(Rn::Rn1pBasis) = [ _nsym(Rn) ]

indexrange(Rn::Rn1pBasis) = NamedTuple{(_nsym(Rn),)}((1:length(Rn),))

isadmissible(b, basis::Rn1pBasis) = (1 <= _n(b, basis) <= length(basis))

degree(b, basis::Rn1pBasis) = _n(b, basis) - 1

get_index(basis::Rn1pBasis, b) = _n(b, basis)

rand_radial(basis::Rn1pBasis) = rand_radial(basis.R)

# ---------------------------  Evaluation code
#

alloc_B(basis::Rn1pBasis, args...) = alloc_B(basis.R)

alloc_dB(basis::Rn1pBasis, args...) =
      zeros( SVector{3, fltype(basis.R)}, length(basis) )

alloc_temp(basis::Rn1pBasis, args...) = alloc_temp(basis.R)

alloc_temp_d(basis::Rn1pBasis, args...) =
      (
      # alloc_temp_d(basis.R)...,
      dRdr = zeros(fltype(basis.R), length(basis.R)),
      )

evaluate!(B, tmp, basis::Rn1pBasis, X::AbstractState) =
      evaluate!(B, tmp, basis.R, norm(_rr(X, basis)))

function evaluate_d!(dB, tmpd, basis::Rn1pBasis, X::AbstractState)
   rr = _rr(X, basis)
   r = norm(rr)
   r̂ = rr / r
   evaluate_d!(tmpd.dRdr, tmpd, basis.R, r)
   dB[:] .= Ref(r̂) .* tmpd.dRdr
   return dB
end

function evaluate_ed!(B, dB, tmpd, basis::Rn1pBasis, X::AbstractState)
   rr = _rr(X, basis)
   r = norm(rr)
   r̂ = rr / r
   evaluate!(B, tmpd, basis.R, norm(rr))
   evaluate_d!(tmpd.dRdr, tmpd, basis.R, r)
   dB[:] .= Ref(r̂) .* tmpd.dRdr
   return nothing
end


text: ## no code, others need to add


text:
export color_palette

const color_palette = [colorant"firebrick";
                colorant"orangered2";
                    colorant"goldenrod1";
             colorant"skyblue1";
                 colorant"slateblue";
                colorant"maroon3";
             colorant"grey70";
                 colorant"dodgerblue4";
              colorant"seagreen4"]


text:
function gencode_dir()
    dir = joinpath(opengene_datadir(), "gencode")
    if !isdir(dir)
        mkpath(dir)
    end
    return dir
end
const gencode_releases = Dict(
    "grch37" => Dict(
        "localfile" => joinpath(gencode_dir(), "gencode.v19.annotation.gtf.gz"),
        "source" => "ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz"
    ),
    "hg19" => Dict(
        "localfile" => joinpath(gencode_dir(), "gencode.v19.annotation.gtf.gz"),
        "source" => "ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz"
    ),
    "grch38" => Dict(
        "localfile" => joinpath(gencode_dir(), "gencode.v24.annotation.gtf.gz"),
        "source" => "ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_24/gencode.v24.annotation.gtf.gz"
    )
)


text:
###############################################################################################################
########################################## Basic resolve functions ############################################
###############################################################################################################

"""
	resolve!(t::Tree, S::SplitList; conflict=:fail, usemask=false, tau=0.)

Add splits in `S` to `t` by introducing internal nodes. New nodes are assigned a time `tau` (`0` by default).
If `conflict != :ignore`, will fail if a split `s` in `S` is not compatible with `t`. Otherwise, silently skip the conflicting splits.
"""
function resolve!(
	t::Tree{T}, S::SplitList;
	conflict=:fail, usemask=false, tau=0., safe=false
) where T
	# Label for created nodes
	label_i = parse(Int64, TreeTools.create_label(t, "RESOLVED")[10:end])
	#
	tsplits = SplitList(t)
	for (i,s) in enumerate(S)
		if !safe && !in(s, tsplits; usemask)
			if iscompatible(s, tsplits; usemask)
				if usemask
					roots = TreeTools.blca([t.lleaves[x] for x in leaves(S,i)]...)
				else
					roots = TreeTools.blca([t.lleaves[x] for x in leaves(S,i)]...)
				end
				R = lca(roots)
				# Creating a new node with `roots` as children and `r` as ancestor.
				nr = TreeNode(T(), label="RESOLVED_$(label_i)")
				label_i += 1
				for r in roots
					prunenode!(r)
					graftnode!(nr,r)
				end
				graftnode!(R, nr, tau=tau)
				push!(tsplits.splits, s)
			elseif conflict != :ignore
				error("Tried to resolve tree with an incompatible split.")
			end
		end
	end
	node2tree!(t, t.root)
end

"""
	resolve(trees::Dict{T, <:Tree}, splits::Dict{T, <:SplitList}; kwargs...) where T

Resolve `trees[s]` with splits in `splits[s]` by calling `resolve!`. Return dictionary of resolved trees. `trees` and `splits` must share keys. This is meant to be used for dictionaries of trees/splits indexed by flu segments.
"""
function resolve(trees::Dict{T, <:Tree}, splits::Dict{T, <:SplitList}; kwargs...) where T
	resolved_trees = Dict(k => copy(t) for (k,t) in trees)
	for (s,S) in splits
		resolve!(resolved_trees[s], S; kwargs...)
	end
	return resolved_trees
end

"""
	resolve!(S1::SplitList, t1::Tree, S2::SplitList)

Add splits of `S2` in `S1` if they resolve `t1`.
"""
function resolve!(S1new, S1::SplitList, t1::Tree, S2::SplitList)
	c = 0
	for s2 in S2
		r1 = t1.lleaves[S2.leaves[s2.dat[1]]]
		for i in 2:length(s2.dat)
			r1 = lca(r1, t1.lleaves[S2.leaves[s2.dat[i]]])
		end
		#r1 = lca(t1, S2.leaves[s2.dat]) # Ancestor of nodes in s2 in t1
		s1 = S1.splitmap[r1.label]
		if s1 != s2 && !in(s2, S1new) && arecompatible(s1, s2)
			# Consider the set of splits just below r1 that are subsplits of s2
			# If I join those, I should get exactly s2
			# Otherwise, can't use s2 to resolve r1
			stmp = Split(0)
			for n in r1.child
				if n.isleaf
					i = findfirst(==(n.label), S2.leaves)
					if in(i, s2.dat)
						TreeTools.joinsplits!(stmp, Split([i]))
					end
				else
					if TreeTools.is_sub_split(S1.splitmap[n.label], s2)
						TreeTools.joinsplits!(stmp, S1.splitmap[n.label])
					end
				end
			end

			if stmp == s2
				push!(S1.splits, s2)
				push!(S1new.splits, s2)
				c += 1
			end
		end
	end

	return c
end

function resolve!(S1::SplitList, S2::SplitList, t1::Tree, t2::Tree)
	nit = 0
	nitmax = 20
	flag = true
	S1new = SplitList(S1.leaves)
	S2new = SplitList(S2.leaves)
	while flag && nit < nitmax
		flag = false
		c = resolve!(S1new, S1, t1, S2)
		c != 0 && (flag = true)
		c = resolve!(S2new, S2, t2, S1)
		c != 0 && (flag = true)
		nit += 1
	end

	if nit == nitmax
		@warn "Maximum number of iterations reached"
	end

	return [S1new, S2new]
end

###############################################################################################################
################################### Resolve trees using eachother splits ######################################
###############################################################################################################

"""
	resolve!(t1::Tree, t2::Tree; tau=0.)

Resolve `t1` using splits of `t2` and inversely. Every split of `t2` a tree that is compatible with `t1` is introduced in `t1` with branch length `tau` (and inversely). Return new splits in each tree.
"""
function resolve!(t1::Tree, t2::Tree; tau=0.)
	S = [SplitList(t) for t in (t1,t2)]
	Snew = resolve!(S[1], S[2], t1, t2)
	for (t, s) in zip((t1,t2), S)
		resolve!(t, s, conflict=:fail, usemask=false, tau=tau)
	end

	return Snew
end

###############################################################################################################
############################## Resolve trees using inferred compatible clades #################################
###############################################################################################################

"""
	resolve!(t1::Tree, t2::Tree, MCCs; tau = 0.)

Resolve `t1` using `t2` and inversely using the list of MCCs.
New branches have a length `tau`.
Return the list of resolved splits in each tree.
"""
function resolve!(t1::Tree, t2::Tree, MCCs; tau = 0.)
	resolvable_splits = TreeKnit.new_splits(MCCs, t1, t2)
	resolve!(t1, resolvable_splits[1]; conflict=:fail, tau)
	resolve!(t2, resolvable_splits[2]; conflict=:fail, tau)

	return resolvable_splits
end

"""
	resolve_from_mccs!(infer_mccs::Function, trees::Vararg{Tree}; verbose=false, kwargs...)
	resolve_from_mccs!(infer_mccs::Function, trees::Dict; verbose=false, kwargs...)

Resolve `trees` using pairwise MCC **inference**.
The `infer_mccs` function must take a `Dict{<:Any, Tree}` as input.
`kwargs...` are fed to `infer_mccs`.
Return the set of compatible splits introduced.

This is done in four steps:
1. Compute MCCs for pairs of trees while resolving them
2. Find all new splits in each tree that are introduced by the found MCCs
3. Select only the splits compatible with all others
4. Introduce those in trees

These four steps are **iterated**, while new compatible splits are found.
When no new compatible split is found, compute final MCCs without resolving.
"""
function resolve_from_mccs!(infer_mccs::Function, trees::Dict{<:Any, <:Tree}; verbose=false)
	maxit = 10
	#
	verbose && println("It. 1/$(maxit)")
	newsplits = _resolve_from_mccs!(infer_mccs, trees, verbose=verbose)
	#
	nS = newsplits
	it = 1
	while !mapreduce(isempty, *, values(nS), init=true) && it < maxit
		verbose && println("It. $(it+1)/$(maxit)")
		nS = _resolve_from_mccs!(infer_mccs, trees, verbose=verbose)
		for (s,S) in nS
			append!(newsplits[s].splits, S.splits)
		end
		it += 1
	end
	return newsplits
end

function _resolve_from_mccs!(infer_mccs::Function, trees::Dict{<:Any, <:Tree}; verbose=false)
	verbose && println("--- Resolving with MCCs... ---\n")
	# 1. `MCCs[u,v]` corresponds to `trees[u]` and `trees[v]`
	MCCs = _computeMCCs(infer_mccs, trees)
	# 2. `resolvable_splits[s]` is an array of `SplitList` objects that can be introduced in `trees[s]` from each other tree.
	resolvable_splits = TreeKnit.new_splits(trees, MCCs)
	if verbose
		for (s,t) in trees
			Y = union(resolvable_splits[s]...)
			println("Tree $s: $(length(Y)) resolvable splits (potentially incompatible).")
		end
		println()
	end
	# 3. `cmpt_splits[s]` is a `SplitList` object containing a selection of splits
	# we can choose between several methods for this choice.
	# cmpt_splits = TreeKnit.max_clique_splits(resolvable_splits; verbose)
	cmpt_splits = TreeKnit.compatible_splits(resolvable_splits)
	if verbose
		for (s,t) in trees
			println("Tree $s: $(length(cmpt_splits[s])) compatible resolvable splits.")
		end
		println()
	end
	# 4.
	for (s,S) in cmpt_splits
		resolve!(trees[s], S, conflict=:fail)
		TreeTools.check_tree(trees[s])
	end
	#
	return cmpt_splits
end

function resolve_from_mccs!(
	infer_mccs::Function, trees::Vararg{Tree};
	verbose=false, kwargs...
)
	resolve_from_mccs!(
		infer_mccs, Dict(i=>t for (i,t) in enumerate(trees));
		verbose=verbose, kwargs...
	)
end

"""
    compatible_splits(rS::Dict{<:Any, Array{SplitList,1}})

Find splits in `rS` that are compatible with all other splits. Return a dictionary with the same keys as `rS`.
"""
function compatible_splits(rS::Dict{<:Any, Array{SplitList,1}})
    cs = Dict{Any, SplitList}() # Compatible splits (output)
    for (seg, aS) in rS # aS: new splits for segment seg from all other trees
        for (i,S) in enumerate(aS) # S: new splits in seg from tree i
            if !haskey(cs, seg)
                cs[seg] = SplitList(S.leaves)
            end
            for s in S
                # Check if s is compatible with aS[j] for all `j!=i`
                cmpt = true
                for j in Iterators.filter(!=(i), 1:length(aS))
                    if !iscompatible(s, aS[j], usemask=false)
                        cmpt = false
                        break
                    end
                end
                cmpt && push!(cs[seg].splits, s)
            end
        end
        cs[seg] = unique(cs[seg])
    end
    return cs
end

# function max_clique_splits(nS::Dict; verbose=false)
# 	cS = Dict{Any, SplitList}()
# 	for (seg, S) in nS
# 		verbose && println(seg)
# 		cS[seg] = max_clique_splits(S; verbose)
# 	end
# 	return cS
# end
# function max_clique_splits(nS; verbose=false)
# 	S = union(nS...)
# 	verbose && println("Finding max clique among $(length(S)) splits.")
# 	g = build_compat_graph(S)
# 	all_cliques = LightGraphs.maximal_cliques(g)
# 	verbose && println("Found $(length(all_cliques)) cliques.")
# 	if isempty(all_cliques)
# 		return S
# 	else
# 		max_clique = sort(all_cliques, by=x->length(x), rev=true)[1]
# 		# Delete splits not in the max clique
# 		todel = findall(!in(max_clique), 1:length(S))
# 		deleteat!(S.splits, todel)
# 		return S
# 	end
# end

# function build_compat_graph(S::SplitList)
# 	L = length(S)
# 	g = Graph(L)
# 	for i in 1:L, j in (i+1):L
# 		if arecompatible(S[i], S[j])
# 			add_edge!(g, i, j)
# 		end
# 	end
# 	return g
# end

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
#=
using Statistics: std #We have some problem in the calculation of the std

function BOSSfit(X,  y,
    threshold::Float64=0.92,
    n_parameter_samples::Int64=250,
    max_ensemble_size::Int64=500,
    min_window::Int64=10,
    word_lengths::Array=[16, 14, 12, 10, 8],
    weights,
    max_win_len_prop::Int64= 1,
    norm_options::Array=[true, false],
    #norm = false    use normalise instead
    alphabet_size::Int64=4,
    save_words::Bool=true, # for BOSSIndividual and false for SAF
    random_state=nothing,
    levels::Int64=1,
    normalise_dft::Bool=true,
    #bigrams = false   As it is always false in BOSS
    remove_repeat_words::Bool=false,
    level_weights::Array=[1, 2, 4, 16, 32, 64, 128]
               )
    n_classes = unique(y) #Take it directly from the interface,
    n_instances, series_length = size(X)
    max_window_searches = series_length ÷ 4
    max_window = series_length * max_win_len_prop
    win_inc = (max_window - min_window) ÷ max_window_searches
    word_length = word_lengths[1] # Try to use one word for it
    # Try making a struct with following
    classifiers = []
    classifiers_word_len = []
    classifiers_accuracy = []

    max_acc = -1
    min_max_acc = -1

    for i = 1:length(norm_options)
        for win_size = min_window:win_inc:(max_window + 1)
            # window_size = win_size   norm = normalise  Try to replace it with one word
            inverse_sqrt_win_size = 1 / sqrt(win_size)
            breakpoints = _mcb(X, win_size, alphabet_size, word_length, norm_options[i], inverse_sqrt_win_size, normalise_dft)
            bags, words_lists = transform(X, word_length, win_size, norm_options[i], inverse_sqrt_win_size,
                                                                 series_length, breakpoints, save_words, normalise_dft)

            best_classifier_for_win_size = temp = bags
            best_acc_for_win_size = -1
            best_word_len = word_length
            for j=1:length(word_lengths)
                if j > 1
                    bags = _shorten_bags(word_lengths[j], words_lists)  # Fix function
                end
                # Make Array or somthing for accuracy/ may be change the variable
                accuracy = _individual_train_acc(bags, y, train_size, best_acc_for_win_size)
                if accuracy >= best_acc_for_win_size
                    best_acc_for_win_size = accuracy
                    best_classifier_for_win_size = bags
                    best_word_len = word_lengths[j]
                end
            end
            if some_logic
                #def _clean(self):
                #    self.transformer.words = None
                #    self.transformer.save_words = False
                #def _set_word_len(self, word_len):
                #    self.word_length = word_len
                #    self.transformer.word_length = word_len
                push!(classifiers, best_classifier_for_win_size)
                push!(classifiers_word_len, best_word_len)
                push!(classifiers_accuracy, best_acc_for_win_size)
                if best_acc_for_win_size > max_acc

                end
                if length(classifiers) > max_ensemble_size

                end
            end
            # saving the class if acc > 92 of the best one _include_in_ensemble
        end
    end
end

function _discrete_fourier_transform(series, word_length, norm, inverse_sqrt_win_size, normalise_dft)
    _length = length(series)
    output_length = ceil(Int, word_length/2)
    start = norm ? 1 : 0

    _std = 1

    if normalise_dft
        s = std(series)  # some issue with std statement ??
        if s!=0
            _std = s
        end
    end

    dft = [ sum( [ [  series[n] * cos(2 * π * (n-1) * i / _length),
                     -series[n] * sin(2 * π * (n-1) * i / _length) ]
                                                   for n=1:_length ] ) # sum end
                                                   for i = start:(start + output_length-1) ] #end of dft

    dft = vcat(dft...)

    if normalise_dft
        dft = dft * (inverse_sqrt_win_size/ _std)
    end
    return dft
end

# self._create_word(), self._add_to_pyramid(),  self.levels,  self._add_to_bag(), self.bigrams (are not used ),
# self.save_words, self.words.append, self.window_size, self.word_length,

function transform(X, word_length, window_size, norm, inverse_sqrt_win_size, series_length, breakpoints, save_words, normalise_dft)    # What is diff b\w transform and _shorten_bags
    # check_is_fitted()
    # X = check_X(X, enforce_univariate=True)
    # X = tabularize(X, return_array=True)
    # bags = pd.DataFrame()
    dim = []
    levels = 1
    norm = norm # norm used in _discrete_fourier_transform and _mft while calcuation
    remove_repeat_words = false
    words_list = zeros(Float64, size(X)[1], size(dfts)[1])
    for i = 1:size(X)[1]
        dfts = _mft(X[i, :], word_length, window_size, norm, inverse_sqrt_win_size, series_length, normalise_dft)
        bag = Dict() # {} This is how they use the ditionary in the python
        last_word = -1
        repeat_words = 0
        words = zeros(Float64, size(dfts)[1])
        for window = 1:size(dfts)[1]
            # word = _create_word(dfts[window, :], breakpoints, alphabet_size, word_length)
            @inbounds for i = 1:word_length
                @inbounds for bp = 1:alphabet_size
                    if dfts[window, :] <= breakpoints[i][bp]
                        words[window] = (words[window] << 2) | bp
                        break
                    end
                end
            end
            # append!(words, word)
            # levels
            repeat_word = levels > 1 ? (_add_to_pyramid(bag, words[window], last_word, (window - (repeat_words ÷2)),
                          remove_repeat_words, series_length, level_weights, levels) : _add_to_bag(bag, words[window], last_word, remove_repeat_words))

            if repeat_word               # chek if order of end or if else got messed up
                repeat_words += 1
            else
                last_word = words[window]
                repeat_words = 0
            end
        end
       # bigrams removed
        if save_words
            words_list[i, :] =  words    #  self.words is equl to words_list
        end
        append!(dim, bag)
    end
    return dim, words_list
end

function _shorten_bags(word_len, words_lists)      # Check the logic changes w.r.t transform function
    dim = []
    w_instance, w_length = size(words_lists)
    for i = 1:w_instance
        bag = Dict()
        last_word = -1
        repeat_words = 0
        new_words = zeros(Int64, w_length)
        for j = 1:w_length
            new_word_len = shorten(16 - word_len)
            new_words[j] = (words_lists[i, j] >> (2*new_word_len))    # This might be wrong check corner cases
            repeat_word = levels > 1 ? (_add_to_pyramid(bag, new_words[j], last_word, (window - (repeat_words ÷2)),
                          remove_repeat_words, series_length, level_weights, levels) : _add_to_bag(bag, new_words[j], last_word, remove_repeat_words))
            if repeat_word
                repeat_words += 1
            else
                last_word = new_words[j]
                repeat_words = 0
            end
        end
        append!(dim, bag)
    end
    return dim
end

function _mft(series, word_length, window_size, norm, inverse_sqrt_win_size, series_length, normalise_dft)  # Relation b\w transform and mft while runing main loop
    start_offset = norm ? 2 : 0
    lengthl = (word_length + word_length) ÷ 2

    phis = [  [cos(2 * π * (-((i * 2) + start_offset) / 2) / window_size),
               -sin(2 * π * (-((i * 2) + start_offset) / 2) / window_size) ]
                                          for i = 0:((lengthl÷2)-1) ] #end of dft
    phis = vcat(phis...)
    endl = max(1, series_length - window_size + 1)
    stds = _calc_incremental_mean_std(series, endl, window_size)
    transformed = zeros(Float64, endl, lengthl)
    mft_data = []  # Find way to do catching
    @inbounds for i = 0:endl-1
        if i > 0
            @inbounds for n = 1:2:lengthl
                real = mft_data[n] + series[i + window_size] - series[i]
                imag = mft_data[n + 1]
                mft_data[n] = real * phis[n] - imag * phis[n + 1]
                mft_data[n + 1] = real * phis[n + 1] + phis[n] * imag
            end
        else
            mft_data = _discrete_fourier_transform(series[1:window_size], word_length,
                                               norm, inverse_sqrt_win_size, normalise_dft) #normalise =false
        end
        normalising_factor = ( (1 / ( stds[i + 1] > 0 ? stds[i + 1] : 1) ) * inverse_sqrt_win_size)
        transformed[i + 1, :] = mft_data * normalising_factor
    end
    return transformed
end

function _calc_incremental_mean_std(series, endl, window_size)
    means = zeros(endl)
    stds = zeros(endl)
    window = series[1:window_size]
    series_sum = sum(window)
    square_sum = sum(window.*window)

    r_window_length = 1 / window_size
    means[1] = series_sum * r_window_length
    buf = square_sum * r_window_length - means[1]^2
    stds[1] = buf > 0 ? sqrt(buf) : 0

    @inbounds for w =1:(endl-1)
        series_sum += series[w + window_size] - series[w]
        means[w + 1] = series_sum * r_window_length
        square_sum += series[w + window_size]^2 - series[w]^2
        buf = square_sum * r_window_length - means[w + 1]^2
        stds[w + 1] = buf > 0 ? sqrt(buf) : 0
    end
    return stds
end
# self.series_length  self.window_size, self.n_instances, self.alphabet_size,  self.word_length

function _mcb(X, window_size, alphabet_size, word_length, norm, inverse_sqrt_win_size, normalise_dft)
    n_instances, series_length = size(X)
    num_windows_per_inst = ceil(series_length / window_size)
    dft = zeros(Float64, n_instances, num_windows_per_inst, word_length)
    @inbounds for k=1:n_instances
        @inbounds for i=1:num_windows_per_inst-1
            dft[k, i, :] = _discrete_fourier_transform(view(X, k, ((i-1)*window_size + 1):(i*window_size)),
                                                word_length, norm, inverse_sqrt_win_size, normalise_dft)
        end
        dft[k, end, :] = _discrete_fourier_transform(view(X, k, (series_length-window_size + 1):series_length),
                                                word_length, norm, inverse_sqrt_win_size, normalise_dft)
    end
    total_num_windows = n_instances * num_windows_per_inst
    breakpoints = zeros(Float64, word_length, alphabet_size)
    @inbounds for letter in range(word_length)                        # sortring algo, making a array and finding the aproproate ramge to place
        sort!(dft[((letter-1)*total_num_windows + 1):(letter*total_num_windows)])
        bin_index = 0
        target_bin_depth = total_num_windows / alphabet_size
        @inbounds for bp in range(alphabet_size - 1)
            bin_index += target_bin_depth
            breakpoints[letter][bp] = dft[int(bin_index + (letter-1)*total_num_windows)]  # Check indexing here
        end
        breakpoints[letter][alphabet_size] = maxintfloat(Float64)   # check this out in the system
    return breakpoints
end

# self.window_size, self.window_size
#=
function _mcb_dft(series, num_windows_per_inst, window_size, series_length,
                        word_length, norm, inverse_sqrt_win_size, normalise)
    # Splits individual time series into windows and returns the DFT for
    # each
    Z = zeros(Float64, num_windows_per_inst, window_size)
    @inbounds for i=1:num_windows_per_inst-1
        Z[i, :] = _discrete_fourier_transform(series[((i-1)*window_size + 1):(i*window_size)],
                                            word_length, norm, inverse_sqrt_win_size, normalise)
    end
    Z[end, :] = _discrete_fourier_transform(series[ (series_length-window_size + 1):series_length],
                                            word_length, norm, inverse_sqrt_win_size, normalise)
    return Z # find the way to have two indiaces for he same loops @inoans may help
end
=#
# self.word_length, self.alphabet_size, self.breakpoints

#=
function _create_word(dft, word)                      # check out this function
    word = 0.0
    @inbounds for i = 1:word_length
        @inbounds for bp = 1:alphabet_size
            if dft[i] <= breakpoints[i][bp]
                word = (word << 2) | letter
                break
            end
        end
    end
    return word
end
=#
# self.remove_repeat_words

function _add_to_bag(bag, word, last_word, remove_repeat_words)
    if remove_repeat_words && word == last_word         #word.word
        return false
    end
    bag[word] = get(bag, word, 0) + 1   # check out this get fucntion
    return true
end
# self.levels, self.remove_repeat_words, self.window_size,  self.series_length, self.level_weights

function _add_to_pyramid(bag, word, last_word, window_ind, remove_repeat_words, series_length, level_weights, levels)
    if remove_repeat_words && word == last_word  # check if python `and` and `&&` from the julia are same
        return false
    end
    start = 0
    for i = 0:(levels-1)
        num_quadrants = 2^i                # Check out this function pow
        quadrant_size = series_length / num_quadrants
        pos = window_ind + (window_size ÷ 2)
        quadrant = start + (pos / quadrant_size)
        bag[(word, quadrant)] = (get(bag, (word, quadrant), 0) + level_weights[i+1])
        start += num_quadrants
    end
    return true
end

function _train_predict(train_num, bags)
    best_dist = maxintfloat(Float64)
    nn = nothing                                             # Check on what type nn should be?
    for n =  1:length(bags)
        if n == train_num
            continue
        end
        dist = boss_distance(bags[train_num], bags[n], best_dist)
        if dist < best_dist
            best_dist = dist
            nn = class_vals[n]                        # Make class_vals to loop over BOSSIndividual
        end
    end
    return nn
end

function _individual_train_acc(bags, y, train_size, lowest_acc)
    correct = 0
    required_correct = floor(Int64, lowest_acc * train_size)
    for i = 1:train_size
        if correct + train_size - i < required_correct     # Check the boundary case i or i-1 ?
            return -1
        end
        c = _train_predict(i, bags)
        if c == y[i]
            correct += 1
        end
    end
    return correct / train_size
end

function boss_distance(first, second, best_dist=maxintfloat(Float64))
    dist = 0
    if typeof(first) <: Dict
        for (word, val_a) in first
            val_b = second[word]         # check the cal again
            dist += (val_a - val_b) * (val_a - val_b)
            if dist >= best_dist
                return maxintfloat(Float64)
            end
        end
    else
        dist = sum( [ first[n] == 0 ? 0 : ( (first[n] - second[n]) * (first[n] - second[n]) ) # Check this code again
                                         for n=1:length(first) ] )
    end
    return dist
end

"""
julia> linespace(2,4,5, Array{Int})
5-element Array{Int64,1}:
 2
 2
 3
 3
 4
"""
linespace(start,stop,_length, _type) =  convert(_type, floor.(LinRange(start,stop,_length)))
=#


text:
module TestUtils
    using Compat, Compat.Test, DataFrames, StatsBase, Compat.Random
    using Compat: @warn
    import DataFrames: identifier

    @test identifier("%_B*_\tC*") == :_B_C_
    @test identifier("2a") == :x2a
    @test identifier("!") == :x!
    @test identifier("\t_*") == :_
    @test identifier("begin") == :_begin
    @test identifier("end") == :_end

    @test DataFrames.make_unique([:x, :x, :x_1, :x2], makeunique=true) == [:x, :x_2, :x_1, :x2]
    # TODO uncomment this line after deprecation period when makeunique=false throws error
    # @test_throws ArgumentError DataFrames.make_unique([:x, :x, :x_1, :x2], makeunique=false)
    @test DataFrames.make_unique([:x, :x_1, :x2], makeunique=false) == [:x, :x_1, :x2]

    # Check that reserved words are up to date

    f = "$(Compat.Sys.BINDIR)/../../src/julia-parser.scm"
    if isfile(f)
        r1 = r"define initial-reserved-words '\(([^)]+)"
        r2 = r"define \(parse-block s(?: \([^)]+\))?\)\s+\(parse-Nary s (?:parse-eq '\([^(]*|down '\([^)]+\) '[^']+ ')\(([^)]+)"
        body = read(f, String)
        m1, m2 = match(r1, body), match(r2, body)
        if m1 == nothing || m2 == nothing
            error("Unable to extract keywords from 'julia-parser.scm'.")
        else
            s = replace(m1.captures[1]*" "*m2.captures[1], r";;.*?\n", "")
            rw = Set(split(s, r"\W+"))
            @test rw == DataFrames.RESERVED_WORDS
        end
    else
        @warn("Unable to validate reserved words against parser. ",
              "Expected if Julia was not built from source.")
    end

    @test DataFrames.countmissing([1:3;]) == 0

    data = Vector{Union{Float64, Missing}}(rand(20))
    @test DataFrames.countmissing(data) == 0
    data[sample(1:20, 11, replace=false)] = missing
    @test DataFrames.countmissing(data) == 11
    data[1:end] = missing
    @test DataFrames.countmissing(data) == 20

    pdata = Vector{Union{Int, Missing}}(sample(1:5, 20))
    @test DataFrames.countmissing(pdata) == 0
    pdata[sample(1:20, 11, replace=false)] = missing
    @test DataFrames.countmissing(pdata) == 11
    pdata[1:end] = missing
    @test DataFrames.countmissing(pdata) == 20

    funs = [mean, sum, var, x -> sum(x)]
    if string(funs[end]) == "(anonymous function)" # Julia < 0.5
        @test DataFrames._fnames(funs) == ["mean", "sum", "var", "λ1"]
    else
        @test DataFrames._fnames(funs) == ["mean", "sum", "var", string(funs[end])]
    end

    @testset "describe" begin
        io = IOBuffer()
        df = DataFrame(Any[collect(1:4), Vector{Union{Int, Missing}}(2:5),
                           CategoricalArray(3:6),
                           CategoricalArray{Union{Int, Missing}}(4:7)],
                       [:arr, :missingarr, :cat, :missingcat])
        describe(io, df)
        DRT = CategoricalArrays.DefaultRefType
        str =
            """
            arr
            Summary Stats:
            Mean:           2.500000
            Minimum:        1.000000
            1st Quartile:   1.750000
            Median:         2.500000
            3rd Quartile:   3.250000
            Maximum:        4.000000
            Length:         4
            Type:           $Int

            missingarr
            Summary Stats:
            Mean:           3.500000
            Minimum:        2.000000
            1st Quartile:   2.750000
            Median:         3.500000
            3rd Quartile:   4.250000
            Maximum:        5.000000
            Length:         4
            Type:           $(Union{Missing,Int})
            Number Missing: 0
            % Missing:      0.000000

            cat
            Summary Stats:
            Length:         4
            Type:           $(CategoricalValue{Int,DRT})
            Number Unique:  4

            missingcat
            Summary Stats:
            Length:         4
            Type:           $(Union{Missing, CategoricalValue{Int,DRT}})
            Number Unique:  4
            Number Missing: 0
            % Missing:      0.000000

            """
            out = String(take!(io))
            @test out == str
    end

    @testset "describe" begin
        io = IOBuffer()
        df = DataFrame(Any[collect(1:4), collect(Union{Int, Missing}, 2:5),
                           CategoricalArray(3:6),
                           CategoricalArray{Union{Int, Missing}}(4:7)],
                       [:arr, :missingarr, :cat, :missingcat])
        describe(io, df)
        str =
            """
            arr
            Summary Stats:
            Mean:           2.500000
            Minimum:        1.000000
            1st Quartile:   1.750000
            Median:         2.500000
            3rd Quartile:   3.250000
            Maximum:        4.000000
            Length:         4
            Type:           $Int

            missingarr
            Summary Stats:
            Mean:           3.500000
            Minimum:        2.000000
            1st Quartile:   2.750000
            Median:         3.500000
            3rd Quartile:   4.250000
            Maximum:        5.000000
            Length:         4
            Type:           $(Union{Missing, Int})
            Number Missing: 0
            % Missing:      0.000000

            cat
            Summary Stats:
            Length:         4
            Type:           $(CategoricalValue{Int,CategoricalArrays.DefaultRefType})
            Number Unique:  4

            missingcat
            Summary Stats:
            Length:         4
            Type:           $(Union{Missing, CategoricalValue{Int,CategoricalArrays.DefaultRefType}})
            Number Unique:  4
            Number Missing: 0
            % Missing:      0.000000

            """
        out = String(take!(io))
        @test out == str
    end
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

module Gaston

export closefigure, closeall, clearfigure, figure, plot, histogram, imagesc,
    surf, printfigure, gaston_demo, gaston_tests, set_terminal, set_filename,
    set_default_legend, set_default_plotstyle, set_default_color,
    set_default_marker, set_default_linewidth, set_default_pointsize,
    set_default_title, set_default_xlabel, set_default_ylabel,
    set_default_zlabel, set_default_box, set_default_axis, set_print_color,
    set_print_fontface, set_print_fontsize, set_print_fontscale,
    set_print_linewidth, set_print_size, gnuplot_send,
    addcoords, addconf, llplot, CurveConf, AxesConf, meshgrid

# before doing anything else, verify gnuplot is present on this system
if !success(`which gnuplot`)
    error("Gaston cannot be loaded: gnuplot is not available on this system.")
end
if readchomp(`gnuplot --version`)[1:11] != "gnuplot 4.6"
    println("Warning: Gaston has only been tested on gnuplot version 4.6")
end

# load files
include("gaston_types.jl")
include("gaston_aux.jl")
include("gaston_lowlvl.jl")
include("gaston_midlvl.jl")
include("gaston_hilvl.jl")
include("gaston_config.jl")
include("gaston_demo.jl")
include("gaston_test.jl")

# set up global variables
# global variable that stores gnuplot's state
gnuplot_state = GnuplotState(false,0,0,string("/tmp/gaston-",ENV["USER"],
"-",randstring(5),"/"),[])
# when gnuplot_state goes out of scope, close the pipe
finalizer(gnuplot_state,gnuplot_exit)

# global variable that stores Gaston's configuration
gaston_config = GastonConfig()

end