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.
=#

#using Debug

"""
Try to hoist allocations outside the loop if possible.
"""
function hoistAllocation(ast::Array{Any,1}, lives, domLoop::DomLoops, state :: expr_state)
    # Only allocations that are not aliased can be safely hoisted.
    # Note that we must rule out simple re-assignment in alias analysis to be conservative about object uniqueness
    # (instead of just variable uniqueness).
    body = CompilerTools.LambdaHandling.getBody(ast, CompilerTools.LambdaHandling.getReturnType(state.LambdaVarInfo))
    uniqSet = AliasAnalysis.from_lambda(state.LambdaVarInfo, body, lives, pir_alias_cb, nothing; noReAssign = true)
    @dprintln(3, "HA: uniqSet = ", uniqSet)
    for l in domLoop.loops
        @dprintln(3, "HA: loop from block ", l.head, " to ", l.back_edge)
        headBlk = lives.cfg.basic_blocks[ l.head ]
        tailBlk = lives.cfg.basic_blocks[ l.back_edge ]
        if length(headBlk.preds) != 2
            continue
        end
        preBlk = nothing
        for blk in headBlk.preds
            if blk.label != tailBlk.label
                preBlk = blk
                break
            end
        end

        #if (is(preBlk, nothing) || length(preBlk.statements) == 0) continue end
        if is(preBlk, nothing) continue end
        tls = lives.basic_blocks[ preBlk ]

        # sometimes the preBlk has no statements
        # in this case we go to preBlk's previous block to find the previous statement of the current loop (for allocations to be inserted)
        while length(preBlk.statements)==0
            if length(preBlk.preds)==1
                preBlk = next(preBlk.preds,start(preBlk.preds))[1]
            end
        end
        if length(preBlk.statements)==0 continue end
        preHead = preBlk.statements[end].index

        head = headBlk.statements[1].index
        tail = tailBlk.statements[1].index
        @dprintln(3, "HA: line before head is ", ast[preHead-1])
        # Is iterating through statement indices this way safe?
        for i = head:tail
            if isAssignmentNode(ast[i]) && isAllocation(ast[i].args[2])
                @dprintln(3, "HA: found allocation at line ", i, ": ", ast[i])
                lhs = ast[i].args[1]
                rhs = ast[i].args[2]
                lhs = toLHSVar(lhs)
                if in(lhs, uniqSet) && (haskey(state.array_length_correlation, lhs))
                    c = state.array_length_correlation[lhs]
                    for (d, v) in state.symbol_array_correlation
                        if v == c
                            ok = true

                            for j = 1:length(d)
                                if !(isa(d[j],Int) || in(d[j], tls.live_out))
                                    ok = false
                                    break
                                end
                            end
                            @dprintln(3, "HA: found correlation dimension ", d, " ", ok, " ", length(rhs.args)-6)
                            if ok && length(rhs.args) - 6 == 2 * length(d) # dimension must match
                                rhs.args = rhs.args[1:6]
                                for s in d
                                    push!(rhs.args, s)
                                    push!(rhs.args, 0)
                                end
                                @dprintln(3, "HA: hoist ", ast[i], " out of loop before line ", head)
                                ast = [ ast[1:preHead-1]; ast[i]; ast[preHead:i-1]; ast[i+1:end] ]
                                break
                            end
                        end
                    end
                end
            end
        end
    end
    return ast
end

function isDeadCall(rhs::Expr, live_out)
    if isCall(rhs)
        fun = getCallFunction(rhs)
        args = getCallArguments(rhs)
        if in(fun, CompilerTools.LivenessAnalysis.wellknown_all_unmodified)
            @dprintln(3, rhs)
            return true
        elseif in(fun, CompilerTools.LivenessAnalysis.wellknown_only_first_modified) &&
                !in(toLHSVar(args[1]), live_out)
            return true
        end
    end
    return false
end

function isDeadCall(rhs::ANY, live_out)
    return false
end

type DictInfo
    live_info
    expr
end

"""
State for the remove_no_deps and insert_no_deps_beginning phases.
"""
type RemoveNoDepsState
    lives             :: CompilerTools.LivenessAnalysis.BlockLiveness
    top_level_no_deps :: Array{Any,1}
    hoistable_scalars :: Set{LHSVar}
    dict_sym          :: Dict{LHSVar, DictInfo}
    change            :: Bool

    function RemoveNoDepsState(l, hs)
        new(l, Any[], hs, Dict{LHSVar, DictInfo}(), false)
    end
end

"""
Works with remove_no_deps below to move statements with no dependencies to the beginning of the AST.
"""
function insert_no_deps_beginning(node, data :: RemoveNoDepsState, top_level_number, is_top_level, read)
    if is_top_level && top_level_number == 1
        return [data.top_level_no_deps; node]
    end
    nothing
end

"""
# This routine gathers up nodes that do not use
# any variable and removes them from the AST into top_level_no_deps.  This works in conjunction with
# insert_no_deps_beginning above to move these statements with no dependencies to the beginning of the AST
# where they can't prevent fusion.
"""
function remove_no_deps(node :: Expr, data :: RemoveNoDepsState, top_level_number, is_top_level, read)
    @dprintln(3,"remove_no_deps starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(3,"remove_no_deps node = ", node, " type = ", typeof(node))
    @dprintln(3,"node.head: ", node.head)
    head = node.head

    if is_top_level
        @dprintln(3,"remove_no_deps is_top_level")

        if head==:gotoifnot
            # Empty the state at the end or begining of a basic block
            data.dict_sym = Dict{LHSVar,DictInfo}()
        end

        live_info = CompilerTools.LivenessAnalysis.find_top_number(top_level_number, data.lives)
        # Remove line number statements.
        if head == :line
            return CompilerTools.AstWalker.ASTWALK_REMOVE
        end
        if live_info == nothing
            @dprintln(3,"remove_no_deps no live_info")
        else
            @dprintln(3,"remove_no_deps live_info = ", live_info)
            @dprintln(3,"remove_no_deps live_info.use = ", live_info.use)

            if isa(node, Number) || isa(node, RHSVar)
                @dprintln(3,"Eliminating dead node: ", node)
                return CompilerTools.AstWalker.ASTWALK_REMOVE
            elseif isAssignmentNode(node)
                @dprintln(3,"Is an assignment node.")
                lhs = node.args[1]
                @dprintln(4,lhs)
                rhs = node.args[2]
                @dprintln(4,rhs)

                if isa(rhs, Expr) && (is(rhs.head, :parfor) || is(rhs.head, :mmap!))
                    # Always keep parfor assignment in order to work with fusion
                    @dprintln(3, "keep assignment due to parfor or mmap! node")
                    return node
                end
                if isa(lhs, RHSVar)
                    lhs_sym = toLHSVar(lhs)
                    @dprintln(3,"remove_no_deps found assignment with lhs symbol ", lhs, " ", rhs, " typeof(rhs) = ", typeof(rhs))
                    # Remove a dead store
                    if !in(lhs_sym, live_info.live_out)
                        data.change = true
                        @dprintln(3,"remove_no_deps lhs is NOT live out")
                        if hasNoSideEffects(rhs) || isDeadCall(rhs, live_info.live_out)
                            @dprintln(3,"Eliminating dead assignment. lhs = ", lhs, " rhs = ", rhs)
                            return CompilerTools.AstWalker.ASTWALK_REMOVE
                        else
                            # Just eliminate the assignment but keep the rhs
                            @dprintln(3,"Eliminating dead variable but keeping rhs, dead = ", lhs_sym)
                            return rhs
                        end
                    else
                        @dprintln(3,"remove_no_deps lhs is live out")
                        if isa(rhs, RHSVar)
                            rhs_sym = toLHSVar(rhs)
                            @dprintln(3,"remove_no_deps rhs is symbol ", rhs_sym)
                            if !in(rhs_sym, live_info.live_out)
                                @dprintln(3,"remove_no_deps rhs is NOT live out")
                                if haskey(data.dict_sym, rhs_sym)
                                    di = data.dict_sym[rhs_sym]
                                    di_live = di.live_info
                                    prev_expr = di.expr

                                    if !in(lhs_sym, di_live.live_out)
                                        prev_expr.args[1] = lhs_sym
                                        delete!(data.dict_sym, rhs_sym)
                                        data.dict_sym[lhs_sym] = DictInfo(di_live, prev_expr)
                                        @dprintln(3,"Lhs is live but rhs is not so substituting rhs for lhs ", lhs_sym, " => ", rhs_sym)
                                        @dprintln(3,"New expr = ", prev_expr)
                                        return CompilerTools.AstWalker.ASTWALK_REMOVE
                                    else
                                        delete!(data.dict_sym, rhs_sym)
                                        @dprintln(3,"Lhs is live but rhs is not.  However, lhs is read between def of rhs and current statement so not substituting.")
                                    end
                                end
                            else
                                @dprintln(3,"Lhs and rhs are live so forgetting assignment ", lhs_sym, " ", rhs_sym)
                                delete!(data.dict_sym, rhs_sym)
                            end
                        else
                            data.dict_sym[lhs_sym] = DictInfo(live_info, node)
                            @dprintln(3,"Remembering assignment for symbol ", lhs_sym, " ", rhs)
                        end
                    end
                end
            else
                @dprintln(3,"Not an assignment node.")
            end

            for j = live_info.use
                delete!(data.dict_sym, j)
            end

            # Here we try to determine which scalar assigns can be hoisted to the beginning of the function.
            #
            # If this statement defines some variable.
            if !isempty(live_info.def)
                @dprintln(3, "Checking if the statement is hoistable.")
                @dprintln(3, "Previous hoistables = ", data.hoistable_scalars)
                # Assume that hoisting is safe until proven otherwise.
                dep_only_on_parameter = true
                # Look at all the variables on which this statement depends.
                # If any of them are not a hoistable scalar then we can't hoist the current scalar definition.
                for i in live_info.use
                    if !in(i, data.hoistable_scalars)
                        @dprintln(3, "Could not hoist because the statement depends on :", i)
                        dep_only_on_parameter = false
                        break
                    end
                end

                # See if there are any calls with side-effects that could prevent moving.
                sews = SideEffectWalkState()
                ParallelAccelerator.ParallelIR.AstWalk(node, hasNoSideEffectWalk, sews)
                if sews.hasSideEffect
                    dep_only_on_parameter = false
                end

                if dep_only_on_parameter
                    # If this statement is defined in more than one place then it isn't hoistable.
                    for i in live_info.def
                        @dprintln(3,"Checking if ", i, " is multiply defined.")
                        @dprintln(4,"data.lives = ", data.lives)
                        if CompilerTools.LivenessAnalysis.countSymbolDefs(i, data.lives) > 1
                            @dprintln(3, "Could not hoist because the function has multiple definitions of: ", i)
                            dep_only_on_parameter = false
                            break
                        end
                    end

                    if dep_only_on_parameter
                        @dprintln(3,"remove_no_deps removing ", node, " because it only depends on hoistable scalars.")
                        push!(data.top_level_no_deps, node)
                        # If the defs in this statement are hoistable then other statements which depend on them may also be hoistable.
                        for i in live_info.def
                            push!(data.hoistable_scalars, i)
                        end
                        return CompilerTools.AstWalker.ASTWALK_REMOVE
                    end
                end
            end
        end
    end

    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function remove_no_deps(node::Union{LabelNode,GotoNode}, data :: RemoveNoDepsState, top_level_number, is_top_level, read)
    if is_top_level
        # Empty the state at the end or begining of a basic block
        data.dict_sym = Dict{LHSVar,DictInfo}()
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function remove_no_deps(node :: LineNumberNode, data :: RemoveNoDepsState, top_level_number, is_top_level, read)
    if is_top_level
        # remove line number nodes
        return CompilerTools.AstWalker.ASTWALK_REMOVE
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function remove_no_deps(node::ANY, data :: RemoveNoDepsState, top_level_number, is_top_level, read)
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

"""
Empty statements can be added to the AST by some passes in ParallelIR.
This pass over the statements of the :body excludes such "nothing" statements from the new :body.
"""
function removeNothingStmts(args :: Array{Any,1}, state)
    newBody = Any[]
    for i = 1:length(args)
        if args[i] != nothing
            push!(newBody, args[i])
        end
    end
    return newBody
end

"""
Holds liveness information for the remove_dead AstWalk phase.
"""
type RemoveDeadState
    lives :: CompilerTools.LivenessAnalysis.BlockLiveness
end

"""
An AstWalk callback that uses liveness information in "data" to remove dead stores.
"""
function remove_dead(node, data :: RemoveDeadState, top_level_number, is_top_level, read)
    @dprintln(3,"remove_dead starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(3,"remove_dead node = ", node, " type = ", typeof(node))
    if typeof(node) == Expr
        @dprintln(3,"node.head = ", node.head)
    end
    ntype = typeof(node)

    if is_top_level
        @dprintln(3,"remove_dead is_top_level")
        live_info = CompilerTools.LivenessAnalysis.find_top_number(top_level_number, data.lives)
        if live_info != nothing
            @dprintln(3,"remove_dead live_info = ", live_info)
            @dprintln(3,"remove_dead live_info.use = ", live_info.use)

            if isAssignmentNode(node)
                @dprintln(3,"Is an assignment node.")
                lhs = node.args[1]
                @dprintln(4,lhs)
                rhs = node.args[2]
                @dprintln(4,rhs)

                if isa(lhs,RHSVar)
                    lhs_sym = toLHSVar(lhs)
                    @dprintln(3,"remove_dead found assignment with lhs symbol ", lhs, " ", rhs, " typeof(rhs) = ", typeof(rhs))
                    # Remove a dead store
                    if !in(lhs_sym, live_info.live_out)
                        @dprintln(3,"remove_dead lhs is NOT live out")
                        if hasNoSideEffects(rhs) || isDeadCall(rhs, live_info.live_out)
                            @dprintln(3,"Eliminating dead assignment. lhs = ", lhs, " rhs = ", rhs)
                            return CompilerTools.AstWalker.ASTWALK_REMOVE
                        else
                            # Just eliminate the assignment but keep the rhs
                            @dprintln(3,"Eliminating dead variable but keeping rhs, dead = ", lhs_sym, " rhs = ", rhs)
                            return rhs
                        end
                    end
                end
            elseif isInvoke(node)
                @dprintln(3,"isInvoke. head = ", node.head, " type = ", typeof(node.args[2]), " name = ", node.args[2].name)
                if hasNoSideEffects(node) || isDeadCall(node, live_info.live_out)
                    @dprintln(3,"Eliminating dead call. node = ", node)
                    return CompilerTools.AstWalker.ASTWALK_REMOVE
                end
            elseif isCall(node)
                @dprintln(3,"isCall. head = ", node.head, " type = ", typeof(node.args[1]), " name = ", node.args[1].name)
                if hasNoSideEffects(node) || isDeadCall(node, live_info.live_out)
                    @dprintln(3,"Eliminating dead call. node = ", node)
                    return CompilerTools.AstWalker.ASTWALK_REMOVE
                end
            end
        end
    end

    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

"""
State to aide in the transpose propagation phase.
"""
type TransposePropagateState
    lives  :: CompilerTools.LivenessAnalysis.BlockLiveness
    transpose_map :: Dict{LHSVar, LHSVar} # transposed output -> matrix in

    function TransposePropagateState(l)
        new(l, Dict{LHSVar, LHSVar}())
    end
end

function transpose_propagate(node :: ANY, data :: TransposePropagateState, top_level_number, is_top_level, read)
    @dprintln(3,"transpose_propagate starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(3,"transpose_propagate node = ", node, " type = ", typeof(node))
    if typeof(node) == Expr
        @dprintln(3,"node.head = ", node.head)
    end
    ntype = typeof(node)

    if is_top_level
        @dprintln(3,"transpose_propagate is_top_level")
        live_info = CompilerTools.LivenessAnalysis.find_top_number(top_level_number, data.lives)

        if live_info != nothing
            # Remove matrices from data.transpose_map if either original or transposed matrix is modified by this statement.
            # For each symbol modified by this statement...
            for def in live_info.def
                @dprintln(4,"Symbol ", def, " is modifed by current statement.")
                # For each transpose map we currently have recorded.
                for mat in data.transpose_map
                    @dprintln(4,"Current mat in data.transpose_map = ", mat)
                    # If original or transposed matrix is modified by the statement.
                    if def == mat[1] || def==mat[2]
                    #@bp
                        @dprintln(3,"transposed or original matrix is modified so removing ", mat," from data.transpose_map.")
                        # Then remove the lhs = rhs entry from copies.
                        delete!(data.transpose_map, mat[1])
                    end
                end
            end
        end

        if isa(node, LabelNode) || isa(node, GotoNode) || (isa(node, Expr) && is(node.head, :gotoifnot))
            # Only transpose propagate within a basic block.  this is now a new basic block.
            empty!(data.transpose_map)
        elseif isAssignmentNode(node) && isCall(node.args[2])
            @dprintln(3,"Is an assignment call node.")
            lhs = toLHSVar(node.args[1])
            rhs = node.args[2]
            func = getCallFunction(rhs)
            if func==GlobalRef(Base,:transpose!)
                @dprintln(3,"transpose_propagate transpose! found.")
                args = getCallArguments(rhs)
                original_matrix = toLHSVar(args[2])
                transpose_var1 = toLHSVar(args[1])
                transpose_var2 = lhs
                data.transpose_map[transpose_var1] = original_matrix
                data.transpose_map[transpose_var2] = original_matrix
            elseif func==GlobalRef(Base,:transpose)
                @dprintln(3,"transpose_propagate transpose found.")
                args = getCallArguments(rhs)
                original_matrix = toLHSVar(args[1])
                transpose_var = lhs
                data.transpose_map[transpose_var] = original_matrix
            elseif func==GlobalRef(Base.LinAlg,:gemm_wrapper!)
                @dprintln(3,"transpose_propagate GEMM found.")
                args = getCallArguments(rhs)
                A = toLHSVar(args[4])
                if haskey(data.transpose_map, A)
                    args[4] = data.transpose_map[A]
                    args[2] = 'T'
                    @dprintln(3,"transpose_propagate GEMM replace transpose arg 1.")
                end
                B = toLHSVar(args[5])
                if haskey(data.transpose_map, B)
                    args[5] = data.transpose_map[B]
                    args[3] = 'T'
                    @dprintln(3,"transpose_propagate GEMM replace transpose arg 2.")
                end
                rhs.args = rhs.head == :invoke ? [ rhs.args[1:2]; args ] : [ rhs.args[1]; args ]
            elseif func==GlobalRef(Base.LinAlg,:gemv!)
                args = getCallArguments(rhs)
                A = toLHSVar(args[3])
                if haskey(data.transpose_map, A)
                    args[3] = data.transpose_map[A]
                    args[2] = 'T'
                end
                rhs.args = rhs.head == :invoke ? [ rhs.args[1:2]; args ] : [ rhs.args[1]; args ]
            # replace arraysize() calls to the transposed matrix with original
            elseif isBaseFunc(func, :arraysize)
                args = getCallArguments(rhs)
                if haskey(data.transpose_map, args[1])
                    args[1] = data.transpose_map[args[1]]
                    if args[2] ==1
                        args[2] = 2
                    elseif args[2] ==2
                        args[2] = 1
                    else
                        throw("transpose_propagate matrix dim error")
                    end
                end
                rhs.args = rhs.head == :invoke ? [ rhs.args[1:2]; args ] : [ rhs.args[1]; args ]
            end
            return node
        end
    end

    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

"""
State to aide in the copy propagation phase.
"""
type CopyPropagateState
    lives  :: CompilerTools.LivenessAnalysis.BlockLiveness
    copies :: Dict{LHSVar, Union{LHSVar,Number}}
    # if ISASSIGNEDONCE flag is set for a variable, its safe to keep it across block boundaries
    safe_copies :: Dict{LHSVar, Union{LHSVar,Number}}
    linfo

    function CopyPropagateState(l, c,s,li)
        new(l,c,s,li)
    end
end

"""
In each basic block, if there is a "copy" (i.e., something of the form "a = b") then put
that in copies as copies[a] = b.  Then, later in the basic block if you see the symbol
"a" then replace it with "b".  Note that this is not SSA so "a" may be written again
and if it is then it must be removed from copies.
"""
function copy_propagate(node :: ANY, data :: CopyPropagateState, top_level_number, is_top_level, read)
    @dprintln(3,"copy_propagate starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(3,"copy_propagate node = ", node, " type = ", typeof(node))
    if typeof(node) == Expr
        @dprintln(3,"node.head = ", node.head)
    end
    ntype = typeof(node)

    if is_top_level
        @dprintln(3,"copy_propagate is_top_level")
        live_info = CompilerTools.LivenessAnalysis.find_top_number(top_level_number, data.lives)

        if live_info != nothing
            # Remove elements from data.copies if the original RHS is modified by this statement.
            # For each symbol modified by this statement...
            for def in live_info.def
                @dprintln(4,"Symbol ", def, " is modifed by current statement.")
                # For each copy we currently have recorded.
                for copy in data.copies
                    @dprintln(4,"Current entry in data.copies = ", copy)
                    # If the rhs of the copy is modified by the statement.
                    if def == copy[2]
                        @dprintln(3,"RHS of data.copies is modified so removing ", copy," from data.copies.")
                        # Then remove the lhs = rhs entry from copies.
                        delete!(data.copies, copy[1])
                    elseif def == copy[1]
                        # LHS is def.  We can maintain the mapping if RHS is dead.
                        if in(copy[2], live_info.live_out)
                            @dprintln(3,"LHS of data.copies is modified and RHS is live so removing ", copy," from data.copies.")
                            # Then remove the lhs = rhs entry from copies.
                            delete!(data.copies, copy[1])
                        end
                    end
                end
            end
        end

        if isa(node, LabelNode) || isa(node, GotoNode) || (isa(node, Expr) && is(node.head, :gotoifnot))
            # Only copy propagate within a basic block.  this is now a new basic block.
            # if ISASSIGNEDONCE flag is set for a variable, its safe to keep it across block boundaries
            data.copies = copy(data.safe_copies)
        elseif isAssignmentNode(node)
            @dprintln(3,"Is an assignment node.")
            # ignore LambdaInfo nodes generated by domain-ir that are essentially dead nodes here
            # TODO: should these nodes be traversed here recursively?
            if isa(node.args[2],LambdaInfo)
                return node
            end
            lhs = AstWalk(node.args[1], copy_propagate, data)
            @dprintln(4,lhs)
            rhs = node.args[2] = AstWalk(node.args[2], copy_propagate, data)
            @dprintln(4,rhs)
            # sometimes lhs can already be replaced with a constant
            if !isa(lhs, RHSVar)
                return node
            end
            node.args[1] = lhs
            if isa(rhs, RHSVar) || (isa(rhs, Number) && !isa(rhs,Complex)) # TODO: fix complex number case
                lhs = toLHSVar(lhs)
                rhs = toLHSVarOrNum(rhs)
                desc = CompilerTools.LambdaHandling.getDesc(lhs, data.linfo)
                if desc & ISASSIGNEDBYINNERFUNCTION != ISASSIGNEDBYINNERFUNCTION
                    @dprintln(3,"Creating copy, lhs = ", lhs, " rhs = ", rhs)
                    # Record that the left-hand side is a copy of the right-hand side.
                    data.copies[lhs] = rhs
                    if (desc & ISASSIGNEDONCE == ISASSIGNEDONCE) &&
                        (isa(rhs, Number) || CompilerTools.LambdaHandling.getDesc(rhs, data.linfo) & ISASSIGNEDONCE == ISASSIGNEDONCE)
                        @dprintln(3,"Creating safe copy, lhs = ", lhs, " rhs = ", rhs)
                        #@bp
                        data.safe_copies[lhs] = rhs
                    end
                end
            end
            return node
        end
    end

    return copy_propagate_helper(node, data, top_level_number, is_top_level, read)
end

function copy_propagate_helper(node::Union{Symbol,RHSVar},
                               data::CopyPropagateState,
                               top_level_number,
                               is_top_level,
                               read)

    lhsVar = toLHSVar(node)
    if haskey(data.copies, lhsVar)
        @dprintln(3,"Replacing ", lhsVar, " with ", data.copies[lhsVar])
        tmp_node = data.copies[lhsVar]
        return isa(tmp_node, Symbol) ? toRHSVar(tmp_node, getType(node, data.linfo), data.linfo) : tmp_node
    end

    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function copy_propagate_helper(node::DomainLambda,
                               data::CopyPropagateState,
                               top_level_number,
                               is_top_level,
                               read)

    @dprintln(3,"Found DomainLambda in copy_propagate, dl = ", node)
    intersection_dict = Dict{LHSVar,Any}()

    # all copies of escaping_defs should be deleted, since there is no guarantee that their values
    # remain the same at when DomainLambda is actually called.
    for v in CompilerTools.LambdaHandling.getEscapingVariables(node.linfo)
        if haskey(data.copies, v) && !haskey(data.safe_copies, v)
            @dprintln(3, "Found escaping_defs for ", v, " remove it from data.copies")
            delete!(data.copies, v)
            @dprintln(3, "data.copies = ", data.copies)
        end
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function copy_propagate_helper(node::ANY,
                               data::CopyPropagateState,
                               top_level_number,
                               is_top_level,
                               read)
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function create_equivalence_classes_assignment(lhs::RHSVar, rhs::RHSVar, state)
    rhs = toLHSVar(rhs)
    lhs = toLHSVar(lhs)

    rhs_corr = getOrAddArrayCorrelation(rhs, state)
    @dprintln(3,"assignment correlation lhs = ", lhs, " type = ", typeof(lhs))
    # if an array has correlation already, there might be a case of multiple assignments
    # in this case, try to make sure sizes are the same or assign a new negative value otherwise
    if haskey(state.array_length_correlation, lhs)
        prev_corr = state.array_length_correlation[lhs]
        prev_size = []
        rhs_size = []
        for (d, v) in state.symbol_array_correlation
            if v==prev_corr
                prev_size = d
            end
            if v==rhs_corr
                rhs_size = d
            end
        end
        if prev_size==[] || rhs_size==[] || prev_size!=rhs_size
            # can't make sure sizes are always equal, assign negative correlation to lhs
            state.array_length_correlation[lhs] = getNegativeCorrelation(state)
            @dprintln(3, "multiple assignment detected, negative correlation assigned for ", lhs)
        end
    else
        lhs_corr = getOrAddArrayCorrelation(toLHSVar(lhs), state)
        merge_correlations(state, lhs_corr, rhs_corr)
        @dprintln(3,"Correlations after assignment merge into lhs")
        print_correlations(3, state)
    end

    CompilerTools.AstWalker.ASTWALK_RECURSE
end

function create_equivalence_classes_assignment(lhs, rhs::Expr, state)
    @dprintln(4,lhs)
    @dprintln(4,rhs)

    if rhs.head == :assertEqShape
        # assertEqShape lets us know that the array mentioned in the assertEqShape node must have the same shape.
        @dprintln(3,"Creating array length assignment from assertEqShape")
        from_assertEqShape(rhs, state)
    elseif rhs.head == :alloc
        # Here an array on the left-hand side is being created from size specification on the right-hand side.
        # Map those array sizes to the corresponding array equivalence class.
        sizes = Any[ x for x in rhs.args[2]]
        n = length(sizes)
        assert(n >= 1 && n <= 3)
        @dprintln(3, "Detected :alloc array allocation. dims = ", sizes)
        checkAndAddSymbolCorrelation(lhs, state, sizes)
    elseif isCall(rhs)
        @dprintln(3, "Detected call rhs in from_assignment.")
        @dprintln(3, "from_assignment call, arg1 = ", rhs.args[1])
        if length(rhs.args) > 1
            @dprintln(3, " arg2 = ", rhs.args[2])
        end
        fun = getCallFunction(rhs)
        args = getCallArguments(rhs)
        if isBaseFunc(fun, :ccall)
            # Same as :alloc above.  Detect an array allocation call and map the specified array sizes to an array equivalence class.
            if args[1] == QuoteNode(:jl_alloc_array_1d)
                dim1 = args[6]
                @dprintln(3, "Detected 1D array allocation. dim1 = ", dim1, " type = ", typeof(dim1))
                checkAndAddSymbolCorrelation(lhs, state, Any[dim1])
            elseif args[1] == QuoteNode(:jl_alloc_array_2d)
                dim1 = args[6]
                dim2 = args[8]
                @dprintln(3, "Detected 2D array allocation. dim1 = ", dim1, " dim2 = ", dim2)
                checkAndAddSymbolCorrelation(lhs, state, Any[dim1, dim2])
            elseif args[1] == QuoteNode(:jl_alloc_array_3d)
                dim1 = args[6]
                dim2 = args[8]
                dim3 = args[10]
                @dprintln(3, "Detected 2D array allocation. dim1 = ", dim1, " dim2 = ", dim2, " dim3 = ", dim3)
                checkAndAddSymbolCorrelation(lhs, state, Any[dim1, dim2, dim3])
            end
        elseif  isBaseFunc(fun, :arraylen)
            # This is the other direction.  Takes an array and extract dimensional information that maps to the array's equivalence class.
            array_param = args[1]                  # length takes one param, which is the array
            assert(isa(array_param, RHSVar))
            array_param_type = CompilerTools.LambdaHandling.getType(array_param, state.LambdaVarInfo) # get its type
            if ndims(array_param_type) == 1            # can only associate when number of dimensions is 1
                dim_symbols = Any[toLHSVar(lhs)]
                @dprintln(3,"Adding symbol correlation from arraylen, name = ", array_param, " dims = ", dim_symbols)
                checkAndAddSymbolCorrelation(toLHSVar(array_param), state, dim_symbols)
            end
        elseif isBaseFunc(fun, :arraysize)
            # This is the other direction.  Takes an array and extract dimensional information that maps to the array's equivalence class.
            if length(args) == 1
                array_param = args[1]                  # length takes one param, which is the array
                assert(isa(array_param, TypedVar))         # should be a TypedVar
                array_param_type = getType(array_param, state.LambdaVarInfo)  # get its type
                array_dims = ndims(array_param_type)
                dim_symbols = Any[]
                for dim_i = 1:array_dims
                    push!(dim_symbols, lhs[dim_i])
                end
                lhsVar = toLHSVar(args[1])
                @dprintln(3,"Adding symbol correlation from arraysize, name = ", lhsVar, " dims = ", dim_symbols)
                checkAndAddSymbolCorrelation(lhsVar, state, dim_symbols)
            elseif length(args) == 2
                @dprintln(1,"Can't establish symbol to array length correlations yet in the case where dimensions are extracted individually.")
            else
                throw(string("arraysize AST node didn't have 2 or 3 arguments."))
            end
        elseif isBaseFunc(fun, :reshape)
            # rhs.args[2] is the array to be reshaped, lhs is the result, rhs.args[3] is a tuple with new shape
            if haskey(state.tuple_table, args[2])
                checkAndAddSymbolCorrelation(lhs, state, state.tuple_table[args[2]])
            end
        elseif isBaseFunc(fun, :tuple)
            ok = true
            for s in args
                if !(isa(s,TypedVar) || isa(s,Int))
                    ok = false
                end
            end
            if ok
                state.tuple_table[lhs]=args[1:end]
            end
        end
    elseif rhs.head == :mmap! || rhs.head == :mmap || rhs.head == :map! || rhs.head == :map
        # Arguments to these domain operations implicit assert that equality of sizes so add/merge equivalence classes for the arrays to this operation.
        rhs_corr = extractArrayEquivalencies(rhs, state)
        @dprintln(3,"lhs = ", lhs, " type = ", typeof(lhs))
        if rhs_corr != nothing && isa(lhs, RHSVar)
            lhs = toLHSVar(lhs)
            # if an array has correlation already, there might be a case of multiple assignments
            # in this case, try to make sure sizes are the same or assign a new negative value otherwise
            if haskey(state.array_length_correlation, lhs)
                prev_corr = state.array_length_correlation[lhs]
                prev_size = []
                rhs_size = []
                for (d, v) in state.symbol_array_correlation
                    if v==prev_corr
                        prev_size = d
                    end
                    if v==rhs_corr
                        rhs_size = d
                    end
                end
                if prev_size==[] || rhs_size==[] || prev_size!=rhs_size
                    # can't make sure sizes are always equal, assign negative correlation to lhs
                    state.array_length_correlation[lhs] = getNegativeCorrelation(state)
                    @dprintln(3, "multiple assignment detected, negative correlation assigned for ", lhs)
                end
            else
                lhs_corr = getOrAddArrayCorrelation(toLHSVar(lhs), state)
                merge_correlations(state, lhs_corr, rhs_corr)
                @dprintln(3,"Correlations after map merge into lhs")
                print_correlations(3, state)
            end
        end
    end
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function create_equivalence_classes_assignment(lhs, rhs::ANY, state)
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function getNegativeCorrelation(state)
    state.multi_correlation -= 1
    return state.multi_correlation
end

function print_correlations(level, state)
    if !isempty(state.array_length_correlation)
        dprintln(level,"array_length_correlations = ", state.array_length_correlation)
    end
    if !isempty(state.symbol_array_correlation)
        dprintln(level,"symbol_array_correlations = ", state.symbol_array_correlation)
    end
    if !isempty(state.range_correlation)
        dprintln(level,"range_correlations = ")
        for i in state.range_correlation
            dprint(level, "    ")
            for j in i[1]
                if isa(j, RangeData)
                    dprint(level, j.exprs, " ")
                else
                    dprint(level, j, " ")
                end
            end
            dprintln(level, " => ", i[2])
        end
    end
end

"""
AstWalk callback to determine the array equivalence classes.
"""
function create_equivalence_classes(node :: Expr, state :: expr_state, top_level_number :: Int64, is_top_level :: Bool, read :: Bool)
    @dprintln(3,"create_equivalence_classes starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(3,"create_equivalence_classes node = ", node, " type = ", typeof(node))
    @dprintln(3,"node.head: ", node.head)
    print_correlations(3, state)

    if node.head == :lambda
        save_LambdaVarInfo  = state.LambdaVarInfo
        linfo, body = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(node)
        state.LambdaVarInfo = linfo
        AstWalk(body, create_equivalence_classes, state)
        state.LambdaVarInfo = save_LambdaVarInfo
        return node
    end

    # We can only extract array equivalences from top-level statements.
    if is_top_level
        @dprintln(3,"create_equivalence_classes is_top_level")

        if isAssignmentNode(node)
            # Here the node is an assignment.
            @dprintln(3,"Is an assignment node.")
            # return value here since this function can replace arraysize() calls
            return create_equivalence_classes_assignment(toLHSVar(node.args[1]), node.args[2], state)
        else
            if node.head == :mmap! || node.head == :mmap || node.head == :map! || node.head == :map
                extractArrayEquivalencies(node, state)
            end
        end
    end

    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function create_equivalence_classes(node :: ANY, state :: expr_state, top_level_number :: Int64, is_top_level :: Bool, read :: Bool)
    @dprintln(3,"create_equivalence_classes starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(3,"create_equivalence_classes node = ", node, " type = ", typeof(node))
    @dprintln(3,"Not an assignment or expr node.")
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

"""
Given an array whose name is in "x", allocate a new equivalence class for this array.
"""
function addUnknownArray(x :: LHSVar, state :: expr_state)
    @dprintln(3, "addUnknownArray x = ", x, " next = ", state.next_eq_class)
    m = state.next_eq_class
    state.next_eq_class += 1
    state.array_length_correlation[x] = m + 1
end

"""
Given an array of RangeExprs describing loop nest ranges, allocate a new equivalence class for this range.
"""
function addUnknownRange(x :: Array{DimensionSelector,1}, state :: expr_state)
    m = state.next_eq_class
    state.next_eq_class += 1
    state.range_correlation[x] = m + 1
end

"""
If we somehow determine that two sets of correlations are actually the same length then merge one into the other.
"""
function merge_correlations(state, unchanging, eliminate)
    if unchanging < 0 || eliminate < 0
        @dprintln(3,"merge_correlations will not merge because ", unchanging, " and/or ", eliminate, " represents an array that is multiply defined within the function.")
        return unchanging
    end

    # For each array in the dictionary.
    for i in state.array_length_correlation
        # If it is in the "eliminate" class...
        if i[2] == eliminate
            # ...move it to the "unchanging" class.
            state.array_length_correlation[i[1]] = unchanging
        end
    end
    # The symbol_array_correlation shares the equivalence class space so
    # do the same re-numbering here.
    for i in state.symbol_array_correlation
        if i[2] == eliminate
            state.symbol_array_correlation[i[1]] = unchanging
        end
    end
    # The range_correlation shares the equivalence class space so
    # do the same re-numbering here.
    for i in state.range_correlation
        if i[2] == eliminate
            state.range_correlation[i[1]] = unchanging
        end
    end

    return unchanging
end

"""
If we somehow determine that two arrays must be the same length then
get the equivalence classes for the two arrays and merge those equivalence classes together.
"""
function add_merge_correlations(old_sym :: LHSVar, new_sym :: LHSVar, state :: expr_state)
    @dprintln(3, "add_merge_correlations ", old_sym, " ", new_sym)
    print_correlations(3, state)
    old_corr = getOrAddArrayCorrelation(old_sym, state)
    new_corr = getOrAddArrayCorrelation(new_sym, state)
    ret = merge_correlations(state, old_corr, new_corr)
    @dprintln(3, "add_merge_correlations post")
    print_correlations(3, state)

    return ret
end

"""
Return a correlation set for an array.  If the array was not previously added then add it and return it.
"""
function getOrAddArrayCorrelation(x :: LHSVar, state :: expr_state)
    if !haskey(state.array_length_correlation, x)
        @dprintln(3,"Correlation for array not found = ", x)
        addUnknownArray(x, state)
    end
    state.array_length_correlation[x]
end

"""
"node" is a domainIR node.  Take the arrays used in this node, create an array equivalence for them if they
don't already have one and make sure they all share one equivalence class.
"""
function extractArrayEquivalencies(node :: Expr, state)
    input_args = node.args

    # Make sure we get what we expect from domain IR.
    # There should be two entries in the array, another array of input array symbols and a DomainLambda type
    if(length(input_args) < 2)
        throw(string("extractArrayEquivalencies input_args length should be at least 2 but is ", length(input_args)))
    end

    # First arg is an array of input arrays to the mmap!
    input_arrays = input_args[1]
    len_input_arrays = length(input_arrays)
    @dprintln(2,"Number of input arrays: ", len_input_arrays)
    @dprintln(3,"input_arrays =  ", input_arrays)
    assert(len_input_arrays > 0)

    # Second arg is a DomainLambda
    ftype = typeof(input_args[2])
    @dprintln(2,"extractArrayEquivalencies function = ",input_args[2])
    if(ftype != DomainLambda)
        throw(string("extractArrayEquivalencies second input_args should be a DomainLambda but is of type ", typeof(input_args[2])))
    end

#    if !isa(input_arrays[1], RHSVar)
#        @dprintln(1, "extractArrayEquivalencies input_arrays[1] is not RHSVar")
#        return nothing
#    end

    inputInfo = InputInfo[]
    for i = 1 : length(input_arrays)
        push!(inputInfo, get_mmap_input_info(input_arrays[i], state))
    end
#    num_dim_inputs = findSelectedDimensions(inputInfo, state)
    @dprintln(3, "inputInfo = ", inputInfo)

    main_length_correlation = getCorrelation(inputInfo[1], state)
    # Get the correlation set of the first input array.
    #main_length_correlation = getOrAddArrayCorrelation(toLHSVar(input_arrays[1]), state)

    # Make sure each input array is a TypedVar
    # Also, create indexed versions of those symbols for the loop body
    for i = 2:length(inputInfo)
        @dprintln(3,"extractArrayEquivalencies input_array[i] = ", input_arrays[i], " type = ", typeof(input_arrays[i]))
        this_correlation = getCorrelation(inputInfo[i], state)
        # Verify that all the inputs are the same size by verifying they are in the same correlation set.
        if this_correlation != main_length_correlation
            merge_correlations(state, main_length_correlation, this_correlation)
        end
    end

    @dprintln(3,"extractArrayEq result")
    print_correlations(3, state)
    return main_length_correlation
end

"""
Make sure all the dimensions are TypedVars or constants.
Make sure each dimension variable is assigned to only once in the function.
Extract just the dimension variables names into dim_names and then register the correlation from lhs to those dimension names.
"""
function checkAndAddSymbolCorrelation(lhs :: LHSVar, state, dim_array)
    dim_names = Union{RHSVar,Int}[]

    for i = 1:length(dim_array)
        # constant sizes are either TypedVars, Symbols or Ints, TODO: expand to GenSyms that are constant
        if !(isa(dim_array[i],RHSVar) || isa(dim_array[i], Int))
            @dprintln(3, "checkAndAddSymbolCorrelation dim not Int or RHSVar ", dim_array[i])
            return false
        end
        dim_array[i] = toLHSVar(dim_array[i])
        desc = 0
        if isa(dim_array[i],Symbol) && !(CompilerTools.LambdaHandling.getType(dim_array[i], state.LambdaVarInfo)<:Int)
            @dprintln(3, "checkAndAddSymbolCorrelation dim symbol not Int ", dim_array[i])
            throw(string("Dimension not an Int"))
        end
        if !isa(dim_array[i],Int)
            desc = CompilerTools.LambdaHandling.getDesc(dim_array[i], state.LambdaVarInfo)
        end
        # FIXME: description of input parameters not changed in function is always 0?
        if !isa(dim_array[i],Int) && ((desc & ISASSIGNED == ISASSIGNED) && !(desc & ISASSIGNEDONCE == ISASSIGNEDONCE))
            @dprintln(3, "checkAndAddSymbolCorrelation dim not Int or assigned once ", dim_array[i])
            return false
        end
        push!(dim_names, dim_array[i])
    end

    @dprintln(3, "Will establish array length correlation for const size lhs = ", lhs, " dims = ", dim_names)
    getOrAddSymbolCorrelation(lhs, state, dim_names)
    return true
end

"""
Gets (or adds if absent) the range correlation for the given array of RangeExprs.
"""
function getOrAddRangeCorrelation(array, ranges :: Array{DimensionSelector,1}, state :: expr_state)
    @dprintln(3, "getOrAddRangeCorrelation for ", array, " with ranges = ", ranges)
    if print_times
        @dprintln(3, "with hash = ", hash(ranges))
    end
    print_correlations(3, state)

    # We can't match on array of RangeExprs so we flatten to Array of Any
    all_mask = true
    for i = 1:length(ranges)
        all_mask = all_mask & isa(ranges[i], MaskSelector)
    end

    if !haskey(state.range_correlation, ranges)
        @dprintln(3,"Exact match for correlation for range not found = ", ranges)
        # Look for an equivalent but non-exact range in the dictionary.
        nonExactCorrelation = nonExactRangeSearch(ranges, state.range_correlation)
        if nonExactCorrelation == nothing
            @dprintln(3, "No non-exact match so adding new range")
            range_corr = addUnknownRange(ranges, state)
            # If all the dimensions are selected based on masks then the iteration space
            # is that of the entire array and so we can establish a correlation between the
            # DimensionSelector and the whole array.
            if all_mask
                masked_array_corr = getOrAddArrayCorrelation(toLHSVar(array), state)
                @dprintln(3, "All dimension selectors are masks so establishing correlation to main array ", masked_array_corr)
                range_corr = merge_correlations(state, masked_array_corr, range_corr)

                if length(ranges) == 1
                    print_correlations(3, state)
                    mask_correlation = getCorrelation(ranges[1].value, state)

                    @dprintln(3, "Range length is 1 so establishing correlation between range ", range_corr, " and the mask ", ranges[1].value, " with correlation ", mask_correlation)
                    range_corr = merge_correlations(state, mask_correlation, range_corr)
                end
            end
        else
            # Found an equivalent range.
            @dprintln(3, "Adding non-exact range match to class ", nonExactCorrelation)
            state.range_correlation[ranges] = nonExactCorrelation
        end
        @dprintln(3, "getOrAddRangeCorrelation final correlations")
        print_correlations(3, state)
    end
    state.range_correlation[ranges]
end

"""
A new array is being created with an explicit size specification in dims.
"""
function getOrAddSymbolCorrelation(array :: LHSVar, state :: expr_state, dims :: Array{Union{RHSVar,Int},1})
    if !haskey(state.symbol_array_correlation, dims)
        # We haven't yet seen this combination of dims used to create an array.
        @dprintln(3,"Correlation for symbol set not found, dims = ", dims)
        if haskey(state.array_length_correlation, array)
            return state.symbol_array_correlation[dims] = state.array_length_correlation[array]
        else
            # Create a new array correlation number for this array and associate that number with the dim sizes.
            return state.symbol_array_correlation[dims] = addUnknownArray(array, state)
        end
    else
        @dprintln(3,"Correlation for symbol set found, dims = ", dims)
        # We have previously seen this combination of dim sizes used to create an array so give the new
        # array the same array length correlation number as the previous one.
        return state.array_length_correlation[array] = state.symbol_array_correlation[dims]
    end
end

"""
Replace arraysize() calls for arrays with known constant sizes.
Constant size is Int constants, as well as assigned once variables which are
in symbol_array_correlation. Variables should be assigned before current statement, however.
"""
function replaceConstArraysizes(node :: Expr, state::expr_state, top_level_number :: Int64, is_top_level :: Bool, read :: Bool)
    @dprintln(4,"replaceConstArraysizes starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(4,"replaceConstArraysizes node = ", node, " type = ", typeof(node))
    @dprintln(4,"node.head: ", node.head)
    print_correlations(3, state)

    # TODO: handle arraylen similarly

    live_info = CompilerTools.LivenessAnalysis.find_top_number(top_level_number, state.block_lives)

    if isCall(node) && isBaseFunc(getCallFunction(node), :arraysize)
        # replace arraysize calls when size is known and constant
        args = getCallArguments(node)
        arr = toLHSVar(args[1])
        if isa(args[2],Int) && haskey(state.array_length_correlation, arr)
            arr_class = state.array_length_correlation[arr]
            for (d, v) in state.symbol_array_correlation
                if v==arr_class
                    res = d[args[2]]
                    # only replace when the size is constant or a valid live variable
                    # check def since a symbol correlation might be defined with current arraysize() in reverse direction
                    if isIntType(res) || ( live_info!=nothing && in(res, live_info.live_in) && !in(res,live_info.def) )
                        @dprintln(3, "arraysize() replaced: ", node," res ",res)
                        return res
                    end
                end
            end
        end
    end

    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function replaceConstArraysizes(node :: ANY, state :: expr_state, top_level_number :: Int64, is_top_level :: Bool, read :: Bool)
    @dprintln(4,"replaceConstArraysizes starting top_level_number = ", top_level_number, " is_top = ", is_top_level)
    @dprintln(4,"replaceConstArraysizes node = ", node, " type = ", typeof(node))
    @dprintln(4,"Not an expr node.")
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

"""
Implements one of the main ParallelIR passes to remove assertEqShape AST nodes from the body if they are statically known to be in the same equivalence class.
"""
function removeAssertEqShape(args :: Array{Any,1}, state)
    newBody = Any[]
    for i = 1:length(args)
        # Add the current statement to the new body unless the statement is an assertEqShape Expr and the array in the assertEqShape are known to be the same size.
        if !(typeof(args[i]) == Expr && args[i].head == :assertEqShape && from_assertEqShape(args[i], state))
            push!(newBody, args[i])
        end
    end
    return newBody
end

"""
Create array equivalences from an assertEqShape AST node.
There are two arrays in the args to assertEqShape.
"""
function from_assertEqShape(node::Expr, state)
    @dprintln(3,"from_assertEqShape ", node)
    a1 = node.args[1]    # first array to compare
    a2 = node.args[2]    # second array to compare
    a1_corr = getOrAddArrayCorrelation(toLHSVar(a1), state)  # get the length set of the first array
    a2_corr = getOrAddArrayCorrelation(toLHSVar(a2), state)  # get the length set of the second array
    if a1_corr == a2_corr
        # If they are the same then return an empty array so that the statement is eliminated.
        @dprintln(3,"assertEqShape statically verified and eliminated for ", a1, " and ", a2)
        return true
    else
        @dprintln(3,"a1 = ", a1, " ", a1_corr, " a2 = ", a2, " ", a2_corr, " correlations")
        print_correlations(3, state)
        # If assertEqShape is called on e.g., inputs, then we can't statically eliminate the assignment
        # but if the assert doesn't fire then we do thereafter know that the arrays are in the same length set.
        merge_correlations(state, a1_corr, a2_corr)
        @dprintln(3,"assertEqShape NOT statically verified.  Merge correlations")
        print_correlations(3, state)
        return false
    end
end


text:
import BlockArrays

"""
    FieldVector

A `FieldVector` is a wrapper around one or more `Field`s that acts like vector
of the underlying arrays.

It is similar in spirit to [`ArrayPartition` from
RecursiveArrayTools.jl](https://github.com/SciML/RecursiveArrayTools.jl#arraypartition),
but allows referring to fields by name.

# Constructors

    FieldVector(;name1=field1, name2=field2, ...)

Construct a `FieldVector`, wrapping `field1, field2, ...` using the names
`name1, name2, ...`.
"""
struct FieldVector{T, M} <: BlockArrays.AbstractBlockVector{T}
    values::M
end
FieldVector{T}(values::M) where {T, M} = FieldVector{T, M}(values)

"""
    Fields.ScalarWrapper(val) <: AbstractArray{T,0}

This is a wrapper around scalar values that allows them to be mutated as part of
a FieldVector. A call `getproperty` on a `FieldVector` with this component will
return a scalar, instead of the boxed object.
"""
mutable struct ScalarWrapper{T} <: AbstractArray{T, 0}
    val::T
end
Base.size(::ScalarWrapper) = ()
Base.getindex(s::ScalarWrapper, args...) = s.val
Base.setindex!(s::ScalarWrapper, value, args...) = s.val = value
Base.similar(s::ScalarWrapper) = ScalarWrapper(s.val)

"""
    Fields.wrap(x)

Construct a mutable wrapper around `x`. This can be extended for new types
(especially immutable ones).
"""
wrap(x) = x
wrap(x::Real) = ScalarWrapper(x)
wrap(x::NamedTuple) = FieldVector(; pairs(x)...)

"""
    Fields.unwrap(x::T)

This is called when calling `getproperty` on a `FieldVector` property of element
type `T`.
"""
unwrap(x) = x
unwrap(x::ScalarWrapper) = x[]

function FieldVector(; kwargs...)
    values = map(wrap, NamedTuple(kwargs))
    T = promote_type(
        map(RecursiveArrayTools.recursive_bottom_eltype, values)...,
    )
    return FieldVector{T}(values)
end

_values(fv::FieldVector) = getfield(fv, :values)

"""
    backing_array(x)

The `AbstractArray` that is backs an object `x`, allowing it to be treated as a
component of a `FieldVector`.
"""
backing_array(x) = x
backing_array(x::Field) = parent(x)

Base.propertynames(fv::FieldVector) = propertynames(_values(fv))
@inline function Base.getproperty(fv::FieldVector, name::Symbol)
    unwrap(getfield(_values(fv), name))
end

@inline function Base.setproperty!(fv::FieldVector, name::Symbol, value)
    x = getfield(_values(fv), name)
    x .= value
end

BlockArrays.blockaxes(fv::FieldVector) =
    (BlockArrays.BlockRange(1:length(_values(fv))),)
Base.axes(fv::FieldVector) =
    (BlockArrays.blockedrange(map(length ∘ backing_array, Tuple(_values(fv)))),)

Base.getindex(fv::FieldVector, block::BlockArrays.Block{1}) =
    backing_array(_values(fv)[block.n...])
function Base.getindex(fv::FieldVector, bidx::BlockArrays.BlockIndex{1})
    X = fv[BlockArrays.block(bidx)]
    X[bidx.α...]
end
Base.getindex(fv::FieldVector, i::Integer) =
    getindex(fv, BlockArrays.findblockindex(axes(fv, 1), i))

function Base.setindex!(fv::FieldVector, val, bidx::BlockArrays.BlockIndex{1})
    X = fv[BlockArrays.block(bidx)]
    X[bidx.α...] = val
end
Base.setindex!(fv::FieldVector, val, i::Integer) =
    setindex!(fv, val, BlockArrays.findblockindex(axes(fv, 1), i))

Base.similar(fv::FieldVector{T}) where {T} =
    FieldVector{T}(map(similar, _values(fv)))
Base.similar(fv::FieldVector{T}, ::Type{T}) where {T} =
    FieldVector{T}(map(similar, _values(fv)))
_similar(x, ::Type{T}) where {T} = similar(x, T)
_similar(x::Field, ::Type{T}) where {T} =
    Field(DataLayouts.replace_basetype(field_values(x), T), axes(x))
Base.similar(fv::FieldVector{T}, ::Type{T′}) where {T, T′} =
    FieldVector{T′}(map(x -> _similar(x, T′), _values(fv)))

Base.copy(fv::FieldVector{T}) where {T} = FieldVector{T}(map(copy, _values(fv)))
Base.zero(fv::FieldVector{T}) where {T} = FieldVector{T}(map(zero, _values(fv)))

# Ensure that Field axes are preserved when FieldVectors are deepcopied and when
# objects that contain FieldVectors are deepcopied.
_deepcopy_internal(x, stackdict::IdDict) = Base.deepcopy_internal(x, stackdict)
_deepcopy_internal(x::Field, stackdict::IdDict) =
    Field(Base.deepcopy_internal(field_values(x), stackdict), axes(x))
Base.deepcopy_internal(fv::FieldVector{T}, stackdict::IdDict) where {T} =
    FieldVector{T}(map(x -> _deepcopy_internal(x, stackdict), _values(fv)))

struct FieldVectorStyle <: Base.Broadcast.AbstractArrayStyle{1} end

Base.Broadcast.BroadcastStyle(::Type{<:FieldVector}) = FieldVectorStyle()

Base.Broadcast.BroadcastStyle(
    fs::FieldVectorStyle,
    as::Base.Broadcast.DefaultArrayStyle{0},
) = fs
Base.Broadcast.BroadcastStyle(
    fs::FieldVectorStyle,
    as::Base.Broadcast.AbstractArrayStyle{0},
) = fs
Base.Broadcast.BroadcastStyle(
    fs::FieldVectorStyle,
    as::Base.Broadcast.DefaultArrayStyle,
) = as
Base.Broadcast.BroadcastStyle(
    fs::FieldVectorStyle,
    as::Base.Broadcast.AbstractArrayStyle,
) = as

function Base.similar(
    bc::Base.Broadcast.Broadcasted{FieldVectorStyle},
    ::Type{T},
) where {T}
    for arg in bc.args
        if arg isa FieldVector ||
           arg isa Base.Broadcast.Broadcasted{FieldVectorStyle}
            return similar(arg, T)
        end
    end
    error("Cannot construct FieldVector")
end

Base.mapreduce(f, op, fv::FieldVector) =
    mapreduce(x -> mapreduce(f, op, backing_array(x)), op, _values(fv))

Base.any(f, fv::FieldVector) = any(x -> any(f, backing_array(x)), _values(fv))
Base.any(f::Function, fv::FieldVector) = # avoid ambiguities
    any(x -> any(f, backing_array(x)), _values(fv))
Base.any(fv::FieldVector) = any(identity, A)

Base.all(f, fv::FieldVector) = all(x -> all(f, backing_array(x)), _values(fv))
Base.all(f::Function, fv::FieldVector) =
    all(x -> all(f, backing_array(x)), _values(fv))
Base.all(fv::FieldVector) = all(identity, fv)

# TODO: figure out a better way to handle these
# https://github.com/JuliaArrays/BlockArrays.jl/issues/185
LinearAlgebra.ldiv!(
    x::FieldVector,
    A::LinearAlgebra.QRCompactWY,
    b::FieldVector,
) = x .= LinearAlgebra.ldiv!(A, Vector(b))
LinearAlgebra.ldiv!(A::LinearAlgebra.QRCompactWY, x::FieldVector) =
    x .= LinearAlgebra.ldiv!(A, Vector(x))

LinearAlgebra.ldiv!(x::FieldVector, A::LinearAlgebra.LU, b::FieldVector) =
    x .= LinearAlgebra.ldiv!(A, Vector(b))

LinearAlgebra.ldiv!(A::LinearAlgebra.LU, x::FieldVector) =
    x .= LinearAlgebra.ldiv!(A, Vector(x))

ext: jl


text:
"""
    manufacturer(x)
    manufacturer!(x, val)

The manufacturer of the equipment that produced the composite instances.
"""
@defprop Manufacturer{:manufacturer}::String

"""
    manufacturer_model_name(x)
    manufacturer_model_name!(x, val)

The manufacturer's model name of the equipment that produced the composite instances.
"""
@defprop ManufacturerModelName{:manufacturer_model_name}::String

"""
    device_serial_number(x)
    device_serial_number!(x, val)

The serial number of the equipment that produced the composite instances. A
pseudonym can also be used to prevent the equipment from being identifiable,
so long as each pseudonym is unique within the dataset.
"""
@defprop DeviceSerialNumber{:device_serial_number}::String

"""
    software_versions(x)
    software_versions!(x, val)

The manufacturer’s designation of software version of the equipment that produced the composite instances.
"""
@defprop SoftwareVersions{:software_versions}::String

"""
    station_name(x)
    station_name!(x, val)

Institution defined name of the machine that produced the composite instances.
"""
@defprop StationName{:station_name}::String

"""
    HardwareMetadata

Metadata structure for general MRI sequence information.

## Supported Properties
$(description_list(manufacturer, manufacturer_model_name, device_serial_number))

## Examples
```jldoctest
julia> using NeuroCore

julia> using NeuroCore.NeuroMetadata

julia> m = HardwareMetadata("a", "b", "c");

julia> m.device_serial_number
"c"

julia> m.manufacturer_model_name
"b"

julia> m.manufacturer
"a"
```
"""
struct HardwareMetadata{M} <: AbstractPropertyList{M}
    manufacturer::String
    model::String
    serial_number::String
    extension::M
end

FieldProperties.proplist(m::HardwareMetadata) = getfield(m, :extension)

function HardwareMetadata(manufacturer, model, serial_number; kwargs...)
    return HardwareMetadata(manufacturer, model, serial_number, PropertyList(; kwargs...))
end

@properties HardwareMetadata begin
    manufacturer(self) => :manufacturer
    manufacturer_model_name(self) => :model
    device_serial_number(self) => :serial_number
    Any(self) => :extension
    Any!(self, val) => :extension
end


text:
# Input/output of tensors

"""
    read(file) -> C, rho

Return the elasticity matrix `C` and density `rho` from the ecs-format file
`file` (which usually has the extension `.ecs`).  This should be in the format:

    1 1 C11
    . . ...
    i j Cij
    . . ...
    6 6 C66
    7 7 rho

Comments are permitted at any point, preceded by a '#' character.

In this version of the code, all elastic constants must be specified, even if
they are 0.
"""
function read(file::AbstractString)
    isfile(file) || error("file \"$file\" does not exist.")
    d = readdlm(file)
    size(d, 1) == 22 || error("file \"$file\" is not in expected format")
    C = zero(EC)
    rho = 0
    for l in 1:22
        i = Int(d[l,1])
        j = Int(d[l,2])
        if i == j == 7
            rho = d[l,3]
        else
            1 <= i <= j <= 6 || error("unexpected indices '$i' and '$j' in file \"$file\"")
            C[i,j] = d[l,3]
            C[j,i] = d[l,3]
        end
    end
    C, rho
end

"""
    write(C, rho, file, comment="")

Write the density-normalised elastic constants `C` and density `rho` to the file
`file` in the following format:

    1 1 C11
    . . ...
    i j Cij
    . . ...
    6 6 C66
    7 7 rho

If supplied, `comment` should be a string of one or more lines to write at the
end of the file to describe the constants.  By default, the following is written:

    # Saved by <user> on <hostname> on <date> using <Julia version>

Such files usually have a `.ecs` file extension.
"""
function Base.write(C::Union{EC,Array{<:Number,2}}, rho::Number, file, comment::AbstractString="")
    is_6x6(C) || throw(ArgumentError("elastic constants must be an EC or 6x6 matrix"))
    is_stable(C) && is_symm(C) ||
        warn("elastic constants are not in the right form.  "
              * "May be asymmetric or unstable.")
    isdir(dirname(file)) || error("directory \"$(dirname(file))\" does not exist")
    open(file, "w") do f
        for i = 1:6, j=i:6
            @printf(f, "%d %d %12.6e\n", i, j, C[i,j]*rho)
        end
        @printf(f, "7 7 %9.3f\n", rho)
        if comment != ""
            if comment[1] != '#'
                comment = "# " * comment
            end
            write(f, replace(chomp(comment), "\n"=>"\n# ") * "\n")
        end
        user = ENV["USER"]
        hostname = readchomp(`hostname`)
        write(f, "# Saved by user $user on $hostname on $(now()) using CIJ.jl on Julia $VERSION\n")
    end
    nothing
end


text:
module Shape
    abstract type T end
    area(shape::T) = throw(MethodError(area, shape))
    combined_area(a::T, b::T) = area(a) + area(b)
end

module Circle
    import ..Shape

    struct T <: Shape.T
        diameter::Float64
    end
    radius(c::T) = c.diameter / 2
    Shape.area(c::T) = π * radius(c) ^ 2
end

module AbstractRectangle
    import ..Shape

    abstract type T <: Shape.T end
    width(rectangle::T) = throw(MethodError(width, rectangle))
    height(rectangle::T) = throw(MethodError(width, rectangle))
    Shape.area(r::T) = width(r) * height(r)
end

module Rectangle
    import ..AbstractRectangle

    struct T <: AbstractRectangle.T
        width::Float64
        height::Float64
    end
    AbstractRectangle.width(r::T) = r.width
    AbstractRectangle.height(r::T) = r.height
end

module Square
    import ..AbstractRectangle

    struct T <: AbstractRectangle.T
        length::Float64
    end
    AbstractRectangle.width(s::T) = s.length
    AbstractRectangle.height(s::T) = s.length
end

c = Circle.T(3)
s = Square.T(3)
r = Rectangle.T(3, 2)

@show Shape.combined_area(c, s)
@show Shape.combined_area(s, r)

max_stars_repo_name: ninjaaron/dispatching-design-patterns

max_issues_repo_name: ninjaaron/dispatching-design-patterns

max_forks_repo_name: ninjaaron/dispatching-design-patterns


text:
"""
    new_trace = elliptical_slice(trace, addr, mu, cov)

Apply an elliptical slice sampling update to a given random choice with a multivariate normal prior.

Also takes the mean vector and covariance matrix of the prior.

[Reference URL](http://proceedings.mlr.press/v9/murray10a/murray10a.pdf)
"""
function elliptical_slice(trace, addr, mu, cov)
    args = get_args(trace)
    argdiffs = map((_) -> NoChange(), args)

    # sample nu
    nu = mvnormal(zeros(length(mu)), cov)

    # sample u
    u = uniform(0, 1)

    # sample initial theta and define bracket
    theta = uniform(0, 2*pi)
    theta_min, theta_max = (theta - 2*pi, theta)

    # previous value
    f = trace[addr] .- mu

    new_f = f * cos(theta) + nu * sin(theta)
    new_trace, weight = update(trace, args, argdiffs, choicemap((addr, new_f .+ mu)))
    while weight <= log(u)
        if theta < 0
            theta_min = theta
        else
            theta_max = theta
        end
        theta = uniform(theta_min, theta_max)
        new_f = f * cos(theta) + nu * sin(theta)
        new_trace, weight = update(trace, args, argdiffs, choicemap((addr, new_f .+ mu)))
    end
    return new_trace
end

export elliptical_slice


text:
@testset "Slices" begin
    xs = Slice{Int,10}()

    @test_throws BoundsError xs[0]
    @test_throws BoundsError xs[11]
    @test_throws BoundsError xs[0] = 1
    @test_throws BoundsError xs[11] = 1
    @test_throws BoundsError pop!(xs)
    @test_throws BoundsError insert!(xs, 0, 1)

    for x in 1:10
        push!(xs, x)
    end

    @test_throws BoundsError push!(xs, 11)
    @test_throws BoundsError insert!(xs, 11, 1)
    @test_throws BoundsError resize!(xs, 11)
    @test_throws BoundsError splice!(xs, 0)
    @test_throws BoundsError splice!(xs, 11)
end

# Abuse low-level functions to test cases that wouldn't otherwise occur
@testset "Low Level" begin
    # findidx should return 0 when called on an empty node
    K = Int
    V = SegmentValue{Int,Int}
    B = 32
    x = Segment{Int}(1, 1)
    node = InternalNode{K,V,B}()
    @test SegmentTrees.findidx(node, x) === 0
    @test SegmentTrees.firstfrom(node, 1) == (node, 0)
    node = LeafNode{K,V,B}()
    @test SegmentTrees.findidx(node, x) == 0

    result = SegmentTrees.Intersection{K,V,B}()
    SegmentTrees.firstintersection!(node, x, nothing, result, SegmentTrees.true_cmp)
    @test result.index == 0

    @test SegmentTrees.firstfrom(node, 1) == (node, 0)

    push!(node, SegmentValue{Int,Int}(1, 1, 1))
    @test SegmentTrees.firstfrom(node, 2) == (node, 0)

    # test that delete! still works on a contrived tree with one internal and
    # one leaf node
    t = SegmentBTree{K,V,B}()
    t.root = InternalNode{K,V,B}()
    push!(t.root.children, LeafNode{K,V,B}())
    push!(t.root.maxends, 1)
    push!(t.root.keys, x)
    push!(t.root.children[1], SegmentValue{Int,Int}(1, 1, 1))
    t.root.maxend = 1
    t.root.children[1].maxend = 1
    delete!(t, (1, 1))
    @test isa(t.root, LeafNode)

    ## test that the right thing happens if you delete the last key in a non-root
    ## leaf node (which can't actually happen)
    node = LeafNode{Int,SegmentValue{Int},32}()
    push!(node, SegmentValue{Int,Int}(1, 1, 1))
    @test SegmentTrees.findidx(node, x) == 1
end


text:
## WARNING: this is an attempt to wrap cuDNNs dropout implementation
##          it doesn't work and should be considered work in progress

## Dropout descriptor

mutable struct DOD
    ptr::Ptr{Void}
    states::CuArray{UInt8,1}
end

# TODO: what should be passed for states? what size should it have?
function DOD(dropout::Float32; handle=cudnnhandle(), states=C_NULL,
             seed=UInt64(42))
    d = Ref{Cptr}(0)
    @cuda(cudnn, cudnnCreateDropoutDescriptor, (Ref{Cptr},), d)

    statesSizeInBytes = Ref{Csize_t}(0)
    @cuda(cudnn, cudnnDropoutGetStatesSize, (Cptr, Ref{Csize_t}), handle, statesSizeInBytes)

    states = CuArray{UInt8,1}(statesSizeInBytes[])
    # TODO: states will be garbage collected and thus its memory will be reclained
    # should we return it in DOD?
    @cuda(cudnn, cudnnSetDropoutDescriptor,
          (Ref{Cptr}, Cptr, Cfloat, Cptr, Cuint, Culong),
          d, handle, dropout, states.ptr, statesSizeInBytes[], seed)
    dod = DOD(d[], states)
    finalizer(dod, x->@cuda(cudnn, cudnnDestroyDropoutDescriptor, (Cptr,), x.ptr))
    # finalizer(dod, x -> cudnnDestroyDropoutDescriptor(x.ptr))
    return dod
end

# doesn't work on my machine
function get_dropout_desc(dod::DOD)
    dropout = Ref{Cfloat}(0)
    states = Ref{Cptr}(0)
    seed = Ref{Culonglong}(0)
    # TODO: cudnnGetDropoutDescriptor cannot be found. Wrong library version?
    @cuda(cudnn, cudnnGetDropoutDescriptor,
          (Cptr, Cptr, Ptr{Cfloat}, Ptr{Cptr}, Ptr{Culong}),
          d, handle, dropout, states, seed)
end

function  get_reserve_space_size(x::CuArray{T}) where T
    td = TD(x)
    sz = Ref{Csize_t}(0)
    @cuda(cudnn, cudnnDropoutGetReserveSpaceSize, (Cptr, Ref{Csize_t}), td.ptr, sz)
    return sz[]
end

function  dropout_get_states_size(handle)
    sz = Ref{Csize_t}(0)
    @cuda(cudnn, cudnnDropoutGetStatesSize, (Cptr, Ref{Csize_t}), handle, sz)
    return sz[]
end

function dropout_forward!(y::CuArray{T}, x::CuArray{T}, dropout::Float64;
                          handle=cudnnhandle(), rs=nothing) where T
    dod = DOD(Float32(dropout))
    xtd = TD(x)
    ytd = TD(y)
    if rs == nothing
        rs_sz = get_reserve_space_size(x)
        rs = CuArray(zeros(UInt8, rs_sz))
        # rs = Ptr{Void}(0)
    else
        rs_sz = length(rs)
    end
    @cuda(cudnn, cudnnDropoutForward,
          #  (Cptr, Cptr, Cptr, Cptr, Cptr, Cptr, Cptr, Csize_t),
          (Cptr, Cptr, Cptr, Cptr, Cptr, Cptr, Cptr, Csize_t),
          handle, dod.ptr, xtd.ptr, x.ptr, ytd.ptr, y.ptr, rs.ptr, rs_sz)

end

function main()
    x = CuArray(randn(Float32, 5, 4, 3, 2))
    y = similar(x)
    dropout = 0.5
    handle = cudnnhandle()
    rs = nothing

    dropoutDesc = DOD(Float32(dropout))
    xdesc = TD(x)
    ydesc = TD(y)
    x = x
    y = y

    reserveSpace = C_NULL
    reserveSizeInBytes = Cint(get_reserve_space_size(x))

    cudnnDropoutForward(handle, dropoutDesc.ptr, xdesc, x.ptr, ydesc, y.ptr, reserveSpace,
                        reserveSizeInBytes)
end


text:
module Sorting

  importall Base

  export mergeSortPartialBy!, mergeSortBy!, mergeSortPartial!, mergeSort!,
         mergeSortPartialBy, mergeSortBy, mergeSortPartial, mergeSort,
         mergeSortByFunc!, mergeSortByFunc,
         removeDuplicates!, union!,
         removeDuplicates, union

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

  function mergePartialBy!{T}( islessf, arr :: Array{T}, a :: Int, m :: Int, b :: Int )
    aa = Array{T}(b-a+1)
    i = a
    j = m+1
    k = 1
    while i <= m && j <= b
      if islessf(arr[i], arr[j])
        aa[k] = arr[i]
        i += 1
      else
        aa[k] = arr[j]
        j += 1
      end
      k += 1
    end
    while i <= m
      aa[k] = arr[i]
      i += 1
      k += 1
    end
    while j <= b
      aa[k] = arr[j]
      j += 1
      k += 1
    end
    for x = a:b
      arr[x] = aa[x-a+1]
    end
    return arr
  end

  function mergeSortPartialBy!{T}( islessf, arr :: Array{T}, a :: Int, b :: Int )
    if a == b
      return arr
    end
    if b-a == 1
      if islessf(arr[b], arr[a])
        c = arr[a]
        arr[a] = arr[b]
        arr[b] = c
      end
      return arr
    end
    m = div(b-a+1,2) + a - 1
    mergeSortPartialBy!(islessf, arr, a, m)
    mergeSortPartialBy!(islessf, arr, m+1, b)
    return mergePartialBy!(islessf, arr, a, m, b)
  end

  function mergeSortBy!{T}( islessf, arr :: Array{T} )
    if length(arr) == 1
      return arr
    end
    return mergeSortPartialBy!( islessf, arr, Int(1), length(arr) )
  end

  function mergeSortByFunc!{T}( f, arr :: Array{T} )
    return mergeSortBy!( (x,y) -> f(x) < f(y), arr )
  end

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

  # assumes the array is sorted already
  function removeDuplicates!{T}( arr :: Array{T} )
    n = length(arr)
    if n==0
      return arr
    end
    j = 0
    for i = 1:(n-1)
      if arr[i] != arr[i+1]
        j += 1
        arr[j] = arr[i]
      end
    end
    j += 1
    arr[j] = arr[n]
    # j is number of unique values
    resize!(arr,j)
    return arr
  end

  function union!{T}( arr :: Array{T} )
    return removeDuplicates!( mergeSort!( arr ) )
  end

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

  # assumes the array is sorted already
  function removeDuplicates{T}( arr :: Array{T} )
    return removeDuplicates!( copy(arr) )
  end

  function union{T}( arr :: Array{T} )
    return union!( copy(arr) )
  end

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

  function mergePartial!{T}( arr :: Array{T}, a :: Int, m :: Int, b :: Int )
    return mergePartialBy!( isless, arr, a, m, b )
  end

  function mergeSortPartial!{T}( arr :: Array{T}, a :: Int, b :: Int )
    return mergeSortPartialBy!( isless, arr, a, b )
  end

  function mergeSort!{T}( arr :: Array{T} )
    return mergeSortBy!( isless, arr )
  end

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

  function mergePartialBy{T}( islessf, arr :: Array{T}, a :: Int, m :: Int, b :: Int )
    return mergePartialBy!( islessf, copy(arr), a, m, b )
  end

  function mergeSortPartialBy{T}( islessf, arr :: Array{T}, a :: Int, b :: Int )
    return mergeSortPartialBy!( islessf, copy(arr), a, b )
  end

  function mergeSortBy{T}( islessf, arr :: Array{T} )
    return mergeSortBy!( islessf, copy(arr) )
  end

  function mergeSortByFunc{T}( f, arr :: Array{T} )
    return mergeSortByFunc!( f, copy(arr) )
  end

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

  function mergePartial{T}( arr :: Array{T}, a :: Int, m :: Int, b :: Int )
    return mergePartial!( copy(arr), a, m, b )
  end
  function mergeSortPartial{T}( arr :: Array{T}, a :: Int, b :: Int )
    return mergeSortPartial!( copy(arr), a, b )
  end
  function mergeSort{T}( arr :: Array{T} )
    return mergeSort!( copy(arr) )
  end

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

end

max_stars_repo_name: clauderichard/RubikSearch

max_issues_repo_name: clauderichard/RubikSearch

max_forks_repo_name: clauderichard/RubikSearch


text:
# line styles - Plots
styles = filter((s->begin
                s in Plots.supported_styles()
            end), [:solid, :dash, :dot, :dashdot, :dashdotdot])
styles = reshape(styles, 1, length(styles))
n = length(styles)
y = cumsum(randn(20, n), dims = 1)
Plots.plot(y, line = (5, styles), label = map(string, styles), legendtitle = "linestyle")

# line styles - makie
using CairoMakie

f = Figure()
Axis(f[1, 1])

xs = 0:0.01:10
ys = 0.5 .* sin.(xs)

for (i, lw) in enumerate([1, 2, 3])
    lines!(xs, ys .- i/6, linestyle = nothing, linewidth = lw)
    lines!(xs, ys .- i/6 .- 1, linestyle = :dash, linewidth = lw)
    lines!(xs, ys .- i/6 .- 2, linestyle = :dot, linewidth = lw)
    lines!(xs, ys .- i/6 .- 3, linestyle = :dashdot, linewidth = lw)
    lines!(xs, ys .- i/6 .- 4, linestyle = :dashdotdot, linewidth = lw)
    lines!(xs, ys .- i/6 .- 5, linestyle = [0.5, 1.0, 1.5, 2.5], linewidth = lw)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
divisors(number) = [i for i in 1:(number - 1) if iszero(number % i)]

function validint(number)
    if number < 1
        throw(DomainError("Number needs to be a positive integer"))
    end
    return true
end

isperfect(number) = validint(number) && sum(divisors(number)) == number

isabundant(number) = validint(number) && sum(divisors(number)) > number

isdeficient(number) = validint(number) && sum(divisors(number)) < number

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using CSV, DataFrames, Plots, StatsPlots, Missings, CategoricalArrays
using Flux
using MLDataUtils
using PyCall

data = CSV.read("julia/cleaned.csv")

### Analyze the participants by career stage
career_group = groupby(data, :Position);
career_breakdown = sort(combine(career_group, nrow), :x1);
career_breakdown[!, :Index] = [1, 3, 4, 2];
sort!(career_breakdown, :Index);

career_breakdown.Position

# plot the data
gr(size=(600, 500))
chart = bar(
    career_breakdown[!, :x1],
    color = "#b85231",
    leg = false,
    lw = 1.0,
    bar_width = 0.3,
    background_color = "#f9efde",
    titlefont = font(14, "Open Sans"),
)
title!("Participant breakdown")
ylabel!("Counts", font = font(12, "Open Sans"))
xlabel!("Career stage", font = font(12, "Open Sans"))
xticks!([1, 2, 3, 4], career_breakdown[!, :Position])
savefig("julia/career_plot.png")

# Generate a corner plot for everything
mat_rep = Matrix(data[[:Lat, :Long, :HotDogEncoding, :PosEncoding]]);

gr(size = (600, 500))
corrplot(
    mat_rep,
    label = ["Latitude", "Longitude", "Hot Dog?", "Career Step"],
    grid = false,
    compact = true,
)

# Get raw stats out; how likely each career stage is to be affirmative
true_stats = combine(career_group, :HotDogEncoding => sum);
true_stats[!,:nentries] = combine(career_group, nrow)[:x1]
true_stats[!,:percent] = true_stats[:HotDogEncoding_sum] ./ true_stats[:nentries];
true_stats[!,:Index] = [2,4,3,1];
sort!(true_stats, :Index)

gr(size=(600, 500))
chart = bar(
    true_stats[!, :percent],
    color = "#b85231",
    leg = false,
    lw = 1.0,
    bar_width = 0.3,
    background_color = "#f9efde",
    titlefont = font(14, "Open Sans"),
)
hline!([1.], lw=1.5)
title!("Sandwich vs. Career stage")
ylabel!("Fraction of Yes", font = font(12, "Open Sans"))
xlabel!("Career stage", font = font(12, "Open Sans"))
xticks!([1, 2, 3, 4], career_breakdown[!, :Position])
annotate!(1.3, 0.95, "Correct answer", font=font(1, "Open Sans"))
savefig("julia/career_sandwich.png")

for col in [:Position, :Hotdog]
    categorical!(data, col);
end

gr(size=(600, 500))
@df data violin(:Position, :HotDogEncoding, alpha=0.7, background_color = "#f9efde", color="#b85231", marker=(0.2,"#b85231",stroke(0)), leg=false)
@df data dotplot!(:Position, :HotDogEncoding, marker=(:black,stroke(0)), leg=false, alpha=0.4)
yticks!([0., 1.], ["No", "Yes"])
ylabel!("Response")
savefig("julia/career_sandwich.png")

yaynay = combine(groupby(data, :HotDogEncoding), nrow);

gr(size=(600, 500))
chart = bar(
    yaynay[!, :x1],
    color = "#b85231",
    leg = false,
    lw = 1.0,
    bar_width = 0.3,
    background_color = "#f9efde",
    titlefont = font(14, "Open Sans"),
)
title!("Are hotdogs sandwiches?")
ylabel!("Counts", font = font(12, "Open Sans"))
xlabel!("Response", font = font(12, "Open Sans"))
xticks!([1, 2], ["Yes", "No"])
annotate!(1., 70., "Correct answer", font=font(1, "Open Sans"))
savefig("julia/hotdog_responses.png")

US = filter(row -> row.Long < -20., data)
EU = filter(row -> row.Long > -20., data)

US_yay = sum(US.HotDogEncoding) / nrow(US)
EU_yay = sum(EU.HotDogEncoding) / nrow(EU)

df = DataFrame([[US_yay, EU_yay]])

gr(size=(600, 500))
chart = bar(
    df[:,:x1],
    color = "#b85231",
    leg = false,
    lw = 1.0,
    bar_width = 0.3,
    background_color = "#f9efde",
    titlefont = font(14, "Open Sans"),
)
# title!("Europeans more likely to be right")
ylabel!("Fraction of Yes", font = font(12, "Open Sans"))
xlabel!("Side of the World", font = font(12, "Open Sans"))
xticks!([1, 2], ["Negative longitude", "Positive longitude"])
annotate!(1.3, 0.95, "Correct answer", font=font(1, "Open Sans"))
hline!([1.], lw=1.5)
savefig("julia/split_world.png")

length(unique(data.Org))


text:
abstract type AccelerationStyle end

struct UnknownStyle <: AccelerationStyle end

struct NoAccelerationStyle <: AccelerationStyle end

struct QuasiNewtonStyle <: AccelerationStyle end

struct NesterovStyle <: AccelerationStyle end

acceleration_style(::Type{<:Any}) = UnknownStyle()


text:
using LightGraphs
using GraphPlot

function visible_graph(x)
    N = length(x)
    g = SimpleGraph(N)
    for i = 1:N
        for j = i+1:N
            cnt = 0
            for k = i+1:j-1
                if x[k]>=x[j]+(x[i]-x[j])*(j-k)/(j-i)
                    cnt = 1
                    break
                end
            end
            if cnt == 0
                add_edge!(g,i,j)
            end
        end
    end
    return g
end

max_stars_repo_name: ryusei-nyaan/visiblegraph

max_issues_repo_name: ryusei-nyaan/visiblegraph

max_forks_repo_name: ryusei-nyaan/visiblegraph


text:
@testset "Country.jl" begin
    filename = "data\\test.country"
    dc = DDR2import.Country.read(filename)
    @test dc["AG"].name == "Solomon Islands"
    @test dc["EB"].member
    @test !dc["ZY"].member
end


text:
using MixedModels
using Test

using MixedModels: dataset, likelihoodratiotest
using GLM: ProbitLink

@testset "likelihoodratio test" begin
    slp = dataset(:sleepstudy);
    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp);
    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp);
    lrt = likelihoodratiotest(fm0,fm1);

    @test [deviance(fm0), deviance(fm1)] == lrt.deviance
    @test deviance(fm0) - deviance(fm1) == only(lrt.tests.deviancediff)
    @test only(lrt.tests.dofdiff) == 1
    @test sum(map(length,lrt.tests)) == 3
    @test sum(map(length,lrt.pvalues)) == 1
    @test sum(map(length,lrt.models)) == 4
    @test length(lrt.formulae) == 2
    show(IOBuffer(),lrt);
    @test :pvalues in propertynames(lrt)

    # mix of REML and ML
    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp, REML=true);
    @test_throws ArgumentError likelihoodratiotest(fm0,fm1)

    # differing FE with REML
    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp, REML=true);

    @test_throws ArgumentError likelihoodratiotest(fm0,fm1)

    contra = MixedModels.dataset(:contra);
    gm0 = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urbdist)), contra, Bernoulli(), fast=true);
    gm1 = fit(MixedModel, @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urbdist)), contra, Bernoulli(), fast=true);
    lrt = likelihoodratiotest(gm0,gm1);
    @test [deviance(gm0), deviance(gm1)] == lrt.deviance
    @test deviance(gm0) - deviance(gm1) == only(lrt.tests.deviancediff)
    @test first(lrt.tests.dofdiff) == 1
    @test sum(length, lrt.tests) == 3
    @test sum(length, lrt.pvalues) == 1
    @test sum(length, lrt.models) == 4
    @test length(lrt.formulae) == 2

    # mismatched links
    gm_probit = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urbdist)), contra, Bernoulli(), ProbitLink(), fast=true);
    @test_throws ArgumentError likelihoodratiotest(gm0,gm_probit)

    # mismatched families
    gm_poisson = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urbdist)), contra, Poisson(), fast=true);
    @test_throws ArgumentError MixedModels.likelihoodratiotest(gm0,gm_poisson)
end


text:
# ENV["JULIA_DEBUG"] = "all"

include("debugger.jl")

let
    local pipenames = VSCodeDebugger.clean_up_ARGS_in_launch_mode()
    try
        VSCodeDebugger.startdebug(pipenames[1])
    catch err
        VSCodeDebugger.global_err_handler(err, catch_backtrace(), pipenames[2])
    end
end

max_stars_repo_name: cloudgiant/julia-vscode

max_issues_repo_name: cloudgiant/julia-vscode

max_forks_repo_name: cloudgiant/julia-vscode


text:
# Special functions
import Base.Math.@horner

# See:
#   Martin Maechler (2012) "Accurately Computing log(1 − exp(− |a|))"
#   http://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf

# log(1-exp(x))
# NOTE: different than Maechler (2012), no negation inside parantheses
log1mexp(x::Real) = x >= loghalf ? log(-expm1(x)) : log1p(-exp(x))
# log(1+exp(x))
log1pexp(x::Real) = x > zero(x) ? x + log1p(exp(-x)) : log1p(exp(x))
log1pexp(x::Float64) = x <= 18.0 ? log1p(exp(x)) : x <= 33.3 ? x + exp(-x) : x
log1pexp(x::Float32) = x <= 9f0 ? log1p(exp(x)) : x <= 16f0 ? x + exp(-x) : x
log1pexp(x::Integer) = log1pexp(float(x))

log2mexp(x::Real) = log1p(-expm1(x))
# log(exp(x)-1)
logexpm1(x::Real) = x + log1mexp(-x)
logexpm1(x::Float64) = x <= 18.0 ? log(expm1(x)) : x <= 33.3 ? x - exp(-x) : x
logexpm1(x::Float32) = x <= 9f0 ? log(expm1(x)) : x <= 16f0 ? x - exp(-x) : x
logexpm1(x::Integer) = logexpm1(float(x))

# log(1+x^2)
log1psq(x::FloatingPoint) = (ax = abs(x); ax < maxintfloat(x) ? log1p(ax*ax) : 2*log(ax))

# log(1+x)-x
# accurate ~2ulps for -0.227 < x < 0.315
function log1pmx_kernel(x::Float64)
    r = x/(x+2.0)
    t = r*r
    w = @horner(t,
                6.66666666666666667e-1, # 2/3
                4.00000000000000000e-1, # 2/5
                2.85714285714285714e-1, # 2/7
                2.22222222222222222e-1, # 2/9
                1.81818181818181818e-1, # 2/11
                1.53846153846153846e-1, # 2/13
                1.33333333333333333e-1, # 2/15
                1.17647058823529412e-1) # 2/17
    hxsq = 0.5*x*x
    r*(hxsq+w*t)-hxsq
end

# use naive calculation or range reduction outside kernel range.
# accurate ~2ulps for all x
function log1pmx(x::Float64)
    if !(-0.7 < x < 0.9)
        return log1p(x) - x
    elseif x > 0.315
        u = (x-0.5)/1.5
        return log1pmx_kernel(u) - 9.45348918918356180e-2 - 0.5*u
    elseif x > -0.227
        return log1pmx_kernel(x)
    elseif x > -0.4
        u = (x+0.25)/0.75
        return log1pmx_kernel(u) - 3.76820724517809274e-2 + 0.25*u
    elseif x > -0.6
        u = (x+0.5)*2.0
        return log1pmx_kernel(u) - 1.93147180559945309e-1 + 0.5*u
    else
        u = (x+0.625)/0.375
        return log1pmx_kernel(u) - 3.55829253011726237e-1 + 0.625*u
    end
end

# log(x) - x + 1
function logmxp1(x::Float64)
    if x <= 0.3
        return (log(x) + 1.0) - x
    elseif x <= 0.4
        u = (x-0.375)/0.375
        return log1pmx_kernel(u) - 3.55829253011726237e-1 + 0.625*u
    elseif x <= 0.6
        u = 2.0*(x-0.5)
        return log1pmx_kernel(u) - 1.93147180559945309e-1 + 0.5*u
    else
        return log1pmx(x-1.0)
    end
end

φ(z::Real) = exp(-0.5*z*z)/√2π
logφ(z::Real) = -0.5*(z*z + log2π)

Φ(z::Real) = 0.5*erfc(-z/√2)
Φc(z::Real) = 0.5*erfc(z/√2)
logΦ(z::Real) = z < -1.0 ? log(0.5*erfcx(-z/√2)) - 0.5*z*z : log1p(-0.5*erfc(z/√2))
logΦc(z::Real) = z > 1.0 ? log(0.5*erfcx(z/√2)) - 0.5*z*z : log1p(-0.5*erfc(-z/√2))

# Rational approximations for the inverse cdf, from:
#   Wichura, M.J. (1988) Algorithm AS 241: The Percentage Points of the Normal Distribution
#   Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 37, No. 3, pp. 477-484
for (fn,arg) in ((:Φinv,:p),(:logΦinv,:logp))
    @eval begin
        function $fn($arg::Float32)
            if $(fn == :Φinv)
                q = p - 0.5f0
            else
                q = exp(logp) - 0.5f0
            end
            if abs(q) <= 0.425f0
                r = 0.180625f0 - q*q
                return q * @horner(r,
                                   3.38713_27179f0,
                                   5.04342_71938f1,
                                   1.59291_13202f2,
                                   5.91093_74720f1,
                                   ) /
                @horner(r,
                        1.0f0,
                        1.78951_69469f1,
                        7.87577_57664f1,
                        6.71875_63600f1)
            else
                if $(fn == :Φinv)
                    if p <= 0f0
                        return p == 0f0 ? -inf(Float32) : nan(Float32)
                    elseif p >= 1f0
                        return p == 1f0 ? inf(Float32) : nan(Float32)
                    end
                    r = sqrt(q < 0f0 ? -log(p) : -log1p(-p))
                else
                    if logp == -Inf
                        return -inf(Float32)
                    elseif logp >= 0f0
                        return logp == 0f0 ? inf(Float32) : nan(Float32)
                    end
                    r = sqrt(qf0 < 0 ? -logp : -log1mexp(logp))
                end
                if r < 5.0f0
                    r -= 1.6f0
                    z = @horner(r,
                                1.42343_72777f0,
                                2.75681_53900f0,
                                1.30672_84816f0,
                                1.70238_21103f-1) /
                    @horner(r,
                            1.0f0,
                            7.37001_64250f-1,
                            1.20211_32975f-1)
                else
                    r -= 5.0f0
                    z = @horner(r,
                                6.65790_51150f0,
                                3.08122_63860f0,
                                4.28682_94337f-1,
                                1.73372_03997f-2) /
                    @horner(r,
                            1.0f0,
                            2.41978_94225f-1,
                            1.22582_02635f-2)
                end
                return copysign(z,q)
            end
        end

        function $fn($arg::Real)
            if $(fn == :Φinv)
                q = p - 0.5
            else
                q = exp(logp) - 0.5
            end
            if abs(q) <= 0.425
                r = 0.180625 - q*q
                return q * @horner(r,
                                   3.38713_28727_96366_6080e0,
                                   1.33141_66789_17843_7745e2,
                                   1.97159_09503_06551_4427e3,
                                   1.37316_93765_50946_1125e4,
                                   4.59219_53931_54987_1457e4,
                                   6.72657_70927_00870_0853e4,
                                   3.34305_75583_58812_8105e4,
                                   2.50908_09287_30122_6727e3) /
                @horner(r,
                        1.0,
                        4.23133_30701_60091_1252e1,
                        6.87187_00749_20579_0830e2,
                        5.39419_60214_24751_1077e3,
                        2.12137_94301_58659_5867e4,
                        3.93078_95800_09271_0610e4,
                        2.87290_85735_72194_2674e4,
                        5.22649_52788_52854_5610e3)
            else
                if $(fn == :Φinv)
                    if p <= 0.0
                        return p == 0.0 ? -inf(Float64) : nan(Float64)
                    elseif p >= 1.0
                        return p == 1.0 ? inf(Float64) : nan(Float64)
                    end
                    r = sqrt(q < 0 ? -log(p) : -log1p(-p))
                else
                    if logp == -Inf
                        return -inf(Float64)
                    elseif logp >= 0.0
                        return logp == 0.0 ? inf(Float64) : nan(Float64)
                    end
                    r = sqrt(q < 0 ? -logp : -log1mexp(logp))
                end
                if r < 5.0
                    r -= 1.6
                    z = @horner(r,
                                1.42343_71107_49683_57734e0,
                                4.63033_78461_56545_29590e0,
                                5.76949_72214_60691_40550e0,
                                3.64784_83247_63204_60504e0,
                                1.27045_82524_52368_38258e0,
                                2.41780_72517_74506_11770e-1,
                                2.27238_44989_26918_45833e-2,
                                7.74545_01427_83414_07640e-4) /
                    @horner(r,
                            1.0,
                            2.05319_16266_37758_82187e0,
                            1.67638_48301_83803_84940e0,
                            6.89767_33498_51000_04550e-1,
                            1.48103_97642_74800_74590e-1,
                            1.51986_66563_61645_71966e-2,
                            5.47593_80849_95344_94600e-4,
                            1.05075_00716_44416_84324e-9)
                else
                    r -= 5.0
                    z = @horner(r,
                                6.65790_46435_01103_77720e0,
                                5.46378_49111_64114_36990e0,
                                1.78482_65399_17291_33580e0,
                                2.96560_57182_85048_91230e-1,
                                2.65321_89526_57612_30930e-2,
                                1.24266_09473_88078_43860e-3,
                                2.71155_55687_43487_57815e-5,
                                2.01033_43992_92288_13265e-7) /
                    @horner(r,
                            1.0,
                            5.99832_20655_58879_37690e-1,
                            1.36929_88092_27358_05310e-1,
                            1.48753_61290_85061_48525e-2,
                            7.86869_13114_56132_59100e-4,
                            1.84631_83175_10054_68180e-5,
                            1.42151_17583_16445_88870e-7,
                            2.04426_31033_89939_78564e-15)
                end
                return copysign(z,q)
            end
        end
    end
end

# Multidimensional gamma / partial gamma function
function lpgamma(p::Int, a::Float64)
    res::Float64 = p * (p - 1.0) / 4.0 * log(pi)
    for ii in 1:p
        res += lgamma(a + (1.0 - ii) / 2.0)
    end
    return res
end

# Remainder term after Stirling's approximation to the log-gamma function
# lstirling(x) = lgamma(x) + x - (x-0.5)*log(x) - 0.5*log2π
#              = 1/(12x) - 1/(360x^3) + 1/(1260x^5) + ...
# Asymptotic expansion from:
#   Temme, N. (1996) Special functions: An introduction to the classical
#   functions of mathematical physics, Wiley, New York, ISBN: 0-471-11313-1,
#   Chapter 3.6, pp 61-65.
# Relative error of approximation is bounded by
#   (174611/125400 x^-19) / (1/12 x^-1 - 1/360 x^-3)
# which is < 1/2 ulp for x >= 10.0
# total numeric error appears to be < 2 ulps
#
lstirling_asym(x::Integer) = lstirling_asym(float(x))

function lstirling_asym(x::Float64)
    t = 1.0/(x*x)
    @horner(t,
             8.33333333333333333e-2, #  1/12 x^-1
            -2.77777777777777778e-3, # -1/360 x^-3
             7.93650793650793651e-4, #  1/1260 x^-5
            -5.95238095238095238e-4, # -1/1680 x^-7
             8.41750841750841751e-4, #  1/1188 x^-9
            -1.91752691752691753e-3, # -691/360360 x^-11
             6.41025641025641026e-3, #  1/156 x^-13
            -2.95506535947712418e-2, # -3617/122400 x^-15
             1.79644372368830573e-1)/x #  43867/244188 x^-17
end

function lstirling_asym(x::Float32)
    t = 1f0/(x*x)
    @horner(t,
             8.333333333333f-2, #  1/12 x^-1
            -2.777777777777f-3, # -1/360 x^-3
             7.936507936508f-4, #  1/1260 x^-5
            -5.952380952381f-4, # -1/1680 x^-7
             8.417508417508f-4)/x #  1/1188 x^-9
end