text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module LivenessAnalysis

import ..DebugMsg
DebugMsg.init()

using CompilerTools
using CompilerTools.Helper
import CompilerTools.CFGs
using CompilerTools.LambdaHandling
using Core: Box, IntrinsicFunction

import Base.show

include("function-descriptions.jl")

export BlockLiveness, find_bb_for_statement, show

type Access
    sym
    read
end

"""
Liveness information for a TopLevelStatement in the CFG.
Contains a pointer to the corresponding CFG TopLevelStatement.
Contains def, use, live_in, and live_out for the current statement.
"""
type TopLevelStatement
    tls :: CFGs.TopLevelStatement

    def      :: Set{LHSVar}
    use      :: Set{LHSVar}
    live_in  :: Set{LHSVar}
    live_out :: Set{LHSVar}

    TopLevelStatement(t) = new(t, Set{LHSVar}(), Set{LHSVar}(), Set{LHSVar}(), Set{LHSVar}())
end

"""
Overload of Base.show to pretty-print a LivenessAnalysis.TopLevelStatement.
"""
function show(io::IO, tls::TopLevelStatement)
    print(io, "TLS ", tls.tls.index)

    print(io, " Def = (")
    for i in tls.def
      print(io, i, " ")
    end
    print(io, ") ")

    print(io, "Use = (")
    for i in tls.use
      print(io, i, " ")
    end
    print(io, ") ")

    print(io, "LiveIn = (")
    for i in tls.live_in
      print(io, i, " ")
    end
    print(io, ") ")

    print(io, "LiveOut = (")
    for i in tls.live_out
      print(io, i, " ")
    end
    print(io, ") ")

    println(io, "Expr: ", tls.tls.expr)
end

"""
Sometimes if new AST nodes are introduced then we need to ask for their def and use set as a whole
and then incorporate that into our liveness analysis directly.
"""
type AccessSummary
    def
    use
end

"""
Liveness information for a BasicBlock.
Contains a pointer to the corresponding CFG BasicBlock.
Contains def, use, live_in, and live_out for this basic block.
Contains an array of liveness information for the top level statements in this block.
"""
type BasicBlock
    cfgbb :: CFGs.BasicBlock

    def      :: Set{LHSVar}
    use      :: Set{LHSVar}
    live_in  :: Set{LHSVar}
    live_out :: Set{LHSVar}

    statements :: Array{TopLevelStatement,1}

    BasicBlock(bb) = new(bb, Set{LHSVar}(), Set{LHSVar}(), Set{LHSVar}(), Set{LHSVar}(), TopLevelStatement[])
end

"""
Overload of Base.show to pretty-print a LivenessAnalysis.BasicBlock.
"""
function show(io::IO, bb::BasicBlock)
    print(io, bb.cfgbb)
    print(io," Defs(")
    for j in bb.def
        print(io, " ", j)
    end
    print(io," ) Uses(")
    for j in bb.use
        print(io, " ", j)
    end
    print(io," ) LiveIn(")
    for j in bb.live_in
        print(io, " ", j)
    end
    print(io," ) LiveOut(")
    for j in bb.live_out
        print(io, " ", j)
    end
    println(io,")")

    tls = bb.statements
    if length(tls) == 0
        println(io,"Basic block without any statements.")
    end
    for j = 1:length(tls)
        if(DEBUG_LVL >= 5)
            print(io, "    ",tls[j].tls.index, "  ", tls[j].tls.expr)
            print(io,"  Defs(")
            for k in tls[j].def
                print(io, " ", k)
            end
            print(io," ) Uses(")
            for k in tls[j].use
                print(io, " ", k)
            end
            print(io," ) LiveIn(")
            for k in tls[j].live_in
                print(io, " ", k)
            end
            print(io," ) LiveOut(")
            for k in tls[j].live_out
                print(io, " ", k)
            end
            print(io," )")
            println(io)
        end
    end
end

"""
Called when AST traversal finds some Symbol or GenSym "sym" in a basic block "bb".
"read" is true if the symbol is being used and false if it is being defined.
"""
function add_access(bb, sym, read)
    @dprintln(3,"add_access ", sym, " ", read, " ", typeof(bb), " ")
    if bb == nothing    # If not in a basic block this is dead code so ignore.
        return nothing
    end

    assert(length(bb.statements) != 0)
    tls = bb.statements[end]    # Get the statement to which we will add access information.
    @dprintln(3, "tls = ", tls)
    write = !read

    # If this synbol is already in the statement then it is already in the basic block as a whole.
    # Therefore nothing would change if we did the rest of the function and we can shortcut quit here.
    if in(sym, tls.def)
        @dprintln(3, "sym already in tls.def")
        return nothing
    end

    # If the first use is a read then it goes in tls.use.
    # If there is a write after a read then it will be in tls.use and tls.def.
    # If there is a read after a write (which can happen for basic blocks) then we
    # ignore the read from the basic block perspective.

    # Handle access modifications at the statement level.
    if in(sym, tls.use)
        if write
            @dprintln(3, "sym already in tls.use so adding to def")
            # Handle write after read for a statement.
            push!(tls.def, sym)
        end
    elseif read
        if in(sym, tls.def)
            throw(string("Found a read after a write at the statement level in liveness analysis."))
        end
        @dprintln(3, "adding sym to tls.use ", sym)
        # Record the first use of this sym in this statement.
        push!(tls.use, sym)
    else # must be a write
        @dprintln(3, "adding sym to tls.def ", sym)
        # Record the first def of this sym in this statement.
        push!(tls.def, sym)
    end

    @dprintln(3, "tls after = ", tls)

    # Handle access modifications at the basic block level.
    if in(sym, bb.use)
        if write
            @dprintln(3, "sym already in bb.use so adding to def")
            # Handle write after read for a  basic block.
            push!(bb.def, sym)
        end
    elseif read
        if !in(sym, bb.def)
            @dprintln(3, "adding sym to bb.use ", sym)
            push!(bb.use, sym)
        end
    else # must be a write
        @dprintln(3, "adding sym to bb.def ", sym)
        push!(bb.def, sym)
    end

    @dprintln(4, "bb after = ", bb)

    nothing
end

"""
Holds the state during the AST traversal.
cfg = the control flow graph from the CFGs module.
map = our own map of CFG basic blocks to our own basic block information with liveness in it.
cur_bb = the current basic block in which we are processing AST nodes.
read = whether the sub-tree we are currently processing is being read or written.
ref_params = those arguments to the function that are passed by reference.
"""
type expr_state
    cfg :: CFGs.CFG
    map :: Dict{CFGs.BasicBlock, BasicBlock}
    cur_bb
    read
    ref_params :: Array{LHSVar, 1}
    params_not_modified :: Dict{Tuple{Any,Array{DataType,1}}, Array{Int64,1}} # Store function/signature mapping to an array whose entries corresponding to whether that argument passed to that function can be modified.
    params_not_modified_cb # a callback function to ask if some unknown function has unmodified args
    li :: Union{Void, LambdaVarInfo}

    function expr_state(cfg, no_mod, no_mod_cb)
        new(cfg, Dict{CFGs.BasicBlock, BasicBlock}(), nothing, true, LHSVar[], no_mod, no_mod_cb, nothing)
    end
end

"""
The main return type from LivenessAnalysis.
Contains a dictionary that maps CFG basic block to liveness basic blocks.
Also contains the corresponding CFG.
"""
type BlockLiveness
    basic_blocks :: Dict{CFGs.BasicBlock, BasicBlock}
    cfg :: CFGs.CFG

    function BlockLiveness(bb, cfg)
      new(bb, cfg)
    end
end

function getBasicBlockFromBlockNumber(block_num, bl :: BlockLiveness)
    cfg_bb = bl.cfg.basic_blocks[block_num]
    return bl.basic_blocks[cfg_bb]
end

"""
Helper function to say if the given LHSVar is in some basic block in "bl" for the given field (which is :def, :use, :live_in, or :live_out).
"""
function is_internal(x :: LHSVar, bl :: BlockLiveness, field)
    for entry in bl.basic_blocks
        if in(x, getfield(entry[2], field))
            return true
        end
    end

    return false
end

"""
Returns true if LHSVar x is part of the live_in field in some basic block in bl.
"""
function is_livein(x :: LHSVar, bl :: BlockLiveness)
    return is_internal(x, bl, :live_in)
end

"""
Returns true if LHSVar x is part of the live_out field in some basic block in bl.
"""
function is_liveout(x :: LHSVar, bl :: BlockLiveness)
    return is_internal(x, bl, :live_out)
end

"""
Returns true if LHSVar x is part of the def field in some basic block in bl.
"""
function is_def(x :: LHSVar, bl :: BlockLiveness)
    return is_internal(x, bl, :def)
end

"""
Returns true if LHSVar x is part of the use field in some basic block in bl.
"""
function is_use(x :: LHSVar, bl :: BlockLiveness)
    return is_internal(x, bl, :use)
end

"""
The live_in, live_out, def, and use routines are all effectively the same but just extract a different field name.
Here we extract this common behavior where x can be a liveness or CFG basic block or a liveness or CFG statement.
bl is BlockLiveness type as returned by a previous LivenessAnalysis.
field is the name of the field requested.
"""
function get_info_internal(x::Union{BasicBlock,TopLevelStatement}, bl::BlockLiveness, field)
    return getfield(x, field)
end

function get_info_internal(x::CFGs.BasicBlock, bl::BlockLiveness, field)
    bb = bl.basic_blocks[x]
    return getfield(bb, field)
end

function get_info_internal(x::CFGs.TopLevelStatement, bl::BlockLiveness, field)
    for i in bl.basic_blocks
        for j in i[2].statements
            if x == j.tls
                return getfield(j, field)
            end
        end
    end
    throw(string("Couldn't find liveness statement corresponding to cfg statement. "))
end

function get_info_internal(x::Any, bl::BlockLiveness, field)
    throw(string("get_info_internal called with non-BB and non-TopLevelStatement input."))
end

"""
Get the live_in information for "x" where x can be a liveness or CFG basic block or a liveness or CFG statement.
"""
function live_in(x, bl :: BlockLiveness)
    get_info_internal(x, bl, :live_in)
end

"""
Get the live_out information for "x" where x can be a liveness or CFG basic block or a liveness or CFG statement.
"""
function live_out(x, bl :: BlockLiveness)
    get_info_internal(x, bl, :live_out)
end

"""
Get the def information for "x" where x can be a liveness or CFG basic block or a liveness or CFG statement.
"""
function def(x, bl :: BlockLiveness)
    get_info_internal(x, bl, :def)
end

"""
Get the use information for "x" where x can be a liveness or CFG basic block or a liveness or CFG statement.
"""
function use(x, bl :: BlockLiveness)
    get_info_internal(x, bl, :use)
end

"""
Overload of Base.show to pretty-print BlockLiveness type.
"""
function show(io::IO, bl::BlockLiveness)
    println(io)
    body_order = CFGs.getBbBodyOrder(bl.cfg)
    for i = 1:length(body_order)
      cfgbb = bl.cfg.basic_blocks[body_order[i]]
      if !haskey(bl.basic_blocks, cfgbb)
        println(io,"Could not find LivenessAnalysis basic block for CFG basic block = ", cfgbb)
      else
        bb = bl.basic_blocks[bl.cfg.basic_blocks[body_order[i]]]
        show(io, bb)
        println(io)
      end
    end
end

"""
Query if the symbol in argument "x" is defined in live_info which can be a BasicBlock or TopLevelStatement.
"""
function isDef(x :: LHSVar, live_info)
  in(x, live_info.def)
end

"""
Search for a statement with the given top-level number in the liveness information.
Returns a LivenessAnalysis.TopLevelStatement having that top-level number or "nothing" if no such statement could be found.
"""
function find_top_number(top_number::Int, bl::BlockLiveness)
  # Liveness information stored in blocks so scan each block.
  for bb in bl.basic_blocks
    stmts = bb[2].statements
    # Scan each statement in this block for a matching statement number.
    for j = 1:length(stmts)
      if stmts[j].tls.index == top_number
        return stmts[j]
      end
    end
  end
  nothing
end

"""
Search for a basic block containing a statement with the given top-level number in the liveness information.
Returns a basic block label of a block having that top-level number or "nothing" if no such statement could be found.
"""
# Search for a statement with the given number in the liveness information.
function find_bb_for_statement(top_number::Int, bl::BlockLiveness)
  # Liveness information stored in blocks so scan each block.
  for bb in bl.basic_blocks
    stmts = bb[2].statements
    # Scan each statement in this block for a matching statement number.
    for j = 1:length(stmts)
      if stmts[j].tls.index == top_number
        return bb[1].label
      end
    end
  end

  @dprintln(3,"Didn't find statement top_number in basic_blocks.")
  nothing
end

"""
Clear the live_in and live_out data corresponding to all basic blocks and statements and then recompute liveness information.
"""
function recompute_live_ranges(state, dfn; array_params_live_out=true)
    for bb in state.basic_blocks
        empty!(bb.live_in)
        empty!(bb.live_out)
        for s in bb.statements
          empty!(s.live_in)
          empty!(s.live_out)
        end
    end

    compute_live_ranges(state, dfn, array_params_live_out)

    nothing
end

"""
Compute the live_in and live_out information for each basic block and statement.
"""
function compute_live_ranges(state :: expr_state, dfn, array_params_live_out)
    found_change = true
    bbs = state.cfg.basic_blocks

    # Compute the inter-block live ranges first.
    while found_change
        # Iterate until quiescence.
        found_change = false

        @dprintln(4, "Starting iteration of while loop in compute_live_ranges.")

        # For each basic block in reverse depth-first order.
        # This order optimizes the convergence time but isn't necessary for correctness.
        for i = length(dfn):-1:1
            bb_index = dfn[i]
            bb = state.map[bbs[bb_index]]   # get the basic block
            @dprintln(4, "Working on block ", bb_index)

            # add escaping variables to accum
            accum = state.li == nothing ? Set{LHSVar}() : Set{LHSVar}(getEscapingVariablesAsLHSVar(state.li))

            if bb_index == -2
              # Special case for final block.
              @dprintln(3,"Final block live_out = ", accum)
              # Nothing is live out of the final block.
              if array_params_live_out
                  accum = Set{LHSVar}(state.ref_params)
              end
            else
              # The live_out of any non-final block is the union of the live_in of every successor block.
              for j in bb.cfgbb.succs
                  @dprintln(4, bb_index, " adding from ", j.label, " live_in = ", state.map[j].live_in)
                  accum = union(accum, state.map[j].live_in)
              end
            end

            bb.live_out = accum
            @dprintln(4, bb_index, " bb.live_out = ", bb.live_out, " bb.use = ", bb.use, " bb.def = ", bb.def)

            old_size = length(bb.live_in)

            # The live_in to this block is the union of things used in this block
            # with the set of things used by successors and not defined in this block.
            # Note that for basic blocks, we do not create a "use" if the first "use"
            # is after a "def".
            bb.live_in = union(bb.use, setdiff(bb.live_out, bb.def))
            @dprintln(4, bb_index, " bb.live_in = ", bb.live_in)

            new_size = length(bb.live_in)
            if new_size != old_size
                found_change = true
            end
        end
    end

    # Compute the intra-block live ranges using the inter-block live ranges.
    # For each basic block.
    for i = 1:length(dfn)
        bb = state.map[bbs[dfn[i]]]

        tls = bb.statements

        # The lives at the end of the last statement of this block is the inter-block live_out.
        cur_live_out = bb.live_out

        # Go through the statements in reverse order.
        for j = length(tls):-1:1
            tls[j].live_out = cur_live_out;
            tls[j].live_in  = union(tls[j].use, setdiff(tls[j].live_out, tls[j].def))
            cur_live_out    = tls[j].live_in
        end
    end

    nothing
end

"""
Dump a bunch of debugging information about BlockLiveness.
"""
function dump_bb(bl :: BlockLiveness)
    if DEBUG_LVL >= 4
      f = open("bbs.dot","w")
      println(f, "digraph bbs {")
    end

    body_order = CFGs.getBbBodyOrder(bl.cfg)
    if DEBUG_LVL >= 3
      println("body_order = ", body_order)
    end

    for i = 1:length(body_order)
        bb = bl.basic_blocks[bl.cfg.basic_blocks[body_order[i]]]
        @dprint(2,bb)

        if DEBUG_LVL >= 4
            for j in bb.cfgbb.succs
                println(f, bb.cfgbb.label, " -> ", j.label, ";")
            end
        end
    end

    if DEBUG_LVL >= 4
      println(f, "}")
      close(f)
    end
end

"""
Walk through an array of expressions.
Just recursively call from_expr for each expression in the array.
"""
function from_exprs(ast, depth :: Int64, state :: expr_state, callback :: Function, cbdata :: ANY)
    # sequence of expressions
    # ast = [ expr, ... ]
    local len = length(ast)
    for i = 1:len
      @dprintln(2,"Processing ast #",i," depth=",depth)
      @dprintln(3,"ast[", i, "] = ", ast[i])
      from_expr(ast[i], depth, state, callback, cbdata)
    end
    nothing
end

"""
Walk through an assignment expression.
"""
function from_assignment(ast :: Array{Any,1}, depth :: Int64, state :: expr_state, callback :: Function, cbdata :: ANY)
    # :(=) assignment
    # ast = [ ... ]
    assert(length(ast) == 2)
    local lhs = ast[1]
    local rhs = ast[2]
    @dprintln(3,"liveness from_assignment lhs = ", lhs, " rhs = ", rhs)
    # Process the right-hand side of the assignment unless it is a lambda.
    if isa(rhs, Expr) && rhs.head == :lambda
      # skip handling rhs lambdas
    else
      from_expr(rhs, depth, state, callback, cbdata)
    end
    @dprintln(3,"liveness from_assignment handling lhs")
    # Handle the left-hand side of the assignment which is being written.
    read_save = state.read
    state.read = false
    from_expr(lhs, depth, state, callback, cbdata)
    state.read = read_save
    @dprintln(3,"liveness from_assignment done handling lhs")
end

"""
Add an entry the dictionary of which arguments can be modified by which functions.
"""
function addUnmodifiedParams(func, signature :: Array{DataType,1}, unmodifieds, state :: expr_state)
    state.params_not_modified[(func, signature)] = unmodifieds
end

"""
Get the type of some AST node.
"""
function typeOfOpr(x::Expr, li :: LambdaVarInfo)
    @dprintln(3,"starting typeOfOpr, type = Expr")
    return typeOfOpr_fixType(x.typ)
end

function typeOfOpr(x::Symbol, li :: LambdaVarInfo)
    @dprintln(3,"starting typeOfOpr, type = LHSVar, x = ", x)
    return isLocalVariable(x, li) ? typeOfOpr_fixType(getType(x, li)) : Void
end

function typeOfOpr(x::RHSVar, li :: LambdaVarInfo)
    @dprintln(3,"starting typeOfOpr, x = ", x)
    typ1 = getType(toLHSVar(x), li)
    @dprintln(3,"typ1 = ", typ1)
    if isa(x, TypedVar) && x.typ != typ1
        @dprintln(2, "typeOfOpr x.typ and lambda type different")
        @dprintln(2, "x = ", x, " typ1 = ", typ1)
        @dprintln(2, "li = ", li)
        if (x.typ <: typ1) || (typ1 === Box)
            typ1 = x.typ
        elseif (typ1 <: x.typ) || (x.typ === Box)
        else
            throw(string("typeOf Opr ", x, " is incompatible with its type in lambda ", typ1))
        end
    end
    assert(isa(typ1, Type))
    return typeOfOpr_fixType(typ1)
end

function typeOfOpr(x::GlobalRef, li :: LambdaVarInfo)
    @dprintln(3,"starting typeOfOpr, type = GlobalRef ", x)
    if x.name == :end
        return Int64
    else
        return typeOfOpr_fixType(typeof(eval(x)))
    end
end

function typeOfOpr(x::SimpleVector, li :: LambdaVarInfo)
    @dprintln(3,"starting typeOfOpr, type = SimpleVector")
    svec_types = [ typeOfOpr(x[i], li) for i = 1:length(x) ]
    return Tuple{svec_types...}
end

function typeOfOpr(x::Any, li :: LambdaVarInfo)
    @dprintln(3,"starting typeOfOpr, type = ",typeof(x))
    return typeOfOpr_fixType(typeof(x))
end

function typeOfOpr_fixType(ret::DataType)
    return ret
end

function typeOfOpr_fixType(ret::Type)
    return ret
end

if VERSION >= v"0.6.0-pre"
function typeOfOpr_fixType(ret::Core.TypeofBottom)
    return ret
end
end

function typeOfOpr_fixType(ret::Union)
    return Tuple{ret.types...}
end

"""
Returns true if a parameter is passed by reference.
isbits types are not passed by ref but everything else is (is this always true..any exceptions?)
"""
function isPassedByRef(x, state :: expr_state)
    if isa(x, Tuple)
        return true
    elseif isbits(x)
        return false
    else
        return true
    end
end

function showNoModDict(dict)
    for i in dict
      try
        @dprintln(4, "(", i[1][1], ",", i[1][2], ") => ", i[2])
      catch
        targs = i[1][2]
        assert(isa(targs, Tuple))
        println("EXCEPTION: type = ", typeof(targs))
        for j = 1:length(targs)
           println(j, " = ", typeof(targs[j]))
           println(targs[j])
        end
      end
    end
end

# If true, will assume that functions without "!" don't update their arguments.
use_inplace_naming_convention = false
function set_use_inplace_naming_convention()
    global use_inplace_naming_convention = true
end

wellknown_all_unmodified = Set{Any}()
wellknown_only_first_modified = Set{Any}()

function __init__()
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(./)), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(.*)), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(.+)), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(.-)), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(/)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(*)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(+)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(-)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(<=)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(<)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(>=)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:(>)),  force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:size), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:maximum), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:minimum), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:max), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:min), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:mean), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:ctranspose), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base.LinAlg,:norm), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:Ac_mul_B), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:Ac_mul_Bc), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:box), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:arraylen), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:arraysize), force = true))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base,:transpose), force = true))
#  push!(wellknown_all_unmodified, eval(TopNode(:(!))))
  push!(wellknown_all_unmodified, Base.resolve(GlobalRef(Base.LinAlg,:typeof), force = true))

    push!(wellknown_only_first_modified, Base.resolve(GlobalRef(Base.LinAlg,:gemm_wrapper!), force = true))
    push!(wellknown_only_first_modified, Base.resolve(GlobalRef(Base.LinAlg,:gemv!), force = true))
    push!(wellknown_only_first_modified, Base.resolve(GlobalRef(Base,:transpose!), force = true))
end

"""
For a given function and signature, return which parameters can be modified by the function.
If we have cached this information previously then return that, else cache the information for some
well-known functions or default to presuming that all arguments could be modified.
"""
function getUnmodifiedArgs(func :: ANY, args, arg_type_tuple :: Array{DataType,1}, state :: expr_state)
  @dprintln(3,"getUnmodifiedArgs func = ", func, " type = ", typeof(func))
  @dprintln(3,"getUnmodifiedArgs args = ", args)
  @dprintln(3,"getUnmodifiedArgs arg_type_tuple = ", arg_type_tuple)
  @dprintln(3,"getUnmodifiedArgs len(args) = ", length(arg_type_tuple))
  showNoModDict(state.params_not_modified)

  default_result = Int64[(isPassedByRef(x, state) ? 0 : 1) for x in arg_type_tuple]
  if length(default_result) == 0
    return default_result
  end

  if typeof(func) == GlobalRef
    func = Base.resolve(func, force=true)
    @dprintln(3,"getUnmodifiedArgs func = ", func, " type = ", typeof(func))
  elseif isTopNode(func)
    func = eval(func)
    @dprintln(3,"getUnmodifiedArgs func = ", func, " type = ", typeof(func))
  elseif typeof(func) == Expr
    return default_result
  end

  # We are seeing Symbol's getting here as well due to incomplete name resolution.  Once this is
  # fixed then maybe we re-enable this assertion as a sanity check.
#  assert(typeof(func) == Function || typeof(func) == IntrinsicFunction)

  fs = (func, arg_type_tuple)
  if haskey(state.params_not_modified, fs)
    res = state.params_not_modified[fs]
    assert(length(res) == length(args))
    @dprintln(3,"function already in params_not_modified so returning previously computed value")
    return res
  end

  for i in state.params_not_modified
    (f1, t1) = i[1]
    @dprintln(3,"f1 = ", f1, " t1 = ", t1, " f1type = ", typeof(f1), " len(t1) = ", length(t1))
    if func == f1
      @dprintln(3,"function matches")
      if Tuple{arg_type_tuple...} <: Tuple{t1...}
        res = i[2]
        assert(length(res) == length(args))
        addUnmodifiedParams(func, arg_type_tuple, res, state)
        @dprintln(3,"exact match not found but sub-type match found")
        return res
      end
    end
  end

  if in(func, wellknown_all_unmodified)
    @dprintln(3,"Well-known function known not to modify args.")
    addUnmodifiedParams(func, arg_type_tuple, ones(Int64, length(args)), state)
  elseif in(func, wellknown_only_first_modified)
    # only first argument is def
    local params_res = ones(Int64, length(args))
    params_res[1] = 0
    addUnmodifiedParams(func, arg_type_tuple, params_res, state)
  else
    if isBaseFunc(func, :tuple)
      @dprintln(3,"Detected tuple in getUnmodifiedArgs so returning that no arguments are modified.")
      addUnmodifiedParams(func, arg_type_tuple, [1 for x in arg_type_tuple], state)
      return state.params_not_modified[fs]
    end

    @dprintln(3,"is func generic => ", isa(func,Function))
    if use_inplace_naming_convention && isa(func,Function) && !in('!', string(Base.function_name(func)))
      @dprintln(3,"using naming convention that function has no ! so it doesn't modify anything in place.")
      addUnmodifiedParams(func, arg_type_tuple, [1 for x in arg_type_tuple], state)
    else
      if state.params_not_modified_cb != nothing && typeof(func) == GlobalRef
        res = state.params_not_modified_cb(func, arg_type_tuple)
        if res == nothing
          @dprintln(3,"fallback to args passed by ref as modified.")
          addUnmodifiedParams(func, arg_type_tuple, default_result, state)
        else
          @dprintln(3,"callback returned unmodified information.")
          addUnmodifiedParams(func, arg_type_tuple, res, state)
        end
      else
        @dprintln(3,"fallback to args passed by ref as modified.")
        addUnmodifiedParams(func, arg_type_tuple, default_result, state)
      end
    end
  end

  return state.params_not_modified[fs]
end

"""
Walk through a call expression.
"""
function from_foreigncall(ast :: Expr , depth :: Int64, state :: expr_state, callback :: Function, cbdata :: ANY)
  local fun  = getCallFunction(ast)
  local args = getCallArguments(ast)
  @dprintln(2,"from_foreigncall fun = ", fun, " typeof fun = ", typeof(fun), " args = ", args)

  # Form the signature of the call in an array.
  arg_type_array = DataType[]
  # For each argument to the call.
  for i = 1:length(args)
    @dprintln(3, "arg ", i, " = ", args[i], " typeof arg = ", typeof(args[i]))
    # Make sure the type of the argument resolves to a DataType.
    too = typeOfOpr(args[i], state.li)
    if !isa(too, DataType)
      @dprintln(0, "arg type = ", too, " tootype = ", typeof(too))
    end
    push!(arg_type_array, typeOfOpr(args[i], state.li))
  end
  @dprintln(3, "arg_type_array = ", arg_type_array)
  #arg_type_tuple = Tuple{arg_type_array...}
  # See which arguments to the function can be modified by the function.

  # We can do better/tighter analysis if we know which arguments to some call are modified.
  # bits type arguments are never modified.
  # non-bits types argument may be modified unless we know otherwise.
  # We prime a dictionary with some common, well-known functions and record that their non-bits type
  # arguments are not modified.
  # getUnmodifiedArgs returns an array to indicate which args to the call may be modified.
  unmodified_args = getUnmodifiedArgs(fun, args, arg_type_array, state)
  assert(length(unmodified_args) == length(args))
  @dprintln(3,"unmodified_args = ", unmodified_args)

  # For each argument.
  for i = 1:length(args)
    argtyp = typeOfOpr(args[i], state.li)
    @dprintln(2,"cur arg = ", args[i], " type = ", argtyp, " state.read = ", state.read)
    read_cache = state.read

    state.read = true
    # We can always potentially read first.
    from_expr(args[i], depth+1, state, callback, cbdata)
    if unmodified_args[i] == 0
      # The argument could be modified so treat it as a "def".
      state.read = false
      from_expr(args[i], depth+1, state, callback, cbdata)
    end
    state.read = read_cache
  end
end

"""
Walk through a call expression.
"""
function from_call(ast :: Expr , depth :: Int64, state :: expr_state, callback :: Function, cbdata :: ANY)
  local fun  = getCallFunction(ast)
  local args = getCallArguments(ast)
  @dprintln(2,"from_call fun = ", fun, " typeof fun = ", typeof(fun), " args = ", args)

  # Form the signature of the call in an array.
  arg_type_array = DataType[]
  # For each argument to the call.
  for i = 1:length(args)
    @dprintln(3, "arg ", i, " = ", args[i], " typeof arg = ", typeof(args[i]))
    # Make sure the type of the argument resolves to a DataType.
    too = typeOfOpr(args[i], state.li)
    if !isa(too, DataType)
      @dprintln(0, "arg type = ", too, " tootype = ", typeof(too))
    end
    push!(arg_type_array, too)
  end
  @dprintln(3, "arg_type_array = ", arg_type_array)
  #arg_type_tuple = Tuple{arg_type_array...}
  # See which arguments to the function can be modified by the function.

  # We can do better/tighter analysis if we know which arguments to some call are modified.
  # bits type arguments are never modified.
  # non-bits types argument may be modified unless we know otherwise.
  # We prime a dictionary with some common, well-known functions and record that their non-bits type
  # arguments are not modified.
  # getUnmodifiedArgs returns an array to indicate which args to the call may be modified.
  unmodified_args = getUnmodifiedArgs(fun, args, arg_type_array, state)
  assert(length(unmodified_args) == length(args))
  @dprintln(3,"unmodified_args = ", unmodified_args)

  # symbols don't need to be translated
  if typeof(fun) != Symbol
      from_expr(fun, depth, state, callback, cbdata)
  end

  # For each argument.
  for i = 1:length(args)
    argtyp = typeOfOpr(args[i], state.li)
    @dprintln(2,"cur arg = ", args[i], " type = ", argtyp, " state.read = ", state.read)
    read_cache = state.read

    state.read = true
    # We can always potentially read first.
    from_expr(args[i], depth+1, state, callback, cbdata)
    if unmodified_args[i] == 0
      # The argument could be modified so treat it as a "def".
      state.read = false
      from_expr(args[i], depth+1, state, callback, cbdata)
    end
    state.read = read_cache
  end
end

"""
The default callback that processes no non-standard Julia AST nodes.
"""
function not_handled(a,b)
  nothing
end

"""
Count the number of times that the symbol in "s" is defined in all the basic blocks.
"""
function countSymbolDefs(s, lives)
  @dprintln(3,"countSymbolDefs: ", s)
  count = 0
  for (j,bb) in lives.basic_blocks
    @dprintln(3,"Examining block ", j.label)
    for stmt in bb.statements
      if in(s, stmt.def)
          @dprintln(3, "Found symbol defined in block ", j.label, " in statement: ", stmt)
          count += 1
      end
    end
  end
  return count
end

"""
ENTRY point to liveness analysis.
You must pass a lambda either as two parts, LambdaVarInfo and body, or as Expr (with :lambda head).
If you have non-standard AST nodes, you may pass a callback that will be given a chance to process the non-standard node first.
"""
function from_lambda(LambdaVarInfo :: LambdaVarInfo, body::ANY, callback=not_handled, cbdata :: ANY = nothing, no_mod=Dict{Tuple{Any,Array{DataType,1}}, Array{Int64,1}}(); no_mod_cb = nothing, array_params_live_out=true)
  if isa(body, Array)
    body = CompilerTools.LambdaHandling.getBody(body, getReturnType(LambdaVarInfo))
  end
  @assert (isa(body, Expr)) "Expect body to be Expr, but got " * string(typeof(body))
  @dprintln(3, "from_lambda body = ", body)
  #@dprintln(3,"liveness from_expr no_mod = ", no_mod)
  cfg = CFGs.from_lambda(body)      # Create the CFG from this lambda Expr.
  live_res = expr_state(cfg, no_mod, no_mod_cb)
  # Just to process the lambda and extract what the ref_params are.
  live_res.li = LambdaVarInfo
  live_res.ref_params = CompilerTools.LambdaHandling.getRefParams(live_res.li)
  # Process the body of the function via the CFG.
  res = fromCFG(live_res, cfg, callback, cbdata, array_params_live_out=array_params_live_out)
  input_params = CompilerTools.LambdaHandling.getInputParameters(live_res.li)
  ip_lhsvar = [CompilerTools.LambdaHandling.toLHSVar(x, live_res.li) for x in input_params]
  escape_lhsvar = CompilerTools.LambdaHandling.getEscapingVariablesAsLHSVar(live_res.li)
  start_live_in = getBasicBlockFromBlockNumber(CompilerTools.CFGs.CFG_ENTRY_BLOCK, res).live_in
  combined = union(ip_lhsvar, escape_lhsvar)
  use_before_init = setdiff(start_live_in, combined)
  @dprintln(2, "input_params = ", input_params)
  @dprintln(2, "ip_lhsvar = ", ip_lhsvar)
  @dprintln(2, "escape_lhsvar = ", escape_lhsvar)
  @dprintln(2, "combined = ", combined)
  @dprintln(2, "start_live_in = ", start_live_in)
  @dprintln(2, "use_before_init = ", use_before_init)
  if !isempty(use_before_init)
     @dprintln(1, "WARNING: use of variables before initialization = ", use_before_init)
     @dprintln(1, "LambdaVarInfo = ", live_res.li)
#     throw(string("WARNING: use of variables before initialization = ", use_before_init))
  end
  res
end

if VERSION >= v"0.6.0-pre"
import ..LambdaHandling.LambdaInfo
end

function from_lambda(lambda::Union{Expr,LambdaInfo}, callback=not_handled, cbdata :: ANY = nothing, no_mod=Dict{Tuple{Any,Array{DataType,1}}, Array{Int64,1}}(); no_mod_cb = nothing, array_params_live_out=true)
    linfo, body = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(lambda)
    from_lambda(linfo, body, callback, cbdata, no_mod; no_mod_cb = no_mod_cb, array_params_live_out = array_params_live_out)
end

"""
Extract liveness information from the CFG.
"""
function fromCFG(live_res, cfg :: CFGs.CFG, callback :: Function, cbdata :: ANY; array_params_live_out=true)
  @dprintln(2,"fromCFG")
  CFGs.dump_bb(cfg)   # Dump debugging information if set_debug_level is high enough.

  # For each basic block.
  for bb in cfg.basic_blocks
    live_res.map[bb[2]] = BasicBlock(bb[2])
    live_res.cur_bb = live_res.map[bb[2]]

    @dprintln(4,"fromCFG cur_bb = ", live_res.cur_bb)
    # For each statement in each block.
    for i = 1:length(bb[2].statements)
       cur_stmt = bb[2].statements[i]
       # Add this statement to our list of statements in the current LivenessAnalysis.BasicBlock.
       push!(live_res.cur_bb.statements, TopLevelStatement(cur_stmt))
       @dprintln(3,"fromCFG stmt = ", cur_stmt.expr)
       # Process the statement looking for def and use.
       from_expr(cur_stmt.expr, 1, live_res, callback, cbdata)
    end
  end

  @dprintln(4,"fromCFG live_res before live ranges = ", live_res)
  # Compute live_in and live_out for basic blocks and statements.
  compute_live_ranges(live_res, cfg.depth_first_numbering, array_params_live_out)
  @dprintln(2,"Dumping basic block info from_expr.")
  ret = BlockLiveness(live_res.map, cfg)
  dump_bb(ret)
  return ret
end

"""
Process a return Expr node which is just a recursive processing of all of its args.
"""
function from_return(args, depth :: Int64, state :: expr_state, callback :: Function, cbdata :: ANY)
    @dprintln(2,"Expr return: ")
    from_exprs(args, depth, state, callback, cbdata)
    nothing
end

"""
Process a gotoifnot which is just a recursive processing of its first arg which is the conditional.
"""
function from_if(args, depth :: Int64, state :: expr_state, callback :: Function, cbdata :: ANY)
    # The structure of the if node is an array of length 2.
    assert(length(args) == 2)
    # The first index is the conditional.
    if_clause  = args[1]

    # Process the conditional as part of the current basic block.
    from_expr(if_clause, depth, state, callback, cbdata)
    nothing
end

"""
Generic routine for how to walk most AST node types.
"""
function from_expr(ast::LambdaInfo,
                   depth::Int64,
                   state::expr_state,
                   callback::Function,
                   cbdata::ANY)
    # skip processing LambdaInfo
    @dprintln(3, "LivenessAnalysis: skip processing nested LambdaInfo")
    return nothing
end

function from_expr(ast::ANY,
                   depth::Int64,
                   state::expr_state,
                   callback::Function,
                   cbdata::ANY)
    @dprintln(2,"from_expr ast = ", ast, " depth = ",depth," ", " asttyp = ", typeof(ast), " state.read = ", state.read)

    handled = callback(ast, cbdata)
    @dprintln(3,"callback handled ", handled)
    if handled != nothing
        if length(handled) > 0
            @dprintln(3,"Processing expression from callback for ", ast)
            @dprintln(3,handled)
            from_exprs(handled, depth+1, state, callback, cbdata)
            @dprintln(3,"Done processing expression from callback.")
        end
        return nothing
    end

    from_expr_helper(ast, depth, state, callback, cbdata)
end

function from_expr_helper(ast::Union{Tuple,SimpleVector},
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    # N.B. This also handles the empty tuple correctly.

    for i = 1:length(ast)
        from_expr(ast[i], depth, state, callback, cbdata)
    end
end

function from_expr_helper(ast::Expr,
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    # If you have an assignment with an Expr as its left-hand side then you get here with "read = false"
    # but that doesn't  mean the whole Expr is written.  In fact, none of it is written so we set read
    # back to true and then restore in the incoming read value at the end.
    read_save  = state.read
    state.read = true

    @dprint(2,"Expr ")
    local head = ast.head
    local args = ast.args
    local typ  = ast.typ
    @dprintln(2,head, " ", args)
    if head == :lambda
        # from_lambda(ast, depth, state, callback, cbdata)
        @dprintln(0, "nested :lambda skipped")
    elseif head == :body
        @dprintln(0,":body found in from_expr")
        throw(string(":body found in from_expr"))
    elseif head == :(=)
        from_assignment(args, depth, state, callback, cbdata)
    elseif head == :return
        from_return(args, depth, state, callback, cbdata)
    elseif head == :invoke || head == :call || head == :call1
        from_call(ast, depth, state, callback, cbdata)
    elseif head == :foreigncall
        from_foreigncall(ast, depth, state, callback, cbdata)
    elseif head == :gotoifnot
        from_if(args,depth,state, callback, cbdata)
    elseif head == :line || head == :inbounds || head == :boundscheck || head == :meta || head == :type_goto || head == :static_parameter
        # Intentionally do nothing.
    elseif head == :copyast
        @dprintln(2,"copyast type")
        # Intentionally do nothing.
    elseif head == :alloc
        from_exprs(args[2], depth+1, state, callback, cbdata)
    elseif head == :assert || head == :select || head == :ranges || head == :range || head == :static_typeof ||
        head == :tomask || head == :arraysize || head == :copy || head == :new ||
        head == :tuple || head == :getindex || head == :quote || head == Symbol("'") ||
        head == :(&) || head == :(~) || head == :(|) || head == :($) || head == :(>>>) || head == :(>>) || head == :(<<)
        from_exprs(args, depth+1, state, callback, cbdata)
    elseif head == :(.)
        # skip handling fields of a type
        # ISSUE: will this cause precision issue, or correctness issue? I guess it is precision?
    elseif head == :llvmcall
        # Intentionally do nothing.
    elseif head == :simdloop
        # Intentionally do nothing.
    else
        throw(string("from_expr: unknown Expr head :", head))
    end
    state.read = read_save
end

function from_expr_helper(ast::LabelNode,
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    assert(false)
    # from_label(ast.label, state, callback, cbdata)
end

function from_expr_helper(ast::Union{GotoNode,LineNumberNode,TopNode,Module},
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    # Intentionally do nothing.
end

function from_expr_helper(ast::Union{Symbol,RHSVar},
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    # addStatement(top_level, state, ast)
    @dprintln(2, "ast = ", ast, "::", typeof(ast))
    if isVariableDefined(ast, state.li)
        v = toLHSVar(ast, state.li)
        add_access(state.cur_bb, v, state.read)
    end
end

function from_expr_helper(ast::Union{DataType,AbstractString,NewvarNode,Void},
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    # addStatement(top_level, state, ast)
end

function from_expr_helper(ast::QuoteNode,
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    # addStatement(top_level, state, ast)
    local value = ast.value
    #TODO: fields: value
    @dprintln(3,"QuoteNode type ",typeof(value))
end

function from_expr_helper(ast::AccessSummary,
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    @dprintln(3, "Incorporating AccessSummary")
    for i in ast.use
        add_access(state.cur_bb, i, true)
    end
    for i in ast.def
        add_access(state.cur_bb, i, false)
    end
end

function from_expr_helper(ast::ANY,
                          depth::Int64,
                          state::expr_state,
                          callback::Function,
                          cbdata::ANY)
    asttyp = typeof(ast)

    if isdefined(:GetfieldNode) && asttyp == GetfieldNode  # GetfieldNode = value + name
        # addStatement(top_level, state, ast)
        @dprintln(3,"GetfieldNode type ",typeof(ast.value), " ", ast)
    elseif isdefined(:GlobalRef) && asttyp == GlobalRef
        # addStatement(top_level, state, ast)
        @dprintln(3,"GlobalRef type ",typeof(ast.mod), " ", ast)  # GlobalRef = mod + name
    elseif isbits(asttyp)
        # addStatement(top_level, state, ast)
        # skip
    #elseif asttyp.name == Array.name
    elseif asttyp <: Array
        # addStatement(top_level, state, ast)
        @dprintln(3,"Handling case of AST node that is an array. ast = ", ast, " typeof(ast) = ", asttyp)
        for i = 1:length(ast)
            from_expr(ast[i], depth, state, callback, cbdata)
        end
    elseif asttyp == Symbol
        # skip
    else
        throw(string("from_expr: unknown AST type :", asttyp, " ", ast))
    end
end

end


text:
import Manifolds: is_decorator_manifold

struct NotImplementedOperation <: AbstractGroupOperation end

struct NotImplementedManifold <: Manifold end

struct NotImplementedGroupDecorator{M} <: Manifold
    manifold::M
end
is_decorator_manifold(::NotImplementedGroupDecorator) = Val(true)

"""
    test_group(
        G,
        g_pts::AbstractVector,
        v_pts::AbstractVector = [];
        atol = 1e-10,
        test_mutating = true,
        test_diff = false,
        diff_convs = [(), (LeftAction(),), (RightAction(),)],
    )

Tests general properties of the group `G`, given at least three different points
elements of it (contained in `g_pts`).
Optionally, specify `test_diff` to test differentials of translation, using `v_pts`, which
must contain at least one tangent vector at `g_pts[1]`, and the direction conventions
specified in `diff_convs`.
"""
function test_group(
    G,
    g_pts::AbstractVector,
    v_pts::AbstractVector = [];
    atol = 1e-10,
    test_mutating = true,
    test_diff = false,
    diff_convs = [(), (LeftAction(),), (RightAction(),)],
)
    e = Identity(G)

    @testset "Basic group properties" begin
        @testset "Closed" begin
            for g1 in g_pts, g2 in g_pts
                g3 = compose(G, g1, g2)
                @test is_manifold_point(G, g3, true; atol = atol)
            end
        end

        @testset "Associative" begin
            g12_3 = compose(G, compose(G, g_pts[1], g_pts[2]), g_pts[3])
            g1_23 = compose(G, g_pts[1], compose(G, g_pts[2], g_pts[3]))
            @test isapprox(G, g12_3, g1_23; atol = atol)

            test_mutating && @testset "mutating" begin
                g12, g23, g12_3, g1_23 = allocate.(repeat([g_pts[1]], 4))
                @test compose!(G, g12, g_pts[1], g_pts[2]) === g12
                @test compose!(G, g23, g_pts[2], g_pts[3]) === g23
                @test compose!(G, g12_3, g12, g_pts[3]) === g12_3
                @test compose!(G, g1_23, g_pts[1], g23) === g1_23
                @test isapprox(G, g12_3, g1_23; atol = atol)
            end
        end

        @testset "Identity" begin
            @test isapprox(G, e, e)
            @test identity(G, e) === e
            @test compose(G, e, e) === e
            @test copyto!(e, e) === e

            for g in g_pts
                @test isapprox(G, compose(G, g, e), g)
                @test isapprox(G, compose(G, e, g), g)

                ge = identity(G, g)
                @test isapprox(G, compose(G, g, ge), g)
                @test isapprox(G, compose(G, ge, g), g)
            end

            test_mutating && @testset "mutating" begin
                for g in g_pts
                    h = allocate(g)
                    @test compose!(G, h, g, e) === h
                    @test isapprox(G, h, g)
                    h = allocate(g)
                    @test compose!(G, h, e, g) === h
                    @test isapprox(G, h, g)

                    ge = allocate(g)
                    @test identity!(G, ge, e) === ge
                    @test isapprox(G, compose(G, g, ge), g)
                    @test isapprox(G, compose(G, ge, g), g)

                    ge = allocate(g)
                    @test compose!(G, ge, e, e) === ge
                    @test isapprox(G, ge, e)
                end
            end
        end

        @testset "Inverse" begin
            for g in g_pts
                ginv = inv(G, g)
                @test isapprox(G, compose(G, g, ginv), e; atol = atol)
                @test isapprox(G, compose(G, ginv, g), e; atol = atol)
                @test isapprox(G, e, compose(G, g, ginv); atol = atol)
                @test isapprox(G, e, compose(G, ginv, g); atol = atol)
                @test inv(G, e) === e

                test_mutating && @testset "mutating" begin
                    ginv = allocate(g)
                    @test inv!(G, ginv, g) === ginv
                    @test isapprox(G, compose(G, g, ginv), e; atol = atol)
                    @test isapprox(G, compose(G, ginv, g), e; atol = atol)

                    @test inv(G, e) === e
                    geinv = allocate(g)
                    @test inv!(G, geinv, e) === geinv
                    @test isapprox(G, geinv, e; atol = atol)
                end
            end
        end
    end

    @testset "translation" begin
        convs = ((), (LeftAction(),), (RightAction(),))

        @test isapprox(G, translate(G, g_pts[1], g_pts[2]), compose(G, g_pts[1], g_pts[2]); atol = atol)
        @test isapprox(G, translate(G, g_pts[1], g_pts[2], LeftAction()), compose(G, g_pts[1], g_pts[2]); atol = atol)
        @test isapprox(G, translate(G, g_pts[1], g_pts[2], RightAction()), compose(G, g_pts[2], g_pts[1]); atol = atol)

        for conv in convs
            @test isapprox(G, inverse_translate(G, g_pts[1], translate(G, g_pts[1], g_pts[2], conv...), conv...), g_pts[2]; atol = atol)
            @test isapprox(G, translate(G, g_pts[1], inverse_translate(G, g_pts[1], g_pts[2], conv...), conv...), g_pts[2]; atol = atol)
        end

        test_mutating && @testset "mutating" begin
            for conv in convs
                g = allocate(g_pts[1])
                @test translate!(G, g, g_pts[1], g_pts[2], conv...) === g
                @test isapprox(G, g, translate(G, g_pts[1], g_pts[2], conv...); atol = atol)

                g = translate(G, g_pts[1], g_pts[2], conv...)
                g2 = allocate(g)
                @test inverse_translate!(G, g2, g_pts[1], g, conv...) === g2
                @test isapprox(G, g2, g_pts[2]; atol = atol)

                g = inverse_translate(G, g_pts[1], g_pts[2], conv...)
                g2 = allocate(g)
                @test translate!(G, g2, g_pts[1], g, conv...) === g2
                @test isapprox(G, g2, g_pts[2]; atol = atol)
            end
        end
    end

    test_diff && @testset "translation differential" begin
        v = v_pts[1]
        g21 = compose(G, g_pts[2], g_pts[1])
        g12 = compose(G, g_pts[1], g_pts[2])
        @test isapprox(G, g12, translate_diff(G, g_pts[2], g_pts[1], v), translate_diff(G, g_pts[2], g_pts[1], v, LeftAction()); atol = atol)
        @test is_tangent_vector(G, g12, translate_diff(G, g_pts[2], g_pts[1], v, LeftAction()), true; atol = atol)
        RightAction() in diff_convs && @test is_tangent_vector(G, g21, translate_diff(G, g_pts[2], g_pts[1], v, RightAction()), true; atol = atol)

        for conv in diff_convs
            g2g1 = translate(G, g_pts[2], g_pts[1], conv...)
            g2invg1 = inverse_translate(G, g_pts[2], g_pts[1], conv...)
            @test isapprox(G, g_pts[1], inverse_translate_diff(G, g_pts[2], g2g1, translate_diff(G, g_pts[2], g_pts[1], v, conv...), conv...), v; atol = atol)
            @test isapprox(G, g_pts[1], translate_diff(G, g_pts[2], g2invg1, inverse_translate_diff(G, g_pts[2], g_pts[1], v, conv...), conv...), v; atol = atol)
        end

        test_mutating && @testset "mutating" begin
            for conv in diff_convs
                g2g1 = translate(G, g_pts[2], g_pts[1], conv...)
                g2invg1 = inverse_translate(G, g_pts[2], g_pts[1], conv...)
                vout = allocate(v)
                @test translate_diff!(G, vout, g_pts[2], g_pts[1], v, conv...) === vout
                @test isapprox(G, g2g1, vout, translate_diff(G, g_pts[2], g_pts[1], v, conv...); atol = atol)

                vout = translate_diff(G, g_pts[2], g_pts[1], v, conv...)
                vout2 = allocate(vout)
                @test inverse_translate_diff!(G, vout2, g_pts[2], g2g1, vout, conv...) === vout2
                @test isapprox(G, g_pts[1], vout2, v; atol = atol)

                vout = inverse_translate_diff(G, g_pts[2], g_pts[1], v, conv...)
                vout2 = allocate(vout)
                @test translate_diff!(G, vout2, g_pts[2], g2invg1, vout, conv...) === vout2
                @test isapprox(G, g_pts[1], vout2, v; atol = atol)
            end
        end
    end
end

"""
    test_action(
        A::AbstractGroupAction,
        a_pts::AbstractVector,
        m_pts::AbstractVector;
        atol = 1e-10,
        test_optimal_alignment = false,
        test_mutating = true,
        test_diff = false,
    )

Tests general properties of the action `A`, given at least three different points
that lie on it (contained in `a_pts`) and three different point that lie
on the manifold it acts upon (contained in `m_pts`).
"""
function test_action(
        A::AbstractGroupAction,
        a_pts::AbstractVector,
        m_pts::AbstractVector,
        v_pts = [];
        atol = 1e-10,
        test_optimal_alignment = false,
        test_mutating = true,
        test_diff = false,
    )

    G = base_group(A)
    M = g_manifold(A)
    e = Identity(G)

    @testset "Basic action properties" begin
        @testset "Direction" begin
            Aswitch = switch_direction(A)
            if isa(A, AbstractGroupAction{LeftAction})
                @test direction(A) === LeftAction()
                @test isa(Aswitch, AbstractGroupAction{RightAction})
                @test direction(Aswitch) === RightAction()
            else
                @test direction(A) === RightAction()
                @test isa(Aswitch, AbstractGroupAction{LeftAction})
                @test direction(Aswitch) === LeftAction()
            end
        end

        @testset "Closed" begin
            @testset "over actions" begin
                for a1 in a_pts, a2 in a_pts
                    a3 = compose(A, a1, a2)
                    @test is_manifold_point(G, a3, true; atol = atol)
                end
            end
            @testset "over g-manifold" begin
                for a in a_pts, m in m_pts
                    @test is_manifold_point(M, apply(A, a, m), true; atol = atol)
                    @test is_manifold_point(M, inverse_apply(A, a, m), true; atol = atol)
                end
            end
        end

        @testset "Associative" begin
            a12 = compose(A, a_pts[1], a_pts[2])
            a23 = compose(A, a_pts[2], a_pts[3])

            @testset "over compose" begin
                a12_a3 = compose(A, a12, a_pts[3])
                a1_a23 = compose(A, a_pts[1], a23)
                @test isapprox(G, a12_a3, a1_a23)
            end

            @testset "over apply" begin
                for m in m_pts
                    a12_a3_m = apply(A, a12, apply(A, a_pts[3], m))
                    a1_a23_m = apply(A, a_pts[1], apply(A, a23, m))
                    @test isapprox(M, a12_a3_m, a1_a23_m)
                end
            end

            test_mutating && @testset "mutating" begin
                a12, a23, a12_3, a1_23 = allocate.(repeat([a_pts[1]], 4))
                @test compose!(A, a12, a_pts[1], a_pts[2]) === a12
                @test compose!(A, a23, a_pts[2], a_pts[3]) === a23
                @test compose!(A, a12_3, a12, a_pts[3]) === a12_3
                @test compose!(A, a1_23, a_pts[1], a23) === a1_23
                @test isapprox(G, a12_3, a1_23)

                for m in m_pts
                    a12_a3_m, a1_a23_m = allocate(m), allocate(m)
                    @test apply!(A, a12_a3_m, a12, apply(A, a_pts[3], m)) === a12_a3_m
                    @test apply!(A, a1_a23_m, a_pts[1], apply(A, a23, m)) === a1_a23_m
                    @test isapprox(M, a12_a3_m, a1_a23_m)
                end
            end
        end

        @testset "Identity" begin
            @test compose(A, e, e) === e

            for a in a_pts
                @test isapprox(G, compose(A, a, e), a)
                @test isapprox(G, compose(A, e, a), a)

                ge = identity(G, a)
                @test isapprox(G, compose(A, a, ge), a)
                @test isapprox(G, compose(A, ge, a), a)

                for m in m_pts
                    @test isapprox(M, apply(A, e, m), m)
                    @test isapprox(M, apply(A, ge, m), m)
                    @test isapprox(M, inverse_apply(A, e, m), m)
                    @test isapprox(M, inverse_apply(A, ge, m), m)
                end
            end

            test_mutating && @testset "mutating" begin
                for a in a_pts
                    h = allocate(a)
                    @test compose!(A, h, a, e) === h
                    @test isapprox(G, h, a)
                    h = allocate(a)
                    @test compose!(A, h, e, a) === h
                    @test isapprox(G, h, a)

                    ge = identity(G, a)
                    @test isapprox(G, compose(A, a, ge), a)
                    @test isapprox(G, compose(A, ge, a), a)

                    ge = identity(G, a)
                    @test compose!(A, ge, e, e) === ge
                    @test isapprox(G, ge, e)

                    for m in m_pts
                        em = allocate(m)
                        @test apply!(A, em, e, m) === em
                        @test isapprox(M, em, m)
                        em = allocate(m)
                        @test apply!(A, em, ge, m) === em
                        @test isapprox(M, em, m)
                        em = allocate(m)
                        @test inverse_apply!(A, em, e, m) === em
                        @test isapprox(M, em, m)
                        em = allocate(m)
                        @test inverse_apply!(A, em, ge, m) === em
                        @test isapprox(M, em, m)
                    end
                end
            end
        end

        @testset "Inverse" begin
            for a in a_pts
                ainv = inv(G, a)
                @test isapprox(G, compose(A, a, ainv), e)
                @test isapprox(G, compose(A, ainv, a), e)
                @test isapprox(G, e, compose(A, a, ainv))
                @test isapprox(G, e, compose(A, ainv, a))

                for m in m_pts
                    @test isapprox(M, apply(A, a, m), inverse_apply(A, ainv, m))
                    @test isapprox(M, apply(A, ainv, m), inverse_apply(A, a, m))
                end
            end
        end
    end

    test_diff && @testset "apply differential" begin
        for (m, v) in zip(m_pts, v_pts)
            for a in a_pts
                am, av = apply(A, a, m), apply_diff(A, a, m, v)
                ainvm, ainvv = inverse_apply(A, a, m), inverse_apply_diff(A, a, m, v)
                @test is_tangent_vector(M, am, av, true; atol = atol)
                @test is_tangent_vector(M, ainvm, ainvv, true; atol = atol)
            end

            a12 = compose(A, a_pts[1], a_pts[2])
            a2m = apply(A, a_pts[2], m)
            a12v = apply_diff(A, a12, m, v)
            a2v = apply_diff(A, a_pts[2], m, v)
            @test isapprox(M, a2m, apply_diff(A, a_pts[1], a2m, a2v), a12v; atol = atol)

            @test isapprox(M, m, apply_diff(A, e, m, v), v; atol = atol)
            @test isapprox(M, m, inverse_apply_diff(A, e, m, v), v; atol = atol)
        end

        test_mutating && @testset "mutating" begin
            for (m, v) in zip(m_pts, v_pts)
                for a in a_pts
                    am = apply(A, a, m)
                    av = allocate(v)
                    @test apply_diff!(A, av, a, m, v) === av
                    ainvm = inverse_apply(A, a, m)
                    ainvv = allocate(v)
                    @test inverse_apply_diff!(A, ainvv, a, m, v) === ainvv
                    @test is_tangent_vector(M, am, av, true; atol = atol)
                    @test is_tangent_vector(M, ainvm, ainvv, true; atol = atol)
                end

                a12 = compose(A, a_pts[1], a_pts[2])
                a2m = apply(A, a_pts[2], m)
                a12m = apply(A, a12, m)
                a12v, a2v, a1_a2v = allocate(v), allocate(v), allocate(v)
                @test apply_diff!(A, a12v, a12, m, v) === a12v
                @test apply_diff!(A, a2v, a_pts[2], m, v) === a2v
                @test apply_diff!(A, a1_a2v, a_pts[1], a2m, a2v) === a1_a2v
                @test isapprox(M, a12m, a1_a2v, a12v; atol = atol)

                ev = allocate(v)
                @test apply_diff!(A, ev, e, m, v) === ev
                @test isapprox(G, m, ev, v; atol = atol)
                ev = allocate(v)
                @test inverse_apply_diff!(A, ev, e, m, v) === ev
                @test isapprox(G, m, ev, v; atol = atol)
            end
        end
    end

    test_optimal_alignment && @testset "Center of orbit" begin
        act = center_of_orbit(A, [m_pts[1]], m_pts[2])
        act2 = center_of_orbit(A, [m_pts[1]], m_pts[2], GradientDescentEstimation())
        act_opt = optimal_alignment(A, m_pts[2], m_pts[1])
        @test isapprox(G, act, act_opt; atol = atol)
        @test isapprox(G, act2, act_opt; atol = atol)

        test_mutating && @testset "mutating" begin
            act_opt2 = allocate(act_opt)
            optimal_alignment!(A, act_opt2, m_pts[2], m_pts[1])
            @test isapprox(G, act_opt, act_opt2; atol = atol)
        end
    end
end

ext: jl


text:
function born_approximation(output::Ts, input::Ts, iflag::Ti;
         irz::Vector, irx::Vector, location_flag="index", dir_sourceside="NULL", normalization_flag=true, mute_index=0, fidiff::TdParams) where {Ts<:String, Ti<:Int64}

    # clean the space
    rm(output, force=true, recursive=true)

    if iflag == 1
       born_approximation_forward(output, input, irz, irx, dir_sourceside, fidiff;
                                  location_flag=location_flag, normalization_flag=normalization_flag, mute_index=mute_index)

    elseif iflag == 2
       born_approximation_adjoint(output, input, dir_sourceside, fidiff;
                                  normalization_flag=normalization_flag, mute_index=mute_index)

    else
       error("undefined operator behavior")
    end

    return nothing
end

"""
   perform y = ax + by for multiple recordings
"""
function recordings_axpby!(a, dir_x::Ts, b, dir_y::Ts) where {Ts<:String}

    # determine number of files
    file_name = readdir(dir_x)
    nx        = length(file_name)
    if file_name[1] == ".DS_Store" # in case of MacOS
       nx = nx - 1
    end

    file_name = readdir(dir_y)
    ny        = length(file_name)
    if file_name[1] == ".DS_Store" # in case of MacOS
       ny = ny - 1
    end

    nx == ny || error("number of files doesn't match")

    r = 0.0
    for i = 1 : nx
        file_name = join(["recordings_" "$i" ".bin"])

        path_x    = joinpath(dir_x, file_name)
        path_y    = joinpath(dir_y, file_name)

        rec_x     = read_recordings(path_x)
        rec_y     = read_recordings(path_y)

        for i2 = 1 : rec_x.nr
            for i1 = 1 : rec_x.nt
                rec_y.p[i1,i2] = a*rec_x.p[i1,i2] + b*rec_y.p[i1,i2]
            end
        end

        write_recordings(path_y, rec_y)
    end

    return nothing
end

"""
   perform y = ax + by for variables in model space
"""
function image_axpby!(a, path_x::Ts, b, path_y::Ts) where {Ts<:String}

    (hdr1, x) = read_RSdata(path_x)
    (hdr2, y) = read_RSdata(path_y)

    for i = 1 : length(x)
        y[i] = a*x[i] + b*y[i]
    end

    write_RSdata(path_y, hdr2, y)

    return nothing
end

function recordings_norm(dir::Ts) where {Tv<:Number, Ts<:String}

    # determine number of files
    file_name = readdir(dir)
    n         = length(file_name)
    if file_name[1] == ".DS_Store" # in case of MacOS
       n = n - 1
    end

    r = 0.0
    for i = 1 : n
        file_name = join(["recordings_" "$i" ".bin"])
        path      = joinpath(dir, file_name)
        rec       = read_recordings(path)
        tmp       = norm(rec.p)
        r        += tmp * tmp
    end

    return sqrt(r)
end


text:
# # Tutorial: missing value handling
#
# This tutorial illustrates the usage of Transducers.jl by stepping
# through various handling of missing values.

using Transducers
using LiterateTest                                                     #src
using Test                                                             #src

# ## Dot product
#
# Here is a simple way to compute a dot product using
# [`MapSplat`](@ref):

@test begin
    zip(1:3, 10:2:14) |> MapSplat(*) |> sum
end == 76

# Let's see what it does step by step.  First we create a "printer"
# transducer using the following function (see [`Map`](@ref)):

xf_printer(label) = Map() do x
    println(label, ": ", x)
    return x  # just return it as-is
end
nothing  # hide

# This transducer just pass-through the input while printing its value
# (prefixed by a `label`).  Let's sandwich the previous `MapSplat(*)`
# with it:

@test begin
    zip(1:3, 10:2:14) |>
    xf_printer(" input") |>
    MapSplat(*) |>
    xf_printer("output") |>
    sum
end == 76
flush(stdout)  # hide

# You can see that the input tuple `(1, 10)` is splatted into function
# `*` by `MapSplat` which then outputs `10`.  This is repeated for all
# inputs.

# Perhaps unfortunately, this way of computing a dot product
# propagates any missing values contained in the input arrays to the
# result (which may actually be desired in certain cases).

xs = [1, missing, 3, 2]
ys = [10, 14, missing, 12]
zip(xs, ys) |> MapSplat(*) |> sum

# However, it is very simple to ignore any missing values using
# [`OfType`](@ref):

@test begin
    zip(xs, ys) |> OfType(Tuple{Vararg{Number}}) |> MapSplat(*) |> sum
end == 34

# Here, `Tuple{Vararg{Number}}` is a type that matches with a tuple of
# any length with numbers.  It does not match with a tuple if it has a
# `missing`.

@assert (1, 0.5) isa Tuple{Vararg{Number}}
@assert (1, 0.5, 2im) isa Tuple{Vararg{Number}}
@assert !((1, missing) isa Tuple{Vararg{Number}})

# The part `... |> OfType(Tuple{Vararg{Number}}) |> MapSplat(*)` can
# be factored out using [`opcompose`](@ref ∘):

xf_mdot = opcompose(OfType(Tuple{Vararg{Number}}), MapSplat(*));

# or equivalently (in Julia ≥ 1.5):
#
# ```julia
# OfType(Tuple{Vararg{Number}}) ⨟ MapSplat(*)
# ```

# The transducer `xf_mdot` can be used where previously
# `OfType(Tuple{Vararg{Number}}) |> MapSplat(*)` was used:

@test begin
    zip(xs, ys) |> xf_mdot |> sum
end == 34

# ## Covariance
#
# Transducer `xf_mdot` above can also be used to compute the
# covariance.  First, we need the number of pairs of elements in `xs`
# and `ys` that _both_ of them are not `missing`:

nonmissings = zip(xs, ys) |> Map(x -> x isa Tuple{Vararg{Number}}) |> count
@test nonmissings == 2

# Finally, we have to pre-process the input to `xf_mdot` by
# subtracting the average.  It's easy to do this with `Map`:

using Statistics: mean

function xf_demean(xs, ys)
    xmean = mean(skipmissing(xs))
    ymean = mean(skipmissing(ys))
    return Map(((x, y),) -> (x - xmean, y - ymean))
end
nothing  # hide

# We can then compute the covariance by combining `xf_demean` and
# `xf_mdot`:

@test begin
    s = zip(xs, ys) |> xf_demean(xs, ys) |> xf_mdot |> sum
    s / nonmissings
end == 1.0

# In `xf_demean`, the averages of the vectors `xs` and `ys` are
# computed separately.  It is also easy to compute the averages of the
# elements where both `xs` and `ys` are non-`missing`:

function xf_demean2(xs, ys)
    n, xsum, ysum =
        zip(xs, ys) |>
        OfType(Tuple{Vararg{Number}}) |>
        Map(((x, y),) -> (1, x, y)) |>
        Broadcasting() |>
        sum
    xmean = xsum / n
    ymean = ysum / n
    return Map(((x, y),) -> (x - xmean, y - ymean))
end

if VERSION >= v"1.1"                                                   #src
    @test begin
        s = zip(xs, ys) |> xf_demean2(xs, ys) |> xf_mdot |> sum
        s / nonmissings
    end == 0.5
end                                                                    #src

# In `xf_demean2`, we used [`Broadcasting`](@ref) transducer to
# broadcast elements of the tuple `(1, x, y)` over the reducing
# function of `sum` (i.e., `+`).

# ### Advanced: `TeeRF` and `ProductRF`
#
# Alternatively, it can also be computed using by combining
# [`foldxl`](@ref), [`ProductRF`](@ref), [`TeeRF`](@ref), and
# [`DataTools.inc1`](https://juliafolds.github.io/DataTools.jl/dev/)
# (see [below](@ref tutorial-findminmax) for how `TeeRF` and
# `ProductRF` work):

using DataTools: inc1

function xf_demean3(xs, ys)
    n, (xsum, ysum) =
        zip(xs, ys) |>
        OfType(Tuple{Vararg{Number}}) |>
        foldxl(TeeRF(inc1, ProductRF(+, +)))
    xmean = xsum / n
    ymean = ysum / n
    return Map(((x, y),) -> (x - xmean, y - ymean))
end

@test begin
    s = zip(xs, ys) |> xf_demean3(xs, ys) |> xf_mdot |> sum
    s / nonmissings
end == 0.5

# ## Vectorized reduction
#
# `reduce`, `sum`, etc. in `Base` support `dims` argument.
# Transducers.jl does not support this argument as of writing.
# However, this can easily be emulated using `eachcol`, `eachrow`, or
# `eachslice` iterators in `Base`.

xs = [
    0       missing 1       2
    3       4       5       missing
    missing 6       7       missing
]

if VERSION < v"1.1"
    using Compat: eachcol
end
@test begin
    eachcol(xs) |> Broadcasting() |> NotA(Missing) |> sum
end == [3, 12, 13]

# Here, we use [`NotA`](@ref) transducer that filters out `missing`
# values:

@test begin
    [1, 2, missing, 3] |> NotA(Missing) |> collect
end == [1, 2, 3]

# Above computation returns the sum over each row without taking into
# account the relationship within a column.  Another possibly useful
# reduction is the sum of the columns with no missing values.  This
# can easily be done by filtering before:

@test begin
    eachcol(xs) |> Filter(x -> !any(ismissing, x)) |> Broadcasting() |> sum
end == [1, 5, 7]

# ## `findmax` and `findmin`
#
# Another useful operation is `findmax`/`findmin`.  Using `Filter`,
# `missing` values can be filtered out by

filtered_pairs = [1, 3, missing, 0] |> pairs |> Filter(!(ismissing ∘ last))
@test begin
    collect(filtered_pairs)
end == [1 => 1, 2 => 3, 4 => 0]

# These key-value pairs can be accumulated by the following reducing
# step function:

##                     ,--- current state
##                     |
##                     |              ,-- input
##                     |              |
function findmax_step((argmax, max), (index, value))
    argmax, max = value > max ? (index, value) : (argmax, max)
    return argmax => max
    ##        \
    ##         \__ next state
end

@test begin
    foldxl(findmax_step, filtered_pairs)
end == (2 => 3)

# or equivalently

@test begin
    [1, 3, missing, 0] |> pairs |> Filter(!(ismissing ∘ last)) |> foldxl(findmax_step)
end == (2 => 3)

# [`foldxl`](@ref) is like [`foldl`](@refe) but always uses
# Transducers.jl's extended fold protocol.  It also has the unary
# curried method `foldxl(rf)` defined as `xs -> foldxl(rf, xs)`.  It
# is handy to use in the piping context as in the latter example.

# [`DataTools.rightif`](https://juliafolds.github.io/DataTools.jl/dev/)
# can be used for defining `findmax`/`findmin`-like functions on the
# fly:

using DataTools: rightif

@test begin
    foldxl(rightif(<, last), filtered_pairs)
end === Pair{Int64,Union{Missing, Int}}(2, 3)

# ### Side note: why `Pair{Int64,Union{Missing,Int}}`?
#
# The result type just above using `rightif` is
# `Pair{Int64,Union{Missing,Int}}`:

@test begin
    typeof(foldxl(rightif(<, last), filtered_pairs))
end === Pair{Int64,Union{Missing,Int}}

# This is because that's the element type of `pairs([1, 3, missing,
# 0])` and `rightif` does not re-construct the input `Pair` like
# `findmax_step`:

@test begin
    [1, 3, missing, 0] |> pairs |> first |> typeof
end === Pair{Int64,Union{Missing, Int}}

# We can avoid this by pre-processing the input with `MapSplat(Pair)`:

@test begin
    foldxl(rightif(<, last), filtered_pairs |> MapSplat(Pair))
end === (2 => 3)

# ### `findmin`

# Similarly, we can define `findmin` with

function findmin_step((argmin, min), (index, value))
    argmin, min = value < min ? (index, value) : (argmin, min)
    return argmin => min
end

@test begin
    foldxl(findmin_step, filtered_pairs)
end == (4 => 0)

# and

@test begin
    foldxl(rightif(>, last), filtered_pairs)
end == (4 => 0)

# ## [Extrema (`findminmax`)](@id tutorial-findminmax)
#
# To compute `findmax` and `findmax` in a single sweep, we can use
# [`TeeRF`](@ref) to "fan out" the input stream to multiple reducing
# step functions:

@test begin
    foldxl(TeeRF(findmin_step, findmax_step), filtered_pairs)
end == (4 => 0, 2 => 3)

# or equivalently

@test begin
    foldxl(TeeRF(rightif(>, last), rightif(<, last)), filtered_pairs)
end == (4 => 0, 2 => 3)

# In general, multiple folds on a same collection
#
# ```julia
# a₁ = foldxl(rf₁, xs)
# a₂ = foldxl(rf₂, xs)
# ...
# aₙ = foldxl(rfₙ, xs)
# ```
#
# can be fused into a single fold using `TeeRF`
#
# ```julia
# a₁, a₂, ..., aₙ = foldxl(TeeRF(rf₁, rf₂, ..., rfₙ), xs)
# ```
#
# provided that the input collection `xs` and the reducing functions
# `rf₁`, `rf₂`, ..., and `rfₙ` are not stateful.

# ### More fusing by transforming reducing functions

# In the above computation, we have a [reducing (step) function](@ref
# glossary-rf)

rf = TeeRF(rightif(>, last), rightif(<, last));

# a [transducer](@ref glossary-transducer)

xf = Filter(!(ismissing ∘ last));

# and an iterable

xs = pairs([1, 3, missing, 0]);

# In Transducers.jl, a transducer acts as iterator transformation
# `xf(xs)` as well as reducing function transformation `xf'(rf)`.
# Thus, the following calls are equivalent:

#src `xf(xs)` impossible before https://github.com/JuliaLang/julia/pull/31916
@dedent if VERSION >= v"1.3"
    @assert foldxl(rf, xf, xs) == (4 => 0, 2 => 3)
    @assert foldxl(rf, xf(xs)) == (4 => 0, 2 => 3)
    @assert foldxl(xf'(rf), xs) == (4 => 0, 2 => 3)
end

# By exploiting this equality, we can fuse more computations by moving
# the transformation on the side of reducing function.  For example,
# we can compute non-missing extrema and count missings at the same
# time:

@test begin
    [1, 3, missing, 0] |>
    pairs |>
    MapSplat(Pair) |>  # avoid `Pair{Int64,Union{Missing, Int}}`
    foldxl(TeeRF(
        Map(ismissing ∘ last)'(+),  # count number of missings
        Filter(!(ismissing ∘ last))'(TeeRF(
            rightif(>, last),  # find non-missing minimum
            rightif(<, last),  # find non-missing maximum
        )),
    ))
end == (1, (4 => 0, 2 => 3))

# Using `ProductRF`, we can compute `findmax` of individual and
# `zip`ped items at the same time

@test begin
    zip(
        pairs([1, 3, missing, 0]),  # produces k1 => v1
        pairs([4, missing, 6, 5]),  # produces k2 => v2
    ) |>
    Map() do ((k1, v1), (k2, v2))
        (k1 => v1, k2 => v2, (k1, k2) => (v1, v2))
    end |>
    foldxl(ProductRF(
        Filter(!(ismissing ∘ last))'(rightif(<, last)),  # max of k1 => v1
        Filter(!(ismissing ∘ last))'(rightif(<, last)),  # max of k2 => v2
        Filter(((_, (v1, v2)),) -> v1 !== missing && v2 !== missing)'(
            rightif(<, last)  # max (k1, k2) => (v1, v2)
        ),
    ))
end == (2 => 3, 3 => 6, (1, 1) => (1, 4))

# `ProductRF` is like `TeeRF` but acts on the input that is already a
# tuple.  That is to say, given a collection of `n`-tuple `xs`,
# multiple folds on a same collection
#
# ```julia
# a₁ = foldxl(rf₁, (x[1] for x in xs))
# a₂ = foldxl(rf₂, (x[2] for x in xs))
# ...
# aₙ = foldxl(rfₙ, (x[n] for x in xs))
# ```
#
# can be fused into a single fold using `ProductRF`
#
# ```julia
# a₁, a₂, ..., aₙ = foldxl(ProductRF(rf₁, rf₂, ..., rfₙ), xs)
# ```
#
# provided that the input collection `xs` and the reducing functions
# `rf₁`, `rf₂`, ..., and `rfₙ` are not stateful.

# ## Early termination

# `Base`'s `maximum` reports the maximum to be `missing` when it
# receives a container with a `missing`:

@test begin
    maximum([1, 2, missing, 3])
end === missing

# We can obtain the same behavior by using `isless` instead of `>` in
# `findmax_step′`:

function findmax_step′((argmax, max), (index, value))
    argmax, max = isless(max, value) ? (index, value) : (argmax, max)
    return argmax => max
end

@test begin
    foldl(findmax_step′, pairs([1, 2, missing, 3]))
end === (3 => missing)

@testset "foldl(findmax_step′, pairs(...))" begin
    f(xs) = foldl(findmax_step′, pairs(xs))
    @test f([1, 2, 3]) === (3 => 3)
    @test f([1, NaN, 3]) === (2 => NaN)
    @test f([1, missing, 3, missing]) === (2 => missing)
end

# `foldl(findmax_step′, ...)` does not stop even after it observed
# `missing`.  We can easily add early termination by using
# [`ReduceIf`](@ref):

@test begin
    [1, 2, missing, 3] |> pairs |> ReduceIf(ismissing ∘ last) |> foldxl(findmax_step′)
end === (3 => missing)

# Note that [`ReduceIf(f)`](@ref ReduceIf) is not same as
# [`TakeWhile(!f)`](@ref TakeWhile):

@test begin
    [1, 2, missing, 3] |> pairs |> TakeWhile(!(ismissing ∘ last)) |> foldxl(findmax_step′)
end === (2 => 2)

# That is to say, `TakeWhile` does not evaluate the inner reducing
# step function with the item that triggers the early termination.
# That's why we need `ReduceIf` here.

# ### `findmin` with `missing` and `NaN`

# Unfortunately, we do not have `min`-compatible "total" ordering in
# `Base`.  Thus, we need to create a function that special-cases
# `missing` and `NaN`:

function isgreater(x, y)
    xisnan = x != x
    xisnan isa Bool || return false  # x is missing
    yisnan = y != y
    yisnan isa Bool || return true  # y is missing
    xisnan && return false
    yisnan && return true
    return isless(y, x)
end

@assert isgreater(2, 1)
@assert isgreater(1, missing)
@assert isgreater(NaN, missing)
@assert isgreater(1, NaN)
@assert !isgreater(1, 2)
@assert !isgreater(missing, 1)
@assert !isgreater(missing, NaN)
@assert !isgreater(NaN, 1)

# Using `isgreater` instead of `<`, we can define `findmin_step′` like
# `findmax_step′`:

function findmin_step′((argmin, min), (index, value))
    argmin, min = isgreater(min, value) ? (index, value) : (argmin, min)
    return argmin => min
end

@test begin
    foldl(findmax_step′, pairs([1, 2, missing, 3]))
end === (3 => missing)

@testset "foldl(findmin_step′, pairs(...))" begin
    f(xs) = foldl(findmin_step′, pairs(xs))
    @test f([1, 2, 3]) === (1 => 1)
    @test f([1, NaN, 3]) === (2 => NaN)
    @test f([1, missing, 3, missing]) === (2 => missing)
end

# ### Extrema (`findminmax`) with early termination

# As before, we can fuse `findmin_step′` and `findmax_step′` using
# `TeeRF`.  This can also be composed with `ReduceIf`:

@test begin
    [1, 2, 3, 0, 2] |>
    pairs |>
    ReduceIf(ismissing ∘ last) |>
    foldxl(TeeRF(findmin_step′, findmax_step′))
end === ((4 => 0), (3 => 3))
#-

@test begin
    [1, 2, missing, 0, 2] |>
    pairs |>
    ReduceIf(ismissing ∘ last) |>
    foldxl(TeeRF(findmin_step′, findmax_step′))
end === ((3 => missing), (3 => missing))

# ## Ad-hoc imputation

# Using [`Scan`](@ref), it is straightforward to fill `missing` items
# with the last non-`missing` item (last observation carried forward):

@test begin
    [1, 3, missing, 0, 2, missing, missing] |>
    pairs |>
    MapSplat(Pair) |>
    Scan() do last, (k, v)
        ismissing(v) ? last : k => v
    end |>
    collect
end == [1 => 1, 2 => 3, 2 => 3, 4 => 0, 5 => 2, 5 => 2, 5 => 2]

# `rightif` can also be used:

@test begin
    [1, 3, missing, 0, 2, missing, missing] |>
    pairs |>
    MapSplat(Pair) |>
    Scan(rightif(!(ismissing ∘ right), last)) |>
    collect
end == [1 => 1, 2 => 3, 2 => 3, 4 => 0, 5 => 2, 5 => 2, 5 => 2]
nothing  # hide

# Note that the output still may contain `missing` if the first item
# is `missing`:

@test begin
    [missing, 1, missing] |>
    pairs |>
    MapSplat(Pair) |>
    Scan(rightif(!(ismissing ∘ right), last)) |>
    collect
end |> isequal([1 => missing, 2 => 1, 2 => 1])

# This can be worked around by specifying `init` argument for `Scan`:

@test begin
    [missing, 1, missing] |>
    pairs |>
    MapSplat(Pair) |>
    Scan(rightif(!(ismissing ∘ right), last), 0 => 0) |>
    collect
end == [0 => 0, 2 => 1, 2 => 1]


text:
raw"""
# Conjugate Gradient Method

## References

- [wiki]: https://en.wikipedia.org/wiki/Conjugate_gradient_method
- [painless]: Shewchuk, Jonathan Richard. "An introduction to the conjugate gradient method without
  the agonizing pain." (1994): 1.

The derivation in <https://en.wikipedia.org/wiki/Derivation_of_the_conjugate_gradient_method> is
rather complicated and unintuitive, so it's better to go through [painless].
"""
module ConjugateGradientMethod

using LinearAlgebra: norm, I, (⋅)

export conjugate_gradient, conjugate_directions

#=

At first "conjugate" sounds confusing, but it becomes clear after realizing it is actually "$p_i$
and $p_j$ (i ≠ j) are conjugate **w.r.t. $A$**" (while $r_i^T r_j = 0$ for i ≠ j). So the whole
iteration is indeed doing Gram-Schmidt process to find $A$-orthogonal search directions.

Recall some definitions:

    x_* = \sum _{i=1}^{n} α_i p_i  ...mutually conjugate vectors, w.r.t. $A$
    α_i = <p_i, r_i> / <p_i, p_i>_A

    r_i = b - A x_i  ...the residual, negative gradient of $f(x) = 1/2 x^T A x - b^T x$ at $x_i$
    e_i = x_i - x_*  ...the error

Some notes taken from [painless]:

- "In fact, if the search vectors are constructed by conjugation of the axial unit vectors,
  Conjugate Directions becomes equivalent to performing Gaussian elimination."
- In Conjugate Directions, $e_{i + 1}$ is $A$-orthogonal to $p_{i}$. This is equivalent to finding
  the minimum point along $p_i$.
- Conjugate Directions chooses value from $e_i = e_0 + P_i$ that minimizes $<e_0, e_0>_A$, where
  $P_i$ is `span {p_0, p_1, ..., p_{i - 1}}`. (See [painless] 7.3)
- Moreover, not only $<p_{i - 1}, e_i>_A = 0$, but for $i < j$, $<p_i, e_j>_A = 0$ also holds. As
  $p_i^T A e_j = 0$ and $r_j = -A e_j$, we know $p_i$ and $r_j$ are orthogonal.
- Finally, CG is simply a special case of Conjugate Directions where search directions are
  constructed by conjugation of the residuals. This choice makes $r_i^T r_j = 0$ and $p_i^T r_i =
  r_i^T r_i$ and $r_i^T A d_i = d_i^T A d_i$ and $r_{i + 1}$ being $A$-orthogonal to $P_i$, etc.,
  many nice properties, so there is no need to store all the previous search directions for doing
  Gram-Schmidt process.

Recapitulate some derivations of equations:

    α_i = <p_i, r_i> / <p_i, p_i>_A    (CD, given $p_i$, minimization)
        = <r_i, r_i> / <p_i, p_i>_A    (CG, u_i = r_i, c.f. [painless] Eq. 42)

    β_ij = -<u_i, p_j>_A / <p_j, p_j>_A    (CD, Gram-Schmidt, c.f. [painless] Eq. 37)

    β_i := β_{i}{i - 1}
        = -<r_{i + 1}, p_i>_A / <p_i, p_i>_A    (CG, Gram-Schmidt in Krylov subspace)
        = <r_{i + 1}, r_{i + 1}> / <r_i, r_i>    (CD, c.f. [painless] Eq. 43, definition of α_i)

=#

# Just a proof-of-concept. VERY slow.
#
# <https://en.wikipedia.org/wiki/Derivation_of_the_conjugate_gradient_method#The_conjugate_direction_method>
function conjugate_directions(
    A::AbstractMatrix{T},
    b::Vector{T},
)::Vector{T} where T<:AbstractFloat
    n = length(b)
    x = zeros(T, n)
    if norm(b) < eps(T)
        return x
    end

    U = Matrix(T.(I(n)))
    P = zeros(T, n, n)
    u(k::Int) = @view U[:, k]
    p(k::Int) = @view P[:, k]

    r = b
    for i in 1:n
        # Gram-Schmidt orthogonalization, directly translated from its definition.
        t = u(i)
        for j = 1:(i - 1)
            Ap = A * p(j)
            β = -(u(i)' * Ap) / (p(j)' * Ap)
            t += β * p(j)
        end
        p(i) .= t
        Ap = A * p(i)
        α = (p(i)' * r) / (p(i)' * Ap)
        x += α * p(i)
        r -= α * Ap
    end
    x
end

function conjugate_gradient(
    A::AbstractMatrix{T},
    b::Vector{T},
)::Vector{T} where T<:AbstractFloat
    n = length(b)
    x = zeros(T, n)
    if norm(b) < eps(T)
        return x
    end

    r = b - A * x
    p = r
    rr = r ⋅ r

    MAX_ITERS = n * 10
    for i in 1:MAX_ITERS
        Ap = A * p
        α = rr / (p ⋅ Ap)
        x += α * p
        r -= α * Ap
        rr_new = r ⋅ r
        if sqrt(rr_new) < eps(T)
            break
        end
        β = rr_new / rr
        p = r + β * p
        rr = rr_new
    end
    x
end

# TODO: Preconditioning.

end  # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
module Util
using ..DSP: xcorr
import Base: *
import LinearAlgebra.BLAS
using LinearAlgebra: mul!
using FFTW

export  hilbert,
        fftintype,
        fftouttype,
        fftabs2type,
        nextfastfft,
        dB,
        dBa,
        pow2db,
        amp2db,
        db2pow,
        db2amp,
        rms,
        rmsfft,
        meanfreq,
        unsafe_dot,
        shiftin!,
        finddelay,
        shiftsignal,
        shiftsignal!,
        alignsignals,
        alignsignals!

function hilbert(x::StridedVector{T}) where T<:FFTW.fftwReal
# Return the Hilbert transform of x (a real signal).
# Code inspired by Scipy's implementation, which is under BSD license.
    N = length(x)
    X = zeros(Complex{T}, N)
    p = plan_rfft(x)
    mul!(view(X, 1:(N >> 1)+1), p, x)
    for i = 2:div(N, 2)+isodd(N)
        @inbounds X[i] *= 2.0
    end
    return ifft!(X)
end
hilbert(x::AbstractVector{T}) where {T<:Real} = hilbert(convert(Vector{fftintype(T)}, x))

"""
    hilbert(x)

Computes the analytic representation of x, ``x_a = x + j
\\hat{x}``, where ``\\hat{x}`` is the Hilbert transform of x,
along the first dimension of x.
"""
function hilbert(x::AbstractArray{T}) where T<:Real
    N = size(x, 1)
    xc = Vector{fftintype(T)}(undef, N)
    X = Vector{fftouttype(T)}(undef, N)
    out = similar(x, fftouttype(T))

    p1 = plan_rfft(xc)
    Xsub = view(X, 1:(N >> 1)+1)
    p2 = plan_bfft!(X)

    normalization = 1/N
    off = 1
    for i = 1:Base.trailingsize(x, 2)
        copyto!(xc, 1, x, off, N)

        # fft
        fill!(X, 0)
        mul!(Xsub, p1, xc)

        # scale real part
        for i = 2:div(N, 2)+isodd(N)
            @inbounds X[i] *= 2.0
        end

        # ifft
        mul!(X, p2, X)

        # scale and copy to output
        @simd for j = 1:N
            @inbounds out[off+j-1] = X[j]*normalization
        end

        off += N
    end

    out
end

## FFT TYPES

# Get the input element type of FFT for a given type
fftintype(::Type{T}) where {T<:FFTW.fftwNumber} = T
fftintype(::Type{T}) where {T<:Real} = Float64
fftintype(::Type{T}) where {T<:Complex} = ComplexF64

# Get the return element type of FFT for a given type
fftouttype(::Type{T}) where {T<:FFTW.fftwComplex} = T
fftouttype(::Type{T}) where {T<:FFTW.fftwReal} = Complex{T}
fftouttype(::Type{T}) where {T<:Union{Real,Complex}} = ComplexF64

# Get the real part of the return element type of FFT for a given type
fftabs2type(::Type{Complex{T}}) where {T<:FFTW.fftwReal} = T
fftabs2type(::Type{T}) where {T<:FFTW.fftwReal} = T
fftabs2type(::Type{T}) where {T<:Union{Real,Complex}} = Float64

# Get next fast FFT size for a given signal length
const FAST_FFT_SIZES = [2, 3, 5, 7]
"""
    nextfastfft(n)

Return the closest product of 2, 3, 5, and 7 greater than or equal
to `n`. FFTW contains optimized kernels for these sizes and
computes Fourier transforms of input that is a product of these
sizes faster than for input of other sizes.
"""
nextfastfft(n) = nextprod(FAST_FFT_SIZES, n)
nextfastfft(ns::Tuple) = nextfastfft.(ns)
nextfastfft(ns...) = nextfastfft(ns)

## COMMON DSP TOOLS

struct dBconvert end
struct dBaconvert end
const dB = dBconvert()
const dBa = dBaconvert()
# for using e.g. 3dB or -3dBa
*(a::Real, ::dBconvert) = db2pow(a)
*(a::Real, ::dBaconvert) = db2amp(a)

"""
    db2pow(a)

Convert dB to a power ratio. This function call also be called
using `a*dB`, i.e. `3dB == db2pow(3)`. The inverse of `pow2db`.
"""
db2pow(a::Real) = 10^(a/10)

"""
    db2amp(a)

Convert dB to an amplitude ratio. This function call also be called
using `a*dBa`, i.e. `3dBa == db2amp(3)`. The inverse of `amp2db`.
"""
db2amp(a::Real) = 10^(a/20)

"""
    pow2db(a)

Convert a power ratio to dB (decibel), or ``10\\log_{10}(a)``.
The inverse of `db2pow`.
"""
pow2db(a::Real) = 10*log10(a)

"""
    amp2db(a)

Convert an amplitude ratio to dB (decibel), or ``20
\\log_{10}(a)=10\\log_{10}(a^2)``. The inverse of `db2amp`.
"""
amp2db(a::Real) = 20*log10(a)

"""
    rms(s)

Return the root mean square of signal `s`.
"""
rms(s::AbstractArray{T}) where {T<:Number} = sqrt(sum(abs2, s)/length(s))

"""
    rmsfft(f)

Return the root mean square of signal `s` given the FFT transform
`f = fft(s)`. Equivalent to `rms(ifft(f))`.
"""
rmsfft(f::AbstractArray{T}) where {T<:Complex} = sqrt(sum(abs2, f))/length(f)

"""
    meanfreq(x, fs)

Calculate the mean power frequency of `x` with a sampling frequency of `fs`, defined as:
```math
MPF = \\frac{\\sum_{i=1}^{F} f_i X_i^2 }{\\sum_{i=0}^{F} X_i^2 } Hz
```
where ``F`` is the Nyquist frequency, and ``X`` is the power spectral density.
"""
function meanfreq(x::AbstractVector{<:Real}, fs=2*π)
    pxx = abs2.(rfft(x))

    len = length(x)
    npoints = fld(len,2)
    freqrg = fs/len.*(0:(npoints))

    mf = sum(pxx.*freqrg)./sum(pxx)
    return mf
end

# Computes the dot product of a single column of a, specified by aColumnIdx, with the vector b.
# The number of elements used in the dot product determined by the size(A)[1].
# Note: bIdx is the last element of b used in the dot product.
function unsafe_dot(a::AbstractMatrix, aColIdx::Integer, b::AbstractVector, bLastIdx::Integer)
    aLen     = size(a, 1)
    bBaseIdx = bLastIdx - aLen
    dotprod  = a[1, aColIdx] * b[ bBaseIdx + 1]
    @simd for i in 2:aLen
        @inbounds dotprod += a[i, aColIdx] * b[bBaseIdx + i]
    end

    return dotprod
end

@inline function unsafe_dot(a::Matrix{T}, aColIdx::Integer, b::Vector{T}, bLastIdx::Integer) where T<:BLAS.BlasReal
    BLAS.dot(size(a, 1), pointer(a, size(a, 1)*(aColIdx-1) + 1), 1, pointer(b, bLastIdx - size(a, 1) + 1), 1)
end

function unsafe_dot(a::AbstractMatrix, aColIdx::Integer, b::AbstractVector{T}, c::AbstractVector{T}, cLastIdx::Integer) where T
    aLen = size(a, 1)
    bLen = length(b)
    bLen == aLen-1  || error( "length(b) must equal to length(a)[1] - 1" )
    cLastIdx < aLen || error( "cLastIdx but be < length(a)")

    dotprod = a[1, aColIdx] * b[cLastIdx]
    @simd for i in 2:aLen-cLastIdx
        @inbounds dotprod += a[i, aColIdx] * b[i+cLastIdx-1]
    end
    @simd for i in 1:cLastIdx
        @inbounds dotprod += a[aLen-cLastIdx+i, aColIdx] * c[i]
    end

    return dotprod
end

function unsafe_dot(a::T, b::AbstractArray, bLastIdx::Integer) where T
    aLen     = length(a)
    bBaseIdx = bLastIdx - aLen
    @inbounds dotprod  = a[1] * b[bBaseIdx + 1]
    @simd for i in 2:aLen
        @inbounds dotprod += a[i] * b[bBaseIdx + i]
    end

    return dotprod
end

@inline function unsafe_dot(a::Vector{T}, b::Array{T}, bLastIdx::Integer) where T<:BLAS.BlasReal
    BLAS.dot(length(a), pointer(a), 1, pointer(b, bLastIdx - length(a) + 1), 1)
end

function unsafe_dot(a::AbstractVector, b::AbstractVector{T}, c::AbstractVector{T}, cLastIdx::Integer) where T
    aLen    = length(a)
    dotprod = zero(a[1]*b[1])
    @simd for i in 1:aLen-cLastIdx
        @inbounds dotprod += a[i] * b[i+cLastIdx-1]
    end
    @simd for i in 1:cLastIdx
        @inbounds dotprod += a[aLen-cLastIdx+i] * c[i]
    end

    return dotprod
end

# Shifts b into the end a.
# julia> DSP.Util.shiftin!( [1,2,3,4], [5, 6])
# 4-element Array{Int64,1}:
#  3
#  4
#  5
#  6
function shiftin!(a::AbstractVector{T}, b::AbstractVector{T}) where T
    aLen = length(a)
    bLen = length(b)

    if bLen >= aLen
        copyto!(a, 1, b, bLen - aLen + 1, aLen)
    else

        for i in 1:aLen-bLen
            @inbounds a[i] = a[i+bLen]
        end
        bIdx = 1
        for i in aLen-bLen+1:aLen
            @inbounds a[i] = b[bIdx]
            bIdx += 1
        end
    end

    return a
end

# DELAY FINDING UTILITIES

"""
    finddelay(x, y)

Estimate the delay of x with respect to y by locating the peak of their
cross-correlation.

The output delay will be positive when x is delayed with respect y, negative if
advanced, 0 otherwise.

# Example
```jldoctest
julia> finddelay([0, 0, 1, 2, 3], [1, 2, 3])
2

julia> finddelay([1, 2, 3], [0, 0, 1, 2, 3])
-2
```
"""
function finddelay(x::AbstractVector{<: Real}, y::AbstractVector{<: Real})
    s = xcorr(y, x, padmode=:none)
    max_corr = maximum(abs, s)
    max_idxs = findall(x -> abs(x) == max_corr, s)

    center_idx = length(x)
    # Delay is position of peak cross-correlation relative to center.
    # If the maximum cross-correlation is not unique, use the position
    # closest to the center.
    d_ind = argmin(abs.(center_idx .- max_idxs))
    d = center_idx - max_idxs[d_ind]
end

"""
    shiftsignal!(x, s)

Mutating version of shiftsignals(): shift x of s samples and fill the spaces
with zeros in-place.

See also [`shiftsignal`](@ref).
"""
function shiftsignal!(x::AbstractVector, s::Integer)
    l = length(x)
    if abs(s) > l
        error("The absolute value of s must not be greater than the length of x")
    end
    if s > 0
        x[s + 1:l] = x[1:l - s]
        x[1:s] .= 0
    elseif s < 0
        x[1:l + s] = x[1 - s:l]
        x[l + s + 1:l] .= 0
    end
    x
end

"""
    shiftsignal(x, s)

Shift elements of signal x in time by a given amount s of samples and fill
the spaces with zeros. For circular shifting, use circshift.

# Example
```jldoctest
julia> shiftsignal([1, 2, 3], 2)
3-element Array{Int64,1}:
 0
 0
 1

julia> shiftsignal([1, 2, 3], -2)
3-element Array{Int64,1}:
 3
 0
 0
```

See also [`shiftsignal!`](@ref).
"""
shiftsignal(x::AbstractVector, s::Integer) = shiftsignal!(copy(x), s)

"""
    alignsignals!(x, y)

Mutating version of alignsignals(): time align x to y in-place.

See also [`alignsignals`](@ref).
"""
function alignsignals!(x, y)
    d = finddelay(x, y)
    x = shiftsignal!(x, -d)
    x, d
end

"""
    alignsignals(x, y)

Use finddelay() and shiftsignal() to time align x to y. Also return the delay
of x with respect to y.

# Example
```jldoctest
julia> alignsignals([0, 0, 1, 2, 3], [1, 2, 3])
([1, 2, 3, 0, 0], 2)

julia> alignsignals([1, 2, 3], [0, 0, 1, 2, 3])
([0, 0, 1], -2)
```

See also [`alignsignals!`](@ref).
"""
alignsignals(x, y) = alignsignals!(copy(x), y)

end # end module definition


text:
using Test
using SafeTestsets

@safetestset "Polynomials" begin include("polynomials_test.jl") end

@safetestset "HyperRectangle" begin include("hyperrectangle_test.jl") end

@safetestset "Tensor lagrange interpolant" begin include("tensorlaginterp_test.jl") end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
module ImperoHooks

using Impero
using ClimateMachine

using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Elements
import ClimateMachine.Mesh.Elements: baryweights
using GaussQuadrature
using Base.Threads
using Makie, GLMakie, AbstractPlotting
using ImageTransformations, Colors
using AbstractPlotting.MakieLayout, LinearAlgebra

include("kernels.jl")
include("grid.jl")
include("calculus.jl")
include("utils.jl")
include("permutations.jl")
include("find_element.jl")
include("lagrange_interpolation.jl")
include("fields.jl")
include("gridhelper.jl")

end # module

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
@testset "solve_karger.jl" begin

end