text:
# setup and activate project environment
using Pkg
Pkg.activate(@__DIR__)

# setup conda environment
Pkg.add("Conda")
using Conda
Conda.add_channel("conda-forge")
Conda.add("python=3.7.8")
Conda.add("h5py=2.10.0")
Conda.add("pillow=8.0.1")
Conda.add("numpy=1.17")
Conda.add("matplotlib=3.1.2")
Conda.add("tensorflow-gpu=1.15")
Conda.add("tqdm=4.53.0")
Conda.pip_interop(true)
Conda.pip("install", "keras-adf==19.1.0")

# setup pycall and make it use the conda env
ENV["PYTHON"] = ""  # empty path defaults to local conda
Pkg.add("PyCall")

max_stars_repo_name: jmaces/fw-rde

max_issues_repo_name: ZIB-IOL/fw-rde

max_forks_repo_name: ZIB-IOL/fw-rde


text:
#
# Structure that contains the information about the solute and solvent molecules
#

struct Selection

  natoms :: Int # Total number of atoms
  nmols :: Int # Number of molecules
  natomspermol :: Int # Number of atoms per molecule

  index :: Vector{Int} # Index of each atom in the full vector of coordinates
  imol :: Vector{Int} # index of the molecule to which each atom belongs

  names :: Vector{String} # Types of the atoms, to be used in the atom-contributions

end

# Initialize providing the file name, and calling by default PDBTools.select

function Selection( file :: String, selection :: String;
                    nmols :: Int = 0, natomspermol :: Int = 0 )
  sel = PDBTools.readPDB(file,selection)
  return Selection( sel, nmols = nmols, natomspermol = natomspermol )
end

# If the input is a vector of PDBTools.Atom types, load the index and types

function Selection( atoms :: Vector{PDBTools.Atom}; nmols :: Int = 0, natomspermol :: Int = 0)
  indexes = [ at.index for at in atoms ]
  names = [ at.name for at in atoms ]
  return Selection( indexes, names, nmols=nmols, natomspermol=natomspermol )
end

# If no names are provided, just repeat the indexes

function Selection( indexes :: Vector{Int}; nmols :: Int = 0, natomspermol :: Int = 0)
  names = [ "$(index[i])" for i in 1:length(index) ]
  return Selection( indexes, names, nmols=nmols, natomspermol=natomspermol )
end

# Function to initialize the structures

function Selection( indexes :: Vector{Int}, names :: Vector{String};
                    nmols :: Int = 0, natomspermol :: Int = 0)

  if nmols == 0 && natomspermol == 0
    error("Set nmols or natomspermol when defining a selection.")
  end

  natoms = length(indexes)
  if natoms == 0
    error(" Vector of atom indexes is empty. ")
  end
  if length(names) == 0
    error(" Vector of atom names was explicitly provided but is empty. ")
  end
  if length(names) != natoms
    error(" The vector of atom indexes has a different number of elements than the vector of atom names. ")
  end
  if nmols != 0
    if natoms%nmols != 0
      error(" Number of atoms in selection must be a multiple of nmols.")
    end
    natomspermol = round(Int,natoms/nmols)
  else
   if natoms%natomspermol != 0
     error(" Number of atoms in selection must be a multiple of natomspermols.")
   end
   nmols = round(Int,natoms/natomspermol)
  end

  # Setting the vector that contains the index of the molecule of each atom

  imol = Vector{Int}(undef,natoms)
  iat = 0
  for j in 1:nmols
    for i in 1:natomspermol
      iat = iat + 1
      imol[iat] = j
    end
  end

  return Selection(natoms,nmols,natomspermol,indexes,imol,names[1:natomspermol])

end

import Base.show
function Base.show( io :: IO, s :: Selection )
  if s.nmols == 1
    mol = "molecule"
  else
    mol = "molecules"
  end
  if s.natoms == 1
    at = "atom"
  else
    at = "atoms"
  end
  println(" Selection of $(s.natoms) $at belonging to $(s.nmols) $mol. ")
end

ext: jl


text:
# Julia wrapper for header: /Users/jgoldfar/Documents/misc/julia/librsync/gen/../deps/usr/include/librsync.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function rs_trace_set_level(level::rs_loglevel)
    ccall((:rs_trace_set_level, librsync), Void, (rs_loglevel,), level)
end

# TODO: Find out what this is meant to be
const rs_trace_fn_t = Void
function rs_trace_to(arg1)
    ccall((:rs_trace_to, librsync), Void, (Ptr{rs_trace_fn_t},), arg1)
end

function rs_trace_stderr(level::rs_loglevel, msg)
    ccall((:rs_trace_stderr, librsync), Void, (rs_loglevel, Cstring), level, msg)
end

function rs_supports_trace()
    ccall((:rs_supports_trace, librsync), Cint, ())
end

function rs_hexify(to_buf, from_buf, from_len::Cint)
    ccall((:rs_hexify, librsync), Void, (Cstring, Ptr{Void}, Cint), to_buf, from_buf, from_len)
end

function rs_unbase64(s)
    ccall((:rs_unbase64, librsync), Csize_t, (Cstring,), s)
end

function rs_base64(buf, n::Cint, out)
    ccall((:rs_base64, librsync), Void, (Ptr{Cuchar}, Cint, Cstring), buf, n, out)
end

function rs_strerror(r::rs_result)
    ccall((:rs_strerror, librsync), Cstring, (rs_result,), r)
end

function rs_mdfour(out, _in, arg1::Csize_t)
    ccall((:rs_mdfour, librsync), Void, (Ptr{Cuchar}, Ptr{Void}, Csize_t), out, _in, arg1)
end

function rs_mdfour_begin(md)
    ccall((:rs_mdfour_begin, librsync), Void, (Ptr{rs_mdfour_t},), md)
end

function rs_mdfour_update(md, in_void, n::Csize_t)
    ccall((:rs_mdfour_update, librsync), Void, (Ptr{rs_mdfour_t}, Ptr{Void}, Csize_t), md, in_void, n)
end

function rs_mdfour_result(md, out)
    ccall((:rs_mdfour_result, librsync), Void, (Ptr{rs_mdfour_t}, Ptr{Cuchar}), md, out)
end

function rs_format_stats(stats, buf, size::Csize_t)
    ccall((:rs_format_stats, librsync), Cstring, (Ptr{rs_stats_t}, Cstring, Csize_t), stats, buf, size)
end

function rs_log_stats(stats)
    ccall((:rs_log_stats, librsync), Cint, (Ptr{rs_stats_t},), stats)
end

function rs_free_sumset(arg1)
    ccall((:rs_free_sumset, librsync), Void, (Ptr{rs_signature_t},), arg1)
end

function rs_sumset_dump(arg1)
    ccall((:rs_sumset_dump, librsync), Void, (Ptr{rs_signature_t},), arg1)
end

function rs_job_iter(job, buffers)
    ccall((:rs_job_iter, librsync), rs_result, (Ptr{rs_job_t}, Ptr{rs_buffers_t}), job, buffers)
end

# TODO: Find out what the below is meant to be
const rs_driven_cb = Any
function rs_job_drive(job, buf, in_cb::rs_driven_cb, in_opaque, out_cb::rs_driven_cb, out_opaque)
    ccall((:rs_job_drive, librsync), rs_result, (Ptr{rs_job_t}, Ptr{rs_buffers_t}, rs_driven_cb, Ptr{Void}, rs_driven_cb, Ptr{Void}), job, buf, in_cb, in_opaque, out_cb, out_opaque)
end

function rs_job_statistics(job)
    ccall((:rs_job_statistics, librsync), Ptr{rs_stats_t}, (Ptr{rs_job_t},), job)
end

function rs_job_free(arg1)
    ccall((:rs_job_free, librsync), rs_result, (Ptr{rs_job_t},), arg1)
end

function rs_sig_begin(new_block_len::Csize_t, strong_sum_len::Csize_t, sig_magic::rs_magic_number)
    ccall((:rs_sig_begin, librsync), Ptr{rs_job_t}, (Csize_t, Csize_t, rs_magic_number), new_block_len, strong_sum_len, sig_magic)
end

function rs_delta_begin(arg1)
    ccall((:rs_delta_begin, librsync), Ptr{rs_job_t}, (Ptr{rs_signature_t},), arg1)
end

function rs_loadsig_begin(arg1)
    ccall((:rs_loadsig_begin, librsync), Ptr{rs_job_t}, (Ptr{Ptr{rs_signature_t}},), arg1)
end

function rs_build_hash_table(sums)
    ccall((:rs_build_hash_table, librsync), rs_result, (Ptr{rs_signature_t},), sums)
end

# TODO: Find out what the below is meant to be
const rs_copy_cb = Void
function rs_patch_begin(copy_cb, copy_arg)
    ccall((:rs_patch_begin, librsync), Ptr{rs_job_t}, (Ptr{rs_copy_cb}, Ptr{Void}), copy_cb, copy_arg)
end

# TODO: Unify this type with the RawFD
const FILE = Void
function rs_sig_file(old_file, sig_file, block_len, strong_len, sig_magic, stats)
    ccall((:rs_sig_file, librsync), rs_result, (Ptr{FILE}, Ptr{FILE}, Csize_t, Csize_t, ANONYMOUS_2#=rs_magic_number=#, Ptr{rs_stats_t}), old_file, sig_file, block_len, strong_len, sig_magic, stats)
end

function rs_loadsig_file(sig_file, sumset, stats)
    ccall((:rs_loadsig_file, librsync), rs_result, (Ptr{FILE}, Ptr{Ptr{rs_signature_t}}, Ptr{rs_stats_t}), sig_file, sumset, stats)
end

function rs_file_copy_cb(arg, pos::rs_long_t, len, buf)
    ccall((:rs_file_copy_cb, librsync), rs_result, (Ptr{Void}, rs_long_t, Ptr{Csize_t}, Ptr{Ptr{Void}}), arg, pos, len, buf)
end

function rs_delta_file(arg1, new_file, delta_file, arg2)
    ccall((:rs_delta_file, librsync), rs_result, (Ptr{rs_signature_t}, Ptr{FILE}, Ptr{FILE}, Ptr{rs_stats_t}), arg1, new_file, delta_file, arg2)
end

function rs_patch_file(basis_file, delta_file, new_file, arg1)
    ccall((:rs_patch_file, librsync), rs_result, (Ptr{FILE}, Ptr{FILE}, Ptr{FILE}, Ptr{rs_stats_t}), basis_file, delta_file, new_file, arg1)
end


text:
# Macros:
#   reordering loops---done  (add @simd?)
#   tiling---done
#   parallelization
#   loop unrolling macro?
#   automatic bounds-checking of inputs? the infrastructure is basically in place already
#   writing generated code to an executable file

# add?: @test_looporder_simd and @loophoist_simd
#
# Parallelization:
#    @thread n for n = 1:size(A,2), ...
#    @thread (m,n) for n = 1:size(A,2), ...
# An issue will be allocation of temporaries in tiles; how to express this, if @thread
# is supposed to a separate step? Make need @tile_thread? Or, any allocations that come
# before the for loop could be interpreted as needing to be done once for each thread.

const assignment_heads = (:(=), :(+=), :(-=), :(*=), :(/=))

# Debugging: show useful information. Recognized settings:
#    :inferredbounds    Shows the sizes and offsets of temporaries
#    :tileindex         Shows the current outer tile variable values
#    :bodyindex         Shows (outer, inner) iteration variables in the body
#    :preindex          Shows (outer, inner) iteration variables in the pre-expression(s)
const KT_DEBUG = Dict{Symbol, Any}()
function setdebug(sym::Symbol, val)
    global KT_DEBUG
    KT_DEBUG[sym] = val
end

### loop ordering & hoisting
# In a loop nest like this:
# for n = 1:N
#     for m = 1:M
#         for k = 1:K
#             C[k, m, n] = A[m,n]*B[k]
#         end
#     end
# end
# the expression A[m, n] might be replaced by a stack-allocated variable `sym` assigned
# prior to the k loop.
type HoistInfo
    sym::Symbol
    loopstage::Int    # how many nestings deep at which the indexes are constant
    isassignedto::Bool
end

macro test_looporder(ex)
    loopVars, loopRanges = extract_loopvars(ex)
    nvars = length(loopVars)
    body = ex.args[2]
    exout = Expr(:block, Any[])
    for lperm in permutations(1:nvars)
        loopvars = loopVars[lperm]
        loopranges = loopRanges[lperm]
        newbody = gen_hoisted(body, loopvars, loopranges)
        push!(exout.args, :(println($loopvars); @time $newbody))
    end
    exout
end

macro loophoist(ex)
    loopvars, loopranges = extract_loopvars(ex)
    body = ex.args[2]
    gen_hoisted(body, loopvars, loopranges)
end

function extract_loopvars(ex::Expr)
    if ex.head != :for
        error("Must be a for expression")
    end
    loopVars = Symbol[]
    loopRanges = Any[]
    if ex.args[1].head == :(=)
        push!(loopVars, ex.args[1].args[1])
        push!(loopRanges, ex.args[1].args[2])
    elseif ex.args[1].head == :block
        for a in ex.args[1].args
            a.head == :(=) || error("Malformed loop-index block")
            push!(loopVars, a.args[1])
            push!(loopRanges, a.args[2])
        end
    else
        error("Cannot parse loop variables")
    end
    loopVars, loopRanges
end

function gen_hoisted(body, loopvars, loopranges)
    nvars = length(loopvars)
    hoisted, newbody = hoist(body, loopvars)
    # Ordering hoists by loopstage
    hoistex = collect(keys(hoisted))
    vals = collect(values(hoisted))
    hoistsym = [v.sym for v in vals]
    hoiststage = [v.loopstage for v in vals]
    hoistassignedto = [v.isassignedto for v in vals]
    p = sortperm(hoiststage)
    hoistex = hoistex[p]
    hoistsym = hoistsym[p]
    hoiststage = hoiststage[p]
    ## Code generation
    newbody = esc(newbody)
    k = length(hoiststage)
    for ls = nvars:-1:1
        newbody = quote
            for $(esc(loopvars[ls])) = $(esc(loopranges[ls]))
                $newbody
            end
        end
        while k > 0 && hoiststage[k] == ls-1
            reassign = hoistassignedto[k] ? (:($(esc(hoistex[k])) = $(esc(hoistsym[k])))) : :nothing
            newbody = quote
                $(esc(hoistsym[k])) = $(esc(hoistex[k]))
                $newbody
                $reassign
            end
            k -= 1
        end
    end
    newbody
end

hoist(body::Expr, varnames) = hoist!(Dict{Expr, HoistInfo}(), copy(body), varnames)
function hoist!(hoisted::Dict{Expr, HoistInfo}, exex::Expr, varnames)
    # This examines the arguments rather than the head expression, because we need to be able to replace
    # expressions inside the "container." One place where this will have trouble is if the top-level expression
    # is a single :ref expression, but that doesn't seem likely/useful in a real computation.
    isassignment = in(exex.head, assignment_heads)
    for k = 1:length(exex.args)
        isassignedto = isassignment && k == 1
        if isa(exex.args[k], Expr)
            ex::Expr = exex.args[k]
            # Check inside ex
            hoist!(hoisted, ex, varnames)
            # Now check ex itself
            if ex.head == :ref
                ls = 0  # loopstage for this reference
                for i = 2:length(ex.args)
                    ls = max(ls, loopstage_indexes(ex.args[i], varnames))
                end
                if ls < length(varnames)
                    # We can hoist. But first check whether this expression is already in the list
                    if haskey(hoisted, ex)
                        hinfo = hoisted[ex]
                        s = hinfo.sym
                        if isassignedto
                            hinfo.isassignedto = isassignedto
                        end
                    else
                        s = ktgensym()
                        hoisted[ex] = HoistInfo(s, ls, isassignedto)
                    end
                    exex.args[k] = s
                end
            end
        end
    end
    hoisted, exex
end

loopstage_indexes(s::Symbol, varnames) = indexin_scalar(s, varnames)
function loopstage_indexes(s::Expr, varnames)
    ls = 0
    for a in s.args
        ls = max(ls, loopstage_indexes(a, varnames))
    end
    ls
end
loopstage_indexes(arg, varnames) = 0

# A gensym-variant that generates names that can be used to write kernels to a parsable file
# (The `##32199` symbols produced by gensym() are not parsable)
let counter = 0
global ktgensym
global ktgensym_counter_reset
ktgensym_counter_reset() = (counter = 0)
function ktgensym()
    counter += 1
    symbol(string("_kt_", counter))
end
end
ktgensym(s::Symbol) = symbol("_kt_"*string(s))

# like indexin but first argument must be a scalar
function indexin_scalar(a, b)
    i = 1
    for bb in b
        if a == bb
            return i
        end
        i += 1
    end
    0
end

### Tiling ###

macro tile(tilesettings, args...)
    _tile(tilesettings, args...)
end
function _tile(tilesettings, args...)
    if length(args) == 2
        pre = args[1]
        body = args[2]
    elseif length(args) == 1
        pre = nothing
        body = args[1]
    end
    tilevars = Symbol[]  # these are the variables we will be tiling over. Might be a subset of all loopvars.
    outervars = Symbol[] # the names used for the outer variables (the tile indexes)
    tilesize = Int[]     # the tile size (i.e., the step size for each outer variable)
    k = 1
    (isa(tilesettings, Expr) && tilesettings.head == :tuple) || error("First argument must be a tuple-expression")
    while k < length(tilesettings.args)
        isa(tilesettings.args[k], Symbol) || error("error parsing tilesettings")
        push!(tilevars, tilesettings.args[k])
        k += 1
        if isa(tilesettings.args[k], Int)
            tv = symbol("_kt_outer_"*string(tilevars[end]))
        elseif isa(tilesettings.args[k], Symbol)
            tv = tilesettings.args[k]
            k += 1
        else
            error("error parsing tilesettings")
        end
        push!(outervars, tv)
        push!(tilesize, tilesettings.args[k])
        k += 1
    end
    # Create variables to hold the value of the tilesize.
    # We could have hard-coded these sizes, but using a variable makes it easier to
    # write code to loop over different settings (as in @test_tile).
    tilesizesym = [symbol("_kt_tsz_"*string(tv)) for tv in tilevars]
    ex = gen_tiled(tilevars, outervars, tilesizesym, pre, body)
    # Set the tilesize variables to those specified in the tuple provided by the user.
    settilesize = Expr(:block, Any[:(const $(esc(tilesizesym[i])) = $(tilesize[i])) for i = 1:length(tilesize)]...)
    return quote
        $settilesize
        $ex
    end
end

function gen_tiled(tilevars, outervars, tilesizesym, pre, body::Expr)
    initblocks = Any[]  # Statements to allocate temporaries used on a per-tile basis
    preblocks = Any[]   # Statements that run at the beginning of working on each tile
    innervars = [symbol("_kt_inner_"*string(tv)) for tv in tilevars]
    # Look for the initial "for" statement. It might be buried inside @loophoist
    body = copy(body)
    bodymod = body  # will point to the head of the "real" body expression
    if bodymod.head == :block
        bodymod = Base.is_linenumber(bodymod.args[1]) ? bodymod.args[2] : bodymod.args[1]
    end
    if is_macrocall(bodymod, :KernelTools, symbol("@loophoist"))
        bodymod = bodymod.args[2]
    end
    if bodymod.head != :for
        error("@tile requires that the body expression starts with a for loop; got ", bodymod)
    end
    # Extract all loop variables and their ranges (even those that we're not tiling over)
    looprangeexpr = bodymod.args[1]
    looprangedict = Dict{Symbol,Any}()    # might be an expression, :(1:5), or a symbol, :inner_x
    if looprangeexpr.head == :block
        # A multidimensional loop
        for ex in looprangeexpr.args
            looprangedict[ex.args[1]::Symbol] = copy(ex.args[2])
        end
    else
        # A single-variable loop
        looprangedict[looprangeexpr.args[1]::Symbol] = copy(looprangeexpr.args[2])
    end
    # Check that the tiled variables are a subset of the loop variables
    for tv in tilevars
        if !haskey(looprangedict, tv)
            error("Tile variables must be loop variables; did not find ", tv)
        end
    end

    # Create the outer loop nest and modify the inner loop nest to iterate over a tile
    outerlooprangeexprs = Expr[]
    innerlooprangeexprs = bodymod.args[1].head == :block ? bodymod.args[1].args : [bodymod.args[1]]
    outerorder = Symbol[]
    for i = 1:length(innerlooprangeexprs)
        ex = innerlooprangeexprs[i]
        lv = ex.args[1]
        ind = indexin_scalar(lv, tilevars)
        if ind > 0
            rng = looprangedict[lv]
            push!(outerlooprangeexprs, tilerange(outervars[ind], tilesizesym[ind], rng))
            push!(outerorder, lv)
            innerlooprangeexprs[i] = :($(innervars[ind]) = 0:min($(tilesizesym[ind])-1,last($rng)-$(outervars[ind])))
        end
    end
    if bodymod.args[1].head != :block
        bodymod.args[1] = innerlooprangeexprs[1]
    end
    p = sortperm(indexin(outerorder, tilevars))  # put in order specified by user
    outerlooprangeexprs = outerlooprangeexprs[p]

    # Handle any tilewise temporaries
    if pre != nothing
        ## Infer the sizes of the temporaries
        # Extract all references, and determine which temporaries will be needed
        tmpnames = assignments(pre)  # these are the temporary arrays we'll need
        refs = refexprs(pre, skiplhs=true)
        refexprs!(refs, body)        # refs contains exprs for all :ref expressions, other than setindex! statements in pre
        # Gather all indexes used for a given array
        # If A is used as A[x,y] and A[x+1,y], then
        #   indexes_by_array[:A] == Vector{Any}[[:x,:y], [:(x+1),:y]]
        indexes_by_array = Dict{Symbol,Vector{Vector{Any}}}()
        for r in refs
            arraysym = r.args[1]
            indexes = r.args[2:end]
            if !haskey(indexes_by_array, arraysym)
                indexes_by_array[arraysym] = Vector{Any}[]
            end
            push!(indexes_by_array[arraysym], indexes)
        end
        # From the set of indexing operations, compute bounds for each temporary.
        # This includes what are usually called "ghost cells" beyond the edges of tiles, so
        # that expressions like (tmparray[x-1] + tmparray[x+1]) will be in-bounds.
        # Temporaries will be indexed with x = innervar + offset, where innervar starts at 0.
        # The offset compensates for the smallest index in the stencil.
        offsetsyms = Any[]    # Names of variables holding the offsets for each temporary
        sizesyms = Any[]      # Names of variables holding the sizes of each temporary
        lastsyms = Any[]      # Names of variables holding the largest value of (outer+inner) along each axis
                              # at which the temporary needs to be initialized.
                              # (if tmp[x] = A[x] and out[x] = tmp[x]+tmp[x+1], we need to evaluate
                              #  tmp at one higher x than the range used for out)
        boundsblock = Any[]   # Statements needed to compute the bounds
        for (arrayname, indexes) in indexes_by_array
            if !in(arrayname, tmpnames)
                continue
            end
            osyms, ssyms, lsyms = constructbounds!(boundsblock, arrayname, indexes, tilevars, tilesizesym, looprangedict)
            push!(offsetsyms, osyms)
            push!(sizesyms, ssyms)
            push!(lastsyms, lsyms)
        end
        # Annotate the bounds computation so future analysis will know what this is about
        push!(initblocks, Expr(:block, esc(:(KT_BOUNDSCOMP_BLOCK = 1)), boundsblock...))
        # Allocate the temporaries. Now that we have the sizes, the only problem left is the types.
        # We'll compute one element of each temporary and use typeof.
        allocblock = Any[]
        for (i,Asym) in enumerate(tmpnames)
            # Find the first statement in pre that defines each temporary
            ex, _ = assignmentexpr(pre, Asym)
            # Change all indexes on the RHS to offsets. This is a safe place to evaluate it
            outfirst = [:(first($(looprangedict[tv]))) for tv in tilevars]
            ex1 = tileindex(ex.args[2], tilevars, outfirst, zeros(Int, length(tilevars)), tmpnames, offsetsyms[i])
            push!(allocblock, esc(:($Asym = Array(typeof($ex1), $(sizesyms[i]...)))))
        end
        push!(initblocks, Expr(:block, esc(:(KT_ALLOC_BLOCK = 1)), allocblock...))
        # Now we build code that will run inside the tile loops
        # Initialize each temporary
        for (i,Asym) in enumerate(tmpnames)
            # Here's a design problem: what about temporaries that take more than one line to compute?
            ex = assignmentblock(pre, Asym)  # find the expression for calculating the chosen temporary
            exa, _ = assignmentexpr(ex, Asym)
            osyms = offsetsyms[i]
            ssyms = sizesyms[i]
            lsyms = lastsyms[i]
            tv = exa.args[1].args[2:end]  # the local tilevars used for initializing this temporary
            ind = indexin(tv, tilevars)
            any(ind.==0) && error("Pre expressions must use the same indexing variables as the overall loop")
            ov, iv = outervars[ind], innervars[ind]
            loopranges = [esc(:($(iv[d]) = 1-$(osyms[d]):min($(ssyms[d])-$(osyms[d]), $(lsyms[d])-$(ov[d])))) for d = 1:length(iv)]
            indswap = indexin(outerorder, tv)
            indswap = indswap[indswap .> 0]
            loopranges = loopranges[indswap]
            initbody = tileindex(esc(ex), tv, ov, iv, tmpnames, offsetsyms)
            if get(KT_DEBUG, :preindex, false)
                initbody = Expr(:block, esc(Expr(:macrocall, symbol("@show"), Expr(:tuple, iv...))), initbody)
                push!(preblocks, :(println("Initializing ", $(Expr(:quote, Asym)))))
            end
            push!(preblocks, Expr(:for, Expr(:block, loopranges...), initbody))
        end
        # Replace the body indexing to be tile-specific
        bodymod.args[2] = tileindex(bodymod.args[2], tilevars, outervars, innervars, tmpnames, offsetsyms)
    else
        bodymod.args[2] = tileindex(bodymod.args[2], tilevars, outervars, innervars, Symbol[], [])
    end
    # Create the final expression
    block = Any[]
    if !isempty(initblocks)
        push!(block, initblocks...)
    end
    if get(KT_DEBUG, :bodyindex, false)
        unshift!(bodymod.args[2].args, Expr(:macrocall, symbol("@show"), Expr(:tuple, innervars...)))
    end
    body = esc(body)
    if get(KT_DEBUG, :tileindex, false)
        unshift!(preblocks, esc(Expr(:macrocall, symbol("@show"), Expr(:tuple, outervars...))))
    end
    if !isempty(preblocks)
        body = Expr(:block, Expr(:block, esc(:(KT_TILEWISE_BLOCK = 1)), preblocks...), body)
    end
    push!(block, Expr(:for, Expr(:block, outerlooprangeexprs...), body))
    Expr(:block, block...)
end

# On output, stmnts will have the expressions needed to compute the size parameters (offset, size, last) of the given
# temporary array.
# See more detailed description of the outputs above in the call to constructbounds!
function constructbounds!(stmnts::Vector{Any}, arrayname, indexes, tilevars, tilesizesym, looprangedict)
    nd = length(indexes[1])
    if !all(x->length(x) == nd, indexes)
        error("indexes are not all of the same dimensionality: $indexes")
    end
    # Parse all indexing operations on the temporary
    indexexprs = [Any[] for d = 1:nd]  # Indexing expressions evaluated at first&last indexes of a tile
    lastexprs = [Any[] for d = 1:nd]    # Used for the final tile along an axis
    for I in indexes
        for d = 1:nd
            ind = AffineIndex(I[d])
            tvindex = indexin_scalar(ind.sym, tilevars)
            rng = looprangedict[ind.sym]
            if tvindex != 0
                # when the inner tile variable is at min (0)
                push!(indexexprs[d], ind.offset)
                # when the inner tile variable is at max (tilesizesym-1)
                push!(indexexprs[d], :($(ind.coeff) * ($(tilesizesym[tvindex])-1) + $(ind.offset)))
            else
                # It's not a variable we're tiling over, so use the whole range
                push!(indexexprs[d], :(first($rng)))
                push!(indexexprs[d], :(last($rng)))
            end
            push!(lastexprs[d], :($(ind.coeff) * last($rng) + $(ind.offset)))
        end
    end
    # Add statements that compute the offset and size expressions
    offsetsyms = Array(Symbol, nd)
    sizesyms = Array(Symbol, nd)
    lastsyms = Array(Symbol, nd)
    toshow = get(KT_DEBUG, :inferredbounds, false)
    for d = 1:nd
        tag = string(arrayname)*"_"*string(d)
        # The offset is computed in terms of the minimum value of all indexing operations
        offsetsym = symbol("_kt_offset_"*tag)
        expr = Expr(:call, :min, indexexprs[d]...)
        push!(stmnts, esc(:(const $offsetsym = 1 - $expr)))
        sizesym = symbol("_kt_size_"*tag)
        lastsym = symbol("_kt_last_"*tag)
        # The size is computed in terms of the maximum value of all indexing operations
        expr = copy(expr)
        expr.args[1] = :max
        push!(stmnts, esc(:(const $sizesym = $offsetsym + $expr)))
        # The lastvalue is the maximum of all statements at total loop edges
        expr = length(lastexprs[d]) > 1 ? Expr(:call, :max, lastexprs[d]...) : lastexprs[d][1]
        push!(stmnts, esc(:(const $lastsym = $expr)))
        offsetsyms[d] = offsetsym
        sizesyms[d] = sizesym
        lastsyms[d] = lastsym
        if toshow
            push!(stmnts, esc(:(@show $offsetsym)))
            push!(stmnts, esc(:(@show $sizesym)))
            push!(stmnts, esc(:(@show $lastsym)))
        end
    end
    offsetsyms, sizesyms, lastsyms
end

# Get the symbol associated with assignment statements, e.g., in :(A[5] = 7) return :A
assignments(ex::Expr) = (asgn = Symbol[]; assignments!(asgn, ex))

function assignments!(asgn, ex::Expr)
    if in(ex.head, assignment_heads)
        if isa(ex.args[1], Expr)
            a = ex.args[1]::Expr
            if a.head == :ref
                push!(asgn, a.args[1]::Symbol)
            end
        end
    else
        for i = 1:length(ex.args)
            assignments!(asgn, ex.args[i])
        end
    end
    asgn
end
assignments!(asgn, a) = asgn

# Get the expression that first assigns a particular array symbol
function assignmentblock(ex::Expr, Asym::Symbol)
    if ex.head == :block
        args = ex.args
        for a in args
            _, found = assignmentexpr(a)
            if found
                return a
            end
        end
        error("Not found")
    end
    ex
end

function assignmentexpr(ex::Expr, Asym::Symbol)
    if in(ex.head, assignment_heads)
        if isa(ex.args[1], Expr)
            a = ex.args[1]::Expr
            if a.head == :ref && a.args[1] == Asym
                return ex, true
            end
        end
    end
    for i = 1:length(ex.args)
        if isa(ex.args[i], Expr)
            newex, found = assignmentexpr(ex.args[i]::Expr, Asym)
            if found
                return newex, true
            end
        end
    end
    return ex, false
end

# Get all reference expressions, i.e., :(A[i,j]). Optionally skip those that occur on the
# lhs of an assignment.
refexprs(ex::Expr; skiplhs=false) = (refs = Expr[]; refexprs!(refs, ex, skiplhs=skiplhs))

function refexprs!(refs, ex::Expr; skiplhs=false)
    if in(ex.head, assignment_heads)
        if isa(ex.args[1], Expr)
            a = ex.args[1]::Expr
            if !skiplhs && a.head == :ref
                push!(refs, a)
            end
        end
        refexprs!(refs, ex.args[2], skiplhs=skiplhs)
    else
        if ex.head == :ref
            push!(refs, ex)
        end
        for i = 1:length(ex.args)
            refexprs!(refs, ex.args[i], skiplhs=skiplhs)
        end
    end
    refs
end
refexprs!(refs, s; skiplhs=false) = refs

# Convert untiled indexing statements to tiled indexing statements
# For arrays that cover the whole domain, i -> i_outer + i_inner
# For temporary arrays defined just on a tile, i -> i_inner + array_i_offset
# Note that the pattern-matching is based on tilevars and tmpnames, and everything
# else is based on the position (i.e., index) within outervars/innervars/offsetsyms.
function tileindex(ex::Expr, tilevars, outervars, innervars, tmpnames, offsetsyms)
    domainexprs = [:($(outervars[i]) + $(innervars[i]))  for i = 1:length(tilevars)]
    exout = tileindex!(copy(ex), tilevars, domainexprs, innervars, tmpnames, offsetsyms)
    exout
end

function tileindex!(ex::Expr, tilevars, domainexprs, innervars, tmpnames, offsetsyms)
    if ex.head == :ref
        sym = ex.args[1]
        ind = indexin_scalar(sym, tmpnames)
        if ind > 0
            osyms = offsetsyms[ind]
            for i = 2:length(ex.args)
                ex.args[i] = replaceindex!(ex.args[i], tilevars, innervars, osyms[i-1])
            end
        else
            for i = 2:length(ex.args)
                ex.args[i] = replaceindex!(ex.args[i], tilevars, domainexprs)
            end
        end
    end
    for i = 1:length(ex.args)
        if isa(ex.args[i], Expr)
            tileindex!(ex.args[i]::Expr, tilevars, domainexprs, innervars, tmpnames, offsetsyms)
        elseif isa(ex.args[i], Symbol)
            s = ex.args[i]::Symbol
            ind = indexin_scalar(s, tilevars)
            ex.args[i] = ind > 0 ? domainexprs[ind] : s
        end
    end
    ex
end

replaceindex!(arg, tilevars, replacevars) = arg
function replaceindex!(s::Symbol, tilevars, replacevars)
    ind = indexin_scalar(s, tilevars)
    if ind > 0
        return replacevars[ind]
    end
    s
end
function replaceindex!(ex::Expr, tilevars, replacevars)
    for i = 1:length(ex.args)
        ex.args[i] = replaceindex!(ex.args[i], tilevars, replacevars)
    end
    ex
end
replaceindex!(arg, tilevars, replacevars, offset::Symbol) = arg
function replaceindex!(s::Symbol, tilevars, replacevars, offset::Symbol)
    ind = indexin_scalar(s, tilevars)
    if ind > 0
        return :($(replacevars[ind]) + $offset)
    end
#     :($s + $offset)
    s
end
function replaceindex!(ex::Expr, tilevars, replacevars, offset::Symbol)
    for i = 1:length(ex.args)
        ex.args[i] = replaceindex!(ex.args[i], tilevars, replacevars, offset)
    end
    ex
end

tilerange(outersym, tilesizesym, rng::Symbol) = esc(:($outersym = first($rng):$tilesizesym:last($rng)))
function tilerange(outersym, tilesizesym, rng::Expr)
    rng.head == :(:) || error("Must be a range expression")
    length(rng.args) == 2 || error("Sorry, step size not yet supported")
    esc(:($outersym = first($rng):$tilesizesym:last($rng)))
end

is_macrocall(a, modulename, sym::Symbol) = false
function is_macrocall(ex::Expr, modulename, sym::Symbol)
    if ex.head == :macrocall
        ex.args[1] == sym && return true
        if isa(ex.args[1], Expr)
            ex = ex.args[1]::Expr
            return ex.head == :. && ex.args[1] == modulename &&
                ((isa(ex.args[2], QuoteNode) && (ex.args[2]::QuoteNode).value == sym) ||
                 (isa(ex.args[2], Expr) && (ex.args[2]::Expr).head == :quote && (ex.args[2]::Expr).args[1] == sym))
        end
    end
    false
end

# Constant-folding for parsing indexes
# Parse indexing expressions like A[2k+1] to extract the coefficient (2), symbol (:k), and offset (1)
# Used for bounds-inference on temporary arrays. Each index must involve at most one symbol.
immutable AffineIndex
    sym::Symbol
    coeff::Int
    offset::Int
end
AffineIndex() = AffineIndex(:nothing, 0, 0)
AffineIndex(s::Symbol) = AffineIndex(s, 1, 0)
AffineIndex(i::Int) = AffineIndex(:nothing, 0, i)
function AffineIndex(ex::Expr)
    if ex.head == :call
        f = ex.args[1]
        if f == :+
            return +(AffineIndex(ex.args[2]), ex.args[3:end]...)
        elseif f == :-
            return -(AffineIndex(ex.args[2]), ex.args[3:end]...)
        elseif f == :*
            return *(AffineIndex(ex.args[2]), ex.args[3:end]...)
        else
            error("Operation $f not recognized")
        end
    else
        error("Expression appears not to be algebraic")
    end
end

(+)(a::AffineIndex, b::Int) = AffineIndex(a.sym, a.coeff, a.offset+b)
(-)(a::AffineIndex, b::Int) = AffineIndex(a.sym, a.coeff, a.offset-b)
(*)(a::AffineIndex, b::Int) = AffineIndex(a.sym, b*a.coeff, a.offset)
(+)(a::AffineIndex, b::Symbol) = a.sym == :nothing ? AffineIndex(b, a.coeff+1, a.offset) : error("cannot evaluate $a + $b")
(-)(a::AffineIndex, b::Symbol) = a.sym == :nothing ? AffineIndex(b, a.coeff-1, a.offset) : error("cannot evaluate $a - $b")
(*)(a::AffineIndex, b::Symbol) = a.sym == :nothing ? AffineIndex(b, a.offset, 0) : error("cannot evaluate $a * $b")
function (+)(a::AffineIndex, b::AffineIndex)
    sym = (a.sym == b.sym) ? a.sym : (a.sym == :nothing ? b.sym : (b.sym == :nothing ? a.sym : error("cannot combine affine expressions with different symbols")))
    AffineIndex(sym, a.coeff+b.coeff, a.offset+b.offset)
end
function (-)(a::AffineIndex, b::AffineIndex)
    sym = (a.sym == b.sym) ? a.sym : (a.sym == :nothing ? b.sym : (b.sym == :nothing ? a.sym : error("cannot combine affine expressions with different symbols")))
    AffineIndex(sym, a.coeff-b.coeff, a.offset-b.offset)
end
function (*)(a::AffineIndex, b::AffineIndex)
    sym = a.sym == :nothing ? b.sym : (b.sym == :nothing ? a.sym : error("cannot combine affine expressions with different symbols"))
    AffineIndex(sym, a.coeff*b.offset + a.offset*b.coeff, a.offset*b.offset)
end
(+)(a::AffineIndex, b::Expr) = a + AffineIndex(b)
(-)(a::AffineIndex, b::Expr) = a - AffineIndex(b)
(*)(a::AffineIndex, b::Expr) = a * AffineIndex(b)
(+)(a::AffineIndex, b...) = a + (+)(AffineIndex(b[1]), b[2:end]...)
(*)(a::AffineIndex, b...) = a * (*)(AffineIndex(b[1]), b[2:end]...)


text:
export getFaceMassMatrix

getFaceMassMatrix(mesh::TetraMesh, sigma::Array{Float64,1}) =
  getFaceMassMatrix(mesh, [repmat(sigma,1,3) zeros(length(sigma),3)])

function getFaceMassMatrix(mesh::TetraMesh, Sigma::Array{Float64,2})

	nc = mesh.nc
  nf = mesh.nf
  nn = mesh.nn

	M = zeros(4,4,nc)
   	N = round(Int128,nn)
   	globalFaceNumbering = arrayToUniqueInt(mesh.faces,N)
   	localFaceNumbering  = zeros(Int64,nc,4)

	for i=1:nc
		Ti = vec(mesh.Tetras[i,:])
	  	Xi = mesh.Points[Ti,:]

		M[:,:,i] = getLocalFaceMassMatrix(Sigma[i,:],Xi)

		localFaces = [
      Ti[1] Ti[2] Ti[3]
      Ti[1] Ti[2] Ti[4]
      Ti[1] Ti[3] Ti[4]
      Ti[2] Ti[3] Ti[4]]
		localFaceNumbering[i,:] = arrayToUniqueInt(localFaces,N)

	end
	faceInd = reshape(indexin(localFaceNumbering,globalFaceNumbering),nc,4)

	ii = zeros(Int64,16*nc)
	jj = zeros(Int64,16*nc)
  # mm = zeros(Float64,16*nc)

	cnt = 1
	for i=1:nc
		for i2=1:4
			for i1=1:4
				ii[cnt] = faceInd[i,i1]
				jj[cnt] = faceInd[i,i2]
        # mm[cnt] = M[i1,i2,i]
				cnt    += 1
			end
		end
	end
  mm = vec(M)

	return sparse(ii,jj,mm,nf,nf)
end

function getLocalFaceMassMatrix(sigma::Array{Float64},X::Array{Float64})
#              * X4
#             /|\
#            / | \
#           /  |  \
#          /   |   \
#         /    |    \
#     t6 /     |     \t5
#       /      |      \
#      /       |t3     \
#     /        |        \
#    /         |         \
# X3/_________ |t4_______ \ X2
#    .         |         .
#      .       |       .
#        .     |     . t1
#     t2   .   |   .
#            . * .
#             X1

	phi = [ zeros(3) for i = 1:4, j = 1:4 ]
	s = (1 - 1 / sqrt(5)) / 4
	t = 1 - 3 * s

	b1  = X[1, :]
	b2  = X[2, :]
	b3  = X[3, :]
	b4  = X[4, :]
  b21 = b2 - b1
  b31 = b3 - b1
  b41 = b4 - b1
	B  = [b21 b31 b41]

  A = getFaceArea(X[2, :], X[3, :], X[4, :])		  # face 234 - opposite vertex b1
	phi[1,4] =  A * B * [s;s;s]		            # phi_4
	phi[2,4] =  A * B * [t;s;s]
	phi[3,4] =  A * B * [s;t;s]
	phi[4,4] =  A * B * [s;s;t]

  A = getFaceArea(X[1, :], X[3, :], X[4, :])	# face 134 - opposite vertex b2
	phi[1,3] = -A * (B * [s;s;s] - b21)	# phi_3
	phi[2,3] = -A * (B * [t;s;s] - b21)
	phi[3,3] = -A * (B * [s;t;s] - b21)
	phi[4,3] = -A * (B * [s;s;t] - b21)

  A = getFaceArea(X[1, :], X[2, :], X[4, :])	# face 124 - opposite vertex b3
	phi[1,2] =  A * (B * [s;s;s] - b31)	# phi_2
	phi[2,2] =  A * (B * [t;s;s] - b31)
	phi[3,2] =  A * (B * [s;t;s] - b31)
	phi[4,2] =  A * (B * [s;s;t] - b31)

  A = getFaceArea(X[1, :], X[2, :], X[3, :])		  # face 123 - opposite vertex b4
	phi[1,1] = -A * (B * [s;s;s] - b41)	# phi_1
	phi[2,1] = -A * (B * [t;s;s] - b41)
	phi[3,1] = -A * (B * [s;t;s] - b41)
	phi[4,1] = -A * (B * [s;s;t] - b41)

	Sigma = [sigma[1] sigma[4]  sigma[5];
	         sigma[4] sigma[2]  sigma[6];
		       sigma[5] sigma[6]  sigma[3]]

   	M = zeros(4,4)
        for q = 1:4
	 for k = 1:4
          vk = Sigma * phi[q,k]
 	  for j = 1:4
            M[j,k] += dot(phi[q,j], vk)
          end
	 end
	end

	v = getCellVolume(X[1,:], X[2,:], X[3,:], X[4,:])
	return (1 / (36 * v)) * M
end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
using QuantumCumulants
using SymbolicUtils
using Test

@testset "fock" begin

# Test
hf = FockSpace(:c)

a = Destroy(hf,:a)
ad = a'

@test !isequal(hash(a), hash(ad))
b = Destroy(hf,:b)
@test !isequal(hash(a), hash(b))

@test isequal(a,ad')
@test isequal(simplify(a+a),2*a)
@test isequal(simplify(a/2 + 0.5*a),a)
@test isequal(a*a' , 1+a'*a)
@test isequal(simplify(a*a' + 1) , 2 + a'*a)

@test isequal(simplify(-1*(a'*a + 1)*a + a) , -1*a'*a^2)
@test isequal(simplify(a'*a*a - a*a'*a) , -1*a)

# Single mode
ωc = 0.1313
H = ωc*a'*a
da = simplify(1.0im*(H*a - a*H))
@test isequal(da , (0.0-1.0im)*ωc*a)

# Test substitute by numbers
@syms x
@test iszero(substitute(x*a, Dict(x=>0)))
@test isequal(substitute(x*a, Dict(x=>1)), a)
@test iszero(substitute(x*a, Dict(a=>0)))

# Test substitute by syms
@syms y
@test isequal(substitute(x*a, Dict(x=>y)), y*a)
@test isequal(substitute(x*a, Dict(a=>y)), x*y)
@test isequal(substitute(x*(a+a'), Dict(x => y)), y*(a + a'))

end # testset


text:
@testset "LibGEOS invalid geometry" begin
    # LibGEOS shouldn't crash but error out on invalid geometry

    # Self intersecting polygon
    polygon = LibGEOS._readgeom("POLYGON((0 0, 10 10, 0 10, 10 0, 0 0))")
    @test !LibGEOS.isValid(polygon)

    # Hole outside of base
    polygon = LibGEOS._readgeom("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),
        (15 15, 15 20, 20 20, 20 15, 15 15))")
    @test !LibGEOS.isValid(polygon)
end


text:
type AminoAcid
    name
    length
    content
end

type Alignment
    seqX # string
    seqY # string
    startPos   # (x, y)
    endPos
    score
end

type Group
    startCoord
    endCoord
    members
    memberNum
end

aminoAcids = Set()
push!(aminoAcids,
"A",
"R",
"N",
"D",
"B",
"C",
"E",
"Q",
"Z",
"G",
"H",
"I",
"L",
"K",
"M",
"F",
"P",
"S",
"T",
"W",
"Y",
"V")

function loadData(dataPath)
    data = open(dataPath)
    lines = readlines(data)

    aminoAcidlist = []

    for i = 1:2:length(lines)
        lineEles = split(lines[i], ",")
        contend = split(lines[1+i], "")

        for j = contend
            if !checkElement(j, aminoAcids)
                print("Wrong element: ")
                print(j)
                print("\n")
                return
            end
        end

        aminoAcid = AminoAcid(
            split(lineEles[1])[2],
            split(lineEles[2])[2],
            split(lines[1+i], "")
        )
        push!(aminoAcidlist, aminoAcid)
    end

    return aminoAcidlist
end

function checkElement(target, ele)
    for i = ele
        if i == target
            return true
        end
    end
    return false
end

#=
A   ALA
R   ARG
N   ASN
D   ASP
C   CYS
Q   GLN
E   GLU
G   GLY
H   HIS
I   ILE
L   LEU
K   LYS
M   MET
F   PHE
P   PRO
S   SER
T   THR
W   TRP
Y   TYR
V   VAL

A   1
R   2
N   3
D   4
C   5
Q   6
E   7
G   8
H   9
I   10
L   11
K   12
M   13
F   14
P   15
S   16
T   17
W   18
Y   19
V   20

4   -1  -2  -2  0   -1  -1  0   -2  -1  -1  -1  -1  -2  -1  1   0   -3  -2  0
-1  5   0   -2  -3  1   0   -2  0   -3  -2  2   -1  -3  -2  -1  -1  -3  -2  -3
-2  0   6   1   -3  0   0   0   1   -3  -3  0   -2  -3  -2  1   0   -4  -2  -3
-2  -2  1   6   -3  0   2   -1  -1  -3  -4  -1  -3  -3  -1  0   -1  -4  -3  -3
0   -3  -3  -3  9   -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1
-1  1   0   0   -3  5   2   -2  0   -3  -2  1   0   -3  -1  0   -1  -2  -1  -2
-1  0   0   2   -4  2   5   -2  0   -3  -3  1   -2  -3  -1  0   -1  -3  -2  -2
0   -2  0   -1  -3  -2  -2  6   -2  -4  -4  -2  -3  -3  -2  0   -2  -2  -3  -3
-2  0   1   -1  -3  0   0   -2  8   -3  -3  -1  -2  -1  -2  -1  -2  -2  2   -3
-1  -3  -3  -3  -1  -3  -3  -4  -3  4   2   -3  1   0   -3  -2  -1  -3  -1  3
-1  -2  -3  -4  -1  -2  -3  -4  -3  2   4   -2  2   0   -3  -2  -1  -2  -1  1
-1  2   0   -1  -3  1   1   -2  -1  -3  -2  5   -1  -3  -1  0   -1  -3  -2  -2
-1  -1  -2  -3  -1  0   -2  -3  -2  1   2   -1  5   0   -2  -1  -1  -1  -1  -1
-2  -3  -3  -3  -2  -3  -3  -3  -1  0   0   -3  0   6   -4  -2  -2  1   3   1
-1  -2  -2  -1  -3  -1  -1  -2  -2  -3  -3  -1  -2  -4  7   -1  -1  -4  -3  -2
1   -1  1   0   -1  0   0   0   -1  -2  -2  0   -1  -2  -1  4   1   -3  -2  -2
0   -1  0   -1  -1  -1  -1  -2  -2  -1  -1  -1  -1  -2  -1  1   5   -2  -2  0
-3  -3  -4  -4  -2  -2  -3  -2  -2  -3  -2  -3  -1  1   -4  -3  -2  11  2   -3
-2  -2  -2  -3  -2  -1  -2  -3  2   -1  -1  -2  -1  3   -3  -2  -2  2   7   -1
0   -3  -3  -3  -1  -2  -2  -3  -3  3   1   -2  -1  1   -2  -2  0   -3  -1  4

=#

using DataFrames

function getBlosum()
    dataPath = joinpath(@__DIR__, "..", "data")
    blosumPath = joinpath(dataPath, "ezBlosum.csv")
    blosum = readtable(blosumPath, separator = ',')
    return blosum
end

blosum = getBlosum()

blosumIndex = Dict(
    "A" => 1,
    "R" => 2,
    "N" => 3,
    "D" => 4,
    "C" => 5,
    "Q" => 6,
    "E" => 7,
    "G" => 8,
    "H" => 9,
    "I" => 10,
    "L" => 11,
    "K" => 12,
    "M" => 13,
    "F" => 14,
    "P" => 15,
    "S" => 16,
    "T" => 17,
    "W" => 18,
    "Y" => 19,
    "V" => 20
)

function getMappingScore(col, row, blosum)
    col = uppercase(col)
    row = uppercase(row)

    col = blosumIndex[col]
    row = blosumIndex[row]

    return blosum[col][row]
end

#=
This function searches for the maximum score left in the score matrix. If
it is >= the score threshold, then it will report the alignment corresponding
to it and return true, else false is returned.
@param matrix Score matrix
@param path Path matrix
@param xLen Length of seqX + 1
@param yLen Length of seqY + 1
@param alignmentList List of alignments to which we will append new alignment
@param threshold Lowest allowable score
@return True if alignment has sufficient score and is successfully added to
        list, else false
=#
function getAlignment(matrix, path, xLen, yLen, alignmentList, threshold, width, xAminoAcid, yAminoAcid)
    maxScore = 0
    col = 0
    row = 0
    ## First get the maxScore and its coordinates
    for r = 2:yLen
        for c = 2:xLen
            if matrix[r, c] > maxScore
                maxScore = matrix[r, c]
                col = c
                row = r
            end
        end
    end
    ## Check if maxScore >= threshold
    if maxScore < threshold
        return false ## score too low
    end
    ## Now retrieve the path
    seqX = ""
    seqY = ""
    x1 = 0
    x2 = col
    lastX = col
    y1 = 0
    y2 = row
    lastY = row

    #zero out the region around the maxScore
    for w in 0:width
        if row - w >= 2 && col - w >= 2
            matrix[row - w, col - w] = 0
        end
        if row - w >= 2 && col + w <= xLen
            matrix[row - w, col + w] = 0
        end
        if row + w <= yLen && col - w >= 2
            matrix[row + w, col - w] = 0
        end
        if row + w <= yLen && col + w <= xLen
            matrix[row + w, col + w] = 0
        end
    end
    pathVal = path[row, col] ## init path location
    while pathVal != 0
        ## First zero out the score element and nearby left/right region
        matrix[row, col] = 0
        for w in 1:width
            if col - w >= 2
                matrix[row, col - w] = 0
            end
            if col + w <= xLen
                matrix[row, col + w] = 0
            end
        end
        lastX = row
        lastY = col
        if pathVal == 1  # match/mismatch
            row -= 1
            col -= 1
            seqX = string(xAminoAcid.content[col], seqX)
            seqY = string(yAminoAcid.content[row], seqY)
        elseif pathVal == 2
            row -= 1
            seqX = string("-", seqX)
            seqY = string(yAminoAcid.content[row], seqY)
        else ## pathVal == 3
            col -= 1
            seqX = string(xAminoAcid.content[col], seqX)
            seqY = string("-", seqY)
        end
        pathVal = path[row, col]
    end
    x1 = lastX
    y1 = lastY
    ## Add this alignment to the list of alignments
    alignment = Alignment(
        seqX,
        seqY,
        (x1, y1),
        (y2, x2),
        maxScore
    )
    push!(alignmentList, alignment)
    return true
end

#=
This function returns a filtered list of alignments. If alignments are
overlapping (i.e. have the same start position), then only the one with
the highest score will be added.
@param alignmentList List of alignments
@return filtered list of non-overlapping alignments
=#
function filterAlignments(alignmentList)
    # sort the alignments first
    sort!(alignmentList, by = x -> x.startPos)
    filtList = []
    numFilt = 1
    i = length(alignmentList) - 1
    # push / pop the very first alignment
    push!(filtList, pop!(alignmentList))
    while i > 0
        if alignmentList[i].startPos == filtList[numFilt].startPos
            if alignmentList[i].score > filtList[numFilt].score
                pop!(filtList)
                push!(filtList,pop!(alignmentList))
            else
                pop!(alignmentList)
            end
        else # new start position
            push!(filtList,pop!(alignmentList))
            numFilt += 1
        end
        i -= 1
    end
    return filtList
end

#=
 = alignmentList is a list of Alignment object
 =#

function mergeAlignment(alignmentList)
    seqXGroups = [] # list of group x object
    seqYGroups = [] # list of group y object
    for alignment = alignmentList
        pushX = false
        for xGroup in seqXGroups
            if isCoordSimilar((xGroup.startCoord, xGroup.endCoord), (alignment.startPos[2], alignment.endPos[2]))
                push!(xGroup.members, alignment)
                xGroup.startCoord = (xGroup.startCoord * xGroup.memberNum + alignment.startPos[2]) / (xGroup.memberNum + 1)
                xGroup.endCoord = (xGroup.endCoord * xGroup.memberNum + alignment.endPos[2]) / (xGroup.memberNum + 1)
                xGroup.memberNum += 1
                pushX = true
            end
        end

        if pushX == false
            newSet = Set()
            push!(newSet, alignment)
            push!(
                seqXGroups,
                Group(
                    alignment.startPos[2],
                    alignment.endPos[2],
                    newSet,
                    1
                )
            )
        end

        pushY = false
        for yGroup in seqYGroups
            if isCoordSimilar((yGroup.startCoord, yGroup.endCoord), (alignment.startPos[1], alignment.endPos[1]))
                push!(yGroup.members, alignment)
                yGroup.startCoord = (yGroup.startCoord * yGroup.memberNum + alignment.startPos[1]) / (yGroup.memberNum + 1)
                yGroup.endCoord = (yGroup.endCoord * yGroup.memberNum + alignment.endPos[1]) / (yGroup.memberNum + 1)
                yGroup.memberNum += 1
                pushY = true
            end
        end

        if pushY == false
            newSet = Set()
            push!(newSet, alignment)
            push!(
                seqYGroups,
                Group(
                    alignment.startPos[1],
                    alignment.endPos[1],
                    newSet,
                    1
                )
            )
        end
    end
    return (seqXGroups, seqYGroups)
end

#=
    xCoord1 = (x1, x2) / (y1, y2)
    xCoord2 = (x1, x2) / (y1, y2)

    e.g.
    |xCoord1(x1) - xCoord2(x1)| <= distance (usually be 1)
    |xCoord1(x2) - xCoord2(x2)| <= distance (usually be 1)
=#
distance = 2
function isCoordSimilar(xCoord1, xCoord2)
    if (abs(xCoord1[1] - xCoord2[1]) <= distance) && (abs(xCoord1[2] - xCoord2[2]) <= distance)
        return true
    else
        return false
    end
end

function plotAlignmentProfile(alignmentListFilt, path, scoreMatrix)
    mark = 0
    for a in alignmentListFilt
        mark += 1
        row = a.endPos[1]
        col = a.endPos[2]
        score = []
        while path[row, col] != 0
            unshift!(score, scoreMatrix[row, col])
            if path[row, col] == 1
                row -= 1
                col -= 1
            elseif path[row, col] == 2 # indel from left
                col -= 1
            elseif path[row, col] == 3 # indel from right
                row -= 1
            end
        end
        graphTitle = string("Alignment", dec(mark), " Score Profile")

        alignmentScore = plot(x = collect(1:length(score)), y = score, Guide.xlabel("Alignment Position"), Guide.ylabel("Score"), Guide.title(graphTitle))
        fileName = string("Alignment", dec(mark))
        fileName = string(fileName, "ScoreProfile")
        fileName = string(fileName, ".svg")
        draw(SVG(fileName, 3inch, 3inch), alignmentScore)
    end
end

max_stars_repo_name: yguoogo/VPSLA

max_issues_repo_name: yguoogo/VPSLA

max_forks_repo_name: yguoogo/VPSLA


text:
"""
Tensor{StoreT,IndsT}

A plain old tensor (with order independent
interface and no assumption of labels)
"""
struct Tensor{ElT,N,StoreT<:TensorStorage,IndsT} <: AbstractArray{ElT,N}
  storage::StoreT
  inds::IndsT

  """
      Tensor{ElT,N,StoreT,IndsT}(inds, store::StorageType)

  Internal constructor for creating a Tensor from the
  storage and indices.

  The Tensor is a view of the tensor storage.

  For normal usage, use the Tensor(store::TensorStorage, inds)
  and tensor(store::TensorStorage, inds) constructors.
  """
  function Tensor{ElT,N,StoreT,IndsT}(
    ::AllowAlias, storage::TensorStorage, inds::Tuple
  ) where {ElT,N,StoreT<:TensorStorage,IndsT}
    return new{ElT,N,StoreT,IndsT}(storage, inds)
  end
end

function Tensor{ElT,N,StoreT,IndsT}(
  ::NeverAlias, storage::TensorStorage, inds
) where {ElT,N,StoreT<:TensorStorage,IndsT}
  return Tensor{ElT,N,StoreT,IndsT}(AllowAlias(), copy(storage), inds)
end

# Allow the storage and indices to be input in opposite ordering
(T::Type{<:Tensor})(as::AliasStyle, inds, storage::TensorStorage) = T(as, storage, inds)

"""
    Tensor(storage::TensorStorage, inds)

Construct a Tensor from a tensor storage and indices.
The Tensor holds a copy of the storage data.

The indices `inds` will be converted to a `Tuple`.
"""
function Tensor(as::AliasStyle, storage::TensorStorage, inds::Tuple)
  return Tensor{eltype(storage),length(inds),typeof(storage),typeof(inds)}(
    as, storage, inds
  )
end

# Automatically convert to Tuple if the indices are not a Tuple
# already (like a Vector). In the future this may be lifted
# to allow for very large tensor orders in which case Tuple
# operations may become too slow.
function Tensor(as::AliasStyle, storage::TensorStorage, inds)
  return Tensor(as, storage, Tuple(inds))
end

tensor(args...; kwargs...) = Tensor(AllowAlias(), args...; kwargs...)
Tensor(storage::TensorStorage, inds) = Tensor(NeverAlias(), storage, inds)
Tensor(inds, storage::TensorStorage) = Tensor(storage, inds)

storage(T::Tensor) = T.storage

# TODO: deprecate
store(T::Tensor) = storage(T)

data(T::Tensor) = data(storage(T))

datatype(T::Tensor) = datatype(storage(T))

indstype(::Type{<:Tensor{<:Any,<:Any,<:Any,IndsT}}) where {IndsT} = IndsT
indstype(T::Tensor) = indstype(typeof(T))

storagetype(::Type{<:Tensor{<:Any,<:Any,StoreT}}) where {StoreT} = StoreT
storagetype(T::Tensor) = storagetype(typeof(T))

# TODO: deprecate
storetype(args...) = storagetype(args...)

inds(T::Tensor) = T.inds

ind(T::Tensor, j::Integer) = inds(T)[j]

eachindex(T::Tensor) = CartesianIndices(dims(inds(T)))

eachblock(T::Tensor) = eachblock(inds(T))

eachdiagblock(T::Tensor) = eachdiagblock(inds(T))

eltype(::Tensor{ElT}) where {ElT} = ElT

strides(T::Tensor) = dim_to_strides(inds(T))

setstorage(T, nstore) = tensor(nstore, inds(T))

setinds(T, ninds) = tensor(storage(T), ninds)

#
# Generic Tensor functions
#

# The size is obtained from the indices
dims(T::Tensor) = dims(inds(T))
dim(T::Tensor) = dim(inds(T))
dim(T::Tensor, i::Int) = dim(inds(T), i)
maxdim(T::Tensor) = maxdim(inds(T))
mindim(T::Tensor) = mindim(inds(T))
diaglength(T::Tensor) = mindim(T)
size(T::Tensor) = dims(T)
size(T::Tensor, i::Int) = dim(T, i)

# Needed for passing Tensor{T,2} to BLAS/LAPACK
# TODO: maybe this should only be for DenseTensor?
function unsafe_convert(::Type{Ptr{ElT}}, T::Tensor{ElT}) where {ElT}
  return unsafe_convert(Ptr{ElT}, storage(T))
end

copy(T::Tensor) = setstorage(T, copy(storage(T)))

copyto!(R::Tensor, T::Tensor) = (copyto!(storage(R), storage(T)); R)

complex(T::Tensor) = setstorage(T, complex(storage(T)))

Base.real(T::Tensor) = setstorage(T, real(storage(T)))

Base.imag(T::Tensor) = setstorage(T, imag(storage(T)))

# Define Base.similar in terms of NDTensors.similar
Base.similar(T::Tensor, args...) = similar(T, args...)

#
# Necessary to overload since the generic fallbacks are
# slow
#

LinearAlgebra.norm(T::Tensor) = norm(storage(T))

conj(as::AliasStyle, T::Tensor) = setstorage(T, conj(as, storage(T)))
conj(T::Tensor) = conj(AllowAlias(), T)

randn!!(T::Tensor) = (randn!(T); T)
Random.randn!(T::Tensor) = (randn!(storage(T)); T)

LinearAlgebra.rmul!(T::Tensor, α::Number) = (rmul!(storage(T), α); T)
scale!(T::Tensor, α::Number) = rmul!(storage(T), α)

fill!!(T::Tensor, α::Number) = fill!(T, α)
fill!(T::Tensor, α::Number) = (fill!(storage(T), α); T)

-(T::Tensor) = setstorage(T, -storage(T))

#function similar(::Type{<:Tensor{ElT,N,StoreT}},dims) where {ElT,N,StoreT}
#  return tensor(similar(StoreT,dim(dims)),dims)
#end

# TODO: make sure these are implemented correctly
#similar(T::Type{<:Tensor},::Type{S}) where {S} = tensor(similar(storage(T),S),inds(T))
#similar(T::Type{<:Tensor},::Type{S},dims) where {S} = tensor(similar(storage(T),S),dims)

similar(T::Tensor) = setstorage(T, similar(storage(T)))

# TODO: for BlockSparse, this needs to include the offsets
# TODO: for Diag, the storage is not just the total dimension
similar(T::Tensor, dims::Tuple) = _similar_from_dims(T, dims)
function similar(::Type{TensorT}, dims::Tuple) where {TensorT<:Tensor}
  return _similar_from_dims(TensorT, dims)
end
function similar(::Type{TensorT}, dims::Tuple{}) where {TensorT<:Tensor}
  return _similar_from_dims(TensorT, dims)
end

# To handle method ambiguity with AbstractArray
#similar(T::Tensor,dims::Dims) = _similar_from_dims(T,dims)

similar(T::Tensor, ::Type{S}) where {S} = setstorage(T, similar(storage(T), S))

function similar(::Type{TensorT}, ::Type{S}, dims) where {TensorT<:Tensor,S<:Number}
  return _similar_from_dims(TensorT, S, dims)
end

similar(T::Tensor, ::Type{S}, dims) where {S<:Number} = _similar_from_dims(T, S, dims)

# To handle method ambiguity with AbstractArray
similar(T::Tensor, ::Type{S}, dims::Dims) where {S<:Number} = _similar_from_dims(T, S, dims)

_similar_from_dims(T::Tensor, dims) = tensor(similar(storage(T), dim(dims)), dims)

function _similar_from_dims(T::Tensor, ::Type{S}, dims) where {S<:Number}
  return tensor(similar(storage(T), S, dim(dims)), dims)
end

function _similar_from_dims(::Type{TensorT}, dims) where {TensorT<:Tensor}
  return _similar_from_dims(TensorT, eltype(TensorT), dims)
end

function _similar_from_dims(
  ::Type{TensorT}, ::Type{S}, dims
) where {TensorT<:Tensor,S<:Number}
  return tensor(similar(storagetype(TensorT), S, dim(dims)), dims)
end

function convert(
  ::Type{<:Tensor{<:Number,N,StoreR,Inds}}, T::Tensor{<:Number,N,<:Any,Inds}
) where {N,Inds,StoreR}
  return setstorage(T, convert(StoreR, storage(T)))
end

function zeros(TensorT::Type{<:Tensor{ElT,N,StoreT}}, inds) where {ElT,N,StoreT}
  return error("zeros(::Type{$TensorT}, inds) not implemented yet")
end

function promote_rule(
  ::Type{<:Tensor{ElT1,N1,StoreT1,IndsT1}}, ::Type{<:Tensor{ElT2,N2,StoreT2,IndsT2}}
) where {ElT1,ElT2,N1,N2,StoreT1,StoreT2,IndsT1,IndsT2}
  StoreR = promote_type(StoreT1, StoreT2)
  ElR = eltype(StoreR)
  return Tensor{ElR,N3,StoreR,IndsR} where {N3,IndsR}
end

function promote_rule(
  ::Type{<:Tensor{ElT1,N,StoreT1,Inds}}, ::Type{<:Tensor{ElT2,N,StoreT2,Inds}}
) where {ElT1,ElT2,N,StoreT1,StoreT2,Inds}
  StoreR = promote_type(StoreT1, StoreT2)
  ElR = eltype(StoreR)
  return Tensor{ElR,N,StoreR,Inds}
end

# Convert the tensor type to the closest dense
# type
function dense(::Type{<:Tensor{ElT,NT,StoreT,IndsT}}) where {ElT,NT,StoreT,IndsT}
  return Tensor{ElT,NT,dense(StoreT),IndsT}
end

dense(T::Tensor) = setstorage(T, dense(storage(T)))

function similartype(
  ::Type{<:Tensor{ElT,<:Any,StoreT,<:Any}}, ::Type{IndsR}
) where {ElT,StoreT,IndsR}
  return Tensor{ElT,length(IndsR),StoreT,IndsR}
end

function similartype(
  ::Type{<:Tensor{ElT,<:Any,StoreT,<:Any}}, ::Type{IndsR}
) where {ElT,StoreT,IndsR<:NTuple{NR}} where {NR}
  return Tensor{ElT,NR,StoreT,IndsR}
end

# Convert to Array, avoiding copying if possible
array(T::Tensor) = array(dense(T))
matrix(T::Tensor{<:Number,2}) = array(T)
vector(T::Tensor{<:Number,1}) = array(T)

isempty(T::Tensor) = isempty(storage(T))

#
# Helper functions for BlockSparse-type storage
#

"""
nzblocks(T::Tensor)

Return a vector of the non-zero blocks of the BlockSparseTensor.
"""
nzblocks(T::Tensor) = nzblocks(storage(T))

eachnzblock(T::Tensor) = eachnzblock(storage(T))

blockoffsets(T::Tensor) = blockoffsets(storage(T))
nnzblocks(T::Tensor) = nnzblocks(storage(T))
nnz(T::Tensor) = nnz(storage(T))
nblocks(T::Tensor) = nblocks(inds(T))
blockdims(T::Tensor, block) = blockdims(inds(T), block)
blockdim(T::Tensor, block) = blockdim(inds(T), block)

"""
offset(T::Tensor, block::Block)

Get the linear offset in the data storage for the specified block.
If the specified block is not non-zero structurally, return nothing.

offset(T::Tensor,pos::Int)

Get the offset of the block at position pos
in the block-offsets list.
"""
offset(T::Tensor, block) = offset(storage(T), block)

"""
isblocknz(T::Tensor,
          block::Block)

Check if the specified block is non-zero
"""
isblocknz(T::Tensor, block) = isblocknz(storage(T), block)

function blockstart(T::Tensor{<:Number,N}, block) where {N}
  start_index = @MVector ones(Int, N)
  for j in 1:N
    ind_j = ind(T, j)
    for block_j in 1:(block[j] - 1)
      start_index[j] += blockdim(ind_j, block_j)
    end
  end
  return Tuple(start_index)
end

function blockend(T::Tensor{<:Number,N}, block) where {N}
  end_index = @MVector zeros(Int, N)
  for j in 1:N
    ind_j = ind(T, j)
    for block_j in 1:block[j]
      end_index[j] += blockdim(ind_j, block_j)
    end
  end
  return Tuple(end_index)
end

#
# Some generic getindex and setindex! functionality
#

@propagate_inbounds @inline setindex!!(T::Tensor, x, I...) = setindex!(T, x, I...)

insertblock!!(T::Tensor, block) = insertblock!(T, block)

"""
getdiagindex

Get the specified value on the diagonal
"""
function getdiagindex(T::Tensor{<:Number,N}, ind::Int) where {N}
  return getindex(T, CartesianIndex(ntuple(_ -> ind, Val(N))))
end

"""
setdiagindex!

Set the specified value on the diagonal
"""
function setdiagindex!(T::Tensor{<:Number,N}, val, ind::Int) where {N}
  setindex!(T, val, CartesianIndex(ntuple(_ -> ind, Val(N))))
  return T
end

#
# Some generic contraction functionality
#

function zero_contraction_output(
  T1::TensorT1, T2::TensorT2, indsR::IndsR
) where {TensorT1<:Tensor,TensorT2<:Tensor,IndsR}
  return zeros(contraction_output_type(TensorT1, TensorT2, IndsR), indsR)
end

#
# Broadcasting
#

BroadcastStyle(::Type{T}) where {T<:Tensor} = Broadcast.ArrayStyle{T}()

function Base.similar(
  bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{T}}, ::Type{ElT}
) where {T<:Tensor,ElT}
  A = find_tensor(bc)
  return similar(A, ElT)
end

"`A = find_tensor(As)` returns the first Tensor among the arguments."
find_tensor(bc::Broadcast.Broadcasted) = find_tensor(bc.args)
find_tensor(args::Tuple) = find_tensor(find_tensor(args[1]), Base.tail(args))
find_tensor(x) = x
find_tensor(a::Tensor, rest) = a
find_tensor(::Any, rest) = find_tensor(rest)

function summary(io::IO, T::Tensor)
  for (dim, ind) in enumerate(inds(T))
    println(io, "Dim $dim: ", ind)
  end
  println(io, typeof(storage(T)))
  return println(io, " ", Base.dims2string(dims(T)))
end

#
# Printing
#

print_tensor(io::IO, T::Tensor) = Base.print_array(io, T)
print_tensor(io::IO, T::Tensor{<:Number,1}) = Base.print_array(io, reshape(T, (dim(T), 1)))


text:
@testset "trace" begin
    @testset "basic usage" begin
        global_tracer_provider(
            TracerProvider(
                span_processor = CompositSpanProcessor(
                    SimpleSpanProcessor(InMemoryExporter()),
                ),
            ),
        )

        tracer = Tracer()

        with_span(Span("test_async", tracer)) do
            @sync for i in 1:5
                @async with_span(Span("asyn_span_$i", tracer)) do
                    current_span()["my_id"] = i
                end
            end
        end

        p = global_tracer_provider()
        sp = p.span_processor.span_processors[1]
        @test length(sp.span_exporter.pool) == 6

        push!(p, SimpleSpanProcessor(ConsoleExporter()))

        with_span(Span("test", tracer)) do
            println("hello world!")
        end

        force_flush!(p)
        shut_down!(p)
        @test length(sp.span_exporter.pool) == 0

        s = Span("test", tracer)
        @test span_context(s) === INVALID_SPAN_CONTEXT
    end

    @testset "sampling" begin
        p = TracerProvider(
            span_processor = SimpleSpanProcessor(InMemoryExporter()),
            sampler = ALWAYS_OFF,
        )

        t = Tracer(; provider = p)

        s = Span("test", t)
        @test span_context(s) === INVALID_SPAN_CONTEXT

        p = TracerProvider(
            span_processor = SimpleSpanProcessor(InMemoryExporter()),
            sampler = DEFAULT_OFF,
        )

        t = Tracer(; provider = p)

        s = Span("test", t)
        @test span_context(s) === INVALID_SPAN_CONTEXT

        p = TracerProvider(
            span_processor = SimpleSpanProcessor(InMemoryExporter()),
            sampler = TraceIdRatioBased(0.5),
            id_generator = RandomIdGenerator(; rng = MersenneTwister(123)),
        )

        t = Tracer(; provider = p)
        n_no_op_spans = 0
        for _ in 1:1_000
            s = Span("test", t)
            if span_context(s) === INVALID_SPAN_CONTEXT
                n_no_op_spans += 1
            end
        end
        @test isapprox(n_no_op_spans / 1_000, 0.5; atol = 0.1)
    end
end


text:
#
# filters.jl --
#
# Implementation of basic filters.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2018, Éric Thiébaut.
#

"""
```julia
localmean(A, B)
```

yields the local mean of `A` in a neighborhood defined by `B`.  The result is
an array similar to `A`.

The in-place version is:

```julia
localmean!(dst, A, B) -> dst
```

See also [`localfilter!`](@ref).

"""
localmean(A::AbstractArray{T,N}, B=3) where {T,N} =
    localmean(A, Neighborhood{N}(B))

localmean(A::AbstractArray{T}, B::RectangularBox) where {T} =
    localmean!(similar(Array{float(T)}, axes(A)), A, B)

localmean(A::AbstractArray{T}, B::Kernel{Bool}) where {T} =
    localmean!(similar(Array{float(T)}, axes(A)), A, B)

localmean(A::AbstractArray{T}, B::Kernel{K}) where {T,K} =
    localmean!(similar(Array{float(promote_type(T,K))}, axes(A)), A, B)

function localmean!(dst::AbstractArray{Td,N},
                    A::AbstractArray{Ts,N}, B=3) where {Td,Ts,N}
    localmean!(dst, A, Neighborhood{N}(B))
end

@doc @doc(localmean) localmean!

function localmean!(dst::AbstractArray{Td,N},
                    A::AbstractArray{Ts,N},
                    B::RectangularBox{N}) where {Td,Ts,N}
    @assert axes(dst) == axes(A)
    T = _typeofsum(Ts)
    localfilter!(dst, A, B,
                 (a)     -> (zero(T), 0),
                 (v,a,b) -> (v[1] + a, v[2] + 1),
                 (d,i,v) -> _store!(d, i, v[1]/v[2]))
end

function localmean!(dst::AbstractArray{Td,N},
                    A::AbstractArray{Ts,N},
                    B::Kernel{Bool,N}) where {Td,Ts,N}
    @assert axes(dst) == axes(A)
    T = _typeofsum(Ts)
    localfilter!(dst, A, B,
                 (a)     -> (zero(T), 0),
                 (v,a,b) -> b ? (v[1] + a, v[2] + 1) : v,
                 (d,i,v) -> _store!(d, i, v[1]/v[2]))
end

function localmean!(dst::AbstractArray{Td,N},
                    A::AbstractArray{Ts,N},
                    B::Kernel{Tk,N}) where {Td,Ts,Tk,N}
    @assert axes(dst) == axes(A)
    T = _typeofsum(promote_type(Ts, Tk))
    localfilter!(dst, A, B,
                 (a)     -> (zero(T), zero(T)),
                 (v,a,b) -> (v[1] + a*b, v[2] + b),
                 (d,i,v) -> _store!(d, i, v[1]/v[2]))
end

"""

`_store!(arr, idx, val)` stores value `val` in array `arr` at index `idx`,
taking care of rounding `val` if it is of floating-point type while the
elements of `arr` are integers.

"""
@inline function _store!(arr::AbstractArray{T}, idx,
                         val::AbstractFloat) where {T<:Integer}
    arr[idx] = round(T, val)
end

@inline function _store!(arr::AbstractArray, idx, val)
    arr[idx] = val
end

"""

`_typeofsum(T)` yields a numerical type suitable for storing a sum of elements
of type `T`.

"""
_typeofsum(::Type{T}) where {T} = T
_typeofsum(::Type{T}) where {T<:Integer} =
    (sizeof(T) < sizeof(Int) ? widen(T) : T)

"""
```julia
convolve(A, B)
```

yields the convolution of `A` by the support of the neighborhood defined by
`B` of by the kernel `B` if it is an instance of `LocalFilters.Kernel` with
numerical coefficients.  The result is an array similar to `A`.

The in-place version is:

```julia
convolve!(dst, A, B) -> dst
```

See also [`localfilter!`](@ref).

"""
convolve(A::AbstractArray{T,N}, B=3) where {T,N} =
    convolve(A, Neighborhood{N}(B))

convolve(A::AbstractArray{T}, B::RectangularBox) where {T} =
    convolve!(similar(Array{_typeofsum(T)}, axes(A)), A, B)

convolve(A::AbstractArray{T}, B::Kernel{Bool}) where {T} =
    convolve!(similar(Array{_typeofsum(T)}, axes(A)), A, B)

convolve(A::AbstractArray{T}, B::Kernel{K}) where {T,K} =
    convolve!(similar(Array{_typeofsum(promote_type(T,K))}, axes(A)), A, B)

function convolve!(dst::AbstractArray{Td,N},
                   A::AbstractArray{Ts,N}, B=3) where {Td,Ts,N}
    convolve!(dst, A, Neighborhood{N}(B))
end

@doc @doc(convolve) convolve!

function convolve!(dst::AbstractArray{Td,N},
                   A::AbstractArray{Ts,N},
                   B::RectangularBox{N}) where {Td,Ts,N}
    @assert axes(dst) == axes(A)
    T = _typeofsum(Ts)
    localfilter!(dst, A, B,
                 (a)     -> zero(T),
                 (v,a,b) -> v + a,
                 _store!)
end

function convolve!(dst::AbstractArray{Td,N},
                   A::AbstractArray{Ts,N},
                   B::Kernel{Bool,N}) where {Td,Ts,N}
    @assert axes(dst) == axes(A)
    T = _typeofsum(Ts)
    localfilter!(dst, A, B,
                 (a)     -> zero(T),
                 (v,a,b) -> b ? v + a : v,
                 _store!)
end

function convolve!(dst::AbstractArray{Td,N},
                   A::AbstractArray{Ts,N},
                   B::Kernel{Tk,N}) where {Td,Ts,Tk,N}
    @assert axes(dst) == axes(A)
    T = _typeofsum(promote_type(Ts, Tk))
    localfilter!(dst, A, B,
                 (a)     -> zero(T),
                 (v,a,b) -> v + a*b,
                 _store!)
end

"""
# General local filters

A local filtering operation can be performed by calling:

```julia
localfilter!(dst, A, B, initial, update, store) -> dst
```

where `dst` is the destination, `A` is the source, `B` defines the
neighborhood, `initial`, `update` and `store` are three functions whose
purposes are explained by the following pseudo-code to implement the local
operation:

```
for i ∈ Sup(A)
    v = initial(A[i])
    for j ∈ Sup(A) and i - j ∈ Sup(B)
        v = update(v, A[j], B[i-j])
    end
    store(dst, i, v)
end
```

where `A` `Sup(A)` yields the support of `A` (that is the set of indices in
`A`) and likely `Sup(B)` for `B`.

For instance, to compute a local minimum (that is, an *erosion*):

```julia
localfilter!(dst, A, B,
             (a)     -> typemax(T),
             (v,a,b) -> min(v,a),
             (d,i,v) -> d[i] = v)
```

**Important:** Because the result of an elementary operation can be something
else than just a scalar, the loop(s) in `localfilter!` are performed without
bound checking of the destination and it is the caller's responsability to
insure that the destination have the correct sizes.

"""
function localfilter!(dst, A::AbstractArray{T,N}, B, initial::Function,
                      update::Function, store::Function) where {T,N}
    # Notes: The signature of this method is intentionally as little
    #        specialized as possible to avoid confusing the dispatcher.  The
    #        purpose of this method is just to convert `B ` into a neighborhood
    #        suitable for `A`.
    localfilter!(dst, A, Neighborhood{N}(B), initial, update, store)
end

#
# With `dst` the destination, `A` the source, and `B` the structuring
# element the pseudo-code to implement a local operation writes:
#
#     for i ∈ Sup(A)
#         v = initial(A[i])
#         for k ∈ Sup(B) and j = i - k ∈ Sup(A)
#             v = update(v, A[j], B[k])
#         end
#         dst[i] = final(v)
#     end
#
# where `A` `Sup(A)` yields the support of `A` (that is the set of indices in
# `A`) and likely `Sup(B)` for `B`.  Note that, in this example, destination
# `dst` and source `A` must have the same support.
#
# Equivalent form:
#
#     for i ∈ Sup(A)
#         v = initial(A[i])
#         for j ∈ Sup(A) and k = i - j ∈ Sup(B)
#             v = update(v, A[j], B[k])
#         end
#         dst[i] = final(v)
#     end
#
# In the second form, the bounds for `j` are:
#
#    imin ≤ j ≤ imax   and   kmin ≤ k = i - j ≤ kmax
#
# where `imin = initialindex(A)` and `imax = finalindex(A)` are the bounds for
# `A` while `kmin = initialindex(B)` and `kmax = finalindex(B)` are the bounds
# for `B`.  The above constraints are identical to:
#
#    max(imin, i - kmax) ≤ j ≤ min(imax, i - kmin)
#
function localfilter!(dst,
                      A::AbstractArray{T,N},
                      B::RectangularBox{N},
                      initial::Function,
                      update::Function,
                      store::Function) where {T,N}
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    @inbounds for i in R
        v = initial(A[i])
        @simd for j in cartesianregion(max(imin, i - kmax),
                                       min(imax, i - kmin))
            v = update(v, A[j], true)
        end
        store(dst, i, v)
    end
    return dst
end

#
# Pseudo-code for a local operation on `A` in a neighborhood `B` is:
#
#     for i ∈ Sup(A)
#         v = initial(A[i])
#         for j ∈ Sup(A) and i - j ∈ Sup(B)
#             v = update(v, A[j], kernel[i-j+off])
#         end
#         dst[i] = final(v)
#     end
#
# where `off` is the offset; the bounds for `j` are:
#
#    imin ≤ j ≤ imax   and   kmin ≤ i - j ≤ kmax
#
# where `imin` and `imax` are the bounds for `A` while `kmin` and `kmax` are
# the bounds for `B`.  The above constraints are identical to:
#
#    max(imin, i - kmax) ≤ j ≤ min(imax, i - kmin)
#

function localfilter!(dst,
                      A::AbstractArray{T,N},
                      B::Kernel{K,N},
                      initial::Function,
                      update::Function,
                      store::Function) where {T,K,N}
    R = cartesianregion(A)
    imin, imax = limits(R)
    kmin, kmax = limits(B)
    ker, off = coefs(B), offset(B)
    @inbounds for i in R
        v = initial(A[i])
        k = i + off
        @simd for j in cartesianregion(max(imin, i - kmax),
                                       min(imax, i - kmin))
            v = update(v, A[j], ker[k-j])
        end
        store(dst, i, v)
    end
    return dst
end


text:
include("./MyPlots/target-plot.jl")
import .TargetPlot as tp

tp.plot()

max_stars_repo_name: applied-math-coding/heuristic

max_issues_repo_name: applied-math-coding/heuristic

max_forks_repo_name: applied-math-coding/heuristic


text:
using DiffEqBase: ODEProblem, RODEProblem

function make_tangent_prob(phase_prob::RODEProblem, tangent_dynamics, u0)
    dim_phase = length(phase_prob.u0)
    rand_prototype = if phase_prob.rand_prototype === nothing
        similar(phase_prob.u0, dim_phase)
    else
        phase_prob.rand_prototype
    end

    return remake(
        phase_prob;
        f = tangent_dynamics,
        u0 = u0,
        rand_prototype = rand_prototype)
end

function with_output_noise(ode::ODEProblem,
                           g::Union{AbstractVector, Number};
                           kwargs...)
    function f(du, u, p, t, W)
        ode.f(du, u, p, t)
        du .+= g .* W
    end

    return RODEProblem(f, ode.u0, ode.tspan, ode.p; kwargs...)
end

function with_input_noise(ode::ODEProblem,
                          g::Union{AbstractVector, Number};
                          kwargs...)
    function f(du, u, p, t, W)
        u .+= g .* W
        ode.f(du, u, p, t)
    end

    return RODEProblem(f, ode.u0, ode.tspan, ode.p; kwargs...)
end


text:
module Commands
using ..StructTypes, ..Messages, ..Events
import ..Optional, ..FriendId, ..GroupOrFriendId, ..GroupId, ..TimeStamp, ..DurationSeconds
import ..MessageId, ..EventId, ..DeviceKind, ..MessageChains.MessageChain, ..EventOrMessage

abstract type AbstractCommand end
StructTypes.StructType(::Type{<:AbstractCommand}) = StructTypes.Struct()

module CommandMethods
@enum CommandMethod begin
    GET
    POST
    UPLOAD
end
end
const CommandMethod = CommandMethods.CommandMethod

command(cmd::AbstractCommand) = nameof(typeof(cmd))
subcommand(::AbstractCommand) = nothing

Base.@kwdef struct RESTful{T}
    code::Int
    msg::String
    data::Optional{T}
end
StructTypes.StructType(::Type{<:RESTful}) = StructTypes.Struct()

Base.@kwdef struct RESTfulErrorMessage{T}
    code::Int
    errorMessage::String
    data::Optional{T}
end
StructTypes.StructType(::Type{<:RESTfulErrorMessage}) = StructTypes.Struct()

StructTypes.@Struct struct VersionResponse
    version::VersionNumber
end

struct about <: AbstractCommand end
method(::about) = CommandMethods.GET
response_type(::about) = RESTful{VersionResponse}
response_type_compat(::about) = RESTfulErrorMessage{VersionResponse}

struct messageFromId <: AbstractCommand
    id::Int
end
method(::messageFromId) = CommandMethods.GET
response_type(::messageFromId) = RESTful{MessageChain}
response_type_compat(::messageFromId) = RESTfulErrorMessage{MessageChain}

abstract type AbstractListCommand <: AbstractCommand end
method(::AbstractListCommand) = CommandMethods.GET

Base.@kwdef struct friendList <: AbstractListCommand end
response_type(::friendList) = RESTful{Vector{Friend}}
response_type_compat(::friendList) = Vector{Friend}

Base.@kwdef struct groupList <: AbstractListCommand end
response_type(::groupList) = RESTful{Vector{Group}}
response_type_compat(::groupList) = Vector{Group}

Base.@kwdef struct memberList <: AbstractListCommand end
response_type(::memberList) = RESTful{Vector{Member}}
response_type_compat(::memberList) = Vector{Member}

module Sexes
@enum Sex begin
    UNKNOWN
    MALE
    FEMALE
end
end
const Sex = Sexes.Sex

StructTypes.@Struct struct Profile <: AbstractCommand
    nickname::String
    email::String
    age::Int
    level::Int
    sign::String
    sex::Sex
end

abstract type AbstractGetProfileCommand <: AbstractCommand end
method(::AbstractGetProfileCommand) = CommandMethods.GET
response_type(::AbstractGetProfileCommand) = Profile

Base.@kwdef struct botProfile <: AbstractGetProfileCommand end

Base.@kwdef struct friendProfile <: AbstractGetProfileCommand
    target::FriendId
end

Base.@kwdef struct memberProfile <: AbstractGetProfileCommand
    target::GroupId   # group id
    memberId::FriendId # member qq
end

StructTypes.@Struct struct MessageIdResponse
    code::Int
    msg::String
    messageId::Optional{MessageId}
end

abstract type AbstractMessagingCommand <: AbstractCommand end
StructTypes.names(::Type{<:AbstractMessagingCommand}) = ((:quoteId, :quote),)
StructTypes.omitempties(::Type{<:AbstractMessagingCommand}) = (:target, :qq, :group, :quoteId)
method(::AbstractMessagingCommand) = CommandMethods.POST
response_type(::AbstractMessagingCommand) = MessageIdResponse

Base.@kwdef struct sendFriendMessage <: AbstractMessagingCommand
    target::Optional{FriendId} = nothing
    qq::Optional{FriendId} = nothing # Choose one between target and qq
    quoteId::Optional{MessageId} = nothing
    messageChain::MessageChain
end

Base.@kwdef struct sendGroupMessage <: AbstractMessagingCommand
    target::Optional{GroupId} = nothing
    group::Optional{GroupId} = nothing # Choose one between target and group
    quoteId::Optional{MessageId} = nothing
    messageChain::MessageChain
end

Base.@kwdef struct sendTempMessage <: AbstractMessagingCommand
    qq::FriendId
    group::GroupId
    quoteId::Optional{MessageId} = nothing
    messageChain::MessageChain
end

Base.@kwdef struct sendNudge <: AbstractMessagingCommand
    target::FriendId  # target qq
    subject::GroupOrFriendId # group or friend id
    kind::NudgeKind
end
response_type(::sendNudge) = RESTful{Nothing}

Base.@kwdef struct recall <: AbstractMessagingCommand
    target::MessageId
end
response_type(::recall) = RESTful{Nothing}

abstract type AbstractFileCommand <: AbstractCommand end
StructTypes.omitempties(::Type{<:AbstractFileCommand}) = (:path, :target, :qq, :group, :offset, :size)
method(::AbstractFileCommand) = CommandMethods.POST

StructTypes.@Struct struct DownloadInfo
    sha1::String
    md5::String
    downloadTimes::Int
    uploaderId::FriendId
    uploadTime::TimeStamp
    lastModifyTime::TimeStamp
    url::String
end

StructTypes.@Struct struct FileInfo
    name::String
    id::Optional{String} # Empty for root
    path::String
    parent::Optional{FileInfo}
    contact::Union{Group,Friend}
    isFile::Bool
    isDirectory::Bool
    size::Int
    downloadInfo::Optional{DownloadInfo}
end

Base.@kwdef struct file_list <: AbstractFileCommand
    id::String # folder id
    path::Optional{String} = nothing # folder name
    target::Optional{GroupOrFriendId} = nothing
    group::Optional{GroupId} = nothing
    qq::Optional{FriendId} = nothing
    withDownloadInfo::Bool = false
    offset::Optional{Int} = 1 # page number
    size::Optional{Int} = 10  # page size
end
response_type(::file_list) = RESTful{Vector{FileInfo}}
method(::file_list) = CommandMethods.GET

Base.@kwdef struct file_info <: AbstractFileCommand
    id::String
    path::Optional{String} = nothing
    target::Optional{GroupOrFriendId} = nothing
    group::Optional{GroupId} = nothing
    qq::Optional{FriendId} = nothing
    withDownloadInfo::Bool = false
end
response_type(::file_info) = RESTful{FileInfo}
method(::file_info) = CommandMethods.GET

Base.@kwdef struct file_mkdir <: AbstractFileCommand
    id::String
    path::Optional{String} = nothing
    target::Optional{GroupOrFriendId} = nothing
    group::Optional{GroupId} = nothing
    qq::Optional{FriendId} = nothing
    directoryName::String
end
response_type(::file_mkdir) = RESTful{FileInfo}

Base.@kwdef struct file_delete <: AbstractFileCommand
    id::String
    path::Optional{String} = nothing
    target::Optional{GroupOrFriendId} = nothing
    group::Optional{GroupId} = nothing
    qq::Optional{FriendId} = nothing
end
response_type(::file_delete) = RESTful{Nothing}

Base.@kwdef struct file_move <: AbstractFileCommand
    id::String
    path::Optional{String} = nothing
    target::Optional{GroupOrFriendId} = nothing
    group::Optional{GroupId} = nothing
    qq::Optional{FriendId} = nothing
    moveTo::Optional{String} # destination folder id
    moveToPath::Optional{String}
end
response_type(::file_move) = RESTful{Nothing}

Base.@kwdef struct file_rename <: AbstractFileCommand
    id::String
    path::Optional{String} = nothing
    target::Optional{GroupOrFriendId} = nothing
    group::Optional{GroupId} = nothing
    qq::Optional{FriendId} = nothing
    renameTo::String
end
response_type(::file_rename) = RESTful{Nothing}

abstract type AbstractUploadCommand <: AbstractCommand end
method(::AbstractUploadCommand) = CommandMethods.UPLOAD

module UploadTypes
@enum UploadType begin
    FRIEND
    GROUP
    TEMP
end
end
const UploadType = UploadTypes.UploadType

Base.string(x::UploadType) = x == UploadTypes.FRIEND ? "friend" : x == UploadTypes.GROUP ? "group" : "temp"

StructTypes.@Struct struct ImageIdResponse
    imageId::String
    url::String
end

Base.@kwdef struct uploadImage <: AbstractUploadCommand
    type::UploadType
    img::Any
end
response_type(::uploadImage) = ImageIdResponse

StructTypes.@Struct struct VoiceIdResponse
    voiceId::String
    url::String
end

Base.@kwdef struct uploadVoice <: AbstractUploadCommand
    type::UploadType
    voice::Any
end
response_type(::uploadVoice) = VoiceIdResponse

Base.@kwdef struct file_upload <: AbstractUploadCommand
    type::UploadType
    target::GroupOrFriendId
    path::String
    file::Any # What ever that can be used by HTTP.Form
end
response_type(::file_upload) = RESTful{FileInfo}

Base.@kwdef struct deleteFriend <: AbstractCommand
    target::FriendId
end
method(::deleteFriend) = CommandMethods.POST
response_type(::deleteFriend) = RESTful{Nothing}

abstract type AbstractGroupManagementCommand <: AbstractCommand end
method(::AbstractGroupManagementCommand) = CommandMethods.POST
response_type(::AbstractGroupManagementCommand) = RESTful{Nothing}

Base.@kwdef struct mute <: AbstractGroupManagementCommand
    target::GroupId
    memberId::FriendId
    time::TimeStamp # seconds
end

Base.@kwdef struct unmute <: AbstractGroupManagementCommand
    target::GroupId
    memberId::FriendId
end

Base.@kwdef struct kick <: AbstractGroupManagementCommand
    target::GroupId
    memberId::FriendId
    msg::Optional{String} = nothing
end

Base.@kwdef struct quit <: AbstractGroupManagementCommand
    target::GroupId
end

Base.@kwdef struct muteAll <: AbstractGroupManagementCommand
    target::GroupId
end

Base.@kwdef struct unmuteAll <: AbstractGroupManagementCommand
    target::GroupId
end

Base.@kwdef struct setEssence <: AbstractGroupManagementCommand
    target::MessageId
end

Base.@kwdef struct memberAdmin <: AbstractGroupManagementCommand
    target::GroupId
    memberId::FriendId
    assign::Bool
end

Base.@kwdef struct GroupConfig
    name::Optional{String} = nothing # Group name
    announcement::Optional{String} = nothing
    confessTalk::Optional{Bool} = nothing
    allowMemberInvite::Optional{Bool} = nothing
    autoApprove::Optional{Bool} = nothing
    anonymousChat::Optional{Bool} = nothing
end
StructTypes.StructType(::Type{GroupConfig}) = StructTypes.Struct()

Base.@kwdef struct getGroupConfig <: AbstractGroupManagementCommand
    target::GroupId
end
command(::getGroupConfig) = :groupConfig
subcommand(::getGroupConfig) = :get
method(::getGroupConfig) = CommandMethods.GET
response_type(::getGroupConfig) = GroupConfig

Base.@kwdef struct updateGroupConfig <: AbstractGroupManagementCommand
    target::GroupId
    config::GroupConfig
end
command(::updateGroupConfig) = :groupConfig
subcommand(::updateGroupConfig) = :update
method(::updateGroupConfig) = CommandMethods.POST
response_type(::updateGroupConfig) = RESTful{Nothing}

Base.@kwdef struct getMemberInfo <: AbstractGroupManagementCommand
    target::GroupId
    memberId::FriendId
end
command(::getMemberInfo) = :memberInfo
subcommand(::getMemberInfo) = :get
method(::getMemberInfo) = CommandMethods.GET
response_type(::getMemberInfo) = Member

Base.@kwdef struct updateMemberInfo <: AbstractGroupManagementCommand
    target::GroupId
    memberId::FriendId
end
command(::updateMemberInfo) = :memberInfo
subcommand(::updateMemberInfo) = :update
method(::updateMemberInfo) = CommandMethods.POST
response_type(::updateMemberInfo) = RESTful{Nothing}

abstract type AbstractAnnouncementCommand <: AbstractGroupManagementCommand end
method(::AbstractAnnouncementCommand) = CommandMethods.POST

"Unique announcement identifyer"
const AnnouncementId = String

StructTypes.@Struct struct Announcement
    group::Group
    content::String
    senderId::FriendId
    fid::AnnouncementId
    allConfirmed::Bool
    confirmedMembersCount::Int
    publicationTime::TimeStamp
end

Base.@kwdef struct anno_list <: AbstractAnnouncementCommand
    id::GroupId
    offset:Int = 0
    size::Int = 10
end
response_type(::anno_list) = RESTful{Vector{Announcement}}

Base.@kwdef struct anno_publish <: AbstractAnnouncementCommand
    target::GroupId
    content::String
    pinned::Bool = false
end
response_type(::anno_publish) = RESTful{Vector{announcement}}

Base.@kwdef struct anno_delete <: AbstractAnnouncementCommand
    id::GroupId
    fid::AnnouncementId
end
response_type(::anno_delete) = RESTful{Nothing}

abstract type AbstractRequestCommand <: AbstractCommand end
method(::AbstractRequestCommand) = CommandMethods.POST
response_type(::AbstractRequestCommand) = RESTful{Nothing}

module NewFriendOperations
const ACCEPT = 0
const REFUSE = 1
const REFUSE_AND_BLACKLIST = 2
end

Base.@kwdef struct resp_newFriendRequestEvent <: AbstractRequestCommand
    eventId::EventId
    fromId::FriendId
    groupId::GroupId # 0 if not from a group
    operate::Int
    message::String
end

module MemberJoinOperations
const ACCEPT = 0
const REFUSE = 1
const IGNORE = 2
const REFUSE_AND_BLACKLIST = 3
const IGNORE_AND_BLACKLIST = 4
end

Base.@kwdef struct resp_memberJoinRequestEvent <: AbstractRequestCommand
    eventId::EventId
    fromId::FriendId
    groupId::GroupId
    operate::Int
    message::String
end

module BotInviteJoinGroupOperations
const ACCEPT = 0
const REFUSE = 1
end

Base.@kwdef struct resp_botInvitedJoinGroupRequestEvent <: AbstractRequestCommand
    eventId::EventId
    fromId::FriendId
    groupId::GroupId
    operate::Int
    message::String
end

StructTypes.@Struct struct CompatMessageChainResponse
    code::Int
    data::Vector{EventOrMessage}
end

abstract type AbstractGetMessageCommand <: AbstractCommand end
method(::AbstractGetMessageCommand) = CommandMethods.GET
response_type(::AbstractGetMessageCommand) = RESTful{Vector{EventOrMessage}}
response_type_compat(::AbstractGetMessageCommand) = CompatMessageChainResponse

struct fetchMessage <: AbstractGetMessageCommand
    count::Int
end

struct fetchLatestMessage <: AbstractGetMessageCommand
    count::Int
end

struct peekMessage <: AbstractGetMessageCommand
    count::Int
end

struct peekLatestMessage <: AbstractGetMessageCommand
    count::Int
end

struct countMessage <: AbstractGetMessageCommand end
response_type(::countMessage) = RESTful{Int}
response_type_compat(::countMessage) = RESTfulErrorMessage{Int}

end


text:
module LibCUFFT
using Compat

import Base.Sys: WORD_SIZE

include("cufft_h.jl")
import CUDArt.rt.cudaStream_t

if is_windows()
    const libcufft = Libdl.find_library([string("cufft", WORD_SIZE, "_75"),
      string("cufft", WORD_SIZE, "_70"), string("cufft", WORD_SIZE, "_65")],
      [joinpath(ENV["CUDA_PATH"], "bin")])
else
    const libcufft = Libdl.find_library(["libcufft", "cufft"], ["/usr/lib", "/usr/local/cuda"])
end

isempty(libcufft) && error("Cannot load libcufft")

function checkerror(code::cufftResult)
    if code == CUFFT_SUCCESS
        return nothing
    end
    println("An error was triggered:")
    Base.show_backtrace(STDOUT, backtrace())
    throw(string(Symbol(code)))
end

function cufftPlan1d(plan, nx, _type, batch)
  checkerror(ccall( (:cufftPlan1d, libcufft), cufftResult, (Ptr{cufftHandle}, Cint, cufftType, Cint), plan, nx, _type, batch))
end
function cufftPlan2d(plan, nx, ny, _type)
  checkerror(ccall( (:cufftPlan2d, libcufft), cufftResult, (Ptr{cufftHandle}, Cint, Cint, cufftType), plan, nx, ny, _type))
end
function cufftPlan3d(plan, nx, ny, nz, _type)
  checkerror(ccall( (:cufftPlan3d, libcufft), cufftResult, (Ptr{cufftHandle}, Cint, Cint, Cint, cufftType), plan, nx, ny, nz, _type))
end
function cufftPlanMany(plan, rank, n, inembed, istride, idist, onembed, ostride, odist, _type, batch)
  checkerror(ccall( (:cufftPlanMany, libcufft), cufftResult, (Ptr{cufftHandle}, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Ptr{Cint}, Cint, Cint, cufftType, Cint), plan, rank, n, inembed, istride, idist, onembed, ostride, odist, _type, batch))
end
function cufftDestroy(plan)
  checkerror(ccall( (:cufftDestroy, libcufft), cufftResult, (cufftHandle,), plan))
end
function cufftExecC2C(plan, idata, odata, direction)
  checkerror(ccall( (:cufftExecC2C, libcufft), cufftResult, (cufftHandle, Ptr{cufftComplex}, Ptr{cufftComplex}, Cint), plan, idata, odata, direction))
end
function cufftExecR2C(plan, idata, odata)
  checkerror(ccall( (:cufftExecR2C, libcufft), cufftResult, (cufftHandle, Ptr{cufftReal}, Ptr{cufftComplex}), plan, idata, odata))
end
function cufftExecC2R(plan, idata, odata)
  checkerror(ccall( (:cufftExecC2R, libcufft), cufftResult, (cufftHandle, Ptr{cufftComplex}, Ptr{cufftReal}), plan, idata, odata))
end
function cufftExecZ2Z(plan, idata, odata, direction)
  checkerror(ccall( (:cufftExecZ2Z, libcufft), cufftResult, (cufftHandle, Ptr{cufftDoubleComplex}, Ptr{cufftDoubleComplex}, Cint), plan, idata, odata, direction))
end
function cufftExecD2Z(plan, idata, odata)
  checkerror(ccall( (:cufftExecD2Z, libcufft), cufftResult, (cufftHandle, Ptr{cufftDoubleReal}, Ptr{cufftDoubleComplex}), plan, idata, odata))
end
function cufftExecZ2D(plan, idata, odata)
  checkerror(ccall( (:cufftExecZ2D, libcufft), cufftResult, (cufftHandle, Ptr{cufftDoubleComplex}, Ptr{cufftDoubleReal}), plan, idata, odata))
end
function cufftSetStream(plan, stream)
  checkerror(ccall( (:cufftSetStream, libcufft), cufftResult, (cufftHandle, cudaStream_t), plan, stream))
end
function cufftSetCompatibilityMode(plan, mode)
  checkerror(ccall( (:cufftSetCompatibilityMode, libcufft), cufftResult, (cufftHandle, cufftCompatibility), plan, mode))
end
function cufftGetVersion(version)
  checkerror(ccall( (:cufftGetVersion, libcufft), cufftResult, (Ptr{Cint},), version))
end

end


text:
function γn(xs::Vector{Float64}, s::Int64)
  n::Int64 = length(xs)
  sa::Int64 = abs(s)
  xbar::Float64 = mean(xs)
  acc::Float64 = 0.0
  for i in 1:n-sa
    @inbounds acc += (xs[i]-xbar)*(xs[i+sa]-xbar)
  end
  return acc/n
end

function estimateSV(xs::Vector{Float64}, b::Int64, ws::Vector{Float64})
  v::Float64 = γn(xs, 0)
  for s in 1:b
    v += 2*ws[s]*γn(xs, s)
  end
  return v
end

function _wsSimpleTruncation(b::Int64)
  return ones(Float64, b)
end

function _wsBlackmanTukey(b::Int64, a::Float64)
  ws::Vector{Float64} = Vector{Float64}(b)
  for k in 1:b
    ws[k] = 1 - 2*a + 2*a*cos(π*k/b)
  end
  return ws
end

function _wsTukeyHanning(b::Int64)
  return _wsBlackmanTukey(b, 0.25)
end

function _wsParzen(b::Int64, q::Int64)
  ws::Vector{Float64} = Vector{Float64}(undef, b)
  for k in 1:b
    ws[k] = 1 - (k/b)^q
  end
  return ws
end

function _wsModifiedBartlett(b::Int64)
  return _wsParzen(b, 1)
end

function estimateSV(xs::Vector{Float64}, b::Int64, name::Symbol)
  name == :ModifiedBartlett && return estimateSV(xs, b, _wsModifiedBartlett(b))
end

function estimateSV(xs::Vector{Float64}, name::Symbol)
  b::Int64 = floor(Int64, sqrt(length(xs)))
  name == :ModifiedBartlett && return estimateSV(xs, b, _wsModifiedBartlett(b))
end