text:
@testset "utils tests" begin
    @testset "knn_search" begin
        @inferred knn_search(rand(5, 100), 5, Euclidean())
        @inferred knn_search(rand(5, 100), 5, Euclidean(), Val(:approximate))
        @testset "pairwise tests" begin
            data = [0. 0. 0.; 0. 1.5 2.]
            true_knns = [2 3 2; 3 1 1]
            true_dists = [1.5 .5 .5; 2. 1.5 2.]
            knns, dists = knn_search(data, 2, Euclidean(), Val(:pairwise))
            @test knns == true_knns
            @test dists == true_dists
        end

        @testset "precomputed tests" begin
            dist_mat = [0. 2. 1.;
                        2. 0. 3.;
                        1. 3. 0.]
            true_knns = [3 1 1;
                         2 3 2]
            true_dists = [1 2 1;
                          2 3 3]
            knns, dists = knn_search(dist_mat, 2, :precomputed)
            @test knns == true_knns
            @test dists == true_dists

        end

        @testset "query data tests" begin
            X = [0.0 -2.0 0.0; 0.0 -1.0 2.0]
            Q = [-2 0; 0 -1]
            true_knns = [2 1; 1 2]
            true_dists = [1.0 1.0; 2.0 2.0]
            knns = [3 1 1; 2 3 2]
            dists = [2.0 2.236 2.0; 2.236 3.606 3.606]
            ret_knns, ret_dists = knn_search(X, Q, 2, Euclidean(), knns, dists)
            @test ret_knns == true_knns
            @test isapprox(ret_dists, true_dists, atol=1e-4)
        end
    end

    @testset "combine_fuzzy_sets tests" begin
        A = [1.0 0.1; 0.4 1.0]
        union_res = [1.0 0.46; 0.46 1.0]
        res = combine_fuzzy_sets(A, 1.0)
        @test isapprox(res, union_res)
        inter_res = [1.0 0.04; 0.04 1.0]
        res = combine_fuzzy_sets(A, 0.0)
        @test isapprox(res, inter_res)

    end

    @testset "fit_ab tests" begin
        @test all((1, 2) .== fit_ab(-1, 0, 1, 2))
    end
end


text:
include_dependency("../../dd/construction.jl")
include_dependency("../../dd/modification.jl")
include_dependency("../../dd/computation.jl")

using JuMP, Clp, CPLEX

"""
Generate CGLP cuts.

# Input
- `dd::DecisionDiagram`: DD.
- `fractional_point::Array{Float64}`: The point to be separated (as a vector).
- `tilim::Float64`: Time limit in seconds (optional).

# Output
- `obj`: The optimal value of the CGLP.
- `coefficients`: The coefficient vector of the cut.
- `rhs`: The right-hand-side of the cut.
- `status`: The status of the CGLP.
"""
function CGLP(dd::DecisionDiagram, fractional_point::Array{Float64, 1}; tilim::Float64 = 100000)::Tuple{Float64, Array{Float64, 1}, Float64, Symbol}

    n::Int = nvars(dd)     #the number of arc layers of dd (number of variables)
    @assert(n == length(fractional_point))       #makes sure the size of the fractional point matches the dimension of variables represented by dd

    node_num::Int = nnodes(dd)
    arc_num::Int = narcs(dd)

    C_layersize::Array{Int, 1} = Array{Int, 1}(  0)     # compute the cumulative number of nodes in previous layers of the indexed layer
    push!(C_layersize, 0)
    for i::Int = 2::n+1
        push!(C_layersize, C_layersize[i-1] + get_node_layer_size(dd, i-1))
    end

    # Preparing an optimization model
    m = Model(solver = CplexSolver(CPX_PARAM_TILIM = tilim))

    # Defining variables
    # NOTE: Jump does not support Array{Array} of variables in the definition
    # There are two options:
    # 1: Use anonymous variable definition which has a few restrictions (cannot get index for a copy model)
    # 2: To define "flattend out" theta variables, by numbering them manually
    @variable(m, theta_plus[1:node_num] >= 0)
    @variable(m, theta_minus[1:node_num] >= 0)
    @variable(m, gamma_plus[1:n] >= 0)
    @variable(m, gamma_minus[1:n] >= 0)

    # Fixing the theta of the source node at zero
    setupperbound(theta_plus[1], 0)
    setupperbound(theta_minus[1], 0)

    # Defining objective function
    @objective(m, Max, sum((gamma_plus[i] - gamma_minus[i])*fractional_point[i] for i in 1:n) - theta_plus[node_num] + theta_minus[node_num])

    # Adding projection cone constraints without accessing the inner graph information
    for i::Int=1:n
        for j::int=1:get_arc_layer_size(dd, i)
            @constraint(m, theta_plus[get_arc_tail(dd, i, j) + C_layersize[i]] - theta_minus[get_arc_tail(dd, i, j) + C_layersize[i]] - theta_plus[get_arc_head(dd, i, j) + C_layersize[i+1]] + theta_minus[get_arc_head(dd, i, j) + C_layersize[i+1]] + (gamma_plus[i] - gamma_minus[i])*get_arc_label(dd, i, j) <= 0)
        end
    end

    # Adding normalization constraints
    @constraint(m, sum(gamma_plus[i] + gamma_minus[i] for i in 1:n) <= 1)

    # Solving the model
    status::Symbol = solve(m)

    # Getting the output
    obj::Float64 = getobjectivevalue(m)  #gets the optimal value

    println("--------------")
    println("status: ", status)
    println("CGLP obj: ", obj)

    coefficients::Array{FLoat64, 1} = zeros(n)   #coefficient of the CGLP inequality
    for i::Int=1:n
        coefficients[i] = getvalue(gamma_plus[i]) - getvalue(gamma_minus[i])
    end
    rhs::Float64 = getvalue(theta_plus[node_num]) - getvalue(theta_minus[node_num])  #gets the rhs value for the CGLP inequality

    return obj, coefficients, rhs, status
end

"""
Generate cuts via the subgradient method.

# Input
- `dd::DecisionDiagram`: DD.
- `fractional_point::Array{Float64}`: The point to be separated (as a vector).
- `step_rule::Int64`: Step size rule (optional).
- `starting_point::Array{Float64}`: Starting point (optional).
- `tilim::Float64`: Time limit in seconds (optional).

# Output
- `best_obj`: The optimal value of the subgradient method at termination.
- `best_coefficients`: The coefficient vector of the best inequality at termination.
- `best_rhs`: The right-hand-side of the best inequality at termination.
"""
function subgradient(dd::DecisionDiagram, fractional_point::Array{Float64, 1}; step_rule::Int64 = 0, starting_point::Array{Float64, 1} = zeros(length(fractional_point)), tilim::Float64 = 100000.0)::Tuple{Float64, Array{Float64, 1}, Float64}

    n::Int = nvars(dd)     #the number of arc layers of dd
    @assert(n == length(fractional_point))       #makes sure the size of the fractional point matches the dimension of variables represented by dd

    # Choosing the step rule
    # 0. constant step size: \rho_t = c for some c>0
    # 1. constraint step length: \rho_t = c/||g_t||, where c>0 and g_t is the subgradient at step t
    # 2. diminishing step rule: \rho_t = c/sqrt(t), where c>0
    # 3. square summable diminishing step rule: \rho_t = c/(a+t), where c>0 and a>=0
    if step_rule < 0 || step_rule > 3
        step_rule = 0
        info("Mismatch in the step size rule. The default rule has been used.")
    end

    c0::Int = 1
    a0::Int = 0
    if step_rule == 0
        rho = (t::Int64,norm_g_t::Float64) -> c0
    elseif step_rule == 1
        rho = (t::Int64,norm_g_t::Float64) -> c0/norm_g_t
    elseif step_rule == 2
        rho = (t::Int64,norm_g_t::Float64) -> c0/sqrt(t)
    elseif step_rule == 3
        rho = (t::Int64,norm_g_t::Float64) -> c0/(a0 + t)
    end

    # initialization
    best_obj::Float64 = 0            # best objecive found so far
    best_coefficients::Array{Float64, 1} = zeros(n)      # coefficient vector corresponding to the best objective
    best_rhs::Float64 = 0            # right-hand-side value corresponding to the best objective
    delta::Float64 = 0               # the violation value (objective value)
    subgrad::Array{Float64, 1} = ones(n)       # subgradient vector
    iter::Int = 1                # iteration number
    stop_criterion::Array{Int, 1} = [0, 1]                # stopping criterion
    if norm(starting_point) > 0         # normalizing the starting point and making a copy
        gamma::Array{Float64, 1} = starting_point/norm(starting_point)
    else
        gamma = copy(starting_point)
    end

    # stopping criteria for the entire algorithm
    iter_max::Int = 20      # the maximum number of iterations to execute the algorithm
    tolerance::Float64 = 0.0005    # the tolerance for relative error
    tol_num::Int = 2         # the number of previous objective values (including the current one) wrt which the tolerance is computed
    obj_history::Array{Float64, 1} = Array{Float64, 1}(  tol_num)    # stores the previous (improved) objective values for tolerance check

    # function to define what stopping criterion must be used:
    # 0: time
    # 1: iteration number
    # 2: concavity error: obtained from the first order concavity inequality over the unit ball normalization constraint: at iteration t, the objective error is <= norm of subgradient - violation value
    # 3: objective improvement: the relative difference between k consequtive improvements
    function stop_rule(st::Array{Int64})::Bool   # multiple stopping criteria can be passed as an array
        p = true
        for v::Int in st
            if v == 0         # compares time elapsed with given time limit
                clock_end::Float64 = time()
                time_elapsed::Float64 = float(clock_end - clock_start)
                p = p && (time_elapsed <= tilim)
            end
            if v == 1       # uses the number of iteration criterion
                p = p && (iter <= iter_max)
            end
            if v == 2       # uses the concavity error criterion
                p = p && ((norm(subgrad) - delta)/max(norm(subgrad), 1e-9) > tolerance)
            end
            if v == 3       # uses the objective improvement criterion (ONLY considers positive violation values)
                if obj_history[end] != best_obj   # if the objective value is updated, the obj_history vector must be updated as well and the criterion is checked
                    push!(obj_history, best_obj)
                    popfirst!(obj_history)
                    p = p && ((best_obj - obj_history[1])/best_obj > tolerance)
                end
            end
        end
        return p    # is false when the stopping criterion is met
    end

    # initializing time parameter
    clock_start::Float64 = time()

    # Main algorithm
    while stop_rule(stop_criterion)

        # computing the longest path with the objective function defined by the current direction
        lp1::Array{Float64, 1}, lpv1::Float64 = longest_path(dd, gamma)

        # computing the violatation value of the current inequality at the given fractional point
        delta = gamma'*fractional_point - lpv1

        if delta > best_obj                 # since best_obj >= 0 by construction, this means that a cut is detected
            if best_obj <= 0                # i.e., the current inequality is the first that cuts off the fractional point
                iter_max = iter + 20                   # reduce the number of iterations after detecting a cut
                stop_criterion = [0, 1]             # changes the stop rule to the one that considers both iteration number and concavity error
            end
            best_obj = delta
            best_coefficients = gamma
            best_rhs = lpv1
        end

        # computing the subgradient at current iteation
        subgrad = fractional_point - lp1

        # updating the coefficient vector
        phi::Array{Float64, 1} = gamma + rho(iter,norm(subgrad))*subgrad

        # projecting the updated vector onto the unit ball normalization cosntraint
        #if norm(phi) > 1
        #    gamma = phi./norm(phi)
        #else
            gamma = phi
        #end

        iter += 1

    end

    return best_obj, best_coefficients, best_rhs

end

ext: jl

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
using ApproxFun, MultivariateOrthogonalPolynomials, Plots, BlockArrays
    import MultivariateOrthogonalPolynomials: Vec, DirichletTriangle

# Neumann

d = Triangle()
S = DirichletTriangle{1,1,1}(d)
Dx = Derivative(S,[1,0])
Dy = Derivative(S,[0,1])
B₁ = I : S → Legendre(Vec(0,0) .. Vec(0,1))
B₂ = I : S → Legendre(Vec(0,0) .. Vec(1,0))
B₃ = I : S → Legendre(Vec(0,1) .. Vec(1,0))

L = [0 B₁ 0;
 0  0  B₂;
 0  B₃ B₃;
 Dx -eye(S) 0;
  Dy 0 -eye(S);
   0.1eye(S) Dx Dy]

f = Fun((x,y) -> exp(x*cos(y)), d)

u = \(L , [0;0; 0; 0;0; f]; tolerance=1E-5)
plot(u[1])
\

h=0.01;
 (u[1](0.5,0.5)-u[1](0.5-h,0.5-h))/h
plot(Dx*u[1])
plot(Dy*u[1])
plot((Dx+Dy)*u[1])
plot(u[3])

rhs = Fun([0;0; 0; 0;0; f], rangespace(L))

A = Matrix(L[Block.(1:10), Block.(1:10)])

u = svdfact(A) \ pad(rhs.coefficients, size(A,1))

Laplacian()*Fun(domainspace(L),u)[1] - f |> coefficients

Dx*Fun(domainspace(L),u)[1]  - Fun(domainspace(L),u)[2] |> coefficients

Dx*Fun(domainspace(L),u)[2] + Dy*Fun(domainspace(L),u)[3] - f |> coefficients

plot(Fun(domainspace(L),u)[1])

u
[Dx eye(S)]

L = []

[Dx[

# D to N
d = Triangle() , Triangle(Vec(1,1),Vec(0,1),Vec(1,0))
S = DirichletTriangle{1,1,1}.(d)
B₁₁ = I : S[1] → Legendre(Vec(0,0) .. Vec(0,1))
B₁₂ = I : S[1] → Legendre(Vec(0,0) .. Vec(1,0))
B₂₁ = I : S[2] → Legendre(Vec(1,1) .. Vec(0,1))
B₂₂ = I : S[2] → Legendre(Vec(1,1) .. Vec(1,0))


text:
function projection_points(gp::GPE, border::B, maxdist::Float64) where {B<:BorderType}
    X∂ = Array{Float64}(undef, 2, gp.nobs)
    distances = Vector{Float64}(undef, gp.nobs)
    for i in 1:gp.nobs
        # obtain coordinates for treatment point
        x,y = gp.x[:,i]
        point = LibGEOS.Point(x,y)
        # projection onto border (as distance along border)
        distances[i] = LibGEOS.distance(border, point)
        proj_point = nearestPoints(border, point)[1]
        # get border point from distance
        # proj_point = interpolate(border, proj_dist)
        # get coordinates from point
        proj_x, proj_y = coordinates(proj_point)
        X∂[1,i] = proj_x
        X∂[2,i] = proj_y
    end
    return X∂[:, distances .<= maxdist]
end

function proj_estimator(gpT::GPE, gpC::GPE, border::B, maxdist::Float64) where {B<:BorderType}
    X∂_treat = projection_points(gpT, border, maxdist)
    X∂_ctrol = projection_points(gpC, border, maxdist)
    X∂ = [X∂_treat X∂_ctrol]

    μpost, Σpost = cliff_face(gpT, gpC, X∂)
    return unweighted_mean(μpost, Σpost)
end

function data_hull(gpT::GPE, gpC::GPE)
    # Obtain a convex hull containing all the data.
    X_multi_treat = LibGEOS.MultiPoint([gpT.x[:, i] for i in 1:gpT.nobs])
    X_multi_ctrol = LibGEOS.MultiPoint([gpC.x[:, i] for i in 1:gpC.nobs])
    convexhull_treat = LibGEOS.convexhull(X_multi_treat)
    convexhull_ctrol = LibGEOS.convexhull(X_multi_ctrol)
    hull = LibGEOS.union(convexhull_treat, convexhull_ctrol)
    return hull
end

function infinite_proj_sentinels(gpT::GPE, gpC::GPE, border::B,
                                    region::MultiPolygon,
                                    maxdist::Float64, gridspace::Float64;
                                    density = ((s1,s2) -> 1.0)) where {B<:BorderType}
    # Add a buffer around the border, returns a "border area".
    #=buffer_polygon = LibGEOS.buffer(border, buffer)=#

    # Get the minimum and maximum of X1 and X2 within the region
    env = envelope(region)
    env_coord = coordinates(env)[1]
    X1_min = minimum(xcoord(p) for p in env_coord)
    X1_max = maximum(xcoord(p) for p in env_coord)
    X2_min = minimum(ycoord(p) for p in env_coord)
    X2_max = maximum(ycoord(p) for p in env_coord)

    # Obtain grid of points that covers the region
    X1_grid = X1_min:gridspace:X1_max
    X2_grid = X2_min:gridspace:X2_max

    border_envelope = envelope(border)

    projected_weights = Dict{Vector{Float64}, Float64}()
    for (s1,s2) in Iterators.product(X1_grid, X2_grid)
        # convert to a point obejct
        p = LibGEOS.Point(s1,s2)
        # Only keep points that are both within `maxdist` of
        # the border, and are in the convex hull of the data.
        if LibGEOS.distance(p, border_envelope) > maxdist
            continue
        elseif LibGEOS.distance(p, border) > maxdist
            continue
        elseif !LibGEOS.within(p, region)
            # Note: the `LibGEOS.jl` point-in-polygon implementation is
            # stupidly slow. The `ClimateTools.jl` package provides a simple
            # `inpoly(point, polygon)` function that could be adapted and
            # is about 10 times faster that `LibGEOS.within`
            continue
        end
        # compute the population density at this location
        ρ = density(s1,s2)

        # project the point onto the border
        projected = LibGEOS.nearestPoints(p, border)[2]
        # obtain the coordinates of the projected point
        proj_coords = GeoInterface.coordinates(projected)[1:2]
        if haskey(projected_weights, proj_coords)
            projected_weights[proj_coords] += ρ
        else
            # initialize
            projected_weights[proj_coords] = ρ
        end
    end

    unique_projected_points = collect(keys(projected_weights))
    X∂_projected = [[p[1] for p in unique_projected_points]';
                    [p[2] for p in unique_projected_points]']
    # And the counts as the weights.
    weights = collect(values(projected_weights))
    # Note: assumes keys() and values() are in the same order
    return X∂_projected, weights
end

function infinite_proj_estim(gpT::GPE, gpC::GPE, border::B,
                                    region::MultiPolygon,
                                    maxdist::Float64, gridspace::Float64;
                                    density = ((s1,s2) -> 1.0)) where {B<:BorderType}
    X∂_projected, weights = infinite_proj_sentinels(
            gpT, gpC, border, region, maxdist, gridspace; density=density
            )

    # Obtain the cliff face for these sentinels...
    μpost, Σpost = cliff_face(gpT, gpC, X∂_projected)
    # ... and apply a weighted mean estimator.
    return weighted_mean(μpost, Σpost, weights)
end


text:
using Test
using ManifoldsBase
using ManifoldsBase: AbstractNumbers, ℝ, ℂ, NestedReplacingPowerRepresentation
using StaticArrays

power_array_wrapper(::Type{NestedPowerRepresentation}, ::Int) = identity
power_array_wrapper(::Type{NestedReplacingPowerRepresentation}, i::Int) = SVector{i}

struct TestArrayRepresentation <: AbstractPowerRepresentation end

@testset "Power Manifold" begin

    @testset "Power Manifold with a test representation" begin
        M = ManifoldsBase.DefaultManifold(3)
        N = PowerManifold(M, TestArrayRepresentation(), 2)
        O = PowerManifold(N, TestArrayRepresentation(), 3) # joins instead of nesting.
        @test repr(O) ==
              "PowerManifold(DefaultManifold(3; field = ℝ), TestArrayRepresentation(), 2, 3)"
        p = zeros(6)
        X = zeros(6)
        @test ManifoldsBase.check_power_size(N, p) === nothing
        @test ManifoldsBase.check_power_size(O, p) isa DomainError
        @test ManifoldsBase.check_power_size(N, p, X) === nothing
        @test ManifoldsBase.check_power_size(O, p, X) isa DomainError
    end

    @testset "PowerManifold and allocation with empty representation size" begin
        M = ManifoldsBase.DefaultManifold()
        N = PowerManifold(M, NestedPowerRepresentation(), 2)
        p = [1, 1]
        X = [2, 2]
        # check - though only because this function exists for avoiding ambiguities.
        cm = ManifoldsBase.allocate_result(N, get_coordinates, p, X, DefaultBasis())
        @test size(X) == size(cm)
    end

    for PowerRepr in [NestedPowerRepresentation, NestedReplacingPowerRepresentation]
        @testset "PowerManifold with $(PowerRepr)" begin
            M = ManifoldsBase.DefaultManifold(3)
            for pow_size in ((2,), (2, 3))
                pow_cart = CartesianIndices(pow_size)
                N = PowerManifold(M, PowerRepr(), pow_size...)
                wrapper_2 = power_array_wrapper(PowerRepr, 2)
                wrapper_3 = power_array_wrapper(PowerRepr, 3)

                p = wrapper_3.([i[1] == 1 ? zeros(3) : ones(3) for i in pow_cart])
                q = wrapper_3.([i[1] == 1 ? ones(3) : zeros(3) for i in pow_cart])

                pow_string = join(pow_size, ", ")

                @testset "Constructors" begin
                    @test repr(N) ==
                          "PowerManifold(DefaultManifold(3; field = ℝ), $(PowerRepr)(), $pow_string)"
                    # add to type
                    @test repr(N^3) ==
                          "PowerManifold(DefaultManifold(3; field = ℝ), $(PowerRepr)(), $pow_string, 3)"
                    # add to type
                    @test repr(PowerManifold(N, 3)) ==
                          "PowerManifold(DefaultManifold(3; field = ℝ), $(PowerRepr)(), $pow_string, 3)"
                    # nest
                    @test repr(PowerManifold(N, PowerRepr(), 3)) ==
                          "PowerManifold(PowerManifold(DefaultManifold(3; field = ℝ), $(PowerRepr)(), $pow_string), $(PowerRepr)(), 3)"
                end

                @test N^3 == PowerManifold(M, PowerRepr(), pow_size..., 3)

                if pow_size == (2,)

                    @test ManifoldsBase.get_iterator(N^3) ==
                          Base.product(Base.OneTo(2), Base.OneTo(3))
                    @testset "point/tangent checks" begin
                        pE1 = [wrapper_3(zeros(3)), wrapper_2(ones(2))] # one component wrong
                        pE2 = wrapper_2.([zeros(2), ones(2)]) # both wrong
                        @test is_point(N, p)
                        @test is_point(N, p, true)
                        @test !is_point(N, pE1)
                        @test !is_point(N, pE2)
                        @test_throws ComponentManifoldError is_point(N, pE1, true)
                        @test_throws CompositeManifoldError is_point(N, pE2, true)
                        # tangent - test base
                        @test is_vector(N, p, p)
                        @test is_vector(N, p, p, true)
                        @test !is_vector(N, pE1, p)
                        @test !is_vector(N, pE2, p)
                        # tangents - with proper base
                        @test is_vector(N, p, p, true)
                        @test !is_vector(N, p, pE1)
                        @test !is_vector(N, p, pE2)
                        @test_throws ComponentManifoldError is_vector(N, p, pE1, true)
                        @test_throws CompositeManifoldError is_vector(N, p, pE2, true)
                    end
                end
                @testset "specific functions" begin
                    @test distance(N, p, q) == sqrt(sum(distance.(Ref(M), p, q) .^ 2))
                    @test exp(N, p, q) == p .+ q

                    @test retract(N, p, q) == p .+ q
                    @test retract(N, p, q, ExponentialRetraction()) == p .+ q
                    r = allocate(p)
                    @test retract!(N, r, p, q, ExponentialRetraction()) == p .+ q
                    @test r == p .+ q
                    @test inverse_retract(N, p, r) == q
                    @test inverse_retract(N, p, r, LogarithmicInverseRetraction()) == q
                    X = allocate(p)
                    @test inverse_retract!(N, X, p, r, LogarithmicInverseRetraction()) == q
                    @test X == q

                    if pow_size == (2,)
                        @test ManifoldsBase.get_iterator(N) == Base.OneTo(2)
                    else
                        @test ManifoldsBase.get_iterator(N) ==
                              Base.product(Base.OneTo(2), Base.OneTo(3))
                    end
                    @test injectivity_radius(N) == injectivity_radius(M)
                    @test injectivity_radius(N, ExponentialRetraction()) ==
                          injectivity_radius(M)
                    @test injectivity_radius(N, p) == injectivity_radius(M, p)
                    p2 = allocate(p)
                    copyto!(N, p2, p)
                    @test !(p === p2)
                    @test p == p2
                    q2 = allocate(q)
                    copyto!(N, q2, p, q)
                    @test !(q === q2)
                    @test q == q2
                    @test inner(N, p, q, q) == sum(inner.(Ref(M), p, q, q))
                    @test isapprox(N, p, q) == (all(isapprox.(Ref(M), p, q)))
                    @test isapprox(N, p, p) == (all(isapprox.(Ref(M), p, p)))
                    @test isapprox(N, p, q, p) == (all(isapprox.(Ref(M), p, q, p)))
                    @test isapprox(N, p, p, p) == (all(isapprox.(Ref(M), p, p, p)))
                    @test log(N, p, q) == q .- p
                    @test inverse_retract(N, p, q) == q .- p
                    @test manifold_dimension(N) == prod(pow_size) * manifold_dimension(M)
                    @test mid_point(N, p, q) == mid_point.(Ref(M), p, q)
                    @test sqrt(inner(N, p, q, q)) ≈ norm(N, p, q)
                    @test project(N, p) == p
                    @test project(N, p, q) == q
                    @test power_dimensions(N) == pow_size
                    @test power_dimensions(N^3) == (pow_size..., 3)
                    m = ParallelTransport()
                    @test vector_transport_to(N, p, p, q) == p
                    @test vector_transport_to(N, p, p, q, m) == p
                    @test parallel_transport_to(N, p, p, q) == p
                    q2 = [zeros(3) for _ in pow_cart]
                    vector_transport_to!(N, q2, p, q, p)
                    @test q2 == q
                    q2 = [zeros(3) for _ in pow_cart]
                    parallel_transport_to!(N, q2, p, q, p)
                    @test q2 == q
                    @test vector_transport_direction(N, p, p, q) == p
                    @test vector_transport_direction(N, p, p, q, m) == p
                    @test parallel_transport_direction(N, p, p, q) == p
                    q2 = [zeros(3) for _ in pow_cart]
                    vector_transport_direction!(N, q2, p, q, p)
                    @test q2 == q
                    q2 = [zeros(3) for _ in pow_cart]
                    parallel_transport_direction!(N, q2, p, q, p)
                    @test q2 == q
                    if pow_size == (2,)
                        @test p[N, 1] == p[1]
                        p[N, 1] = 2 .* ones(3)
                        @test p[N, 1] == 2 .* ones(3)
                        @test p[N, [2, 1]] == [p[2], p[1]]
                        if PowerRepr == NestedPowerRepresentation
                            @test view(p, N, 2) isa SubArray
                            @test view(p, N, 2) == p[2]
                        end
                    end
                end
                @testset "Basis, coordinates & vector" begin
                    v = get_coordinates(N, p, q, DefaultBasis())
                    @test v == reduce(vcat, q)
                    @test get_vector(N, p, v, DefaultBasis()) == q
                    B = get_basis(N, p, DefaultBasis())
                    @test get_coordinates(N, p, q, B) == v
                    # the method tested below should not be used but it prevents ambiguities from occurring
                    # and the test is here to make coverage happy
                    @test ManifoldsBase.allocate_result(N, get_coordinates, p, q, B) isa
                          Vector
                    v2 = similar(v)
                    get_coordinates!(N, v2, p, q, B)
                    @test v2 == v
                    @test get_coordinates(N, p, q, DefaultOrthonormalBasis()) == v
                    v3 = similar(v2)
                    @test get_coordinates!(N, v3, p, q, DefaultOrthonormalBasis()) == v
                    @test v3 == v
                    @test B.data.bases[1].data == get_basis(M, p[1], DefaultBasis()).data
                    @test B.data.bases[2].data == get_basis(M, p[2], DefaultBasis()).data
                    @test get_vector(N, p, v, B) == q
                    q2 = similar.(q)
                    @test get_vector!(N, q2, p, v, B) == q
                    @test q2 == q
                    q3 = similar.(q)
                    @test get_vector(N, p, v, DefaultOrthonormalBasis()) == q
                    @test get_vector!(N, q3, p, v, DefaultOrthonormalBasis()) == q
                    @test q3 == q
                    @test zero_vector(N, p) == [zeros(3) for _ in pow_cart]
                    B2 = DiagonalizingOrthonormalBasis([ones(3) for _ in pow_cart])
                    B3 = get_basis(N, p, B2)
                    if pow_size == (2,)
                        @test sprint(show, "text/plain", B) ==
                              """$(DefaultBasis()) for a power manifold
          Basis for component (1,):
          $(sprint(show, "text/plain", B.data.bases[1]))
          Basis for component (2,):
          $(sprint(show, "text/plain", B.data.bases[2]))
          """
                    end

                    bv = get_vectors(N, p, B)
                    @test length(bv) == number_of_coordinates(N, B)
                    for i in 1:number_of_coordinates(N, B)
                        @test bv[i] == get_vector(
                            N,
                            p,
                            ManifoldsBase._euclidean_basis_vector(v, i),
                            B,
                        )
                    end
                end
                @testset "zero argument allocation" begin
                    @test size(ManifoldsBase.allocate_result(N, rand)) == size(p)
                end
            end

            @testset "Zero index manifold" begin
                Mzero = ManifoldsBase.DefaultManifold()
                N = PowerManifold(Mzero, PowerRepr(), 3)
                p = [1.0, 2.0, 3.0]

                @test p[N, 1] == 1.0
                @test zero_vector(N, p) == zero(p)
            end
        end
    end
end


text:
let
    f(x) = x[1]^2
    function g!(x, out)
        out[1] = 2.*x[1]
    end
    function h!(x, out)
        out[1,1] = 2.
    end
    for Optimizer in (AcceleratedGradientDescent,
                      GradientDescent, ConjugateGradient, LBFGS, BFGS,
                      MomentumGradientDescent)

        res = optimize(f, g!, [0.], Optimizer())
        @assert isapprox(res.minimum[1], 0.)
    end

    for Optimizer in (Newton, NewtonTrustRegion)
        res = optimize(f, g!, h!, [0.], Optimizer())
        @assert isapprox(res.minimum[1], 0.)
    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -12.331138821460632 <= q <= 11.287379993959185)

                     add_NL_constraint(m, :(softplus(0.16893740747119512 + 0.361142465903193*softplus(0.3759145594542668 + 0.23452173305395663*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.21914525995597334*softplus(-0.017859443401365915 + -0.7776638835685663*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.7287835676612944*softplus(0.6615910141517634 + 0.18731902906517428*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + softplus(-0.9123133183860701 + 0.273513749523993*softplus(0.3759145594542668 + 0.23452173305395663*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.052599029592083735*softplus(-0.017859443401365915 + -0.7776638835685663*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.2642412570560455*softplus(0.6615910141517634 + 0.18731902906517428*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + softplus(0.5514504616591211 + 0.27732859795053555*softplus(0.3759145594542668 + 0.23452173305395663*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + -0.08702067567597149*softplus(-0.017859443401365915 + -0.7776638835685663*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.40713382739201665*softplus(0.6615910141517634 + 0.18731902906517428*softplus(-0.8645963938709578 + -0.1467408759040536*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*softplus(-0.9609056399980664 + 0.06060702618594416*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*softplus(-0.9973666264873402 + 0.1126936245130743*softplus(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*softplus(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*softplus(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
###
#   transforms concs and dirality
###
import NumericExtensions.negate!, NumericExtensions.divide!, NumericExtensions.add!, NumericExtensions.exp!

function flow(Pot :: Matrix)
    d1, d2 = size(Pot)

    Flow = Array(Number8, d1, d2)

    for j in 1:d2
        for i in 1:d1

            west  = j == 1  ? d2 : j-1
            east  = j == d2 ? 1  : j+1
            north = i == d1 ? 1  : i+1
            south = i == 1  ? d1 : i-1

            ###
            # Calculate the probability of flow out of a cell based on the pot differnece

            p   = Pot[i    ,   j]    ;
            psw = Pot[south,west] - p;
            psj = Pot[south,j   ] - p;
            pse = Pot[south,east] - p;
            pie = Pot[i    ,east] - p;
            pne = Pot[north,east] - p;
            pnj = Pot[north,j   ] - p;
            pnw = Pot[north,west] - p;
            piw = Pot[i    ,west] - p;

            ge = [psw, psj, pse, pie, pne, pnj, pnw, piw]
            ge2 = [psw, psj, pse, pie, pne, pnj, pnw, piw]
            ge = negate!(divide!(ge2, add!(negate!(exp!(ge)), 1))) #

            ###
            # Remove NaN
            ###

            ge_wo_na = Number8([isnan(x) ? 1.0 / 8.0 : x/8.0 for x in ge])

            ###
            # Normalize by eight.
            ###

            Flow[i,j] = ge_wo_na
        end
    end
    return Flow
end

function diffusionJl!(Conc :: Matrix, pot :: Matrix, p_move :: Matrix) #concentration and direction
    d1, d2 = size(Conc)

    Flow = flow(pot)

    for j in 1:d2
        for i in 1:d1

            west  = j == 1  ? d2 : j-1
            east  = j == d2 ? 1  : j+1
            north = i == d1 ? 1  : i+1
            south = i == 1  ? d1 : i-1

            ###
            # Get the flow out of cell ij
            ###

            outflow =  Conc[i,j] * sum(Flow[i,j])

            ###
            # Calculate the inflow based on the outflow from other cells into this on.
            ###

            inflow  =   Conc[south, west] * Flow[south, west].s4 +
                        Conc[south, j   ] * Flow[south, j   ].s5 +
                        Conc[south, east] * Flow[south, east].s6 +
                        Conc[i    , east] * Flow[i    , east].s7 +
                        Conc[north, east] * Flow[north, east].s0 +
                        Conc[north, j   ] * Flow[north, j   ].s1 +
                        Conc[north, west] * Flow[north, west].s2 +
                        Conc[i    , west] * Flow[i    , west].s3 ;

            ###
            # Inflow - outflow = change
            ###

            p_move[i,j] = inflow - outflow
        end
    end
end


text:
using DeduplicateFiles
using Base.Test

# A temporary directory to do tests in
const testdir = mktempdir();

# Some filepaths to use in tests
const testpaths = joinpath.([testdir],["dir1" "dir2"],
    ["file1", "file2", "file3", "file4", "file5", "file6", "file7", "file8"])

const dir1 = dirname(testpaths[1,1]);
const dir2 = dirname(testpaths[1,2]);

# Utility function that creates a file containing a string.
function strfile(str, pth)
  open(pth,"w") do f
    write(f,str)
  end
end

function testproc(list)
  return [list]
end

try
  # Generate a few directories, files and symlinks to test deduplication on
  mkdir(dir1)
  mkdir(dir2)
  strfile(b"dup 1", testpaths[1,1])
  strfile(b"dup 1", testpaths[1,2])
  strfile(b"dup 1", testpaths[2,1])
  strfile(b"nodup", testpaths[2,2])
  strfile(b"dup two", testpaths[3,1])
  strfile(b"dup two", testpaths[3,2])
  strfile(b"unique", testpaths[4,1])
  strfile(b"unique2", testpaths[4,2])
  strfile(b"MAHAVFT (CRC32 collision)", testpaths[5,1])
  strfile(b"VJM (CRC32 collision)", testpaths[5,2])
  strfile(b"KSETVMW (CRC32 collision 2)", testpaths[6,1])
  strfile(b"XNGMFOX (CRC32 collision 2)", testpaths[6,2])
  symlink(testpaths[1,2], testpaths[7,1])
  symlink(testpaths[1,2], testpaths[7,2])
  symlink(testpaths[6,2], testpaths[8,1])
  symlink(dirname(testpaths[1,1]), testpaths[8,2])

  @assert(!islink(testpaths[1,1]))
  @assert(islink(testpaths[7,1]))

  #run(`ls -lR $testdir`) # list the test files and links that were just created.

  # tests
  @test(same_file(testpaths[1], testpaths[1]) === true)
  @test(same_file(testpaths[1], testpaths[2]) === false)
  @test(same_file(testpaths[1], testpaths[3]) === false)

  @test(identical_files(testpaths[1], testpaths[2]) === true)
  @test(identical_files(testpaths[1], testpaths[4]) === false)
  @test(identical_files(testpaths[1], testpaths[2], buflen=1) === true)
  @test(identical_files(testpaths[1], testpaths[4], buflen=1) === false)
  @test(identical_files(testpaths[7], testpaths[8], buflen=1) === false)

  @test(delete_duplicate_file(delete=testpaths[1], keep=testpaths[2]) === true)
  @test(!isfile(testpaths[1]))
  strfile("dup 1", testpaths[1]) # re-create deleted file

  @test(delete_duplicate_file(delete=testpaths[1], keep=testpaths[4]) === false)
  @test(isfile(testpaths[1]))

  @test(!islink(testpaths[1]))
  @test(islink(testpaths[15]))

  filelist = collect(indexfiles(testdir))
  @test(length(filelist)==12)
  @test(length(indexfiles(testdir,follow_symlinks=true))==12)

  procresult = DeduplicateFiles.processdups(x->x.stat.size, testproc, filelist)
  #show(procresult)

  @test(length(deduplicate_files([dir1, dir2], (x,y)->(x.start==dir1 && y.start==dir2), dry_run=true, verbose=true))==3)
  @test(isfile(testpaths[1,1]) && isfile(testpaths[1,2]) &&
        isfile(testpaths[2,1]) && isfile(testpaths[2,2]) &&
        isfile(testpaths[3,1]) && isfile(testpaths[3,2]) &&
        isfile(testpaths[4,1]) && isfile(testpaths[4,2]) &&
        isfile(testpaths[5,1]) && isfile(testpaths[5,2]) &&
        isfile(testpaths[6,1]) && isfile(testpaths[6,2]) )

  @test(length(deduplicate_files([dir1, dir2], (x,y)->(x.start==dir1 && y.start==dir2)))==3)
  @test(!isfile(testpaths[1,1]) && isfile(testpaths[1,2]) &&
        !isfile(testpaths[2,1]) && isfile(testpaths[2,2]) &&
        !isfile(testpaths[3,1]) && isfile(testpaths[3,2]) &&
        isfile(testpaths[4,1]) && isfile(testpaths[4,2]) &&
        isfile(testpaths[5,1]) && isfile(testpaths[5,2]) &&
        isfile(testpaths[6,1]) && isfile(testpaths[6,2]) )

  @test(length(deduplicate_files([dir1, dir2], (x,y)->(x.start==dir2 && y.start==dir1)))==0)
  strfile("dup 1", testpaths[1,1])
  strfile("dup 1", testpaths[2,1])
  strfile("dup two", testpaths[3,1])

  lst = deduplicate_files([dir1, dir2], (x,y)->(x.start==dir2 && y.start==dir1))
  @test(length(lst)==2)
  @test(DeduplicateFiles.totalsize(lst)==5+7)
  @test(isfile(testpaths[1,1]) && !isfile(testpaths[1,2]) &&
        isfile(testpaths[2,1]) && isfile(testpaths[2,2]) &&
        isfile(testpaths[3,1]) && !isfile(testpaths[3,2]) &&
        isfile(testpaths[4,1]) && isfile(testpaths[4,2]) &&
        isfile(testpaths[5,1]) && isfile(testpaths[5,2]) &&
        isfile(testpaths[6,1]) && isfile(testpaths[6,2]) )
  strfile("dup 1", testpaths[1,2])
  strfile("dup two", testpaths[3,2])

  # TODO: Try using "mount -t bind" and/or creating hard-linked directories
  #       producing different paths to the same file.

finally
  rm(testdir, force=true, recursive=true)
end


text:
@testset "Combinations" begin

    # length of the construct
    n_rand = rand(2:5)
    #input modules
    mt_String_array = groupmod(["1", "2", "3", "4","5","6","7","8","9"])
    test_com = Combination(mt_String_array.m,n_rand)

    @testset "setup" begin
        @test length(test_com) == binomial(length(mt_String_array),n_rand)
        @test size(test_com) == (binomial(length(mt_String_array),n_rand),1)
        @test length(collect(test_com)) == length(test_com)
        # all unique
        @test length(Set(collect(test_com))) == length(test_com)
    end

    # check if getindex equals the iterator
    # chekc if all combination have unique elements
    @testset "getindex" begin
        for (a,b) in zip(test_com,1:length(test_com))
            @test isequal(a,test_com[b])
            @test allunique(a) == true
        end
    end

    test_collect = collect(test_com)
    @testset "test_collect" begin
        @test length(test_collect) == binomial(9,n_rand)
        @test eltype(test_collect) == UnorderedConstruct{Mod{String}}
    end

    test_com_1 = Combination(mt_String_array.m,1)
    test_len = _len1(test_com_1,[2])
    @testset "test_collect" begin
        @test typeof(test_len[1])== UnorderedConstruct{Mod{String}}
        @test test_len[2] == [3]
    end

end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

include("define_save.jl") #trajLoad

function groupattributes{T<:AbstractString}(files::Vector{T}, getlabel::Function, getx::Function, gety::Function)

  outdict = Dict{ASCIIString, Vector{(Float64, Vector{Float64})}}()

  M = Array(Any, length(files), 3) #number of files by L,x,y

  for (i, f) in enumerate(files)

    d = trajLoad(f)
    M[i, :] = [getlabel(d), getx(d), gety(d)] #each row is L, x, y
  end

  for (L, L_inds) in groupbycol(M, 1)

    Mx = M[L_inds, :]

    xyvecs = map(groupbycol(Mx, 2)) do tup #vector of (x, yvec)
      x, xinds = tup

      return (Float64(x), Float64(Mx[xinds, 3])) #convert to float for plotting
    end

    sort!(xyvecs, by = v -> v[1]) #sort the x's of easier plotting

    outdict[L] = xyvecs
  end

  return outdict
end

function groupbycol(M::Array{Any, 2}, col::Int64)
  #Looks at all the entries in column 'col' of M.  For each unique entry, return a list of row indices

  labels = unique(M[:, col])
  inds = map(l -> find(x -> x == l, M[:, col]), labels)

  return zip(labels, inds)
end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using Expect
using Test

# like @test, but stops execution of the testset
macro require(expr)
    quote
        value = $(esc(expr))
        if value
            res = Test.Pass(:test, nothing, nothing, nothing)
            Test.record(Test.get_testset(), res)
        else
            throw(Main.Base.AssertionError(string($(Expr(:inert, expr)))))
        end
    end
end


text:
############################################################################################
# Import External PackagesJK
using Test
using Random: Random, AbstractRNG, seed!
using Statistics

############################################################################################
# Import Baytes Packages
using BaytesCore

import BaytesCore:
    get_result,
    result!,
    ResamplingMethod

############################################################################################
# Include Files
include("TestHelper.jl")

############################################################################################
# Run Tests
@testset "All tests" begin
    include("test-Algorithm.jl")
    include("test-Buffer.jl")
    include("test-Core.jl")
    include("test-Sampling.jl")
    include("test-Tempering.jl")
end


text:
## Constants
go = 9.80665;
d2r = pi/180;
r2d = 180/pi;

REARTH = 6370987.308;         # Mean earth radius - m
WEII3 = 7.292115e-5;          # Angular rotation of earth - rad/s
GO = 9.80665;                 # Gravity Constant
AGRAV = 9.80675445;           # Standard value of gravity acceleration - m/s^2
GCONST = 6.673e-11;                # Universal gravitational constant - Nm^2/kg^2
EARTH_MASS = 5.973332e24;     # Mass of the earth - kg
MU = 3.9860044e14;            # Gravitational parameter=G*EARTH_MASS - m^3/s^2
C20 = -4.8416685e-4;          # Second degree zonal gravitational coefficient - ND
FLATTENING = 3.33528106e-3;   # Flattening of the Earth (WGS84) - ND
SMAJOR_AXIS = 6378137.;       # Semi-major axis of Earth's ellipsoid (WGS84) - m
GW_CLONG = 0.0;               # Greenwich celestial longitude at start of simulation - rad
RGAS = 287.053;               # Ideal gas constant - J/(K*kg)=N*m/(K*kg)
KBOLTZ = 1.38e-23;            # Boltzmann's constant - Ws/K
R2D = 57.29577951309314;      # Radians to Degrees
D2R = 0.01745329251994;       # Degrees to Radiancs

#numerical constants
PI = 3.141592653589793;       # Circumference of unit diameter circle
EPS = 1.e-10;                 # Machine precision error (type double)
SMALL = 1e-7;                 # Small real number
ILARGE = 9999;                # Large integer number
LARGE = 1e10;                 # Large real number (type double)

#conversion factors
M2F = 3.280834;               # Conversion factor m->ft
M2NMI = 5.399568e-4;          # Conversion factor  m->nm

# WGS84 Constants
J2 = -sqrt(5)*C20;
aa = 6378137.0;
finv = 298.257223563;
ff = 1.0/finv;
bb = aa*( 1.0 - ff );
e2 = 1.0 - ( 1.0 - ff )^2;

wgs_aa = aa;
wgs_finv = finv;
wgs_ff = ff;
wgs_bb = bb;
wgs_e2 = e2;

###########################################################################
## Transformation about Axis 1
function TR1(x)
  Tr1 = [
         1.0    0.0      0.0
         0.0    cos(x) sin(x)
         0.0   -sin(x) cos(x)
  ];

  return Tr1;
end

###########################################################################
## Transformation about Axis 2
function TR2(x)
  Tr2 = [
         cos(x)    0.0      -sin(x)
         0.0       1.0       0.0
         sin(x)    0.0       cos(x)
  ];

  return Tr2;
end

###########################################################################
## Transformation about Axis 3
function TR3(x)
  Tr3 = [
         cos(x)    sin(x)     0.0
         -sin(x)   cos(x)     0.0
         0.0       0.0        1.0
  ];

  return Tr3;
end

###########################################################################
## Transformation from NED to Body
function TR_BN(roll, pitch, heading)

  TRBN = TR1(roll)*TR2(pitch)*TR3(heading);
  return TRBN;

end

###########################################################################
## Transformation from ECEF to NED
function TR_NE(latitude, longitude)
  clat = cos(latitude);
  slat = sin(latitude);
  clon = cos(longitude);
  slon = sin(longitude);

  TRNE = [
          -slat*clon   -slat*slon   clat
          -slon         clon        0.0
          -clat*clon   -clat*slon   -slat
  ];
  return TRNE;
end

###########################################################################
## Transformation from ECI to ECEF
function TR_EI(time)
  TREI = TR3(WEII3*time);
  return TREI;
end

###########################################################################
## Transformation from Body to Wind Skid To Turn
function TR_WB(alpha, beta)
  TRWB = TR3(beta)*TR2(-alpha);
  return TRWB;
end

###########################################################################
## Transformation from Body to Wind Skid To Turn
function TR_VN(gam, azi)
    TRVN = zeros(3,3);
    TRVN[1,1] = cos(gam)*cos(azi);
    TRVN[1,2] = cos(gam)*sin(azi);
    TRVN[1,3] = -sin(gam);
    TRVN[2,1] = -sin(azi);
    TRVN[2,2] = cos(azi);
    TRVN[2,3] = 0.0;
    TRVN[3,1] = sin(gam)*cos(azi);
    TRVN[3,2] = sin(gam)*sin(azi);
    TRVN[3,3] = cos(gam);

  return TRVN;
end

###########################################################################
## Transformation from Body to Wind Skid To Turn
function TR_BV( alfa,  phii)
  ca = cos(alfa);
  sa = sin(alfa);
  cp = cos(phii);
  sp = sin(phii);

  TBV = zeros(3,3);
  TBV[1,1] = ca;
  TBV[1,2] = sa*sp;
  TBV[1,3] = -sa*cp;
  TBV[2,1] = 0.0;
  TBV[2,2] = cp;
  TBV[2,3] = sp;
  TBV[3,1] = sa;
  TBV[3,2] = -ca*sp;
  TBV[3,3] = cp*ca;

  return TBV;
end
###########################################################################
# Atmosphere Function from Zipfel's book
function Atmosphere(alt, vmag)

  # Universal Air Constant
  Rair = 287.058;

  # Atmosphere Parameters
  if( alt < 11000.0 )
    Temp = 288.15 - 0.0065*alt;
    Press = 101325.0*(Temp/288.15)^5.2559;
  else
    Temp = 216.0;
    Press = 22630.0*exp(-0.00015769*(alt - 11000.0));
  end
  rho = Press/(Rair*Temp);
  aa = sqrt(1.4*Rair*Temp);

  # Mach and Dynamic Pressure
  Mach = vmag/aa;
  qbar = 0.5*rho*vmag*vmag;

  return [Press, Temp, rho, aa, Mach, qbar];

end

###########################################################################
function ECEF2LLH(rm_ecef)

  X = rm_ecef[1];
  Y = rm_ecef[2];
  Z = rm_ecef[3];

  # Longitude
  lon = atan(Y,X);

  # Eccentricity
  ee = e2*(aa/bb)^2;

  p = sqrt( X.*X + Y.*Y);
  r = sqrt( p.*p + Z.*Z);
  u = atan( bb.*Z.*( 1 + ee.*bb./r)./(aa.*p));

  # Latitude
  lat = atan( (Z + ee.*bb.*sin(u).^3)./( p - e2.*aa.*cos(u).^3) );

  v = aa./sqrt( 1.0 - e2.*sin(lat).^2 );

  # Altitude
  alt = p.*cos(lat) + Z.*sin(lat) - aa*aa./v;

  return [lat lon alt];
end

###########################################################################
function LLH2ECEF(lat, lon, alt)

  vv = aa/sqrt( 1.0 - e2*sin(lat)*sin(lat) );
  xx = ( vv + alt )*cos(lat)*cos(lon);
  yy = ( vv + alt )*cos(lat)*sin(lon);
  zz = ( vv*( 1.0 - e2 ) + alt)*sin(lat);

  rm_ecef = [xx;yy;zz];
  return rm_ecef;
end

###########################################################################
function WGS84_GRAVITY(rm_ecef)
  rmag = norm(rm_ecef);
  px = rm_ecef[1];
  py = rm_ecef[2];
  pz = rm_ecef[3];
  tmp1 = pz/rmag;

  vec = [
          ( 1.0 + 1.5*J2*(aa/rmag)^2*(1 - 5*tmp1^2) )*px/rmag
          ( 1.0 + 1.5*J2*(aa/rmag)^2*(1 - 5*tmp1^2) )*py/rmag
          ( 1.0 + 1.5*J2*(aa/rmag)^2*(3 - 5*tmp1^2) )*pz/rmag
  ];

  # Raw Gravity Vector
  GRAV_ECEF = -MU/( rmag^2 )*vec

  # Earth Rotation Vector
  wei_ecef = [0.0;0.0;WEII3];

  # Effective Gravity Vector
  grav_ecef = GRAV_ECEF - cross( wei_ecef , cross( wei_ecef, rm_ecef ) );

  return grav_ecef;
end

###########################################################################
function QuatInit(roll, pitch, heading)
  q0 = cos(heading/2)*cos(pitch/2)*cos(roll/2) + sin(heading/2)*sin(pitch/2)*sin(roll/2);
  q1 = cos(heading/2)*cos(pitch/2)*sin(roll/2) - sin(heading/2)*sin(pitch/2)*cos(roll/2);
  q2 = cos(heading/2)*sin(pitch/2)*cos(roll/2) + sin(heading/2)*cos(pitch/2)*sin(roll/2);
  q3 = sin(heading/2)*cos(pitch/2)*cos(roll/2) - cos(heading/2)*sin(pitch/2)*sin(roll/2);

  quat = [q0;q1;q2;q3];

  # Normalization
  quat = quat/norm(quat);

  return quat;
end
###########################################################################
function QuatInit( Tr_bi::Array{Float64,2} )
  t11 = Tr_bi[1,1];
  t12 = Tr_bi[1,2];
  t13 = Tr_bi[1,3];
  t23 = Tr_bi[2,3];
  t33 = Tr_bi[3,3];

  pitch = asin(-t13);
  heading = atan(t12,t11);
  roll = atan(t23,t33);

  q0 = cos(heading/2)*cos(pitch/2)*cos(roll/2) + sin(heading/2)*sin(pitch/2)*sin(roll/2);
  q1 = cos(heading/2)*cos(pitch/2)*sin(roll/2) - sin(heading/2)*sin(pitch/2)*cos(roll/2);
  q2 = cos(heading/2)*sin(pitch/2)*cos(roll/2) + sin(heading/2)*cos(pitch/2)*sin(roll/2);
  q3 = sin(heading/2)*cos(pitch/2)*cos(roll/2) - cos(heading/2)*sin(pitch/2)*sin(roll/2);

  quat = [q0;q1;q2;q3];

  # Normalization
  quat = quat/norm(quat);

  return quat;
end

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

function QuatDerivative(quat::Array{Float64,1}, wb_body::Array{Float64,1} )
  pp = wb_body[1];
  qq = wb_body[2];
  rr = wb_body[3];

  lam = 1.0 - norm(quat)^2;

  BB = [
        0   -pp   -qq    -rr
        pp   0     rr    -qq
        qq  -rr    0      pp
        rr   qq   -pp     0
  ];

  kgain = 100.0;
  quatdot = 0.5*BB*quat + kgain*lam*quat;

  return quatdot;
end

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

function QuatToDCM(quat)

  q0 = quat[1];
  q1 = quat[2];
  q2 = quat[3];
  q3 = quat[4];

  Tr_bi = zeros(3,3);
  Tr_bi[1,1] = q0^2 + q1^2 - q2^2 - q3^2;
  Tr_bi[2,2] = q0^2 - q1^2 + q2^2 - q3^2;
  Tr_bi[3,3] = q0^2 - q1^2 - q2^2 + q3^2;

  Tr_bi[2,1] = 2.0*(q1*q2 - q0*q3);
  Tr_bi[1,2] = 2.0*(q1*q2 + q0*q3);

  Tr_bi[3,1] = 2.0*(q1*q3 + q0*q2);
  Tr_bi[1,3] = 2.0*(q1*q3 - q0*q2);

  Tr_bi[3,2] = 2.0*(q2*q3 - q0*q1);
  Tr_bi[2,3] = 2.0*(q2*q3 + q0*q1);

  return Tr_bi';
end

###########################################################################
function ComputeEuler( Tr_bi::Array{Float64,2} )

  t11 = Tr_bi[1,1];
  t12 = Tr_bi[1,2];
  t13 = Tr_bi[1,3];
  t23 = Tr_bi[2,3];
  t33 = Tr_bi[3,3];

  pitch = asin(-t13);
  heading = atan(t12,t11);
  roll = atan(t23,t33);

  euler = [roll; pitch; heading];

  return euler;
end

###########################################################################
function Pause()
  print("Program Paused:  Press any key to Continue, or press q to quit: ");
  xx = readline();
  if xx == "q\r\n"
    error("Terminated Program");
  end
  return
end

###########################################################################
function CloseAll()

  for i = 1:500
    close();
  end

  return
end

###########################################################################
function Bound(xx::Float64, xmin::Float64, xmax::Float64)
  yy = xx;

  if (xx > xmax)
    yy = xmax;
  elseif(xx < xmin)
    yy = xmin;
  else
    yy = xx;
  end

  return yy;
end


text:
module MiniFastai

using Flux
using Flux.Data: MNIST, DataLoader
using Flux: onehotbatch, onecold, logitcrossentropy
using Flux.Optimise: update!, Descent
using Parameters: @with_kw
using Statistics: mean, std
using Zygote: pullback, @adjoint
using CUDA

include("data.jl")
include("callback.jl")
include("cnn.jl")
include("learner.jl")
include("hook.jl")
include("training.jl")
include("utils.jl")

export my_f

export get_processed_data_MNIST
export AbstractDatabunch
export Databunch
export normalize
export normalize_imgs

export AbstractCallback
export TrainEvalCallback
export AvgStatsCallback

export get_model

export AbstractLearner
export Learner

export get_model_hooks

export accuracy

export fit!

end


text:
#=
Define a discrete state space for the car MDP problem
=#

function n_states(mdp::CarMDP)
    n_ego = n_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    routes = get_car_routes(mdp.env)
    n_cars = 0
    for route in routes
        n_cars += n_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
    end
    return n_ego*(n_cars + 1) # do not forget absent state
end

function POMDPs.states(mdp::CarMDP)
    state_space = CarMDPState[]
    ego_states = get_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    routes = get_car_routes(mdp.env)
    for route in routes
        for ego in ego_states
            car_states = get_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
            for car in car_states
                crash =  is_colliding(Vehicle(ego, mdp.ego_type, 1), Vehicle(car, mdp.car_type, 2))
                # enumerate all possible routes
                lane = get_lane(mdp.env.roadway, car)
                push!(state_space, CarMDPState(crash, ego, car, SVector{2, LaneTag}(route[1], route[end])))
            end
        end
    end
    for ego in ego_states
        # add absent states
        push!(state_space, CarMDPState(false, ego, get_off_the_grid(mdp), SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0, 0))))
    end
    return state_space
end

function POMDPs.stateindex(mdp::CarMDP, s::CarMDPState)
    n_ego = n_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    routes = get_car_routes(mdp.env)
    # step 1: find ego_index
    ego_i = ego_stateindex(mdp.env, s.ego, mdp.pos_res, mdp.vel_res)
    # handle off the grid case
    if s.car.posG == mdp.off_grid
        si = 0
        for route in routes
            si += n_ego * n_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
        end
        si += ego_i
    else
        n_car = n_car_states(mdp.env, find_route(mdp.env, s.route), mdp.pos_res, mdp.vel_res)
        # step 2: find route_index
        route_i = 0
        for (i, route) in enumerate(routes)
            if s.route[end] == route[end] && s.route[1] == route[1]
                route_i = i
            end
        end
        # step 3: find car_index in car states
        car_i = car_stateindex(mdp.env, s.car, find_route(mdp.env, s.route), mdp.pos_res, mdp.vel_res)
        # linear/cartesian magic
        si = LinearIndices((n_car, n_ego))[car_i, ego_i]

        for i=2:route_i
            size_r = n_ego * n_car_states(mdp.env, routes[i-1], mdp.pos_res, mdp.vel_res)
            si += size_r
        end
    end
    return si
end

function POMDPs.initialstate(mdp::CarMDP)
    ego = initial_ego_state(mdp)
    init_car_states, init_car_routes = initial_car_state_distribution(mdp)
    init_states = Vector{CarMDPState}(length(init_car_states))
    for i=1:length(init_car_states)
        car = init_car_states[i]
        route = init_car_routes[i]
        crash = is_colliding(Vehicle(ego, mdp.ego_type, 1), Vehicle(car, mdp.car_type, 2))
        init_states[i] = CarMDPState(crash, ego, car, SVector{2, LaneTag}(route[1], route[end]))
    end
    # uniform
    probs = ones(length(init_car_states))
    normalize!(probs, 1)
    return SparseCat(init_states, probs)
end

function initial_ego_state(mdp::CarMDP)
    lanes = get_ego_route(mdp.env)
    posF = Frenet(mdp.env.roadway[lanes[1]], mdp.ego_start)
    v0 = 0.
    return VehicleState(posF, mdp.env.roadway, v0)
end

function initial_car_state_distribution(mdp::CarMDP)
    routes = get_car_routes(mdp.env)
    init_car_routes = []
    init_car_states = VehicleState[]
    for route in routes
        car_states = get_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
        for cs in car_states
            push!(init_car_states, cs)
            push!(init_car_routes, route)
        end
    end
    push!(init_car_states, get_off_the_grid(mdp))
    push!(init_car_routes, SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0,0)))
    return init_car_states, init_car_routes
end

function initial_car_state(mdp::CarMDP, rng::AbstractRNG)
    init_car_states, init_car_routes = initial_car_state_distribution(mdp)
    return rand(rng, init_car_states), rand(rng, init_car_routes)
end

function car_starting_states(mdp::CarMDP, min_speed::Float64 = 6.0)
    # list of car starting states
    routes = get_car_routes(mdp.env)
    v_space = min_speed:mdp.vel_res:mdp.env.params.speed_limit
    N_states = 0
    for route in routes
        N_states += length(v_space)
    end
    car_states = Vector{VehicleState}(N_states)
    start_routes = Vector{SVector{2, LaneTag}}(N_states)
    i = 1
    for route in routes
        lane = route[1]
        for v in v_space
            car_states[i] = VehicleState(Frenet(mdp.env.roadway[lane], 0.), mdp.env.roadway, v)
            start_routes[i] = SVector{2, LaneTag}(route[1], route[end])
            i += 1
        end
    end
    return car_states, start_routes
end


text:
# TODO
#	- Implment broadcasting
#
#	- Consider generalizing this to static bit arrays of any dimension.
#	- Also, consider incorporating this into the BitArray framework. This would
#		require a reworking BitArray; possibly parameterize by backing type?
#		e.g.	BitArray{N} = AbstractBitArray{N, Array{N,UInt64}}
#				SBitArray{S,N} = AbstractBitArray{N, SArray{S,N,UInt64}}
#				MBitArray{S,N} = AbstractBitArray{N, MArray{S,N,UInt64}}
#		then SBitVector{L} = SBitArray{Tuple{L}, 1}, etc.
"""
	BitStrings (module)

Dense-packed, stack-allocated bit strings.
(Basically a merging of SVector and BitVector.)
"""
module BitStrings

#using BaseExtensions: allequal
using StaticArrays
import Base: convert, promote_type, promote_rule
import Base: string, print, show, display
import Base: length, size, checkbounds, ndims, axes, getindex, setindex, setindex!, iterate
import Base: vcat
import Base: map, map!		#, bit_map!
import Base.Broadcast: broadcastable, BroadcastStyle, Broadcasted, broadcasted
import Base: count, sum
import Base: +, -, *, /
import LinearAlgebra: dot

export AbstractBitString, SBitString, MBitString, parity, hammdist, setindex

# TODO:	Support for other binary element types
# TODO:	Support multidimensional arrays?

## From BaseExtensions

const Iterable = Union{Tuple, AbstractArray, UnitRange, Base.Generator}

function allequal(it::Iterable)
	if length(it) <= 1
		return true
	end
	v = first(it)
	for v_ in it
		v_ === v || return false
	end
	return true
end

## OK, Now the actual stuff

abstract type AbstractBitString{L,N} <: AbstractVector{Bool} end

struct SBitString{L,N} <: AbstractBitString{L,N}
	chunks::SVector{N,UInt64}
end

struct MBitString{L,N} <: AbstractBitString{L,N}
	chunks::MVector{N,UInt64}
	MBitString{L,N}(ch) where {L,N} = new{L,N}(ch)
	MBitString{L,N}(::UndefInitializer) where {L,N} = new{L,N}(MVector{N,UInt64}(undef))
end

# These determine the appropriate constructor type for operations
promote_type(::Type{<:SBitString}) = SBitString
promote_type(::Type{<:MBitString}) = MBitString
# These are needed to strip the {L,N} parameters
promote_type(::Type{T}, ::Type{T}) where {T<:SBitString{L,N}} where {L,N} = SBitString
promote_type(::Type{T}, ::Type{T}) where {T<:MBitString{L,N}} where {L,N} = MBitString
promote_rule(::Type{<:SBitString}, ::Type{<:MBitString}) = MBitString

SBitString(b::Bool) = SBitString{1,1}((UInt64(b),))
MBitString(b::Bool) = MBitString{1,1}((UInt64(b),))

convert(::Type{SBitString}, b::Bool) = SBitString(b)
convert(::Type{MBitString}, b::Bool) = MBitString(b)

# AbstractBitString(bs::AbstractBitString{L,N}) where {L,N} = AbstractBitString{L,N}(bs.chunks)
#
#
# # Construction from integer
# AbstractBitString(i::Integer) = AbstractBitString{8*sizeof(typeof(i)),1}(SVector{1}(UInt64(i)))
#

# Construction from BitVector (mre efficient than for general iterables)
function (::Type{T})(b::BitVector) where T<:AbstractBitString
	T{length(b), length(b.chunks)}(b.chunks)
end

# Construction from any iterable
function (::Type{T})(A) where T<:AbstractBitString
	chunks = compute_chunks(A)
	T{length(A), length(chunks)}(chunks)
end

convert(::Type{SBitString{L,N}}, v::SVector{N,UInt64}) where {L,N} = SBitString{L,N}(v)
convert(::Type{MBitString{L,N}}, v::MVector{N,UInt64}) where {L,N} = MBitString{L,N}(v)

convert(::Type{T}, bs::T) where {T<:AbstractBitString} = bs
convert(::Type{SBitString}, bs::AbstractBitString{L,N}) where {L,N} = SBitString{L,N}(bs.chunks)
convert(::Type{SBitString{L,N}}, bs::AbstractBitString{L,N}) where {L,N} = SBitString{L,N}(bs.chunks)
convert(::Type{MBitString}, bs::AbstractBitString{L,N}) where {L,N} = MBitString{L,N}(bs.chunks)
convert(::Type{MBitString{L,N}}, bs::AbstractBitString{L,N}) where {L,N} = MBitString{L,N}(bs.chunks)

nchunks(len::Int) = Base._div64(len + 63)
nchunks(b::AbstractBitString{L,N}) where {L,N} = N

# Compute chunks from an iterable
function compute_chunks(A)
	len = length(A)
	nch = nchunks(len)

	# empty vector
	nch == 0 && return MVector{0,UInt64}()

	chunks = MVector{nch,UInt64}(undef)
	itr = iterate(A)
	@inbounds begin
		for ich = 1:nch-1
			ch = UInt64(0)
			for j = 0:63
				ch |= (UInt64(convert(Bool, itr[1])) << j)
				itr = iterate(A, itr[2])
			end
			chunks[ich] = ch
		end
		ch = UInt64(0)
		for j = 0:Base._mod64(len-1)
			ch |= (UInt64(convert(Bool, itr[1])) << j)
			itr = iterate(A, itr[2])
		end
		chunks[nch] = ch
	end
	return chunks
end

length(b::AbstractBitString{L}) where {L} = L
size(b::AbstractBitString{L}) where {L} = (L,)
ndims(b::AbstractBitString)  = 1
ndims(::Type{T}) where {T<:AbstractBitString} = 1
axes(b::AbstractBitString{L}) where {L} = (Base.OneTo(L),)

# # Conversion to BitVector
# convert(::Type{BitVector}, bs::AbstractBitString) = BitVector(bs)
# convert(::Type{AbstractBitString}, bv::BitVector) = AbstractBitString(bv)
#
# function BitVector(bs::AbstractBitString)
# 	bv = BitVector(undef, length(bs))
# 	bv.chunks = Vector(bs.chunks)
# 	bv
# end

"""
	string(bs::AbstractBitString)

Create a string representation of a `AbstractBitString`.
Bits are shown as 0,1 in order from left to right, with an underscore every 8 bits.
"""
function string(bs::AbstractBitString)
	mv = MVector{length(bs) + max(0,(length(bs)-1))>>3, Char}(undef)
	j = 1
	for i in 1:length(bs)
		mv[j] = bs[i] + '0'
		j += 1
		if mod(i, 8) == 0 && i < length(bs)
			mv[j] = '_'
			j += 1
		end
	end
	string(mv...)
end

display(bs::AbstractBitString) = show(bs)
show(io::IO, bs::SBitString) = print(io, string(bs))
show(io::IO, bs::MBitString) = print(io, "(", string(bs), ")")
#printstyled(io, string((bs.bvec .+ '0')...); bold=true)

##  Indexing

checkbounds(b::AbstractBitString{L}, I...) where {L} = Base.checkbounds_indices(Bool, (Base.OneTo(L),), I) || Base.throw_boundserror(b, I)

@inline function get_chunk_index(i::Integer)
	i1, i2 = Base.get_chunks_id(i)
	msk = UInt64(1) << i2
	return (i1, msk)
end

# getindex - entry point
@inline function getindex(b::AbstractBitString, i)
   @boundscheck checkbounds(b, i)
	_getindex(b, i)
end

# implementations (no bounds checking)
@inline _getindex(b::AbstractBitString, i::CartesianIndex{1}) = _getindex(b, i[1])
@inline _getindex(b::AbstractBitString, i::Integer) = _getindex(b, get_chunk_index(i))

# unsafe getindex using precomputed chunk index and mask
@inline function _getindex(b::AbstractBitString, (ich,msk)::Tuple{Int64, UInt64})
	# @info "_getindex (chunk): ich = $ich, msk = $msk"
	@inbounds r = (b.chunks[ich] & msk) != 0
	return r
end

# fallback / index by an iterable
function _getindex(b::T, itr) where T<:AbstractBitString
	 v = [_getindex(b, get_chunk_index(i)) for i in itr]
	 promote_type(T)(v)
end

function setindex!(bs::MBitString, val, i)
	@boundscheck checkbounds(bs, i)
	_setindex!(bs, val, i)
end

@inline function _setindex!(bs::MBitString, val::Bool, i::Integer)
	i1, i2 = Base.get_chunks_id(i)
	msk = ~(UInt64(1) << i2)
	@inbounds bs.chunks[i1] = (bs.chunks[i1] & msk) | (val << i2)
end

function _setindex!(bs::MBitString, val, itr)
	for (iv,ib) in enumerate(itr)
		_setindex!(bs, val[iv], ib)
	end
end

# function setindex(bs::AbstractBitString, val, i)
# 	@boundscheck checkbounds(bs, i)
# 	_setindex(bs, val, i)
# end
#
#
# @inline function _setindex(bs::AbstractBitString{L,N}, val::Bool, i::Integer) where {L,N}
# 	temp = MVector{N,UInt64}(bs.chunks)
# 	i1, i2 = Base.get_chunks_id(i)
# 	msk = ~(UInt64(1) << i2)
# 	@inbounds temp[i1] = (temp[i1] & msk) | (val << i2)
# 	AbstractBitString{L,N}(SVector{N,UInt64}(temp))
# end
#
#
#
# @inline function _setindex(bs::AbstractBitString{L,N}, val::AbstractVector{Bool}, itr) where {L,N}
# 	temp = MVector{N,UInt64}(bs.chunks)
# 	for (iv,ib) in enumerate(itr)
# 		i1, i2 = Base.get_chunks_id(ib)
# 		msk = ~(UInt64(1) << i2)
# 		@inbounds temp[i1] = (temp[i1] & msk) | (val[iv] << i2)
# 	end
# 	AbstractBitString{L,N}(SVector{N,UInt64}(temp))
# end
#
#
#

##  Iteration

function iterate(bs::AbstractBitString, i::Int=0)
    i >= length(bs) && return nothing
    (bs.chunks[Base._div64(i)+1] & (UInt64(1)<<Base._mod64(i)) != 0, i+1)
end

vcat(a::AbstractBitString) = a

@inline function vcat(a::AbstractBitString, b::AbstractBitString)
	lena = length(a)
	lenb = length(b)
	nch = nchunks(lena + lenb)
	na = length(a.chunks)
	nb = length(b.chunks)
	chunks = MVector{nch, UInt64}(undef)
	@inbounds for i = 1:na
		chunks[i] = a.chunks[i]
	end
	off = Base._mod64(length(a))
	@inbounds for j = 1:nb
		chunks[na + j-1]  |= b.chunks[j] << off
		if na + j <= nch
			chunks[na + j] = b.chunks[j] >> (64-off)
		end
	end
	return promote_type(typeof(a), typeof(b)){lena+lenb, nch}(chunks)
end
#
# totallength() = 0
# totallength(a::AbstractBitString) = length(a)
# totallength(a::AbstractBitString, others::AbstractBitString...) = length(a) + totallength(others...)
#
# # Without inlining it is SLOW
# @inline function vcat(args::AbstractBitString...)
# 	totlen = totallength(args...)
# 	nch = nchunks(totlen)
# 	chunks = MVector{nch, UInt64}(undef)
# 	a = args[1]
# 	for i = 1:length(a.chunks)
# 		chunks[i] = a.chunks[i]
# 	end
# 	len = length(a)
#
# 	for a in Base.tail(args)
# 		ich = nchunks(len)
# 		off = Base._mod64(len)
# 		for j = 1:length(a.chunks)
# 			chunks[ich - 1 + j]  |= a.chunks[j] << off
# 			if ich + j <= nch
# 				chunks[ich + j] = a.chunks[j] >> (64-off)
# 			end
# 		end
# 		len += length(a)
# 	end
# 	return AbstractBitString{totlen, nch}(chunks)
# end

# Basic arithmetic operations -- result in a non-dense array
-(b::AbstractBitString) = (-).(b)
*(b::AbstractBitString, x::Number) = b .* x
*(x::Number, b::AbstractBitString) = x .* b
/(b::AbstractBitString, x::Number) = b ./ x

# Efficient versions of map and map! for boolean functions.
map(::Union{typeof(~), typeof(!)}, a::AbstractBitString) = bit_map(~, a)
map(::Union{typeof(&), typeof(*), typeof(min)}, a::AbstractBitString, b::AbstractBitString) = bit_map(&, a, b)
map(::Union{typeof(|), typeof(max)}, a::AbstractBitString, b::AbstractBitString) = bit_map(|, a, b)
map(::Union{typeof(⊻), typeof(!=)}, a::AbstractBitString, b::AbstractBitString) = bit_map(⊻, a, b)
map(::typeof(*), a::AbstractBitString, b::AbstractBitString) = bit_map(*, a, b)
map(::typeof(==), a::AbstractBitString, b::AbstractBitString) = bit_map((x,y) -> ~xor(x,y), a, b)
map(::typeof(^), a::AbstractBitString, b::AbstractBitString) = bit_map((x,y) -> x | ~y, a, b)
map(::typeof(>), a::AbstractBitString, b::AbstractBitString) = bit_map((x,y) -> x & ~y, a, b)
map(::typeof(>=), a::AbstractBitString, b::AbstractBitString) = bit_map((x,y) -> x | ~y, a, b)
map(::typeof(<), a::AbstractBitString, b::AbstractBitString) = bit_map((x,y) -> y & ~x, a, b)
map(::typeof(<=), a::AbstractBitString, b::AbstractBitString) = bit_map((x,y) -> y | ~x, a, b)
map(::typeof(min), args...) where T = bit_map(&, args...)
map(::typeof(max), args...) = bit_map(|, args...)

# 1-ary functions
@inline function bit_map(f::F, a::AbstractBitString{L,N}) where {L,N,F}
	temp = MVector{N,UInt64}(undef)
#	isempty(A) && return AbstractBitString{0,1}()
	for i = 1:N
   	temp[i] = f(a.chunks[i])
   end
   temp[N] &= Base._msk_end(L)
   promote_type(typeof(a))(temp)
end

# 2-ary functions
function bit_map(f::F, a::AbstractBitString{L,N},  b::AbstractBitString{L,N}) where {F,L,N}
	#	isempty(A) && return AbstractBitString{0,1}()
	temp = MVector{N,UInt64}(undef)
	for i = 1:N
   	temp[i] = f(a.chunks[i], b.chunks[i])
   end
   temp[N] &= Base._msk_end(L)
   promote_type(typeof(a), typeof(b)){L,N}(temp)
end

# n-ary functions
@inline function bit_map(f::F, args::AbstractBitString{L,N}...) where {F,L,N}
	allequal(length(arg) for arg in args) || throw(DimensionMismatch("sizes of all arguments must match"))
	temp = MVector{N,UInt64}(undef)
	for i = 1:N
   	temp[i] = f((arg.chunks[i] for arg in args)...)
   end
   temp[N] &= Base._msk_end(L)
   AbstractBitString{L,N}(SVector{N}(temp))
end
#

map!(::Union{typeof(~), typeof(!)}, dest::MBitString, a::AbstractBitString) = bit_map!(~, dest, a)
map!(::Union{typeof(&), typeof(*), typeof(min)}, a::AbstractBitString, b::AbstractBitString) = bit_map!(&, dest, a, b)
map!(::Union{typeof(|), typeof(max)}, a::AbstractBitString, b::AbstractBitString) = bit_map!(|, dest, a, b)
map!(::typeof(⊻), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!(⊻, dest, a, b)
map!(::typeof(*), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!(*, dest,a, b)
map!(::typeof(!=), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!(xor, dest, a, b)
map!(::typeof(==), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!((x,y) -> ~xor(x,y), dest, a, b)
map!(::typeof(^), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!((x,y) -> x | ~y, dest, a, b)
map!(::typeof(>), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!((x,y) -> x & ~y, dest, a, b)
map!(::typeof(>=), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!((x,y) -> x | ~y, dest, a, b)
map!(::typeof(<), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!((x,y) -> y & ~x, dest, a, b)
map!(::typeof(<=), dest::MBitString, a::AbstractBitString, b::AbstractBitString) = bit_map!((x,y) -> y | ~x, dest, a, b)
map!(::typeof(min), dest::MBitString, args::AbstractBitString...) where T = bit_map!(&, dest, args...)
map!(::typeof(max), dest::MBitString, args::AbstractBitString...) = bit_map!(|, dest, args...)

@inline function bit_map!(f, dest, args...)
	checklengths(args...)
	unsafe_bit_map!(f, dest, args...)
end

checklengths(a) = nothing

@inline function checklengths(a, b)
	length(a) == length(b) || throw(DimensionMismatch("sizes of A and B must match"))
	nothing
end

@inline function checklengths(args...)
	allequal(length(arg) for arg in args) || throw(DimensionMismatch("all arguments must have the same length"))
	nothing
end

# 1-ary functions
function unsafe_bit_map!(f::F, dest::MBitString{L,N}, A::AbstractBitString{L,N}) where {F,L,N}
	destc = dest.chunks
	for i = 1:N
   	destc[i] = f(A.chunks[i])
   end
	destc[N] &= Base._msk_end(L)
	dest
end

# 2-ary functions
@inline function unsafe_bit_map!(f::F, dest::MBitString{L,N}, A::AbstractBitString{L,N}, B::AbstractBitString{L,N}) where {L,N,F}
	destc = dest.chunks
	@inbounds for i = 1:N
   	destc[i] = f(A.chunks[i], B.chunks[i])
   end
   @inbounds destc[N] &= Base._msk_end(L)
	dest
end

# n-ary functions
function unsafe_bit_map!(f::F, dest::MBitString{L,N}, args::AbstractBitString{L,N}...) where {L,N,F}
	destc = dest.chunks
	@inbounds for i = 1:N
   	destc[i] = f((arg.chunks[i] for arg in args)...)
   end
	@inbounds destc[N] &= Base._msk_end(L)
	dest
end

# ## Broadcasting
#
# # broadcastable(b::AbstractBitString) = b
# #
# # # associative bitwise functions
# # const BitwiseFun = Union{typeof(&), typeof(|), typeof(xor), typeof(~)}
# #
# #
# # # Create a custom style so that we can dispatch to our custom implementation
# # struct BitStringStyle{L,N} <: Broadcast.AbstractArrayStyle{1} end
# # BroadcastStyle(::Type{AbstractBitString{L,N}}) where {L,N} = BitStringStyle{L,N}()
# # BroadcastStyle(::BitStringStyle, t::Broadcast.DefaultArrayStyle{N}) where {N} = Broadcast.DefaultArrayStyle{max(1,N)}()
# #
# # # fallback -- default to ArrayStyle
# # broadcasted(bc::BitStringStyle, f::F, args...) where {F} = broadcasted(Broadcast.DefaultArrayStyle{1}(), f, args...)
# #
# # broadcasted(bc::BitStringStyle, f::BitwiseFun, args...) = bs_broadcast(f, args...)
# # broadcasted(bc::BitStringStyle, ::typeof(min), args...) = bs_broadcast(&, args...)
# # broadcasted(bc::BitStringStyle, ::typeof(max), args...) = bs_broadcast(|, args...)
# # broadcasted(bc::BitStringStyle, ::typeof(!=), a, b) = bs_broadcast(xor, a, b)
# # broadcasted(bc::BitStringStyle, ::typeof(==), a, b) = bs_broadcast((x,y) -> ~xor(x,y), a, b)
# # broadcasted(bc::BitStringStyle, ::typeof(^),  a, b) = bs_broadcast((x,y) -> x | ~y, a, b)
# # broadcasted(bc::BitStringStyle, ::typeof(>), a, b) = bs_broadcast((x,y) -> x & ~y, a, b)
# # broadcasted(bc::BitStringStyle, ::typeof(>=), a, b) = bs_broadcast((x,y) -> x | ~y, a, b)
# # broadcasted(bc::BitStringStyle, ::typeof(<), a, b) = bs_broadcast((x,y) -> y & ~x, a, b)
# # broadcasted(bc::BitStringStyle, ::typeof(<=), a, b) = bs_broadcast((x,y) -> y | ~x, a, b)
# #
# # # Custom implementation.  For N=1, should take about 2ns (2μs per 1000)
# # function bs_broadcast(f::F, a::AbstractBitString{L,N}, b::AbstractBitString{L,N}) where {F,L,N}
# # 	mv = MVector{N,UInt64}(undef)
# # 	for i in 1:N
# # 		mv[i] = f(a.chunks[i], b.chunks[i])
# # 	end
# # 	AbstractBitString{L,N}(Tuple(mv))
# # end

## Linear algebra

@inline function dot(x::AbstractBitString, y::AbstractBitString)
    # simplest way to mimic Array dot behavior
    length(x) == length(y) || throw(DimensionMismatch())
    s = 0
    xc = x.chunks
    yc = y.chunks
    @inbounds for i = 1:length(xc)
        s += count_ones(xc[i] & yc[i])
    end
    s
end

# Just like dot, bit with xor instead of &
function hammdist(a::AbstractBitString, b::AbstractBitString)
	# simplest way to mimic Array dot behavior
	length(x) == length(y) || throw(DimensionMismatch())
	s = 0
	xc = x.chunks
	yc = y.chunks
	@inbounds for i = 1:length(xc)
		 s += count_ones(xc[i] ⊻ yc[i])
	end
	s
end

sum(b::AbstractBitString) = count(b)

function count(b::AbstractBitString)
	s = 0
	chk = b.chunks
	@inbounds for i = 1:length(chk)
		 s += count_ones(chk[i])
	end
	s
end

function parity(b::AbstractBitString)
	s = false
	chk = b.chunks
	@inbounds for i = 1:length(chk)
		 s ⊻= isodd(count_ones(chk[i]))
	end
	s
end

end	# module


text:
# Example of the circuit example of Figure 2 of the Macia & Sole (2009) paper.
#include("../../CGP.jl/src/CGP.jl")
using Main.CGP
#include("../../information_theory/src/entropy.jl")
using Combinatorics
using Test

# According to Macia & Sole's definition of Z as the number of interacting units,
#  you would think that Z = 5.  However, their Figure 2 example clearly states
#  that X consists of gates 1, 2, and 3 (which are int3, int4, and int5 below).
#  Thus, Z = 3 even though gates 4 and 5 are computation nodes.
# As of March 13, I am using n from Tononi rather than Z from Macia & Sole
# This means results for degeneracy and complexity may not match Macia and Sole

# Runs without errors on 4/17/21
function test_macia_circuit()
  p0 = Parameters( numinputs=2, numoutputs=2, numinteriors=5, numlevelsback=8 )
  default_funcs( p0.numinputs )  # need to call default_funcs() to define global vars ONE and Ones
  # Circuit of Figure 2 of the Macia & Sole paper
  println("Figure 2 circuit")
  Ones = Main.CGP.construct_ones(p0.numinputs)[p0.numinputs]
  One() = Ones
  funcs = [ Func(Main.CGP.Nand, 2, "NAND") ]
  in1 = InputNode(1)
  in2 = InputNode(2)
  int3 = InteriorNode(NAND, [1, 2])
  int4 = InteriorNode(NAND, [3, 2])
  int5 = InteriorNode(NAND, [3, 2])
  int6 = InteriorNode(NAND, [4, 1])
  int7 = InteriorNode(NAND, [2, 5])
  out8 = OutputNode(6)
  out9 = OutputNode(7)
  c0 = Chromosome(p0,[in1,in2],[int3,int4,int5,int6,int7],[out8,out9],0.0,0.0)
  context = construct_contexts(p0.numinputs)[p0.numinputs]
  execute_chromosome(c0,context)
  println("Macia Sole Figure 2 chromosome: ")
  print_chromosome(c0)
  println("num_mutate_locations: ",num_mutate_locations(c0,funcs))
  @assert number_active(c0) == number_active_old(c0)
  nv = node_values(c0)
  println("nv: ",nv)

  # Node labels are from Macia Figure 2
  I1 = context[1]
  I2 = context[2]
  G1 = Nand(I1,I2)
  G2 = Nand(G1,I2)
  G3 = Nand(G1,I2)
  G4 = Nand(G2,I1)
  G5 = Nand(I2,G3)
  @test nv == ([I1,I2],[G1,G2,G3,G4,G5],[G4,G5])
  IN = nv[1];
  XX = nv[2]   # 5 gates, while Macia uses only the first 3 gates for X
  X = XX[1:3]  # 3 gates, what Macia calls X
  O = nv[3]
  @test nv == (IN,XX,O)  # Failed on 8/15/21 because node caching in evaluate_node() is commented out
  println("Verifying the results in comparison to Tables 3 and 4 of Macia and Sole")
  gbX = get_bits( X, p0.numinputs )
  gbO = get_bits( O, p0.numinputs )
  X2_1 = X[1:2]
  Xh2_1 = X[3:3]
  gbX2_1 = get_bits( X2_1, p0.numinputs )
  gbXh2_1 = get_bits( Xh2_1, p0.numinputs )
  println("(gbX,gbO): ",(gbX,gbO))
  miX_O = mutual_information(gbX,gbO,base=exp(1))
  jeX_O = joint_entropy(gbX,gbO,base=exp(1))
  print("je(X,O): ",jeX_O)
  println("   mi(X,O): ",miX_O)
  println("(gbX2_1,gbO): ",(gbX2_1,gbO))
  jeX2_1_O = joint_entropy(gbX2_1,gbO,base=exp(1))
  miX2_1_O = mutual_information(gbX2_1,gbO,base=exp(1))
  println("je(X2_1,O): ",jeX2_1_O,"  mi(X2_1,O): ",miX2_1_O)
  jeXh2_1_O = joint_entropy(gbXh2_1,gbO,base=exp(1))
  miXh2_1_O = mutual_information(gbXh2_1,gbO,base=exp(1))
  println("je(Xh2_1,O): ",jeXh2_1_O,"  mi(Xh2_1,O): ",miXh2_1_O)
  println("je(X2_1,O) + je(Xh2_1,O) - je(X,O): ",jeX2_1_O + jeXh2_1_O - jeX_O)
  println("mi(X2_1,O) + mi(Xh2_1,O) - mi(X,O): ",miX2_1_O + miXh2_1_O - miX_O)
  #mi0 = map(x->(x,mutual_information(x,gbO)),map(x->get_bits(x,p0.numinputs),[s for s in combinations(X,2)]))
  #println("mi0: ",mi0)

  # Simplified and modified version of the above circuit with only one output gate
  println("Alternate simple chromosome: ")
  p1 = Parameters(1, 4, 0.05, 0.0, 2, 1, 2, 4, 8)
  in1 = InputNode(1)
  in2 = InputNode(2)
  int3 = InteriorNode(NAND, [1, 2])  # (in1,in2)
  int4 = InteriorNode(NAND, [3, 2])  # (int3,in2)
  int5 = InteriorNode(NAND, [1, 4])  # (in1,int4)
  int6 = InteriorNode(NAND, [5, 3])  # (int5,int3)
  out7 = OutputNode(6)
  c1 = Chromosome(p1,[in1,in2],[int3,int4,int5,int6],[out7],0.0,0.0)
  context = construct_contexts(p1.numinputs)[p1.numinputs]
  execute_chromosome(c1,context)
  #println("c1: ",c1)
  #print_chromosome(c1)
  println("num_mutate_locations: ",num_mutate_locations(c1,funcs))
  @assert number_active(c1) == number_active_old(c1)
  nv = node_values(c1)
  println("nv: ",nv)

  I1 = context[1]
  I2 = context[2]
  INT3 = Nand(I1,I2)
  INT4 = Nand(INT3,I2)
  INT5 = Nand(I1,INT4)
  INT6 = Nand(INT5,INT3)
  OUT7 = INT6
  @test nv == ([I1,I2],[INT3,INT4,INT5,INT6],[OUT7])
  IN = nv[1];
  X = nv[2]
  O = nv[3]
  @test nv == (IN,X,O)
  gbX = get_bits( X, p1.numinputs )
  gbO = get_bits( O, p1.numinputs )
  println("(gbX,gbO): ",(gbX,gbO))
  println("mi(X,O): ",mutual_information(gbX,gbO))
  mi1 = map(x->(x,mutual_information(x,gbO)),map(x->get_bits(x,p1.numinputs),[s for s in combinations(X,2)]))
  println("mi1: ",mi1)

  # Duplicate the above circuit with each part going to a separate output
  p2 = Parameters(1, 4, 0.05, 0.0, 2, 2, 2, 8, 16)
  print_parameters(p2)
  in1 = InputNode(1)
  in2 = InputNode(2)
  int3 = InteriorNode(NAND, [1, 2])  # (in1,in2)
  int4 = InteriorNode(NAND, [3, 2])  # (int3,in2)
  int5 = InteriorNode(NAND, [1, 4])  # (in1,int4)
  int6 = InteriorNode(NAND, [5, 3])  # (int5,int3)
  int7 = InteriorNode(NAND, [1, 2])  # (in1,in2)
  int8 = InteriorNode(NAND, [7, 2])  # (int7,in2)
  int9 = InteriorNode(NAND, [1, 8])  # (in1,int8)
  int10 = InteriorNode(NAND, [9, 7])  # (int9,int7)
  out11 = OutputNode(6)
  out12 = OutputNode(10)
  c2 = Chromosome(p2,[in1,in2],[int3,int4,int5,int6,int7,int8,int9,int10],[out11,out12],0.0,0.0)
  context = construct_contexts(p2.numinputs)[p2.numinputs]
  execute_chromosome(c2,context)
  #println("c2: ",c2)
  print_chromosome(c2)
  println("num_mutate_locations: ",num_mutate_locations(c2,funcs))
  @assert number_active(c2) == number_active_old(c2)
  nv = node_values(c2)
  println("nv: ",nv)

  I1 = context[1]
  I2 = context[2]
  INT3 = Nand(I1,I2)
  INT4 = Nand(INT3,I2)
  INT5 = Nand(I1,INT4)
  INT6 = Nand(INT5,INT3)
  INT7 = Nand(I1,I2)
  INT8 = Nand(INT7,I2)
  INT9 = Nand(I1,INT8)
  INT10 = Nand(INT9,INT7)
  OUT11 = INT6
  OUT12 = INT10
  @test nv == ([I1,I2],[INT3,INT4,INT5,INT6,INT7,INT8,INT9,INT10],[OUT11,OUT12])
  IN = nv[1];
  X = nv[2]
  O = nv[3]
  @test nv == (IN,X,O)
  gbX = get_bits( X, p2.numinputs )
  gbO = get_bits( O, p2.numinputs )
  gbX2_1 = get_bits(X2_1, p2.numinputs )
  println("(gbX,gbO): ",(gbX,gbO))
  println("mi(X,O): ",mutual_information(gbX,gbO))
  mi2 = map(x->(x,mutual_information(x,gbO)),map(x->get_bits(x,p2.numinputs),[s for s in combinations(X,2)]))
  println("mi2: ",mi2)
end