text:
# finds the longest number of observations in arguments like
# x::Svector{2}, y::Vector{SVector{2}}, z::SVector{3}
# (in this case, length of y)
# slow and fails shamelessly if only SVectors are supplied
function _getsize(arrays::Union{AbstractArray,<:SVector}...)
# function _getsize(arrays::Union{AbstractVector,Adjoint,<:SVector}...)
    vectors_length = Int[length(a) for a in arrays if ~isa(a, SVector)]
    @assert length(vectors_length)>0
    maximum(vectors_length)
end

@inline _getobs(x::AbstractArray, i::Int) = x[i]
# @inline _getobs(x::Union{AbstractVector,Adjoint}, i::Int) = x[i]
@inline _getobs(x::SVector, i::Int) = x

# constructs expression SVector(a,b,c) from [:a,:b,:c]
_sym_sarray(v::AbstractVector{Symbol}) = Expr(:call,:SVector, v...)
_sym_sarray(v::AbstractMatrix{Symbol}) = Expr(:call, :(SMatrix{$(size(v)...)}), v...)

list_syms(eq::Expr) = list_symbols(eq).parameters
list_syms(eq::Symbol) = [eq]
list_syms(eq::Number) = Symbol[]

diff_symbol(k::Symbol, j::Symbol) = Symbol("∂", k, "_∂", j)

function add_derivatives(all_eqs, gradient_components)

    new_eqs = OrderedDict{Symbol, Union{Expr, Symbol, Int64, Float64}}()
    for g in gradient_components
        k = diff_symbol(g,g)
        new_eqs[k] = 1
    end

    for (sym, eq) in all_eqs

        new_eqs[sym] = eq
        # symbols the equation depends on
        deps = Dolang.list_parameters(eq)

        for θ in gradient_components

            deq_terms = []

            for a in deps

                # compute dy/da * da/dθ

                if a == θ
                    push!(deq_terms, ( deriv(eq, a), 1) )
                elseif !(a in keys(new_eqs))  # I should keep a list of defined variables
                    # nothing to do
                else
                    da_dθ = diff_symbol(a, θ)
                    if da_dθ in keys(new_eqs) # otherwise it is 0
                        dy_da = diff_symbol(sym, a)
                        if !(dy_da in keys(new_eqs))
                            new_eqs[dy_da] = deriv(eq, a)
                        end
                        push!(deq_terms, (dy_da, da_dθ) )
                    end
                end
            end
            if length(deq_terms)>0
                if length(deq_terms)==1
                    rhs = Expr(:call, :*, deq_terms[1]...)
                else
                    rhs = Expr(:call, :+, [Expr(:call, :*, a, b) for (a,b) in deq_terms]...)
                end
                new_eqs[diff_symbol(sym, θ)] = rhs
            end
        end
    end

    return new_eqs

end

# function add_derivatives(dd::OrderedDict, jac_args::Vector{Symbol})
#     # given a list of equations:
#     # x=p+a
#     # y=p+x
#     # and a list of symbols to differetiate with
#     # [:a, :p]
#     # produces a new list of equations with derivatives:
#     # x=p+a
#     # ∂x_∂a = 1
#     # ∂x_∂p = 1
#     # y=p+x
#     # ∂y_∂a = ..
#     # ∂y_∂p = ...

#     diff_eqs = OrderedDict{Symbol, SymExpr}()
#     for (var, eq) in dd
#         diff_eqs[var] = eq
#         deps = list_syms(eq) # list of variables eq depends on

#         for k in [kk for kk in deps if !(kk in jac_args)]
#             for l in jac_args
#                 dv = diff_symbol(var, k)
#                 ddv = diff_symbol(k, l)
#                 cdv = diff_symbol(var, l)
#                 if haskey(diff_eqs, ddv)
#                     if !haskey(diff_eqs, dv)
#                         deq = Dolang.deriv(eq, k)
#                         diff_eqs[dv] = deq
#                     end
#                     if cdv in keys(diff_eqs)
#                         diff_eqs[cdv] = :($(diff_eqs[cdv])+$dv*$ddv)
#                     else
#                         diff_eqs[cdv] = :($dv*$ddv)
#                     end
#                 end
#             end
#         end
#         for k in [kk for kk in deps if kk in jac_args]
#             dv = diff_symbol(var, k)
#             deq = Dolang.deriv(eq, k)
#             diff_eqs[dv] = deq
#         end

#     end

#     return diff_eqs

# end

function clean_unused(code)

    lines = code.args[2].args

    uses = Vector{Symbol}[]
    defs = Symbol[]

    for l in lines
        if l.head==:(=)
        rhs = l.args[2]
        tt = []
        MacroTools.postwalk(x-> ( ((x isa(Symbol)) && string(x)[end]=='_' ) ? push!(tt, x) : nothing ) , rhs)
        push!(uses, tt)
        push!(defs, l.args[1])
        end

    end

    newlines = [:(return res_)]
    useful = [ :(res_) ]
    for j=1:length(lines)-1
        i = length(lines) - j
        if occursin("oo", string(defs[i])) || (defs[i] ∈ useful)
            push!(newlines, lines[i])
            append!(useful, uses[i])
        end
    end

    newcode = deepcopy(code)
    reverse!(newlines)
    newcode.args[2].args = newlines
    return newcode

end

"""
Create a non allocating kernel from the function factory.

`fff`: assumed to be a `FlatFunctionFactory` object with empty preamble.
`diff`: index of variables to differentiate with or list of indices of variables positions.

The generated kernel looks like (diff=[0, 1])
```
function myfun(x::SVector{1, Float64}, y::SVector{3, Float64}, z::SVector{2, Float64}, p::SVector{1, Float64})
    _a_m1_ = x[1]
    _a__0_ = y[1]
    _b__0_ = y[2]
    _c__0_ = y[3]
    _c__1_ = z[1]
    _d__1_ = z[2]
    _u_ = p[1]
    _foo__0_ = log(_a__0_) + _b__0_ / (_a_m1_ / (1 - _c__0_))
    _bar__0_ = _c__1_ + _u_ * _d__1_
    d__foo__0__d__a_m1_ = (-(1 / (1 - _c__0_)) * _b__0_) / (_a_m1_ / (1 - _c__0_)) ^ 2
    d__bar__0__d__a_m1_ = 0
    oo_0_ = SVector(_foo__0_, _bar__0_)
    oo_1_ = SMatrix{2, 1}(d__foo__0__d__a_m1_, d__bar__0__d__a_m1_)
    res_ = (oo_0_, oo_1_)
    return res_
end
```

If diff is a scalar, the result of the kernel is a static vector (or a static matrix).
If diff is a list, the result is a tuple.
"""
function gen_kernel(fff::FlatFunctionFactory, diff::Vector{Int}; funname=fff.funname, arguments=fff.arguments)

    targets = [keys(fff.equations)...]
    equations = [values(fff.equations)...]
    # names of symbols to output
    output_names = []
    for d in diff
        if d == 0
            push!(output_names, targets)
        else
            diff_args = collect(values(fff.arguments))[d]
            p = length(targets)
            q = length(diff_args)
            mat = Matrix{Symbol}(undef, p, q)
            for i in 1:p
                for j in 1:q
                    mat[i, j] = diff_symbol(targets[i], diff_args[j])
                end
            end
            push!(output_names, mat)
        end
    end

    argnames = collect(keys(arguments))

    all_eqs = cat(values(fff.preamble)..., equations, dims=1)
    all_args = cat(values(fff.arguments)..., dims=1)
    if maximum(diff)>0
        jac_args = cat([collect(values(fff.arguments))[i] for i in diff if i!=0]..., dims=1)
    else
        jac_args = []
    end
    # jac_args = Symbol.(jac_args) # strange type of output can by Any[]
    jac_args = Symbol[Symbol(e) for e in jac_args]

    # concatenate preamble and equations (doesn't make much sense...)
    dd = OrderedDict()
    for (k, v) in (fff.preamble)
        dd[k] = v
    end
    for (target, eq) in zip(targets, equations)
        dd[target] = eq
    end
    # compute all equations to write
    diff_eqs = add_derivatives(dd, jac_args)
    for out in output_names
        for k in out
            if !(haskey(diff_eqs, k))
                diff_eqs[k] = 0.0
            end
        end
    end

    diff_eqs = reorder_triangular_block(diff_eqs)

    # create function block
    code = []
    for (k, args) in zip(argnames, values(arguments))
        for (i, a) in enumerate(args)
            push!(code, :($a = ($k)[$i]))
        end
    end
    for (k, v) in diff_eqs
        push!(code, :($k=$v))
    end

    return_args = []
    for (d, names) in enumerate(output_names)
        outname = Symbol("oo_", d, "_")
        push!(code, :($outname = $(_sym_sarray(names))))
        push!(return_args, outname)
    end

    # this is to make inserting the resulting code in another function easier
    push!(code, :(res_ = $(Expr(:tuple, return_args...))))
    push!(code, :(return res_))

    # now we construct the function
    typed_args = [:($k::SVector{$(length(v)), Float64}) for (k, v) in arguments]
    fun_args = Expr(:call, funname, typed_args...)

    ncode = Expr(:function, fun_args, Expr(:block, code...))

    return clean_unused(ncode)

end

# NOTE: for small arrays (probably <= 10 elements) the splatting below is
#       faster than `reinterpret(SVector{length(v), Float64}, v)`
to_SA(v::AbstractVector{Float64}) = SVector(v...)
to_SA(v::AbstractMatrix{Float64}) = copy(reshape(reinterpret(SVector{size(v, 2),Float64}, vec(copy(v'))), (size(v, 1),)))

function from_SA(v::AbstractVector{SMatrix{p,q,Float64,k}}) where p where q where k
    copy(permutedims(reshape(reinterpret(Float64, vec(v)), (p, q, size(v, 1))), [3, 1, 2]))
end
from_SA(v::AbstractVector{SVector{d,Float64}}) where d  = copy(reshape(reinterpret(Float64, vec(v)), (d, size(v, 1)))')
from_SA(v::SVector{d,Float64}) where d = Vector(v)

# that one is a bit risky as from_SA(to_SA(mat)) != mat
from_SA(v::SMatrix{p,q,Float64,k}) where p where q where k = Matrix(v)
from_SA(v::Tuple) = tuple([from_SA(e) for e in v]...)

"""
Creates a vectorized function, which can accept points or list-of-points as
arguments.

`fff`: assumed to be a `FlatFunctionFactory` object with empty preamble.
`diff`: index of variables to differentiate with or list of indices of variables positions.
`out`: optional preallocated output

If at least one of the arguments is a list of points, the result is a list of
points (or a tuple mad of lists of points). In this case preallocated
structures can be passed as `out`.
"""
function gen_gufun(fff::FlatFunctionFactory, to_diff::Union{Array{Int}, Int};
    funname=fff.funname)

    if to_diff isa Int
        diff = [to_diff]::Array{Int}
    else
        diff = to_diff::Array{Int}
    end

    # Generate the code for the kernel
    # change arguments in the kernel to avoid clashes
    arguments_ = OrderedDict{Symbol,Vector{Symbol}}(
        (Symbol(k, "_"), v) for (k, v) in fff.arguments
    )
    kernel_code = gen_kernel(fff, diff; funname=:kernel, arguments=arguments_)
    # remove function def and return statement
    kernel_code_stripped = kernel_code.args[2].args[1:end-1]

    args = [keys(fff.arguments)...]
    args_scalar = [keys(arguments_)...]
    args_out = [Symbol("out_", i) for i in 1:length(diff)]
    args_length = [length(v) for v in values(fff.arguments)]

    targets = [keys(fff.equations)...]

    out_types = []
    p = length(targets)
    for d in diff
        if d == 0
            push!(out_types, :(SVector{$p,Float64}))
        else
            q = args_length[d]
            push!(out_types, :(SMatrix{$p,$q,Float64,$(p*q)}))
        end
    end

    code = quote
        function $funname($(args...),out=nothing)

            @fastmath begin
            # @inline $kernel_code # unpure...
            # revisit in the future

            ### if all arguments are 1d.vectors we do vector things
            # if they are static
            if (($(args...)),) isa Tuple{$([:(SVector) for i=1:length(args)]...)}
                # ret = kernel($(args...))
                # return $(to_diff isa Int? :(ret[1]) :  :(ret) )
                $([:($_a = $a) for (_a, a) in zip(args_scalar, args)]...)
                $(kernel_code_stripped...)
                return $(to_diff isa Int ? :(res_[1]) :  :(res_) )
            end

            # if all arguments are array vectors
            if (($(args...)),) isa Tuple{$([:(AbstractVector{Float64}) for i=1:length(args)]...)}
                 # oo = kernel( $( [ :(SVector( $(e)...)) for e in args]...) )
                 $([:($_a = $a) for (_a, a) in zip(args_scalar, args)]...)
                 $(kernel_code_stripped...)
                 ret = ( $( [:(Array(res_[$i])) for i=1:length(args_out)]...),)
                 return $(to_diff isa Int ? :(ret[1]) :  :(ret) )
             end
            # if arguments are array vectors and one of them is not a vector (plausibly a matrix then...)
            hackish = ( (($(args...)),) isa Tuple{$([:(AbstractArray{Float64}) for i=1:length(args)]...)} )
            if hackish
                $([:($a = Dolang.to_SA($a)) for a in args]...)
            end

            N = Dolang._getsize($(args...))::Int

            if isa(out, Nothing)
                $([:($a=zeros($t, N)) for (a, t) in zip(args_out, out_types)]...)
            else
                $([:($a=out[$i]::Vector{$t}) for (i,(a, t)) in enumerate(zip(args_out, out_types))]...)
            end

            @inbounds @simd for n=1:N
                $([:($a_ = Dolang._getobs($a, n)) for (a_, a) in zip(args_scalar, args)]...)
                # res_ = kernel($(args_scalar...))
                $(kernel_code_stripped...)
                $([:($a[n] = res_[$i]) for (i, a) in enumerate(args_out)]...)
            end

            ret = ($(args_out...),)

            if hackish
                ret = Dolang.from_SA(ret)
            end

            return $(to_diff isa Int ? :(ret[1]) :  :(ret) )
        end

        end
    end

    return code

end

# We now produce the code for generated functions.
# The generated functions have a signature like:
# genfun((Val(1),Val(2)), x, y, z, p, out=)
# where the first argument denotes the arguments to differentiate with (0, is zero-order diff)

_get_nums(::Union{Val{n}, Type{Val{n}}, Type{Type{Val{n}}}}) where n = n
_get_nums(t::Type{<:Tuple}) = [_get_nums(i) for i in getfield(t, 3)]
_get_oorders(x::Array{Int,0}) = x[1]
_get_oorders(x::Union{Tuple,<:Array{Int}}) = x

function gen_generated_kernel(fff::FlatFunctionFactory)

    funname = fff.funname
    meta_code = quote
        @generated function $funname(orders, x, y, z, p, out=nothing)
            fff = $(fff) # this is amazing !
            oorders = Dolang._get_nums(orders) # convert into tuples
            code = Dolang.gen_kernel(fff, oorders)
            code.args[2]
        end
    end
end

function gen_generated_gufun(fff::FlatFunctionFactory; funname=fff.funname, dispatch=Nothing)

    args = collect(keys(fff.arguments))
    dispatch_argtype = Meta.parse(string(dispatch))
    dispatch_arg = dispatch == Nothing ? [] : [:(::$(dispatch_argtype))]
    meta_code = quote
        # basic fun for compat
        @generated function $funname($(dispatch_arg...), $(args...), out=nothing)
            fff = $(fff) # this is amazing !
            code = Dolang.gen_gufun(fff, 0)
            code.args[2].args[2]
        end
        @generated function $funname($(dispatch_arg...), orders::Tuple, $(args...), out=nothing)
            fff = $(fff) # this is amazing !
            oorders = Dolang._get_nums(orders) # convert into tuples
            code = Dolang.gen_gufun(fff, oorders)
            code.args[2].args[2]
        end

        @generated function $funname($(dispatch_arg...), orders::Union{Val,Type{<:Val}}, $(args...), out=nothing)
            fff = $(fff) # this is amazing !
            oorders = _get_oorders(collect(Dolang._get_nums(orders)))
            code = Dolang.gen_gufun(fff, oorders)
            code.args[2].args[2]
        end
    end
end

max_issues_repo_name: EconForge/Dolang


text:
# feedforward
export forward!

function forward!(v::voidLayer)
end
function forward!(d::dataLayer)
end
function forward!(h::hiddenLayer)
  rᵢ=h.bottom[1].rₒ
  for i=2:length(h.bottom)
  rᵢ=[rᵢ;h.bottom[i].rₒ]
  end
  h.rᵢ=rᵢ

  h.I=h.W*h.rᵢ.+h.b
  h.rₒ=h.A(h.I)
end

function forward!(p::predictLayer)
  rᵢ=p.bottom[1].rₒ
  for i=2:length(p.bottom)
  rᵢ=[rᵢ;p.bottom[i].rₒ]
  end
  p.rᵢ=rᵢ

  p.rₒ=p.Predict(p.rᵢ)
end
function forward!(l::lossLayer)
  rᵢ=l.bottom[1].rₒ
  for i=2:length(l.bottom)
  rᵢ=[rᵢ;l.bottom[i].rₒ]
  end
  l.rᵢ=rᵢ

  rₒ=l.top[1].rᵢ
  for i=2:length(l.bottom)
  rₒ=[rₒ;l.top[i].rᵢ]
  end
  l.rₒ=rₒ
  l.loss=l.Cost(l.rₒ, l.rᵢ)

end

function forward!(net::Array{Layer})
  for idx in order
    forward!(net[idx])
  end
end

ext: jl


text:
#=
# Create a recursive ordering
# Each state maps to a configuration of states with a rule and ordering
#
# For example, 2D Hilbert has 4 states: n, ], [, u labeled here as s1,s2,s3,s4
#
#       s1 s1
# s1 -> s2 s3  This state, s1, corresponds to rule [s2,s3,s1,s1] and ordering [1,3,4,2]
#
# Morton is the simplest with one rule [s1,s2,s3,s4] and ordering [1,2,3,4]
=#
export get_recursive_order, reorder_grid_recursive!
export HILBERT_ORDERING, MORTON_ORDERING

struct RecursiveOrdering
    states::Int                 # Number of possible states
    rules::Array{Array{Int,1}}  # Rules for each state mapping sectors to states
    orders::Array{Array{Int,1}} # Ordering for each state
end

HILBERT_ORDERING = "hilbert ordering";
MORTON_ORDERING = "morton ordering";

# Use this function to build an ordering for a specified grid size
function get_recursive_order(type, dim, griddim)
    nnodes = 1;
    for i=1:length(griddim)
        nnodes = nnodes*griddim[i];
    end

    if dim == 1
        # This doesn't make sense, to just keep regular ordering
        return 1:nnodes;

    elseif dim == 2
        # Find the smallest N such that 2^N >= max(griddim)
        maxgriddim = max(griddim[1],griddim[2]);
        N=0;
        while maxgriddim > (2^N)
            N = N+1;
        end
        twoN = 2^N;

        if type == HILBERT_ORDERING
            rorder_rules = hilbert_ordering_rules(2); # The RecursiveOrdering object
            rorder = build_ordering_2d(rorder_rules, N, false); # The ordering (false -> inverse)
        elseif type == MORTON_ORDERING
            rorder_rules = morton_ordering_rules(2); # The RecursiveOrdering object
            rorder = build_ordering_2d(rorder_rules, N, false); # The ordering (false -> inverse)
        end

        # Truncate ordering to match griddim size
        subrorder = zeros(Int,nnodes);
        subind = 1;
        for i=1:length(rorder)
            # Get coordinates in 2^N grid
            hind = rorder[i];
            hx = mod(hind-1,twoN) + 1;
            hy = Int(floor((hind-1) / (twoN))) + 1;

            if hx <= griddim[1] && hy <= griddim[2]
                gind = hx + griddim[1]*(hy-1);
                subrorder[subind] = gind;
                subind += 1;
            end
        end
        return subrorder;

    elseif dim == 3
        # Find the smallest N such that 2^N >= max(griddim)
        maxgriddim = max(griddim[1],max(griddim[2],griddim[3]));
        N=0;
        while maxgriddim > (2^N)
            N = N+1;
        end
        twoN = 2^N;

        if type == HILBERT_ORDERING
            rorder_rules = hilbert_ordering_rules(3); # The RecursiveOrdering object
            rorder = build_ordering_3d(rorder_rules, N, false); # The ordering (false -> inverse)
        elseif type == MORTON_ORDERING
            rorder_rules = morton_ordering_rules(3); # The RecursiveOrdering object
            rorder = build_ordering_3d(rorder_rules, N, false); # The ordering (false -> inverse)
        end

        # Truncate ordering to match griddim size
        subrorder = zeros(Int,nnodes);
        subind = 1;
        for i=1:length(rorder)
            # Get coordinates in 2^N grid
            hind = rorder[i];
            hx = mod(hind-1,twoN) + 1;
            hy = Int(floor(mod(hind-1,twoN*twoN) / (twoN))) + 1;
            hz = Int(floor((hind-1) / (twoN*twoN))) + 1;

            if hx <= griddim[1] && hy <= griddim[2] && hz <= griddim[3]
                gind = hx + griddim[1]*((hy-1) + griddim[2]*(hz-1));
                subrorder[subind] = gind;
                subind += 1;
            end
        end
        return subrorder;
    end
end

# Use this to reorder the nodes in a given grid.
function reorder_grid_recursive!(grid, griddim, type)
    dim = size(grid.allnodes,1);

    # get the ordering
    rordering = get_recursive_order(type, dim, griddim);
    # invert it for transfering the grid
    rordering = invert_ordering(rordering);

    reorder_grid_nodes!(grid, rordering);
    return grid;
end

#####################################################################################
# ordering rule sets

# Build a Morton Ordering
function morton_ordering_rules(dim)
    if dim == 1
        return RecursiveOrdering(1,[[1]],[[1]]);
    elseif dim == 2
        return RecursiveOrdering(1,[[1,1,1,1]], [[1,2,3,4]]);
    elseif dim == 3
        return RecursiveOrdering(1,[[1,1,1,1,1,1,1,1]], [[1,2,3,4,5,6,7,8]]);
    end
end

# Build a Hilbert Ordering
function hilbert_ordering_rules(dim)
    if dim == 1
        return RecursiveOrdering(1,[[1]], [[1]]);
    elseif dim == 2
        return RecursiveOrdering(4,
                                [[2,3,1,1], [1,2,4,2], [3,1,3,4], [4,4,2,3]], # state indices
                                [[1,4,2,3], [1,2,4,3], [3,4,2,1], [3,2,4,1]]);# indexing order
    elseif dim == 3
        # Oh boy rotations for each state
        rots = [[],
                [1,2],
                [1,2,1,2],
                [2],
                [1,2,2],
                [3],
                [1,1,1],
                [1,1,2],
                [3,3,3],
                [2,3],
                [2,1,1,1],
                [3,3],
                [2,2,2],
                [1],
                [1,1,3],
                [1,2,2,2],
                [1,3],
                [2,2],
                [1,3,3,3],
                [3,2,2,1],
                [1,1],
                [3,3,2],
                [2,2,1],
                [2,2,3]];
        first = [1,4,2,3,8,5,7,6];
        orders = [];
        for i=1:24
            push!(orders, rotations(first, rots[i]));
        end

        # state rules
        rules = [[2,12,3,3,20,12,17,17],
                [3,1,11,21,16,1,16,21],
                [1,18,19,19,2,10,2,10],
                [], [],[], [],[],[],  # These ones aren't used, so I'll just leave them empty
                [18,3,12,11,18,20,12,20],
                [17,17,21,12,2,10,2,10],
                [11,11,1,10,19,19,1,16],
                [], [],[],
                [18,2,12,2,18,17,12,19],
                [20,16,20,16,1,18,11,11],
                [21,10,3,3,21,16,17,17],
                [20,16,20,16,3,3,21,12],
                [10,1,10,21,17,1,19,21],
                [11,11,2,18,19,19,20,18],
                [], [],[]];

        return RecursiveOrdering(24, rules, orders);
    end
end

function rotation_3d(a,axis)
    r = [5 6 1 2 7 8 3 4 ; 2 6 4 8 1 5 3 7 ; 2 4 1 3 6 8 5 7];
    b = zeros(Int,length(a));
    for i=1:length(a)
        b[r[axis,i]] = a[i];
    end
    return b;
end

function rotations(a,rots)
    b = a;
    for i=1:length(rots)
        b = rotation_3d(b,rots[i]);
    end
    return b;
end
############################################################################################

# Build the ordering for a (2^N)^D grid
function build_ordering_2d(ord, N, invert = true)
    nnodes = (2^N)^2;
    lexorder = zeros(Int,nnodes,2);
    twoN = 2^N;
    for j=1:twoN
        for i=1:twoN
            ni = i + twoN*(j-1);
            lexorder[ni,1] = i;
            lexorder[ni,2] = j;
        end
    end

    bbox = zeros(4);
    bbox[1] = 1;
    bbox[3] = 1;
    bbox[2] = twoN;
    bbox[4] = twoN;

    center = [(bbox[1]+bbox[2])/2, (bbox[3]+bbox[4])/2];
    tmpbbox = zeros(4);
    tmpcenter = zeros(2);

    result = zeros(Int64, nnodes);
    for ni=1:nnodes
        vertex = lexorder[ni,:];
        for bi=1:4
            tmpbbox[bi] = bbox[bi];
        end
        tmpcenter[1] = center[1];
        tmpcenter[2] = center[2];

        # Find the index for this node
        index = Int64(0);
        state = 1;
        for i=1:N
            stupid = index;
            stupid = state;
            index = index<<2;
            # Which quadrant is it in?
            quad = 0;
            if vertex[1] > tmpcenter[1]
                quad = quad+1;
            end
            if vertex[2] > tmpcenter[2]
                quad = quad+2;
            end
            quad += 1; # 1-based index
            #println("ni="*string(ni)*" level="*string(i)*" quad="*string(quad));

            # Add this quad's ordering position to index
            index += ord.orders[state][quad] - 1;

            # Update the state to that quad's state according to rule for this one.
            state = ord.rules[state][quad];

            # Shrink the box to one quad
            for j=1:2
                if vertex[j] > tmpcenter[j]
                    tmpbbox[2*j-1] = tmpcenter[j];
                else
                    tmpbbox[2*j] = tmpcenter[j];
                end
            end
            tmpcenter = [(tmpbbox[1]+tmpbbox[2])/2, (tmpbbox[3]+tmpbbox[4])/2];
        end

        # add one for 1-based index
        index += 1;
        #println("index for ni="*string(ni)*" is "*string(index));
        if invert
            result[ni] = index;
        else
            result[index] = ni;
        end
    end

    return result;
end

# Build the ordering for a (2^N)^D grid
function build_ordering_3d(ord, N, invert = true)
    nnodes = (2^N)^3;
    lexorder = zeros(Int,nnodes,3);
    twoN = 2^N;
    for k=1:twoN
        for j=1:twoN
            for i=1:twoN
                ni = i + twoN*((j-1) + twoN*(k-1));
                lexorder[ni,1] = i;
                lexorder[ni,2] = j;
                lexorder[ni,3] = k;
            end
        end
    end

    bbox = zeros(6);
    bbox[1] = 1;
    bbox[3] = 1;
    bbox[5] = 1;
    bbox[2] = twoN;
    bbox[4] = twoN;
    bbox[6] = twoN;

    center = [(bbox[1]+bbox[2])/2, (bbox[3]+bbox[4])/2, (bbox[5]+bbox[6])/2];
    tmpbbox = zeros(6);
    tmpcenter = zeros(3);

    result = zeros(Int64, nnodes);
    for ni=1:nnodes
        vertex = lexorder[ni,:];
        for bi=1:6
            tmpbbox[bi] = bbox[bi];
        end
        tmpcenter[1] = center[1];
        tmpcenter[2] = center[2];
        tmpcenter[3] = center[3];

        # Find the index for this node
        index = Int64(0);
        state = 1;
        for i=1:N
            stupid = index;
            stupid = state;
            index = index<<3;
            # Which octant is it in?
            octant = 0;
            if vertex[1] > tmpcenter[1]
                octant = octant+1;
            end
            if vertex[2] > tmpcenter[2]
                octant = octant+2;
            end
            if vertex[3] > tmpcenter[3]
                octant = octant+4;
            end
            octant += 1; # 1-based index
            #println("ni="*string(ni)*" level="*string(i)*" oct="*string(octant));

            # Add this octant's ordering position to index
            index += ord.orders[state][octant] - 1;

            # Update the state to that octant's state according to rule for this one.
            state = ord.rules[state][octant];

            # Shrink the box to one octant
            for j=1:3
                if vertex[j] > tmpcenter[j]
                    tmpbbox[2*j-1] = tmpcenter[j];
                else
                    tmpbbox[2*j] = tmpcenter[j];
                end
            end
            tmpcenter = [(tmpbbox[1]+tmpbbox[2])/2, (tmpbbox[3]+tmpbbox[4])/2, (tmpbbox[5]+tmpbbox[6])/2];
        end

        # add one for 1-based index
        index += 1;
        #println("index for ni="*string(ni)*" is "*string(index));
        if invert
            result[ni] = index;
        else
            result[index] = ni;
        end
    end

    return result;
end

function invert_ordering(ord)
    iord = zeros(Int,length(ord));
    for i=1:length(ord)
        iord[ord[i]] = i;
    end
    return iord;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
function bus_kernel(
    baseMVA, nbus, gen_start, line_start,
    FrStart, FrIdx, ToStart, ToIdx, GenStart, GenIdx,
    Pd, Qd, u, v, l, rho, YshR, YshI
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))
    if I <= nbus
        common_wi = 0.0
        common_ti = 0.0
        inv_rhosum_pij_ji = 0.0
        inv_rhosum_qij_ji = 0.0
        rhosum_wi_ij_ji = 0.0
        rhosum_ti_ij_ji = 0.0

        @inbounds begin
            for k=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 8*(FrIdx[k]-1)
                common_wi += l[pij_idx+4] + rho[pij_idx+4]*u[pij_idx+4]
                common_ti += l[pij_idx+6] + rho[pij_idx+6]*u[pij_idx+6]
                inv_rhosum_pij_ji += 1.0 / rho[pij_idx]
                inv_rhosum_qij_ji += 1.0 / rho[pij_idx+1]
                rhosum_wi_ij_ji += rho[pij_idx+4]
                rhosum_ti_ij_ji += rho[pij_idx+6]
            end
            for k=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 8*(ToIdx[k]-1)
                common_wi += l[pij_idx+5] + rho[pij_idx+5]*u[pij_idx+5]
                common_ti += l[pij_idx+7] + rho[pij_idx+7]*u[pij_idx+7]
                inv_rhosum_pij_ji += 1.0 / rho[pij_idx+2]
                inv_rhosum_qij_ji += 1.0 / rho[pij_idx+3]
                rhosum_wi_ij_ji += rho[pij_idx+5]
                rhosum_ti_ij_ji += rho[pij_idx+7]
            end
        end

        common_wi /= rhosum_wi_ij_ji

        rhs1 = 0.0
        rhs2 = 0.0
        inv_rhosum_pg = 0.0
        inv_rhosum_qg = 0.0

        @inbounds begin
            for g=GenStart[I]:GenStart[I+1]-1
                pg_idx = gen_start + 2*(GenIdx[g]-1)
                rhs1 += u[pg_idx] + (l[pg_idx]/rho[pg_idx])
                rhs2 += u[pg_idx+1] + (l[pg_idx+1]/rho[pg_idx+1])
                inv_rhosum_pg += 1.0 / rho[pg_idx]
                inv_rhosum_qg += 1.0 / rho[pg_idx+1]
            end

            rhs1 -= (Pd[I] / baseMVA)
            rhs2 -= (Qd[I] / baseMVA)

            for k=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 8*(FrIdx[k]-1)
                rhs1 -= u[pij_idx] + (l[pij_idx]/rho[pij_idx])
                rhs2 -= u[pij_idx+1] + (l[pij_idx+1]/rho[pij_idx+1])
            end

            for k=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 8*(ToIdx[k]-1)
                rhs1 -= u[pij_idx+2] + (l[pij_idx+2]/rho[pij_idx+2])
                rhs2 -= u[pij_idx+3] + (l[pij_idx+3]/rho[pij_idx+3])
            end

            rhs1 -= YshR[I]*common_wi
            rhs2 += YshI[I]*common_wi

            A11 = (inv_rhosum_pg + inv_rhosum_pij_ji) + (YshR[I]^2 / rhosum_wi_ij_ji)
            A12 = -YshR[I]*(YshI[I] / rhosum_wi_ij_ji)
            A21 = A12
            A22 = (inv_rhosum_qg + inv_rhosum_qij_ji) + (YshI[I]^2 / rhosum_wi_ij_ji)
            mu2 = (rhs2 - (A21/A11)*rhs1) / (A22 - (A21/A11)*A12)
            mu1 = (rhs1 - A12*mu2) / A11
            #mu = A \ [rhs1 ; rhs2]
            wi = common_wi + ( (YshR[I]*mu1 - YshI[I]*mu2) / rhosum_wi_ij_ji )
            ti = common_ti / rhosum_ti_ij_ji

            for k=GenStart[I]:GenStart[I+1]-1
                pg_idx = gen_start + 2*(GenIdx[k]-1)
                v[pg_idx] = u[pg_idx] + (l[pg_idx] - mu1) / rho[pg_idx]
                v[pg_idx+1] = u[pg_idx+1] + (l[pg_idx+1] - mu2) / rho[pg_idx+1]
            end
            for j=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 8*(FrIdx[j]-1)
                v[pij_idx] = u[pij_idx] + (l[pij_idx] + mu1) / rho[pij_idx]
                v[pij_idx+1] = u[pij_idx+1] + (l[pij_idx+1] + mu2) / rho[pij_idx+1]
                v[pij_idx+4] = wi
                v[pij_idx+6] = ti
            end
            for j=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 8*(ToIdx[j]-1)
                v[pij_idx+2] = u[pij_idx+2] + (l[pij_idx+2] + mu1) / rho[pij_idx+2]
                v[pij_idx+3] = u[pij_idx+3] + (l[pij_idx+3] + mu2) / rho[pij_idx+3]
                v[pij_idx+5] = wi
                v[pij_idx+7] = ti
            end
        end
    end
end

function bus_kernel_cpu(
    baseMVA, nbus, gen_start, line_start,
    FrStart, FrIdx, ToStart, ToIdx, GenStart, GenIdx,
    Pd, Qd, u, v, l, rho, YshR, YshI
)
    Threads.@threads for I=1:nbus
        common_wi = 0.0
        common_ti = 0.0
        inv_rhosum_pij_ji = 0.0
        inv_rhosum_qij_ji = 0.0
        rhosum_wi_ij_ji = 0.0
        rhosum_ti_ij_ji = 0.0

        @inbounds begin
            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 8*(FrIdx[k]-1)
                    common_wi += l[pij_idx+4] + rho[pij_idx+4]*u[pij_idx+4]
                    common_ti += l[pij_idx+6] + rho[pij_idx+6]*u[pij_idx+6]
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+1]
                    rhosum_wi_ij_ji += rho[pij_idx+4]
                    rhosum_ti_ij_ji += rho[pij_idx+6]
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 8*(ToIdx[k]-1)
                    common_wi += l[pij_idx+5] + rho[pij_idx+5]*u[pij_idx+5]
                    common_ti += l[pij_idx+7] + rho[pij_idx+7]*u[pij_idx+7]
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx+2]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+3]
                    rhosum_wi_ij_ji += rho[pij_idx+5]
                    rhosum_ti_ij_ji += rho[pij_idx+7]
                end
            end
        end

        common_wi /= rhosum_wi_ij_ji

        rhs1 = 0.0
        rhs2 = 0.0
        inv_rhosum_pg = 0.0
        inv_rhosum_qg = 0.0

        @inbounds begin
            if GenStart[I] < GenStart[I+1]
                for g=GenStart[I]:GenStart[I+1]-1
                    pg_idx = gen_start + 2*(GenIdx[g]-1)
                    rhs1 += u[pg_idx] + (l[pg_idx]/rho[pg_idx])
                    rhs2 += u[pg_idx+1] + (l[pg_idx+1]/rho[pg_idx+1])
                    inv_rhosum_pg += 1.0 / rho[pg_idx]
                    inv_rhosum_qg += 1.0 / rho[pg_idx+1]
                end
            end

            rhs1 -= (Pd[I] / baseMVA)
            rhs2 -= (Qd[I] / baseMVA)

            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 8*(FrIdx[k]-1)
                    rhs1 -= u[pij_idx] + (l[pij_idx]/rho[pij_idx])
                    rhs2 -= u[pij_idx+1] + (l[pij_idx+1]/rho[pij_idx+1])
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 8*(ToIdx[k]-1)
                    rhs1 -= u[pij_idx+2] + (l[pij_idx+2]/rho[pij_idx+2])
                    rhs2 -= u[pij_idx+3] + (l[pij_idx+3]/rho[pij_idx+3])
                end
            end

            rhs1 -= YshR[I]*common_wi
            rhs2 += YshI[I]*common_wi

            A11 = (inv_rhosum_pg + inv_rhosum_pij_ji) + (YshR[I]^2 / rhosum_wi_ij_ji)
            A12 = -YshR[I]*(YshI[I] / rhosum_wi_ij_ji)
            A21 = A12
            A22 = (inv_rhosum_qg + inv_rhosum_qij_ji) + (YshI[I]^2 / rhosum_wi_ij_ji)
            mu2 = (rhs2 - (A21/A11)*rhs1) / (A22 - (A21/A11)*A12)
            mu1 = (rhs1 - A12*mu2) / A11
            #mu = A \ [rhs1 ; rhs2]
            wi = common_wi + ( (YshR[I]*mu1 - YshI[I]*mu2) / rhosum_wi_ij_ji )
            ti = common_ti / rhosum_ti_ij_ji

            for k=GenStart[I]:GenStart[I+1]-1
                pg_idx = gen_start + 2*(GenIdx[k]-1)
                v[pg_idx] = u[pg_idx] + (l[pg_idx] - mu1) / rho[pg_idx]
                v[pg_idx+1] = u[pg_idx+1] + (l[pg_idx+1] - mu2) / rho[pg_idx+1]
            end
            for j=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 8*(FrIdx[j]-1)
                v[pij_idx] = u[pij_idx] + (l[pij_idx] + mu1) / rho[pij_idx]
                v[pij_idx+1] = u[pij_idx+1] + (l[pij_idx+1] + mu2) / rho[pij_idx+1]
                v[pij_idx+4] = wi
                v[pij_idx+6] = ti
            end
            for j=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 8*(ToIdx[j]-1)
                v[pij_idx+2] = u[pij_idx+2] + (l[pij_idx+2] + mu1) / rho[pij_idx+2]
                v[pij_idx+3] = u[pij_idx+3] + (l[pij_idx+3] + mu2) / rho[pij_idx+3]
                v[pij_idx+5] = wi
                v[pij_idx+7] = ti
            end
        end
    end
end

function bus_kernel_two_level(
    baseMVA, nbus, gen_start, line_start, bus_start,
    FrStart, FrIdx, ToStart, ToIdx, GenStart, GenIdx,
    Pd, Qd, v, xbar, z, l, rho, YshR, YshI
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))
    if I <= nbus
        inv_rhosum_pij_ji = 0.0
        inv_rhosum_qij_ji = 0.0

        @inbounds begin
            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+1]
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx+2]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+3]
                end
            end
        end

        rhs1 = 0.0
        rhs2 = 0.0
        inv_rhosum_pg = 0.0
        inv_rhosum_qg = 0.0

        @inbounds begin
            if GenStart[I] < GenStart[I+1]
                for g=GenStart[I]:GenStart[I+1]-1
                    pg_idx = gen_start + 2*(GenIdx[g]-1)
                    rhs1 += xbar[pg_idx] - z[pg_idx] - (l[pg_idx] / rho[pg_idx])
                    rhs2 += xbar[pg_idx+1] - z[pg_idx+1] - (l[pg_idx+1] / rho[pg_idx+1])
                    inv_rhosum_pg += 1.0 / rho[pg_idx]
                    inv_rhosum_qg += 1.0 / rho[pg_idx+1]
                end
            end

            rhs1 -= (Pd[I] / baseMVA)
            rhs2 -= (Qd[I] / baseMVA)

            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    rhs1 += (-xbar[pij_idx] + z[pij_idx]) + (l[pij_idx] / rho[pij_idx])
                    rhs2 += (-xbar[pij_idx+1] + z[pij_idx+1]) + (l[pij_idx+1] / rho[pij_idx+1])
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    rhs1 += (-xbar[pij_idx+2] + z[pij_idx+2]) + (l[pij_idx+2] / rho[pij_idx+2])
                    rhs2 += (-xbar[pij_idx+3] + z[pij_idx+3]) + (l[pij_idx+3] / rho[pij_idx+3])
                end
            end

            vi_idx = bus_start + 2*(I-1)
            common_wi = -xbar[vi_idx] + z[vi_idx] + (l[vi_idx] / rho[vi_idx])
            rhs1 += YshR[I]*common_wi
            rhs2 -= YshI[I]*common_wi

            A11 = (inv_rhosum_pg + inv_rhosum_pij_ji) + (YshR[I]^2 / rho[vi_idx])
            A12 = -YshR[I]*(YshI[I] / rho[vi_idx])
            A21 = A12
            A22 = (inv_rhosum_qg + inv_rhosum_qij_ji) + (YshI[I]^2 / rho[vi_idx])
            mu2 = (rhs2 - (A21/A11)*rhs1) / (A22 - (A21/A11)*A12)
            mu1 = (rhs1 - A12*mu2) / A11
            #mu = A \ [rhs1 ; rhs2]

            for k=GenStart[I]:GenStart[I+1]-1
                pg_idx = gen_start + 2*(GenIdx[k]-1)
                v[pg_idx] = xbar[pg_idx] - z[pg_idx] - ((l[pg_idx] + mu1) / rho[pg_idx])
                v[pg_idx+1] = xbar[pg_idx+1] - z[pg_idx+1] - ((l[pg_idx+1] + mu2) / rho[pg_idx+1])
            end
            for j=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 4*(FrIdx[j]-1)
                v[pij_idx] = xbar[pij_idx] - z[pij_idx] - ((l[pij_idx] - mu1) / rho[pij_idx])
                v[pij_idx+1] = xbar[pij_idx+1] - z[pij_idx+1] - ((l[pij_idx+1] - mu2) / rho[pij_idx+1])
            end
            for j=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 4*(ToIdx[j]-1)
                v[pij_idx+2] = xbar[pij_idx+2] - z[pij_idx+2] - ((l[pij_idx+2] - mu1) / rho[pij_idx+2])
                v[pij_idx+3] = xbar[pij_idx+3] - z[pij_idx+3] - ((l[pij_idx+3] - mu2) / rho[pij_idx+3])
            end

            vi_idx = bus_start + 2*(I-1)
            v[vi_idx] = xbar[vi_idx] - z[vi_idx] - ((l[vi_idx] - YshR[I]*mu1 + YshI[I]*mu2) / rho[vi_idx])
            v[vi_idx+1] = xbar[vi_idx+1] - z[vi_idx+1] - (l[vi_idx+1] / rho[vi_idx+1])
        end
    end

    return
end

function bus_kernel_two_level(
    baseMVA, nbus, gen_start, line_start, bus_start,
    FrStart, FrIdx, ToStart, ToIdx, GenStart, GenIdx,
    Pd, Qd, v, xbar, z, l, rho, YshR, YshI,
    slack, slack_rho
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))
    if I <= nbus
        inv_rhosum_pij_ji = 0.0
        inv_rhosum_qij_ji = 0.0

        @inbounds begin
            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+1]
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx+2]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+3]
                end
            end
        end

        rhs1 = 0.0
        rhs2 = 0.0
        inv_rhosum_pg = 0.0
        inv_rhosum_qg = 0.0

        @inbounds begin
            if GenStart[I] < GenStart[I+1]
                for g=GenStart[I]:GenStart[I+1]-1
                    pg_idx = gen_start + 2*(GenIdx[g]-1)
                    rhs1 += xbar[pg_idx] - z[pg_idx] - (l[pg_idx] / rho[pg_idx])
                    rhs2 += xbar[pg_idx+1] - z[pg_idx+1] - (l[pg_idx+1] / rho[pg_idx+1])
                    inv_rhosum_pg += 1.0 / rho[pg_idx]
                    inv_rhosum_qg += 1.0 / rho[pg_idx+1]
                end
            end

            rhs1 -= (Pd[I] / baseMVA)
            rhs2 -= (Qd[I] / baseMVA)

            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    rhs1 += (-xbar[pij_idx] + z[pij_idx]) + (l[pij_idx] / rho[pij_idx])
                    rhs2 += (-xbar[pij_idx+1] + z[pij_idx+1]) + (l[pij_idx+1] / rho[pij_idx+1])
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    rhs1 += (-xbar[pij_idx+2] + z[pij_idx+2]) + (l[pij_idx+2] / rho[pij_idx+2])
                    rhs2 += (-xbar[pij_idx+3] + z[pij_idx+3]) + (l[pij_idx+3] / rho[pij_idx+3])
                end
            end

            vi_idx = bus_start + 2*(I-1)
            common_wi = -xbar[vi_idx] + z[vi_idx] + (l[vi_idx] / rho[vi_idx])
            rhs1 += YshR[I]*common_wi
            rhs2 -= YshI[I]*common_wi

            A11 = (inv_rhosum_pg + inv_rhosum_pij_ji) + (YshR[I]^2 / rho[vi_idx]) + (1.0/slack_rho)
            A12 = -YshR[I]*(YshI[I] / rho[vi_idx])
            A21 = A12
            A22 = (inv_rhosum_qg + inv_rhosum_qij_ji) + (YshI[I]^2 / rho[vi_idx]) + (1.0/slack_rho)
            mu2 = (rhs2 - (A21/A11)*rhs1) / (A22 - (A21/A11)*A12)
            mu1 = (rhs1 - A12*mu2) / A11
            #mu = A \ [rhs1 ; rhs2]

            for k=GenStart[I]:GenStart[I+1]-1
                pg_idx = gen_start + 2*(GenIdx[k]-1)
                v[pg_idx] = xbar[pg_idx] - z[pg_idx] - ((l[pg_idx] + mu1) / rho[pg_idx])
                v[pg_idx+1] = xbar[pg_idx+1] - z[pg_idx+1] - ((l[pg_idx+1] + mu2) / rho[pg_idx+1])
            end
            for j=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 4*(FrIdx[j]-1)
                v[pij_idx] = xbar[pij_idx] - z[pij_idx] - ((l[pij_idx] - mu1) / rho[pij_idx])
                v[pij_idx+1] = xbar[pij_idx+1] - z[pij_idx+1] - ((l[pij_idx+1] - mu2) / rho[pij_idx+1])
            end
            for j=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 4*(ToIdx[j]-1)
                v[pij_idx+2] = xbar[pij_idx+2] - z[pij_idx+2] - ((l[pij_idx+2] - mu1) / rho[pij_idx+2])
                v[pij_idx+3] = xbar[pij_idx+3] - z[pij_idx+3] - ((l[pij_idx+3] - mu2) / rho[pij_idx+3])
            end

            vi_idx = bus_start + 2*(I-1)
            v[vi_idx] = xbar[vi_idx] - z[vi_idx] - ((l[vi_idx] - YshR[I]*mu1 + YshI[I]*mu2) / rho[vi_idx])
            v[vi_idx+1] = xbar[vi_idx+1] - z[vi_idx+1] - (l[vi_idx+1] / rho[vi_idx+1])

            slack[2*I-1] = (-mu1) / slack_rho
            slack[2*I] = (-mu2) / slack_rho
        end
    end

    return
end

function bus_kernel_two_level_cpu(
    baseMVA, nbus, gen_start, line_start, bus_start,
    FrStart, FrIdx, ToStart, ToIdx, GenStart, GenIdx,
    Pd, Qd, v, xbar, z, l, rho, YshR, YshI
)
    for I=1:nbus
        inv_rhosum_pij_ji = 0.0
        inv_rhosum_qij_ji = 0.0

        @inbounds begin
            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+1]
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx+2]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+3]
                end
            end
        end

        rhs1 = 0.0
        rhs2 = 0.0
        inv_rhosum_pg = 0.0
        inv_rhosum_qg = 0.0

        @inbounds begin
            if GenStart[I] < GenStart[I+1]
                for g=GenStart[I]:GenStart[I+1]-1
                    pg_idx = gen_start + 2*(GenIdx[g]-1)
                    rhs1 += xbar[pg_idx] - z[pg_idx] - (l[pg_idx] / rho[pg_idx])
                    rhs2 += xbar[pg_idx+1] - z[pg_idx+1] - (l[pg_idx+1] / rho[pg_idx+1])
                    inv_rhosum_pg += 1.0 / rho[pg_idx]
                    inv_rhosum_qg += 1.0 / rho[pg_idx+1]
                end
            end

            rhs1 -= (Pd[I] / baseMVA)
            rhs2 -= (Qd[I] / baseMVA)

            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    rhs1 += (-xbar[pij_idx] + z[pij_idx]) + (l[pij_idx] / rho[pij_idx])
                    rhs2 += (-xbar[pij_idx+1] + z[pij_idx+1]) + (l[pij_idx+1] / rho[pij_idx+1])
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    rhs1 += (-xbar[pij_idx+2] + z[pij_idx+2]) + (l[pij_idx+2] / rho[pij_idx+2])
                    rhs2 += (-xbar[pij_idx+3] + z[pij_idx+3]) + (l[pij_idx+3] / rho[pij_idx+3])
                end
            end

            vi_idx = bus_start + 2*(I-1)
            common_wi = -xbar[vi_idx] + z[vi_idx] + (l[vi_idx] / rho[vi_idx])
            rhs1 += YshR[I]*common_wi
            rhs2 -= YshI[I]*common_wi

            A11 = (inv_rhosum_pg + inv_rhosum_pij_ji) + (YshR[I]^2 / rho[vi_idx])
            A12 = -YshR[I]*(YshI[I] / rho[vi_idx])
            A21 = A12
            A22 = (inv_rhosum_qg + inv_rhosum_qij_ji) + (YshI[I]^2 / rho[vi_idx])
            mu2 = (rhs2 - (A21/A11)*rhs1) / (A22 - (A21/A11)*A12)
            mu1 = (rhs1 - A12*mu2) / A11
            #mu = A \ [rhs1 ; rhs2]

            for k=GenStart[I]:GenStart[I+1]-1
                pg_idx = gen_start + 2*(GenIdx[k]-1)
                v[pg_idx] = xbar[pg_idx] - z[pg_idx] - ((l[pg_idx] + mu1) / rho[pg_idx])
                v[pg_idx+1] = xbar[pg_idx+1] - z[pg_idx+1] - ((l[pg_idx+1] + mu2) / rho[pg_idx+1])
            end
            for j=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 4*(FrIdx[j]-1)
                v[pij_idx] = xbar[pij_idx] - z[pij_idx] - ((l[pij_idx] - mu1) / rho[pij_idx])
                v[pij_idx+1] = xbar[pij_idx+1] - z[pij_idx+1] - ((l[pij_idx+1] - mu2) / rho[pij_idx+1])
            end
            for j=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 4*(ToIdx[j]-1)
                v[pij_idx+2] = xbar[pij_idx+2] - z[pij_idx+2] - ((l[pij_idx+2] - mu1) / rho[pij_idx+2])
                v[pij_idx+3] = xbar[pij_idx+3] - z[pij_idx+3] - ((l[pij_idx+3] - mu2) / rho[pij_idx+3])
            end

            vi_idx = bus_start + 2*(I-1)
            v[vi_idx] = xbar[vi_idx] - z[vi_idx] - ((l[vi_idx] - YshR[I]*mu1 + YshI[I]*mu2) / rho[vi_idx])
            v[vi_idx+1] = xbar[vi_idx+1] - z[vi_idx+1] - (l[vi_idx+1] / rho[vi_idx+1])
        end
    end
end

function bus_kernel_two_level_cpu(
    baseMVA, nbus, gen_start, line_start, bus_start,
    FrStart, FrIdx, ToStart, ToIdx, GenStart, GenIdx,
    Pd, Qd, v, xbar, z, l, rho, YshR, YshI,
    slack, slack_rho
)
    for I=1:nbus
        inv_rhosum_pij_ji = 0.0
        inv_rhosum_qij_ji = 0.0

        @inbounds begin
            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+1]
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    inv_rhosum_pij_ji += 1.0 / rho[pij_idx+2]
                    inv_rhosum_qij_ji += 1.0 / rho[pij_idx+3]
                end
            end
        end

        rhs1 = 0.0
        rhs2 = 0.0
        inv_rhosum_pg = 0.0
        inv_rhosum_qg = 0.0

        @inbounds begin
            if GenStart[I] < GenStart[I+1]
                for g=GenStart[I]:GenStart[I+1]-1
                    pg_idx = gen_start + 2*(GenIdx[g]-1)
                    rhs1 += xbar[pg_idx] - z[pg_idx] - (l[pg_idx] / rho[pg_idx])
                    rhs2 += xbar[pg_idx+1] - z[pg_idx+1] - (l[pg_idx+1] / rho[pg_idx+1])
                    inv_rhosum_pg += 1.0 / rho[pg_idx]
                    inv_rhosum_qg += 1.0 / rho[pg_idx+1]
                end
            end

            rhs1 -= (Pd[I] / baseMVA)
            rhs2 -= (Qd[I] / baseMVA)

            if FrStart[I] < FrStart[I+1]
                for k=FrStart[I]:FrStart[I+1]-1
                    pij_idx = line_start + 4*(FrIdx[k]-1)
                    rhs1 += (-xbar[pij_idx] + z[pij_idx]) + (l[pij_idx] / rho[pij_idx])
                    rhs2 += (-xbar[pij_idx+1] + z[pij_idx+1]) + (l[pij_idx+1] / rho[pij_idx+1])
                end
            end

            if ToStart[I] < ToStart[I+1]
                for k=ToStart[I]:ToStart[I+1]-1
                    pij_idx = line_start + 4*(ToIdx[k]-1)
                    rhs1 += (-xbar[pij_idx+2] + z[pij_idx+2]) + (l[pij_idx+2] / rho[pij_idx+2])
                    rhs2 += (-xbar[pij_idx+3] + z[pij_idx+3]) + (l[pij_idx+3] / rho[pij_idx+3])
                end
            end

            vi_idx = bus_start + 2*(I-1)
            common_wi = -xbar[vi_idx] + z[vi_idx] + (l[vi_idx] / rho[vi_idx])
            rhs1 += YshR[I]*common_wi
            rhs2 -= YshI[I]*common_wi

            A11 = (inv_rhosum_pg + inv_rhosum_pij_ji) + (YshR[I]^2 / rho[vi_idx]) + (1.0/slack_rho)
            A12 = -YshR[I]*(YshI[I] / rho[vi_idx])
            A21 = A12
            A22 = (inv_rhosum_qg + inv_rhosum_qij_ji) + (YshI[I]^2 / rho[vi_idx]) + (1.0/slack_rho)
            mu2 = (rhs2 - (A21/A11)*rhs1) / (A22 - (A21/A11)*A12)
            mu1 = (rhs1 - A12*mu2) / A11
            #mu = A \ [rhs1 ; rhs2]

            for k=GenStart[I]:GenStart[I+1]-1
                pg_idx = gen_start + 2*(GenIdx[k]-1)
                v[pg_idx] = xbar[pg_idx] - z[pg_idx] - ((l[pg_idx] + mu1) / rho[pg_idx])
                v[pg_idx+1] = xbar[pg_idx+1] - z[pg_idx+1] - ((l[pg_idx+1] + mu2) / rho[pg_idx+1])
            end
            for j=FrStart[I]:FrStart[I+1]-1
                pij_idx = line_start + 4*(FrIdx[j]-1)
                v[pij_idx] = xbar[pij_idx] - z[pij_idx] - ((l[pij_idx] - mu1) / rho[pij_idx])
                v[pij_idx+1] = xbar[pij_idx+1] - z[pij_idx+1] - ((l[pij_idx+1] - mu2) / rho[pij_idx+1])
            end
            for j=ToStart[I]:ToStart[I+1]-1
                pij_idx = line_start + 4*(ToIdx[j]-1)
                v[pij_idx+2] = xbar[pij_idx+2] - z[pij_idx+2] - ((l[pij_idx+2] - mu1) / rho[pij_idx+2])
                v[pij_idx+3] = xbar[pij_idx+3] - z[pij_idx+3] - ((l[pij_idx+3] - mu2) / rho[pij_idx+3])
            end

            vi_idx = bus_start + 2*(I-1)
            v[vi_idx] = xbar[vi_idx] - z[vi_idx] - ((l[vi_idx] - YshR[I]*mu1 + YshI[I]*mu2) / rho[vi_idx])
            v[vi_idx+1] = xbar[vi_idx+1] - z[vi_idx+1] - (l[vi_idx+1] / rho[vi_idx+1])

            slack[2*I-1] = (-mu1) / slack_rho
            slack[2*I] = (-mu2) / slack_rho
        end
    end
end


text:
using Statistics
using StatsBase: sample
using Base.Threads: @threads
using BenchmarkTools
using Revise
using EvoTrees

n_obs = Int(1e6)
n_vars = 100
n_bins = 255
𝑖 = collect(1:n_obs);
𝑗 = collect(1:n_vars);
δ = rand(n_obs);
δ² = rand(n_obs);

hist_δ = zeros(n_bins, n_vars);
hist_δ² = zeros(n_bins, n_vars);
X_bin = rand(UInt8, n_obs, n_vars);

# split row ids into left and right based on best split condition
function update_set_1(set, best, x_bin)
    left = similar(set)
    right = similar(set)
    left_count = 0
    right_count = 0
    @inbounds for i in set
        if x_bin[i] <= best
            left_count += 1
            left[left_count] = i
        else
            right_count += 1
            right[right_count] = i
        end
    end
    resize!(left, left_count)
    resize!(right, right_count)
    return left, right
end

@time update_set_1(𝑖, 16, X_bin[:,1]);
@btime update_set_1($𝑖, 16, $X_bin[:,1]);
@btime update_set_1($𝑖, 64, $X_bin[:,1]);
@btime update_set_1($𝑖, 128, $X_bin[:,1]);
@btime update_set_1($𝑖, 240, $X_bin[:,1]);

# add a leaf id update - to indicate to which leaf the set is associated
function update_set_2!(leaf_vec::Vector{T}, set, best_feat, best_cond, x_bin, depth::T) where {T}
    @inbounds for i in set
        left_id = leaf_vec[i] + 2^depth
        right_id = left_id + 1
        x_bin[i, best_feat[leaf_vec[i]]] <= best_cond[leaf_vec[i]] ? leaf_vec[i] = left_id : leaf_vec[i] = right_id
    end
end

leaf_vec = ones(UInt16, n_obs);
leaf_id = 0
depth = UInt16(1)
depth = 1
best_feat = UInt16.(sample(1:100, 100000))
best_cond = rand(UInt16, 100000);

@time update_set_2!(leaf_vec, 𝑖, best_feat, best_cond, X_bin, depth);
@btime update_set_2!($leaf_vec, $𝑖, $best_feat, $best_cond, $X_bin, $depth);
Int.(leaf_vec)


text:
module LinearSolvers

using CUDA
using KernelAbstractions
using IterativeSolvers
using Krylov
using LinearAlgebra
using Printf
using SparseArrays

import ..ExaPF: xnorm, csclsvqr!
import Base: show

export bicgstab, list_solvers
export DirectSolver, BICGSTAB, EigenBICGSTAB, KrylovBICGSTAB
export get_transpose

@enum(
    SolveStatus,
    Unsolved,
    MaxIterations,
    NotANumber,
    Converged,
    Diverged,
)

include("preconditioners.jl")
include("bicgstab.jl")
include("bicgstab_eigen.jl")

abstract type AbstractLinearSolver end
abstract type AbstractIterativeLinearSolver <: AbstractLinearSolver end

"""
    list_solvers(::KernelAbstractions.Device)

List linear solvers available on current device. Currently,
only `CPU()` and `CUDADevice()` are supported.

"""
function list_solvers end

get_transpose(::AbstractLinearSolver, M::AbstractMatrix) = transpose(M)

"""
    ldiv!(solver, y, J, x)

* `solver::AbstractLinearSolver`: linear solver to solve the system
* `y::AbstractVector`: Solution
* `J::AbstractMatrix`: Input matrix
* `x::AbstractVector`: RHS

Solve the linear system `J * y = Fx

"""
function ldiv! end

"""
    DirectSolver <: AbstractLinearSolver

Solve linear system ``A * x = y`` with direct linear algebra.

* On the CPU, `DirectSolver` uses UMFPACK to solve the linear system
* On CUDA GPU, `DirectSolver` redirects the resolution to the method `CUSOLVER.csrlsvqr`

"""
struct DirectSolver <: AbstractLinearSolver end
DirectSolver(precond) = DirectSolver()
function ldiv!(::DirectSolver,
    y::Vector, J::AbstractMatrix, x::Vector,
)
    y .= J \ x
    return 0
end
function ldiv!(::DirectSolver,
    y::CuVector, J::CUDA.CUSPARSE.CuSparseMatrixCSR, x::CuVector,
)
    CUSOLVER.csrlsvqr!(J, x, y, 1e-8, one(Cint), 'O')
    return 0
end
function ldiv!(::DirectSolver,
    y::CuVector, J::CUDA.CUSPARSE.CuSparseMatrixCSC, x::CuVector,
)
    csclsvqr!(J, x, y, 1e-8, one(Cint), 'O')
    return 0
end
get_transpose(::DirectSolver, M::CUDA.CUSPARSE.CuSparseMatrixCSR) = CuSparseMatrixCSC(M)

function update!(solver::AbstractIterativeLinearSolver, J)
    update(J, solver.precond)
end

"""
    BICGSTAB <: AbstractIterativeLinearSolver
    BICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)

Custom BICGSTAB implementation to solve iteratively the linear system
``A * x = y``.
"""
struct BICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    maxiter::Int
    tol::Float64
    verbose::Bool
end
BICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false) = BICGSTAB(precond, maxiter, tol, verbose)
function ldiv!(solver::BICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P

    y[:], n_iters, status = bicgstab(J, x, P, y; maxiter=solver.maxiter,
                                        verbose=solver.verbose, tol=solver.tol)
    if status != Converged
        @warn("BICGSTAB failed to converge. Final status is $(status)")
    end
    return n_iters
end

"""
    EigenBICGSTAB <: AbstractIterativeLinearSolver
    EigenBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false)

Julia's port of Eigen's BICGSTAB to solve iteratively the linear system
``A * x = y``.
"""
struct EigenBICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    maxiter::Int
    tol::Float64
    verbose::Bool
end
EigenBICGSTAB(precond; maxiter=2_000, tol=1e-8, verbose=false) = EigenBICGSTAB(precond, maxiter, tol, verbose)
function ldiv!(solver::EigenBICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P

    y[:], n_iters, status = bicgstab_eigen(J, x, P, y; maxiter=solver.maxiter,
                                            verbose=solver.verbose, tol=solver.tol)
    if status != Converged
        error("EigenBICGSTAB failed to converge. Final status is $(status)")
    end

    return n_iters
end

struct RefBICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    verbose::Bool
end
RefBICGSTAB(precond; verbose=true) = RefBICGSTAB(precond, verbose)
function ldiv!(solver::RefBICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    y[:], history = IterativeSolvers.bicgstabl(P*J, P*x, log=solver.verbose)
    return history.iters
end

struct RefGMRES <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    restart::Int
    verbose::Bool
end
RefGMRES(precond; restart=4, verbose=true) = RefGMRES(precond, restart, verbose)
function ldiv!(solver::RefGMRES,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    y[:], history = IterativeSolvers.gmres(P*J, P*x, restart=solver.restart, log=solver.verbose)
    return history.iters
end

struct DQGMRES <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    memory::Int
    verbose::Bool
end
DQGMRES(precond; memory=4, verbose=false) = DQGMRES(precond, memory, verbose)
function ldiv!(solver::DQGMRES,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    (y[:], status) = Krylov.dqgmres(J, x, N=P, memory=solver.memory)
    return length(status.residuals)
end

"""
    KrylovBICGSTAB <: AbstractIterativeLinearSolver
    KrylovBICGSTAB(precond; verbose=false, rtol=1e-10, atol=1e-10)

Wrap `Krylov.jl` BICGSTAB algorithm to solve iteratively the linear system
``A * x = y``.
"""
struct KrylovBICGSTAB <: AbstractIterativeLinearSolver
    precond::AbstractPreconditioner
    verbose::Bool
    atol::Float64
    rtol::Float64
end
KrylovBICGSTAB(precond; verbose=false, rtol=1e-10, atol=1e-10) = KrylovBICGSTAB(precond, verbose, atol, rtol)
function ldiv!(solver::KrylovBICGSTAB,
    y::AbstractVector, J::AbstractMatrix, x::AbstractVector,
)
    P = solver.precond.P
    (y[:], status) = Krylov.bicgstab(J, x, N=P,
                                     atol=solver.atol,
                                     rtol=solver.rtol,
                                     verbose=solver.verbose)
    return length(status.residuals)
end

list_solvers(::CPU) = [RefBICGSTAB, RefGMRES, DQGMRES, BICGSTAB, EigenBICGSTAB, DirectSolver, KrylovBICGSTAB]
list_solvers(::CUDADevice) = [BICGSTAB, DQGMRES, EigenBICGSTAB, DirectSolver, KrylovBICGSTAB]

end


text:
__precompile__(true)
module NEPA
# Nonparametric Efficiency and Productivity Analysis (NEPA)

using MathProgBase
using Clp

# Export functions and types
export RS, CRS, VRS, NIRS, NDRS
export AbstractDEA, AbstractDataEnvelopment, getdata
export DEAData, getnrdmu, getiodim, getindexes, setindexes!, size, linearindexing, getindex
export DEAResult, geteff, getpeers, getsx, getsy, show
export Convex, FreeDisposal, DDF
export DEA, DEA_CRS, DEA_VRS, DEA_NIRS, DEA_NDRS
export FDH, FDH_CRS, FDH_VRS, FDH_NIRS, FDH_NDRS
export SBM
export WACM
export Luenberger,TEI,TC
export HMB,TEI,TC,SEC
export orderm
export Hyperbolic

# Source files
include("RTS.jl")
include("AbstractDEA.jl")
include("DEAData.jl")
include("DEAResult.jl")
include("DDF.jl")
include("DEA.jl")
include("FDH.jl")
include("SBM.jl")
include("WACM.jl")
include("Luenberger.jl")
include("HMB.jl")
include("OrderM.jl")
include("Hyperbolic.jl")

end

max_stars_repo_name: kepiej/NEPA


text:
# Not currently working because of Julia threading problems
type Thread
    ptr::Ptr{Void}

    function Thread(ptr::Ptr{Void})
        t = new(ptr)
        finalizer(t, destroy)
        t
    end
end

function async_send(func::Ptr{Void})
    ccall(:uv_async_send, Cint, (Ptr{Void},), func)
    C_NULL
end

function Thread(callback::Function)
    callback_c = Base.SingleAsyncWork(data -> callback())

    # async_send(func::Ptr{Void}) = (ccall(:uv_async_send, Cint, (Ptr{Void},), func); C_NULL)
    const c_async_send = cfunction(async_send, Ptr{Void}, (Ptr{Void},))

    Thread(ccall((:sfThread_create, libcsfml_system), Ptr{Void}, (Ptr{Void}, Ptr{Void},), c_async_send, callback_c.handle))
end

function destroy(thread::Thread)
    ccall((:sfThread_destroy, libcsfml_system), Void, (Ptr{Void},), thread.ptr)
end

function launch(thread::Thread)
    ccall((:sfThread_launch, libcsfml_system), Void, (Ptr{Void},), thread.ptr)
end

function wait(thread::Thread)
    println("Wait")
    ccall((:sfThread_wait, libcsfml_system), Void, (Ptr{Void},), thread.ptr)
end

function terminate(thread::Thread)
    ccall((:sfThread_terminate, libcsfml_system), Void, (Ptr{Void},), thread.ptr)
end

# real_callback(data) = print(data)

# function real_callback_handle(obj::Any)
#   work = Base.SingleAsyncWork(real_callback)
#   work.handle
# end

# function callback_other_thread(handle::Ptr{Void})
#   ccall(:uv_async_send, Cint, (Ptr{Void},), handle)
# end

# const callback_other_thread_c = cfunction(callback_other_thread, Cint, (Ptr{Void},))


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

# TODO: QuadGK is used to "hack" normalised Jacobi Polynomials.
#       to make sure they are orthonormal, not just orthogonal,
#       but this should be done properly...
#       or - indeed - move all this to a more general Pair potential
#       with arbitrary orthogonality - measure

using QuadGK

import JuLIP: evaluate,
              evaluate!,
              evaluate_d!,
              read_dict,
              write_dict

import JuLIP.MLIPs: alloc_B, alloc_dB, IPBasis

import Base.==

export Jacobi, Chebyshev

Chebyshev(N, args...; kwargs...) = Jacobi(0.5, 0.5, N, args...; kwargs...)

"""
`Jacobi{T} : ` represents the basis of Jacobi polynomials
parameterised by α, β up to some fixed maximum degree. Recall that Jacobi
polynomials are orthogonal on [-1,1] w.r.t. the weight
w(x) = (1-x)^α (1+x)^β.

### Constructor:
```
Jacobi(α, β, N)   # N = max degree
```

### Evaluate basis and derivatives:
```
x = 2*(rand() - 0.5)
P = zeros(N)
evaluate!(P, J, x)
dP = zeros(N)
evaluate_d!(P, dP, J, x)   # evaluates both P, dP
```

### Notes

`Jacobi(...)` precomputes the recursion coefficients using arbitrary
precision arithmetic, then stores them as `Vector{Float64}`. The recursion
is then given by
```
P_{n} = (A[n] * x + B[n]) * P_{n-1} + C[n] * P_{n-2}
```
"""
struct Jacobi{T} <: IPBasis
   α::T
   β::T
   A::Vector{T}
   B::Vector{T}
   C::Vector{T}
   nrm::Vector{T}
end

==(J1::Jacobi, J2::Jacobi) = (
      (J1.α == J2.α) && (J1.β == J2.β) && (length(J1) == length(J2))
   )

write_dict(J::Jacobi{T}) where {T} =
      Dict( "__id__" => "ACE_Jacobi",
            "T" => write_dict(T),
            "alpha" => J.alpha,
            "beta" => J.beta,
            "N" => length(J.A),
            "normalise" => !isempty(J.nrm) )

read_dict(::Val{SHIPs_Jacobi}, D::Dict) = read_dict(Val{ACE_Jacobi}(), D)

read_dict(::Val{:ACE_Jacobi}, D::Dict; T = read_dict(D["T"])) =
      Jacobi(T(D["alpha"]), T(D["beta"]), D["N"], T = T;
             normalise = D["normalise"])

function Jacobi(α, β, N, T=Float64; normalise=true)
   # precompute the recursion coefficients
   A = zeros(T, N)
   B = zeros(T, N)
   C = zeros(T, N)
   for n = 2:N
      c1 = big(2*n*(n+α+β)*(2*n+α+β-2))
      c2 = big(2*n+α+β-1)
      A[n] = T( big(2*n+α+β)*big(2*n+α+β-2)*c2 / c1 )
      B[n] = T( big(α^2 - β^2) * c2 / c1 )
      C[n] = T( big(-2*(n+α-1)*(n+β-1)*(2n+α+β)) / c1 )
   end
   J = Jacobi(T(α), T(β), A, B, C, T[])
   if normalise
      integrand = x -> evaluate(J, x).^2 * ((1-x)^α * (1+x)^β)
      nrm2 = quadgk(integrand, -1.0, 1.0)[1]
      J = Jacobi(T(α), T(β), A, B, C, nrm2.^(-0.5))
   end
   return J
end

Base.length(J::Jacobi) = maxdegree(J) + 1
maxdegree(J::Jacobi) = length(J.A)
alloc_B(J::Jacobi{T}, args...) where {T} = zeros(T, length(J))
alloc_dB(J::Jacobi{T}, args...) where {T} = zeros(T, length(J))

function evaluate!(P::AbstractVector, tmp, J::Jacobi, x)
   N = maxdegree(J)
   @assert length(P) >= N + 1
   @assert N+1 <= length(P)
   α, β = J.α, J.β
   P[1] = 1
   if N > 0
      #       (A * x + B) * P[1]
      #       A = 0.5 * (α+β+2)
      #       B = (α+1) - 0.5 * (α+β+2)
      P[2] = (α+1) + 0.5 * (α+β+2) * (x-1)
      if N > 1
         # 3-pt recursion
         @inbounds for n = 2:N
            P[n+1] = (J.A[n] * x + J.B[n]) * P[n] + J.C[n] * P[n-1]
         end
      end
   end
   # normalise
   if !isempty(J.nrm)  # if we want an orthonormal basis
      P .= P .* J.nrm
   end
   return P
end

function evaluate_d!(P::AbstractVector, dP::AbstractVector, tmp,
                    J::Jacobi, x::Number)
   N = maxdegree(J)
   @assert length(P) >= N+1
   @assert length(dP) >= N+1
   @assert N+1 <= min(length(P), length(dP))
   α, β = J.α, J.β
   P[1] = 1
   dP[1] = 0
   if N > 0
      P[2] = (α+1) + 0.5 * (α+β+2) * (x-1)
      dP[2] = 0.5 * (α+β+2)
      if N > 1
         @inbounds for n = 2:N
            c1 = J.A[n] * x + J.B[n]
            c2 = J.C[n]
            P[n+1] = c1 * P[n] + c2 * P[n-1]
            dP[n+1] = J.A[n] * P[n] + c1 * dP[n] + J.C[n] * dP[n-1]
         end
      end
   end
   if !isempty(J.nrm)  # if we want an orthonormal basis
      P .= P .* J.nrm
      dP .= dP .* J.nrm
   end
   # return P, dP
   return dP
end


text:
#!/usr/bin/env julia

using TimerOutputs, GLMakie, FileIO

N=4096
NPx = N
NPy = convert(Int64,round(3*NPx/4))		# 3/4 aspect ratio
const xmin	= -2.0
const xmax 	=  1.0
const ymin 	= -1.25
const ymax	=  1.25

dx = (xmax-xmin)/NPx
dy = (ymax-ymin)/NPy

xs = [ xmin + (i-1)*dx for i in 1:NPx]
ys = [ ymin + (i-1)*dy for i in 1:NPy]

to = TimerOutput()
out_array = zeros(Float64,NPx,NPy)

function mb_iterate(i,j,dx,dy,xmin,ymin,Niter)
	z = c =  xmin + (i-1) * dx + im*(ymin + (j-1) * dy)
	j = 0.0
	for i in 1:Niter
		j += 1.0
		z_prime = z*z + c
		if abs2(z_prime) <= 4.0
			z = z_prime
		else
			z = j + 0*im
			break
		end
	end
	return abs(z)
end

# compile mb_iterate with 1 call, then time it on the array.
x = mb_iterate(1,1,dx,dy,xmin,ymin,60)
print("x = $x\n") # to defeat being optimized away

@timeit to "run iterations for loops" @inbounds Threads.@threads for i in 1:NPx
	@simd for j in 1:NPy
		out_array[i,j] = mb_iterate(i,j,dx,dy,xmin,ymin,60)
	end
end

jrange = 1:NPy
irange = 1:NPx
@timeit to "run iterations for loop + broadcast" Threads.@threads for j in 1:NPy
		out_array[irange,j] = mb_iterate.(irange,j,dx,dy,xmin,ymin,30)
end

@timeit to "create figure and save" FileIO.save("mbs.png",
					heatmap(xs,ys,out_array,colormap = Reverse(:deep));
					resolution = size(out_array))

show(to)
print("\n")

max_stars_repo_name: joelandman/mbs

max_issues_repo_name: joelandman/mbs

max_forks_repo_name: joelandman/mbs


text:
import Base: filter, map, reduce

"""
    map(f, d::DTable) -> DTable

Applies `f` to each row of `d`.
The applied function needs to return a `Tables.Row` compatible object (e.g. `NamedTuple`).

# Examples
```julia
julia> d = DTable((a = [1, 2, 3], b = [1, 1, 1]), 2);

julia> m = map(x -> (r = x.a + x.b,), d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(m)
(r = [2, 3, 4],)

julia> m = map(x -> (r1 = x.a + x.b, r2 = x.a - x.b), d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(m)
(r1 = [2, 3, 4], r2 = [0, 1, 2])
```
"""
function map(f, d::DTable)
    chunk_wrap = (_chunk, _f) -> begin
        if isnonempty(_chunk)
            m = TableOperations.map(_f, _chunk)
            Tables.materializer(_chunk)(m)
        else
            _chunk
        end
    end
    chunks = map(c -> Dagger.spawn(chunk_wrap, c, f), d.chunks)
    DTable(chunks, d.tabletype)
end

"""
    reduce(f, d::DTable; cols=nothing, [init]) -> NamedTuple

Reduces `d` using function `f` applied on all columns of the DTable.

By providing the kwarg `cols` as a `Vector{Symbol}` object it's possible
to restrict the reduction to the specified columns.
The reduced values are provided in a NamedTuple under names of reduced columns.

For the `init` kwarg please refer to `Base.reduce` documentation,
as it follows the same principles.

# Examples
```julia
julia> d = DTable((a = [1, 2, 3], b = [1, 1, 1]), 2);

julia> r1 = reduce(+, d)
EagerThunk (running)

julia> fetch(r1)
(a = 6, b = 3)

julia> r2 = reduce(*, d, cols=[:a])
EagerThunk (running)

julia> fetch(r2)
(a = 6,)
```
"""
function reduce(f, d::DTable; cols=nothing::Union{Nothing, Vector{Symbol}}, init=Base._InitialValue())
    # TODO replace this with checking the colnames in schema, once schema handling gets introduced
    if length(d.chunks) > 0
        columns = isnothing(cols) ? Tables.columnnames(Tables.columns(_retrieve(d.chunks[1]))) : cols
    else
        return Dagger.@spawn NamedTuple()
    end
    col_in_chunk_reduce = (_f, _c, _init, _chunk) -> reduce(_f, Tables.getcolumn(_chunk, _c); init=deepcopy(_init))

    chunk_reduce = (_f, _chunk, _cols, _init) -> begin
        # TODO: potential speedup enabled by commented code below by reducing the columns in parallel
        v = [col_in_chunk_reduce(_f, c, _init, _chunk) for c in _cols]
        (; zip(_cols, v)...)

        # TODO: uncomment and define a good threshold for parallelization when this get's resolved
        # https://github.com/JuliaParallel/Dagger.jl/issues/267
        # This piece of code (else option) below is causing the issue above
        # when reduce is repeatedly executed or @btime is used.
        # if length(_cols) <= 1
        #     v = [col_in_chunk_reduce(_f, c, _init, _chunk) for c in _cols]
        # else
        #     values = [Dagger.spawn(col_in_chunk_reduce, _f, c, _init, _chunk) for c in _cols]
        #     v = fetch.(values)
        # end
        # (; zip(_cols, v)...)
    end
    chunk_reduce_results = [Dagger.@spawn chunk_reduce(f, c, columns, deepcopy(init)) for c in d.chunks]

    construct_single_column = (_col, _chunk_results...) -> getindex.(_chunk_results, _col)
    result_columns = [Dagger.@spawn construct_single_column(c, chunk_reduce_results...) for c in columns]

    reduce_result_column = (_f, _c, _init) -> reduce(_f, _c; init=_init)
    reduce_chunks = [Dagger.@spawn reduce_result_column(f, c, deepcopy(init)) for c in result_columns]

    construct_result = (_cols, _vals...) -> (; zip(_cols, _vals)...)
    Dagger.@spawn construct_result(columns, reduce_chunks...)
end

"""
    filter(f, d::DTable) -> DTable

Filter `d` using `f`.
Returns a filtered `DTable` that can be processed further.

# Examples
```julia
julia> d = DTable((a = [1, 2, 3], b = [1, 1, 1]), 2);

julia> f = filter(x -> x.a < 3, d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(f)
(a = [1, 2], b = [1, 1])

julia> f = filter(x -> (x.a < 3) .& (x.b > 0), d)
DTable with 2 partitions
Tabletype: NamedTuple

julia> fetch(f)
(a = [1, 2], b = [1, 1])
```
"""
function filter(f, d::DTable)
    chunk_wrap = (_chunk, _f) -> begin
        m = TableOperations.filter(_f, _chunk)
        Tables.materializer(_chunk)(m)
    end
    DTable(map(c -> Dagger.spawn(chunk_wrap, c, f), d.chunks), d.tabletype)
end


text:
@testset "Test stdout redirector" begin
    log_file, io = mktemp()
    close(io)
    try
        IS.open_file_logger(log_file) do file_logger
            multi_logger = IS.MultiLogger([file_logger], IS.LogEventTracker())
            count = 3
            Logging.with_logger(multi_logger) do
                IS.redirect_stdout_to_log() do
                    for i in 1:count
                        message = "hello"
                        println("$message $i\n$message $i")
                        sleep(0.01)
                    end
                end
            end
            events = IS.get_log_events(multi_logger.tracker, Logging.Info)
            @test length(events) == 1
            for event in events
                # There is one event per complete line.
                @test event.count == count * 2
            end
        end
    finally
        rm(log_file)
    end
end


text:
export TOC

struct TOC
  token::String
  TOC() = new("")
  TOC(tok::String) = new(tok)
end

function process_source!(toc::TOC,
                         source::AbstractString,
                         fileinfo::FileInfo;
                         args...)
  toc_re = Regex(toc.token,"is")
  if occursin(toc_re,source)
    output = ""
    toc_html = "\n\n\n
\n"
    toc_html *= "Table of Contents

\n"
    lev = 1
    sec_re = r"\n(#+)(.*)"
    count = 1
    pos = 1
    for m in eachmatch(sec_re,source)
      nlev = length(m.captures[1])
      output *= source[pos:m.offset-1]
      if nlev > 1
        output *= " \n"
        name = strip(convert(String,m.captures[2]))
        name = replace(name,r"\\cite{.*?}" => "")
        name = replace(name,r"\\onlinecite{.*?}" => "")
        for n in 1:nlev toc_html *= "  " end
        if nlev == lev+1
          toc_html *= "
"
        elseif nlev == lev-1
          toc_html *= "
"
        end
        toc_html *= "
$name
\n"
        lev = nlev
        count += 1
      end
      output *= m.match
      pos = m.offset+length(m.match)
    end
    #if has_refs
    #  toc_html *= "
References
\n"
    #end
    toc_html *= "

\n\n\n"
    output *= source[pos:end]
    #println(toc_html)
    return replace(output,toc_re => toc_html)
  end
  return source
end


text:
mutable struct AdaptiveGrid{dim,T <: Real} <: JuAFEM.AbstractGrid{dim}
    rootcells::Vector{Element}

    cells::Vector{Element} # leafcells - don't really need it, we can just collect all leaves that are masters among the cells
    faces::Vector{Element} # leaffaces - don't really need it, we can just collect all leaves that are masters among the cells
    nodes::Vector{Node{dim,T}} # actually don't need that either

    cellsets::Dict{String,Set{Element}}
    nodesets::Dict{String,Set{Node{dim,T}}}
    facesets::Dict{String,Set{Element}} # careful about master-slave-subtleties
end

function update_leafmesh!(grid::AdaptiveGrid)
    leafcells = collect_children(grid.rootcells)
    leaffacecandidates = collect(Iterators.flatten((x -> x.boundaries).(leafcells)))
    leaffaces = filter(x -> x.ismaster, leaffacecandidates)
    nodecandidates = collect(Iterators.flatten(vertices.(leafcells)))
    nodes = unique!(x -> x, nodecandidates)
    grid.cells = leafcells
    grid.faces = leaffaces
    grid.nodes = nodes
    return grid
end

function JuAFEM.getcells(grid::AdaptiveGrid)
    return (x -> eltocell(x, grid.nodes)).(grid.cells)
end

function JuAFEM.getcells(grid::AdaptiveGrid, v::Union{Int,Vector{Int}})
    if typeof(v) <: Integer return eltocell(grid.cells[v], grid.nodes) end
    return (x -> eltocell(x, grid.nodes)).(grid.cells[v])
end

function JuAFEM.getcells(grid::AdaptiveGrid, set::String)
    leaves = collect_children(grid.cellsets[set])
    return (x -> eltocell(x, grid.nodes)).(leaves)
end

function JuAFEM.getnodes(grid::AdaptiveGrid, set::String)
    return collect(grid.nodesets[set])
end

function JuAFEM.getnodeset(grid::AdaptiveGrid, set::String)
    return Set((x -> findfirst(a -> a == x, grid.nodes)).(grid.nodesets[set]))
end

function JuAFEM.getnodesets(grid::AdaptiveGrid)
    n = Dict{String,Set{Int}}()
    for d in grid.nodesets
        n[d[1]] = getnodeset(grid, d[1])
    end
    return n
end

function JuAFEM.n_faces_per_cell(grid::AdaptiveGrid)
    return nfaces(grid.cells[1])
end

function JuAFEM.getcellset(grid::AdaptiveGrid, set::String)
    s = Set(findall(x -> ultimate_ancestor(x) ∈ grid.cellsets[set], grid.cells))
    return s
end

function JuAFEM.getcellsets(grid::AdaptiveGrid)
    c = Dict{String,Set{Int}}()
    for d in grid.cellsets
        c[d[1]] = getcellset(grid, d[1])
    end
    return c
end

function getfaceidx(grid::AdaptiveGrid, f::Element)
    i = findfirst(x -> x == f.boundaryof, grid.cells)
    j = findfirst(x -> x == f, f.boundaryof.boundaries)
    return (i, j)
end

function JuAFEM.getfaceset(grid::AdaptiveGrid, set::String) # not clear if this is correct
    leaffaces = collect_children(grid.facesets[set])
    faces = Set(get_master.(leaffaces))
    S = (x -> getfaceidx(grid, x)).(faces)
    return Set(S)
end

function JuAFEM.getfacesets(grid::AdaptiveGrid)
    f = Dict{String,Set{Tuple{Int,Int}}}()
    for d in grid.facesets
        f[d[1]] = getfaceset(grid, d[1])
    end
    return f
end

# use the ones in grid.jl of JuAFEM?
_check_setname(dict, name) = haskey(dict, name) && throw(ArgumentError("there already exists a set with the name: $name"))
_warn_emptyset(set) = length(set) == 0 && @warn("no entities added to set")

function JuAFEM.addcellset!(grid::AdaptiveGrid, name::String, cellid::Union{Set{Int},Vector{Int}})
    _check_setname(grid.cellsets,  name)
    cells = Set(cellid)
    _warn_emptyset(cells)
    grid.cellsets[name] = Set((x -> grid.cells[x]).(cells))
    grid
end

# this is probably very suboptimal for our purposes
function JuAFEM.addcellset!(grid::AdaptiveGrid, name::String, f::Function; all::Bool=true)
    _check_setname(grid.cellsets,  name)
    cells = Set{Int}()
    for (i, cell) in enumerate(getcells(grid))
        pass = all
        for node_idx in cell.nodes
            node = grid.nodes[node_idx]
            v = f(node.x)
            all ? (!v && (pass = false; break)) : (v && (pass = true; break))
        end
        pass && push!(cells, i)
    end
    _warn_emptyset(cells)
    grid.cellsets[name] = Set(grid.cells[cells])
    grid
end

function JuAFEM.addfaceset!(grid::AdaptiveGrid, name::String, faceid::Set{Tuple{Int,Int}})
    _check_setname(grid.facesets, name)
    faceset = Set(faceid)
    _warn_emptyset(faceset)
    grid.facesets[name] = (x -> (grid.cells[x[1]]).boundaries[x[2]]).(faceset)
    grid
end

function JuAFEM.addfaceset!(grid::AdaptiveGrid, name::String, f::Function; all::Bool=true)# does this work?
    _check_setname(grid.facesets, name)
    faceset = Set{Tuple{Int,Int}}()
    for (cell_idx, cell) in enumerate(getcells(grid))
        for (face_idx, face) in enumerate(faces(cell))
            pass = all
            for node in face
                v = f(node.x)
                all ? (!v && (pass = false; break)) : (v && (pass = true; break))
            end
            pass && push!(faceset, (cell_idx, face_idx))
        end
    end
    _warn_emptyset(faceset)
    grid.facesets[name] = (x -> (grid.cells[x[1]]).boundaries[x[2]]).(faceset)
    grid
end

function JuAFEM.addnodeset!(grid::AdaptiveGrid, name::String, nodeid::Union{Vector{Int},Set{Int}})
    _check_setname(grid.nodesets, name)
    grid.nodesets[name] = Set((x -> grid.nodes[x]).(Set(nodeid)))
    _warn_emptyset(grid.nodesets[name])
    grid
end

function JuAFEM.addnodeset!(grid::AdaptiveGrid, name::String, f::Function)# this can be simplified
    _check_setname(grid.nodesets, name)
    nodes = Set{Int}()
    for (i, n) in enumerate(getnodes(grid))
        f(n.x) && push!(nodes, i)
    end
    grid.nodesets[name] = Set((x -> grid.nodes[x]).(Set(nodeid)))
    _warn_emptyset(grid.nodesets[name])
    grid
end

getcelltype(grid::AdaptiveGrid) = typeof(grid.cells[1])

# these two should be obsolete

function JuAFEM.getcoordinates!(x::Vector{Vec{dim,T}}, grid::AdaptiveGrid, cell::Int) where {dim,T}
    @inbounds for i in 1:length(x)
        x[i] = grid.cells[cell].nodes[i].x
    end
end

function JuAFEM.getcoordinates(grid::AdaptiveGrid, cell::Int)
    nodeex = grid.cells[cell].nodes[1]
    dim = typeof(nodeex).parameters[1]
    T = typeof(nodeex).parameters[2]
    return (n -> n.x).(grid.cells[cell].nodes)# ::Vector{Vec{dim,T}}
end

function JuAFEM.faces(e::Element, grid::AdaptiveGrid)
    getindex = (n -> findfirst(x -> x == n, grid.nodes))
    res = (f -> getindex.(f)).(faces(e))
    return res
end

function JuAFEM.edges(e::Element, grid::AdaptiveGrid)
    getindex = (n -> findfirst(x -> x == n, grid.nodes))
    res = (b -> getindex.(b)).(edges(e))
    return res
end

function JuAFEM.vertices(e::Element, grid::AdaptiveGrid)
    getindex = (n -> findfirst(x -> x == n, grid.nodes))
    res = getindex.(vertices(e))
    return res
end

function refine!(grid::AdaptiveGrid, cells::Vector{Int},  maxdepth::Int=-1)
    cs = grid.cells[cells]
    for c in cs
        if maxdepth != -1 && getdepth(c) > maxdepth error("too deep, can not split") end
        isotropic_split!(c)
        @debug check(grid)
    end
    grid = update_leafmesh!(grid)
    return grid
end

function refine!(grid::AdaptiveGrid, cells::Vector{Int}, axes::Vector{Vector{Int}},  maxdepth::Int=-1)
    cs = grid.cells[cells]
    for i in 1:length(cs)
        if maxdepth != -1 && getdepth(cs[i]) > maxdepth error("too deep, can not split") end
        split!(cs[i], axes[i])
        @debug check(grid)
    end
    grid = update_leafmesh!(grid)
    return grid
end

function derefine!(grid::AdaptiveGrid, cells::Vector{Int})
    cs = grid.cells[cells]
    par = unique!((x -> x.parent).(cs))
    if nothing ∈ par error("trying to derefine a root leaf, shame on you") end
    for c in par
        if issubset(Set(c.children), Set(cs))
            derefine!(c)
            @debug check(grid)
        else
            error("not all children of a cell marked for derefinement")
        end
    end
    grid = update_leafmesh!(grid)
    return grid
end

function getncinfo(grid::AdaptiveGrid{2})
    output = []
    for f in grid.faces
        ns, s = getinteriornodeinfo(f)
        for i in 1:length(ns)
            nodenr = findfirst(x -> x == ns[i], grid.nodes)
            masterfaceidx = getfaceidx(grid, f)
            slavefaceidxes = (x -> getfaceidx(grid, x)).(s[i])
            push!(output, [nodenr, masterfaceidx, slavefaceidxes])
        end
    end
    return output
end

function getedgeindex(grid::AdaptiveGrid, e::Element)
    i, j = getfaceidx(grid, e.boundaryof)
    k = findfirst(x -> x == e, e.boundaryof.boundaries)
    return (i, j, k)
end

getidx(grid::AdaptiveGrid{3}, e::Element{2,4,4}) = getfaceidx(grid, e)
getidx(grid::AdaptiveGrid{3}, e::Element{1,2,2}) = getedgeindex(grid, e)

function getncinfo(grid::AdaptiveGrid{3})
    output = []
    for f in grid.faces
        for n in getnodeneighborhood(f.boundaryof)
            if (node_in_el(f, n)) && !(n ∈ f.nodes)
                nodenr = findfirst(x -> x == n, grid.nodes)
                minfo = containmentinfo(f, n)
                relevantelements = (x -> containmentinfo(x, n)).(union(collect_children(f.attachedto), [f]))
                relevantelements = filter!(x -> length(x) != 0, relevantelements)
                relevantelements = collect(Iterators.flatten(relevantelements))
                slinfo = filter(x -> n ∈ vertices(x), relevantelements)
                minfo = filter(x -> !(n ∈ vertices(x)), relevantelements)
                push!(output, [nodenr, (x -> getidx(grid, x)).(minfo), (x -> getidx(grid, x)).(slinfo) ])
            end
        end
    end
    return output
end

# for testing purposes
@debug function check(grid::AdaptiveGrid)
    leafcells = collect(Iterators.flatten(collect_children.(grid.rootcells)))
    leaffacecandidates = collect(Iterators.flatten((x -> x.boundaries).(leafcells)))
    for l in leaffacecandidates
        get_master(l)
    end
end


text:
function sgd(X, y, x0, f, g,eps,lambda, maxiter, stepsize,num_data_points)
    @printf("iter \t\t j-pick \t\t fval \t\t\t diff \t\t opt \t\t alp \t\t lsiter \n")

    # safe best try => pk not always descent direction
    f_best = Inf
    x_best = x0
    x_new = x0
    x_k = x0
    oval = Inf
    val_array = zeros(0)
    stop_array = zeros(0)

    for i in 1:maxiter
        # pick idx
        idx = rand(1:num_data_points)

        x_i = X[idx:idx,:]'
        fval = f(x_new)
        gval = g(x_i, y[idx] , x_k, lambda)
        if fval < f_best
            f_best = fval
            x_best = x_new
        end
        stoppingCriteria = norm(gval,Inf)
        append!( val_array, fval )
        append!( stop_array, stoppingCriteria )

        x_new = x_k + stepsize * - gval

        x_k = x_new
        @printf("%i \t\t %i \t\t %f \t\t %f \t %e \n",i,idx,fval,(oval-fval),stoppingCriteria)

        oval = fval
    end
    return (x_k,"iter",val_array,stop_array)
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
"""
Scripts to run Ensemble Kalman Filter experiments with tensor graphical models

Author: Wayne Wang
Last modified: 08/05/2021
"""

dynamic_type = "poisson"
obs_type = "identity"
method_list = ["glasso", "kpca", "kglasso", "teralasso", "sg_palm"]
# method_list = ["sg_palm"]
T = 20
N = 25
px = py = (64, 64)
obs_noise = 0.1
process_noise = 0.1
add_process_noise = false

using SparseArrays
using LinearAlgebra
using SpecialMatrices
using BandedMatrices
using Distributions
using Random
using DataFrames
using MatrixEquations
using PDMats
using TensorToolbox
using Printf
using Plots
using PyPlot
using MIRT: jim
ENV["GKSwstype"] = 100
ENV["GKSwstype"] = "nul"

include("utils.jl")
include("glasso.jl")
include("kglasso.jl")
include("sg_palm.jl")
include("kron_pca.jl")
include("teralasso.jl")
include("utils_sim.jl")
include("enkf.jl")

# generate ground truth data
X, Y = gen_kalmanfilter_data(dynamic_type, obs_type, T, px, py, obs_noise, process_noise, add_process_noise)

# run enkf with tensor graphical models
NRMSEs_list = []
time_list = []
Omegahat_list = []
for method in method_list
    starttime = time()
    ## run enkf
    Xhat, Xhat_bar, Omegahat = enkf(Y,
                            method_str_to_type(method),
                            dynamic_type,
                            obs_type,
                            px,
                            py,
                            N,
                            obs_noise,
                            process_noise,
                            add_process_noise)
    ## timer
    stoptime = time() - starttime
    push!(time_list, stoptime)
    ## compute NRMSEs
    NRMSEs = compute_nrmse(X, Xhat)
    push!(NRMSEs_list, NRMSEs)
    ## store est. precision matrix
    push!(Omegahat_list, Omegahat)
end

# set up plots of nrmses
plt = plot()
xlabel!("Time step")
ylabel!("RMSE")
for method in method_list
    if method == "sg_palm"
        NRMSEs = NRMSEs_list[end]
    elseif method == "teralasso"
        NRMSEs = NRMSEs_list[4]
    elseif method == "kglasso"
        NRMSEs = NRMSEs_list[3]
    elseif method == "glasso"
        NRMSEs = NRMSEs_list[1]
    elseif method == "kpca"
        NRMSEs = NRMSEs_list[2]
    end
    ## plot rmse progression for each method
    plot_nrmse!(NRMSEs, method)
end
display(plt)


text:
export CubicLattice

struct CubicLattice <: Viznet.AbstractLattice
    Nx::Int
    Ny::Int
    Nz::Int
end

Base.length(lt::CubicLattice) = lt.Nx * lt.Ny * lt.Nz

function getbonds(lt::CubicLattice)
    out = Pair{Tuple{Int,Int,Int}, Tuple{Int,Int,Int}}[]
    for k=1:lt.Nz, j=1:lt.Ny, i=1:lt.Nx
        i!=lt.Nx && push!(out, (i,j,k)=>(i+1,j,k))
        j!=lt.Ny && push!(out, (i,j,k)=>(i,j+1,k))
        k!=lt.Nz && push!(out, (i,j,k)=>(i,j,k+1))
    end
    out
end

function getvertices(lt::CubicLattice)
    out = Tuple{Int,Int,Int}[]
    for k=1:lt.Nz, j=1:lt.Ny, i=1:lt.Nx
        push!(out, (i,j,k))
    end
    out
end

function rand_Js(::Type{T}, ::Ferro, lt::CubicLattice) where T
    Dict(b=>one(T) for b in getbonds(lt))
end
function rand_Js(::Type{T}, ::Randpm, lt::CubicLattice) where T
    Dict(b=>T(rand([1,-1])) for b in getbonds(lt))
end
function rand_hs(::Type{T}, ::Zero, lt::CubicLattice) where T
    Dict(b=>zero(T) for b in getvertices(lt))
end

function solve(::Type{TT}, sg::AbstractSpinglass{LT}; usecuda=false) where {LT<:CubicLattice,TT}
    # Yao gates
    lt = sg.lattice
    Lx, Ly, Lz = lt.Nx, lt.Ny, lt.Nz
    LI = LinearIndices((Lx, Ly))
    n = Lx * Ly
    reg = _init_reg(TT, n, Val(usecuda))

    for k=1:Lz
        println("Layer $k/$Lz")
        k!=1 && for j=1:Ly, i=1:Lx
            reg |> put(n, LI[i,j]=>Ghb(bondtensor(TT, sg, (i,j,k-1)=>(i,j,k))))
        end
        for j=1:Ly, i=1:Lx
            reg |> put(n, LI[i,j]=>Gh(vertextensor(TT, sg, (i,j,k))))
        end
        for j=1:Ly, i=1:Lx
            j!=Ly && (reg |> put(n, (LI[i,j],LI[i,j+1])=>Gvb(bondtensor(TT, sg, (i,j,k)=>(i,j+1,k)))))
            i!=Lx && (reg |> put(n, (LI[i,j],LI[i+1,j])=>Gvb(bondtensor(TT, sg, (i,j,k)=>(i+1,j,k)))))
        end
    end
    sum(state(reg))
end

regsize(lt::CubicLattice) = lt.Ny * lt.Nx

function _init_reg(::Type{T}, lt::CubicLattice, usecuda) where T<:TropicalTypes
    _init_reg(T, regsize(lt), usecuda)
end


text:
using DSGE
using Distributions

# Test Parameter type

# UnscaledParameter, fixed=false
α =  parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999), DSGE.SquareRoot(), Normal(0.30, 0.05), fixed=false)
@test isa(α, UnscaledParameter)
@test α.key == :α
@test isa(α.prior.value, Normal)
@test α.prior.value.μ == 0.3
@test α.description == ""
@test α.tex_label == ""
@test isa(α.transform, DSGE.SquareRoot)

# UnscaledParameter, fixed = true
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), DSGE.Untransformed(), Normal(0.30, 0.05), fixed=true)
@test α_fixed.transform_parameterization == (0.1596,0.1596)
@test isa(α_fixed.transform, DSGE.Untransformed)

# UnscaledParameter, fixed = true, transform should be overwritten given fixed
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), DSGE.SquareRoot(), Normal(0.30, 0.05), fixed=true)
@test isa(α_fixed.transform, DSGE.Untransformed)

# Fixed UnscaledParameter, minimal constructor
δ = parameter(:δ, 0.025)
@test δ.fixed
@test δ.transform_parameterization == (0.025, 0.025)
@test δ.valuebounds == (0.025, 0.025)

# Scaled parameter
β = parameter(:β, 0.1402, (1e-5, 10.), (1e-5, 10.), DSGE.Exponential(), GammaAlt(0.25, 0.1), fixed=false,  scaling = x -> (1 + x/100)\1, description="β: Discount rate.", tex_label="\\beta ")
@test isa(β, ScaledParameter)
@test isa(β.prior.value, Gamma)
@test isa(β.transform, DSGE.Exponential)

# Invalid transform
@test_throws UndefVarError α_bad = parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999),
                                              InvalidTransform(), Normal(0.30, 0.05), fixed=false)

# Arithmetic with parameters
@test promote_type(AbstractParameter{Float64}, Float16) == Float64
@test promote_type(AbstractParameter{Float64}, Int8) == Float64
## @test promote_rule(AbstractParameter{Float64}, Float16) == Float64
## @test promote_rule(AbstractParameter{Float64}, Int8) == Float64
@test δ + δ == 0.05
@test δ^2 == 0.025^2
@test -δ == -0.025
@test log(δ) == log(0.025)

# transform_to_real_line and transform_to_model_space
cx = 2 * (α - 1/2)
@test_approx_eq_eps(transform_to_real_line(α), cx / sqrt(1 - cx^2), 0.001)
@test transform_to_real_line(δ) == 0.025

m = Model990()
lastparam = parameter(:p, 0.0)
for θ in m.parameters
    isa(θ, Parameter) && (lastparam = θ)
end
@test isa(lastparam, Parameter)
@test lastparam.value == 0.0181

# transform_to_real_line and transform_to_model_space, acting on the entire parameter vector. they should be inverses!
pvec = m.parameters
vals = transform_to_real_line(pvec)
transform_to_model_space!(m, vals)
@test pvec == m.parameters

# all fixed parameters should be unchanged by both transform_to_real_line and transform_to_model_space
for θ in m.parameters
    if θ.fixed
        @test θ.value == transform_to_real_line(θ)
        @test θ.value == transform_to_model_space(θ, θ.value)
    end
end

# prior
priordensity = exp(prior(m))
@test 0 <= priordensity <= 1

# settings
# settings - boolean, string, and number. adding to model. overwriting. filestrings. testing/not testing.
n_mh_blocks = Setting(:n_mh_blocks, 22) # short constructor
reoptimize = Setting(:reoptimize, false)
vint = Setting(:data_vintage, "REF", true, "vint", "Date of data") # full constructor
@test promote_rule(Setting{Float64}, Float16) == Float64
@test promote_rule(Setting{Bool}, Bool) == Bool
@test promote_rule(Setting{ASCIIString}, AbstractString) == UTF8String
@test convert(Int64, n_mh_blocks) == 22
@test convert(ASCIIString, vint) == "REF"

@test get_setting(m, :n_mh_blocks) == m.settings[:n_mh_blocks].value
m.testing = true
@test get_setting(m, :n_mh_blocks) == m.test_settings[:n_mh_blocks].value
@test DSGE.filestring(m) == "_test"

m.testing = false
m <= Setting(:n_mh_blocks, 5, true, "mhbk", "Number of blocks for Metropolis-Hastings")
@test m.settings[:n_mh_blocks].value == 5
@test ismatch(r"^\s*_mhbk=5_vint=(\d{6})", DSGE.filestring(m))
DSGE.filestring(m, "key=val")
DSGE.filestring(m, ["key=val", "foo=bar"])
m.testing = true

# Overwriting settings
a = gensym() # unlikely to clash
b = gensym()
m <= Setting(a, 0, true, "abcd", "a")
m <= Setting(a, 1)
@test m.test_settings[a].value == 1
@test m.test_settings[a].print == true
@test m.test_settings[a].code == "abcd"
@test m.test_settings[a].description == "a"
m <= Setting(b, 2, false, "", "b")
m <= Setting(b, 3, true, "abcd", "b1")
@test m.test_settings[b].value == 3
@test m.test_settings[b].print == true
@test m.test_settings[b].code == "abcd"
@test m.test_settings[b].description == "b1"

# model paths. all this should work without errors
m.testing = true
addl_strings = ["foo=bar", "hat=head", "city=newyork"]
for fn in [:rawpath, :workpath, :tablespath, :figurespath]
    @eval $(fn)(m, "test")
    @eval $(fn)(m, "test", "temp")
    @eval $(fn)(m, "test", "temp", addl_strings)
end

nothing


text:
include("src/RequestsCache.jl")

import RequestsCache: Session, CachedSession
import RequestsCache: create_query, execute
import RequestsCache: get

session = Session()
#session = CachedSession()
session = CachedSession(cache_name="cache.jld", backend="jld", expire_after=Base.Dates.Day(1))
#println(session)

#prepared_query = create_query("GET", "http://httpbin.org/get", query = Dict("title" => "page1"))
#response = execute(prepared_query; session=session)

response = get(session, "http://httpbin.org/get"; query = Dict("title" => "page1"))

println(readall(response))