text:
using Test, Random, LSHFunctions

@testset "Test intervals with LSHFunctions.@interval" begin
    Random.seed!(0)

    @testset "Construct an interval" begin
        # Create a new interval corresponding to [0,1)
        interval = LSHFunctions.@interval(0 ≤ x < 1)
        @test isa(interval, LSHFunctions.RealInterval)
        @test interval.lower == 0
        @test interval.upper == 1
        @test interval.closed_below == true
        @test interval.closed_above == false

        # Create another interval corresponding to (-Inf,Inf]
        interval = LSHFunctions.@interval(-Inf < x ≤ Inf)
        @test isa(interval, LSHFunctions.RealInterval)
        @test interval.lower == -Inf
        @test interval.upper == Inf
        @test interval.closed_below == false
        @test interval.closed_above == true

        # Construct an interval containing a single point
        interval = LSHFunctions.@interval(0 ≤ x ≤ 0)
        @test interval.lower == interval.upper == 0
        @test interval.closed_below == interval.closed_above == true

        # Construct an interval in reverse order
        @test LSHFunctions.@interval(1 > x ≥ 0) == LSHFunctions.@interval(0 ≤ x < 1)
        @test LSHFunctions.@interval(Inf > x ≥ -Inf) == LSHFunctions.@interval(-Inf ≤ x < Inf)

        # Construct interval using the <= symbols instead of ≤
        @test LSHFunctions.@interval(0 <= x <= 1) == LSHFunctions.@interval(0 ≤ x ≤ 1)
        @test LSHFunctions.@interval(1 >= x >= 0) == LSHFunctions.@interval(1 ≥ x ≥ 0)

        # Construct an interval using variables and expressions for the endpoints
        a, b = -1, 1
        @test LSHFunctions.@interval(a ≤ x ≤ b) == LSHFunctions.@interval(-1 ≤ x ≤ 1)
        @test LSHFunctions.@interval(2a ≤ x ≤ 2b) == LSHFunctions.@interval(-2 ≤ x ≤ 2)

        # Construct an interval that contains a different symbol in the middle
        @test LSHFunctions.@interval(0 ≤ y < 1) == LSHFunctions.@interval(0 ≤ x < 1)
        @test LSHFunctions.@interval(0 ≤ mynum < 1) == LSHFunctions.@interval(0 ≤ x < 1)
    end

    @testset "Incorrect interval construction raises errors" begin
        # Try to construct intervals using invalid expressions
        @test_throws(ErrorException, LSHFunctions.@interval(0))
        @test_throws(ErrorException, LSHFunctions.@interval(0 ≤ x))
        @test_throws(ErrorException, LSHFunctions.@interval(0 ≤ 1 ≤ 2))

        # Try to construct intervals using inequalities pointing in different directions
        @test_throws(ErrorException, LSHFunctions.@interval(0 ≤ x ≥ 1))
        @test_throws(ErrorException, LSHFunctions.@interval(0 > x < 1))

        # Try to construct intervals not on the real line
        @test_throws(MethodError, LSHFunctions.@interval(-1im ≤ x ≤ 1))
        @test_throws(MethodError, LSHFunctions.@interval(0 ≤ x ≤ 1im))
        @test_throws(MethodError, LSHFunctions.@interval(-1im ≤ x ≤ 1im))

        # Try to construct intervals using undefined endpoints
        @test_throws(UndefVarError, LSHFunctions.@interval(a ≤ x ≤ Inf))
        @test_throws(UndefVarError, LSHFunctions.@interval(-Inf ≤ x ≤ b))
    end

    @testset "Test interval methods" begin
        # isempty
        @test isempty(LSHFunctions.@interval(0 ≤ x < -1))
        @test isempty(LSHFunctions.@interval(0 ≤ x < 0))

        # width
        @test LSHFunctions.width(LSHFunctions.@interval(0 ≤ x < 1)) == 1
        @test LSHFunctions.width(LSHFunctions.@interval(0 ≤ x < Inf)) == Inf
        @test LSHFunctions.width(LSHFunctions.@interval(0 ≤ x < -1)) == 0
    end

    @testset "Test membership within an interval" begin
        # Test 1: interval [0,1)
        interval = LSHFunctions.@interval(0 ≤ x < 1)

        @test all(x ∈ interval for x in rand(32))
        @test 1 ∉ interval

        # Test 2: interval (-Inf, Inf]
        interval = LSHFunctions.@interval(-Inf < x ≤ Inf)

        @test all(log(x) ∈ interval for x in rand(32))
        @test Inf ∈ interval
        @test -Inf ∉ interval
    end

    @testset "Test interval intersection" begin
        @test LSHFunctions.@interval(0 < x ≤ 1) ∩ LSHFunctions.@interval(0 ≤ x ≤ 1) ==
              LSHFunctions.@interval(0 < x ≤ 1)
        @test LSHFunctions.@interval(0 < x ≤ 1) ∩ LSHFunctions.@interval(1/2 ≤ x < 1) ==
              LSHFunctions.@interval(1/2 ≤ x < 1)

        # Ensure that the intersect() function is equal to ∩
        @test let success = true
            for ii = 1:128
                a, c = rand(), rand()
                b, d = a + rand(), c + rand()
                interval_1 = LSHFunctions.@interval(a < x ≤ b)
                interval_2 = LSHFunctions.@interval(c ≤ x < d)
                success &= (LSHFunctions.intersect(interval_1, interval_2) == interval_1 ∩ interval_2)
            end
            success
        end

        # Intersect two disjoint intervals to create an empty interval
        @test isempty(LSHFunctions.@interval(0 < x < 1) ∩ LSHFunctions.@interval(1 < x < 2))
    end

    @testset "Test empty intervals" begin
        # Test for equality between empty intervals
        @test LSHFunctions.@interval(0 < x < 0) == LSHFunctions.@interval(1 < x < 1)
        @test LSHFunctions.@interval(0 < x < -1) == LSHFunctions.@interval(0 < x < 0)
        @test LSHFunctions.@interval(0 < x < 0) == LSHFunctions.@interval(0 < x < -1)
    end
end


text:
cd(@__DIR__); include("setups/rgc100.jl");
pyplot(dpi = 200)

##
ϵ = 0.3
for J = 1:2
    Uf = unitary_folding_operator(W, GP; ϵ = ϵ, J = J)
    gplot(W, X; width = 1)
    scatter_gplot!(X; marker = diag(Uf), plotOrder = :l2s, ms = 4)
    plt = plot!(xlims = [-180, 220], xlabel = "x(μm)", ylabel = "y(μm)",
          clims = (0.6, 1), frame = :box, size = (520, 500))
    savefig(plt, "../figs/RGC100_unitary_folding_diag_J$(J).png")
end

ext: jl


text:
export TotalDegreeSolutionIterator, totaldegree_startsystem, totaldegree

"""
    TotalDegreeSolutionIterator(degrees, b)

Given the `Vector`s `degrees` and `b` `TotalDegreeSolutionIterator` enumerates all solutions
of the system
```math
\\begin{align*}
    z_1^{d_1} - b_1 &= 0 \\\\
    z_1^{d_2} - b_2 &= 0 \\\\
    &\\vdots \\\\
    z_n^{d_n} - b_n &= 0 \\\\
\\end{align*}
```
where ``d_i`` is `degrees[i]` and ``b_i`` is `b[i]`.
"""
struct TotalDegreeSolutionIterator{T<:Complex, Iter}
    degrees::Vector{Int}
    b::Vector{T}

    iterator::Iter
end
function TotalDegreeSolutionIterator(degrees::Vector{Int}, b::Vector{T}) where {T<:Complex}
    iterator = Base.Iterators.product(map(d -> 0:d-1, degrees)...)
    TotalDegreeSolutionIterator(degrees, b, iterator)
end

Base.start(iter::TotalDegreeSolutionIterator) = start(iter.iterator)
function Base.next(iter::TotalDegreeSolutionIterator{T}, state) where {T<:Complex}
    indices, nextstate = next(iter.iterator, state)

    value = map(1:length(indices), indices) do i, k
        d = iter.degrees[i]
        convert(T, iter.b[i]^(1/d) * exp(2π*im/d)^k)
    end
    value, nextstate
end
Base.done(iter::TotalDegreeSolutionIterator, state) = done(iter.iterator, state)
Base.length(iter::TotalDegreeSolutionIterator) = length(iter.iterator)
Base.eltype(iter::TotalDegreeSolutionIterator{T}) where {T} = Vector{T}

"""
    totaldegree(H::Type{AbstractPolynomialHomotopy}, F, [unitroots=false])

Construct a  total degree homotopy of type `H` with `F` and an iterator of its solutions.
This is the homotopy with start system
```math
\\begin{align*}
    z_1^{d_1} &- b_1\\\\
    z_1^{d_2} &- b_2\\\\
    &\\vdots \\\\
    z_n^{d_n} &- b_n\\\\
\\end{align*}
```
and target system `F`, where ``d_i`` is the degree of the ``i``-th polynomial of `F`.
If `unitroots=true` then ``b_i=1`` otherwise ``b_i`` is a random
complex number (with real and imaginary part in the unit interval).

## Example
```julia
H, startsolutions = totaldegree(StraightLineHomotopy{Complex128}, [x^2+y+1, x^3*y-2])
```
"""
function totaldegree(H::Type{<:AbstractPolynomialHomotopy{T}}, fs::Vector{<:MP.AbstractPolynomial{U}}; kwargs...) where {T<:Complex, U<:Number}
    F = convert(Vector{FP.Polynomial{T}}, fs)
    totaldegree(H, F; kwargs...)
end
function totaldegree(H::Type{<:AbstractPolynomialHomotopy}, fs::Vector{<:MP.AbstractPolynomial{U}}; kwargs...) where {U<:Number}
    F = convert(Vector{FP.Polynomial{promote_type(Complex128, U)}}, fs)
    totaldegree(H, F; kwargs...)
end

function totaldegree(H::Type{<:AbstractPolynomialHomotopy}, F::Vector{FP.Polynomial{U}}; kwargs...) where {U<:Number}
    G, solutions = totaldegree_startsystem(F; kwargs...)
    H(G, F), solutions
end

"""
    totaldegree_startsystem(F::Vector{FP.Polynomial{<:Complex}}, [unit_roots=false])

Return the system
```math
\\begin{align*}
    z_1^{d_1} &- b_1\\\\
    z_1^{d_2} &- b_2\\\\
    &\\vdots \\\\
    z_n^{d_n} &- b_n\\\\
\\end{align*}
```
where ``d_i`` is the degree of the ``i``-th polynomial of `F` and an iterator of its
solutions.
If `unitroots=true` then ``b_i=1`` otherwise ``b_i`` is a random
complex number (with real and imaginary part in the unit interval).
"""
function totaldegree_startsystem(F::Vector{FP.Polynomial{Complex{T}}}; unitroots=false) where {T<:AbstractFloat}
    vars = FP.variables(F[1])
    degrees = FP.degree.(F)
    if unitroots
        b = ones(Complex{T}, length(degrees))
    else
        b = convert.(Complex{T}, rand(Complex128, length(degrees)))
    end

    n = length(degrees)
    N = length(vars)
    make_homogenous = !all(FP.ishomogenous, F) && length(vars) == length(F) + 1
    if n + 1 == N
        if !all(FP.ishomogenous, F)
            return error("In order to create a total degree start system your input system needs " *
                "to have the same number of variables as number of equations or to have N variables
                and N - 1 homogenous polynomials. Currently your system has " *
                "$(length(vars)) variables and $(length(degrees)) equations and is not homogenous.")
        end

        G = map(2:N, degrees) do i, d
            exponents = zeros(Int, N, 2)
            exponents[i, 1] = d
            exponents[1, 2] = d
            coeffs = zeros(Complex{T}, 2)
            coeffs[1] = one(Complex{T})
            coeffs[2] = -b[i - 1]
            FP.Polynomial(exponents, coeffs, vars)
        end

        startsolutions::Vector{Vector{Complex{T}}} =
            map(TotalDegreeSolutionIterator(degrees, b)) do sol
                unshift!(sol, 1)
            end
        G, startsolutions
    elseif n != N
        return error("In order to create a total degree start system your input system needs " *
            "to have the same number of variables as number of equations. Currently your system has " *
            "$(length(vars)) variables and $(length(degrees)) equations.")
    else
        G = map(1:n, degrees) do i, d
            exponents = zeros(Int, n, 2)
            exponents[i, 1] = d
            coeffs = zeros(Complex{T}, 2)
            coeffs[1] = one(Complex{T})
            coeffs[2] = -b[i]
            FP.Polynomial(exponents, coeffs, vars)
        end
        G, collect(TotalDegreeSolutionIterator(degrees, b))::Vector{Vector{Complex{T}}}
    end
end
function totaldegree_startsystem(F::Vector{FP.Polynomial{T}}; kwargs...) where {T<:Number}
    totaldegree_startsystem(convert.(FP.Polynomial{complex(float(T))}, F); kwargs...)
end


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

"""
Returns true if any variable in the collection "vars" is used in any statement whose top level number is in "top_level_numbers".
We use expr_state "state" to get the block liveness information from which we use "def" and "use" to determine if a variable usage is present.
"""
function isSymbolsUsed(vars :: Dict{LHSVar,LHSVar}, top_level_numbers :: Array{Int,1}, state)
    @dprintln(3,"isSymbolsUsed: vars = ", vars, " typeof(vars) = ", typeof(vars), " top_level_numbers = ", top_level_numbers)
    bl = state.block_lives

    for i in top_level_numbers
        tls = CompilerTools.LivenessAnalysis.find_top_number(i, bl)
        assert(tls != nothing)

        for v in vars
            if in(v[1], tls.def) || in(v[2], tls.def)
                @dprintln(3, "isSymbolsUsed: ", v, " defined in statement ", i)
                return true
            elseif in(v[1], tls.use) || in(v[2], tls.use)
                @dprintln(3, "isSymbolsUsed: ", v, " used in statement ", i)
                return true
            end
        end
    end

    @dprintln(3, "isSymbolsUsed: ", vars, " not used in statements ", top_level_numbers)
    return false
end

function compareIndex(sn1 :: RHSVar, sn2 :: RHSVar)
    return toLHSVar(sn1) == toLHSVar(sn2)
end

function compareIndex(s1::ANY, s2::ANY)
    return s1 == s2
end

function create_merged_output_from_map(output_map, unique_id, state, sym_to_type, loweredAliasMap)
    @dprintln(3,"create_merged_output_from_map, output_map = ", output_map, " sym_to_type = ", sym_to_type)
    # If there are no outputs then return nothing.
    if length(output_map) == 0
        return (nothing, [], true, nothing, [])
    end

    # If there is only one output then all we need is the symbol to return.
    if length(output_map) == 1
        for i in output_map
            new_lhs = toRHSVar(i[1], sym_to_type[i[1]], state.LambdaVarInfo)
            new_rhs = toRHSVar(getAliasMap(loweredAliasMap, i[2]), sym_to_type[i[2]], state.LambdaVarInfo)
            return (new_lhs, [new_lhs], true, [new_rhs])
        end
    end

    lhs_order = RHSVar[]
    rhs_order = RHSVar[]
    for i in output_map
        @dprintln(3,"Working on output ", i, " lhs = ", i[1], " rhs = ", i[2])
        lhs_rhsvar = toRHSVar(i[1], sym_to_type[i[1]], state.LambdaVarInfo)
        push!(lhs_order, lhs_rhsvar)
        @dprintln(3,"lhs_rhsvar = ", lhs_rhsvar)
        rhs_rhsvar = toRHSVar(getAliasMap(loweredAliasMap, i[2]), sym_to_type[i[2]], state.LambdaVarInfo)
        push!(rhs_order, rhs_rhsvar)
        @dprintln(3,"rhs_rhsvar = ", rhs_rhsvar)
    end
    num_map = length(lhs_order)

    # Multiple outputs.

    # First, form the type of the tuple for those multiple outputs.
    tt = Expr(:tuple)
    for i = 1:num_map
        push!(tt.args, CompilerTools.LambdaHandling.getType(rhs_order[i], state.LambdaVarInfo))
    end
    temp_type = eval(tt)

    ( createRetTupleType(lhs_order, unique_id, state), lhs_order, false, rhs_order )
end

"""
Test whether we can fuse the two most recent parfor statements and if so to perform that fusion.
"""
function fuse(body, body_index, cur::Expr, state)
    global fuse_limit
    prev = body[body_index]

    # Handle the debugging case where we want to limit the amount of parfor fusion to a certain number.
    if fuse_limit == 0
        return false
    end
    if fuse_limit > 0
        global fuse_limit = fuse_limit - 1
    end

    @dprintln(2, "Testing if fusion is possible.")
    prev_parfor = getParforNode(prev)
    cur_parfor  = getParforNode(cur)

    sym_to_type = Dict{LHSVar, DataType}()

    @dprintln(2, "prev = ", prev)
    @dprintln(2, "cur = ", cur)
    @dprintln(2, "prev.typ = ", prev.typ)
    @dprintln(2, "cur.typ = ", cur.typ)

    prev_assignment = isAssignmentNode(prev)

    out_correlation = getParforCorrelation(prev_parfor, state)
    if out_correlation == nothing
        @dprintln(3, "Fusion will not happen because out_correlation = nothing")
        return false
    end
    in_correlation  = getParforCorrelation(cur_parfor, state)
    if in_correlation == nothing
        @dprintln(3, "Fusion will not happen because in_correlation = nothing")
        return false
    end
    @dprintln(3,"Fusion correlations ", out_correlation, " ", in_correlation)

    is_prev_multi = isFusionAssignment(prev)
    is_cur_multi  = isFusionAssignment(cur)

    prev_num_dims = length(prev_parfor.loopNests)
    cur_num_dims  = length(cur_parfor.loopNests)

    map_prev_lhs_post, map_prev_lhs_reduction = createMapLhsToParfor(prev, prev_parfor, is_prev_multi, sym_to_type, state)
    map_prev_lhs_all = merge(map_prev_lhs_post, map_prev_lhs_reduction)
    map_cur_lhs_post,  map_cur_lhs_reduction  = createMapLhsToParfor(cur,  cur_parfor,  is_cur_multi,  sym_to_type, state)
    map_cur_lhs_all  = merge(map_cur_lhs_post, map_cur_lhs_reduction)
    prev_output_arrays = collect(values(map_prev_lhs_post))
    prev_output_reduce = collect(values(map_prev_lhs_reduction))
    cur_output_arrays  = collect(values(map_cur_lhs_post))
    cur_output_reduce  = collect(values(map_cur_lhs_reduction))

    merge!(prev_parfor.array_aliases, map_prev_lhs_post)
    assert(length(cur_parfor.array_aliases) == 0)

    loweredAliasMap = createLoweredAliasMap(prev_parfor.array_aliases)

    @dprintln(3, "map_prev_lhs_post = ", map_prev_lhs_post)
    @dprintln(3, "map_prev_lhs_reduction = ", map_prev_lhs_reduction)
    @dprintln(3, "map_cur_lhs_post = ", map_cur_lhs_post)
    @dprintln(3, "map_cur_lhs_reduction = ", map_cur_lhs_reduction)
    @dprintln(3, "sym_to_type = ", sym_to_type)
    reduction_var_used = isSymbolsUsed(map_prev_lhs_reduction, cur_parfor.top_level_number, state)
    @dprintln(3, "reduction_var_used = ", reduction_var_used)
    prev_iei = iterations_equals_inputs(prev_parfor)
    cur_iei  = iterations_equals_inputs(cur_parfor)
    @dprintln(3, "iterations_equals_inputs prev and cur = ", prev_iei, " ", cur_iei)
    @dprintln(3, "loweredAliasMap = ", loweredAliasMap)
    @dprintln(3, "cur_parfor.arrays_read_past_index = ", cur_parfor.arrays_read_past_index)
    arrays_non_simply_indexed_in_cur_that_access_prev_output = intersect(cur_parfor.arrays_read_past_index, prev_output_arrays)
    @dprintln(3, "arrays_non_simply_indexed_in_cur_that_access_prev_output = ", arrays_non_simply_indexed_in_cur_that_access_prev_output)
    @dprintln(3, "prev_parfor.arrays_written_past_index = ", prev_parfor.arrays_written_past_index)
    # Compute which scalars and arrays are ever read or written by the body of the parfor
    prev_rws = CompilerTools.ReadWriteSet.from_exprs(prev_parfor.body, pir_rws_cb, state.LambdaVarInfo, state.LambdaVarInfo)
    # Compute which scalars and arrays are ever read or written by the body of the parfor
    cur_rws = CompilerTools.ReadWriteSet.from_exprs(cur_parfor.body, pir_rws_cb, state.LambdaVarInfo, state.LambdaVarInfo)
    cur_accessed = CompilerTools.ReadWriteSet.getArraysAccessed(cur_rws)
    arrays_non_simply_indexed_in_prev_that_are_read_in_cur = intersect(prev_parfor.arrays_written_past_index, cur_accessed)
    @dprintln(3, "arrays_non_simply_indexed_in_prev_that_are_read_in_cur = ", arrays_non_simply_indexed_in_prev_that_are_read_in_cur)
    if !isempty(arrays_non_simply_indexed_in_cur_that_access_prev_output)
        @dprintln(1, "Fusion won't happen because the second parfor accesses some array created by the first parfor in a non-simple way.")
    end
    if !isempty(arrays_non_simply_indexed_in_prev_that_are_read_in_cur)
        @dprintln(1, "Fusion won't happen because the first parfor wrote to some array index in a non-simple way that the second parfor needs to access.")
    end

    @dprintln(3, "Fusion prev_rws = ", prev_rws)
    @dprintln(3, "Fusion cur_rws = ", cur_rws)
    prev_write_set_array = Set(collect(keys(prev_rws.writeSet.arrays)))
    #cur_read_set_array   = Set(collect(keys(cur_rws.readSet.arrays)))
    #cur_read_with_aliases = getAllAliases(cur_read_set_array, map_prev_lhs_post)
    @dprintln(3, "Fusion prev_write = ", prev_write_set_array)
    #write_in_prev_read_in_cur = intersect(prev_write_with_aliases, cur_read_with_aliases)
    #@dprintln(3, "Fusion intersection = ", write_in_prev_read_in_cur)
    found_array_access_problem = false
    @dprintln(3, "cur_parfor.loopNests = ", cur_parfor.loopNests)
    for cur_read_array in collect(keys(cur_rws.readSet.arrays))
        with_aliases = getAllAliases(Set{LHSVar}([cur_read_array]), map_prev_lhs_post)
        @dprintln(3, "Checking array ", cur_read_array, " with aliases = ", with_aliases)
        if isempty(intersect(prev_write_set_array, with_aliases))
            continue
        end
        @dprintln(3, "Array  ", cur_read_array, " produced by prev parfor.")
        accesses = cur_rws.readSet.arrays[cur_read_array]
        for index_expr in accesses
            @dprintln(3, "Checking usage  ", index_expr)
            for index = 1:length(index_expr)
                lnindex = length(cur_parfor.loopNests) + 1 - index
                if lnindex < 1 || lnindex > length(cur_parfor.loopNests)
                    @dprintln(3, "Found index expression not just using this dimension's loop index variable")
                    found_array_access_problem = true
                    continue
                end
                good_index = cur_parfor.loopNests[lnindex].indexVariable
                @dprintln(3, "Good index here = ", good_index, " cur_index = ", index_expr[index])
                if !compareIndex(good_index, index_expr[index])
                    @dprintln(3, "Found index expression not just using this dimension's loop index variable")
                    found_array_access_problem = true
                end
            end
        end
    end

    @dprintln(3, (prev_iei ,
        cur_iei  ,
        out_correlation == in_correlation ,
        !reduction_var_used ,
        isempty(arrays_non_simply_indexed_in_cur_that_access_prev_output) ,
        isempty(arrays_non_simply_indexed_in_prev_that_are_read_in_cur) ,
        (prev_num_dims == cur_num_dims) ,
        !found_array_access_problem        ))

    if  prev_iei &&
        cur_iei  &&
        out_correlation == in_correlation &&
        !reduction_var_used &&
        isempty(arrays_non_simply_indexed_in_cur_that_access_prev_output) &&
        isempty(arrays_non_simply_indexed_in_prev_that_are_read_in_cur) &&
        (prev_num_dims == cur_num_dims) &&
        !found_array_access_problem

        @dprintln(3, "Fusion will happen here.")

        # Get the top-level statement for the previous parfor.
        prev_stmt_live_first = CompilerTools.LivenessAnalysis.find_top_number(prev_parfor.top_level_number[1], state.block_lives)
        assert(prev_stmt_live_first != nothing)
        @dprintln(2,"Prev parfor first = ", prev_stmt_live_first)
        prev_stmt_live_last = CompilerTools.LivenessAnalysis.find_top_number(prev_parfor.top_level_number[end], state.block_lives)
        assert(prev_stmt_live_last != nothing)
        @dprintln(2,"Prev parfor last = ", prev_stmt_live_last)

        # Get the top-level statement for the current parfor.
        assert(length(cur_parfor.top_level_number) == 1)
        cur_stmt_live  = CompilerTools.LivenessAnalysis.find_top_number(cur_parfor.top_level_number[1],  state.block_lives)
        assert(cur_stmt_live != nothing)
        @dprintln(2,"Cur parfor = ", cur_stmt_live)

        # Get the variables live after the previous parfor.
        live_in_prev = prev_stmt_live_first.live_in
        # def_prev     = prev_stmt_live_first.def
        def_prev = Set{LHSVar}()
        # The "def" for a fused parfor is to a first approximation the union of the "def" of each statement in the parfor.
        # Some variables will be eliminated by fusion so we follow this first approximation up by an intersection with the
        # live_out of the last statement of the previous parfor.
        for prev_parfor_stmt in prev_parfor.top_level_number
            constituent_stmt_live = CompilerTools.LivenessAnalysis.find_top_number(prev_parfor_stmt, state.block_lives)
            def_prev = union(def_prev, constituent_stmt_live.def)
        end
        def_prev = intersect(def_prev, prev_stmt_live_last.live_out)

        @dprintln(2,"live_in_prev = ", live_in_prev, " def_prev = ", def_prev)

        # Get the variables live after the previous parfor.
        live_out_prev = prev_stmt_live_last.live_out
        @dprintln(2,"live_out_prev = ", live_out_prev)

        # Get the live variables into the current parfor.
        live_in_cur   = cur_stmt_live.live_in
        def_cur       = cur_stmt_live.def
        @dprintln(2,"live_in_cur = ", live_in_cur, " def_cur = ", def_cur)

        # Get the variables live after the current parfor.
        live_out_cur  = cur_stmt_live.live_out
        @dprintln(2,"live_out_cur = ", live_out_cur)

        surviving_def = intersect(union(def_prev, def_cur), live_out_cur)

        new_in_prev = setdiff(live_out_prev, live_in_prev)
        new_in_cur  = setdiff(live_out_cur,  live_in_cur)
        @dprintln(3,"new_in_prev = ", new_in_prev)
        @dprintln(3,"new_in_cur = ", new_in_cur)
        @dprintln(3,"surviving_def = ", surviving_def)

        # The things that come in live to cur but leave it dead.
        not_used_after_cur = setdiff(live_out_prev, live_out_cur)
        @dprintln(2,"not_used_after_cur = ", not_used_after_cur)

        live_out_prev_aliases = getAllAliases(live_out_prev, prev_parfor.array_aliases)
        live_out_cur_aliases  = getAllAliases(live_out_cur, prev_parfor.array_aliases)
        @dprintln(2, "live_out_prev_aliases = ", live_out_prev_aliases)
        @dprintln(2, "live_out_cur_aliases  = ", live_out_cur_aliases)
        not_used_after_cur_with_aliases = setdiff(live_out_prev_aliases, live_out_cur_aliases)
        @dprintln(2,"not_used_after_cur_with_aliases = ", not_used_after_cur_with_aliases)

        unique_id = get_unique_num() # prev_parfor.unique_id

        # Output of this parfor are the new things in the current parfor plus the new things in the previous parfor
        # that don't die during the current parfor.
        output_map = Dict{LHSVar, LHSVar}()
        for i in map_prev_lhs_all
            if !in(i[1], not_used_after_cur)
                output_map[i[1]] = i[2]
            end
        end
        for i in map_cur_lhs_all
            output_map[i[1]] = i[2]
        end

        new_aliases = Dict{LHSVar, LHSVar}()
        for i in map_prev_lhs_post
            if !in(i[1], not_used_after_cur)
                new_aliases[i[1]] = i[2]
            end
        end

        @dprintln(3,"output_map = ", output_map)
        @dprintln(3,"new_aliases = ", new_aliases)

        # return code 2 if there is no output in the fused parfor
        # this means the parfor is dead and should be removed
        if length(surviving_def)==0
            @dprintln(1,"No output for the fused parfor so the parfor is dead and is being removed.")
            return 2;
        end

        first_arraylen = getFirstArrayLens(prev_parfor, prev_num_dims, state)
        @dprintln(3,"first_arraylen = ", first_arraylen)

        # Merge each part of the two parfor nodes.

        # loopNests - nothing needs to be done to the loopNests
        # But we use them to establish a mapping between corresponding indices in the two parfors.
        # Then, we use that map to convert indices in the second parfor to the corresponding ones in the first parfor.
        index_map = Dict{LHSVar, LHSVar}()
        assert(length(prev_parfor.loopNests) == length(cur_parfor.loopNests))
        for i = 1:length(prev_parfor.loopNests)
            index_map[toLHSVar(cur_parfor.loopNests[i].indexVariable)] = toLHSVar(prev_parfor.loopNests[i].indexVariable)
        end

        @dprintln(3,"array_aliases before merge ", prev_parfor.array_aliases)
        for i in map_cur_lhs_post
            from = i[1]
            to   = i[2]
            prev_parfor.array_aliases[from] = to
        end
        @dprintln(3,"array_aliases after merge ", prev_parfor.array_aliases)

        # postParFor - can merge everything but the last entry in the postParFor's.
        # The last entries themselves need to be merged extending the tuple if the prev parfor's output stays live and
        # just keeping the cur parfor output if the prev parfor output dies.
        (new_lhs, all_rets, single, output_items) = create_merged_output_from_map(output_map, unique_id, state, sym_to_type, loweredAliasMap)
        @dprintln(3,"new_lhs = ", new_lhs)
        @dprintln(3,"all_rets = ", all_rets)
        @dprintln(3,"single = ", single)
        @dprintln(3,"output_items = ", output_items)
        output_items_set = live_out_cur
        output_items_with_aliases = getAllAliases(output_items_set, prev_parfor.array_aliases)

        @dprintln(3,"output_items_set = ", output_items_set)
        @dprintln(3,"output_items_with_aliases = ", output_items_with_aliases)

        escaping_sets = Set{LHSVar}()
        # Create a dictionary of arrays to the last variable containing the array's value at the current index space.
        save_body = prev_parfor.body
        arrayset_dict = Dict{LHSVar, RHSVar}()
        for i = 1:length(save_body)
            x = save_body[i]
            if isArraysetCall(x)
                # Here we have a call to arrayset.
                array_name = x.args[2]
                value      = x.args[3]
                array_var  = toLHSVar(array_name)
                assert(isa(array_name, RHSVar))
                assert(isa(value, RHSVar))
                array_var = toLHSVar(array_name)
                arrayset_dict[array_var] = value
                if getDesc(array_var, state.LambdaVarInfo) & ISCAPTURED == ISCAPTURED
                    push!(escaping_sets, array_var)
                end
            elseif typeof(x) == Expr && x.head == :(=)
                lhs = x.args[1]
                rhs = x.args[2]
                #assert(isa(lhs, RHSVar))
                if isArrayrefCall(rhs)
                    array_name = rhs.args[2]
                    assert(isa(array_name, RHSVar))
                    array_var = toLHSVar(array_name)
#                    arrayset_dict[array_var] = lhs
                    if getDesc(array_var, state.LambdaVarInfo) & ISCAPTURED == ISCAPTURED
                        push!(escaping_sets, array_var)
                    end
                end
            end
        end
        @dprintln(3,"arrayset_dict = ", arrayset_dict)

        # Extend the arrayset_dict to include the lhs of the prev parfor.
        for i in map_prev_lhs_post
            lhs_sym = i[1]
            rhs_sym = i[2]
            arrayset_dict[lhs_sym] = arrayset_dict[rhs_sym]
            if getDesc(lhs_sym, state.LambdaVarInfo) & ISCAPTURED == ISCAPTURED
                push!(escaping_sets, lhs_sym)
            end
        end
        @dprintln(3,"arrayset_dict = ", arrayset_dict)
        # make sure escaping sets include aliases and prev's output_items
        escaping_sets = union(escaping_sets, not_used_after_cur_with_aliases)
        for (k,v) in prev_parfor.array_aliases
            if in(k, escaping_sets)
                push!(escaping_sets, toLHSVar(v))
            end
        end
        @dprintln(3,"escaping_sets = ", escaping_sets)

        # body - Append cur body to prev body but replace arrayset's in prev with a temp variable
        # and replace arrayref's in cur with the same temp.
        arrays_set_in_cur_body = Set{LHSVar}()
        # Convert the cur_parfor body.
        new_cur_body = map(x -> substitute_cur_body(x, arrayset_dict, index_map, arrays_set_in_cur_body, loweredAliasMap, state), cur_parfor.body)
        arrays_set_in_cur_body_with_aliases = getAllAliases(arrays_set_in_cur_body, prev_parfor.array_aliases)
        @dprintln(3,"arrays_set_in_cur_body = ", arrays_set_in_cur_body)
        @dprintln(3,"arrays_set_in_cur_body_with_aliases = ", arrays_set_in_cur_body_with_aliases)
        combined = union(arrays_set_in_cur_body_with_aliases, not_used_after_cur_with_aliases)
        @dprintln(2,"combined = ", combined)

        prev_parfor.body = Any[]
        for i = 1:length(save_body)
            new_body_line = substitute_arrayset(save_body[i], combined, output_items_with_aliases, escaping_sets)
            @dprintln(3,"new_body_line = ", new_body_line)
            if new_body_line != nothing
                push!(prev_parfor.body, new_body_line)
            end
        end
        append!(prev_parfor.body, new_cur_body)
        @dprintln(2,"New body = ")
        printBody(2, prev_parfor.body)

        # preParFor - Append cur preParFor to prev parParFor but eliminate array creation from
        # prevParFor where the array is in allocs_to_eliminate.
        removed_allocs = Set{LHSVar}()
        @dprintln(3, "preParFor cur prestatements are array input = ", parforArrayInput(prev_parfor))
        filtered_cur_pre = filter(x -> !is_eliminated_arraylen(x), cur_parfor.preParFor)
        prev_parfor.preParFor = [ filter(x -> !is_eliminated_allocation_map(x, output_items_with_aliases, removed_allocs), prev_parfor.preParFor);
                                  map(x -> substitute_arraylen(x, first_arraylen) , filtered_cur_pre)]
#                                  parforArrayInput(prev_parfor) ? map(x -> substitute_arraylen(x, first_arraylen) , filtered_cur_pre) : cur_parfor.preParFor]
        @dprintln(3,"removed_allocs = ", removed_allocs)
        filter!( x -> !is_eliminated_arraysize(x, removed_allocs, prev_parfor.array_aliases), prev_parfor.preParFor)
        @dprintln(2,"New preParFor = ", prev_parfor.preParFor)

        append!(prev_parfor.hoisted, cur_parfor.hoisted)
        @dprintln(2,"New hoisted = ", prev_parfor.hoisted)

        # if allocation of an array is removed, arrayset should be removed as well since the array doesn't exist anymore
        @dprintln(4,"prev_parfor.body before removing dead arrayset: ", prev_parfor.body)
        filter!( x -> !is_dead_arrayset(x, removed_allocs), prev_parfor.body)
        @dprintln(4,"prev_parfor.body after_ removing dead arrayset: ", prev_parfor.body)

        # reductions - a simple append with the caveat that you can't fuse parfor where the first has a reduction that the second one uses
        # need to check this caveat above.
        append!(prev_parfor.reductions, cur_parfor.reductions)
        @dprintln(2,"New reductions = ", prev_parfor.reductions)

        prev_parfor.postParFor = [ prev_parfor.postParFor[1:end-1]; cur_parfor.postParFor[1:end-1]]
        @dprintln(2,"New initial postParFor = ", prev_parfor.postParFor)
        push!(prev_parfor.postParFor, oneIfOnly(output_items))
        @dprintln(2,"New postParFor = ", prev_parfor.postParFor, " typeof(postParFor) = ", typeof(prev_parfor.postParFor), " ", typeof(prev_parfor.postParFor[end]))

        # original_domain_nodes - simple append
        append!(prev_parfor.original_domain_nodes, cur_parfor.original_domain_nodes)
        @dprintln(2,"New domain nodes = ", prev_parfor.original_domain_nodes)

        # top_level_number - what to do here? is this right?
        push!(prev_parfor.top_level_number, cur_parfor.top_level_number[1])

        @dprintln(3,"New lhs = ", new_lhs)
        if prev_assignment
            # The prev parfor was of the form "var = parfor(...)".
            if new_lhs != nothing
                @dprintln(2,"prev was assignment and is staying an assignment")
                # The new lhs is not empty and so this is the normal case where "prev" stays an assignment expression and we update types here and if necessary FusionSentinel.
                prev.args[1] = new_lhs
                prev.typ = getType(new_lhs, state.LambdaVarInfo)
                prev.args[2].typ = prev.typ
                # Strip off a previous FusionSentinel() if it exists in the expression.
                prev.args = prev.args[1:2]
                if !single
                    push!(prev.args, FusionSentinel())
                    append!(prev.args, all_rets)
                    @dprintln(3,"New multiple ret prev args is ", prev.args)
                end
            else
                @dprintln(2,"prev was assignment and is becoming bare")
                # The new lhs is empty and so we need to transform "prev" into an assignment expression.
                body[body_index] = TypedExpr(nothing, :parfor, prev_parfor)
            end
        else
            # The prev parfor was a bare-parfor (not in the context of an assignment).
            if new_lhs != nothing
                @dprintln(2,"prev was bare and is becoming an assignment")
                # The new lhs is not empty so the fused parfor will not be bare and "prev" needs to become an assignment expression.
                body[body_index] = mk_assignment_expr(new_lhs, prev, state)
                prev = body[body_index]
                prev.args[2].typ = CompilerTools.LambdaHandling.getType(new_lhs, state.LambdaVarInfo)
                if !single
                    push!(prev.args, FusionSentinel())
                    append!(prev.args, all_rets)
                    @dprintln(3,"New multiple ret prev args is ", prev.args)
                end
            else
                @dprintln(2,"prev was bare and is staying bare")
            end
        end

        @dprintln(2,"New parfor = ", prev_parfor)

        #throw(string("not finished"))

        return 1
    else
        @dprintln(3, "Fusion could not happen here.")
    end

    return 0

    false
end

"""
Performs the mmap to mmap! phase.
If the arguments of a mmap dies aftewards, and is not aliased, then
we can safely change the mmap to mmap!.
"""
function mmapToMmap!(LambdaVarInfo, body::Expr, lives, uniqSet)
    assert(body.head == :body)
    # For each statement in the body.
    for i =1:length(body.args)
        expr = body.args[i]
        # If the statement is an assignment.
        if isa(expr, Expr) && (expr.head === :(=))
            lhs = toLHSVar(expr.args[1], LambdaVarInfo)
            rhs = expr.args[2]
            # right now assume all
            assert(isa(lhs, RHSVar))
            lhsTyp = CompilerTools.LambdaHandling.getType(lhs, LambdaVarInfo)
            # If the right-hand side is an mmap.
            if isa(rhs, Expr) && (rhs.head === :mmap)
                args = rhs.args[1]
                tls = CompilerTools.LivenessAnalysis.find_top_number(i, lives)
                assert(tls != nothing)
                assert(CompilerTools.LivenessAnalysis.isDef(lhs, tls))
                @dprintln(4, "mmap lhs=", lhs, " args=", args, " live_out = ", tls.live_out)
                reuse = nothing
                j = 0
                # Find some input array to the mmap that is dead after this statement.
                while j < length(args)
                    j = j + 1
                    if isa(args[j], RHSVar)
                        v = toLHSVar(args[j])
                        if isa(v, LHSVar) && !in(v, tls.live_out) && in(v, uniqSet) &&
                            CompilerTools.LambdaHandling.getType(v, LambdaVarInfo) == lhsTyp
                            reuse = v  # Found a dying symbol.
                            break
                        end
                    end
                end
                # If we found a dying array whose space we can reuse.
                if !(reuse === nothing)
                    rhs.head = :mmap!   # Change to mmap!
                    @dprintln(2, "mmapToMMap!: successfully reuse ", reuse, " for ", lhs)
                    if j != 1  # The array to reuse has to be first.  If it isn't already then reorder the args to make it so.
                        # swap j-th and 1st argument
                        rhs.args[1] = DomainIR.arraySwap(rhs.args[1], 1, j)
                        rhs.args[2] = DomainIR.lambdaSwapArg(rhs.args[2], 1, j)
                        @dprintln(3, "mmapToMMap!: after swap, ", lhs, " = ", rhs)
                    end
                end
            end
        end
    end
end

"""
For every basic block, try to push domain IR statements down and non-domain IR statements up so that domain nodes
are next to each other and can be fused.
"""
function maxFusion(bl :: CompilerTools.LivenessAnalysis.BlockLiveness)
    # We will try to optimize the order in each basic block.
    for bb in collect(values(bl.basic_blocks))
        if false
            # One approach to this problem is to create a dependency graph and then use that to calculate valid reorderings
            # that maximize fusion.  We may still want to switch to this at some point but it is more complicated and the
            # simpler approach works for now.

            # Start with a pseudo-statement corresponding to the variables live-in to this basic block.
            livein     = StatementWithDeps(nothing)
            # The variable "last_write" is a dictionary mapping a symbol to the last statement to have written that variable.
            last_write = Dict{Symbol, StatementWithDeps}()
            # For each symbol live-in to this basic block, add it to "last_write".
            for def in bb.live_in
                assert(typeof(def) == Symbol)
                last_write[def] = livein
            end

            # Create the dependency graph.
            for i = 1:length(bb.statements)
                # Get the i'th statement in this basic block.
                stmt    = bb.statements[i]
                # Create an object to associate this statement with its dependencies.
                new_swd = StatementWithDeps(stmt)

                # Establish dependency links from prior statements to this one based on the last statement to have written each symbol used in this statement.
                for use in stmt.use
                    assert(typeof(use) == Symbol)
                    assert(haskey(last_write, use))
                    push!(last_write[use].deps, new_swd)
                end
                # Update last_write with any symbols def in this statement.
                for def in stmt.def
                    assert(typeof(def) == Symbol)
                    last_write[def] = new_swd
                end
            end

            topo_sort = StatementWithDeps[]
            dfsVisit(livein, 1, topo_sort)
        else
            @dprintln(3, "Doing maxFusion in block ", bb)
            # A bubble-sort style of coalescing domain nodes together in the AST.
            earliest_parfor = 1
            found_change = true

            # While the lastest pass over the AST created some change, keep searching for more interchanges that can coalesce domain nodes.
            while found_change
                found_change = false

                # earliest_parfor is an optimization that we don't have to scan every statement in the block but only those statements from
                # the first parfor to the last statement in the block.
                i = earliest_parfor
                # @dprintln(3,"bb.statements = ", bb.statements)
                earliest_parfor = length(bb.statements)

                while i < length(bb.statements)
                    cur  = bb.statements[i]
                    next = bb.statements[i+1]
                    cannot_move_next = mustRemainLastStatementInBlock(next.tls.expr)
                    @dprintln(3,"maxFusion cur = ", cur.tls.expr)
                    @dprintln(3,"maxFusion next = ", next.tls.expr)
                    @dprintln(3,"maxFusion next.use = ", next.use)
                    cur_domain_node  = isDomainNode(cur.tls.expr)
                    next_domain_node = isDomainNode(next.tls.expr)
                    intersection     = intersect(cur.def, next.use)
                    @dprintln(3,"cur_domain_node = ", cur_domain_node, " next_domain_node = ", next_domain_node, " intersection = ", intersection)
                    if cur_domain_node && !cannot_move_next
                        if !next_domain_node
                          if isempty(intersection)
                            # If the current statement is a domain node and the next staterment isn't and we are allowed to move the next node
                            # in the block and the next statement doesn't use anything produced by this statement then we can switch the order of
                            # the current and next statement.
                            @dprintln(3,"bubbling domain node down")
                            (bb.statements[i], bb.statements[i+1]) = (bb.statements[i+1], bb.statements[i])
                            (bb.cfgbb.statements[i], bb.cfgbb.statements[i+1]) = (bb.cfgbb.statements[i+1], bb.cfgbb.statements[i])
                            (bb.cfgbb.statements[i].index, bb.cfgbb.statements[i+1].index) = (bb.cfgbb.statements[i+1].index, bb.cfgbb.statements[i].index)
                            found_change = true
                          else # intersection is not empty, but let's check if it is no longer live
                            if length(intersection) == 1 && isAssignmentNode(cur.tls.expr) && isAssignmentNode(next.tls.expr)
                                lhsvar = toLHSVar(cur.tls.expr.args[1])
                                tmpvar = toLHSVar(next.tls.expr.args[2])
                                if lhsvar == tmpvar && in(tmpvar, intersection) && !in(tmpvar, next.live_out)
                                    @dprintln(3, "next is assignment, and RHS is not live afterwards")
                                    (bb.statements[i], bb.statements[i+1]) = (bb.statements[i+1], bb.statements[i])
                                    (bb.cfgbb.statements[i], bb.cfgbb.statements[i+1]) = (bb.cfgbb.statements[i+1], bb.cfgbb.statements[i])
                                    (bb.cfgbb.statements[i].index, bb.cfgbb.statements[i+1].index) = (bb.cfgbb.statements[i+1].index, bb.cfgbb.statements[i].index)
                                    found_change = true
                                    bb.statements[i+1].def = bb.statements[i].def
                                    bb.statements[i+1].tls.expr.args[1] = bb.statements[i].tls.expr.args[1]
                                    bb.statements[i].tls.expr = nothing
                                    @dprintln(3, "after rewrite, cur = ", bb.statements[i+1])
                                end
                            end
                          end
                        else
                            if i < earliest_parfor
                                earliest_parfor = i
                            end
                        end
                    end
                    i += 1
                end
            end
        end
    end
end


text:
using LightGraphs
using LinearAlgebra: UpperTriangular

"""
UniDirectionalBeta
"""
function beta_graph(beta::Vector; seed::Integer=-1)
    n = length(beta)
    coefM = UpperTriangular(beta' .+ beta)
    g = SimpleGraph(n)
    for i in 1:n
        for j in i:n
            expcoef_ij = exp(coefM[i,j])
            if rand() < (expcoef_ij / (1 + expcoef_ij))
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

"""
BiDirectionalBeta
"""
function beta_graph(alpha::Vector, beta::Vector; seed::Integer=-1)
    n = length(beta)
    @assert length(alpha) == n "alpha and beta must be the same length"
    coefM = alpha' .+ beta
    g = SimpleDiGraph(n)
    for i in 1:n
        for j in 1:n
            expcoef_ij = exp(coefM[i,j])
            if rand() < (expcoef_ij / (1 + expcoef_ij))
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

"""
UniDirectionalSparseBeta
"""
function sparse_beta_graph(beta::Vector, mu::Real; seed::Integer=-1)
    @assert all(beta .>= 0) "beta must be non-negative"
    @assert 0 in beta "beta must contain at least one zero"
    n = length(beta)
    coefM = UpperTriangular(beta' .+ beta)
    g = SimpleGraph(n)
    for i in 1:n
        for j in i:n
            expcoef_ij = exp(mu + coefM[i,j])
            if rand() < (expcoef_ij / (1 + expcoef_ij))
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

max_stars_repo_name: jacobusmmsmit/sparse-beta

max_issues_repo_name: jacobusmmsmit/sparse-beta

max_forks_repo_name: jacobusmmsmit/sparse-beta


text:
function trap_rule(xarr::Vector{T},farr::Vector{T}) where T
    N = length(xarr)

    Δx = xarr[2] - xarr[1]
    quad = zero(T)
    for ti=1:N-1
        quad += farr[ti] + farr[ti+1]
    end

    quad *= 0.5*Δx

    return quad
end

# return J. Compares how well particle field (represented by θ)
# approximates F(y)
function cost(θ,y,F)
    x,s,h = decompose_theta(θ)
    G = get_gauss_field(θ,y)
    f = (G - F).^2
    J = trap_rule(y,f)
    return J
end

function grad_cost(θ,y,F)
    dJdx = compute_dJdx(θ,y,F)
    dJds = compute_dJds(θ,y,F)
    dJdh = compute_dJdh(θ,y,F)
    dJdθ = vcat(dJdx,dJds,dJdh)

    return dJdθ
end

function compute_dJdx(θ,y,F)
    x,s,h = decompose_theta(θ)
    G = get_gauss_field(θ,y)

    N = length(x)
    M = length(y)
    dJdx = zeros(N)

    for ti=1:N
        f = zeros(M)
        for tj=1:M
            f[tj] = 2*(G[tj] - F[tj])*compute_dGdx(x[ti],s[ti],h,y[tj])
        end
        dJdx[ti] = trap_rule(y,f)
    end
    return dJdx
end

function compute_dJds(θ,y,F)
    x,s,h = decompose_theta(θ)
    G = get_gauss_field(θ,y)

    N = length(x)
    M = length(y)
    dJds = zeros(N)

    for ti=1:N
        f = zeros(M)
        for tj=1:M
            f[tj] = 2*(G[tj] - F[tj])*compute_dGds(x[ti],s[ti],h,y[tj])
        end
        dJds[ti] = trap_rule(y,f)
    end
    return dJds
end

function compute_dJdh(θ,y,F)
    x,s,h = decompose_theta(θ)
    G = get_gauss_field(θ,y)

    M = length(y)
    f = zeros(M)
    for tj=1:M
        f[tj] = 2*(G[tj] - F[tj])*compute_dGdh(x,s,h,y[tj])
    end
    dJdh = trap_rule(y,f)
    return dJdh
end

function compute_dGdx(x::T,s::T,h::T,y::T) where T
    dGdx = -2*s*(x - y)*GaussKernel(x-y,h)/h^2
    return dGdx
end

function compute_dGds(x::T,s::T,h::T,y::T) where T
    dGds = GaussKernel(x-y,h)
    return dGds
end

function compute_dGdh(x::Vector{T},s::Vector{T},h::T,y::T) where T
    N = length(x)
    dGdh = zero(T)
    for ti=1:N
        dGdh += s[ti]*(x[ti] - y)^2*GaussKernel(x[ti]-y,h)
    end
    dGdh *= 2/h^3
    return dGdh
end

# computes grad_cost with finite difference method
# to brute-force check the analytic gradient
function grad_cost_FD(θ,y,F)
    n_temp = length(θ)
    N = Int((n_temp-1)/2)

    Δx = 0.01*(maximum(y) -minimum(y))
    Δs = 0.01*maximum(F)
    Δh = Δx
    dJdx = zeros(N)
    dJds = zeros(N)
    dJdh = 0.0

    # compute dJdx
    for ti=1:N
        θtemp = copy(θ)
        θtemp[ti] = θtemp[ti] + Δx
        J1 = cost(θtemp,y,F)
        J0 = cost(θ,y,F)
        dJdx[ti] = (J1 - J0)/Δx
    end

    # compute dJds
    for ti=1:N
        θtemp = copy(θ)
        θtemp[ti+N] = θtemp[ti+N] + Δs
        J1 = cost(θtemp,y,F)
        J0 = cost(θ,y,F)
        dJds[ti] = (J1 - J0)/Δs
    end

    # compute dJdh
    θtemp = copy(θ)
    θtemp[2N+1] = θtemp[2N+1] + Δh
    J1 = cost(θtemp,y,F)
    J0 = cost(θ,y,F)
    dJdh = (J1 - J0)/Δh

    dJdθ = vcat(dJdx,dJds,dJdh)
    return dJdθ
end

function plot_fields(θ,y,F;iter=-1)
    x,s,h = decompose_theta(θ)

    G = get_gauss_field(θ,y)
    p = plot()
    if iter == -1
        plot!(p,y,G,label="G")
    else
        plot!(p,y,G,label="G",title="iteration $(iter)")
    end
    plot!(p,y,F,label="F")
    gui(p)

    nothing
end

function get_gauss_field(θ,y)
    # decompose theta
    x,s,h = decompose_theta(θ)

    # plot gauss field from particles
    d = 1                                   # dimension of data
    M = length(y)                           # number of targets
    N = length(x)                           # number of sources
    ε = 1e-2                                # max tolerance
    W = 1                                   # number of sources being evaluated
    G = zeros(M)                            # output values array
    mydgt!(G,d,M,N,h,ε,x,y,s,W)
    return G
end

function gradient_descent_step!(θ,∇f,α)
    N = length(θ)
    for ti=1:N
        θ[ti] -= α*∇f[ti]
    end
    nothing
end

function gradient_descent(θ0,y,F;NiterMAX=100,TOL=1e-10,pause=false)
    error = 100.0
    α = 0.2

    iter = 1
    θ = θ0
    Jtemp = 100.0
    while iter<NiterMAX && error > TOL
        println("\riteration ",iter,"          ")
        plot_fields(θ,y,F;iter=iter)

        J = cost(θ,y,F)
        error = abs(J - Jtemp)

        #∇J = grad_cost_FD(θ,y,F)
        ∇J = grad_cost(θ,y,F)

        gradient_descent_step!(θ,∇J,α)

        if pause; readline(); end

        Jtemp = J

        iter += 1
    end

    return θ
end

function compose_theta(x,s,h)
    N = length(x)
    θ = zeros(2N+1)
    for ti=1:N
        θ[ti] = x[ti]
        θ[ti+N] = s[ti]
    end
    θ[2N+1] = h
    return θ
end

function decompose_theta(θ)
    n_temp = length(θ)
    N = Int((n_temp-1)/2)
    x = zeros(N)
    s = zeros(N)
    for ti=1:N
        x[ti] = θ[ti]
        s[ti] = θ[ti+N]
    end
    h = θ[2N+1]
    return x,s,h
end


text:
using CUDA
using LogicCircuits

export LayeredParameterCircuit, CuLayeredParameterCircuit
export class_likelihood, class_weights
export one_hot, learn_parameters, update_parameters

#############################################################
############## This is the old implementation ###############
#### Not intended to be used under the current framework ####
#############################################################

# in a parameter circuit
# 1 is true, 2 is false
const TRUE_ID = Int32(1)
const FALSE_ID = Int32(2)

struct LayeredParameterCircuit
    layered_circuit::LayeredBitCircuit
    layered_parameters::Vector{Matrix{Float32}}
end

LayeredParameterCircuit(circuit::LogisticCircuit, nc::Integer, num_features::Integer) = begin
    @assert is⋁gate(circuit)
    decisions::Vector{Vector{Int32}} = Vector{Vector{Int32}}()
    elements::Vector{Vector{Int32}} = Vector{Vector{Int32}}()
    parameters::Vector{Vector{Float32}} = Vector{Vector{Float32}}()
    num_decisions::Int32 = 2 * num_features + 2
    num_elements::Vector{Int32} = Vector{Int32}()
    # num_parameters always equals num_elements

    ensure_layer(i) = begin
        if length(decisions) < i
            # add a new layer
            push!(decisions, Int32[])
            push!(elements, Int32[])
            push!(parameters, Float32[])
            push!(num_elements, 0)
        end
    end

    f_con(n) = LayeredDecisionId(0, istrue(n) ? TRUE_ID : FALSE_ID)
    f_lit(n) = LayeredDecisionId(0,
        ispositive(n) ? Int32(2 + variable(n)) : Int32(2 + num_features + variable(n)))

    f_and(n, cs) = begin
        @assert length(cs) == 2
        LayeredDecisionId[cs[1], cs[2]]
    end
    f_or(n, cs) = begin
        num_decisions += 1
        # determine layer
        layer_id = zero(Int32)
        for c in cs
            if c isa Vector{LayeredDecisionId}
                @assert length(c) == 2
                layer_id = max(layer_id, c[1].layer_id, c[2].layer_id)
            else
                @assert c isa LayeredDecisionId
                layer_id = max(layer_id, c.layer_id)
            end
        end
        layer_id += 1
        ensure_layer(layer_id)
        first_element = num_elements[layer_id] + 1
        foreach(cs, eachrow(n.thetas)) do c, theta
            @assert size(theta)[1] == nc
            append!(parameters[layer_id], theta)
            num_elements[layer_id] += 1
            if c isa Vector{LayeredDecisionId}
                push!(elements[layer_id], c[1].decision_id, c[2].decision_id)
            else
                push!(elements[layer_id], c.decision_id, TRUE_ID)
            end
        end
        push!(decisions[layer_id], num_decisions, first_element, num_elements[layer_id])
        LayeredDecisionId(layer_id, num_decisions)
    end

    foldup_aggregate(circuit, f_con, f_lit, f_and, f_or,
        Union{LayeredDecisionId,Vector{LayeredDecisionId}})

    circuit_layers = map(decisions, elements) do d, e
        Layer(reshape(d, 3, :), reshape(e, 2, :))
    end
    parameter_layers = map(parameters) do p
        reshape(p, nc, :)
    end
    return LayeredParameterCircuit(LayeredBitCircuit(circuit_layers), parameter_layers)
end

struct CuLayeredParameterCircuit
    layered_circuit::CuLayeredBitCircuit
    layered_parameters::Vector{CuMatrix{Float32}}
    CuLayeredParameterCircuit(l::LayeredParameterCircuit) = new(CuLayeredBitCircuit(l.layered_circuit), map(CuMatrix, l.layered_parameters))
end

function class_likelihood(circuit::CuLayeredParameterCircuit, nc::Integer, data::CuMatrix{Float32}, reuse_up=nothing, reuse_down=nothing, reuse_cp=nothing)
    cw, flow, v = class_weights(circuit, nc, data, reuse_up, reuse_down, reuse_cp)
    one = Float32(1.0)
    return @. one / (one + exp(-cw)), flow, v
end

function class_weights(circuit::CuLayeredParameterCircuit, nc::Integer, data::CuMatrix{Float32}, reuse_up=nothing, reuse_down=nothing, reuse_cw=nothing)
    flow, v = compute_flows2(circuit.layered_circuit, data, reuse_up, reuse_down)
    cw = calculate_class_weights(circuit, nc, data, v, flow, reuse_cw)
    return cw, flow, v
end

function calculate_class_weights(circuit::CuLayeredParameterCircuit, nc::Integer, data::CuMatrix{Float32}, v, flow, reuse_cw=nothing)
    ne = num_examples(data)
    cw = if reuse_cw isa CuMatrix{Float32} && size(reuse_cw) == (ne, nc)
        reuse_cw .= zero(Float32)
        reuse_cw
    else
        CUDA.zeros(Float32, ne, nc)
    end

    dec_per_thread = 4
    CUDA.@sync for i = 1:length(circuit.layered_circuit.layers)
        circuit_layer = circuit.layered_circuit.layers[i]
        parameter_layer = circuit.layered_parameters[i]
        ndl = num_decisions(circuit_layer)
        num_threads = balance_threads(ne, ndl / dec_per_thread, 8)
        num_blocks = ceil(Int, ne / num_threads[1]), ceil(Int, ndl / num_threads[2] / dec_per_thread)
        @cuda threads=num_threads blocks=num_blocks calculate_class_weights_layer_kernel_cuda(cw, v, flow, circuit_layer.decisions, circuit_layer.elements, parameter_layer)
    end

    return cw
end

function calculate_class_weights_layer_kernel_cuda(cw, v, flow, decisions, elements, parameters)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    ne, nc = size(cw)
    _, num_decisions = size(decisions)

    for j = index_x:stride_x:ne
        for i = index_y:stride_y:num_decisions
            decision_id = @inbounds decisions[1, i]
            n_up = @inbounds v[j, decision_id]
            if n_up > zero(Float32)
                first_elem = @inbounds decisions[2, i]
                last_elem = @inbounds decisions[3, i]
                n_down = @inbounds flow[j, decision_id]
                for e = first_elem:last_elem
                    e1 = @inbounds elements[1, first_elem]
                    e2 = @inbounds elements[2, first_elem]
                    e_up = @inbounds (v[j, e1] * v[j, e2])
                    edge_flow = e_up / n_up * n_down
                    # following needs to be memory safe
                    for class=1:nc
                        @CUDA.atomic cw[j, class] += edge_flow * parameters[class, e] # atomic is automatically inbounds
                    end
                end
            end
        end
    end

    return nothing
end

function one_hot(labels::Vector, nc::Integer)
    ne = length(labels)
    one_hot_labels = zeros(Float32, ne, nc)
    for (i, j) in enumerate(labels)
        one_hot_labels[i, j + 1] = 1.0
    end
    one_hot_labels
end

function learn_parameters(circuit::CuLayeredParameterCircuit, nc::Integer, data::CuMatrix{Float32}, labels::CuMatrix{Float32}, reuse_up=nothing, reuse_down=nothing, reuse_cp=nothing, num_epochs=20, step_size=0.0001)
    cp, flow, v = class_likelihood(circuit, nc, data, reuse_up, reuse_down, reuse_cp)
    update_parameters(circuit, labels, cp, flow, step_size)
    for _ = 2:num_epochs
        cp, flow, v = class_likelihood(circuit, nc, data, v, flow, cp)
        update_parameters(circuit, labels, cp, v, flow, step_size)
    end
    return nothing
end

function update_parameters(circuit::CuLayeredParameterCircuit, labels, cp, v, flow, step_size=0.0001)
    _, nc = size(labels)
    step_size = Float32(step_size)
    CUDA.@sync for i = 1:length(circuit.layered_circuit.layers)
        circuit_layer = circuit.layered_circuit.layers[i]
        flow_layer = flow[i]
        parameter_layer = circuit.layered_parameters[i]
        ndl = num_decisions(circuit_layer)
        num_threads = balance_threads(ndl, nc, 6)
        num_threads = num_threads[1], num_threads[2],
        num_blocks = ceil(Int, ndl / num_threads[1]), ceil(Int, nc / num_threads[2]), 4
        @cuda threads=num_threads blocks=num_blocks update_parameters_layer_kernel_cuda(labels, cp, flow_layer, circuit_layer.decisions, parameter_layer, step_size)
    end
    return nothing
end

function update_parameters_layer_kernel_cuda(labels, cp, flow, decisions, parameters, step_size)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    index_z = (blockIdx().z - 1) * blockDim().z + threadIdx().z
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    stride_z = blockDim().z * gridDim().z
    ne, nc = size(labels)
    _, num_decisions = size(decisions)

    for class = index_y:stride_y:nc
        for i = index_x:stride_x:num_decisions
            first_elem = @inbounds decisions[2, i]
            last_elem = @inbounds decisions[3, i]
            for e = first_elem:last_elem
                for j = index_z:stride_z:ne
                    edge_flow = e_up / n_up * n_down
                    u = @inbounds edge_flow * (cp[j, class] - labels[j, class]) * step_size
                    # following needs to be memory safe
                    @inbounds parameters[class, e] -= u
                end
            end
        end
    end

    return nothing
end


text:
export direction, state, value

### GPIO and LED ###

"""
Set the direction of digital IO pin to either input (true)
or output (false)
"""
function direction(rp::RedPitaya, pin::String, input::Bool)
  dir = input ? "INP" : "OUTP"
  send(rp,"DIG:PIN:DIR $(pin),$(dir)")
end

function state(rp::RedPitaya, pin::String)
  send(rp,"DIG:PIN? $(pin)")
  return parse(Int, receive(rp)) == 1
end

function state(rp::RedPitaya, pin::String, value::Bool)
  cmd ="DIG:PIN $(pin),$(Int(value))"
  send(rp,cmd)
end

### Analog IO ###

function value(rp::RedPitaya, pin::String)
  send(rp,"ANALOG:PIN? $(pin)")
  return parse(Float64, receive(rp))
end

function value(rp::RedPitaya, pin::String, value::Float64)
  cmd ="ANALOG:PIN $(pin),$(value)"
  send(rp,cmd)
end


text:
"""
    innerloops!(  iz::Int, rb_tar_idx::Int, rb_src_idx::Int, gw_x::Array{T, 2}, gw_y::Array{T, 2}, gw_z::Array{T, 2}, fssrb::PotentialSimulationSetupRB{T, 3, 4, :cartesian},
                                update_even_points::Val{even_points},
                                depletion_handling::Val{depletion_handling_enabled},
                                bulk_is_ptype::Val{_bulk_is_ptype}  )::Nothing where {T, even_points, depletion_handling_enabled, _bulk_is_ptype}

(Vectorized) inner loop for Cartesian coordinates. This function does all the work in the field calculation.
"""
@fastmath function innerloops!( iz::Int, rb_tar_idx::Int, rb_src_idx::Int, gw_x::Array{T, 2}, gw_y::Array{T, 2}, gw_z::Array{T, 2}, fssrb::PotentialSimulationSetupRB{T, 3, 4, :cartesian},
                                update_even_points::Val{even_points},
                                depletion_handling::Val{depletion_handling_enabled},
                                bulk_is_ptype::Val{_bulk_is_ptype},
                                is_weighting_potential::Val{_is_weighting_potential}  )::Nothing where {T, even_points, depletion_handling_enabled, _bulk_is_ptype, _is_weighting_potential}
    @inbounds begin
        inz::Int = iz - 1

        pwwzr::T        = gw_z[1, inz]
        pwwzl::T        = gw_z[2, inz]
        pwΔmpz::T       = gw_z[3, inz]
        Δz_ext_inv_r::T = gw_z[4, inz + 1]
        Δz_ext_inv_l::T = gw_z[4, inz]

        for iy in 2:(size(fssrb.potential, 2) - 1)
            iny::Int = iy - 1
            zyi_is_even::Bool = iseven(iz + iy)
            zyi_is_even_t = zyi_is_even ? Val{true}() : Val{false}()

            pwwyr::T  = gw_y[1, iny]
            pwwyl::T  = gw_y[2, iny]
            pwΔmpy::T = gw_y[3, iny]
            pwΔmpy_pwΔmpz::T = pwΔmpy * pwΔmpz
            Δy_ext_inv_r_pwΔmpz::T  = gw_y[4, iny + 1] * pwΔmpz
            Δy_ext_inv_l_pwΔmpz::T  = gw_y[4, iny]     * pwΔmpz
            Δz_ext_inv_r_pwΔmpy::T = Δz_ext_inv_r * pwΔmpy
            Δz_ext_inv_l_pwΔmpy::T = Δz_ext_inv_l * pwΔmpy

            pwwyr_pwwzr::T = pwwyr * pwwzr
            pwwyr_pwwzl::T = pwwyr * pwwzl
            pwwyl_pwwzr::T = pwwyl * pwwzr
            pwwyl_pwwzl::T = pwwyl * pwwzl

            # for ix in 2:(size(fssrb.potential, 1) - 1)
            @fastmath @inbounds @simd ivdep for ix in 2:(size(fssrb.potential, 1) - 1)
                inx::Int = nidx(ix, update_even_points, zyi_is_even_t)::Int
                ixr::Int = ifelse( zyi_is_even, ix, even_points ? ix - 1 : ix + 1)
                ixr += ifelse(even_points, 1, 0)
                pwwxr::T        = gw_x[1, inx]
                pwwxl::T        = gw_x[2, inx]
                pwΔmpx::T       = gw_x[3, inx]
                Δx_ext_inv_r::T = gw_x[4, inx + 1]
                Δx_ext_inv_l::T = gw_x[4, inx]

                pwwxr_pwwzr::T = pwwxr * pwwzr
                pwwxl_pwwzr::T = pwwxl * pwwzr
                pwwxr_pwwzl::T = pwwxr * pwwzl
                pwwxl_pwwzl::T = pwwxl * pwwzl
                pwwxl_pwwyr::T = pwwxl * pwwyr
                pwwxr_pwwyr::T = pwwxr * pwwyr
                pwwxl_pwwyl::T = pwwxl * pwwyl
                pwwxr_pwwyl::T = pwwxr * pwwyl

                ϵ_rrr::T = fssrb.ϵ[ inx + 1,  iy, iz]
                ϵ_rlr::T = fssrb.ϵ[ inx + 1, iny, iz]
                ϵ_rrl::T = fssrb.ϵ[ inx + 1,  iy, inz ]
                ϵ_rll::T = fssrb.ϵ[ inx + 1, iny, inz ]
                ϵ_lrr::T = fssrb.ϵ[ inx,  iy,  iz ]
                ϵ_llr::T = fssrb.ϵ[ inx, iny,  iz ]
                ϵ_lrl::T = fssrb.ϵ[ inx,  iy, inz ]
                ϵ_lll::T = fssrb.ϵ[ inx, iny, inz ]

                # right weight in r: wrr
                wxr::T = ϵ_rrr * pwwyr_pwwzr
                wxr    = muladd(ϵ_rlr, pwwyl_pwwzr, wxr)
                wxr    = muladd(ϵ_rrl, pwwyr_pwwzl, wxr)
                wxr    = muladd(ϵ_rll, pwwyl_pwwzl, wxr)
                # left weight in r: wrr
                wxl::T = ϵ_lrr * pwwyr_pwwzr
                wxl    = muladd(ϵ_llr, pwwyl_pwwzr, wxl)
                wxl    = muladd(ϵ_lrl, pwwyr_pwwzl, wxl)
                wxl    = muladd(ϵ_lll, pwwyl_pwwzl, wxl)
                # right weight in φ: wφr
                wyr::T = ϵ_lrr * pwwxl_pwwzr
                wyr    = muladd(ϵ_rrr, pwwxr_pwwzr, wyr)
                wyr    = muladd(ϵ_lrl, pwwxl_pwwzl, wyr)
                wyr    = muladd(ϵ_rrl, pwwxr_pwwzl, wyr)
                # left weight in φ: wφl
                wyl::T = ϵ_llr * pwwxl_pwwzr
                wyl    = muladd(ϵ_rlr, pwwxr_pwwzr, wyl)
                wyl    = muladd(ϵ_lll, pwwxl_pwwzl, wyl)
                wyl    = muladd(ϵ_rll, pwwxr_pwwzl, wyl)
                # right weight in z: wzr
                wzr::T = ϵ_rrr * pwwxr_pwwyr
                wzr    = muladd(ϵ_rlr, pwwxr_pwwyl, wzr)
                wzr    = muladd(ϵ_lrr, pwwxl_pwwyr, wzr)
                wzr    = muladd(ϵ_llr, pwwxl_pwwyl, wzr)
                # left weight in z: wzr
                wzl::T = ϵ_rrl * pwwxr_pwwyr
                wzl    = muladd(ϵ_rll, pwwxr_pwwyl, wzl)
                wzl    = muladd(ϵ_lrl, pwwxl_pwwyr, wzl)
                wzl    = muladd(ϵ_lll, pwwxl_pwwyl, wzl)

                wxr *= pwΔmpy_pwΔmpz * Δx_ext_inv_r
                wxl *= pwΔmpy_pwΔmpz * Δx_ext_inv_l
                wyr *= Δy_ext_inv_r_pwΔmpz * pwΔmpx
                wyl *= Δy_ext_inv_l_pwΔmpz * pwΔmpx
                wzr *= Δz_ext_inv_r_pwΔmpy * pwΔmpx
                wzl *= Δz_ext_inv_l_pwΔmpy * pwΔmpx

                vxr::T = fssrb.potential[    ixr,     iy,     iz, rb_src_idx]
                vxl::T = fssrb.potential[ixr - 1,     iy,     iz, rb_src_idx]
                vyr::T = fssrb.potential[     ix, iy + 1,     iz, rb_src_idx]
                vyl::T = fssrb.potential[     ix,    iny,     iz, rb_src_idx]
                vzr::T = fssrb.potential[     ix,     iy, iz + 1, rb_src_idx]
                vzl::T = fssrb.potential[     ix,     iy,    inz, rb_src_idx]

                new_potential::T = _is_weighting_potential ? 0 : fssrb.ρ[ix, iy, iz, rb_tar_idx]
                new_potential = muladd( wxr, vxr, new_potential)
                new_potential = muladd( wxl, vxl, new_potential)
                new_potential = muladd( wyr, vyr, new_potential)
                new_potential = muladd( wyl, vyl, new_potential)
                new_potential = muladd( wzr, vzr, new_potential)
                new_potential = muladd( wzl, vzl, new_potential)

                new_potential *= fssrb.volume_weights[ix, iy, iz, rb_tar_idx]

                old_potential::T = fssrb.potential[ix, iy, iz, rb_tar_idx]

                new_potential -= old_potential
                new_potential = muladd(new_potential, fssrb.sor_const[1], old_potential)

                if depletion_handling_enabled
                    if _bulk_is_ptype # p-type detectors
                        if new_potential < fssrb.minimum_applied_potential
                            # new_potential = fssrb.minimum_applied_potential
                            new_potential -= fssrb.ρ[ix, iy, iz, rb_tar_idx] * fssrb.volume_weights[ix, iy, iz, rb_tar_idx] * fssrb.sor_const[1]
                            if (fssrb.pointtypes[ix, iy, iz, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[ix, iy, iz, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                        else
                            vmin::T = ifelse( vxr <  vxl, vxr,  vxl)
                            vmin    = ifelse( vyr < vmin, vyr, vmin)
                            vmin    = ifelse( vyl < vmin, vyl, vmin)
                            vmin    = ifelse( vzr < vmin, vzr, vmin)
                            vmin    = ifelse( vzl < vmin, vzl, vmin)
                            if new_potential <= vmin # bubble point
                                new_potential -= fssrb.ρ[ix, iy, iz, rb_tar_idx] * fssrb.volume_weights[ix, iy, iz, rb_tar_idx] * fssrb.sor_const[1]
                                if (fssrb.pointtypes[ix, iy, iz, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[ix, iy, iz, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                            else # normal point
                                if (fssrb.pointtypes[ix, iy, iz, rb_tar_idx] & undepleted_bit > 0) fssrb.pointtypes[ix, iy, iz, rb_tar_idx] -= undepleted_bit end # unmark this point
                            end
                        end
                    else # n-type detectors
                        if new_potential > fssrb.maximum_applied_potential
                            new_potential -= fssrb.ρ[ix, iy, iz, rb_tar_idx] * fssrb.volume_weights[ix, iy, iz, rb_tar_idx] * fssrb.sor_const[1]
                            if (fssrb.pointtypes[ix, iy, iz, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[ix, iy, iz, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                        else
                            vmax::T = ifelse( vxr >  vxl, vxr,  vxl)
                            vmax    = ifelse( vyr > vmax, vyr, vmax)
                            vmax    = ifelse( vyl > vmax, vyl, vmax)
                            vmax    = ifelse( vzr > vmax, vzr, vmax)
                            vmax    = ifelse( vzl > vmax, vzl, vmax)
                            if new_potential >= vmax # bubble point
                                new_potential -= fssrb.ρ[ix, iy, iz, rb_tar_idx] * fssrb.volume_weights[ix, iy, iz, rb_tar_idx] * fssrb.sor_const[1]
                                if (fssrb.pointtypes[ix, iy, iz, rb_tar_idx] & undepleted_bit == 0) fssrb.pointtypes[ix, iy, iz, rb_tar_idx] += undepleted_bit end # mark this point as undepleted
                            else # normal point -> unmark
                                if (fssrb.pointtypes[ix, iy, iz, rb_tar_idx] & undepleted_bit > 0) fssrb.pointtypes[ix, iy, iz, rb_tar_idx] -= undepleted_bit end # unmark this point
                            end
                        end
                    end
                end

                fssrb.potential[ix, iy, iz, rb_tar_idx]::T = ifelse(fssrb.pointtypes[ix, iy, iz, rb_tar_idx] & update_bit > 0, new_potential, old_potential)
            end # x loop
        end # y loop
    end # inbounds
end


text:
Base.isfinite(x::Infinite) = false

Base.isinf(x::Infinite) = true

Base.iszero(::Infinite) = false

Base.:(==)(x::Infinite, y::Infinite) = x.signbit == y.signbit
Base.:(==)(x::Infinite, y::Real) = isinf(y) && signbit(x) == signbit(y)
Base.:(==)(x::Real, y::Infinite) = y == x

Base.hash(x::Infinite, h::UInt) = hash(isposinf(x) ? Inf : -Inf, h)

Base.:<(x::Infinite, y::Infinite) = isneginf(x) && isposinf(y)

Base.:≤(x::Infinite, y::Infinite) = !(y < x)

Base.signbit(x::Infinite) = x.signbit

Base.sign(x::Infinite) = isposinf(x) ? 1 : -1

Base.isapprox(x::Infinite, y::Infinite) = x==y


text:
include("helper.jl")

my_tests = [

  "test_factors.jl",
  "test_fast_exponentiation.jl",
  "test_next_prime.jl",
  "test_rsa.jl",

]

for t in my_tests
  include(t)
end


text:
@testset "$(rpad("Runge-Kutta Tableaus",80))" begin

    # explicit midpoint
    a = Rational{Int64}[[0    0]
                        [1//2 0]]
    b = Rational{Int64}[0, 1   ]
    c = Rational{Int64}[0, 1//2]
    o = 2

    @test TableauERK(:explicit_midpoint, o, a, b, c) == getTableauExplicitMidpoint()

    tab_explicit_midpoint1 = @test_nowarn TableauERK(:explicit_midpoint, o, a, b, c)
    tab_explicit_midpoint2 = @test_logs (:warn, r"Initializing TableauDIRK with explicit tableau explicit_midpoint.*") TableauDIRK(:explicit_midpoint, o, a, b, c)
    tab_explicit_midpoint3 = @test_logs (:warn, r"Initializing TableauFIRK with explicit tableau explicit_midpoint.*") TableauFIRK(:explicit_midpoint, o, a, b, c)

    @test tab_explicit_midpoint1 == tab_explicit_midpoint2
    @test tab_explicit_midpoint2 == tab_explicit_midpoint3
    @test tab_explicit_midpoint3 == tab_explicit_midpoint1

    @test !isequal(tab_explicit_midpoint1, tab_explicit_midpoint2)
    @test !isequal(tab_explicit_midpoint2, tab_explicit_midpoint3)
    @test !isequal(tab_explicit_midpoint3, tab_explicit_midpoint1)

    tab_explicit_midpoint1 = TableauERK(:explicit_midpoint, o, a, b, c)
    tab_explicit_midpoint2 = TableauERK(:explicit_midpoint, o, a, b, c)
    @test tab_explicit_midpoint1 == tab_explicit_midpoint2
    @test isequal(tab_explicit_midpoint1, tab_explicit_midpoint2)

    tmp = mktempdir()
    tab_explicit_midpoint1 = TableauERK(:explicit_midpoint, o, a, b, c)
    writeTableauToFile(tmp, tab_explicit_midpoint1)
    tab_explicit_midpoint2 = readTableauERKFromFile(tmp, "explicit_midpoint")
    rm(tmp, recursive=true)
    @test tab_explicit_midpoint1 == tab_explicit_midpoint2
    @test !isequal(tab_explicit_midpoint1, tab_explicit_midpoint2) # fails because data arrays in tab_explicit_midpoint1 are of type Rational{Int64} while data arrays in tab_explicit_midpoint2 are of type Float

    # Heun
    a = [[0.0 0.0]
         [1.0 0.0]]
    b = [0.5, 0.5]
    c = [0.0, 1.0]
    o = 2

    @test_nowarn tab_explicit_heun = TableauERK(:heun, o, a, b, c)
    @test_logs (:warn, r"Initializing TableauDIRK with explicit tableau heun.*") TableauDIRK(:heun, o, a, b, c)
    @test_logs (:warn, r"Initializing TableauFIRK with explicit tableau heun.*") TableauFIRK(:heun, o, a, b, c)

    @test TableauERK(:heun, o, a, b, c) == getTableauHeun()

    tmp = mktempdir()
    tab_explicit_heun1 = TableauERK(:heun, o, a, b, c)
    writeTableauToFile(tmp, tab_explicit_heun1)
    tab_explicit_heun2 = readTableauERKFromFile(tmp, "heun")
    rm(tmp, recursive=true)
    @test tab_explicit_heun1 == tab_explicit_heun2
    @test isequal(tab_explicit_heun1, tab_explicit_heun2)

    # Crouzeix
    fac = 0.5/√3
    a = [[ 0.5+fac 0.0    ]
         [-2.0*fac 0.5+fac]]
    b = [0.5,     0.5    ]
    c = [0.5+fac, 0.5-fac]
    o = 2

    @test TableauDIRK(:crouzeix, o, a, b, c) == getTableauCrouzeix()

    @test_throws AssertionError tab_explicit_crouzeix = TableauERK(:crouzeix, o, a, b, c)
    @test_nowarn tab_explicit_crouzeix = TableauDIRK(:crouzeix, o, a, b, c)
    @test_logs (:warn, r"Initializing TableauFIRK with diagonally implicit tableau crouzeix.*") tab_explicit_crouzeix = TableauFIRK(:crouzeix, o, a, b, c)

    # implicit midpoint
    a = ones(Rational{Int64}, 1, 1) // 2
    b = Rational{Int64}[1   ]
    c = Rational{Int64}[1//2]
    o = 2

    @test TableauFIRK(:implicit_midpoint, o, a, b, c) == getTableauImplicitMidpoint()
    @test TableauFIRK(:implicit_midpoint, o, a, b, c) == getTableauGLRK(1)

    @test_throws AssertionError tab_implicit_midpoint = TableauERK(:implicit_midpoint, o, a, b, c)
    @test_throws AssertionError tab_implicit_midpoint = TableauDIRK(:implicit_midpoint, o, a, b, c)
    tab_implicit_midpoint = TableauFIRK(:implicit_midpoint, o, a, b, c)

    # Gauss-Legendre
    a = Array{Float64}(@dec128 [
         [0.25       0.25-√3/6]
         [0.25+√3/6  0.25     ]
        ])
    b = Array{Float64}(@dec128 [0.5,      0.5     ])
    c = Array{Float64}(@dec128 [0.5-√3/6, 0.5+√3/6])
    o = 4

    @test TableauFIRK(:glrk2, o, a, b, c) == getTableauGLRK(2)

    @test_throws AssertionError tab_explicit_glrk2 = TableauERK(:glrk2, o, a, b, c)
    @test_throws AssertionError tab_explicit_glrk2 = TableauDIRK(:glrk2, o, a, b, c)
    tab_explicit_glrk2 = TableauFIRK(:glrk2, o, a, b, c)

    # instantiate all explicit tableaus
    @test typeof(getTableauExplicitEuler()) <: TableauERK
    @test typeof(getTableauExplicitMidpoint()) <: TableauERK
    @test typeof(getTableauHeun()) <: TableauERK
    @test typeof(getTableauKutta()) <: TableauERK
    @test typeof(getTableauERK4()) <: TableauERK
    @test typeof(getTableauERK438()) <: TableauERK

    # instantiate all diagonally implicit tableaus
    @test typeof(getTableauCrouzeix()) <: TableauDIRK

    # instantiate all fully implicit tableaus
    @test typeof(getTableauImplicitEuler()) <: TableauFIRK
    @test typeof(getTableauImplicitMidpoint()) <: TableauFIRK
    @test typeof(getTableauGLRK(1)) <: TableauFIRK
    @test typeof(getTableauGLRK(2)) <: TableauFIRK

    # instatiate all partitioned tableaus
    @test typeof(getTableauSymplecticEulerA()) <: TableauEPRK
    @test typeof(getTableauSymplecticEulerB()) <: TableauEPRK

    # test instatiation of partioned tableau by composition of two RK tableaus
    @test typeof(TableauEPRK(:PERK4, 4, getTableauERK4().q, getTableauERK4().q)) <: TableauEPRK
    @test TableauEPRK(:PERK4, 4, getTableauERK4().q, getTableauERK4().q) == TableauEPRK(:PERK4, 4, getTableauERK4().q)

    # TODO Add tests for TableauIPRK, TableauSARK, TableauSPARK and TableauGLM.

end


text:
function _h5concat2d(dst, srcs_list, pv...; dims, virtual = true, ka...)
    isfile(dst) && rm(dst)
    @eval GC.gc(true)
    h5open(dst, "w", pv...; ka...) do fid
        layout_map, pos_map = Dict(), Dict()
        prog = Progress(length(srcs_list), desc = "h5concat.config: ")
        verbose = parse(Int, get(ENV, "HDF5_PROGRESS", "1"))
        for (i, srcs) in enumerate(srcs_list)
            for (j, src) in enumerate(srcs)
                fidn = h5open(src, "r", pv...; ka...)
                for c in names(fidn)
                    if isa(fidn[c], HDF5Dataset)
                        vs = VirtualSource(fidn[c])
                        dims′ = [ifelse(d < 0, ndims(fidn[c]) + d + 1, d) for d in dims]
                        if !haskey(layout_map, c)
                            layout_map[c] = VirtualLayout(fill(0, ndims(fidn[c])), eltype(fidn[c]))
                            pos_map[c] = [0, 0]
                        end
                        layout, pos = layout_map[c], pos_map[c]
                        if j == 1
                            pos[1] = dims′[1] == 0 ? 0 : layout.shape[dims′[1]]
                            pos[2] = 0
                        end
                        layout.dtype = promote_type(eltype(fidn[c]), layout.dtype)
                        starts = ntuple(ndims(fidn[c])) do d
                            d == dims′[1] ? pos[1] : d == dims′[2] ? pos[2] : 0
                        end
                        layout[[(1:vs.shape[i]) .+ starts[i] for i in 1:ndims(fidn[c])]...] = vs
                        pos[2] += vs.shape[dims′[2]]
                    elseif !has(fid, c)
                        o_copy(fidn[c], fid, c)
                    end
                end
                close(fidn)
            end
        end
        for c in keys(layout_map)
            d_create_virtual(fid, c, layout_map[c], virtual)
        end
        verbose == 1 && next!(prog)
    end
    return dst
end

function h5concat(dst, srcs, pv...; dims, virtual = true, ka...)
    if length(dims) == 1
        srcs, dims = [srcs], (0, dims)
    end
    _h5concat2d(dst, srcs, pv...; dims = dims, virtual = virtual, ka...)
end

h5concat(srcs, pv...; ka...) = h5concat(randstring() * ".h5", srcs, pv...; ka...)

function h5concat!(h5, name, pattern, pv...; dims, virtual = true, ka...)
    @eval GC.gc(true)
    h5open(h5, "r+", pv...; fclose_degree = H5F_CLOSE_DEFAULT, ka...) do fid
        layout = nothing
        for c in names(fid)
            !occursin(pattern, c) && continue
            !isa(fid[c], HDF5Dataset) && continue
            vs = VirtualSource(fid[c])
            if isnothing(layout)
                layout = VirtualLayout(fill(0, ndims(fid[c])), eltype(fid[c]))
            end
            dims′ = dims < 0 ? ndims(fid[c]) + dims + 1 : dims
            layout.dtype = promote_type(eltype(fid[c]), layout.dtype)
            starts = ntuple(ndims(fid[c])) do d
                d == dims′ ? layout.shape[dims′] : 0
            end
            layout[[(1:vs.shape[i]) .+ starts[i] for i in 1:ndims(fid[c])]...] = vs
        end
        !isnothing(layout) && d_create_virtual(fid, name, layout, virtual)
    end
    return h5
end


text:
module E2ga

import Base.show, Base.+, Base.-, Base.*

export Geometric2

struct Geometric2{T <: Real}
    α::T
    x::T
    y::T
    β::T
end
Geometric2(α::Real, x::Real, y::Real, β::Real) = Geometric2(promote(α, x, y, β)...)
Geometric2(α::Real) = Geometric2(α, zero(α), zero(α), zero(α))

promote_rule(::Type{Geometric2{T}}, ::Type{S}) where {T <: Real,S <: Real} = Geometric2{promote_type(T, S)}
promote_rule(::Type{Geometric2{T}}, ::Type{Geometric2{S}}) where {T <: Real,S <: Real} = Geometric2{promote_type(T, S)}

function +(a::Geometric2, b::Geometric2)
    α = a.α + b.α
    x = a.x + b.x
    y = a.y + b.y
    β = a.β + b.β
    return Geometric2(α, x, y, β)
end

function -(a::Geometric2, b::Geometric2)
    α = a.α - b.α
    x = a.x - b.x
    y = a.y - b.y
    β = a.β - b.β
    return Geometric2(α, x, y, β)
end

function *(a::Geometric2, b::Geometric2)
    α = a.α * b.α + a.x * b.x + a.y * b.y - a.β * b.β
    x = a.α * b.x + a.x * b.α - a.y * b.β + a.β * b.y
    y = a.α * b.y + a.x * b.β + a.y * b.α - a.β * b.x
    β = a.α * b.β + a.x * b.y - a.y * b.x + a.β * b.α
    return Geometric2(α, x, y, β)
end

function show(io::IO, m::Geometric2)
    print(io, "$(m.α) + $(m.x)e1 + $(m.y)e2 + $(m.β)I")
end

end


text:
# function properdivisors(n)
# 	d = Int[]
# 	for i in 1:div(n,2)+1
# 		if n % i == 0
# 			push!(d,i)
# 		end
# 	end
# 	return d
# end

function properdivisors(n)
	d = [1]
	for i in 2:isqrt(n)+1
		if n % i == 0
			push!(d,i)
			other = div(n,i)
			if i != other
				push!(d,other)
			end
		end
	end
	return d
end

function main()
	N = 10000

	s = 0
	for a in 1:N
		da = sum(properdivisors(a))
		b = da
		db = sum(properdivisors(b))
		if a == db && b == da && a != b
			println("$a $b")
			s += a
		end
	end
	println(s)
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
# julia> include("/Users/ian/jd/10-19 Code/12 Julia/12.03 PE/pe1.jl")

println("ans = $(sum([i for i=1:999 if (i%3 == 0) || (i%5 == 0)]))")

max_stars_repo_path: learning exercises/pe/pe1.jl

max_stars_repo_name: ian-double-u/NNfromScratch

max_issues_repo_path: learning exercises/pe/pe1.jl

max_issues_repo_name: ian-double-u/NNfromScratch

max_forks_repo_path: learning exercises/pe/pe1.jl

max_forks_repo_name: ian-double-u/NNfromScratch


text:
@testset "Missing" begin
    @testset "Table equality with missing" begin
        @test ismissing(table([1, 2, missing]) == table([1, 2, missing]))
        @test isequal(table([1,2,missing]), table([1,2,missing]))
        @test ismissing(ndsparse([1], [missing]) == ndsparse([1], [missing]))
        @test isequal(ndsparse([1], [missing]), ndsparse([1], [missing]))
        @test !isequal(ndsparse([2], [missing]), ndsparse([1], [missing]))
    end
    @testset "stack/unstack" begin
        t = table(1:4, [1, missing, 9, 16], [1, 8, 27, missing], names = [:x, :x2, :x3], pkey = :x)
        @test isequal(t, unstack(stack(t)))
    end
end

@testset "dropmissing" begin
    # Missing
    a = table([[rand(Bool) ? missing : rand() for i in 1:30] for i in 1:3]...)

    a2 = dropmissing(a)
    @test all(!ismissing, a2)
    @test all(x -> eltype(x) == Float64, columns(a2))

    a3 = dropmissing(a, 1)
    @test all(!ismissing, select(a3, 1))
    @test eltype(columns(a3)[1]) == Float64
    @test all(x -> Missing <: eltype(x), columns(a3)[2:end])

    # DataValue
    b = table([DataValueArray(rand(30), rand(Bool, 30)) for i in 1:3]...)

    b2 = dropmissing(b)
    @test all(!isna, b2)
    @test all(x -> eltype(x) == Float64, columns(b2))
end
@testset "convertmissing" begin
    a = table([[rand(Bool) ? missing : rand() for i in 1:100] for i in 1:3]..., pkey=1)
    b = ndsparse(a)

    cm = IndexedTables.convertmissing
    @test isequal(a, cm(cm(a, DataValue), Missing))
    @test isequal(b, cm(cm(b, DataValue), Missing))

    c = table([1, missing, 3], DataValueArray(1:3, [true, false, false]))
    @test isequal(cm(c, DataValue), table([1,NA,3],[NA,2,3]))
    @test isequal(cm(c, Missing), table([1,missing,3], [missing,2,3]))
end


text:
using DataFrames

n=10_000_000

# Right now things only work columns of type Array, so
# we need this slighlty cumbersome DataFrame construction
# to prevent DataArray or NullableArray creation
# We are also skipping all Strings because of #14955 (I think)
data_friends = fill(4,n)
data_age = fill(38.2,n)
data_children = fill(2,n)

columns = []
push!(columns, data_friends)
push!(columns, data_age)
push!(columns, data_children)

df = DataFrame(columns, [:friends, :age, :children])