text:
#!/usr/bin/env julia
#=
   Author: Tim Sterne-Weiler, 3/16/2015
   e-mail: tim.sterne.weiler@utoronto.ca

=#

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ## INITIALIZATION  ## ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #

global head = "[aligater table]:"
println(STDERR, "$head Loading Packages..")

using ArgParse
using StatsBase

function parse_cmd()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--gi"
      help = "index for geneid:geneid field in foreground"
      arg_type = Int
      default = 1
    "--ci"
      help = "index for counting field if applicable (non integer fields increment by 1 instead)"
      arg_type = Int
      default = 2
    "--filt","-f"
      help = "pattern filter [syntax= column:regex]"
      arg_type = ASCIIString
      default = "1:."
    "--alpha","-a"
      help = "maximum p-value for output, only applied to the numerator (first element) of --nd if supplied"
      arg_type = Float64
      default = 1.0
    "--delim","-d"
      help = "delimiter for paired index entry"
      arg_type = ASCIIString
      default = ":"
    "--null","-n"
      help = "null value to insert for missing pairs"
      arg_type = ASCIIString
      default = "X"
    "--outsep","-s"
      help = "output separator for table"
      arg_type = ASCIIString
      default = ","
    "--integer","-i"
      help = "only count integers, try to convert if --ci is not an integer"
      action = :store_true
    "--norepeat","-r"
      help = "split names and remove _\\S+"
      action = :store_true
    "--sumOnly"
      action = :store_true
  end
  return parse_args(s)
end

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ### ###  FUNCTIONS ## ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #

include("dictext.jl")

# this accepts a --filt string argument in the form 'column:regexstring'
function parse_colfilt(toparse::ASCIIString)
  cap = match(r"(\d+)\:(.*)", toparse).captures
  @assert( length(cap) == 2 )
  col = parse(Int, cap[1])
  reg = Regex(cap[2])
  col, reg
end #--> Tuple{Integer, Regex}

# if the pairs --TODODOC
function load_interactionfile(io, gInd::Int, cntInd::Int, delim, col, reg; ctype=Float64, spltflag=false)
  cset = Dict{AbstractString,Dict{AbstractString,ctype}}()
  #sizehint!(cset, 100)
  sum = zero(ctype)
  for l in eachline( io )
    s = split(chomp(l), '\t')
    if isa(col, Integer) && isa(reg, Regex) && length(s) >= col
      ismatch(reg, s[col]) || continue
    end
    # make sure we don't get out of bounds error
    (ismatch(Regex(delim), s[gInd])) || continue
    k1,k2 = spltflag ? map( x->split(x, '_')[1], sort( split(s[gInd], delim) ) ) : sort( split(s[gInd], delim) )

    if ismatch(r"sense|transcript|protein|mRNA|sens|int|lncRNA", k1) || ismatch(r"snoRNA|miRNA", k2)
      k1,k2 = k2,k1
    end

    cnum = parse( s[cntInd], raise=false )
    if !haskey(cset, k1)
      cset[k1] = Dict{AbstractString,ctype}()
    end
    indict = cset[k1]
    toinc = isa(cnum, Number) ? convert(ctype, cnum) : one(ctype)
    sum += toinc
    dinc!(indict, k2, toinc)
  end
  cset,sum
end #--> Dict{ASCIIString,Dict{ASCIIString,ctype}},sum::ctype

# tested.
function print_table{K,V}( io, table::Dict{K,Dict{K,V}}; delim=",", null="X" )
 allk = Set{K}()
 for i::K in keys(table), j::K in keys(table[i])
   push!(allk, i)
   push!(allk, j)
 end
 header = "TAB" * delim * join( allk, delim )
 println( header )
 for i::K in allk
   print( io, i )
   if haskey( table, i )
     indict = table[i]
     for j::K in allk
       val = haskey(indict, j) ? string(indict[j]) : null
       print( io, delim * val )
     end
   else
     print( io, repeat( delim * null, length(allk) ) )
   end #if
   print( io, "\n" )
 end
end

## ### ### ### ### ### ### ### ### ### ### ### ### ## #
# ### ### ### ### ### MAIN  ### ### ### ### ### ### ##
## ### ### ### ### ### ### ### ### ### ### ### ### ## #
function main()
  pargs = parse_cmd()

  col,reg = parse_colfilt( pargs["filt"] )
  const mytype = pargs["integer"] ? Int64 : Float64;

  table,sum = load_interactionfile( STDIN, pargs["gi"], pargs["ci"], pargs["delim"], col, reg, ctype=mytype, spltflag=pargs["norepeat"] )
  pargs["sumOnly"] && println(STDOUT, "$(sum)")

  print_table( STDOUT, table; delim=pargs["outsep"], null=pargs["null"] )
end
#####################################################
# main execution here
main()
# eof

max_forks_repo_name: timbitz/aligater


text:
#=
functions to create index points from pointcloud along line defined by start & end point.

profiler generetes DataFrame with XYZ and attribute information from pointcloud,
where X and Y are defined along line

=#

struct Profile{U<:Integer,T<:AbstractFloat}
    pointindex::Vector{U}
    pointlx::Vector{T}
    pointly::Vector{T}
    maxdist::T
    pstart::Tuple{T, T}
    pend::Tuple{T, T}
    len::T
end

# length of profile line segment
linelength(p::Profile) = p.len
pointindex(p::Profile) = p.pointindex
pointlx(p::Profile) = p.pointlx
pointly(p::Profile) = p.pointly

# Base functions for Profile type
Base.length(p::Profile) = length(p.pointindex)

"index points along profileline"
function define_profile(cloud::Cloud, ps::Tuple{T,T}, pe::Tuple{T,T}, maxdist::T;
    pointfilter = nothing) where T <: Real

    # define line and bbox
    x_min, x_max = minmax(ps[1], pe[1])
    y_min, y_max = minmax(ps[2], pe[2])
    x_min -= maxdist
    y_min -= maxdist
    x_max += maxdist
    y_max += maxdist

    # define rotated grid orientation
    transrot = segment_transformation(ps, pe)
    len = hypot(ps[1] - pe[1], ps[2] - pe[2]) #length(ps, pe)

    # loop through Cloud and store points
    n = length(positions(cloud))
    idx = Int32[]  # index of points in Cloud
    lx = Float64[]  # distance along line
    ly = Float64[]  # distance from line
    for i in 1:n
        xi, yi, zi = positions(cloud)[i]
        outside_bbox(x_min, y_min, x_max, y_max, xi, yi) && continue
        (pointfilter != nothing) && (pointfilter(cloud, i) || continue)

        p = SVector(xi, yi)
        p′ = transrot(p)

        if (0.0 <= p′[1] < len) && (abs(p′[2]) <= maxdist)
            push!(idx, i)
            push!(lx, p′[1])
            push!(ly, p′[2])
        end
    end
    np = length(idx)
    # sort points based on distance along line lx
    if np >= 100_000 # RadixSort only faster with many points
        perm = sortperm(lx, alg=RadixSort)
    else
        perm = sortperm(lx)
    end

    Profile(idx[perm], lx[perm], ly[perm], maxdist, ps, pe, len)
end

"create new profile based on pointfilter applied to points in profile"
function filter_profile(cloud::Cloud, prof::Profile, pointfilter)
    n = length(prof)
    idx = Int32[]
    for i in 1:n
        idxi = pointindex(prof)[i]
        (pointfilter != nothing) && (pointfilter(cloud, i) || continue)
        push!(idx, idxi)
    end

    Profile(
        pointindex(prof)[idx],
        pointlx(prof)[idx],
        pointly(prof)[idx],
        prof.maxwidth,
        prof.pstart,
        prof.pend,
        prof.len)
end

# utils
"""Set up a transformation that will map points to a rotated
and translated coordinate system with the origin on the segment
start, and positive x along the segment."""
function segment_transformation(ps, pe)
    θ = -atan2(pe[2]-ps[2], pe[1]-ps[1])
    rot =  SMatrix{2,2}([cos(θ) -sin(θ); sin(θ) cos(θ)])
    trans = SVector(-ps[1], -ps[2])
    transrot(x::SVector{2,Float64}) = *(rot, x+trans)
end

"check if a point falls outside a bounding box"
function outside_bbox(x_min, y_min, x_max, y_max, x, y)
    x < x_min || y < y_min || x > x_max || y > y_max
end

ext: jl


text:
## Docs-specific rendering --------------------------------------------------------------

function writemime(io::IO, mime::MIME"text/plain", docs::Docile.Interface.Docs{:md})
    Markdown.term(io, parsed(docs))
end

function writemime(io::IO, mime::MIME"text/plain", docs::Docile.Interface.Docs{:txt})
    println(io)
    for line in split(parsed(docs), "\n")
        println(io, "  ", line)
    end
end

## General plain-text rendering - REPL --------------------------------------------------

function writemime(io::IO, mime::MIME"text/plain", qr::QueryResults)
    single = length(qr.matches) ≡ 1
    count = 1
    msg = false
    for (score, entries) in reverse!(sort([(a, b) for (a, b) in qr.scores]))
        for entry in entries
            match = qr.matches[entry]
            count, msg = writemime(io, mime, match, count, qr.query.index, single)
        end
    end
    if msg && HELP_MESSAGE[1]
        println(io, """

        Note: To display documentation for a specific entry listed above append
        the corresponding number to your previous query, ie:

            help?> "help" 1""")
        HELP_MESSAGE[1] = false
    end
end

const HELP_MESSAGE = [true]

function writemime(io::IO, mime::MIME"text/plain", m::Match, count, index, single)
    cat = colorize(:blue, "\n[$(category(m.entry))]\n")
    msg = false
    if single
        # Just one entry was found to match the query.
        println(io, cat)
        for object in sort(collect(m.objects); by = x -> string(x))
            print_signature(io, object, m.entry)
            count += 1
        end
        writemime(io, mime, m.entry)
    elseif 0 < index
        # Only show the nth entry from summary output.
        for object in m.objects
            if index == count
                println(io, cat)
                print_signature(io, object, m.entry)
                writemime(io, mime, m.entry)
            end
            count += 1
        end
    else
        # Display a summary of the entries that match a query.
        for object in m.objects
            print(io, colorize(:default, lpad(count, 3) * ": "))
            print_signature(io, object, m.entry)
            count += 1
        end
        msg = true
    end
    count, msg
end

function print_signature(io::IO, object, entry)
    color = isexported(modulename(entry), object) ? :green : :default
    print(io, colorize(color, join(fullname(modulename(entry)), ".")), ".")
    println(io, colorize(:default, writeobj(object, entry)))
end

function writemime(io::IO, mime::MIME"text/plain", entry::Entry, config::Config = Config())
    # Parse docstring into AST and print it out with additional space above and below.
    println(io)
    writemime(io, mime, docs(entry))
    println(io)

    # Print metadata if any is available
    if has_output_metadata(entry, config)
        println(io, colorize(:green, " Details:\n"))
        for m in config.metadata_order
            if haskey(entry.data, m)
                v = entry.data[m]
                if isa(v, Vector)
                    println(io, "\t", m, ":")
                    for line in v
                        if isa(line, NTuple)
                            println(io, "\t\t", colorize(:cyan, string(line[1])), ": ", line[2])
                        else
                            println(io, "\t\t", string(line))
                        end
                    end
                else
                    println(io, "\t", m, ": ", v)
                end
            end
        end
    end
end

function colorize(color::Symbol, text::AbstractString)
    buf = IOBuffer()
    print_with_color(color, buf, text)
    takebuf_string(buf)
end


text:
# Unit tests for analysis/visualization.jl.jl

@testset "analysis/visualization.jl.jl" begin
    # TODO
end


text:
module EnvInfo

using JSON

struct PkgStatusInfo
    name::String
    version::VersionNumber
    revision::String
    is_clean::Bool
end

function status(pkg::AbstractString,
                version::VersionNumber = Pkg.installed(pkg))
    return cd(Pkg.dir(pkg)) do
        PkgStatusInfo(
            pkg,
            version,
            strip(read(`git rev-parse HEAD`, String)),
            strip(read(`git status --short --untracked-files=no`,
                       String)) == "",
        )
    end
end

status(packages::AbstractVector{<: AbstractString}) = map(status, packages)
status() = [status(pkg, v) for (pkg, v) in Pkg.installed()]

function envinfo(packages; include_installed = false)
    info = Dict{String, Any}(
        "status" => status(packages),
    )
    if include_installed
        info["installed"] = Pkg.installed()
    end
    return info
end

function save(filename::AbstractString, packages::AbstractVector;
              assert_clean = false,
              kwargs...)
    info = envinfo(packages; kwargs...)

    if assert_clean
        unclean = []
        for psi in info["status"]
            if ! psi.is_clean
                push!(unclean, psi.name)
            end
        end
        if length(unclean) > 0
            names = join(unclean, ", ")
            error("Package(s) $names are not clean.")
        end
    end

    open(filename, "w") do io
        JSON.print(io, info)
    end
end

# package code goes here

end # module


text:
# XpressMathProgInterface
# Standardized MILP interface
using MathProgBase
using MathProgBase.SolverInterface
const MPB = MathProgBase.SolverInterface

mutable struct XpressMathProgModel <: AbstractLinearQuadraticModel
    inner::Model

    lazycb
    cutcb
    heuristiccb
    infocb

    options
end
function XpressMathProgModel(;options...)

   m = XpressMathProgModel(Model(), nothing, nothing, nothing, nothing, options)

   for (name,value) in options
       setparam!(m.inner, XPRS_CONTROLS_DICT[name], value)
   end
   return m
end

function copy(m::XpressMathProgModel)

    m.lazycb == nothing || Compat.@warn("Callbacks can't be copied, lazy callback ignored")
    m.cutcb == nothing || Compat.@warn("Callbacks can't be copied, cut callback ignored")
    m.heuristiccb == nothing || Compat.@warn("Callbacks can't be copied, heuristic callback ignored")
    m.infocb == nothing || Compat.@warn("Callbacks can't be copied, info callback ignored")

    return XpressMathProgModel(copy(m.inner),
                                nothing,
                                nothing,
                                nothing,
                                nothing,
                                deepcopy(m.options))
end

mutable struct XpressSolver <: AbstractMathProgSolver
    options
end
XpressSolver(;kwargs...) = XpressSolver(kwargs)
MPB.LinearQuadraticModel(s::XpressSolver) = XpressMathProgModel(;s.options...)
MPB.ConicModel(s::XpressSolver) = MPB.LPQPtoConicBridge(LinearQuadraticModel(s))

MPB.supportedcones(m::XpressSolver) = [:Free, :Zero, :NonNeg, :NonPos, :SOC, :SOCRotated]

MPB.loadproblem!(m::XpressMathProgModel, filename:: String) = read_model(m.inner, filename)

function MPB.setparameters!(m::XpressMathProgModel; mpboptions...)
    for (optname, optval) in mpboptions
        if optname == :TimeLimit
            setparam!(m.inner, XPRS_MAXTIME, round(Cint, optval) )
        elseif optname == :Silent
            if optval == true
                setparam!(m.inner, XPRS_OUTPUTLOG, 0)
            end
        else
            error("Unrecognized parameter $optname")
        end
    end
end

MPB.loadproblem!(m::Model, filename::String) = read_model(m.inner, filename)

function MPB.loadproblem!(m::XpressMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
    # throw away old model
    #env = m.inner.env
    #m.inner.finalize_env = false
    free_model(m.inner)

    m.inner = Model()

    for (name,value) in m.options
        setparam!(m.inner, XPRS_CONTROLS_DICT[name], value)
    end

    add_cvars!(m.inner, float(obj), float(collb), float(colub))

    neginf = typemin(eltype(rowlb))
    posinf = typemax(eltype(rowub))

    # check if we have any range constraints
    # to properly support these, we will need to keep track of the
    # slack variables automatically added by xpress.
    rangeconstrs = false#sum((rowlb .!= rowub) .& (rowlb .> neginf) .& (rowub .< posinf))
    for i in eachindex(rowlb)
        if rowlb[i] != rowub[i] && rowlb[i] > neginf && rowub[i] < posinf
            rangeconstrs = true
            break
        end
    end
    if rangeconstrs
        Compat.@warn("Julia Xpress interface doesn't properly support range (two-sided) constraints.")
        add_rangeconstrs!(m.inner, float(A), float(rowlb), float(rowub))
    else
        b = Array{Float64}(undef, length(rowlb))
        senses = Array{Cchar}(undef, length(rowlb))
        for i in 1:length(rowlb)
            if rowlb[i] == rowub[i]
                senses[i] = XPRS_EQ#'='
                b[i] = rowlb[i]
            elseif rowlb[i] > neginf
                senses[i] = XPRS_GEQ#'>'
                b[i] = rowlb[i]
            else
                @assert rowub[i] < posinf
                senses[i] = XPRS_LEQ#'<'
                b[i] = rowub[i]
            end
        end
        add_constrs!(m.inner, float(A), senses, b)
    end

    setsense!(m,sense)
    return nothing
end

MPB.writeproblem(m::XpressMathProgModel, filename:: String) = write_model(m.inner, filename)

MPB.getvarLB(m::XpressMathProgModel)     = get_lb(m.inner)
MPB.setvarLB!(m::XpressMathProgModel, l) = set_lb!(m.inner,l)

MPB.getvarUB(m::XpressMathProgModel)     = get_ub(m.inner)
MPB.setvarUB!(m::XpressMathProgModel, u) = set_ub!(m.inner,u)

function MPB.getconstrLB(m::XpressMathProgModel)
    sense = get_rowtype(m.inner)
    ret   = get_rhs(m.inner)
    for i = 1:num_constrs(m.inner)
        if sense[i] == XPRS_GEQ || sense[i] == XPRS_EQ
            # Do nothing
        else
            # LEQ constraint, so LB is -Inf
            ret[i] = -Inf
        end
     end
     return ret
end

function MPB.getconstrUB(m::XpressMathProgModel)
    sense = get_rowtype(m.inner)
    ret   = get_rhs(m.inner)
    for i = 1:num_constrs(m.inner)
        if sense[i] == XPRS_LEQ || sense[i] == XPRS_EQ
            # Do nothing
        else
            # GEQ constraint, so UB is +Inf
            ret[i] = +Inf
        end
    end
    return ret
end

MPB.setconstrLB!(m::XpressMathProgModel, lb) = set_constrLB!(m.inner,lb)
MPB.setconstrUB!(m::XpressMathProgModel, ub) = set_constrUB!(m.inner,ub)

MPB.getobj(m::XpressMathProgModel)     = get_obj(m.inner)
MPB.setobj!(m::XpressMathProgModel, c) = set_obj!(m.inner,c)

MPB.addvar!(m::XpressMathProgModel, constridx, constrcoef, l, u, objcoef) = add_var!(m.inner, length(constridx), constridx, float(constrcoef), float(objcoef), float(l), float(u))
MPB.addvar!(m::XpressMathProgModel, l, u, objcoef) = add_cvar!(m.inner, float(objcoef), float(l), float(u))

MPB.delvars!(m::XpressMathProgModel, idx) = del_vars!(m.inner, idx)

function MPB.addconstr!(m::XpressMathProgModel, varidx, coef, lb, ub)
    if lb == -Inf
        # <= constraint
        add_constr!(m.inner, varidx, coef, XPRS_LEQ, ub)
    elseif ub == +Inf
        # >= constraint
        add_constr!(m.inner, varidx, coef, XPRS_GEQ, lb)
    elseif lb == ub
        # == constraint
        add_constr!(m.inner, varidx, coef, XPRS_EQ, lb)
    else
        # Range constraint
        error("Adding range constraints not supported yet.")
    end
end

MPB.delconstrs!(m::XpressMathProgModel, idx) = del_constrs!(m.inner, idx)
MPB.changecoeffs!(m::XpressMathProgModel, cidx, vidx, val) = chg_coeffs!(m.inner, cidx, vidx, val)

MPB.getconstrmatrix(m::XpressMathProgModel) = get_constrmatrix(m.inner)

function MPB.setsense!(m::XpressMathProgModel, sense)
    if sense == :Min
        set_sense!(m.inner, :minimize)
    elseif sense == :Max
        set_sense!(m.inner, :maximize)
    else
        error("Unrecognized objective sense $sense")
    end
end

function MPB.getsense(m::XpressMathProgModel)
    if model_sense(m.inner) == :maximize
        return :Max
    else
        return :Min
    end
end

MPB.numvar(m::XpressMathProgModel)    = num_vars(m.inner)
MPB.numconstr(m::XpressMathProgModel) = num_constrs(m.inner)
MPB.numlinconstr(m::XpressMathProgModel) = num_constrs(m.inner) - num_qconstrs(m.inner)
MPB.numquadconstr(m::XpressMathProgModel) = num_qconstrs(m.inner)

MPB.optimize!(m::XpressMathProgModel) = optimize(m.inner)

function MPB.status(m::XpressMathProgModel)
  is_mip(m.inner) ? s = get_mip_status(m.inner) : s = get_lp_status(m.inner)

  if s == :optimal
    return :Optimal
  elseif s == :infeasible
    return :Infeasible
  elseif s == :unbounded
    return :Unbounded
  elseif s == :cutoff
    return :CutOff
  elseif s == :unfinished
    return :Unfinished
  elseif s == :cutoff_in_dual
    return :CutOff_In_Dual
  elseif s == :unsolved
    return :Unsolved
  elseif s == :nonconvex
    return :Nonconvex
  #mip
  elseif s == :mip_optimal
    return :Optimal
  elseif s == :mip_infeasible
    return :Infeasible
  elseif s == :mip_unbounded
    return :Unbounded
  elseif s == :mip_not_loaded
    return :Unfinished
  elseif s == :mip_lp_not_optimal
    return :Unfinished
  elseif s == :mip_lp_optimal
    return :Unfinished
  elseif s == :mip_no_sol_found
    return :Unfinished
  elseif s == :mip_suboptimal
    return :Unfinished
#old
  elseif s == :iteration_limit || s == :node_limit || s == :time_limit || s == :solution_limit
    return :UserLimit
  elseif s == :numeric
    return :Numeric
  elseif s == :suboptimal
    return :Suboptimal # not very useful status
  elseif s == :interrupted # ended by user?
    return :UserLimit
  else
    error("Unrecognized solution status: $s")
  end
end

MPB.getobjval(m::XpressMathProgModel)   = get_objval(m.inner)
MPB.getobjbound(m::XpressMathProgModel) = get_bestbound(m.inner) #get_objbound(m.inner)
MPB.getsolution(m::XpressMathProgModel) = get_solution(m.inner)

function MPB.getconstrsolution(m::XpressMathProgModel)
    rows = num_constrs(m.inner)
    sense = get_rowtype(m.inner)
    rhs   = get_rhs(m.inner)
    slack = get_slack(m.inner)
    ret   = zeros(rows)
    for i = 1:rows
        if sense[i] == XPRS_EQ
            # Must be equal to RHS if feasible
            ret[i] = rhs[i]
        else
            # Slack variable is non-negative for <= constraints, and non-positive for >= constraints
            ret[i] = rhs[i] - slack[i]
        end
    end
    return ret
end

function MPB.getreducedcosts(m::XpressMathProgModel)
    if is_qcp(m.inner) #&& get_int_param(m.inner.env, "QCPDual") == 0
        return get_reducedcost(m.inner)
        #return fill(NaN, num_vars(m.inner))
    else
        return get_reducedcost(m.inner)
    end
end

function MPB.getconstrduals(m::XpressMathProgModel)
    if is_qcp(m.inner) #&& get_int_param(m.inner.env, "QCPDual") == 0

        nlrows = num_linconstrs(m.inner)
        nqrows = num_qconstrs(m.inner)

        lduals = Array{Float64}(undef,  nlrows)

        qrows = get_qrows(m.inner)

        duals = get_dual(m.inner)

        pos_lr = 1
        pos_qr = 1
        for r = 1:length(duals)

            if r != qrows[pos_qr]
                lduals[pos_lr] = duals[r]
                pos_lr += 1
            else
                if pos_qr < nqrows
                    pos_qr += 1
                end
            end
        end

        return lduals
        #return fill(NaN, num_constrs(m.inner))
    else
        return get_dual(m.inner)
    end
end

function MPB.getquadconstrduals(m::XpressMathProgModel)
    if is_qcp(m.inner) #&& get_int_param(m.inner.env, "QCPDual") == 0

        #nlrows = num_linconstrs(m.inner)
        nqrows = num_qconstrs(m.inner)

        qduals = Array{Float64}(undef,  nqrows)

        qrows = get_qrows(m.inner)

        duals = get_dual(m.inner)

        pos_qr = 1
        for r = 1:length(duals)

            if r == qrows[pos_qr]
                qduals[pos_qr] = duals[r]
                pos_qr += 1
            end
            if pos_qr > nqrows
                break
            end
        end

        return qduals
        #return fill(NaN, num_qconstrs(m.inner))
    else
        return Float64[]
    end
end

MPB.getinfeasibilityray(m::XpressMathProgModel) = getdualray(m.inner)
MPB.getunboundedray(m::XpressMathProgModel) = getprimalray(m.inner)

MPB.getbasis(m::XpressMathProgModel) = get_basis(m.inner)

MPB.getrawsolver(m::XpressMathProgModel) = m.inner

const var_type_map = Dict{Cchar,Symbol}(
  convert(Cchar, 'C') => :Cont,
  convert(Cchar, 'B') => :Bin,
  convert(Cchar, 'I') => :Int,
  convert(Cchar, 'S') => :SemiCont,
  convert(Cchar, 'R') => :SemiInt
)

const rev_var_type_map = Dict{Symbol,Cchar}(
  :Cont => Cchar('C'),
  :Bin => Cchar('B'),
  :Int => Cchar('I'),
  :SemiCont => Cchar('S'),
  :SemiInt => Cchar('R')
)

function MPB.setvartype!(m::XpressMathProgModel, vartype::Vector{Symbol})

    # dealing with semi continuous lower bounds
    # see issue #2
    semivars = false
    if :SemiInt in vartype
        semivars = true
    end
    if :SemiCont in vartype
        semivars = true
    end

    # do this to make sure we deal with new columns
    ind = collect( 1:length(vartype) )

    nvartype = map(x->rev_var_type_map[x], vartype)

    lb = semivars ? getvarLB(m) : Float64[] # for semi vars

    chgcoltypes!(m.inner, ind, nvartype )

    if semivars

        semi_col = Int[]
        semi_lb = Float64[]

        for i in 1:length(vartype)
            if vartype[i] == :SemiCont || vartype[i] == :SemiInt
                push!(semi_col, i)
                push!(semi_lb, lb[i])
            end
        end

        chgsemilb!(m.inner, semi_col, semi_lb)
    end

    nothing
end

function MPB.getvartype(m::XpressMathProgModel)
    ret = get_coltype(m.inner)
    return map(x->var_type_map[x], ret)
end

MPB.setwarmstart!(m::XpressMathProgModel, v) = loadbasis(m.inner, v)
MPB.addsos1!(m::XpressMathProgModel, idx, weight) = add_sos!(m.inner, :SOS1, idx, weight)
MPB.addsos2!(m::XpressMathProgModel, idx, weight) = add_sos!(m.inner, :SOS2, idx, weight)

# QCQP

function MPB.setquadobj!(m::XpressMathProgModel, rowidx, colidx, quadval)
    delq!(m.inner)

    # xpress only accept one input per matrix slot
    k = ones(Bool,length(rowidx)) #replicates holder (they are falses)
    for i in 1:length(rowidx), j in (i+1):length(rowidx)
        if ( rowidx[i] == rowidx[j] && colidx[i] == colidx[j] ) || ( rowidx[i] == colidx[j] && colidx[i] == rowidx[j] )
            quadval[i] += quadval[j]
            quadval[j] = 0.0
            rowidx[j] = -1
            colidx[j] = -1
            k[j] = false
        end
    end

    scaledvals = similar(quadval)
    for i in 1:length(rowidx)
        if rowidx[i] == colidx[i]
            # rescale from matrix format to "terms" format
            scaledvals[i] = quadval[i] #/ 2
        else
            scaledvals[i] = quadval[i]#*2
        end
    end
    add_qpterms!(m.inner, rowidx[k], colidx[k], scaledvals[k])

    return nothing
end

function MPB.addquadconstr!(m::XpressMathProgModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    # xpress only accept one input per matrix slot
    k = ones(Bool,length(quadrowidx)) #replicates holder (they are falses)
    for i in 1:length(quadrowidx), j in (i+1):length(quadrowidx)
        if ( quadrowidx[i] == quadrowidx[j] && quadcolidx[i] == quadcolidx[j] ) || ( quadrowidx[i] == quadcolidx[j] && quadcolidx[i] == quadrowidx[j] )
            quadval[i] += quadval[j]
            quadval[j] = 0.0
            #rowidx[j] = -1
            #colidx[j] = -1
            k[j] = false
        end
    end

    for i in 1:length(quadrowidx)
        if quadrowidx[i] != quadcolidx[i]
            quadval[i] = quadval[i]/2
        end
    end

    add_qconstr!(m.inner, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    return nothing
end

MPB.getsolvetime(m::XpressMathProgModel) = m.time
MPB.getnodecount(m::XpressMathProgModel) = get_node_count(m.inner)


text:
function conic_form!(s::Positive, x::AbstractVariable, unique_conic_forms)
    return conic_form!(x >= 0, unique_conic_forms)
end
function conic_form!(s::Negative, x::AbstractVariable, unique_conic_forms)
    return conic_form!(x <= 0, unique_conic_forms)
end


text:
import Oceananigans.TimeSteppers: calculate_pressure_correction!, pressure_correct_velocities!

calculate_pressure_correction!(::HydrostaticFreeSurfaceModel, Δt) = nothing

#####
##### Barotropic pressure correction for models with a free surface
#####

pressure_correct_velocities!(model::HydrostaticFreeSurfaceModel{T, E, A, <:ExplicitFreeSurface}, Δt) where {T, E, A} = nothing

#####
##### Barotropic pressure correction for models with a free surface
#####

function pressure_correct_velocities!(model::HydrostaticFreeSurfaceModel{T, E, A, <:ImplicitFreeSurface}, Δt) where {T, E, A}

    event = launch!(model.architecture, model.grid, :xyz,
                    _barotropic_pressure_correction,
                    model.velocities,
                    model.grid,
                    Δt,
                    model.free_surface.gravitational_acceleration,
                    model.free_surface.η,
                    dependencies = Event(device(model.architecture)))

    return nothing
end

@kernel function _barotropic_pressure_correction(U, grid, Δt, g, η)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U.u[i, j, k] -= g * Δt * ∂xᶠᶜᵃ(i, j, k, grid, η)
        U.v[i, j, k] -= g * Δt * ∂yᶜᶠᵃ(i, j, k, grid, η)
    end
end


text:
using Test
using PolygonInbounds
using .Demo

@testset "PolygonInboundsTests" begin
@testset "argument checks" begin

    @test inpoly2([0.0 0], [0.0 0]) == [false true]
    @test_throws ArgumentError inpoly2([0.0 0], [0.0 0], [1;2])
    @test inpoly2([0.0 0], [-1.0 1; 1 -1]) == [false true]
    @test inpoly2([0.0 0], [-1.0 1; 1 -1], [1; 2]) == [false false]
    @test inpoly2([0.0 0], [-1.0 1; 1 -1], [2; 1]) == [false true]
    @test_throws ArgumentError inpoly2([0.0 0], [-1.0 1; 1 -1], [2 1])
    @test inpoly2([0.0 0], [-1.0 1; 1 -1], [1 1]) == [false false]
    @test_throws ArgumentError inpoly2([0.0 0], [-1.0 1; 1 -1], [2 1; 2 2])
    @test_throws ArgumentError inpoly2([0.0 0], [-1.0 1; 1 -1], [2 2; 1 2])
    @test_throws ArgumentError  inpoly2([0.0 0], [-1.0 1; 1 -1], [1 1 1 2 3 4])
end

let stat = polydemo(4, r=10^5, tol=1e-2, do_plot=false)
    @testset "functionality" begin
        @test 300 <= count(stat[:,2]) < 600
        @test 29000 <= count(stat[:,1]) < 32000
        @test 67000 <= count((!).(stat[:,1])) < 70000
    end
end

@testset "utilities" begin
    @test PolygonInbounds.mean_overlap(0.0, 0.0, 2.0, 8.0, 10.0) == 0
    @test PolygonInbounds.mean_overlap(2.0, 0.0, 2.0, 8.0, 10.0) == 1.5
    @test PolygonInbounds.mean_overlap(8.0, 0.0, 2.0, 8.0, 10.0) == 6
    @test PolygonInbounds.mean_overlap(11.0, 0.0, 2.0, 8.0, 10.0) == 6
end

@testset "means estimations" begin
    r = 100
    nodes = [0.0 0; 3 4; 7 5; 10 0]
    edges = [1 2; 2 3; 3 4; 4 1]
    poly = PolygonInbounds.PolygonMesh(nodes, edges)
    rpts = PolygonInbounds.PointsInbound(Demo.testbox(r, [-1.0 -1.0], [ 11.0 6.0]))
    x = PolygonInbounds.estimate_effort(rpts, [-1.0 -1.0], [11.0 6.0], poly, [0.0 0.0], [10.0 5.0], 1e-6)
    @test length(x) == 4
end

end

# to improve test coverage
for i = 1:5
    polydemo(i)
end


text:
"""
    locate(what, ofs::UInt32, size::UInt64, buffer::String)

offset-oriented high performance search function without matcher.
searches what in codeunits(buffer,ofs+1:ofs+size).

Special treatment of searches for single bytes and byte pairs.
longer searches use a bloom filter algorithm, an improved version
of the algorithm used in Base._searchindex (in Julia Release 1.4).

"""
function locate end

locate(what::Union{AbstractString,AbstractChar},ofs::UInt32, size::UInt64, buffer::String) = locate(BToken(what), ofs, size, where)

function locate(what::UInt8, ofs::UInt32, size::UInt64, buffer::String)
    @boundscheck checkulimit(ofs+size,usize(buffer))
    p = pointer(buffer)
    q = GC.@preserve buffer ccall(:memchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), p+ofs, b, size%Csize_t)
    return q == C_NULL ? typemax(UInt32) : UInt32(q-p)
end

"very simple hash function used for bloom filters in ExactMatcher and locate"
bloomhash(b::UInt8) = UInt64(1) << (b & 0x3f)

"maximum number of bits allowed to be set in bloom filter"
const BLOOMFILTERLIMIT = 48

function locate(what::BToken, ofs::UInt32, size::UInt64, buffer::String)
    @boundscheck checkulimit(ofs+size,usize(buffer))
    sw = usize(what)
    size<sw && return typemax(UInt32)
    if sw<=2
        sw==1 && return @inbounds locate(byte(what,0),ofs,size,buffer)
        sw==0 && return ofs
        # we search for a byte pair.
        b0 = byte(what,0)
        b1 = byte(what,1)
        o2 = 0%UInt32 # o2 is offset within content to search
        pb = pointer(buffer)+ofs
        while size-1-o2 > 0
            # byte search - higher order methods do not pay off by 2-byte-search
            q = GC.@preserve buffer ccall(:memchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), pb+o2, b0, (size-1-o2)%Csize_t)
            if q == C_NULL
                return typemax(UInt32)
            end
            o2 = (q-pb)%UInt32
            if b1 == byte(buffer,ofs+o2+1)
                return ofs+o2
            end
        end
        return typemax(UInt32)
    end
    # standard case: at least 3 bytes to search, buffer size >= pattern size.
    # we use an improved version of Base._searchindex with bloom filter. see https://github.com/rryi/SearchBenchmarks.jl
    #
    # (1) initialization: copy/paste from Exactmatcher initialization
    skip = sw # bytes to skip if last byte matches, any other byte does not match
    last = codeunit(what,sw)
    bloom_mask = bloomhash(last)
    bloom_skip = sw # no. of bytes to skip if byte not in filter
    j = sw
    if sw<=BLOOMFILTERLIMIT # just fillup bloom filter. No degeneration countermeasures necessary
        bloom_skip = sw
        while (j-=1)>=1  # j=n is already processed
            pj = codeunit(what,j)
            if pj == last && skip == sw
                skip = sw - j
            end
            bloom_mask |= bloomhash(pj)
        end
    else
        # long pattern: do not overfill bloom filter which would result in degeneration
        bloom_bits = 1 # no. of bits set in bloom filter if <= 32
        bloom_skip = 1 # we have to count upwards until no. of bytes to skip if byte not in filter
        while (j-=1)>=1  # j=n is already processed
            pj = codeunit(what,j)
            if pj == last && skip == sw
                skip = sw - j
            end
            if bloom_bits <= BLOOMFILTERLIMIT
                # we continue setting bits in bloom filter. BLOOMFILTERLIMIT should be a guess for max(p(bloom_skip)*bloom_skip)
                hash = bloomhash(pj)
                if hash&bloom_mask == 0
                    if bloom_bits < BLOOMFILTERLIMIT
                        # put in bloom filter up to 32 bits
                        bloom_mask |= hash
                    end
                    bloom_bits += 1 # gets >=BLOOMFILTERLIMIT to terminate bloom filter filling
                end
                if bloom_bits <= BLOOMFILTERLIMIT
                    bloom_skip += 1
                end
            else
                # do not put more bits into bloom filter.
                if skip<sw
                    # we have finished adding hashes to bloom filter
                    # and we have determined the skip distance if matching last byte
                    # nothing remains to be done in preprocessing - stop work.
                    break
                end
            end
        end
    end
    skip -= 1 # because in main loop, 1 byte of skip is skipped separately
    #
    # (2) kernel: loop through buffer until found / end of data
    w = size-sw # we already checked for size>=sw
    # w<0 implies: not found
    i = ofs + sw # index of last byte in matching candidate in search
    while i < size # i==size could be a match, is treated after this loop
        if codeunit(buffer,i) == last
            # check candidate
            j = 1
            while j < sw
                if codeunit(buffer,i-sw+j) != codeunit(what,j)
                    break
                end
                j += 1
                # match found?
                if j == sw
                    return i-n
                end
            end
            # no match: skip and test bloom
            i += skip
            if i>=size # main loop finished?
                break
            end
        end
        i += 1
        if bloom_mask & bloomhash(codeunit(buffer,i)) == 0
            i += bloomskip
        end
    end
    if i==size
        # test end match
        j = 1
        while j <= sw
            if codeunit(buffer,i-sw+j) != codeunit(what,j)
                break # not found
            end
            j += 1
        end
        if j > n # means: loop without break , we have a match at the very end
            return i-sw
        end
    end
    # not found
    return typemax(UInt32)
end

# TEST pair performance ohne C code (danach löschen)
function locate2(what::BToken, ofs::UInt32, size::UInt64, buffer::String)
    @boundscheck checkulimit(ofs+size,usize(buffer))
    sw = usize(what)
    size<sw && return typemax(UInt32)
    if sw==2
        # we search for a byte pair.
        b0 = byte(what,0)
        b1 = byte(what,1)
        o2 = 0%UInt32 # o2 is offset within content to search
        pb = pointer(buffer)+ofs
        pend = pb+size-1 # points to last byte in buffer
        p2 = pb
        GC.@preserve buffer while p2<pb
            if unsafe_load(p2) == b0
                if unsafe_load(p2+1) == b1
                    return ofs+(p2-pb)
                end
            end
            p2 += 1
        end # while
        return typemax(UInt32)
    end
    return locate(what,ofs,size,buffer)
end


text:
#### GEInfo #########

"
Auxiliary structure that keeps some generic information for Grammatical Evolution experiments.

# Fields
- `_grammar::Grammar`: grammar that will be used. This includes terminal and non-terminal
    symbols, production rules and the initial symbol.
- `_maxProductions::Int16`: maximum number of productions permitted to generate
    trees.
- `_maxDepth::Int16`: maximum tree depth permitted.
- `_variables::Array{String}`: array that contains the names of the variables of
    the problem.

See also: [`GPExperimentInfo`](@ref), [`Grammar`](@ref)
"
struct GEInfo <: GPExperimentInfo
    _grammar::Grammar
    _maxProductions::Int16
    _maxDepth::Int16
    _variables::Array{String}
end # struct

#### GEREP ##########

"
Structure that is the representantion of an individual in Grammatical Evolution,
which is a tree.

# Fields
- `_symbol::UInt16`: symbol in the grammar.
- `_depth::UInt8`: actual depth of the node in the tree.
- `_children::Array{GERep}`: children of symbol if it is a non terminal.

See also: [`GEGenotype`](@ref)
"
mutable struct GERep
    _symbol::UInt16
    _depth::UInt8
    _children::Array

    GERep(symbol::UInt16) = new(symbol, 0, Array{GERep}(undef, 0))
    GERep(symbol::UInt16, depth::UInt8) = new(symbol, depth, Array{GERep}(undef, 0))
end

"""
    getSymbol(tree::GERep)::UInt16

Returns the symbol stored.
"""
getSymbol(tree::GERep)::UInt16 = tree._symbol

"""
    getDepth(tree::GERep)::UInt8

Returns the depth of the node in the tree.
"""
getDepth(tree::GERep)::UInt8 = tree._depth

"""
    getChildren(tree::GERep)::Array{GERep}

Returns the array of children.
"""
getChildren(tree::GERep)::Array{GERep} = tree._children

"""
    setDepth!(tree::GERep, value)

Sets the depth of the node in the tree.
"""
setDepth!(tree::GERep, value) = tree.depth = depth

"""
    setDepth!(tree::GERep, value)

Sets the children of the node.
"""
setChildren!(tree::GERep, children::Array{GERep}) = tree._children = children

"""
    pushChild!(tree::GERep, child::GERep)

Introduce a child in the array of children.
"""
pushChild!(tree::GERep, child::GERep) = push!(getChildren(tree), child)

"""
    getNChildren(tree::GERep)::Int

Returns the number of children.
"""
getNChildren(tree::GERep)::Int = length(getChildren(tree))

"""
    isLeaf(tree::GERep)::Bool

Checks whether the node is a leaf or not.
"""
isLeaf(tree::GERep)::Bool = getNChildren(tree::GERep) == 0

"""
    getHeights(tree::GERep)

Obtains the heights of the nodes in preorder.
!!! warning
    DONT USE, it is still experimental.
"""
function getHeights(tree::GERep)

    representation = getPreOrder(tree)
    lenRep = length(representation)
    heights = Array{UInt16}(undef, lenRep)
    heights[end] = 0
    stack = [lenRep]
    currIndex = lenRep - 1

    while currIndex > 0

        if isLeaf(representation[currIndex])
            heights[currIndex] = 0
        else
            childrenIndexes = Array{UInt8}(undef, 0)
            for i=1:getNChildren(representation[currIndex])
                push!(childrenIndexes, pop!(stack))
            end
            heights[currIndex] = maximum(heights[childrenIndexes]) + 1
        end

        push!(stack, currIndex)
        currIndex -= 1
    end
    return heights
end

"""
    changeNodes!(tree1::GERep, tree2::GERep)

Interchange the given nodes.
"""
function changeNodes!(tree1::GERep, tree2::GERep)
    tree1._symbol, tree2._symbol = tree2._symbol, tree1._symbol
    tree1._children, tree2._children = tree2._children, tree1._children
    tree1._depth, tree2._depth = tree2._depth, tree1._depth
end

"""
    copyTree(tree::GERep)::GERep

Makes a copy of the tree.
"""
function copyTree(tree::GERep)::GERep

    stack = [tree]
    newTree = GERep(tree._symbol, tree._depth)
    copyStack = [newTree]

    while !isempty(stack)
        node = popfirst!(stack)
        newnode = popfirst!(copyStack)
        if !isLeaf(node)
            children = getChildren(node)
            for i in eachindex(children)
                newChild = GERep(children[i]._symbol, children[i]._depth)
                pushChild!(newnode, newChild)
                push!(stack, children[i])
                push!(copyStack, newChild)
            end
        end
    end
    return newTree
end

"""
    printTree(tree::GERep, gramm::Grammar)

Print the tree formated in the current output.
"""
function printTree(tree::GERep, gramm::Grammar)
    stack = [tree]
    result = []
    depth = 0

    function printRecursive(node, depth)
        for i=1:depth-1
            print("   ")
        end
        print("|--")
        if isLeaf(node)
            println(getRepresentation(gramm, node))
        else
            println(getRepresentation(gramm, node))
            for i=getNChildren(node):-1:1
                printRecursive(getChildren(node)[i], depth + 1)
            end
        end
    end
    println(getRepresentation(gramm, tree))
    for i=getNChildren(tree):-1:1
        printRecursive(getChildren(tree)[i], 1)
    end
    nothing
end

"""
    getPreOrder(tree::GERep, onlyNonTerminal::Bool=false)

Returns the given tree as an array in preorder.
"""
function getPreOrder(tree::GERep, onlyNonTerminal::Bool=false)

    stack = [tree]
    preOrder = Array{GERep}(undef, 0)

    while !isempty(stack)
        node = popfirst!(stack)
        if isLeaf(node)
            if !onlyNonTerminal
                push!(preOrder, node)
            end
        else
            push!(preOrder, node)
            children = getChildren(node)
            nChildren = length(children)
            for i=nChildren:-1:1
                pushfirst!(stack, children[i])
            end
        end
    end
    return preOrder
end

"""
    countNonTerminal(tree::GERep)

Returns the number of non terminals in the tree.
"""
function countNonTerminal(tree::GERep)
    stack = [tree]
    count = 0
    while !isempty(stack)
        node = popfirst!(stack)

        if !isLeaf(node)
            count = count + 1
            children = getChildren(node)
            nChildren = length(children)
            for i=nChildren:-1:1
                pushfirst!(stack, children[i])
            end
        end
    end
    return count
end

"""
    getRepresentation(gramm::Grammar, tree::GERep)

Interface for `getRepresentation(gramm::Grammar, id::UInt16)`
"""
function getRepresentation(gramm::Grammar, tree::GERep)
    getRepresentation(gramm, getSymbol(tree))
end

"Main struct for Grammatical Evolution genotypes used in GP."
struct GEGenotype <: GPGenotype
    _tree::GERep
    _ind::IntegerGenotype{UInt8}
    _selections::Array{UInt8}
    _productions::UInt16
    _valid::Bool
end # struct

"""
    pruneProtein(genotype1, genotype2)::Tuple{IntegerGenotype{UInt8}, Array{UInt8}, IntegerGenotype{UInt8}, Array{UInt8}}

Perfoms a prune dividing the codons in two, the useful part and the introns for each
genotype. Returning -> useful1, introns1, useful2, introns2
"""
function pruneProtein(genotype1, genotype2)::Tuple{IntegerGenotype{UInt8}, Array{UInt8}, IntegerGenotype{UInt8}, Array{UInt8}}
    protein1 = genotype1._ind
    len1 = genotype1._productions
    protein2 = genotype2._ind
    len2 = genotype2._productions

    maxlen = len1
    if len2 > len1
        maxLen = len2
    end
    subProtein1 = IntegerGenotype{UInt8}(protein1._representation[1:maxlen])
    introns1 = protein1._representation[maxlen+1:end]

    subProtein2 = IntegerGenotype{UInt8}(protein2._representation[1:maxlen])
    introns2 = protein2._representation[maxlen+1:end]

    return subProtein1, introns1, subProtein2, introns2
end

"""
    pruneProtein(genotype1)::Tuple{IntegerGenotype{UInt8}, Array{UInt8}}

Perfoms a prune dividing the codons in two, the useful part and the introns for
genotype. Returning -> useful, introns
"""
function pruneProtein(genotype)::Tuple{IntegerGenotype{UInt8}, Array{UInt8}}

    protein = genotype._ind
    len = genotype._productions

    subProtein = IntegerGenotype(protein._representation[1:len])
    introns = protein._representation[len:end]

    return subProtein, introns
end

"""
    copyGenotype(geGenotype::GEGenotype)::GEGenotype

Returns a copy of the genotype.
"""
function copyGenotype(geGenotype::GEGenotype)::GEGenotype
    tree = copyTree(geGenotype._tree)
    ind = copyGenotype(geGenotype._ind)
    selections = copy(geGenotype._selections)
    productions = geGenotype._productions
    valid = geGenotype._valid
    GEGenotype(tree, ind, selections, productions, valid)
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
abstract SuperLUMat

## Modify this for {T<:Union(Float64,Complex128)}
type NCMat <: SuperLUMat                # sparse matrix in compressed-column format
    sm::SuperMatrix
    colptr::Vector{int_t}
    rowind::Vector{int_t}
    nzval::Vector{Cdouble}
end

type DMat <: SuperLUMat                 # dense matrix
    sm::SuperMatrix
    nzval::Vector{Cdouble}
end


text:
using CameraCalibration, Test
@testset "general" begin
    @test true
end
"""
using CSV, DataFrames
imgpoints = Array{Float64, 2}[]
for i = 0:10
    push!(imgpoints, vcat(convert(Array{Float64, 2}, CSV.File(string("./test/imgpoints", i, ".csv"), header = false) |> DataFrame)', ones(1, 35)))
end
"""


text:
export AR_model
"""
    AR_model(df, y, grouping, d, link; penalized = penalized)

Form the autoregressive (AR(1)) model for intercept only regression with the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance parameter for the AR(1) structured covariance.
    One can put true (e.g. `penalized = true`) to add this penalty for numerical stability (default `penalized = false`).
"""
function AR_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaARObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = GaussianCopulaARObs(Y, X)
    end
    gcm = GaussianCopulaARModel(gcs; penalized = penalized);
    return gcm
end

function AR_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaARObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = GLMCopulaARObs(Y, X, d, link)
    end
    gcm = GLMCopulaARModel(gcs; penalized = penalized);
    return gcm
end

function AR_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaARObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        gcs[i] = NBCopulaARObs(Y, X, d, link)
    end
    gcm = NBCopulaARModel(gcs; penalized = penalized);
    return gcm
end

"""
    AR_model(df, y, grouping, covariates, d, link; penalized = penalized)

Form the autoregressive (AR(1)) model for regression with the specified base distribution (d) and link function (link).

# Arguments
- `df`: A named `DataFrame`
- `y`: Ouctcome variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `grouping`: Grouping or Clustering variable name of interest, specified as a `Symbol`.
    This variable name must be present in `df`.
- `covariates`: Covariate names of interest as a vector of `Symbol`s.
    Each variable name must be present in `df`.
- `d`: Base `Distribution` of outcome from `Distributions.jl`.
- `link`: Canonical `Link` function of the base distribution specified in `d`, from `GLM.jl`.

# Optional Arguments
- `penalized`: Boolean to specify whether or not to add an L2 Ridge penalty on the variance parameter for the AR(1) structured covariance.
    One can put true (e.g. `penalized = true`) to add this penalty for numerical stability (default `penalized = false`).
"""
function AR_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Normal, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GaussianCopulaARObs{Float64}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = GaussianCopulaARObs(Y, X)
    end
    gcm = GaussianCopulaARModel(gcs; penalized = penalized);
    return gcm
end

function AR_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:Union{Poisson, Bernoulli}, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{GLMCopulaARObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = GLMCopulaARObs(Y, X, d, link)
    end
    gcm = GLMCopulaARModel(gcs; penalized = penalized);
    return gcm
end

function AR_model(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol,
    covariates::Vector{Symbol},
    d::D,
    link::Link;
    penalized::Bool = false) where {D<:NegativeBinomial, Link}
    groups = unique(df[!, grouping])
    n = length(groups)
    gcs = Vector{NBCopulaARObs{Float64, D, Link}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        di = count(gidx)
        Y = Float64.(df[gidx, y])
        X = ones(di, 1)
        @inbounds for i in 1:length(covariates)
            U = Float64.(df[gidx, covariates[i]])
            X = hcat(X, U)
        end
        gcs[i] = NBCopulaARObs(Y, X, d, link)
    end
    gcm = NBCopulaARModel(gcs; penalized = penalized);
    return gcm
end


text:
import Plots

#Construct isaPlot3d from type AMFMmodel called with a time index using PlotsGR backend
function isaPlot3d_PlotsGR(z::AMFMmodel,t::Vector{Float64};
    a_max = maximum([maximum(abs.(z.comps[k].a.(t))) for k in 1:length(z.comps)]),
    ω_max = maximum([maximum(z.comps[k].ω.(t)) for k in 1:length(z.comps)]),
    ω_min = minimum(append!([minimum(z.comps[k].ω.(t)) for k in 1:length(z.comps)],[0])),
    s_max = maximum([maximum(abs.(real.(z.comps[k](t)))) for k in 1:length(z.comps)]),
    s_min = -s_max,
    t_max = minimum(t),
    t_min = maximum(t),
    FreqUnits = "rad/s",
    )
    if FreqUnits == "rad/s"
         Fnorm = 1
     elseif FreqUnits == "Hz"
         Fnorm = 1/2π
     else
         error("invalid FreqUnits")
     end
     Plots.plot3d(
        xlims = (t_max, t_min),
        ylims = (Fnorm*ω_min, Fnorm*ω_max),
        zlims = (s_min, s_max),
        legend = :false,
        framestyle = :origin,
        xlab = "time",
        ylab = "freq ("*FreqUnits*")",
        zlab = "real",
        camera = (20,80),
        background_color=ISA.cmap[1],
        foreground_color=:white,
        )
    for k in 1:length(z.comps)
        Plots.plot3d!(
            t,
            Fnorm.*z.comps[k].ω.(t),
            real.(z.comps[k](t)),
            c = ISA.cmap[ max.(min.(round.(Int, abs.(z.comps[k].a.(t)) .* 256/a_max ),256),1) ],
            linealpha = max.(min.( abs.(z.comps[k].a.(t)).^(1/2) .* 1/a_max ,1),0)
            )
    end
    Plots.current()
end

#Construct isaPlot3d from Array of Numerical AMFM Componets using PlotsGR backend
function isaPlot3d_PlotsGR(𝚿ₖ::Array{AMFMcompN,1};
    a_max = maximum([maximum(𝚿ₖ[k].a) for k in 1:length(𝚿ₖ)]),
    ω_max = maximum([maximum(𝚿ₖ[k].ω) for k in 1:length(𝚿ₖ)]),
    ω_min = minimum(append!([minimum(𝚿ₖ[k].ω) for k in 1:length(𝚿ₖ)],[0])),
    s_max = maximum([maximum(abs.(𝚿ₖ[k].s)) for k in 1:length(𝚿ₖ)]),
    s_min = -s_max,
    t_max = maximum([maximum(𝚿ₖ[k].t) for k in 1:length(𝚿ₖ)]),
    t_min = minimum([minimum(𝚿ₖ[k].t) for k in 1:length(𝚿ₖ)]),
    FreqUnits = "rad/s",
    )
    if FreqUnits == "rad/s"
         Fnorm = 1
     elseif FreqUnits == "Hz"
         Fnorm = 1/2π
     else
         error("invalid FreqUnits")
     end
    Plots.plot3d(
        xlims = (t_min, t_max),
        ylims = (Fnorm*ω_min, Fnorm*ω_max),
        zlims = (s_min, s_max),
        legend = :false,
        framestyle = :origin,
        xlab = "time",
        ylab = "freq ("*FreqUnits*")",
        zlab = "real",
        camera = (20,80),
        background_color=ISA.cmap[1],
        foreground_color=:white,
        )
    for k in 1:length(𝚿ₖ)
        Plots.plot3d!(
            𝚿ₖ[k].t,
            Fnorm.*𝚿ₖ[k].ω,
            𝚿ₖ[k].s,
            c = ISA.cmap[ max.(min.(round.(Int, 𝚿ₖ[k].a .* 256/a_max ),256),1) ],
            linealpha = max.(min.( 𝚿ₖ[k].a.^(1/2) .* 1/a_max ,1),0)
            )
    end
    Plots.current()
end

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
"""
A calendar with no holidays. But account for weekends.

`isholiday` always returns false for this calendar.
Remember that `isbday` considers that Saturdays and Sundars are not business days.
"""
struct WeekendsOnly <: HolidayCalendar end
isholiday(::WeekendsOnly, dt::Dates.Date) = false

"""
A calendar with no holidays and no weekends.
`bdays` returns the actual days between dates (`dt1 - d10`).
"""
struct NullHolidayCalendar <: HolidayCalendar end

isholiday(::NullHolidayCalendar, dt::Dates.Date) = false
isbday(::NullHolidayCalendar, dt::Dates.Date) = true
bdays(::NullHolidayCalendar, dt0::Dates.Date, dt1::Dates.Date) = dt1 - dt0

include("australia.jl")
include("brazil.jl")
include("canada.jl")
include("target.jl")
include("uk.jl")
include("us.jl")


text:
using Distributions
using OmegaCore
using Test

function testsolution()
  μ = 1 ~ Normal(0, 1)
  y = 2 ~ ((id, ω) -> Normal(μ(ω), 1)(id, ω))
  μc = μ |ᶜ (y ==ₚ 5.0)
  @test_throws OmegaCore.ConditionException μc(defω())
  sol = solution(μc)
end

function testlogpdfsol()
  μ = 1 ~ Normal(0, 1)
  y = 2 ~ ((id, ω) -> Normal(μ(ω), 1)(id, ω))
  μc = μ |ᶜ (y ==ₚ 5.0)
  ω = solution(μc)
  μ_ = μ(ω)
  @test logpdf(ω) == logpdf(Normal(0, 1), μ_) + logpdf(Normal(μ_, 1), 5.0)
end

function test_solution_intervene()
  μ = 1 ~ Normal(0, 1)
  y = 2 ~ ((id, ω) -> Normal(μ(ω), 1)(id, ω))
  μc = μ |ᶜ (y ==ₚ 5.0)
  yi = y |ᵈ (μ => (ω -> 100.0))
  yi2 = y |ᵈ (μ => (ω -> 100))
  ω = solution(μc)
  a, b = (yi(ω), y(ω))
  @test a != b
end

@testset "Solution" begin
  testsolution()
  testlogpdfsol()
  test_solution_intervene()
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
import DataFrames.DataFrame

@testset "fitter.jl" begin

@testset "fitter.jl exceptions" begin
    f_test(x, a, b) = a .* x .+ b
    g_test(x) = x.^2

    bad_data_test = DataFrame(Any[[1:10], [7:16]])

    fitter_test = Fitter(f_test; autoplot=false)

    @test_throws EasyPhys.CannotFitException  Fitter(g_test)

    @test_throws EasyPhys.NoResultsException  fitter_test |> apply_f(5)
    @test_throws EasyPhys.NoResultsException  fitter_test[:a]

    @test_throws EasyPhys.BadDataException    fit!(fitter_test)
    @test_throws EasyPhys.BadDataException    set_data!(
                                                fitter_test, [1, 2, 3],
                                                [4, 5], [1, 2, 3])
    @test_throws EasyPhys.BadDataException    set_data!(
                                                fitter_test, [1, 2, 3],
                                                [4, 5, 6], [1, 2])
    @test_throws EasyPhys.BadDataException    set_data!(
                                                fitter_test, bad_data_test)

    @test_throws KeyError fitter_test[:notakey] = 7
    @test_throws KeyError something = fitter_test[:stillnotakey]

end

tolerance_test = 0.08
outlier_threshold_test = 2

model_test(x, a, b) = a .* exp.(-x .* b)

xdata_test = linspace(0, 10, 100)
eydata_test = 0.01

@testset "fitter.jl `fit!` and friends" begin
    fitter_test = Fitter(model_test) |> set!(autoplot=false)
    show(fitter_test)

    for a in 1:20:81, b in 1:2:11
        ydata_test = model_test(xdata_test, a, b) + 0.01*randn(length(xdata_test))

        set_data!(fitter_test, xdata_test, ydata_test, eydata_test) |> fit!
        χ²_init_test = reduced_χ²(fitter_test, [a, b])
        χ²_worse_test = reduced_χ²(fitter_test)

        @test all(abs.([fitter_test[:a][1], fitter_test[:b][1]] .- [a, b])
                    ./ [a, b]
                    .<= tolerance_test)

        ignore_outliers!(fitter_test, outlier_threshold_test) |> fit!
        χ²_better_test = reduced_χ²(fitter_test)

        @test χ²_init_test >= χ²_worse_test >= χ²_worse_test

        y_true_test = model_test(xdata_test, a, b)
        @test all(apply_f(fitter_test, xdata_test, [a, b]) == y_true_test)

        fix!(fitter_test; a=1.01a) |> fit! |> apply_f(5)
        free!(fitter_test, :a)
    end

    fitter_test[:xmax] = 90
    fit!(fitter_test)
    show(fitter_test)

    show(parameter_covariance(fitter_test))
    show("\n\n")
end

@testset "fitter.jl fixing and freeing parameters" begin
    fitter_test = Fitter(model_test) |> set!(autoplot=false)

    a = 1
    b = 2
    ydata_test = model_test(xdata_test, a, b) + 0.01*randn(length(xdata_test))
    set_data!(fitter_test, xdata_test, ydata_test, eydata_test)

    fitter_test |> fix!(a=1) |> guess!(b=2) |> fit!
    @test fitter_test[:a] == 1
    show(fitter_test)

    guess!(fitter_test; a=3)
    @test fitter_test[:a] == 1

    fitter_test[:a] = 1.1
    fitter_test[:b] = 3
    fitter_test |> fit!

    fix!(fitter_test; b=2)
    @test_throws EasyPhys.CannotFitException fit!(fitter_test)

    fitter_test |> free!(a=1) |> fit!
    fitter_test |> fix!(b=1) |> guess!([2.0]) |> fit!
    show(fitter_test)

    fitter_test |> free!(:b) |> guess!([2.0, 1.0]) |> fit!
    show(fitter_test)
end

end


text:
using CSV, XLSX
using DataFrames

function load_dir(dir)
    res = []
    for file in readdir(dir)
		df = CSV.read(joinpath(dir, file), DataFrame)
		insertcols!(df, 1, :file=>[file])
        push!(res, df)
    end
    res
end

gap_compare_results = [("results/GAP_exp/ssit", "GAP_comp_ssit.xlsx"), ("results/GAP_exp/norm", "GAP_comp_norm.xlsx")]
gurobi_gap_comp = [
	("results/GAP_exp_Gurobi/first", "gurobi_first.xlsx"),
	("results/GAP_exp_Gurobi/second", "gurobi_second.xlsx"),
	("results/GAP_exp_Gurobi/third", "gurobi_third.xlsx")
]
mmkp_res = [
	("results/long_mmkp/ssit", "long_ssit_mmkp.xlsx"),
	("results/long_mmkp/base", "long_base_mmkp.xlsx")
	]

function write_results(results_pairs)
	for (d, n) in results_pairs
		r = load_dir(d)
		df = vcat(r..., cols=:union)
		XLSX.writetable(
				n,
				results = (
					collect(DataFrames.eachcol(df)),
					DataFrames.names(df)),
				overwrite=true)
	end
end

write_results(mmkp_res)

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
module OptimTestProblems

include("optim_tests/multivariate/unconstrained.jl")
include("optim_tests/univariate/bounded.jl")

end # module


text:
using Clapeyron, Test

@printline

struct TestModel <: EoSModel end

function Clapeyron.eos(model::TestModel,V,T,z)
    z_part = 1*z[1] + 2*z[2]+3*z[3]
    v_part = log(V)
    T_part = T^6
    vt_part = V*T*V*T
    return vt_part+T_part+v_part+z_part
end

∂f∂T_analytical(model::TestModel,V,T,z) =6*T^5 +2*T*V*V
∂f∂V_analytical(model::TestModel,V,T,z) =1/V +2*V*T*T
∂2f∂T2_analytical(model::TestModel,V,T,z) =5*6*T^4 +2*V*V
∂2f∂V2_analytical(model::TestModel,V,T,z) =-1/(V*V) +2*T*T
∂2f∂V∂T_analytical(model::TestModel,V,T,z) =4*V*T
∂3f∂V3_analytical(model::TestModel,V,T,z) =2/(V*V*V)
∂3f∂V2∂T_analytical(model::TestModel,V,T,z) = 4*T
∂3f∂V∂T2_analytical(model::TestModel,V,T,z) = 4*V

@testset "differentials" begin
    model = TestModel()
    T = 500*rand()
    V = 10*rand()
    z = rand(3)
    f = Clapeyron.eos(model,V,T,z)
    v = ∂f∂V_analytical(model,V,T,z)
    t = ∂f∂T_analytical(model,V,T,z)
    vv = ∂2f∂V2_analytical(model,V,T,z)
    vt = ∂2f∂V∂T_analytical(model,V,T,z)
    tt = ∂2f∂T2_analytical(model,V,T,z)
    vvv = ∂3f∂V3_analytical(model,V,T,z)

    df = [v,t]
    d2f = [vv vt;vt tt]
    p = -v
    pv = -vv
    pvv = -vvv
    pt = -vt
    dp = [pv,pt]
    pvt = -∂3f∂V2∂T_analytical(model,V,T,z)
    ptt = -∂3f∂V∂T2_analytical(model,V,T,z)
    d2p = [pvv pvt;pvt ptt]

    @testset "first order" begin
        @test Clapeyron.∂f∂T(model,V,T,z) ≈ t
        @test Clapeyron.∂f∂V(model,V,T,z) ≈ v
        df1 = Clapeyron.∂f(model,V,T,z)
        @test df1[2] ≈ f
        @test df1[1][1] ≈  v
        @test df1[1][2] ≈  t
    end

    @testset "second order" begin
        pdp = Clapeyron.p∂p∂V(model,V,T,z)
        @test pdp[1] ≈ -v
        @test pdp[2] ≈ -vv
        h = Clapeyron.f_hess(model,V,T,z)
        @test all(h .≈ d2f)
        ddf = Clapeyron.∂2f(model,V,T,z)
        @test ddf[3] ≈ f
        @test all(ddf[2] .≈ df)
        @test all(ddf[1] .≈ d2f)
    end

    @testset "third order" begin
        dddf = Clapeyron.∂²³f(model,V,T,z)
        @test dddf[1] ≈ -vv
        @test dddf[2] ≈ -vvv

        ddp = Clapeyron.∂2p(model,V,T,z)
        @test ddp[3] ≈ p
        @test all(ddp[2] .≈ dp)
        @test all(ddp[1] .≈ d2p)
    end
end


text:
function simulation(
    seed_number;
    moving = false,
    vaccin = false)

    folder = (moving ? 'M' : '0') * (vaccin ? 'V' : '0')
    Random.seed!(seed_number);
    if folder == "00"
        print(Crayon(foreground = :white), seed_number)
    elseif folder == "0V"
        print(Crayon(foreground = :light_red), seed_number)
    elseif folder == "M0"
        print(Crayon(foreground = :light_blue), seed_number)
    elseif folder == "MV"
        print(Crayon(foreground = :magenta), seed_number)
    end

    ####################################################################

    n_S_ = Int64[]
    n_E_ = Int64[]
    n_I_ = Int64[]
    n_R_ = Int64[]
    n_V_ = Int64[]
    # n_hub_ = Int64[]

    contact_ = Int64[]
    RT_ = Float64[]
    SI_ = Int64[]
    # entropy_ = Float64[]
    # var_ = Float64[]

    n_moved_S_ = Int64[]
    n_moved_E_ = Int64[]
    n_moved_I_ = Int64[]
    n_moved_V_ = Int64[]
    n_moved_R_ = Int64[]

    # S_influx_ = Int64[]
    # I_influx_ = Int64[]
    # E_influx_ = Int64[]
    # R_influx_ = Int64[]
    # V_influx_ = Int64[]

    # S_outflux_ = Int64[]
    # I_outflux_ = Int64[]
    # E_outflux_ = Int64[]
    # R_outflux_ = Int64[]
    # V_outflux_ = Int64[]

    # n_NODE_S_ = DataFrame(zeros(Int64, end_time, N), :auto)
    # n_NODE_I_ = DataFrame(zeros(Int64, end_time, N), :auto)
    # n_NODE_incd_ = DataFrame(zeros(Int64, end_time, N), :auto)

    transmission = DataFrame(
        T = Int64[],
        t = Int64[],
        node_id = [],
        from = [],
        to = [],
        Break = Int64[]
    )
    non_transmission = copy(transmission)

    ####################################################################

    T = 0 # Macro timestep
    state = repeat(['S'], n) # using SEIR model
    host = rand(ID, number_of_host); state[host] .= 'I'

    LATENT = repeat([-1], n)
    RECOVERY = repeat([-1], n)
    RECOVERY[host] .= round.(rand(recovery_period, number_of_host)) .+ 1

    coordinate = rand(Float16, 2, n) # micro location
    LOCATION = rand(NODE_ID, n) # macro location
    for _ in 1:5 LOCATION = rand.(getindex.(Ref(NODE),LOCATION)) end

while sum(state .== 'E') + sum(state .== 'I') > 0
    if T ≥ end_time
        break
    else
        T += 1
    end
    # if T == 38 Random.seed!(seed_number) end

    LATENT   .-= 1
    RECOVERY .-= 1
    bit_LATENT     = (LATENT   .== 0)
    bit_RECOVERY   = (RECOVERY .== 0)
    state[bit_LATENT  ] .= 'I'
    state[bit_RECOVERY] .= 'R'

    bit_S = (state .== 'S'); n_S = count(bit_S); push!(n_S_, n_S)
    bit_E = (state .== 'E'); n_E = count(bit_E); push!(n_E_, n_E)
    bit_I = (state .== 'I'); n_I = count(bit_I); push!(n_I_, n_I)
    bit_R = (state .== 'R'); n_R = count(bit_R); push!(n_R_, n_R)
    bit_V = (state .== 'V'); n_V = count(bit_V); push!(n_V_, n_V)
    # bit_hub = ((LOCATION .== 5) .| (LOCATION .== 4)) ; n_hub = count(bit_hub); push!(n_hub_, n_hub)

    # if T > 0
    #     println("$T: |E: $(n_E_[T]) |I: $(n_I_[T]) |R:$(n_R_[T]) |V:$(n_V_[T])")
    # end

    mobility = σ
    if n_I > θ
        if vaccin state[bit_S .& (rand(n) .< p_V)] .= 'V' end
        if moving mobility = σ*e_σ end
    end
    moved = (rand(n) .< mobility)
    push!(n_moved_S_, count(moved .& bit_S))
    push!(n_moved_E_, count(moved .& bit_E))
    push!(n_moved_I_, count(moved .& bit_I))
    push!(n_moved_R_, count(moved .& bit_R))
    push!(n_moved_V_, count(moved .& bit_V))

    # bit_in = (moved .& bit_hub)
    # push!(S_outflux_, count(bit_in .& bit_S))
    # push!(E_outflux_, count(bit_in .& bit_E))
    # push!(I_outflux_, count(bit_in .& bit_I))
    # push!(R_outflux_, count(bit_in .& bit_R))
    # push!(V_outflux_, count(bit_in .& bit_V))

    # LOCATION[moved] = rand.(NODE[LOCATION[moved]])
    LOCATION[moved] = rand.(getindex.(Ref(NODE),LOCATION[moved]))
    # bit_out = (moved .& .!(bit_hub)) .& ((LOCATION .== 5) .| (LOCATION .== 4))
    # push!(S_influx_, count(bit_out .& bit_S))
    # push!(E_influx_, count(bit_out .& bit_E))
    # push!(I_influx_, count(bit_out .& bit_I))
    # push!(R_influx_, count(bit_out .& bit_R))
    # push!(V_influx_, count(bit_out .& bit_V))

    contact_t = 0
    SI_t = 0
    NODE_I = unique(LOCATION[bit_I])
    for node in NODE_I
        bit_node = (LOCATION .== node); n_micro = count(bit_node)

        bit_micro_S = bit_node .& bit_S
        # bit_micro_E = bit_node .& bit_E
        bit_micro_I = bit_node .& bit_I
        bit_micro_V = bit_node .& bit_V

        ID_S = ID[bit_micro_S]; n_micro_S = count(bit_micro_S)
        # ID_E = ID[bit_micro_E]; n_micro_E = count(bit_micro_E)
        ID_I = ID[bit_micro_I]; n_micro_I = count(bit_micro_I)
        ID_V = ID[bit_micro_V]; n_micro_V = count(bit_micro_V)

        SI_t += n_micro_S*n_micro_I

        # if seed_number != 0 n_NODE_S_[T, node] = n_micro_S end
        # if seed_number != 0 n_NODE_I_[T, node] = n_micro_I end
        for t in 1:12
            coordinate[:,bit_node] = mod.(coordinate[:,bit_node] + rand(brownian, n_micro), 1.0)
            kdtreeI = KDTree(coordinate[:,ID_I])

            if (n_micro_S == 0) continue end # transmission I to S
            in_δ = inrange(kdtreeI, coordinate[:,ID_S], δ)
            contact = length.(in_δ)
            contact_t += sum(contact)

            bit_infected = (rand(n_micro_S) .< (1 .- (1 - β).^contact))
            ID_infected = ID_S[bit_infected]

            n_infected = count(bit_infected)
            if n_infected > 0
                from_id = ID_I[deep_pop!.(shuffle.(in_δ))[bit_infected]]
                append!(transmission, DataFrame(
                    T = T, t = t, node_id = node,
                    from = from_id, to = ID_infected, Break = 0
                    ))
            end
            state[ID_infected] .= 'E'
            LATENT[ID_infected] .= round.(rand(latent_period, n_infected))
            RECOVERY[ID_infected] .= LATENT[ID_infected] + round.(rand(recovery_period, n_infected))

            if !(e_V > 0) continue end # for efficient simulation
            if (n_micro_V == 0) continue end # transmission I to V
            in_δ = inrange(kdtreeI, coordinate[:,ID_V], δ)
            contact = length.(in_δ)
            contact_t += sum(contact)

            bit_infected = (rand(n_micro_V) .< (1 .- (1 - β*e_V).^contact))
            ID_infected = ID_V[bit_infected]

            n_infected = count(bit_infected)
            if n_infected > 0
                from_id = ID_I[deep_pop!.(shuffle.(in_δ))[bit_infected]]
                append!(transmission, DataFrame(
                    T = T, t = t, node_id = node,
                    from = from_id, to = ID_infected, Break = 1
                    ))
            end
            state[ID_infected] .= 'E'
            LATENT[ID_infected] .= round.(rand(latent_period, n_infected))
            RECOVERY[ID_infected] .= LATENT[ID_infected] + round.(rand(recovery_period, n_infected))
        end
        unique!(transmission, :to)

        infector_list = transmission[transmission.T .== T, :from]
        for defendant ∈ ID_I
            if defendant ∉ infector_list
                append!(non_transmission, DataFrame(
                    T = T, t = 0, node_id = node, #x = 0, y = 0,
                    from = defendant, to = 0, Break = 0
                    ))
            end
        end
    end
    push!(contact_, contact_t)
    push!(SI_, SI_t)

    # transmission_node_id = transmission.node_id
    # count_node_incidence = [count(transmission_node_id .== node_id) for node_id ∈ 1:N]
    # n_NODE_incd_[T,:] = count_node_incidence
    # entropy_ = push!(entropy_, entropy(count_node_incidence ./ sum(count_node_incidence), N))
    # var_ = push!(var_, var(count_node_incidence))

    now_I = ID[state .== 'I']
    push!(
        RT_,
        isempty(now_I) ? 0 : nrow(transmission[transmission[:,:from] .∈ Ref(now_I),:]) / length(now_I))
end

print(": $(n_R_[T]), ")
T1 = sufficientN(RT_ .< 1)
peaktime = argmax(n_I_)
peaksize = maximum(n_I_)
# if T == end_time
# if true
seed = lpad(seed_number, 4, '0')

fromby = groupby(transmission, :from)
agent_id = Int64[]
node_id = []
home = Int64[]
away = Int64[]
for agent ∈ fromby
    node_infected = agent.node_id[1]
    push!(agent_id, agent.from[1])
    push!(node_id, node_infected)
    push!(home, count(agent.node_id .== node_infected))
    push!(away, count(agent.node_id .!= node_infected))
end
case = home .+ away
agent_seconarycases = DataFrame(; agent_id, node_id, case, home, away)

# node_incidence = transmission[transmission.to .> 0,:node_id]
# count_node_incidence = [count(node_incidence .== node_id) for node_id ∈ 1:N]
# incidence5 = count_node_incidence[5]; incidence4 = count_node_incidence[4]
# hub_incidence_rate = (incidence5 + incidence4) / sum(count_node_incidence)
# incidence_entropy = entropy(count_node_incidence ./ sum(count_node_incidence), N)
# incidence_var = var(count_node_incidence)

if seed_number != 0 append!(transmission, non_transmission); sort!(transmission, :T) end
time_evolution = DataFrame(;
n_S_, n_E_, n_I_, n_R_, n_V_, contact_, RT_,
# entropy_, var_, SI_,
# n_hub_,
# ,n_moved_S_, n_moved_E_, n_moved_I_, n_moved_R_, n_moved_V_
# ,S_influx_, I_influx_, E_influx_, R_influx_, V_influx_
# ,S_outflux_, I_outflux_, E_outflux_, R_outflux_, V_outflux_
)

# count_node_incidence = [count(node_incidence .== node_id) for node_id ∈ 1:N]
# incidence5 = count_node_incidence[5]
# incidence4 = count_node_incidence[4]
# hub_incidence_rate = (incidence5 + incidence4) / sum(count_node_incidence)
# incidence_entropy = entropy(count_node_incidence ./ sum(count_node_incidence), N)

summary = DataFrame(Tend = T, Rend = n_R_[T], Vend = n_V_[T], peaktime = peaktime, peaksize = peaksize,
    # incidence5 = incidence5, incidence4 = incidence4, HIR = hub_incidence_rate,
    # incidence_entropy = entropy_[T], incidence_var = var_[T],
    T1 = T1, RT1 = n_R_[T1], VT1 = n_V_[T1],
    home = sum(home), away = sum(away))

# if export_type != :XLSX
# CSV.write("./$folder/$seed ndws.csv", n_NODE_S_)
# CSV.write("./$folder/$seed ndwi.csv", n_NODE_I_)
# CSV.write("./$folder/$seed incd.csv", n_NODE_incd)
# CSV.write("./$folder/$seed ndwf.csv", n_NODE_S_ .* n_NODE_I_) # Force of infection
CSV.write("./$folder/$seed trms.csv", transmission)
CSV.write("./$folder/$seed agnt.csv", agent_seconarycases)
CSV.write("./$folder/$seed tevl.csv", time_evolution)
CSV.write("./$folder/$seed smry.csv", summary, bom = true)
# elseif export_type != :CSV
#     XLSX.writetable(
#         "./$folder/$seed.xlsx",
#         time_evolution = ( collect(eachcol(time_evolution)), names(time_evolution) ),
#         transmission = ( collect(eachcol(transmission)), names(transmission) ),
#         config = ( collect(eachcol(config)), names(config) ),
#         summary = ( collect(eachcol(summary)), names(summary) )
#     )
# end
# end
end

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
# Unit tests for judiWeights
# Rafael Orozco (rorozco@gatech.edu)
# July 2021
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

# number of sources/receivers
nsrc = 1
nrec = 120
nx = 4
ny = 4
nt = 10
dt = 2f0
ftol = 1f-6

################################################# test constructors ####################################################

@testset "judiWavefield Unit Tests with $(nsrc) sources" for nsrc=[1, 2]

    # set up judiWeights
    info = Info(nx*ny, nsrc, nt)

    # Extended source weights
    wf = Array{Array{Float32, 3}, 1}(undef, nsrc)
    for j=1:nsrc
        wf[j] = randn(Float32, nt, ny, ny)
    end
    w = judiWavefield(info, dt, wf)
    w1 = similar(w) .+ 1f0

    @test isequal(length(w.data), nsrc)
    @test isequal(length(w.data), info.nsrc)
    @test isequal(w.info.nsrc, nsrc)
    @test isequal(typeof(w.data), Array{Array{Float32, 3}, 1})
    @test isequal(size(w), (nx*ny*nt*nsrc, 1))
    @test isfinite(w)

#################################################### test operations ###################################################

    # conj, transpose, adjoint
    @test isequal(size(w), size(conj(w)))
    @test isequal(reverse(size(w)), size(transpose(w)))
    @test isequal(reverse(size(w)), size(adjoint(w)))

    # +, -, *, /
    @test iszero(norm(2*w - (w + w)))
    @test iszero(norm(w - (w + w)/2))
    @test iszero(norm(1f0 - w1))
    @test isequal(norm(1f0 + w1, 1), 2f0 * norm(w1, 1))
    #test unary operator
    @test iszero(norm((-w) - (-1*w)))

    # vcat
    w_vcat = [w; w]
    @test isequal(length(w_vcat), 2*length(w))
    @test isequal(w_vcat.info.nsrc, 2*nsrc)
    @test isequal(length(w_vcat.data), 2*nsrc)

    # dot, norm, abs
    @test isapprox(norm(w), sqrt(dot(w, w)))
    @test isapprox(abs.(w.data[1]), abs(w).data[1])

    # Test the norm
    d_ones = judiWavefield(info, dt, 2f0 .* ones(Float32, nt, nx, ny))
    @test isapprox(norm(d_ones, 2), sqrt(dt*nt*nx*ny*4*nsrc))
    @test isapprox(norm(d_ones, 1), dt*nt*nx*ny*2*nsrc)
    @test isapprox(norm(d_ones, Inf), 2)

    # vector space axioms
    u = judiWavefield(info, dt, randn(Float32, nt, nx, ny))
    v = judiWavefield(info, dt, randn(Float32, nt, nx, ny))
    w = judiWavefield(info, dt, randn(Float32, nt, nx, ny))
    a = .5f0 + rand(Float32)
    b = .5f0 + rand(Float32)

    @test isapprox(u + (v + w), (u + v) + w; rtol=ftol)
    @test isapprox(u + v, v + u; rtol=ftol)
    @test isapprox(u, u + 0; rtol=ftol)
    @test iszero(norm(u + u*(-1)))
    @test isapprox(-u, (-1f0)*u; rtol=ftol)
    @test isapprox(a .* (b .* u), (a * b) .* u; rtol=ftol)
    @test isapprox(u, u .* 1; rtol=ftol)
    @test isapprox(a .* (u + v), a .* u + a .* v; rtol=1f-5)
    @test isapprox((a + b) .* v, a .* v + b.* v; rtol=1f-5)

    # test fft
    fw = fft_wavefield(w, 1)
    fwf = fft_wavefield(fw, -1)
    @test isapprox(dot(fwf, w), real(dot(fw, fw)); rtol=ftol)
    @test isapprox(fwf, w; rtol=ftol)

end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
#Gravity related stuff

#Additional source terms
function source_terms(hyd::data2d, dt)

    #self-gravity
    #gxflx, gyflx = selfgravity(hyd.rho, hyd.nx, hyd.ny, r32x, r32y)

    #constant gravity

    #x-dir
    dir = 1

    #+
    gxflx, gyflx = ygravity(hyd.rhop[:,:,dir], hyd.nx, hyd.ny)
    hyd.velxp[:,:,dir] += 0.5*gxflx
    hyd.velyp[:,:,dir] += 0.5*gyflx
    hyd.epsp[:,:,dir] += 0.5*((hyd.velxp[:,:,dir] .* gxflx) .+ (hyd.velyp[:,:,dir] .* gyflx))

    #-
    gxflx, gyflx = ygravity(hyd.rhom[:,:,dir], hyd.nx, hyd.ny)
    hyd.velxm[:,:,dir] += 0.5*gxflx
    hyd.velym[:,:,dir] += 0.5*gyflx
    hyd.epsm[:,:,dir] += 0.5*((hyd.velxm[:,:,dir] .* gxflx) .+ (hyd.velym[:,:,dir] .* gyflx))

    #y-dir
    dir = 2

    #+
    gxflx, gyflx = ygravity(hyd.rhop[:,:,dir], hyd.nx, hyd.ny)
    hyd.velxp[:,:,dir] += 0.5*gxflx
    hyd.velyp[:,:,dir] += 0.5*gyflx
    hyd.epsp[:,:,dir] += 0.5*((hyd.velxp[:,:,dir] .* gxflx) .+ (hyd.velyp[:,:,dir] .* gyflx))

    #-
    gxflx, gyflx = ygravity(hyd.rhom[:,:,dir], hyd.nx, hyd.ny)
    hyd.velxm[:,:,dir] += 0.5*gxflx
    hyd.velym[:,:,dir] += 0.5*gyflx
    hyd.epsm[:,:,dir] += 0.5*((hyd.velxm[:,:,dir] .* gxflx) .+ (hyd.velym[:,:,dir] .* gyflx))

    return hyd
end

#Constant y-dir gravity
function ygravity(rho::AbstractMatrix,
                  nx::Int64,
                  ny::Int64;
                  gx=0.0,
                  gy=-1.0)

    gxflux = gx .* rho
    gyflux = gy .* rho

    return gxflux, gyflux
end

#Calculate x(x^2 + y^2)^-3/2 and y(x^2 + y^2)^-3/2 kernel
#for convolution used in self-gravity
function r32kernel(x, y)

    Nx = length(x)
    Ny = length(y)
    dx = x[2] - x[1]
    dy = y[2] - y[1]

    r32x = zeros(2Ny-1, 2Nx-1)
    r32y = zeros(2Ny-1, 2Nx-1)

    r32q = zeros(Ny, Nx)

    #construct 1 quarter of the kernel
    for j = 0:(Ny-1), i = 0:(Nx-1)
        r32q[j+1, i+1] = ((i*dx)^2.0 + (j*dy)^2.0)^(-1.5)
    end
    r32q[1,1] = 0.0

    #x - axis component
    r32x[1:Ny, Nx:2Nx-1] = flipud(r32q)
    r32x[Ny:2Ny-1, Nx:2Nx-1] = r32q

    #multiply by x
    for j = 1:2Ny-1, i = Nx:2Nx-1
        r32x[j, i] = r32x[j, i]*dx*(i-Nx)
    end

    #flip negative values
    r32x[:, 1:Nx-1] = -1.0*fliplr(r32x[:, Nx+1:end])

    #y - axis component
    r32y[Ny:2Ny-1, 1:Nx] = fliplr(r32q)
    r32y[Ny:2Ny-1, Nx:2Nx-1] = r32q

    #multiply by y
    for j = Ny:2Ny-1, i = 1:2Nx-1
        r32y[j, i] = r32y[j, i]*dy*(j-Ny)
    end

    #flip negative values
    r32y[1:Ny-1, :] = -1.0*flipud(r32y[Ny+1:end, :])

    return r32x, r32y
end

#Self-gravity
function selfgravity(rho::AbstractMatrix,
                     nx::Int64,
                     ny::Int64,
                     r32x::AbstractMatrix,
                     r32y::AbstractMatrix)

    xr = nx:(2nx-1)
    yr = ny:(2ny-1)

    gxflux = conv2(rho, r32x)[xr, yr]
    gyflux = conv2(rho, r32y)[xr, yr]

    gxflux = gxflux .* rho
    gyflux = gyflux .* rho

    return gxflux, gyflux
end

max_stars_repo_name: gzahl/hydro

max_issues_repo_name: gzahl/hydro

max_forks_repo_name: gzahl/hydro


text:
using Documenter, SeisNoiseAWS

makedocs(;
    modules=[SeisNoiseAWS],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/dylanmikesell/SeisNoiseAWS.jl/blob/{commit}{path}#L{line}",
    sitename="SeisNoiseAWS.jl",
    authors="Dylan Mikesell, Boise State University",
    assets=String[],
)

deploydocs(;
    repo="github.com/dylanmikesell/SeisNoiseAWS.jl",
)


text:
module TestParseWDI

using Base.Test
using WorldBankDataTd
using Dates

us_gnp_data = {  { "total"=>23,"per_page"=>"25000","pages"=>1,"page"=>1 },
                 { ["date"=>"2012","value"=>"52340","indicator"=>["id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"],"country"=>["id"=>"US","value"=>"United States"],"decimal"=>"0"],
                   ["date"=>"2011","value"=>"50650","indicator"=>["id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"],"country"=>["id"=>"US","value"=>"United States"],"decimal"=>"0"],
                   ["date"=>"2010","value"=>"48960","indicator"=>["id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"],"country"=>["id"=>"US","value"=>"United States"],"decimal"=>"0"],
                   ["date"=>"2009","value"=>"48040","indicator"=>["id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"],"country"=>["id"=>"US","value"=>"United States"],"decimal"=>"0"],
                   ["date"=>"2008","value"=>"49350","indicator"=>["id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"],"country"=>["id"=>"US","value"=>"United States"],"decimal"=>"0"],
                   ["date"=>"2007","value"=>"48640","indicator"=>["id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"],"country"=>["id"=>"US","value"=>"United States"],"decimal"=>"0"],
                   ["date"=>"2006","value"=>"48080","indicator"=>["id"=>"NY.GNP.PCAP.CD","value"=>"GNI per capita, Atlas method (current US\$)"],"country"=>["id"=>"US","value"=>"United States"],"decimal"=>"0"]
                 }
              }

us_gnp = WorldBankDataTd.parse_wdi("NY.GNP.PCAP.CD",us_gnp_data[2],2006,2012)

@test us_gnp[:time] == Date[Date(2012,12,31), Date(2011,12,31),
                            Date(2010,12,31), Date(2009,12,31),
                            Date(2008,12,31), Date(2007,12,31),
                            Date(2006,12,31)]

@test us_gnp[symbol("NY.GNP.PCAP.CD")] == Float64[52340, 50650, 48960, 48040, 49350, 48640, 48080]

end


text:
# include this file @everywhere for local pic processing
using Images #first include Images to avoid race conditions
using LeafAreaIndex

function listrawfiles(dir::AbstractString; raw_ext=LeafAreaIndex.RAW_EXT)
    allfiles = readdir(dir)
    filenames = ASCIIString[]
    filepaths = ASCIIString[]

    for file in allfiles
        if lowercase(splitext(file)[2]) ∈ raw_ext
            push!(filenames, file)
            push!(filepaths, joinpath(dir, file))
        end
    end
    filenames, filepaths
end

function LeafAreaIndex.PolarImage(file::AbstractString, cl::CameraLens)
    img = rawblueread(file)
    size(img) == size(cl) || error("wrong image dimensions")
    LeafAreaIndex.isoverexposed(img) && warn("Overexposed image $file")

    if Images.height(img) > Images.width(img) # Rotate portrait mode.
        img["spatialorder"] = reverse(img["spatialorder"])
    end

    PolarImage(img, cl)
end

# Always use EllipsOpt and langxiang45 methods
type LAIresult{T<:LeafAreaIndex.ThresholdMethod}
	file::ASCIIString
    LAI::Float64
    thresh::Float64
    thresh_method::T
    clump::Float64
    overexp::Bool
end

function getLAI(file::AbstractString, cl::CameraLens)
    polim = PolarImage(file, cl)
    binars = [RidlerCalvard(), EdgeDetection(), MinimumThreshold()]

    oe = LeafAreaIndex.isoverexposed(polim)

    lais = LAIresult[]
    for binar in binars
    	thresh = threshold(polim, binar)
    	Ω = langxiang45(polim, thresh, 0, pi/2)
    	LAIe = inverse(polim, thresh)
    	filename = last(splitdir(file))
    	push!(lais, LAIresult(filename, LAIe / Ω, Float64(thresh), binar, Ω, oe))
	end
	lais
end


text:
    # Functions and methods for cascade computation

    """
    `Cascade.computeSteps(scheme::Cascade.StepwiseDecayScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})`
        ... computes in turn all the requested transition amplitudes as well as PhotoEmission.Line's, AutoIonization.Line's,
            etc. for all pre-specified decay steps of the cascade. When compared with standard computations of these atomic
            processes, however, the amount of output is largely reduced and often just printed into the summary file.
            A set of  data::Cascade.DecayData  is returned.
    """
    function computeSteps(scheme::Cascade.StepwiseDecayScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})
        linesA = AutoIonization.Line[];    linesR = PhotoEmission.Line[];    cOrbitals = Dict{Subshell, Orbital}()
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        nt = 0;   st = 0;   previousMeanEn = 0.
        for  step  in  stepList
            st = st + 1
            nc = length(step.initialMultiplet.levels) * length(step.finalMultiplet.levels)
            println("\n  $st) Perform $(string(step.process)) amplitude computations for up to $nc decay lines (without selection rules): ")
            if  printSummary   println(iostream, "\n* $st) Perform $(string(step.process)) amplitude computations for " *
                                                 "up to $nc decay lines (without selection rules): ")   end

            if      step.process == Basics.Auger()  &&   comp.approach == Cascade.AverageSCA()
                # First determine the `mean' free-electron energy for this Auger block and calculate a common set of continuum orbital
                meanEn = 0.;    NoEn = 0
                for  p = 1:length(step.initialMultiplet.levels),  q = 1:length(step.finalMultiplet.levels)
                    en = step.initialMultiplet.levels[p].energy - step.finalMultiplet.levels[q].energy
                    if  en > 0.1    meanEn = meanEn + en;    NoEn = NoEn + 1    end
                end
                if  NoEn > 0     meanEn = meanEn/ NoEn     else     meanEn = 0.1    end
                if  abs(meanEn - previousMeanEn) / meanEn  < 0.15   # no new continuum orbitals
                    println(">> No new continum orbitals are generated for $(keys(cOrbitals)) and for the energy $meanEn ")
                else
                    previousMeanEn = meanEn;    cOrbitals = Dict{Subshell, Orbital}()
                    for kappa = -step.settings.maxKappa-1:step.settings.maxKappa        if   kappa == 0     continue    end
                        sh           = Subshell(101, kappa);     nrContinuum = Continuum.gridConsistency(meanEn, comp.grid)
                        contSettings = Continuum.Settings(false, nrContinuum);
                        npot         = Nuclear.nuclearPotential(comp.nuclearModel, comp.grid)
                        ## wp1 = compute("radial potential: core-Hartree", grid, wLevel)
                        ## wp2 = compute("radial potential: Hartree-Slater", grid, wLevel)
                        ## wp3 = compute("radial potential: Kohn-Sham", grid, wLevel)
                        ## wp           = Basics.compute("radial potential: Dirac-Fock-Slater", comp.grid, step.finalMultiplet.levels[1].basis)
                        wp           = Basics.computePotentialDFS(comp.grid, step.finalMultiplet.levels[1])
                        pot          = Basics.add(npot, wp)
                        cOrbital, phase, normF  = Continuum.generateOrbitalLocalPotential(meanEn, sh, pot, contSettings)
                        cOrbitals[sh] = cOrbital
                    end
                    println(">> Generate continum orbitals in DFS potential for $(keys(cOrbitals)) and for the energy $meanEn ")
                end

                newLines = AutoIonization.computeLinesFromOrbitals(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                                   step.settings, cOrbitals, output=true, printout=false)
                append!(linesA, newLines);    nt = length(linesA)
            elseif  step.process == Basics.Auger()
                # Compute continuum orbitals independently for all transitions in the given block.
                newLines = AutoIonization.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                              step.settings, output=true, printout=false)
                append!(linesA, newLines);    nt = length(linesA)
            elseif  step.process == Basics.Radiative()
                newLines = PhotoEmission.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.grid,
                                                             step.settings, output=true, printout=false)
                append!(linesR, newLines);    nt = length(linesR)
            else   error("Unsupported atomic process for cascade computations.")
            end
            println("     Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, giving now rise " *
                    "to a total of $nt $(string(step.process)) decay lines." )
            if  printSummary   println(iostream, "\n*    Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, " *
                                                 "giving now rise to a total of $nt $(string(step.process)) decay lines." )   end
        end
        #
        data = Cascade.DecayData(linesR, linesA)
    end

    """
    `Cascade.determineSteps(scheme::Cascade.StepwiseDecayScheme, comp::Cascade.Computation, blockList::Array{Cascade.Block,1})`
        ... determines all step::Cascade.Step's that need to be computed for this decay cascade. It cycles through all processes of the given
            decay scheme and selects all pairs of blocks due to the selected cascade approach. It is checked that the (averaged) energies
            each block or level supports a `decay' within the step. A stepList::Array{Cascade.Step,1} is returned, and for which subsequently
            all required transition amplitudes and rates are computed.
    """
    function determineSteps(scheme::Cascade.StepwiseDecayScheme, comp::Cascade.Computation, blockList::Array{Cascade.Block,1})
        stepList = Cascade.Step[]
        if  comp.approach  in  [Cascade.AverageSCA(), Cascade.SCA()]
            for  a = 1:length(blockList)
                for  b = 1:length(blockList)
                    minEn = 100000.;   maxEn = -100000.
                    for  p = 1:length(blockList[a].multiplet.levels),  q = 1:length(blockList[b].multiplet.levels)
                        minEn = min(minEn, blockList[a].multiplet.levels[p].energy - blockList[b].multiplet.levels[q].energy)
                        maxEn = max(maxEn, blockList[a].multiplet.levels[p].energy - blockList[b].multiplet.levels[q].energy)
                    end
                    for  process  in  comp.scheme.processes
                        if      process == Basics.Radiative()
                            if  a == b   ||   minEn < 0.    continue   end
                            if  blockList[a].NoElectrons == blockList[b].NoElectrons
                                settings = PhotoEmission.Settings([E1], [UseBabushkin], false, false, LineSelection(), 0., 0., 1.0e6)
                                push!( stepList, Cascade.Step(process, settings, blockList[a].confs, blockList[b].confs,
                                                              blockList[a].multiplet, blockList[b].multiplet) )
                            end
                        elseif  process == Basics.Auger()
                            if  a == b   ||   minEn < 0.    continue   end
                            if  blockList[a].NoElectrons == blockList[b].NoElectrons + 1
                                settings = AutoIonization.Settings(false, false, LineSelection(), 0., 1.0e6, 3, CoulombInteraction())
                                push!( stepList, Cascade.Step(process, settings, blockList[a].confs, blockList[b].confs,
                                                              blockList[a].multiplet, blockList[b].multiplet) )
                            end
                        else    error("stop a")
                        end
                    end
                end
            end
            #
        else  error("Unsupported cascade approach.")
        end
        return( stepList )
    end

    """
    `Cascade.perform(scheme::StepwiseDecayScheme, comp::Cascade.Computation)`
        ... to set-up and perform a cascade computation that starts from a given set of initial configurations and proceeds via
            various steps until a given number of electrons has been removed or the decay stops at some stable levels with regard
            to the given atomic processes. The results of all individual steps are printed to screen but nothing is returned
            otherwise.

    `Cascade.perform(scheme::StepwiseDecayScheme, comp::Cascade.Computation; output::Bool=false, outputToFile::Bool=true)`
        ... to perform the same but to return the complete output in a dictionary;  the particular output depends on the type
            and specifications of the cascade but can easily accessed by the keys of this dictionary.
    """
    function perform(scheme::StepwiseDecayScheme, comp::Cascade.Computation; output::Bool=false, outputToFile::Bool=true)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        # Perform the SCF and CI computation for the intial-state multiplets if initial configurations are given
        if  comp.initialConfigs != Configuration[]
            basis      = Basics.performSCF(comp.initialConfigs, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplet  = Basics.performCI(basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplets = [Multiplet("initial states", multiplet.levels)]
        else
            basis      = comp.initialMultiplets[1].levels[1].basis
            multiplets = comp.initialMultiplets
        end
        # Print out initial configurations and levels
        Cascade.displayLevels(stdout, multiplets, sa="initial ")
        if  printSummary   Cascade.displayLevels(iostream, multiplets, sa="initial ")          end
        #
        # Generate subsequent cascade configurations as well as display and group them together
        wa = Cascade.generateConfigurationList(multiplets, comp.scheme.maxElectronLoss, comp.scheme.NoShakeDisplacements)
        wb = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa, sa="decay ")
        #
        # Determine first all configuration 'blocks' and from them the individual steps of the cascade
        wc = Cascade.generateBlocks(comp, wb, basis.orbitals, sa="for the decay cascade:")
        Cascade.displayBlocks(stdout, wc, sa="for the decay cascade ")
        if  printSummary   Cascade.displayBlocks(iostream, wc, sa="for the decay cascade ")    end
        gMultiplets = Multiplet[];     for block in wc  push!(gMultiplets, block.multiplet)    end
        # Determine, modify and compute the transition data for all steps, ie. the PhotoEmission.Line's, the AutoIonization.Line's, etc.
        wd = Cascade.determineSteps(scheme, comp, wc)
        Cascade.displaySteps(stdout, wd, sa="decay ")
        if  printSummary   Cascade.displaySteps(iostream, wd, sa="decay ")    end
        we   = Cascade.modifySteps(wd)
        @time data = Cascade.computeSteps(scheme, comp, we)
        if output
            results = Base.merge( results, Dict("name"                  => comp.name) )
            results = Base.merge( results, Dict("cascade scheme"        => comp.scheme) )
            results = Base.merge( results, Dict("initial multiplets:"   => multiplets) )
            results = Base.merge( results, Dict("generated multiplets:" => gMultiplets) )
            results = Base.merge( results, Dict("decay line data:"      => data) )
            #
            #  Write out the result to file to later continue with simulations on the cascade data
            if outputToFile
                filename = "zzz-cascade-decay-computations-" * string(Dates.now())[1:13] * ".jld"
                println("\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                        "\n   results = JLD.load(''$filename'')    ... to load the results back from file.")
                if  printSummary   println(iostream, "\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                                                     "\n   results = JLD.load(''$filename'')    ... to load the results back from file." )      end
                ## JLD.save(filename, results)
                JLD2.@save filename results
            end
        else    results = nothing
        end

        return( results )
    end


text:
using Documenter, Perceptron

makedocs()


text:
import Slurm: api

type Configuration
  version
  ctl_conf :: Ptr{api.slurm_ctl_conf_t}

  function Configuration()
    version = api.slurm_api_version()
    timestamp = round(Int, Dates.datetime2unix(now()))
    r_conf = Ref{Ptr{api.slurm_ctl_conf_t}}()
    err = api.slurm_load_ctl_conf(timestamp, r_conf)
    # TODO: handle err

    conf = new(version, r_conf[])
    finalizer(conf, delete)
    return conf
  end
end

function delete(c :: Configuration)
  if c.ctl_conf != C_NULL
    api.slurm_free_ctl_conf(c.ctl_conf)
    c.ctl_conf = C_NULL
  end
end

function get(c :: Configuration, field :: Symbol)
  if field ∉ fieldnames(api.slurm_ctl_conf_t)
    error("$field is not a property of Configuration")
  else
    val = getfield(unsafe_load(c.ctl_conf), field)
    if isa(val, Cstring)
      return unsafe_string(val)
    else
      return val
    end
  end
end