text:
export ProgramBuilder, Program, ProgramMemory, prealloc_readonly!, prealloc_unique!, prealloc_stack!, get_global_ptr

struct ArrayDetails{N}
    dims::NTuple{N, Int}
    type::Type

    ArrayDetails(dims::NTuple{N, Int}, type::Type) where N = new{N}(dims, type)
end

mutable struct ProgramBuilder
    device::Device
    readonly_memory::Vector{AbstractArray}
    global_memory::Vector{Pair{Symbol, ArrayDetails}}
    unique_memory::Vector{ArrayDetails}
    current_stack_size::Int

    readonly_size::Int
    unique_memory_size::Int
    global_memory_size::Int
    memory_size::Int
    max_stack_size::Int

    parent_ref

    ProgramBuilder(device::Symbol) = ProgramBuilder(get_device(device))
    ProgramBuilder(device::Device) = new(Function[], device, AbstractArray[], Pair{Symbol, ArrayDetails}[], ArrayDetails[], 0, 0, 0, 0, 0, 0, nothing)

    (p::ProgramBuilder)() = Program(p)
end

function getarray(memory_block, type, offset, dims)
    len = prod(dims) * sizeof(type)
    return (reshape(reinterpret(type, view(memory_block, offset:(offset + len - 1))), dims...), offset + len)
end

const ReadonlyMemType = 1
const UniqueMemType = 2
const GlobalMemType = 3
const StackMemType = 4

struct PrealloctionPointer
    type::Int
    index
    pb::ProgramBuilder

    PrealloctionPointer(type::Int, index, pb::ProgramBuilder) = new(type, index, pb)

    function (p::PrealloctionPointer)()
        (p.pb.parent_ref === nothing) && error("Cannot get array as program has not finished building!")
        (p.type == ReadonlyMemType) && return p.pb.parent_ref.readonly_memory[p.index]
        (p.type == UniqueMemType) && return p.pb.parent_ref.unique_memory[p.index]
        (p.type == GlobalMemType) && return p.pb.parent_ref.global_memory[p.index]
        if p.type == StackMemType
            init = p.index[1]
            range = init:(init + prod(p.index[3]) * sizeof(p.index[2]) - 1)
            return reshape(reinterpret(p.index[2], view(p.pb.parent_ref.stack, range)), p.index[3]...)
        end
        error("Unknown Memory Pointer Type")
    end
end

Base.push!(pb::ProgramBuilder, fun::Function) = push!(pb.programs, fun)

"Get a pointer accessing global memory"
get_global_ptr(pb::ProgramBuilder, sym::Symbol) =  PrealloctionPointer(GlobalMemType, sym, pb)

"Preallocate readonly memory. Will optimize globally in the program to reduce memory"
function prealloc_readonly!(pb::ProgramBuilder, array)
    idx = findfirst([item == array for item in pb.readonly_memory])
    if idx === nothing
        push!(pb.readonly_memory, array)
        pb.readonly_size += length(array) * sizeof(eltype(array))
        idx = length(pb.readonly_memory)
    end
    return PrealloctionPointer(ReadonlyMemType, idx, pb)
end

"Preallocate unique memory that nothing else is allowed to use. Can Read/Write"
function prealloc_unique!(pb::ProgramBuilder, type::Type, dims...)
    push!(pb.unique_memory, ArrayDetails(dims, type))
    pb.unique_memory_size += prod(dims) * sizeof(type)
    return PrealloctionPointer(UniqueMemType, length(pb.unique_memory), pb)
end

"Preallocate global memory that anything can read/write to"
function prealloc_global!(pb::ProgramBuilder, sym::Symbol, type::Type, dims...)
    push!(pb.global_memory, Pair(sym, ArrayDetails(dims, type)))
    pb.global_memory_size += prod(dims) * sizeof(type)
    return PrealloctionPointer(GlobalMemType, sym, pb)
end

"Preallocate stack memory. Will be overwritten constantly. Only use for temporary information in a single threaded environment"
function prealloc_stack!(pb::ProgramBuilder, type::Type, dims...)
    out = PrealloctionPointer(StackMemType, (pb.current_stack_size + 1, type, dims), pb)
    pb.current_stack_size += prod(dims) * sizeof(type)
    return out
end

function Base.flush(pb::ProgramBuilder)
    pb.max_stack_size = max(pb.current_stack_size, pb.max_stack_size)
    pb.current_stack_size = 0
    pb.memory_size = pb.readonly_size + pb.global_memory_size + pb.unique_memory_size
end

struct ProgramMemory{ArrayType}
    readonly_memory::Array{ArrayType}
    unique_memory::Array{ArrayType}
    global_memory::Dict{Symbol, ArrayType}
    stack::ArrayType
    memory_block::ArrayType

    stack_size::Int
    readonly_size::Int
    unique_memory_size::Int
    global_memory_size::Int
    memory_size::Int

    function ProgramMemory(pb::ProgramBuilder)
        atype = arraytype(pb.device)
        offset = 1

        memory_block = alloc(pb.device, UInt8, pb.memory_size + pb.max_stack_size)
        stack = view(memory_block, (pb.memory_size + 1):length(memory_block))
        readonly_memory = Array{atype}(undef, length(pb.readonly_memory))
        unique_memory = Array{atype}(undef, length(pb.unique_memory))
        global_memory = Dict{Symbol, atype}()

        for idx in eachindex(pb.readonly_memory)
            arr = pb.readonly_memory[idx]
            arrmem = reinterpret(UInt8, arr)
            copyto!(memory_block, offset, arrmem)
            res = getarray(memory_block, eltype(arr), offset, size(arr))
            readonly_memory[idx] = res[1]
            offset = res[2]
        end

        for idx in eachindex(pb.unique_memory)
            res = getarray(memory_block, pb.unique_memory[idx].type, offset, pb.unique_memory[idx].dims)
            unique_memory[idx] = res[1]
            offset = res[2]
        end

        for glob in pb.global_memory
            res = getarray(memory_block, glob[2].type, offset, glob[2].dims)
            global_memory[glob[1]] = res[1]
            offset = res[2]
        end

        new{atype}(readonly_memory, unique_memory, global_memory, stack, memory_block, pb.max_stack_size, pb.readonly_size, pb.unique_memory_size, pb.global_memory_size, pb.memory_size + pb.max_stack_size)
    end
end

mutable struct Program{ArrayType}
    device::Device
    mem::ProgramMemory{ArrayType}

    function Program(pb::ProgramBuilder)
        flush(pb)
        mem = ProgramMemory(pb)
        pb.parent_ref = mem
        new{arraytype(pb.device)}(pb.device, mem)
    end

    function (p::Program)(arg)
        return p.program(arg)
    end
end


text:
# Simple example of a Neumann problem:
#
#           y
#             |
#             |     u_y = cos(omega x)
#          Ly ----------------------------------
#             |                                |
#             |                                |
#             |                                |
#     u_x = 0 |    -u_xx - u_yy + u = 0        | u_x = 0
#             |                                |
#             |                                |
#             |                                |
#           -----------------------------------|-----  x
#                         u_y = 0              Lx
#
# Here, omega = 2 n pi / Lx, and the exact solution is
#
#      u = cos(omega x) * cosh(alpha y) / (alpha sinh(alpha Ly)
#
# where alpha^2 = 1 + omega^2.

using FinElt
using FinElt.PlanarPoisson
using Printf

include("params.jl")

const omega = 2 * pi / Lx
const alpha = sqrt(1.0 + omega^2)

function exact_u(x)
    return cos(omega*x[1]) * cosh(alpha*x[2]) / (alpha*sinh(alpha*Ly))
end

function gN(x)
    return cos(omega*x[1])
end

maxerr = zeros(refinements+1)
@printf("%10s  %12s  %8s  %8s\n\n",
        "N", "max error", "rate", "seconds")
for k = 0:refinements
    start = time()
    mesh = read_msh_file("rect$k.msh")
    vp = VariationalProblem(mesh)
    add_bilin_form!(vp, "Omega", grad_dot_grad!)
    add_bilin_form!(vp, "Omega", func_times_func!)
    add_lin_functnl!(vp, "Top", bdry_source_times_func!, gN)
    A, b = assembled_linear_system(vp)
    ufree = A \ b
    uh = complete_soln(ufree, vp)
    u = get_nodal_vals(exact_u, mesh)
    finish = time()
    maxerr[k+1] = maximum(abs.(uh-u))
    N = length(ufree)
    if k == 0
        @printf("%10d  %12.4e\n", N, maxerr[k+1])
    else
        rate = log2(maxerr[k]/maxerr[k+1])
        elapsed = finish - start
        @printf("%10d  %12.4e  %8.4f  %8.4f\n",
                N, maxerr[k+1], rate, elapsed)
    end
end

ext: jl


text:
#
# Date created: 2022-03-07
# Author: aradclif
#
#
############################################################################################
# Experiments in efficiency improvements
using Base.Cartesian

A = rand(4,3,5);
Da = axes(A)
Db = Da
dims = (1,3)
Db = Base.setindex(Db, 1, dims[1])
Db = Base.setindex(Db, 1, dims[2])
B = similar(A, Db);
fill!(B, 0.0);
sum(B)
@macroexpand @nloops 3 i A begin
    (@nref 3 B d -> length(Db[d]) == 1 ? 1 : i_d) =
        +((@nref 3 B d -> length(Db[d]) == 1 ? 1 : i_d), (@nref 3 A i))
end
@macroexpand @nloops 3 i A begin
    (@nref 3 B d -> length(axes(B, d)) == 1 ? 1 : i_d) =
        +((@nref 3 B d -> length(axes(B, d)) == 1 ? 1 : i_d), (@nref 3 A i))
end
@macroexpand @nref 3 B d -> length(axes(B, d)) == 1 ? 1 : i_d
@macroexpand @nexprs 3 d -> (j_d = length(axes(B, d)) == 1 ? 1 : i_d)
@macroexpand @nref 3 A i
Expr(:ref, :A, :(1), :i_2, :(3))
(LoopVectorization.StaticInt(1), LoopVectorization.StaticInt(3), LoopVectorization.StaticInt(1))

@generated function prepare(A, dims::NTuple{N, Int}) where {N}
    quote
        Da = axes(A)
        Db = Da
        @nexprs $N d -> (Db = Base.setindex(Db, Base.OneTo(1), dims[d]))
        Db
    end
end

@generated function vreduce!(f::Function, B::AbstractArray{T, N}, A::AbstractArray{T, N}) where {T, N}
    quote
        @nloops $N i A begin
            (@nref $N B d -> length(axes(B, d)) == 1 ? 1 : i_d) =
                f((@nref $N B d -> length(axes(B, d)) == 1 ? 1 : i_d), (@nref $N A i))
        end
        B
    end
end

# @generated function _vreduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}) where {F, T, N}
#     quote
#         ee = @macroexpand @nloops $N i A begin
#             (@nref $N B d -> length(axes(B, d)) == 1 ? 1 : i_d) =
#                 $(f.instance)((@nref $N B d -> length(axes(B, d)) == 1 ? 1 : i_d), (@nref $N A i))
#         end
#         ea = ee.args[2]
#         return quote
#             $ea
#             return B
#         end
#     end
#     # ee = quote
#     #     @macroexpand @nloops 3 i A begin
#     #     (@nref 3 B d -> length(axes(B, d)) == 1 ? 1 : i_d) =
#     #         +((@nref 3 B d -> length(axes(B, d)) == 1 ? 1 : i_d), (@nref 3 A i))
#     #     end
#     # end
#     # ee.args[2]
# end

function preparedims(B::AbstractArray{T, N}, A::AbstractArray{T, N}) where {T, N}
    Db = axes(B)
    Dᵥ = ntuple(d -> length(Db[d]) == 1 ? StaticInt(1) : Symbol(:i_, d), N)
end
# preparedims(B, A)

# function vvreduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}) where {F, T, N}
#     dims = preparedims(B, A)
#     _vvreduce!(f, B, A, dims)
# end
# @generated function _vvreduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, T, N, D}
#     e = Expr(:ref, :B, ntuple(d -> D.parameters[d] !== StaticInt{1} ? Symbol(:i_, d) : :(1), N)...)
#     quote
#         ex = @macroexpand @nloops $N i A begin
#             $e = $(F.instance)($e, (@nref $N A i))
#         end
#         loops = ex.args[2]
#         return quote
#             # @turbo $loops
#             $loops
#             return B
#         end
#         # @turbo $loops
#         # return B
#     end
# end
# _vvreduce!(+, B, A, preparedims(B, A))
function branches_quote(F, N::Int, D)
    e = Expr(:ref, :B, ntuple(d -> D.parameters[d] !== StaticInt{1} ? Symbol(:i_, d) : :(1), N)...)
    # or:
    # e = Expr(:ref, :B)
    # params = D.parameters
    # for d ∈ 1:N
    #     push!(e.args, params[d] !== StaticInt{1} ? Symbol(:i_, d) : :(1))
    # end
    ex = :(@nloops $N i A begin
               $e = $(F.instance)($e, (@nref $N A i))
           end)
    exm = macroexpand(Main, ex)
    loops = exm.args[2]
    return quote
        @turbo $loops
        return B
    end
end
@generated function vvvreduce!(f::F, B::AbstractArray{T, N}, A::AbstractArray{T, N}, dims::D) where {F, T, N, D}
    branches_quote(F, N, D)
end
function vvvreduce(f, A::AbstractArray{T}; dims::NTuple{N, Int}) where {N} where {T}
    if N == ndims(A)
        cat(_vvvreduce(f, A); dims=N)
    else
        Db = prepare(A, dims)
        # B = similar(A, Db)
        # fill!(B, 0.0)
        B = zeros(T, Db)
        vvvreduce!(f, B, A, preparedims(B, A))
    end
end
# vvvreduce(f, A::AbstractArray{T}; dims::Int) where {T} = vvvreduce(f, A; dims=(dims,))
@generated function _vvvreduce(f::F, A::AbstractArray{T, N}) where {F, T, N}
    quote
        s = zero($T)
        @turbo for i ∈ eachindex(A)
            s = $(F.instance)(s, A[i])
        end
        return s
    end
end
@timev vvvreduce(+, A)
A = rand(4,3,5);

# branches_quote(typeof(+), 3, typeof(preparedims(B, A)))
# dd = preparedims(B, A)
# @timev vvvreduce!(+, B, A, dd)

@benchmark vvvreduce(+, A, dims=dims)
@benchmark sum(A, dims=dims)
@timev vvvreduce(+, A, dims=(1,))
# This causes an error -- however, one can actually perform this faster by just
# reducing using a scalar, then doing cat(s, dims=N).
B = zeros(Float64, prepare(A, (1,)))
tp = typeof(+)
DD = typeof(preparedims(B, A))
@timev branches_quote(tp, 3, DD)
@timev branches_quote2(tp, 3, DD)
eval(ans)
@timev sum(A, dims=(5,));
@timev vvvreduce(+, A, dims=(5,));

@generated function indexexpr(D::NTuple{N, T}) where {N, T}
    ex = quote
        @ntuple $N d -> length(D[d]) == 1 ? Symbol(1) : Symbol("i_$(d)")
    end
    # return :($ex)
    ex
end
Isym = indexexpr(Db2)
@nref 3 B j -> ($(Isym[j]))

Db2 = prepare(A, dims)
B2 = similar(A, Db2);
fill!(B2, 0.0);
vreduce!(+, B2, A);
_vreduce!(+, B2, A)
vreduce6!(+, B2, A)
vreduce2!(+, B2, A);

s = 0.0
ex = quote @nloops 3 i A begin
    s += @nref 3 A i
end
end
ex2 = :(@nloops 3 i A begin
    s += @nref 3 A i
end)

@macroexpand @turbo for i ∈ eachindex(A)
    s = +(s, A[i])
end

LoopVectorization.@_turbo(ex)

ex = quote
    for i ∈ eachindex(A)
        s += A[i]
    end
end
ex2 = :(for i ∈ eachindex(A)
            s += A[i]
        end)

N = 3
static_dims = [1,2,3]
quote
    @turbo $loops
end
quote
    @turbo $loop3
end

loop3 = @macroexpand @nloops 3 i A begin
    (@nref 3 B d -> length(axes(B, d)) == 1 ? 1 : i_d) =
        +((@nref 3 B d -> length(axes(B, d)) == 1 ? 1 : i_d), (@nref 3 A i))
end
ee = loop3.args[2]
quote
    @turbo $ee
end

function simdsum(A::Array{T, 3}) where {T}
    s = zero(T)
    @turbo for k ∈ axes(A, 3)
        nothing
        for j ∈ axes(A, 2)
            for i ∈ axes(A, 1)
                s = (+)(s, A[i, j, k])
            end
        end
    end
    s
end

simdsum(A)


text:
using PropBiology

include("testprop.jl")

max_stars_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/PropBiology


text:
"""
            Objects for SPH mapping.

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

# using Distributed

"""
    mappingParameters(; x_lim::Vector{Float64},
                        y_lim::Vector{Float64},
                        z_lim::Vector{Float64},
                        pixelSideLength::Float64=0.0,
                        Npixels::Int64=0)

Parameter object for sph to grid mapping. Define either `pixelSideLength` of `Npixels`.
"""
struct mappingParameters

    x_lim::Vector{Float64}
    y_lim::Vector{Float64}
    z_lim::Vector{Float64}
    center::Vector{Float64}
    pixelSideLength::Float64
    pixelArea::Float64
    Npixels::Int64
    x_size::Float64
    y_size::Float64
    z_size::Float64
    x::Vector{Float64}
    y::Vector{Float64}
    z::Vector{Float64}

    function mappingParameters(;x_lim::Vector{Float64}   = [-1.0, -1.0],
                                y_lim::Vector{Float64}   = [-1.0, -1.0],
                                z_lim::Vector{Float64}   = [-1.0, -1.0],
                                center::Vector{Float64}  = [-1.0, -1.0, -1.0],
                                x_size::Float64          =  -1.0,
                                y_size::Float64          =  -1.0,
                                z_size::Float64          =  -1.0,
                                pixelSideLength::Float64 =  -1.0,
                                Npixels::Int64           =   0)

        # calculate limits if center position and sizes are given
        if ( x_lim == [-1.0, -1.0] && ( y_lim == [-1.0, -1.0] && z_lim == [-1.0, -1.0] ))

            if (center != [-1.0, -1.0, -1.0]) && ( x_size != -1.0 && (y_size != -1.0 && z_size != -1.0) )

               x_lim = [ center[1] - 0.5x_size, center[1] + 0.5x_size ]
               y_lim = [ center[2] - 0.5y_size, center[2] + 0.5y_size ]
               z_lim = [ center[3] - 0.5z_size, center[3] + 0.5z_size ]

            else
                error("Giving a center position requires extent in x, y and z direction.")
            end
        end

        # calculate side lengths from limits
        if x_size == -1.0
            x_size = x_lim[2] - x_lim[1]
        end
        if y_size == -1.0
            y_size = y_lim[2] - y_lim[1]
        end
        if z_size == -1.0
            z_size = z_lim[2] - z_lim[1]
        end

        # find the maximum extent of the map
        max_size = max(x_size, y_size)

        if (pixelSideLength == -1.0) & (Npixels != 0)
            pixelSideLength = max_size/Npixels
        elseif (pixelSideLength != -1.0) & (Npixels == 0)
            Npixels = floor(Int64, max_size/pixelSideLength)
        else
            error("Please specify pixelSideLenght or number of pixels!")
        end

        pixelArea = pixelSideLength^2

        x = collect(range(x_lim[1]+0.5pixelSideLength, x_lim[2]-0.5pixelSideLength, step = pixelSideLength))
        y = collect(range(y_lim[1]+0.5pixelSideLength, y_lim[2]-0.5pixelSideLength, step = pixelSideLength))
        z = collect(range(z_lim[1]+0.5pixelSideLength, z_lim[2]-0.5pixelSideLength, step = pixelSideLength))

        new(x_lim, y_lim, z_lim,
            center,
            pixelSideLength,
            pixelArea,
            Npixels,
            x_size, y_size, z_size,
            x, y, z)

    end
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
struct UnitCommitment <: AbstractOperationsModel end

function UnitCommitment(system::PSY.System, transmission::Type{S}; optimizer::Union{Nothing,JuMP.OptimizerFactory}=nothing, kwargs...) where {S <: PM.AbstractPowerFormulation}

    devices = Dict{Symbol, PSI.DeviceModel}(:ThermalGenerators => PSI.DeviceModel(PSY.ThermalGen, PSI.ThermalUnitCommitment),
                                            :RenewableGenerators => PSI.DeviceModel(PSY.RenewableGen, PSI.RenewableFullDispatch),
                                            :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))

branches = Dict{Symbol, PSI.DeviceModel}(:Lines => PSI.DeviceModel(PSY.Branch, PSI.SeriesLine))
services = Dict{Symbol, PSI.ServiceModel}(:Reserves => PSI.ServiceModel(PSY.Reserve, PSI.AbstractReservesForm))
    return PowerOperationModel(UnitCommitment,
                                   transmission,
                                    devices,
                                    branches,
                                    services,
                                    system;
                                    optimizer = optimizer, kwargs...)
end


text:
function _split_classifcation_error_loss(y::Vector, X::DataFrame, obs_row_indcs::Vector{Int}, column_indcs::Vector{Int}, weights::Vector)
    best = NO_BEST
    best_val = -Inf
    X = X[obs_row_indcs, :]
    y = y[obs_row_indcs]

    for j in column_indcs
        keep_row = !isna(X[:, j])                # use only non-NA values
        x_obs = convert(Vector, X[keep_row, j])
        y_obs = y[keep_row]
        wgts_obs = weights[keep_row]
        if length(unique(x_obs)) ≥ 100
            domain_j = quantile(x_obs, linspace(0.01, 0.99, 99))
        else
            domain_j = sort(unique(x_obs))
        end

        for thresh in domain_j[2:end]
            cur_split = x_obs .< thresh         # tons of mem-alloc here

            # ensure we have some observations in y
            if (length(y_obs[cur_split]) > 0) && (length(y_obs[!cur_split]) > 0)
                value = _classifcation_error_loss(y_obs[cur_split], wgts_obs[cur_split]) + _classifcation_error_loss(y_obs[!cur_split], wgts_obs[!cur_split])
                if value > best_val
                    best_val = value
                    best = (j, thresh)
                end
            end
        end
    end
    return best
end

# This method is dispatched when weights are omitted. This
# allows us to compute the loss function 5x faster
function _split_classifcation_error_loss(y::Vector, X::DataFrame, obs_row_indcs::Vector{Int}, column_indcs::Vector{Int})
    best = NO_BEST
    best_val = -Inf
    X = X[obs_row_indcs, :]
    y = y[obs_row_indcs]

    for j in column_indcs
        keep_row::BitArray{1} = !isna(X[:, j])
        x_obs = convert(Vector{Float64}, X[keep_row, j])
        y_obs = y[keep_row]

        domain_j = unique(x_obs)
        if length(domain_j) ≥ 100
            x_obs = convert(Array{Float64, 1}, x_obs)          # can't be Array{Any,1} for quantile()
            domain_j = quantile(x_obs, linspace(0.01, 0.99, 99))
        end
        if length(domain_j) > 1
            for thresh in domain_j[2:end]

                cur_split = x_obs .< thresh         # tons of mem-alloc here

                # ensure we have some observations in y
                if (any(cur_split)) && (!all(cur_split))
                    value = _classifcation_error_loss(y_obs[cur_split]) + _classifcation_error_loss(y_obs[!cur_split])
                    if value > best_val
                        best_val = value
                        best = (j, thresh)
                    end
                end
            end
        end
    end
    return best
end

# n = 1000
# p = 20
# X = DataFrame(randn(n, p));
# y = rand([true, false], n);
# wgt = ones(Int, n);

# @code_warntype _split_classifcation_error_loss(y, X, collect(1:n), collect(1:p), wgt)

#
#
#
# @time surrogate_splits(y, X, collect(1:n), collect(1:p), 5, wgt)

#


text:
function hit_box(minx::T, miny::T, minz::T, maxx::T, maxy::T, maxz::T, ray::Ray{T}) where {T<:AbstractFloat}
    minx2v1 = minx - ray.x
    maxx2v1 = maxx - ray.x
    miny2v1 = miny - ray.y
    maxy2v1 = maxy - ray.y
    minz2v1 = minz - ray.z
    maxz2v1 = maxz - ray.z

    tmin = minx2v1 / ray.dx
    tmax = maxx2v1 / ray.dx
    if tmin > tmax
        tmin, tmax = tmax, tmin
    end

    tymin = miny2v1 / ray.dy
    tymax = maxy2v1 / ray.dy
    if tymin > tymax
        tymin, tymax = tymax, tymin
    end

    if (tmin > tymax) || (tymin > tmax)
        return false
    end

    if tymin > tmin
        tmin = tymin
    end

    if tymax < tmax
        tmax = tymax
    end

    tzmin = minz2v1 / ray.dz
    tzmax = maxz2v1 / ray.dz
    if tzmin > tzmax
        tzmin, tzmax = tzmax, tzmin
    end

    if (tmin > tzmax) || (tzmin > tmax)
        return false
    end

    return true
end

function intersect_walls(walls::Walls, ray::Ray)
    r1::Float64 = Inf
    r2::Float64 = Inf
    r3::Float64 = Inf

    xwall::Char = 'o'
    ywall::Char = 'o'
    zwall::Char = 'o'
    wall::Char = 'x'

    if ray.dx > 0.0
        r1 = (walls.east - ray.x) / ray.dx
        xwall = 'e'
    elseif ray.dx < 0.0
        r1 = (walls.west - ray.x) / ray.dx
        xwall = 'w'
    end

    if ray.dy > 0.0
        r2 = (walls.north - ray.y) / ray.dy
        ywall = 'n'
    elseif ray.dy < 0.0
        r2 = (walls.south - ray.y) / ray.dy
        ywall = 's'
    end

    if ray.dz > 0.0
        r3 = (walls.top - ray.z) / ray.dz
        zwall = 't'
    elseif ray.dz < 0.0
        r3 = (walls.bottom - ray.z) / ray.dz
        zwall = 'b'
    end

    r = min(r1, r2, r3)

    if r === r1
        wall = xwall
    elseif r === r2
        wall = ywall
    else
        wall = zwall
    end

    x = ray.x + ray.dx * r
    y = ray.y + ray.dy * r
    z = ray.z + ray.dz * r

    I::Point = Point(x, y, z)

    return I, wall, r
end

function cycle_ray(walls::Walls, ray::Ray)
    # How to type multiple return values?
    I, wall, r = intersect_walls(walls, ray)

    x::Float64 = I.x
    y::Float64 = I.y
    z::Float64 = I.z

    xsame::Bool = false
    ysame::Bool = false
    zsame::Bool = false

    # make these constants?
    small_num::Float64 = 0.000001
    u::Float64 = 100.0

    if I.x <= (walls.west + small_num)
        x = I.x + walls.xsize
    elseif I.x >= (walls.east - small_num)
        x = I.x - walls.xsize
    else
        xsame = true
    end

    if I.y <= (walls.south + small_num)
        y = I.y + walls.ysize
    elseif I.y >= (walls.north - small_num)
        y = I.y - walls.ysize
    else
        ysame = true
    end

    if I.z <= (walls.bottom + small_num)
        z = I.z + walls.zsize
    elseif I.z >= (walls.top - small_num)
        z = I.z - walls.zsize
    else
        zsame = true
    end

    # assert that at least one booleans evaluates to true
    # skipped for now...
    new_ray::Ray = Ray(x, y, z, ray.dx, ray.dy, ray.dz)

    return new_ray
end

function intersect_bvh(bvh::Bvh{T,S,A,E}, ray::Ray{T}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    #NOTE:these two vector take about a third of the time of this function and furthermore the pushes take a long time too...
    # the hit_box line takes surprisingly little, with only 586/4505
    q::Vector = [1]
    leaf_idxs::Vector{S} = []
    while size(q, 1) > 0
        i = pop!(q)
        minx, miny, minz, maxx, maxy, maxz = bvh.xmin[i], bvh.ymin[i], bvh.zmin[i], bvh.xmax[i], bvh.ymax[i], bvh.zmax[i]
        if hit_box(minx, miny, minz, maxx, maxy, maxz, ray)::Bool
            if bvh.left_child[i] > 1 && bvh.right_child[i] > 1
                q = push!(q, bvh.left_child[i], bvh.right_child[i])
            end
            if bvh.bvh2bb[i] > 0
                leaf_idxs = push!(leaf_idxs, bvh.bvh2bb[i])
            end
        end
    end
    return leaf_idxs
end

function solve_quadratic(a::T, b::T, c::T) where {T<:AbstractFloat}
    x0::T = 0.0
    x1::T = 0.0
    discr::T = b*b - 4*a*c
    if discr < 0
        return false, x0, x1
    elseif (discr==0)
        x0 = x1 = -0.5*b/a
    else
        q = b>0 ? -0.5*(b+sqrt(discr)) : -0.5*(b-sqrt(discr))
        x0 = q / a
        x1 = c / q
    end
    if x0 > x1
        x0, x1 = x1, x0
    end
    return true, x0, x1
end

function intersect(balls::Ball{T,S,A,E}, i::S, ray::Ray{T}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray, E<:AbstractArray}
    Lx::T = ray.x - balls.p1x[i]
    Ly::T = ray.y - balls.p1y[i]
    Lz::T = ray.z - balls.p1z[i]
    a::T = ray.dx*ray.dx + ray.dy*ray.dy + ray.dz*ray.dz
    b::T = 2 * (Lx*ray.dx + Ly*ray.dy + Lz*ray.dz)
    c::T = Lx*Lx + Ly*Ly + Lz*Lz - balls.rsq[i]
    hit, t0, t1 = solve_quadratic(a, b, c)
    if !hit
        nullpoint = Point{T}(Inf, Inf, Inf)
        return false, nullpoint, nullpoint, Inf
    end
    if t0 < 0 # if t0 is negative then use t1
	t0 = t1
    end
    if t0 >= 0
	Ix = ray.x + ray.dx * t0
	Iy = ray.y + ray.dy * t0
	Iz = ray.z + ray.dz * t0
	I = Point{T}(Ix, Iy, Iz)
	dx = Ix - balls.p1x[i]
        dy = Iy - balls.p1y[i]
        dz = Iz - balls.p1z[i]
        l = sqrt(dx^2 + dy^2 + dz^2)
        dx /= l
        dy /= l
        dz /= l
        normal::Point = Point{T}(dx, dy, dz)
	return hit, I, normal, t0
    end
    nullpoint = Point{T}(Inf, Inf, Inf)
    return false, nullpoint, nullpoint, Inf
end

function intersect(disks::Disk{T,S,A,E}, i::S, ray::Ray{T}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray, E<:AbstractArray}
    denom::T = disks.nx[i]*ray.dx + disks.ny[i]*ray.dy + disks.nz[i]*ray.dz
    if abs(denom) > SMALL_NUM
        # check point of intersection on the plane...
        p0l0x = disks.p1x[i] - ray.x
        p0l0y = disks.p1y[i] - ray.y
        p0l0z = disks.p1z[i] - ray.z
        t = (p0l0x*disks.nx[i] + p0l0y*disks.ny[i] + p0l0z*disks.nz[i]) / denom
        if t < 0
            nullpoint = Point{T}(Inf, Inf, Inf)
            return false, nullpoint, nullpoint, Inf
        end
        Ix = ray.x + ray.dx * t
        Iy = ray.y + ray.dy * t
        Iz = ray.z + ray.dz * t
        # see if it is on the disk...
        vx = Ix - disks.p1x[i]
        vy = Iy - disks.p1y[i]
        vz = Iz - disks.p1z[i]
        vv = vx*vx + vy*vy + vz*vz
        hit = vv <= disks.rsq[i]
        I = Point{T}(Ix, Iy, Iz)
	normal::Point{T} = Point{T}(disks.nx[i], disks.ny[i], disks.nz[i])
        return hit, I, normal, t
    end
    nullpoint = Point{T}(Inf, Inf, Inf)
    return false, nullpoint, nullpoint, Inf
end

function intersect(cone::Cone{T,S,A,E}, i::S, ray::Ray) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray, E<:AbstractArray}
  # intersect ray with cone/cylinder...
  xoff = cone.m1x[i]
  yoff = cone.m1y[i]
  zoff = cone.m1z[i]
  beta = cone.beta[i]
  gamma = cone.gamma[i]

  tray = transform(ray, -xoff, -yoff, -zoff, -beta, -gamma)

  a::Float64 = tray.dx^2 + tray.dy^2 - tray.dz^2 * cone.aa[i]
  b::Float64 = 2 * tray.x * tray.dx + 2 * tray.y * tray.dy - 2 * cone.r1a[i] * tray.dz - 2 * tray.z * ray.dz * cone.aa[i]
  c::Float64 = tray.x^2 + tray.y^2 - cone.r1sq[i] - 2 * cone.r1a[i] * tray.z - cone.aa[i] * tray.z^2
  d::Float64 = b^2 - (4*a*c)

  t1::Float64 = Inf
  t2::Float64 = Inf
  t::Float64 = Inf

  if d > 0
    t1 = (-1.0*b - sqrt(d))/(2*a)
    t2 = (-1.0*b + sqrt(d))/(2*a)
  end

  # choose the smallest, non-negative solution for t
  if t1 <= 0
    if t2 > 0
      t = t2
    end
  else
    if t2 <= 0
      t = t1
    elseif t1 < t2
      t = t1
    else
      t = t2
    end
  end
  # t1 and t2 can both be negative, hence t can be negative!

  hit::Bool = false
  tI = Point{T}(tray.x + tray.dx * t, tray.y + tray.dy * t, tray.z + tray.dz * t)
  if (t > 0) && (tI.z > 0.0) && (tI.z < cone.length[i])
    hit = true
  end

  I = Point{T}(ray.x + ray.dx * t, ray.y + ray.dy * t, ray.z + ray.dz * t)
  xx::Float64 = tI.x^2
  yy::Float64 = tI.y^2
  rcirc::Float64 = sqrt(xx + yy)
  nx::Float64, ny::Float64, nz::Float64 = normalize(tI.x, tI.y, cone.a[i] * -1.0 * rcirc)
  nx, ny, nz = transform(nx, ny, nz, xoff, yoff, zoff, beta, gamma)
  normal::Point = Point{T}(nx, ny, nz)

  return hit, I, normal, t
end

function intersect(cylinder::Cylinder{T,S,A,E}, i::S, ray::Ray) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray, E<:AbstractArray}
  # intersect ray with cylinder...

  xoff = cylinder.m1x[i]
  yoff = cylinder.m1y[i]
  zoff = cylinder.m1z[i]
  beta = cylinder.beta[i]
  gamma = cylinder.gamma[i]

  tray = transform(ray, -xoff, -yoff, -zoff, -beta, -gamma)

  a::Float64 = tray.dx^2 + tray.dy^2
  b::Float64 = 2 * tray.x * tray.dx + 2 * tray.y * tray.dy
  c::Float64 = tray.x^2 + tray.y^2 - cylinder.rsq[i]
  d::Float64 = b^2 - (4*a*c)

  t1::Float64 = Inf
  t2::Float64 = Inf
  t::Float64 = Inf

  if d > 0
    t1 = (-1.0*b - sqrt(d))/(2*a)
    t2 = (-1.0*b + sqrt(d))/(2*a)
  end

  # choose the smallest, non-negative solution for t
  if t1 <= 0
    if t2 > 0
      t = t2
    end
  else
    if t2 <= 0
      t = t1
    elseif t1 < t2
      t = t1
    else
      t = t2
    end
  end
  # t1 and t2 can both be negative, hence t can be negative!

  hit::Bool = false
  tI = Point{T}(tray.x + tray.dx * t, tray.y + tray.dy * t, tray.z + tray.dz * t)
  if (t > 0) && (tI.z > 0.0) && (tI.z < cylinder.length[i])
    hit = true
  end

  I = Point{T}(ray.x + ray.dx * t, ray.y + ray.dy * t, ray.z + ray.dz * t)
  nx::Float64, ny::Float64, nz::Float64 = normalize(tI.x, tI.y, 0.0)
  nx, ny, nz = transform(nx, ny, nz, xoff, yoff, zoff, beta, gamma)
  normal::Point = Point{T}(nx, ny, nz)

  return hit, I, normal, t
end

function intersect(mesh::Mesh{T,S,A,E}, i::S, ray::Ray{T}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    w0x = ray.x - mesh.v1x[i]
    w0y = ray.y - mesh.v1y[i]
    w0z = ray.z - mesh.v1z[i]
    a = -((mesh.nx[i] * w0x) + (mesh.ny[i] * w0y) + (mesh.nz[i] * w0z))
    b = (ray.dx * mesh.nx[i]) + (ray.dy * mesh.ny[i]) + (ray.dz * mesh.nz[i])

    if abs(b) < SMALL_NUM
        nullPoint = Point{T}(Inf, Inf, Inf)
        return false, nullPoint, nullPoint, Inf # parallel
    end

    r = a / b
    if r < 0.0
        nullPoint = Point{T}(Inf, Inf, Inf)
        return false, nullPoint, nullPoint, r # away
    end

    # point of plane intersection
    x = ray.x + r * ray.dx
    y = ray.y + r * ray.dy
    z = ray.z + r * ray.dz
    I = Point{T}(x, y, z)
    N = Point{T}(mesh.nx[i], mesh.ny[i], mesh.nz[i])

    mx = x - mesh.v1x[i]
    my = y - mesh.v1y[i]
    mz = z - mesh.v1z[i]
    mu = (mx * mesh.ux[i]) + (my * mesh.uy[i]) + (mz * mesh.uz[i])
    mw = (mx * mesh.wx[i]) + (my * mesh.wy[i]) + (mz * mesh.wz[i])

    s = ((mesh.uw[i] * mw) - (mesh.ww[i] * mu)) / mesh.D[i]
    if (s < 0.0) | (s > 1.0)
        return false, I, N, r # outside
    end

    t = ((mesh.uw[i] * mu) - (mesh.uu[i] * mw)) / mesh.D[i]
    if (t < 0.0) | ((s + t) > 1.0)
        return false, I, N, r # outside
    end

    return true, I, N, r # ray intersects triangle
end

function intersect(geometry::Geometry{T,S,A,E}, ray::Ray{T}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    hit::Bool = false
    In::Point{T} = Point{T}(Inf, Inf, Inf)
    eidx::S = 0
    normal::Point{T} = Point{T}(Inf, Inf, Inf)
    rmin::T = Inf

    # intersect the individual elements of this geometry (e.g. facets)...
    eidxs::Array{S,1} = collect(1:size(geometry.mtl, 1))
    for i::S in eidxs
        hit_, In_, N_, r_ = intersect(geometry, i, ray)
        if (hit_ === true) && (r_ < rmin)
            eidx = i
            rmin = r_
            In = In_
            normal = N_
            hit = true
        end
    end
    return hit, In, eidx, normal, rmin
end

function intersect(geometry::Geometry{T,S,A,E}, bvh::Bvh{T,S,A,E}, ray::Ray{T}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    hit::Bool = false
    In::Point{T} = Point{T}(Inf, Inf, Inf)
    eidx::S = 0
    normal::Point{T} = Point{T}(Inf, Inf, Inf)
    rmin::T = Inf

    # intersect mesh facets
    #eidxs = collect(1:size(geometry.mtl, 1))
    eidxs::Array{S,1} = intersect_bvh(bvh, ray)
    for i::S in eidxs
        # and see if it is nearest
        hit_, I_, N_, r_ = intersect(geometry, i, ray)
        if (hit_ === true) && (r_ < rmin)
            eidx = i
            rmin = r_
            In = I_
            normal = N_
            hit = true
        end
    end
    return hit, In, eidx, normal, rmin
end

function intersect_scene(assets::Array{Asset{T,S},1}, geometries::AbstractArray{Geometry{T,S,A,E},1}, geometry_bvhs::AbstractArray{Bvh{T,S,A,E},1}, scene_bvh::Bvh{T,S,A,E}, ray::Ray{T}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    hit::Bool = false
    In::Point{T} = Point{T}(Inf, Inf, Inf)
    nrml::Point{T} = Point{T}(Inf, Inf, Inf)
    tray::Ray{T} = Ray{T}(Inf, Inf, Inf, Inf, Inf, Inf)
    rmin::T = Inf
    instidx::S = 0
    aidx::S = 0 # asset index
    oidx::S = 0 # object index / geometry index
    eidx::S = 0 # element index (of the geometry indexed with oidx)

    #aidxs::Array{S} = collect(1:length(assets))
    aidxs::Array{S} = intersect_bvh(scene_bvh, ray)
    for i::S in aidxs
        oidx_ = assets[i].oidx

        geometry = geometries[oidx_]
        geometry_bvh = geometry_bvhs[oidx_]

        tray_ = transform_inverse(ray, assets[i]) # from world space to object space
        hit_, In_, eidx_, nrml_, rmin_ = intersect(geometry, geometry_bvh, tray_)
        if (hit_ === true) && (rmin_ < rmin)
            hit = true
            x = ray.x + ray.dx * rmin_
            y = ray.y + ray.dy * rmin_
            z = ray.z + ray.dz * rmin_
            In = Point{T}(x, y, z)
            tray = tray_
            nrml = nrml_
            rmin = rmin_
            aidx = i
            oidx = oidx_
            eidx = eidx_
        end
    end
    return hit, In, tray, nrml, rmin, aidx, oidx, eidx
end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
using BesselK, SpecialFunctions, ArbNumerics
include("shared.jl")

function wbesselkxv(v,x)
  try
    return BesselK.adbesselkxv(v,x)
  catch
    return NaN
  end
end

function rbesselkxv(v,x)
  _v = ArbReal(v)
  _x = ArbReal(x)
  xv = _x^_v
  Float64(xv*ArbNumerics.besselk(ArbFloat(v), ArbFloat(x)))
end
abesselkxv(v,x) = (x^v)*SpecialFunctions.besselk(v, x)

const VGRID = range(0.25, 10.0, length=100)
const XGRID = range(0.0,   8.0, length=201)[2:end] # since other impls can't do x=0.

const BASELINE = [rbesselkxv(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const AMOS     = [abesselkxv(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const OURSOL   = [wbesselkxv(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const TOLS_A   = atolfun.(zip(BASELINE, AMOS))
const TOLS_U   = atolfun.(zip(BASELINE, OURSOL))
const TOLS_AU  = rtolfun.(zip(AMOS,     OURSOL))

#=
gnuplot_save_matrix!("../plotdata/atols_amos.csv",      TOLS_A,  VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_ours.csv",      TOLS_U,  VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/rtols_amos_ours.csv", TOLS_AU, VGRID, XGRID)
=#


text:
"""
Given an array of values and associated array of booleans, sort the values into
true / false and move those values to separate arrays
"""

using CUDA

CUDA.allowscalar(false)
"""
For a batch of size `n` what is the lowest index of the batch `i` is in
"""
function batch_floor(idx, n)
    return idx - (idx - 1) % n
end

"""
For a batch of size `n` what is the highest index of the batch `i` is in
"""
function batch_ceil(idx, n)
    return idx + n - 1 - (idx - 1) % n
end

"""
GPU friendly step function
"""
function Θ(i)
    return 1 & (1 <= i)
end

"""
Suppose we are merging two lists of size n, each of which has all falses before
all trues. Together, they will be indexed 1:2n. This is a fast stepwise function
for the destination index of a value at index `x` in the concatenated input,
where `a` is the number of falses in the first half, b = n - a, and false is the
number of falses in the second half.
"""
function step_swap(x, a, b, c)
    return x + Θ(x - a) * b - Θ(x - (a + c)) * (b + c) + Θ(x - (a + b + c)) * c
end

"""
Generalizes `step_swap` for when the floor index is not 1
"""
function batch_step_swap(x, n, a, b, c)
    idx = (x - 1) % n + 1
    return batch_floor(x, n) - 1 + step_swap(idx, a, b, c)
end

"""
For thread `idx` with current value `value`, merge two batches of size `n` and
return the new value this thread takes. `sums` and `swap` and shared mem
"""
function merge_swap_shmem(value, idx, n, sums, swap)
    b = sums[batch_floor(idx, 2 * n)]
    a = n - b
    d = sums[batch_ceil(idx, 2 * n)]
    c = n - d
    swap[idx] = value
    sync_threads()
    sums[idx] = d + b
    return swap[batch_step_swap(idx, 2 * n, a, b, c)]
end

"""
Given an array `flags`, perform a stable sort into `dest` of the indices of
`flags`.
Eg: [True, False, False, True] => [2, 3, 1, 4]
Note, stores indices modulo 1024 (max thread block size). Assumed use is with
`index_move_kernel`
"""
function index_sort_kernel(dest, flags, batch_sums)
    i = threadIdx().x
    idx0 = (blockIdx().x - 1) * blockDim().x + i # for addressing
    idx = i # for use in shared mem swaps
    swap = @cuStaticSharedMem(Int32, 1024)
    sums = @cuStaticSharedMem(Int32, 1024)
    sums[i] = 1 & flags[idx0]
    idx = merge_swap_shmem(idx, i, 1, sums, swap)
    idx = merge_swap_shmem(idx, i, 2, sums, swap)
    idx = merge_swap_shmem(idx, i, 4, sums, swap)
    idx = merge_swap_shmem(idx, i, 8, sums, swap)
    idx = merge_swap_shmem(idx, i, 16, sums, swap)
    idx = merge_swap_shmem(idx, i, 32, sums, swap)
    idx = merge_swap_shmem(idx, i, 64, sums, swap)
    idx = merge_swap_shmem(idx, i, 128, sums, swap)
    idx = merge_swap_shmem(idx, i, 256, sums, swap)
    idx = merge_swap_shmem(idx, i, 512, sums, swap)

    if i == 1
        batch_sums[blockIdx().x] = sums[i]
    end
    dest[idx0] = idx
    return nothing
end

"""
Given an array of `values` reorder them. Denote v0 as initial `values`, v1 as
modified. The result of this kernel is v1[i] = v0[indices[i]].
Makes use of shared memory so that global mem accesses are coalesced.
WARNING: Assumes T is small enough that we won't run out of shmem
"""
function index_move_kernel(values :: CuDeviceArray{T}, indices) where T
    i = threadIdx().x
    idx = (blockIdx().x - 1) * blockDim().x + i
    swap = @cuStaticSharedMem(T, 1024)
    swap[i] = values[idx]
    sync_threads()
    new_value = swap[indices[idx]]
    sync_threads()
    values[idx] = new_value
    return nothing
end

function global_move_kernel(values)

end


text:
function CatStr(s1::Array,sep::AbstractString,s2::Array)
    #Assume s1 and s2 are arrays of String
    #Also sep is a string
    s12 = s1 .* [sep] .* s2
    return s12
end


text:
using ProgressMeter

open("82.jl.out") do f
    ftc = Dict()
    ft = Dict()
    fc = Dict()
    N = 0
    @showprogress 1 "ahaha" for line in readlines(f)
        words = filter(x -> !isempty(x) && x != "N/A", strip.(split(line)))
        word, context = words[1], words[2:end]

        # f(t,c)
        for c in context
            t = (word, c)
            if t in keys(ftc)
                ftc[t] += 1
            else
                ftc[t] = 1
            end
        end

        # f(t,*)
        if word in keys(ft)
            ft[word] += 1
        else
            ft[word] = 1
        end

        # f(*.c)
        for c in context
            if c in keys(fc)
                fc[c] += 1
            else
                fc[c] = 1
            end
        end

        N += 1
    end

    @show N
end


text:
# tests for properties

using ImageGeoms
using FillArrays: Trues, Falses, Zeros, Ones
using Test: @test, @testset, @test_throws, @inferred

# test 2D properties
function image_geom_test2(ig::ImageGeom)
    @test ig.dims isa Tuple
    @test ig.deltas isa Tuple
    @test ig.offsets isa Tuple

    @test ig.dim isa Tuple
    @test ig.x isa AbstractVector
    @test ig.y isa AbstractVector
    @test ig.wx isa Number
    @test ig.wy isa Number
    @test ig.xg isa AbstractArray
    @test ig.yg isa AbstractArray
    @test ig.fovs isa Tuple
    @test ig.np isa Int
#   ig.mask_outline
    @test ig.ones isa Ones
    @test ig.zeros isa Zeros
    @test ig.u isa AbstractVector
    @test ig.v isa AbstractVector
    @test ig.ug isa AbstractArray
    @test ig.vg isa AbstractArray
#   ig.fg

    @test ig.shape(vec(ig.ones)) == ig.ones
    @test ig.embed(ig.ones[ig.mask]) == ig.mask
    @test ig.maskit(ig.ones) == ones(ig.np)
#   @inferred # todo
    @inferred ig.unitv()
    ig.unitv(j=4)
    ig.unitv(i=ntuple(i->1, length(ig.dim)))
    ig.unitv(c=ntuple(i->0, length(ig.dim)))
#= todo-i: why do these fail?
    @inferred image_geom_ellipse(8, 10, 1, 2)
    @inferred ig.circ()
    @inferred ig.plot()
    @inferred ig.unitv()
=#
    ig.circ()
    how = (args...; kw...) -> args[1] # trick to avoid Plots
    ig.plot(how)
    ig.unitv()
#= todo-i:
    @inferred ig.down(2)
    @inferred ig.over(2)
=#
    @test ig.down(2) isa ImageGeom
    @test ig.over(2) isa ImageGeom
#   image_geom_ellipse(8, 10, 1, 2)
    true
end

function image_geom_test2()
    #@inferred
#   image_geom(nx=16, dx=2, offsets=:dsp, mask=:all_but_edge_xy)
#   @test_throws String image_geom(nx=16, dx=1, offsets=:bad)
#   @test_throws String image_geom(nx=16, dx=1, mask=:bad) # mask type
#   @test_throws String image_geom(nx=16, dx=1, mask=trues(2,2)) # mask size
    ig = ImageGeom(MaskAll(); dims=(9,8), deltas=(2,2))
    image_geom_test2(ig) # todo: fails at embed
    ig = ImageGeom(MaskCircle(); dims=(9,8), deltas=(2,2))
    image_geom_test2(ig)
    ig.over(2)
    ig.down(3) # test both even and non-even factors
    ig.help
    true
end

function image_geom_test3(ig::ImageGeom{3})
    @test image_geom_test2(ig)
    ig.wz
    ig.zg
    ig.mask_or
    @inferred ig.expand_nz(2)
#   @inferred cbct(ig)
#   @test_throws String image_geom(nx=1, dx=1, nz=2, offset_z=-1, offsets=:dsp)
    true
end

function image_geom_test3()
    ig = @inferred ImageGeom( dims = (9,8,4) )
    image_geom_test3(ig)
    true
end

@testset "2d" begin
    @test all(ImageGeom((2,), (3,), (0,)).mask)
    @test image_geom_test2()
end

@testset "3d" begin
    @test image_geom_test3()
end


text:
function triggeredby!(o::AbstractObservable, a::AbstractObservable, b::AbstractObservable)
    f = on(t -> setindex!(a, t), o)
    on(t -> setindex!(o, a[], notify = x -> x != f), b)
    o
end

triggeredby(a::AbstractObservable{T}, b::AbstractObservable) where {T} =
    triggeredby!(Observable{T}(a[]), a, b)

"""
`onchange(w::AbstractWidget, change = w[:changes])`

Return a widget that's identical to `w` but only updates on `change`. For a slider it corresponds to releasing it
and for a textbox it corresponds to losing focus.

## Examples

```julia
sld = slider(1:100) |> onchange # update on release
txt = textbox("Write here") |> onchange # update on losing focuse
```
"""
function onchange(w::AbstractWidget{T}, change = w[:changes]) where T
    o = triggeredby(w, change)
    Widget{T}(w, output = o)
end


text:
module OddFrames
#=========
o=o= Heirarchy overview =o=o
|AbstractOddFrame
|_____ OddFrame
|_____ ImmutableOddFrame
|_____ ImageOddFrame
=========#
abstract type AbstractOddFrame end
include("type/frame.jl")
export OddFrame, getindex, setindex, shape
export read_csv
include("methods.jl")
export shape, show
end


text:
using Iterators

function normalize_neighbors!(m, g::Graph, X::AbstractArray, i::Integer)
    Z = sum(m[i,neighbors(g,i), :])
    #@show i, Z
    for j in neighbors(g,i)
        for k in X
            m[i,j,k] /= Z
        end
    end
    return Z
end

function assert_nonzero(val)
    @assert val != 0 "val was zero"
    return val
end

function propogate!{T<:Real, U<:Integer}(m::AbstractArray{T,3},
                                         ϕ::AbstractArray{T,2},
                                         ψ::AbstractArray{T,2}, A, X::AbstractArray{U, 1})
    n = size(A, 1)
    for i in 1:n
        for j in neighbors(A,i)
            if j == i
                continue
            end
            for k in X
                m[i,j,k] = 0
                for l in X
                    w = ϕ[i,l] * ψ[l,k]
                    @assert w != NaN
                    sp = sum([p!=j ? log(assert_nonzero(m[p,i,l])) : 0 for p in neighbors(A, i)])
                    p = exp(sp) + 1e-16
                    if p == 0
                        warn("sp is $sp")
                        warn("p is 0! $i,$j,$k,$l")
                    end
                    @assert p != NaN
                    @assert m[j,i,l] != 0 "$i,$j,$l have m 0"
                    m[i,j,k] += w * p #/ m[j,i,l]
                end
            end
        end
        Z = normalize_neighbors!(m, A, X, i)
    end
    return m
end

function beliefs!{T<:Real}(m::AbstractArray{T}, b, ϕ, A, X)
    for i in 1:size(A,1)
        for k in X
            w = prod(T[m[p,i,k] for p in neighbors(A, i)])
            if w == 0
                w = 1e-16
            end
            b[i,k] = ϕ[i,k] * w
            @assert b[i,k] >= 0 "Unnormalized beliefs are not all positive offending indices ($i, $k)=$(b[i,k])"
        end
        b[i,:] /= sum(b[i,:])
    end
    #@show b
    @assert all( 0.99 .< sum(b,2) .< 1.01) "Beliefs do not form a probability distributions"
    return b
end

type BeliefProblem{T<:Real, U<:Integer, G}
    m::AbstractArray{T,3}
    ϕ::AbstractArray{T,2}
    ψ::AbstractArray{T,2}
    A::G
    X::AbstractArray{U, 1}
end

type Messages{T<:Real, G} <: AbstractArray{T, 3}
    nrows::Int
    ncols::Int
    edges::Dict{Tuple{Int, Int}, Int}
    states::UnitRange
    storage::Array{T, 2}
end

function Messages(g::Graph, states::UnitRange)
    n=nv(g)
    d = Dict{Tuple{Int,Int}, Int}()
    i = 1
    for e in edges(g)
        d[(src(e), dst(e))] = i
        d[(dst(e), src(e))] = i+1
        i += 2
    end
    storage = ones(Float64, (i, length(states)))
    return Messages{Float64, Graph}(n,n,d, states, storage)
end

Base.size(m::Messages) = (m.nrows, m.ncols, length(m.states))
function Base.getindex(m::Messages, I::Vararg{Int})
    pair = I[1:2]
    if pair in keys(m.edges)
        return m.storage[m.edges[pair], I[3]]
    else
        throw(KeyError(pair))
    end
end

function Base.setindex!(m::Messages, v, I::Vararg{Int})
    pair = I[1:2]
    if pair in keys(m.edges)
        return m.storage[m.edges[pair], I[3]] = v
    else
        throw(KeyError(pair))
    end
end

logodds(b::AbstractMatrix) = log.(b[:,1]./b[:,2])

function beliefprop(g::Graph, ϕ, ψ, maxiter)
    n = size(g, 1)
    k = 2
    m = ones(Float64, (n,n,k))
    m = Messages(g, 1:k)
    #@show sort(collect(keys(m.edges)))
    pr = BeliefProblem(m, ϕ, ψ, g, 1:k)
    iterstates = []
    diffs = []
    i = 1
    m = nth(iterate(m->begin msg = propogate!(m, ϕ,ψ, g, 1:k);
            state = [maximum(msg.storage), minimum(msg.storage)]
            diff = i > 1 ? state-iterstates[i-1]:0;
            i+=1
            #@show diff
            push!(iterstates, state)
            push!(diffs, diff)
            #msg.storage = max(msg.storage, 1e-60)
           return msg
        end, m), maxiter)
    b = similar(pr.ϕ)
#     for i in 1:n
#         if sum(pr.ϕ[i,:]) != sum(b[i, :])
#             println("Crud! "*string(pr.ϕ[i,:])*" "*string(b[i,:]))
#         end
#     end
    beliefs!(pr.m, b, pr.ϕ, g, 1:k)
    lodds = logodds(b)
    #plotbeliefs(g, b[:,1])
    #display(bar(lodds, ylabel="Log Odds", xlabel="Vertex ID"))
    return pr, b, lodds
end

function plotbeliefs(g::Graph, b::AbstractVector)
    #nodefillc=map(x-> RGB(1x, 0, 1-x), b./maximum(b))
    #plo = gplot(H, lx, ly, nodelabel=[si2s[v] for v in 1:nv(g)], nodelabelsize=sizes, nodefillc=nodefillc)
    #draw(PNG("plots/postbeliefprop.png", 60cm, 60cm), plo)
    #return plo
    #return gplot(g, layout=(g)-> spring_layout(g; C=1), nodelabel=1:nv(g), nodelabelsize=1, nodefillc=nodefillc)
    #return gplot(g, layout=(g)-> begin pos = (200*begin Z=spring_layout(g); Z[1]end, log2(float(collect(1:nv(g))))); @show typeof(pos); return pos end, nodelabel=1:nv(g), nodelabelsize=1, nodefillc=nodefillc)
    #return gplot(g, layout=(g)-> (exp10(spectral_layout(g)[1]), log(float(vertices(g)))),nodelabel=1:nv(g), nodelabelsize=1, nodefillc=nodefillc)
end

function Psis(ϵ::Real)
    return [1-ϵ ϵ;
            ϵ 1-ϵ]
end

max_stars_repo_name: jpfairbanks/NewsGraphing

max_issues_repo_name: jpfairbanks/NewsGraphing

max_forks_repo_name: jpfairbanks/NewsGraphing


text:
#!/usr/bin/env julia

using SnoopCompile

SnoopCompile.@snoop "gadfly_compiles.csv" begin
    include("runtests.jl")
end

using Gadfly

blacklist = ["__init__", "Image",
             "unpack_pango_attr_list",
             "newsurface",
             "PangoAttr",
             "update_pango_attr",
             "pango_fmt_float",
             "unpack_pango_attr",
             "pango_text_extents",
             "cairo_linecap",
             "cairo_linejoin"]

data = SnoopCompile.read("gadfly_compiles.csv")
pc, discards = SnoopCompile.parcel(data[end:-1:1,2], blacklist=blacklist)
SnoopCompile.write("snoop", pc)


text:
module InfluxDB

using HTTPClient.HTTPC
using JSON

export InfluxDBClient

# package code goes here
type InfluxDBClient
  host::String
  port::Int64
  username::String
  password::String
  database::String
  #ssl::Bool
  #verifyssl
  #udp
  #scheme
  baseurl::String
  headers::Dict{String, String}

  function InfluxDBClient(host::String, port::Int64, username::String, password::String, database::String)
	baseurl = "$host:$port"
	headers = {"Content-Type" => "application/json"}
	new(host, port, username, password, database, baseurl, headers)
  end

end # InfluxDBClient

end # module


text:
__precompile__()

module Coils

using Statistics
using LinearAlgebra
using LightGraphs
using ProgressMeter

export cuboid_system, getedgei, getedge, find_cells, biotsavart,
    biotsavart_cell, system_matrix, find_edgecurrents, simplify_current_graph,
    findsimpleloop, find_all_simpleloops, solve_system, cuboid_poi,
    decompose, decompose_currents, field_loops, μ0, order_cell,
    real_decomposed_currents, save_result

# To be able to write linspace.(starts, stops, lengths)
linspace(start, stop, length) = range(start, stop = stop, length = length)

include("CoilsPlot.jl")

"skipfaces is in the order [x-, x+, y-, y+, z-, z+]"
function cuboid_system(totalsize, ntiles; skipfaces = falses(6))
    gridplanes = linspace.(-totalsize / 2, totalsize / 2, ntiles .+ 1)
    vertex_positions = []

    # construct the vertices
    for z in gridplanes[3], y in gridplanes[2], x in gridplanes[1]
        any(in.([x,y,z], extrema.(gridplanes))) || continue

        # skip the inside of the X- face
        skipfaces[1] &&
            x == minimum(gridplanes[1]) &&
            !(y in extrema(gridplanes[2])) &&
            !(z in extrema(gridplanes[3])) &&
            continue

        # skip the inside of the X+ face
        skipfaces[2] &&
            x == maximum(gridplanes[1]) &&
            !(y in extrema(gridplanes[2])) &&
            !(z in extrema(gridplanes[3])) &&
            continue

        # skip the inside of the Y- face
        skipfaces[3] &&
            y == minimum(gridplanes[2]) &&
            !(x in extrema(gridplanes[1])) &&
            !(z in extrema(gridplanes[3])) &&
            continue

        # skip the inside of the Y+ face
        skipfaces[4] &&
            y == maximum(gridplanes[2]) &&
            !(x in extrema(gridplanes[1])) &&
            !(z in extrema(gridplanes[3])) &&
            continue

        # skip the inside of the Z- face
        skipfaces[5] &&
            z == minimum(gridplanes[3]) &&
            !(x in extrema(gridplanes[1])) &&
            !(y in extrema(gridplanes[2])) &&
            continue

        # skip the inside of the Z+ face
        skipfaces[6] &&
            z == maximum(gridplanes[3]) &&
            !(x in extrema(gridplanes[1])) &&
            !(y in extrema(gridplanes[2])) &&
            continue

        push!(vertex_positions, [x, y, z])
    end

    # add the edges
    g = DiGraph()
    add_vertices!(g, length(vertex_positions))

    tile_size = mean.(diff.(gridplanes))
    for i in eachindex(vertex_positions)
        for j in (i+1):length(vertex_positions)
            # the "mean(tile_size) is not the best, at the moment. Only works for equal tiles
            if abs(norm(vertex_positions[i] .- vertex_positions[j]) - mean(tile_size)) < 0.01
                add_edge!(g, i, j)
            end
        end
    end
    g, vertex_positions
end

"Identify the edge between vertices a and b (in any direction)."
function getedgei(g, a, b)
    for (i, e) in enumerate(edges(g))
        if (dst(e) == a) & (src(e) == b)
            return i
        elseif (dst(e) == b) & (src(e) == a)
            return i
        end
    end
    return 0
end

"Identify the edge between vertices a and b (in any direction)."
function getedge(g, a, b)
    for e in edges(g)
        if (dst(e) == a) & (src(e) == b)
            return e
        elseif (dst(e) == b) & (src(e) == a)
            return e
        end
    end
    return 0
end

"""Find elementary cells in a graph.
Takes a very long time to find large cells (longer then 13-15 edges).
In those cases large cells can be supplied manually, which drastically speeds up the process.
"""
function find_cells(g; maxpathlength = Inf, cells = [])
    # get all cells in the raster
    edgetaken = zeros(Int, length(edges(g)))

    # mark the manually supplied cells
    for cell in cells
        for k in 1:(length(cell) - 1)
            # mark that the edge is part of one additional cell
            edgetaken[getedgei(g, cell[k], cell[k + 1])] += 1
        end
    end

    # an algorithm to decompose a graph into "cells"
    # we traverse the graph breadth-first, so that the first cell we find
    # has to be the shortest one (smallest cell, so elementary)

    @showprogress 1 for i in eachindex(vertices(g))
        #println("Starting vertex: $i")
        paths = [[i]]
        foundnow = 0

        while !isempty(paths)
            # get all cells including that vertex
            # take the last added path
            path = popfirst!(paths)
            length(path) > maxpathlength && break

            for n in [ collect(outneighbors(g, path[end])); collect(inneighbors(g, path[end])) ]
                # try to extend the path
                if edgetaken[getedgei(g, path[end], n)] >= 2
                    # edge already a member of two cells
                    continue
                elseif (length(path) > 2) && (n == path[1]) && !(sort([path; n]) in sort.(cells))
                    # we are back to the starting point
                    # and the path is not a cell yet
                    for k in 1:(length(path) - 1)
                        # mark that the edge is part of one additional cell
                        edgetaken[getedgei(g, path[k], path[k + 1])] += 1
                    end
                    edgetaken[getedgei(g, path[end], n)] += 1

                    # save the cell
                    push!(cells, [path; n])
                    foundnow += 1
                    break
                elseif !(n in path)
                    # we have just extended the path, maybe it will close to be
                    # a cell further on
                    push!(paths, [path; n])
                end
            end

            # a very generous stop condition
            # if the number of cells found for this vertex is at least number of its neighbors
            # typically can stop earlier
            total_neighbors = (length(outneighbors(g, i)) + length(inneighbors(g, i)))
            already_in = count(cell -> i in cell, cells)
            already_in >= total_neighbors && break
        end
    end
    cells
end

"Sine of the angle between vectors"
sinvectors(x1, x2) = norm(x1 × x2) / (norm(x1) * norm(x2))

const μ0 = 4π * 1e-7

function biotsavart(x1, x2, p)
    # n is a normal vector along the wire
    n = normalize(x2 .- x1)

    # d is a vector perpendicular to the wire through the point p
    # ref. https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation
    # ref. https://de.wikipedia.org/wiki/Biot-Savart-Gesetz#Gerader_Linienleiter
    ρ = (x1 .- p) .- ((x1 .- p) ⋅ n) * n

    # In the formula on wikipedia the angles are directed.
    # We need to take it into account.
    # s1 = norm(p .+ ρ .- x2) > norm(x1 .- x2) ? 1 : -1
    # s2 = norm(p .+ ρ .- x1) > norm(x1 .- x2) ? 1 : -1
    # normB = μ0 * abs(s1 * sinvectors(ρ, x2 - p) - s2 * sinvectors(ρ, x1 - p)) / (4π * norm(ρ))

    # check with s whether poit p+ρ lies on the wire
    s = norm(p .+ ρ .- (x1 .+ x2) / 2) > norm(x1 .- x2) / 2 ? -1 : 1
    normB = μ0 * abs(sinvectors(ρ, x2 - p) + s * sinvectors(ρ, x1 - p)) / (4π * norm(ρ))

    B = normalize((ρ × n)) * normB
end

function numbiotsavart(x1, x2, p; n = 1000)
    B = zeros(3)
    dl = (x2 .- x1) / n

    for x in [ x1 .+ m * dl for m in 0.5:n ]
        r = p .- x
        B += μ0 / 4π * (dl × r) / norm(r)^3
    end
    B
end

"Magnetic field in point p produced by ith cell"
function biotsavart_cell(p, vertex_positions, cell)
    vrts = [ vertex_positions[k] for k in cell ]
    sum([ biotsavart(vrts[i], vrts[i + 1], p) for i in 1:length(vrts)-1 ])
end

function field_loops(p, loops, currents, vertex_positions)
    field = zeros(3)
    for (loop, current) in zip(loops, currents)
        for iv in 1 : length(loop)-1
            x1 = vertex_positions[loop[iv]]
            x2 = vertex_positions[loop[iv+1]]
            field .+= biotsavart(x1, x2, p) * current
        end
    end
    field
end

function system_matrix(poi, vertex_positions, cells)
    M = hcat([
        # for the ith cell the proportionality constant for the field it produces
        # in every point of interest. x, y and z are simply concatenated in an xyzxyzxyz way
        vcat([ biotsavart_cell(p, vertex_positions, cell) for p in poi ]...)
        # this is done for each cell
        for cell in cells ]...)
end

function find_edgecurrents(g, cells, optI; bigfloat = false)
    edgecurrents = Dict([
        (edge, bigfloat ? BigFloat(0.0) : 0.0)
        for edge in edges(g) ])

    for edge in edges(g)
        for icell in eachindex(cells)
            # for each cell
            vrts = cells[icell]
            for i in 1:length(vrts)-1
                # for each edge of the cell
                if vrts[i] == src(edge) && vrts[i + 1] == dst(edge)
                    # the cell is aligned with the edge
                    edgecurrents[edge] += optI[icell]
                end
                if vrts[i] == dst(edge) && vrts[i + 1] == src(edge)
                    # the cell is anti-aligned with the edge
                    edgecurrents[edge] -= optI[icell]
                end
            end
        end
    end

    edgecurrents
end

function simplify_current_graph(g, edgecurrents; tolerance = 1e-13)
    gs = DiGraph()
    add_vertices!(gs, nv(g))

    edgecurrentss = empty!(Dict(edgecurrents))

    for edge in edges(g)
        if abs(edgecurrents[edge]) < tolerance
            # skip edges below the machine precision
            continue
        elseif edgecurrents[edge] > 0
            # put this edge on the new graph
            add_edge!(gs, src(edge), dst(edge))
            edgecurrentss[getedge(gs, src(edge), dst(edge))] = abs(edgecurrents[edge])
        elseif edgecurrents[edge] < 0
            # put this edge on the new graph, but invert it first
            add_edge!(gs, dst(edge), src(edge))
            edgecurrentss[getedge(gs, dst(edge), src(edge))] = abs(edgecurrents[edge])
        end
    end

    gs, edgecurrentss
end

function findsimpleloop(g, edgecurrents; maxcurrent = Inf, maxpathlength = Inf,
        verbose = false, tolerance = 1e-13)

    for current in sort(sort(collect(values(edgecurrents))), rev = true)
        current <= nextfloat(maxcurrent) || continue
        verbose && println("    Trying current $current, number of edges: $(length(edges(g)))")
        # start from the highest current and gradually reduce it
        for i in vertices(g)
            # try to find a loop with at least this high current
            paths = [[i]]

            safetycounter = 0
            while !isempty(paths)
                # get all cells including that vertex
                # take the last added path
                path = popfirst!(paths)
                length(path) > maxpathlength && break
                verbose && (safetycounter += 1) % 10_000 == 0 && println("  Pathlength: $(length(path))")

                # only go in the direction of positive current, so outcoming edges
                for n in outneighbors(g, path[end])
                    # add a safety margin for threshold
                    if edgecurrents[getedge(g, path[end], n)] <= (current - tolerance)
                        # cannot go here, too high current
                        continue
                    end
                    # try to extend the path
                    if (length(path) > 2) && (n == path[1])
                        # we are back to the starting point
                        # and the path is not a cell yet
                        # Return the current and the simple loop
                        verbose && println("  Found a simple loop with current $current : $([path; n])")
                        return (current, [path; n])
                    elseif !(n in path)
                        # we have just extended the path, maybe it will close to be
                        # a cell further on
                        push!(paths, [path; n])
                    end
                end
            end
        end
    end

    # could not find a loop
    NaN, []
end

function find_all_simpleloops(g, edgecurrents; tolerance = 1e-10,
        mincurrent = 0, verbose = false)

    simpleloops = []
    simpleloopscurrents = []

    gs, edgecurrentss = simplify_current_graph(g, edgecurrents,
        tolerance = tolerance)

    # Only look for loops with currents not smaller then the current
    # of the most recently found loop.
    # current = Inf
    while length(edges(gs)) > 0
        current, simpleloop = findsimpleloop(gs, edgecurrentss,
            # maxcurrent = current,
            verbose = verbose,
            tolerance = tolerance)
        if isnan(current)
            println("Could not find a loop, stopping.")
            println("   The graph has $(length(edges(g))) edges left,")
            println("   the maximum current left is $(maximum(abs, values(edgecurrentss)))")
            break
        end

        abs(current) < mincurrent && break

        push!(simpleloopscurrents, current)
        push!(simpleloops, simpleloop)

        # subtract the simple loop from the graph
        for i in 1:length(simpleloop)-1
            edge = getedge(gs, simpleloop[i], simpleloop[i+1])

            if src(edge) == simpleloop[i]
                # the edge is aligned with the simple loop
                edgecurrentss[ edge ] -= current
            elseif dst(edge) == simpleloop[i]
                # the edge is anti-aligned with the simple loop
                edgecurrentss[ edge ] += current
            end
        end

        gs, edgecurrentss = simplify_current_graph(gs, edgecurrentss,
            tolerance = tolerance)
    end

    simpleloops, simpleloopscurrents
end

function printlnflush(s)
    println(s)
    flush(stdout)
end

function solve_system(g, vertex_positions, poi, B; initialcells = [],
        verbose = false, tolerance = 1e-10, bigfloat = false, λ = 0,
        mincurrent = 0)

    verbose && printlnflush("Looking for cells in the graph...")
    cells = find_cells(g, cells = initialcells)

    # here the system's matrix is constructed
    verbose && printlnflush("Calculating the system's matrix...")
    M = system_matrix(poi, vertex_positions, cells)

    # the goal field values at each of poi
    # x, y and z are simply concatenated in an xyzxyzxyz way
    Bpoi = vcat(B.(poi)...)

    # calculate the optimal currents for the given goal field
    # \ returns the least-norm solution already (so least current) in case of ambiguity
    # We are using Tikchonov regularization
    verbose && printlnflush("Solving for currents in the cells...")
    optI = [M; λ * I] \ [Bpoi; zeros(size(M,2))]

    verbose && printlnflush("Finding the total currents along the edges...")
    edgecurrents = find_edgecurrents(g, cells, optI, bigfloat = bigfloat)

    verbose && printlnflush("Decomposing the solution into simple loops...")
    simpleloops, simpleloopscurrents = find_all_simpleloops(g, edgecurrents,
        tolerance = tolerance, mincurrent = mincurrent)
end

function cuboid_poi(voi_size = [1.0, 1.0, 1.0], voi_centre = [0.0, 0.0, 0.0],
        n = [9, 9, 9]; filled = true)
    voi_walls_positions = voi_size / 2
    voi_pos = linspace.(-voi_walls_positions, voi_walls_positions, n)
    thr = voi_walls_positions .- mean.(diff.(voi_pos)) / 2

    poi = [ [x, y, z] .+ voi_centre for x in voi_pos[1], y in voi_pos[2], z in voi_pos[3] if
            filled ||
            # faces
            ( ~(-thr[1] ≤ x ≤ thr[1]) && (-thr[2] ≤ y ≤ thr[2]) && (-thr[3] ≤ z ≤ thr[3]) ) ||
            ( ~(-thr[2] ≤ y ≤ thr[2]) && (-thr[1] ≤ x ≤ thr[1]) && (-thr[3] ≤ z ≤ thr[3]) ) ||
            ( ~(-thr[3] ≤ z ≤ thr[3]) && (-thr[1] ≤ x ≤ thr[1]) && (-thr[2] ≤ y ≤ thr[2]) ) ||
            # edges
            ( (-thr[1] ≤ x ≤ thr[1]) && ~(-thr[2] ≤ y ≤ thr[2]) && ~(-thr[3] ≤ z ≤ thr[3]) ) ||
            ( (-thr[2] ≤ y ≤ thr[2]) && ~(-thr[1] ≤ x ≤ thr[1]) && ~(-thr[3] ≤ z ≤ thr[3]) ) ||
            ( (-thr[3] ≤ z ≤ thr[3]) && ~(-thr[1] ≤ x ≤ thr[1]) && ~(-thr[2] ≤ y ≤ thr[2]) ) ||
            # corners
            ( ~(-thr[3] ≤ z ≤ thr[3]) && ~(-thr[1] ≤ x ≤ thr[1]) && ~(-thr[2] ≤ y ≤ thr[2]) ) ]
end

"""Decompose a number into given factors.
decompose(823.2, [100., 10., 1.]) = [8, 2, 3]
if roundlast is true, the last element is rounded:
decompose(84.86, [10., 1., 0.1], roundlast = true) = [8, 4, 9] """
function decompose(x::T, elements::Vector{T}; roundlast = false) where {T <: AbstractFloat}
    coeffs = similar(elements)

    for i in eachindex(elements)
        coeffs[i] = x
        for j in 1 : i-1
            # use prevfloat for numeric reasons (eg. 4.0 ÷ 0.2 = 19.0)
            # but 4.0 ÷ prevfloat(0.2) = 20.0
            coeffs[i] %= prevfloat(elements[j])
        end
        if roundlast && i == length(elements)
            coeffs[i] /= elements[i]
        else
            coeffs[i] ÷= prevfloat(elements[i])
        end
    end

    round.(Int, coeffs)
end

function decompose(x::Real, elements::AbstractArray; kwargs...)
    xt, elementst = promote([x], elements)
    decompose(xt[], elementst; kwargs...)
end

function decompose_currents(currents, elemcurrents)
    [ decompose(c, elemcurrents, roundlast = true) for c in currents ]
end

function real_decomposed_currents(simpleloopscurrents_decomp, elemcurrents)
    [ sum(decomp .* elemcurrents) for decomp in simpleloopscurrents_decomp ]
end

"Given the set of points that would make up a cell, order them into a proper cell."
function order_cell(g, face_unordered)
    face = [ face_unordered[1] ]

    while true
        for n in [inneighbors(g, face[end]); outneighbors(g, face[end])]
            n in face_unordered || continue
            n in face && continue
            push!(face, n)
            break
        end
        length(face) >= length(face_unordered) && break
    end

    push!(face, face[1])
    face
end

function decompose_string(x, coeffs, elements)
    "$(signif(Float64(x), 6)) = $(join([ "$(c)*$e" for (c, e) in zip(coeffs, elements)], " + "))"
end

function save_result(filename, simpleloops, simpleloopscurrents,
                     simpleloopscurrents_decomp, elemcurrents)
    # truncate the contents of the output file
    open(filename, "w") do file
        for i in eachindex(simpleloopscurrents)
            write(file, lpad(i, 4), "  :  ")
            write(file, decompose_string(simpleloopscurrents[i],
                                         simpleloopscurrents_decomp[i],
                                         elemcurrents))
            write(file, "  :  ")
            write(file, repr(simpleloops[i]))
            write(file, "\n")
        end
    end
end

end


text:
using Astro, Test

fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
quiet = length(ARGS) > 0 && ARGS[1] == "-q"

anyerrors = false

push!(LOAD_PATH, "src")

loaddebugflag = false

my_tests = String[
    "calendar-tests.jl",
    "coordinates-test.jl",
    "dynamical-tests.jl",
    "earth-tests.jl",
    "equinox-tests.jl",
    "kepler-tests.jl",
    "moon-tests.jl",
    "nutation-tests.jl",
    "planet-tests.jl",
    "pluto-tests.jl",
    "riseset-tests.jl",
    "sidereal-tests.jl",
    "sun-tests.jl",
    "utils-tests.jl"
]

println("...Running tests:")

passed = String[]
failed = String[]

for my_test in my_tests
    try
        include(my_test)
        println("$(my_test)")
        push!(passed, my_test)
    catch e
        global anyerrors = true
        push!(failed, my_test)
        println("$(my_test)")
        if fatalerrors
            rethrow(e)
        elseif !quiet
            showerror(stdout, e, backtrace())
            println()
        end
    end
end

for tst in passed
    printstyled("PASSED:\t", color=:green)
    println("$tst")
end

for tst in failed
    printstyled("FAILED:\t", color=:red)
    println("$tst")
end


text:
## test_show.jl : test show commands

@info "Testing display commands"

@show Matern(1,1)
@show Exponential(1)
@show SquaredExponential(1)
@show AnisotropicExponential([1000 0; 0 1000])
@show CovarianceFunction(2,Matern(1.0,2))
@show Cholesky()
@show Spectral()
@show KarhunenLoeve(10)
@show CirculantEmbedding()
@show GaussianRandomField(CovarianceFunction(1,Exponential(1)),Cholesky(),0:0.1:1)
@show SeparableCovarianceFunction(Exponential(1),Exponential(1))

@info " Done testing display commands"


text:
function (F::AForcing)(s::AbstractSimulation)

    cPlaneNumber = 1
    dt = get_dt(s)
    Tf = F.Tf
    ts = 10*dt
    omega = (2*π/(Tf*ts))
    alpha = F.α
    eps = Base.eps()
    _kf = F.Kf

    R = F.R
    Em = F.Em
    Zf = F.Zf
    factor = F.factor

    s1 = fill!(F.forcex,0)
    s2 = fill!(F.forcey,0)

    u1 = s.u.c.x
    u2 = s.u.c.y
    u3 = s.u.c.z

    @inbounds for j=YRANGE
        for i=XRANGE
            k = fsqrt(muladd(KX[i],KX[i],KY[j]^2))
            if (0 < k <= _kf)
                Ef = 0.5*(1 + (i>1))*(abs2(u1[i,j,1]) + abs2(u2[i,j,1]))
                R[i,j] += dt*(-2*alpha*omega*R[i,j] - omega*omega*(Ef - Em[i,j]))
                R[i,j] = max(0.0, R[i,j])
                factor[i,j] = fsqrt(R[i,j]/max(Ef,eps))*Zf[i,j]*dt
            end
        end
    end

    umag = 0.0
  #  # -------------------- Horizontal Forcing --------------------
    ff = 0.0
    fv = 0.0

    @inbounds for j=YRANGE
        conjFactX=1.0
        for i=XRANGE
            k = fsqrt(muladd(KX[i],KX[i],KY[j]^2))
            if (0 < k <= _kf)
                s1[i,j,1] = u1[i,j,1]*factor[i,j]
                s2[i,j,1] = u2[i,j,1]*factor[i,j]
                ff += (abs2(s1[i,j,1]) + abs2(s2[i,j,1]))*conjFactX
                fv += (proj(u1[i,j,1],s1[i,j,1]) + proj(u2[i,j,1],s2[i,j,1]))*conjFactX
            end
            conjFactX=2.0
        end
    end
    ff = ff*0.5/dt
    fv = fv/dt
    horizontalPower = ff+fv
    F.hp[] = horizontalPower
    powerFraction = 0.99

    # # -------------------- Vertical Forcing --------------------
    # /*
    #   Force the points with kx=ky=0 and kz=1,2,3,4,5 so that the power
    #   added is (1-powerFraction) of that added to the horizontal */
    ff = 0.0
    fv = 0.0
    umag = 0.0
    #     // get the sum of the existing velocites
    @inbounds for k = 2:6
        umag += fsqrt(abs2(u1[1,1,k]) + abs2(u2[1,1,k])) + fsqrt(abs2(u1[1,1,NZ-k+2]) + abs2(u2[1,1,NZ-k+2]))
    end

    f = dt*sqrt(2.0*horizontalPower*(1.0-powerFraction))/umag

    @inbounds for k=2:6
        s1[1,1,k] = u1[1,1,k]*f
        s2[1,1,k] = u2[1,1,k]*f
        s1[1,1,NZ-k+2] = u1[1,1,NZ-k+2]*f
        s2[1,1,NZ-k+2] = u2[1,1,NZ-k+2]*f
        ff += abs2(s1[1,1,k]) + abs2(s2[1,1,k]) + abs2(s1[1,1,NZ-k+2]) + abs2(s2[1,1,NZ-k+2])
        fv += proj(u1[1,1,k],s1[1,1,k]) + proj(u2[1,1,k],s2[1,1,k]) + proj(u1[1,1,NZ-k+2],s1[1,1,NZ-k+2]) + proj(u2[1,1,NZ-k+2],s2[1,1,NZ-k+2])
    end
    ff = ff*0.5/dt
    fv = fv/dt
    verticalPower = ff+fv
    F.vp[] = verticalPower
    return 1/dt
end


text:
@testset "Convert" begin

    x = variable()
    for (α, β) in ((1/2, 1/4), (1/4, 1/2), (1,2), (2,1))
        P, Q = Laguerre{α}, Laguerre{β}
        for _ in 1:10
            p = P(rand(1:5,5)); q=P(rand(1:5,5))
            @test convert(P, convert(Q, p))(x) ≈ p(x)
            @test convert(Q, convert(P,q))(x) ≈ q(x)
        end
    end

end


text:
# 3.6 profile the economies under the 5 (including benchmark) scenarios
# NOTE: the crowding out effect
idx_plot = (2015:2025) .- env.START_YEAR .+ 1
plot_symbol = :K
tmpexpr = :(
   xlabel("Year"); grid(true);
   xlim([ idx_plot[1] - 2 + env.START_YEAR, idx_plot[end] + env.START_YEAR ]);
)
figure(figsize = (9,5)) # (C+I)/G
      for tmpreform in 1:length(res_IndiInDesc["reformNames"])
         tmpploty = (res_IndiInDesc["Dt"][tmpreform][:C] .+ res_IndiInDesc["Dt"][tmpreform][:I]  ) ./ res_IndiInDesc["Dt"][tmpreform][:G]
         plot( Dt[:Year][idx_plot], 100 .* tmpploty[idx_plot], res_IndiInDesc["plotstyles"][tmpreform] )
      end
      eval(tmpexpr); ylabel("(C+I)/G (%)");
      legend( res_IndiInDesc["plotlegends"] , loc = "best")
# 3.7 save figure
savefig( "./output/CrowdingOutEffect.pdf", format = "pdf" )

#


text:
function f(x) 42 end
function f(x) 1729 end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using DashBootstrapComponents

pagination = html_div([
    html_div("Small"),
    dbc_pagination(max_value = 5, size = "sm"),
    html_div("Default"),
    dbc_pagination(max_value = 5),
    html_div("Large"),
    dbc_pagination(max_value = 5, size = "lg"),
])

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using Test, Unicode
Base.include(@__MODULE__, joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testhelpers", "FakePTYs.jl"))
import .FakePTYs: open_fake_pty

slave, master = open_fake_pty()

CTRL_C = '\x03'

# Script that we want the REPL to execute, here simply st for Pkg REPLMode
test_script1 = """
using ReplMaker

function parse_to_expr(s)
    quote Meta.parse(\$s) end
end

initrepl(parse_to_expr,
         prompt_text="Expr> ",
         prompt_color = :blue,
         start_key=')',
         mode_name="Expr_mode");

) x + 1

"""*CTRL_C

test_script2 = """
using ReplMaker

function parse_to_expr(s)
    quote Meta.parse(\$s) end
end

mode = initrepl(parse_to_expr,
         prompt_text="Expr> ",
         prompt_color = :blue,
         start_key=')',
         mode_name="Expr_mode");

enter_mode!(mode);
x + 1

"""*CTRL_C

function run_test(test_script)
    slave, master = open_fake_pty()
    # Start a julia process
    p = run(`$(Base.julia_cmd()) --history-file=no --startup-file=no`, slave, slave, slave; wait=false)

    # Read until the prompt
    readuntil(master, "julia>", keep=true)
    done = false
    repl_output_buffer = IOBuffer()

    # A task that just keeps reading the output
    @async begin
        while true
            done && break
            write(repl_output_buffer, readavailable(master))
        end
    end

    # Execute our "script"
    for l in split(test_script, '\n'; keepempty=false)
        write(master, l, '\n')
    end

    # Let the REPL exit
    write(master, "exit()\n")
    wait(p)
    done = true

    # Gather the output
    repl_output = String(take!(repl_output_buffer))
    println(repl_output)
    return split(repl_output, '\n'; keepempty=false)
end

# to avoid compilations when tests get started,
# just compile it here.
run_test("using ReplMaker"*CTRL_C)

@testset "test opening REPL modes manually" begin
    out1 = run_test(test_script1);
    @test out1[end-5] == "\e[?2004h\r\e[0K\e[34m\e[1mExpr> \e[0m\e[0m\r\e[6C\r\e[6C\r\e[0K\e[34m\e[1mExpr> \e[0m\e[0m\r\e[6C\r\e[6C^C\r"
end

@testset "test opening REPL modes automatically" begin
    out2 = run_test(test_script2);
    @test out2[end-5] == "\e[?2004h\r\e[0K\e[34m\e[1mExpr> \e[0m\e[0m\r\e[6C\r\e[6C\r\e[0K\e[34m\e[1mExpr> \e[0m\e[0m\r\e[6C\r\e[6C^C\r"
end


text:
# This file is adapted from test/abstractarray.jl from Julia.

# This file is a part of Julia. License is MIT: https://julialang.org/license

A_abs = rand(5,4,3)
A = OpenCV.Mat(A_abs)
@testset "Bounds checking" begin
    @test checkbounds(Bool, A, 1, 1, 1) == true
    @test checkbounds(Bool, A, 5, 4, 3) == true
    @test checkbounds(Bool, A, 0, 1, 1) == false
    @test checkbounds(Bool, A, 1, 0, 1) == false
    @test checkbounds(Bool, A, 1, 1, 0) == false
    @test checkbounds(Bool, A, 6, 4, 3) == false
    @test checkbounds(Bool, A, 5, 5, 3) == false
    @test checkbounds(Bool, A, 5, 4, 4) == false
    @test checkbounds(Bool, A, 1) == true           # linear indexing
    @test checkbounds(Bool, A, 60) == true
    @test checkbounds(Bool, A, 61) == false
    @test checkbounds(Bool, A, 2, 2, 2, 1) == true  # extra indices
    @test checkbounds(Bool, A, 2, 2, 2, 2) == false
    @test checkbounds(Bool, A, 1, 1)  == false
    @test checkbounds(Bool, A, 1, 12) == false
    @test checkbounds(Bool, A, 5, 12) == false
    @test checkbounds(Bool, A, 1, 13) == false
    @test checkbounds(Bool, A, 6, 12) == false
end

@testset "single CartesianIndex" begin
    @test checkbounds(Bool, A, CartesianIndex((1, 1, 1))) == true
    @test checkbounds(Bool, A, CartesianIndex((5, 4, 3))) == true
    @test checkbounds(Bool, A, CartesianIndex((0, 1, 1))) == false
    @test checkbounds(Bool, A, CartesianIndex((1, 0, 1))) == false
    @test checkbounds(Bool, A, CartesianIndex((1, 1, 0))) == false
    @test checkbounds(Bool, A, CartesianIndex((6, 4, 3))) == false
    @test checkbounds(Bool, A, CartesianIndex((5, 5, 3))) == false
    @test checkbounds(Bool, A, CartesianIndex((5, 4, 4))) == false
    @test checkbounds(Bool, A, CartesianIndex((1,))) == false
    @test checkbounds(Bool, A, CartesianIndex((60,))) == false
    @test checkbounds(Bool, A, CartesianIndex((61,))) == false
    @test checkbounds(Bool, A, CartesianIndex((2, 2, 2, 1,))) == true
    @test checkbounds(Bool, A, CartesianIndex((2, 2, 2, 2,))) == false
    @test checkbounds(Bool, A, CartesianIndex((1, 1,)))  == false
    @test checkbounds(Bool, A, CartesianIndex((1, 12,))) == false
    @test checkbounds(Bool, A, CartesianIndex((5, 12,))) == false
    @test checkbounds(Bool, A, CartesianIndex((1, 13,))) == false
    @test checkbounds(Bool, A, CartesianIndex((6, 12,))) == false
end

@testset "mix of CartesianIndex and Int" begin
    @test checkbounds(Bool, A, CartesianIndex((1,)), 1, CartesianIndex((1,))) == true
    @test checkbounds(Bool, A, CartesianIndex((5, 4)), 3)  == true
    @test checkbounds(Bool, A, CartesianIndex((0, 1)), 1)  == false
    @test checkbounds(Bool, A, 1, CartesianIndex((0, 1)))  == false
    @test checkbounds(Bool, A, 1, 1, CartesianIndex((0,))) == false
    @test checkbounds(Bool, A, 6, CartesianIndex((4, 3)))  == false
    @test checkbounds(Bool, A, 5, CartesianIndex((5,)), 3) == false
    @test checkbounds(Bool, A, CartesianIndex((5,)), CartesianIndex((4,)), CartesianIndex((4,)))  == false
end

@testset "vector indices" begin
    @test checkbounds(Bool, A, 1:5, 1:4, 1:3) == true
    @test checkbounds(Bool, A, 0:5, 1:4, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 0:4, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 1:4, 0:3) == false
    @test checkbounds(Bool, A, 1:6, 1:4, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 1:5, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 1:4, 1:4) == false
    @test checkbounds(Bool, A, 1:60) == true
    @test checkbounds(Bool, A, 1:61) == false
    @test checkbounds(Bool, A, 2, 2, 2, 1:1) == true  # extra indices
    @test checkbounds(Bool, A, 2, 2, 2, 1:2) == false
    @test checkbounds(Bool, A, 1:5, 1:4) == false
    @test checkbounds(Bool, A, 1:5, 1:12) == false
    @test checkbounds(Bool, A, 1:5, 1:13) == false
    @test checkbounds(Bool, A, 1:6, 1:12) == false
end

@testset "logical" begin
    @test checkbounds(Bool, A, trues(5), trues(4), trues(3)) == true
    @test checkbounds(Bool, A, trues(6), trues(4), trues(3)) == false
    @test checkbounds(Bool, A, trues(5), trues(5), trues(3)) == false
    @test checkbounds(Bool, A, trues(5), trues(4), trues(4)) == false
    @test checkbounds(Bool, A, trues(60)) == true
    @test checkbounds(Bool, A, trues(61)) == false
    @test checkbounds(Bool, A, 2, 2, 2, trues(1)) == true  # extra indices
    @test checkbounds(Bool, A, 2, 2, 2, trues(2)) == false
    @test checkbounds(Bool, A, trues(5), trues(12)) == false
    @test checkbounds(Bool, A, trues(5), trues(13)) == false
    @test checkbounds(Bool, A, trues(6), trues(12)) == false
    @test checkbounds(Bool, A, trues(5, 4, 3)) == true
    @test checkbounds(Bool, A, trues(5, 4, 2)) == false
    @test checkbounds(Bool, A, trues(5, 12)) == false
    @test checkbounds(Bool, A, trues(1, 5), trues(1, 4, 1), trues(1, 1, 3)) == false
    @test checkbounds(Bool, A, trues(1, 5), trues(1, 4, 1), trues(1, 1, 2)) == false
    @test checkbounds(Bool, A, trues(1, 5), trues(1, 5, 1), trues(1, 1, 3)) == false
    @test checkbounds(Bool, A, trues(1, 5), :, 2) == false
end

@testset "array of CartesianIndex" begin
    @test checkbounds(Bool, A, [CartesianIndex((1, 1, 1))]) == true
    @test checkbounds(Bool, A, [CartesianIndex((5, 4, 3))]) == true
    @test checkbounds(Bool, A, [CartesianIndex((0, 1, 1))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 0, 1))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 1, 0))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((6, 4, 3))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 5, 3))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 4, 4))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 1))], 1) == true
    @test checkbounds(Bool, A, [CartesianIndex((5, 4))], 3) == true
    @test checkbounds(Bool, A, [CartesianIndex((0, 1))], 1) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 0))], 1) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 1))], 0) == false
    @test checkbounds(Bool, A, [CartesianIndex((6, 4))], 3) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 5))], 3) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 4))], 4) == false
end

println("OpenCV.Mat tests passed")


text:
module TestStaticArrays
using Test
using Setfield
using StaticArrays

@testset "StaticArrays" begin
    obj = StaticArrays.@SMatrix [1 2; 3 4]
    @testset for l in [
            (@lens _[2,1]),
            (@lens _[$2,$1]),
        ]
        @test get(obj, l) == 3
        @test set(obj, l, 5) == StaticArrays.@SMatrix [1 2; 5 4]
        @test setindex(obj, 5, 2, 1) == StaticArrays.@SMatrix [1 2; 5 4]
    end

    v = @SVector [1,2,3]
    @test (@set v[1] = 10) === @SVector [10,2,3]
    @test_broken (@set v[1] = π) === @SVector [π,2,3]
end
end


text:
#
#  7 ---9--- 8 ---12-- 9
#
#  |         |         |
#  |         |         |
#  |         |         |
#  10   4    11   5    13
#  |         |         |
#  |         |         |
#  |         |         |
#
#  4 ---2--- 5 ---8--- 6
#
#  |         | \   3   |
#  |         |  \      |
#  |         |   \     |
#  3    1    4    6    7
#  |         |     \   |
#  |         |  2   \  |
#  |         |       \ |
#
#  1 ---1--- 2 ---5--- 3
#

struct GridMock <: Grid{2,2} end

function get_node_coordinates(::GridMock)
  Point{2,Float64}[(0,0),(1,0),(2,0),(0,1),(1,1),(2,1),(0,2),(1,2),(2,2)]
end

function get_cell_nodes(::GridMock)
                          #1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8
  cell_to_vertices_data = [1,2,4,5,2,3,5,3,6,5,4,5,7,8,5,6,8,9]
  cell_to_vertices_ptrs = Int32[1,      5,    8,    11,     15,    19]
  Table(cell_to_vertices_data,cell_to_vertices_ptrs)
end

function get_reffes(::GridMock)
  [QUAD4, TRI3]
end

function get_cell_type(::GridMock)
  quad = Int8(1)
  tri = Int8(2)
  [quad, tri, tri, quad, quad]
end


text:
using Flux
using Flux: gradient, logitcrossentropy, params

"""
    fit(client, weights)

Train the clients with the given weigths and return the updated weights.
"""
function fit(client::Client, weights::Vector{Float32})::Vector{Float32}
    @info "fit on client..."
    load_flatten_model!(client.model, weights)

    loss(x, y) = logitcrossentropy(client.model(x), y)
    epoch_loss = Vector{Float32}(undef, client.epochs)

    @inbounds for i = 1:client.epochs
        local train_loss = 0.0f0

        for batch in client.dataset
            grads = gradient(params(client.model)) do
                train_loss += loss(batch...)
            end
            Flux.update!(client.optimizer, params(client.model), grads)
        end

        epoch_loss[i] = train_loss / length(client.dataset)
    end

    loss_total = sum(epoch_loss) / length(epoch_loss)
    @info "loss = $loss_total"

    return flatten_model(client.model)
end


text:
"""
    ht_p2_p1_veclap_alt

The manufactured-solution colliding flow example from Elman et al 2014. The
Hood-Taylor formulation with quadratic triangles for the velocity and continuous
pressure on linear triangles.

This implementation is an alternative: for the velocity, a single finite element
space with multiple components (2) is used instead of multiple finite element
spaces.

The formulation is the one derived in Donea, Huerta, Introduction to the finite element
method, 1993. Page 486 ff, and Elman, et al., Finite elements and fast
iterative solvers, p. 132. In other words, it is the vector Laplacian version.
"""
module ht_p2_p1_veclap_alt

using LinearAlgebra
using StaticArrays
using MeshCore: retrieve, nrelations, nentities, ir_identity, attribute, VecAttrib
using MeshSteward: T6block, T6toT3
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary, transform
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T6, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs, eldofvals
using Elfel.FEIterators: jacjac, location
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using UnicodePlots

mu = 1.0 # dynamic viscosity
A = 1.0 # half of the length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # Hood-Taylor pair of meshes is needed
    # This mesh will be for the velocities
    vmesh = Mesh()
    attach!(vmesh, T6block(2 * A, 2 * A, N, N), "velocity")
    ir = baseincrel(vmesh)
    transform(ir, x -> x .- A)
    # This mesh will be used for the pressures. Notice that it needs to be
    # "compatible" with the velocity mesh in the sense that they need to share
    # the nodes at the corners of the triangles.
    pmesh = Mesh()
    attach!(pmesh, T6toT3(baseincrel(vmesh, "velocity")), "pressure")
    return vmesh, pmesh
end

function assembleK(ufesp, pfesp, tndof, mu)
    function integrate!(ass, elits, qpits, mu)
        unedof, pnedof = ndofsperel.(elits)
        uedofcomp = edofcompnt(ufesp)
        kuu = LocalMatrixAssembler(unedof, unedof, 0.0)
        kup = LocalMatrixAssembler(unedof, pnedof, 0.0)
        for el in zip(elits...)
            uel, pel = el
            init!(kuu, eldofs(uel), eldofs(uel))
            init!(kup, eldofs(uel), eldofs(pel))
            for qp in zip(qpits...)
                uqp, pqp = qp
                Jac, J = jacjac(uel, uqp)
                JxW = J * weight(uqp)
                gradNu = bfungrad(uqp, Jac)
                Np = bfun(pqp)
                for j in 1:unedof, i in 1:unedof
                    if uedofcomp[i] == uedofcomp[j]
                        kuu[i, j] += (mu * JxW) * (dot(gradNu[i], gradNu[j]))
                    end
                end
                for j in 1:pnedof, i in 1:unedof
                    kup[i, j] += (-JxW * Np[j]) * gradNu[i][uedofcomp[i]]
                end
            end
            assemble!(ass, kuu)
            assemble!(ass, kup)
            assemble!(ass, transpose(kup))
        end
        return ass
    end

    elits = (FEIterator(ufesp), FEIterator(pfesp))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(ufesp, qargs), QPIterator(pfesp, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrate!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function evaluate_pressure_error(pfesp)
    function integrate!(elit, qpit, truep)
        pnedof = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:pnedof
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(pfesp)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(pfesp, qargs)
    return integrate!(elit, qpit, truep)
end

function evaluate_velocity_error(ufesp)
    function integrate!(elit, qpit, trueux, trueuy)
        unedof = ndofsperel(elit)
        uedofcomp = edofcompnt(ufesp)
        E = 0.0
        for el in elit
            udofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Nu = bfun(qp)
                uxt = trueux(location(el, qp)...)
                uyt = trueuy(location(el, qp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:unedof
                    (uedofcomp[j] == 1) && (uxa += (udofvals[j] * Nu[j]))
                    (uedofcomp[j] == 2) && (uya += (udofvals[j] * Nu[j]))
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(ufesp)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(ufesp, qargs)
    return integrate!(elit, qpit, trueux, trueuy)
end

function run(N)
    vmesh, pmesh = genmesh(N)
    # Velocity spaces
    ufesp = FESpace(Float64, vmesh, FEH1_T6(), 2)
    locs = geometry(vmesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(ufesp, 0, i, 1, trueux(locs[i]...))
            setebc!(ufesp, 0, i, 2, trueuy(locs[i]...))
        end
    end
    # Pressure space
    pfesp = FESpace(Float64, pmesh, FEH1_T3(), 1)
    atcenter = vselect(geometry(pmesh); nearestto = [0.0, 0.0])
    setebc!(pfesp, 0, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([ufesp, pfesp])
    tndof = ndofs(ufesp) + ndofs(pfesp)
    tnunk = nunknowns(ufesp) + nunknowns(pfesp)
    # Assemble the coefficient matrix
    K = assembleK(ufesp, pfesp, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [ufesp, pfesp])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([ufesp, pfesp], U)
    # Postprocessing
    makeattribute(pfesp, "p", 1)
    makeattribute(ufesp, "ux", 1)
    makeattribute(ufesp, "uy", 2)
    ep = evaluate_pressure_error(pfesp)
    ev = evaluate_velocity_error(ufesp)
    vtkwrite("ht_p2_p1_veclap_alt-p", baseincrel(pmesh), [(name = "p",), ])
    vtkwrite("ht_p2_p1_veclap_alt-v", baseincrel(vmesh), [(name = "ux",), (name = "uy",)])
    return (ep, ev)
end

end

using .ht_p2_p1_veclap_alt
let
    N = 4
    for loop in 1:5
        ep, ev = ht_p2_p1_veclap_alt.run(N)
        @show N, ep, ev
        N = N * 2
    end
end


text:
"""
Package de fonction pour la minimisation ActifMPCC

liste des fonctions :
LineSearchSolve(ma::ActifMPCCmod.MPCC_actif,xj::Vector,hd::Any;scaling :: Bool = false)
"""

module UnconstrainedMPCCActif

using ActifMPCCmod
using NLPModels
using OutputLSmod

"""
Input :
ma : MPCC_actif
xj : vecteur initial
d : direction précédente en version étendue
"""

function LineSearchSolve(ma::ActifMPCCmod.MPCC_actif,
                         xj::Vector,d::Vector,
                         step::Float64,gradpen::Vector,hg::Float64;
                         scaling :: Bool = false)

 output=0
 d=ActifMPCCmod.redd(ma,d)
 scale=1.0
 beta=ma.beta

 #xj est un vecteur de taille (n x length(bar_w))
 if length(xj) == (ma.n+2*ma.nb_comp)
  xj=vcat(xj[1:ma.n],xj[ma.n+ma.w13c],xj[ma.n+ma.nb_comp+ma.w24c])
 elseif length(xj) != (ma.n+length(ma.w13c)+length(ma.w24c))
  println("Error dimension : UnconstrainedSolve")
  return
 end

 #Choix de la direction :

 gradf=ActifMPCCmod.grad(ma,xj,gradpen)
 gradft=Array{Float64,1}

 #Calcul d'une direction de descente de taille (n + length(bar_w))
 d=ma.direction(ma,gradf,xj,d,beta)

 slope = dot(gradf,d)

 if slope > 0.0  # restart with negative gradient
  d = - gradf
  slope =  dot(gradf,d)
 end

 #Calcul du pas maximum (peut être infinie!)
 stepmax,wmax,wnew = ActifMPCCmod.PasMax(ma,xj,d)

 step,good_grad,ht,nbarmijo,nbwolfe,gradft=ma.linesearch(ma,xj,d,hg,stepmax,scale*slope,step)

 step*=scale

 ols=OutputLSmod.OutputLS(stepmax,step,slope,beta,nbarmijo,nbwolfe)

 xjp=xj+step*d

 sol = ActifMPCCmod.evalx(ma,xjp)
 dsol = ActifMPCCmod.evald(ma,d)

 gradpen=NLPModels.grad(ma.nlp,sol)
 good_grad || (gradft=ActifMPCCmod.grad(ma,xjp,gradpen))

 s=xjp-xj
 y=gradft-gradf

 #MAJ des paramètres du calcul de la direction
 ma=ma.direction(ma,xjp,beta,gradft,gradf,y,d,step)

 #MAJ du scaling
 if scaling
  scale=dot(y,s)/dot(y,y)
 end
 if scale <= 0.0
  scale=1.0
 end

 small_step=norm(xjp-xj,Inf)<=eps(Float64)?true:false #on fait un pas trop petit

 #si alpha=pas maximum alors on met w à jour.
 if stepmax == step
  ActifMPCCmod.setw(ma,wmax)
 else
  wnew = zeros(Bool,0,0) #si on a pas pris le stepmax alors aucune contrainte n'est ajouté
 end
 if nbarmijo >= ma.paramset.ite_max_armijo || true in isnan.(sol) || true in isnan(ht)
  output=1
 elseif nbwolfe==ma.paramset.ite_max_wolfe
  output=2;
  xjp=xj
 end

 return sol,ma.w,dsol,step,wnew,output,small_step,ols,gradpen,ht
end

#end of module
end


text:
# Compare the selection accuracy of KS test and BIC for different sample size
@everywhere using GaussianMixtureTest: gmm, kstest
@everywhere using Distributions:Normal, MixtureModel
@everywhere using Random:seed!

@everywhere function compareBIC(Ctrue::Int, n::Int, b::Int)

    C_max = max(5, (2*Ctrue - 1))
    mu_all = log(1/0.779 - 1)
    if Ctrue == 2
        mu_true = [mu_all - 1.0, mu_all + 0.8]
        wi_true =  [.6, .4]
        sigmas_true = [1.2, .8]
    elseif Ctrue == 3
        mu_true = [mu_all - 2.0, mu_all + 1.0, mu_all + 3.5]
        wi_true = [.3, .4, .3]
        sigmas_true = [1.2, .8, .9]
    end
    m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)

    C_max = max(5, (2*Ctrue - 1))
    bic = fill(Inf, C_max)
    pvalues = fill(1.0, C_max)

    seed!(b)
    x = rand(m, n);
    for C in 1:C_max
        wi, mu, sigmas, ml = gmm(x, C)
        T, pvalues[C] = kstest(x, C)
        bic[C] = -2 * ml + (3 * C - 1) * log(n)
    end
    vcat(bic, pvalues)
end