text:
module EmpiricalLikelihood

using DiffResults
using ForwardDiff
using JuMP
using Optim
using PyCall
using Statistics

include("HigherOrderKernels.jl")

const scipy_spatial = PyNULL()

function __init__()
    copy!(scipy_spatial, pyimport_conda("scipy.spatial", "scipy"))
end

export optimize_el, setup_jump_problem

const Kernel = EpanechnikovKernel{2}

function optimize_el(n_moments::Integer, moment_func, X::Matrix{T}, Y::Matrix{T}, start_param::Vector{T})::Tuple{Vector{T}, T} where T <: AbstractFloat
    # Automatic bandwith selection according to Silverman's rule of thumb.
    # This is not the best we can do, e.g. in the paper it is suggested to use
    # cross-validation a la Newey (1993). However, since according to the paper
    # the bandwith is not very important, we just use something simple.
    # So we take the average bandwith over the bandwith's suggested
    b = mean(density_bandwidth(Kernel, X[:, i]) for i in 1:size(X, 2))
    return optimize_el(n_moments, moment_func, X, Y, start_param, b)
end

function optimize_el(n_moments::Int, moment_func, X::Matrix{T}, Y::Matrix{T}, start_param::Vector{T}, bandwidth::T, trimming::T=.99999)::Tuple{Vector{T}, T} where T <: AbstractFloat
    return optimize_el(convert(UInt, n_moments), moment_func, X, Y, start_param, bandwidth, trimming)
end

function optimize_el(n_moments::UInt, moment_func, X::Matrix{T}, Y::Matrix{T}, start_param::Vector{T}, bandwidth::T, trimming::T=.99999)::Tuple{Vector{T}, T} where T <: AbstractFloat
    weights, not_trimmed = weights_and_trims(n_moments, moment_func, X, Y, start_param, bandwidth, trimming)

    res = Optim.optimize(
        Optim.only_fg!((f, g, x) -> el_objective!(n_moments, moment_func, X, Y, weights, not_trimmed, x, f, g)),
        start_param,
        LBFGS()
    )

    return Optim.minimizer(res), -1 * Optim.minimum(res)
end

function setup_jump_problem(n_moments::Int, moment_func, X::Matrix{T}, Y::Matrix{T}, start_param::Vector{T}, bandwidth::T, trimming::T=.99999) where T <: AbstractFloat
    return setup_jump_problem(convert(UInt, n_moments), moment_func, X, Y, start_param, bandwidth, trimming)
end

function setup_jump_problem(n_moments::UInt, moment_func, X::Matrix{T}, Y::Matrix{T}, start_param::Vector{T}, bandwidth::T, trimming::T=.99999) where T <: AbstractFloat
    weights, not_trimmed = weights_and_trims(n_moments, moment_func, X, Y, start_param, bandwidth, trimming)

    model = Model()
    n_params = length(start_param)
    @variable(model, theta[1:n_params])
    for i in 1:n_params
        set_start_value(theta[i], start_param[i])
    end
    neg_ll(x...) = begin
        x_vec = vcat(x...)
        el_objective!(n_moments, moment_func, X, Y, weights, not_trimmed, x_vec, 1, nothing)
    end
    ∇neg_ll(grad, x...) = begin
        x_vec = vcat(x...)
        el_objective!(n_moments, moment_func, X, Y, weights, not_trimmed, x_vec, nothing, grad)
    end
    register(model, :neg_empirical_likelihood, n_params, neg_ll, ∇neg_ll)
    @NLobjective(model, Min, neg_empirical_likelihood(theta...))

    return model, neg_ll, ∇neg_ll
end

function weights_and_trims(n_moments::UInt, moment_func, X::Matrix{T}, Y::Matrix{T}, start_param::Vector{T}, bandwidth::T, trimming::T)::Tuple{Matrix{T}, Vector{Bool}} where T <: AbstractFloat
    # Calcualte weights
    n_obs, n_conditioning_vars = size(X)
    weights = Matrix{T}(undef, n_obs, n_obs)
    not_trimmed = Vector{Bool}(undef, n_obs)  # false = we will trim
    for i in 1:n_obs
        for j in 1:n_obs
            weights[i, j] = prod(
                density_kernel(Kernel, (X[i, k] - X[j, k]) / bandwidth)
                for k in 1:n_conditioning_vars
            )
        end

        # Determine if the observation will be trimmed.
        # Contrary to what is claimed in the paper, I find it impossible to get the
        # estimator even working without trimming. (Just consider the case where w_{ij} =1,
        # i.e. an observation without "neighbours". In these cases the inner loop is
        # always ill-defined.)
        # Nadaray-Watson estimator of pdf of X at x_i
        hhat = sum(weights[i, :]) / (n_obs * bandwidth^n_conditioning_vars)
        not_trimmed[i] = hhat >= bandwidth^trimming

        # Normalize weights
        weights[i, :] ./= sum(weights[i, :])
    end

    if !any(not_trimmed)
        throw("All observations are trimmed. Decrease the bandwith or increase the trimming parameter.")
    end

    # Check if starting paramater is valid
    moment_vals, _ = moment_vals_and_jac(n_moments, moment_func, X, Y, start_param)
    if !zero_in_hull_for_all(moment_vals, weights, not_trimmed)
        throw("Zero is not in convex hull of moments at starting value. Try a different starting value or increase the bandwidth.")
    end

    return weights, not_trimmed
end

function el_objective!(n_moments::UInt, moment_func, X::Matrix{S}, Y::Matrix{S}, weights::Matrix{S}, not_trimmed::Vector{Bool}, theta::Vector{T}, out, grad::Union{Nothing, AbstractVector{S}})::Union{Nothing, T} where {T <: Real, S <: AbstractFloat}
    moment_vals, moment_jac = moment_vals_and_jac(n_moments, moment_func, X, Y, theta)

    if !zero_in_hull_for_all(moment_vals, weights, not_trimmed)
        if !isnothing(grad)
            grad[:] .= NaN
        end
        if !isnothing(out)
            return -Inf
        end
        return
    end

    # Inner loop
    n_obs = size(X, 1)
    n_params = length(theta)
    λ = Matrix{T}(undef, n_obs, n_moments)
    # Obj contains the negative of the objective at the top of p. 1674, as we will minimize and not maximize
    obj = 0.
    if grad != nothing
        grad[:] = zeros(T, n_params)
    end
    for i in 1:n_obs
        if !not_trimmed[i]
            # This adds 0 to the SER, just go on
            continue
        end

        # Find λ maximizing as in (2.7) in the paper
        weight_i = weights[i, :]
        res = Optim.optimize(
            Optim.only_fgh!((f, g, h, γ) -> λ_obj_fgh(f, g, h, γ, moment_vals, weight_i)),
            # γ -> -1 * λ_obj(γ, moment_vals, weight_i),
            # (g, γ) -> -1 * λ_grad!(g, γ, moment_vals, weight_i),
            zeros(n_moments),
            Newton()
        )
        if !Optim.converged(res)
            throw("Inner loop did not converge")
        end
        λ[i, :] = Optim.minimizer(res)
        # As in the paper, we will not trim
        #obj += λ_obj(λ[i, :], moment_vals, weight_i)
        obj += -1 * Optim.minimum(res)
        if grad != nothing
            # Use envelope theorem to obtain gradient
            grad[:] += sum(weight_i[j] * λ[i, :] * moment_jac[j, :, :] / (1 + transpose(λ[i, :]) * moment_vals[j, :]) for j in 1:n_obs, dims=1)[1, :]
        end
    end

    if out != nothing
        return obj
    end
end

function moment_vals_and_jac(n_moments::UInt, moment_func, X::Matrix{T}, Y::Matrix{T}, theta::Vector{T})::Tuple{Matrix{T}, Array{T, 3}} where T <: AbstractFloat
    n_obs = size(X, 1)
    n_params = length(theta)
    moment_vals = Matrix{T}(undef, n_obs, n_moments)
    moment_jac = Array{T, 3}(undef, n_obs, n_moments, n_params)
    out = DiffResults.JacobianResult(zeros(n_moments), theta)
    for i in 1:n_obs
        ForwardDiff.jacobian!(out, x -> moment_func(X[i, :], Y[i, :], x), theta)
        moment_vals[i, :] = DiffResults.value(out)
        moment_jac[i, :, :] = DiffResults.jacobian(out)
    end
    return moment_vals, moment_jac
end

function zero_in_hull_for_all(moments, weights, not_trimmed)
    n_obs = size(moments, 1)
    for i in 1:n_obs
        # For trimmed observations, we don't care as they will drop from the SER
        if not_trimmed[i]
            pos_weight = weights[i, :] .> 0
            if !zero_in_hull(moments[pos_weight, :])
                return false
            end
        end
    end
    return true
end

function zero_in_hull(moments)
    if size(moments, 2) == 1
        return any(moments .>= 0) && any(moments .<= 0)
    end
    # Based on https://stackoverflow.com/a/16898636/7388096
    hull = scipy_spatial.Delaunay(moments)
    return hull.find_simplex(zeros(size(moments, 2))) .>= 0
end

function λ_obj_fgh(f, grad, hess, γ, moment_vals, weights)
    n_obs = size(moment_vals, 1)
    if grad != nothing
        grad[:] .= 0
    end
    if hess != nothing
        hess[:, :] .= 0
    end
    out = 0.
    for j in 1:n_obs
        if weights[j] > 0
            inlog = 1 + transpose(γ) * moment_vals[j, :]
            if inlog <= 0
                if grad != nothing
                    grad[:] .= NaN
                end
                if hess != nothing
                    hess[:, :] .= NaN
                end
                if f != nothing
                    return Inf
                end
            end
            if f != nothing
                out -= weights[j] * log(inlog)
            end
            if grad != nothing
                grad[:] -= weights[j] * moment_vals[j, :] / inlog
            end
            if hess != nothing
                hess[:, :] -= weights[j] * moment_vals[j, :] * transpose(moment_vals[j, :]) / inlog
            end
        end
    end
    if f != nothing
        return out
    end
end

end


text:
function sodium_init()
    ccall((:sodium_init, libsodium), Cint, ())
end

function crypto_box_sealbytes()
    ccall((:crypto_box_sealbytes, libsodium), Csize_t, ())
end

function crypto_box_seal(c, m, mlen, pk)
    ccall((:crypto_box_seal, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}, Culonglong, Ptr{Cuchar}), c, m, mlen, pk)
end

ext: jl


text:
@testset "Poly: $TPoly" for (TPoly, T) in [(ArbPoly, Arb), (AcbPoly, Acb)]
    @testset "Length and degree" begin
        for (p, l) in [(zero(TPoly), 0), (one(TPoly), 1), (TPoly(T[0, 0, 1, 0]), 3)]
            @test firstindex(p) == 0
            @test size(p) == (l,)
            @test length(p) == size(p, 1) == Arblib.degree(p) + 1 == lastindex(p) + 1 == l
        end
    end

    @testset "Get and set coefficients" begin
        @test eltype(TPoly()) == T

        p = TPoly(T[i for i = 0:10])

        @test firstindex(p) == 0
        @test lastindex(p) == 10

        @test Arblib.coeffs(p) == p[:] == 0:10

        @test all(p[i] == i for i = 0:10)
        @test p[0:2] == 0:2

        @test p[11] == 0
        @test p[12] == 0

        p[3] = T(7)
        @test p[3] == 7
        p[4] = 8
        @test p[4] == 8
        p[5] = π
        @test isequal(p[5], T(π))
        if T == Arb
            p[6] = ArbRefVector(T[9])[1]
        else
            p[6] = AcbRefVector(T[9])[1]
        end
        @test p[6] == 9

        q = TPoly(prec = 512)
        q[0] = π
        @test Arblib.rel_accuracy_bits(q[0]) > 500

        @test_throws BoundsError p[-1]

        p = TPoly(T[i for i = 0:10])
        @test all(Arblib.ref(p, i) == p[i] for i = 0:10)

        Arblib.set!(Arblib.ref(p, 0), 5)
        @test p[0] == 5

        @test_throws BoundsError Arblib.ref(p, -1)
        @test_throws BoundsError Arblib.ref(p, 11)
        Arblib.zero!(Arblib.ref(p, 10))
        Arblib.normalise!(p)
        @test_throws BoundsError Arblib.ref(p, 10)
    end

    @testset "Constructors" begin
        @test TPoly() == TPoly(T(0)) == TPoly(T[0]) == zero(TPoly) == zero(TPoly())
        @test TPoly(T(1)) == TPoly(T[1]) == one(TPoly) == one(TPoly())
        @test Arblib.isx(TPoly(T[0, 1]))
        @test TPoly(T[1, 2, 0]) == TPoly(T[1, 2])
        @test TPoly(5) == TPoly(5.0) == TPoly([5.0]) == TPoly(T(5))
        @test TPoly(TPoly([1, 2])) == TPoly(ArbSeries([1, 2])) == TPoly([1, 2])

        @test precision(TPoly(prec = 64)) == 64
        @test precision(TPoly(T(0), prec = 64)) == 64
        @test precision(TPoly(T[0], prec = 64)) == 64
        @test precision(zero(TPoly(prec = 64))) == 64
        @test precision(one(TPoly(prec = 64))) == 64
        @test precision(TPoly(TPoly([1, 2], prec = 64))) == 64
        @test precision(TPoly(ArbSeries([1, 2], prec = 64))) == 64

        if TPoly == AcbPoly
            @test TPoly(AcbSeries([1, 2])) == TPoly(ArbPoly([1, 2])) == TPoly([1, 2])

            @test precision(TPoly(AcbSeries(prec = 64))) == 64
            @test precision(TPoly(ArbPoly(prec = 64))) == 64
        end

        if TPoly == ArbPoly
            @test Arblib.fromroots(TPoly, []) == one(TPoly)
            @test Arblib.fromroots(TPoly, [1, -1]) == TPoly([-1, 0, 1])
            @test Arblib.fromroots(TPoly, [], []) == one(TPoly)
            @test Arblib.fromroots(TPoly, [0], [im]) == TPoly([0, 1, 0, 1])
            @test Arblib.fromroots(TPoly, [], [im]) == TPoly([1, 0, 1])
            @test Arblib.fromroots(TPoly, [1, -1], []) == TPoly([-1, 0, 1])
        else
            @test Arblib.fromroots(TPoly, []) == one(TPoly)
            @test Arblib.fromroots(TPoly, [0, im]) == TPoly([0, -im, 1])
            @test Arblib.fromroots(TPoly, [0, im, -im]) == TPoly([0, 1, 0, 1])
        end
    end

    @testset "Arithmetic" begin
        p = TPoly([1, 2, 3])
        q = TPoly([2, 3, 4])

        @test p + q == TPoly([3, 5, 7])
        @test p - q == TPoly([-1, -1, -1])
        @test p * q == TPoly([2, 7, 16, 17, 12])

        @test -p == TPoly([-1, -2, -3])

        @test divrem(p, TPoly([0, 1])) == (TPoly([2, 3]), TPoly([1]))

        @test precision(setprecision(p, 80) + setprecision(q, 90)) == 90
        @test precision(setprecision(p, 80) - setprecision(q, 90)) == 90
        @test precision(setprecision(p, 80) * setprecision(q, 90)) == 90
        @test precision(-setprecision(p, 80)) == 80
        @test precision.(divrem(setprecision(p, 80), TPoly([0, 1], prec = 90))) == (90, 90)

        if TPoly == AcbPoly
            q = ArbPoly([2, 3, 4])
            @test p + q == q + p == TPoly([3, 5, 7])
            @test p - q == TPoly([-1, -1, -1])
            @test q - p == TPoly([1, 1, 1])
            @test p * q == q * p == TPoly([2, 7, 16, 17, 12])
        end
    end

    @testset "Scalar arithmetic" begin
        p = TPoly([1, 2, 3])

        @test p + T(2) ==
              T(2) + p ==
              p + 2 ==
              2 + p ==
              p + 2.0 ==
              2.0 + p ==
              TPoly([3, 2, 3])
        @test p - T(2) == p - 2 == p - 2.0 == TPoly([-1, 2, 3])
        @test T(2) - p == 2 - p == 2.0 - p == TPoly([1, -2, -3])
        @test p * T(2) ==
              T(2) * p ==
              p * 2 ==
              2 * p ==
              p * 2.0 ==
              2.0 * p ==
              TPoly([2, 4, 6])
        @test p / T(2) == p / 2 == p / 2.0 == TPoly([0.5, 1, 1.5])

        # Test with zero polynomial
        @test zero(TPoly) + 1 == TPoly(1)
        @test zero(TPoly) - 1 == TPoly(-1)
        @test 1 - zero(TPoly) == TPoly(1)
        @test 1 * zero(TPoly) == TPoly()
        @test zero(TPoly) / 1 == TPoly()

        # Test that the normalisation works
        @test iszero(TPoly(-1) + 1)
        @test iszero(TPoly(1) - 1)
        @test iszero(1 - TPoly(1))

        let p = setprecision(p, 80)
            @test precision(p + T(2)) ==
                  precision(T(2) + p) ==
                  precision(2 + p) ==
                  precision(p + 2) ==
                  80
            @test precision(p - T(2)) == precision(p - 2) == 80
            @test precision(T(2) - p) == precision(2 - p) == 80
            @test precision(p * T(2)) ==
                  precision(T(2) * p) ==
                  precision(2 * p) ==
                  precision(p * 2) ==
                  80
            @test precision(p / T(2)) == precision(p / 2) == 80
        end

        if TPoly == ArbPoly
            @test p + im == p + Acb(im) == AcbPoly([1 + im, 2, 3])
            @test p - im == p - Acb(im) == AcbPoly([1 - im, 2, 3])
            @test im - p == Acb(im) - p == AcbPoly([im - 1, -2, -3])
            @test p * im == p * Acb(im) == AcbPoly([im, 2im, 3im])
            @test p / im == p / Acb(im) == AcbPoly([-im, -2im, -3im])

            let p = setprecision(p, 80)
                @test precision(p + im) == 80
                @test precision(p - im) == 80
                @test precision(im - p) == 80
                @test precision(p * im) == 80
                @test precision(p / im) == 80
            end
        end

        if TPoly == AcbPoly
            @test p + Arb(2) == Arb(2) + p == TPoly([3, 2, 3])
            @test p - Arb(2) == TPoly([-1, 2, 3])
            @test Arb(2) - p == TPoly([1, -2, -3])
            @test p * Arb(2) == Arb(2) * p == TPoly([2, 4, 6])
            @test p / Arb(2) == TPoly([0.5, 1, 1.5])
        end
    end

    @testset "Composition" begin
        p = TPoly([1, 2])
        q = TPoly([2, 3])

        @test Arblib.taylor_shift(p, T(2)) == Arblib.taylor_shift(p, 2) == TPoly([5, 2])
        @test Arblib.compose(p, q) == TPoly([5, 6])

        @test precision(Arblib.taylor_shift(setprecision(p, 80), T(2))) == 80
        @test precision(Arblib.compose(setprecision(p, 80), setprecision(q, 90))) == 90
    end

    @testset "Evaluation" begin
        p = TPoly([1, 2, 3])

        @test p(Arb(2)) == p(2) == p(2.0) == 17
        @test p(Acb(1, 2)) == p(1 + 2im) == p(1.0 + 2.0im) == -6 + 16im

        @test Arblib.evaluate2(p, Arb(2)) == Arblib.evaluate2(p, 2) == (17, 14)
        @test Arblib.evaluate2(p, Acb(1, 2)) ==
              Arblib.evaluate2(p, 1 + 2im) ==
              (-6 + 16im, 8 + 12im)

        @test precision(p(T())) == precision(p(T(prec = 2precision(p)))) == precision(p)
        @test precision(TPoly(prec = 80)(T())) == 80
    end

    @testset "Differentiation and integration" begin
        p = TPoly([1, 2, 3])

        @test Arblib.derivative(p) == TPoly([2, 6])
        @test Arblib.integral(p) == TPoly([0, 1, 1, 1])

        @test precision(Arblib.derivative(p)) == precision(p)
        @test precision(Arblib.integral(p)) == precision(p)
        @test precision(Arblib.derivative(TPoly(prec = 80))) == 80
        @test precision(Arblib.integral(TPoly(prec = 80))) == 80
    end

    @testset "Power methods" begin
        p = TPoly([1, 2, 3])

        @test p^2 == p^Int(2) == p^UInt(2) == TPoly([1, 4, 10, 12, 9])

        @test precision(TPoly([1, 2, 3], prec = 80)^Int(2)) == 80

        @test_throws InexactError p^Int(-1)
    end
end


text:
using Test
include("ExtendableSparseTest.jl")

@test ExtendableSparseTest.test_constructors()
@test ExtendableSparseTest.test_timing()

@time begin

    @test ExtendableSparseTest.test_assembly(m=10,n=10,xnnz=5)
    @test ExtendableSparseTest.test_assembly(m=100,n=100,xnnz=500,nsplice=2)
    @test ExtendableSparseTest.test_assembly(m=1000,n=1000,xnnz=5000,nsplice=3)

    @test ExtendableSparseTest.test_assembly(m=20,n=10,xnnz=5)
    @test ExtendableSparseTest.test_assembly(m=200,n=100,xnnz=500,nsplice=2)
    @test ExtendableSparseTest.test_assembly(m=2000,n=1000,xnnz=5000,nsplice=3)

    @test ExtendableSparseTest.test_assembly(m=10,n=20,xnnz=5)
    @test ExtendableSparseTest.test_assembly(m=100,n=200,xnnz=500,nsplice=2)
    @test ExtendableSparseTest.test_assembly(m=1000,n=2000,xnnz=5000,nsplice=3)

    for irun=1:10
        m=rand((1:10000))
        n=rand((1:10000))
        nnz=rand((1:10000))
        nsplice=rand((1:5))
        @test ExtendableSparseTest.test_assembly(m=m,n=n,xnnz=nnz,nsplice=nsplice)
    end

    for irun=1:10
        m=rand((1:1000))
        n=rand((1:1000))
        d=0.3*rand()
        @test ExtendableSparseTest.test_transient_construction(m=m,n=n,d=d)
    end

    for irun=1:10
        m=rand((1:1000))
        n=rand((1:1000))
        d=0.3*rand()
        @test ExtendableSparseTest.test_addition(m=m,n=n,d=d)
    end

    @test ExtendableSparseTest.test_operations(10)
    @test ExtendableSparseTest.test_operations(100)
    @test ExtendableSparseTest.test_operations(1000)
end


text:
#from https://github.com/vtjnash/GLib.jl/blob/master/src/GLib.jl
# @vtjnash
#Cross-Referenced with: #checked
#https://github.com/JuliaGraphics/Gtk.jl/blob/3bd7caf1eabf0c295d1f092d44c7ac4fccb67a56/src/GLib/GLib.jl
#requires Libffi_jll
#=
✓ Preferences✓ JLLWrappers✓ Libffi_jll✗ FreeType2_jll✗ Libgcrypt_jll✗ Libtiff_jll✓ Preferences✓ JLLWrappers✓ Libffi_jll
✗ FreeType2_jll✗ Libgcrypt_jll✗ Libtiff_jll✗ XML2_jll✗ Fontconfig_jll✗ XSLT_jll✗ Gettext_jll✗ Glib_jll✗ ATK_jll✗ gdk_pixbuf_jll
✗ Cairo_jll✗ HarfBuzz_jll✗ Pango_jll✗ GTK3_jll✗ Librsvg_jll
=also requires GTK3_jll -> errors while building (Pkg.precompile() )
there're some errors:
in expression starting at C:\Users\adamus\.julia\packages\GTK3_jll\FbeBp\src\wrappers\x86_64-w64-mingw32.jl:4
in expression starting at C:\Users\adamus\.julia\packages\GTK3_jll\FbeBp\src\GTK3_jll.jl:2
this file vitally, heavily relies on Glib_jll
=#
module GLib
#=requires
1. ] add Glib_jll [REPL] {an obscure name for a set of various Julia Libraries} #Requires looking into into
2. file , MutableTypes, {gvalues, gerror, glist,gtype.jl} in the same directory
=#
if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

# Import `libgobject` and whatnot
using Glib_jll

if (false)
    function include(x)
        println("including $x")
        @time Base.include(x)
    end
end

import Base: convert, copy, show, size, length, getindex, setindex!, get,
    iterate, eltype, isempty, ndims, stride, strides, popfirst!,
    empty!, append!, reverse!, pushfirst!, pop!, push!, splice!, insert!,
    sigatomic_begin, sigatomic_end, Sys.WORD_SIZE, unsafe_convert, getproperty,
    getindex, setindex!

using Libdl  #--- define Libdl

export GInterface, GType, GObject, GBoxed, @Gtype, @Gabstract, @Giface
#export GEnum, GError, GValue, gvalue, make_gvalue, g_type
export GEnum, GError, GValue, gvalue, make_gvalue, @make_gvalue, g_type
export GList, GSList, glist_iter, _GSList, _GList, gobject_ref, gobject_move_ref
export signal_connect, signal_emit, signal_handler_disconnect
export signal_handler_block, signal_handler_unblock
export g_timeout_add, g_idle_add, @idle_add
# export setproperty!, getproperty WARNING MAJOR NAME CHANGE (files affected )
export set_gtk_property!, get_gtk_property
export GConnectFlags
export @sigatom, cfunction_

#---- 3 cfunctions (new)

cfunction_(@nospecialize(f), r, a::Tuple) = cfunction_(f, r, Tuple{a...})

@generated function cfunction_(f, R::Type{rt}, A::Type{at}) where {rt,at<:Tuple}
    quote
        @cfunction($(Expr(:$, :f)), $rt, ($(at.parameters...),))
    end
end

const gtk_eventloop_f = Ref{Function}()

# local function, handles Symbol and makes UTF8-strings easier
const AbstractStringLike = Union{AbstractString,Symbol}

bytestring(s) = String(s)
bytestring(s::Symbol) = s
bytestring(s::Ptr{UInt8}) = unsafe_string(s)
# bytestring(s::Ptr{UInt8}, own::Bool=false) = unsafe_string(s)

g_malloc(s::Integer) = ccall((:g_malloc, libglib), Ptr{Nothing}, (Csize_t,), s) # include(joinpath("..", "..", "deps", "ext_glib.jl")) #g_malloc(s::Integer) = ccall((:g_malloc, libglib), Ptr{Nothing}, (Csize_t,), s)
g_free(p::Ptr) = ccall((:g_free, libglib), Nothing, (Ptr{Nothing},), p) # ccall((:g_type_init, libgobject), Void, ())  #g_free(p::Ptr) = ccall((:g_free, libglib), Nothing, (Ptr{Nothing},), p)

include("MutableTypes.jl")
using .MutableTypes

include("gerror.jl")
include("glist.jl")
include("signals.jl")
include("gtype.jl")#ERROR

#include("gvalues.jl")#ERROR
#include("gwrap.jl") #Excluded with new releases

export @g_type_delegate
macro g_type_delegate(eq)
    @assert isa(eq, Expr) && eq.head == :(=) && length(eq.args) == 2
    new = eq.args[1]
    real = eq.args[2]
    newleaf = esc(Symbol(string(new, __module__.suffix)))
    realleaf = esc(Symbol(string(real, __module__.suffix)))
    new = esc(new)
    macroreal = QuoteNode(Symbol(string('@', real)))
    quote
        $newleaf = $realleaf
        macro $new(args...)
            Expr(:macrocall, $macroreal, map(esc, args)...)
        end
    end
end

if Base.VERSION >= v"1.4.2"
    precompile(Tuple{typeof(addref),Any})   # time: 0.003988418 #errors out
    # precompile(Tuple{typeof(gc_ref),Any})   # time: 0.002649791
end#ERROR #addref undefined

end # module

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
# ---
# title: 401. Binary Watch
# id: problem401
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Backtracking, Bit Manipulation
# link: <https://leetcode.com/problems/binary-watch/description/>
# hidden: true
# ---
#
# A binary watch has 4 LEDs on the top which represent the **hours** ( **0-11**
# ), and the 6 LEDs on the bottom represent the **minutes** ( **0-59** ).
#
# Each LED represents a zero or one, with the least significant bit on the
# right.
#
# ![](https://upload.wikimedia.org/wikipedia/commons/8/8b/Binary_clock_samui_moon.jpg)
#
# For example, the above binary watch reads "3:25".
#
# Given a non-negative integer _n_ which represents the number of LEDs that are
# currently on, return all possible times the watch could represent.
#
# **Example:**
#
#
#
#     Input: n = 1
#     Return: ["1:00", "2:00", "4:00", "8:00", "0:01", "0:02", "0:04", "0:08", "0:16", "0:32"]
#
# **Note:**
#
#   * The order of output does not matter.
#   * The hour must not contain a leading zero, for example "01:00" is not valid, it should be "1:00".
#   * The minute must be consist of two digits and may contain a leading zero, for example "10:2" is not valid, it should be "10:02".
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#DECISION put this somewhere else?  Its only an internal grouping
mutable struct RowInfo{LID <: Integer}
    graph::CSRGraph{<:Integer, <:Integer, LID}
    localRow::LID
    allocSize::LID
    numEntries::LID
    offset1D::LID
end

#RowInfo object for re-use
const rowInfoSpare = Union{RowInfo, Nothing}[nothing]

"""
    Gets a `RowInfo` object with the given values, reusing an intance if able
"""
@inline function createRowInfo(graph::CSRGraph{<:Integer, <:Integer, LID}, localRow::LID,
        allocSize::LID, numEntries::LID, offset1D::LID)::RowInfo{LID} where {LID <: Integer}
    global rowInfoSpare

    @inbounds begin#if length(rowInfoSpare) > 0
        #rowInfoSpare should always be size 1
        nextVal = rowInfoSpare[1]
        #ensure object pool haves right type
        if nextVal isa RowInfo{LID}
            rowInfoSpare[1] = nothing

            rowInfo::RowInfo{LID} = nextVal
            rowInfo.graph = graph
            rowInfo.localRow = localRow
            rowInfo.allocSize = allocSize
            rowInfo.numEntries = numEntries
            rowInfo.offset1D = offset1D
            return rowInfo
        end
    end
    #couldn't reuse, create new instance
    return RowInfo{LID}(graph, localRow, allocSize, numEntries, offset1D)
end

"""
    Puts the `RowInfo` object back in the object pool.
    After calling this method remove all references to the object.
"""
@inline function recycleRowInfo(rowInfo::RowInfo{T}) where T
    @inbounds rowInfoSpare[1] = rowInfo
    nothing
end

#TODO implement getLocalDiagOffsets(::CSRGraph)
getLocalGraph(graph::CSRGraph) = graph.localGraph

function updateGlobalAllocAndValues(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}, newAllocSize::Integer, rowValues::AbstractArray{Data, 1}) where {Data, GID, PID, LID}

    resize!(graph.globalIndices2D[rowInfo.localRow], newAllocSize)
    resize!(rowValues, newAllocSize)

    nothing
end

function insertIndicesAndValues(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}, newInds::Union{AbstractArray{GID, 1}, AbstractArray{LID, 1}}, oldRowVals::AbstractArray{Data, 1}, newRowVals::AbstractArray{Data, 1}, lg::IndexType) where {Data, GID, PID, LID}
    numNewInds = insertIndices(graph, rowInfo, newInds, lg)
    oldInd = rowInfo.numEntries+1

    oldRowVals[range(oldInd, step=1, length=numNewInds)] = newRowVals[1:numNewInds]
end

function insertIndices(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}, newInds::Union{AbstractArray{GID, 1}, AbstractArray{LID, 1}}, lg::IndexType) where {GID, PID, LID}
    numNewInds = LID(length(newInds))
    if lg == GLOBAL_INDICES
        if isGloballyIndexed(graph)
            numEntries = rowInfo.numEntries
            (gIndPtr, gIndLen) = getGlobalViewPtr(graph, rowInfo)
            @assert gIndLen >= numNewInds+numEntries
            for i in 1:numNewInds
                unsafe_store!(gIndPtr, GID(newInds[i]), numEntries+i)
            end
        else
            lIndView = getLocalView(graph, rowInfo)
            colMap = graph.colMap

            dest = range(rowInfo.numEntries, step=1, length=numNewInds)
            lIndView[dest] = [lid(colMap, GID(id)) for id in newInds]
        end
    elseif lg == LOCAL_INDICES
        if isLocallyIndexed(graph)
            numEntries = rowInfo.numEntries
            (lIndPtr, lIndLen) = getLocalViewPtr(graph, rowInfo)
            @assert gIndLen >= numNewInds+numEntries
            for i in 1:numNewInds
                unsafe_store!(lIndPtr, LID(newInds[i]), numEntries+i)
            end
        else
            @assert(false,"lg=LOCAL_INDICES, isGloballyIndexed(g) not implemented, "
            * "because it doesn't make sense")
        end
    end

    graph.numRowEntries[rowInfo.localRow] += numNewInds
    graph.nodeNumEntries += numNewInds
    setLocallyModified(graph)

    numNewInds
end

function computeGlobalConstants(graph::CSRGraph{GID, PID, LID}) where {
        GID <: Integer, PID <: Integer, LID <: Integer}

    #short circuit if already computed
    graph.haveGlobalConstants && return

    computeLocalConstants(graph)

    commObj = getComm(getRowMap(graph))

    #if graph.haveGlobalConstants == false  #short circuited above
    graph.globalNumEntries, graph.globalNumDiags = sumAll(commObj,
        [GID(graph.nodeNumEntries), GID(graph.nodeNumDiags)])

    graph.globalMaxNumRowEntries = maxAll(commObj, GID(graph.nodeMaxNumRowEntries))
    graph.haveGlobalConstants = true
end

function clearGlobalConstants(graph::CSRGraph)
    graph.globalNumEntries = 0
    graph.globalNumDiags = 0
    graph.globalMaxNumRowEntries = 0
    graph.haveGlobalConstants = false
end

function computeLocalConstants(graph::CSRGraph{GID, PID, LID}) where {
        GID <: Integer, PID <: Integer, LID <: Integer}

    #short circuit if already computed
    graph.haveLocalConstants && return

    #if graph.haveLocalConstants == false  #short circuited above

    graph.upperTriangle = true
    graph.lowerTriangle = true
    graph.nodeMaxNumRowEntries = 0
    graph.nodeNumDiags = 0

    rowMap = graph.rowMap
    colMap = graph.colMap

    #indicesAreAllocated => true
    if  hasRowInfo(graph)
        numLocalRows = numMyElements(rowMap)
        for localRow = LID(1):numLocalRows
            globalRow = gid(rowMap, localRow)
            rowLID = lid(colMap, globalRow)

            rowInfo = getRowInfo(graph, localRow)
            (rowPtr, rowLen) = getLocalViewPtr(graph, rowInfo)

            for i in 1:rowLen
                if unsafe_load(rowPtr) == rowLID
                    graph.nodeNumDiags += 1
                    break
                end
            end

            smallestCol::LID = unsafe_load(rowPtr, 1)
            largestCol::LID  = unsafe_load(rowPtr, rowLen)
            if smallestCol < localRow
                graph.upperTriangle = false
            end
            if localRow < largestCol
                graph.lowerTriangle = false
            end
            graph.nodeMaxNumRowEntries = max(graph.nodeMaxNumRowEntries, rowInfo.numEntries)

            recycleRowInfo(rowInfo)
        end
    end
    graph.haveLocalConstants = true
end

hasRowInfo(graph::CSRGraph) = (getProfileType(graph) != STATIC_PROFILE
                                || length(graph.rowOffsets) != 0)

Base.@propagate_inbounds function getRowInfoFromGlobalRow(graph::CSRGraph{GID, PID, LID},
        row::Integer)::RowInfo{LID} where {GID, PID, LID <: Integer}
    getRowInfo(graph, lid(graph.rowMap, row))
end

@inline function getRowInfo(graph::CSRGraph{GID, PID, LID}, row::LID)::RowInfo{LID} where {GID, PID, LID <: Integer}

    emptyRowInfo = !hasRowInfo(graph)
    @boundscheck emptyRowInfo = emptyRowInfo || !myLID(graph.rowMap, row)
    if emptyRowInfo
        return createRowInfo(graph, row, LID(0), LID(0), LID(1))
    end

    offset1D::LID = 1
    allocSize::LID = 0

    @inbounds if getProfileType(graph) == STATIC_PROFILE
        if length(graph.rowOffsets) != 0
            offset1D  = LID(graph.rowOffsets[row])
            allocSize = LID(graph.rowOffsets[row+1] - graph.rowOffsets[row])
        end
        numEntries = (length(graph.numRowEntries) == 0 ?
            allocSize : LID(graph.numRowEntries[row]))
    else #dynamic profile
        if isLocallyIndexed(graph) && length(graph.localIndices2D) != 0
            allocSize = LID(length(graph.localIndices2D[row]))

        elseif isGloballyIndexed(graph) && length(graph.globalIndices2D) != 0
            allocSize = LID(length(graph.globalIndices2D[row]))
        end
        numEntries = (length(graph.numRowEntries) == 0 ?
            LID(0) : LID(graph.numRowEntries[row]))
    end
    createRowInfo(graph, row, allocSize, numEntries, offset1D)
end

function getLocalView(rowInfo::RowInfo{LID})::AbstractArray{LID, 1} where LID <: Integer
    graph = rowInfo.graph
    if rowInfo.allocSize == 0
        LID[]
    elseif length(graph.localIndices1D) != 0
        start = rowInfo.offset1D
        len = rowInfo.allocSize

        view(graph.localIndices1D, range(start, step=1, length=len))
    elseif length(graph.localIndices2D[rowInfo.localRow]) != 0
        graph.localIndices2D[rowInfo.localRow]
    else
        LID[]
    end
end

function allocateIndices(graph::CSRGraph{GID, <:Integer, LID},
        lg::IndexType, numAllocPerRow::AbstractArray{<:Integer, 1}) where {
        GID <: Integer, LID <: Integer}
    numRows = getLocalNumRows(graph)
    @assert(length(numAllocPerRow) == numRows,
        "numAllocRows has length = $(length(numAllocPerRow)) "
        * "!= numRows = $numRows")
    allocateIndices(graph, lg, numAllocPerRow, i -> numAllocPerRow[i])
end

function allocateIndices(graph::CSRGraph{GID, <:Integer, LID},
        lg::IndexType, numAllocPerRow::Integer) where {
        GID <: Integer, LID <: Integer}

# Manually inlined since the function version is too expensive

    @assert(isLocallyIndexed(graph) == (lg == LOCAL_INDICES),
        "Graph is $(isLocallyIndexed(graph) ? "" : "not ")locally indexed, but lg=$lg")
    @assert(isGloballyIndexed(graph) == (lg == GLOBAL_INDICES),
        "Graph is $(isGloballyIndexed(graph) ? "" : "not ")globally indexed but lg=$lg")

    numRows = getLocalNumRows(graph)

    if getProfileType(graph) == STATIC_PROFILE
        rowPtrs = Array{LID, 1}(undef, numRows + 1)

        computeOffsets(rowPtrs, numAllocPerRow)

        graph.rowOffsets = rowPtrs
        numInds = rowPtrs[numRows+1]

        if lg == LOCAL_INDICES
            graph.localIndices1D = Array{LID, 1}(undef, numInds)
        else
            graph.globalIndices1D = Array{GID, 1}(undef, numInds)
        end
        graph.storageStatus = STORAGE_1D_UNPACKED
    else
        if lg == LOCAL_INDICES
            graph.localIndices2D = Array{Array{LID, 1}, 1}(undef, numRows)
            for row = 1:numRows
                graph.localIndices2D[row] = Array{LID, 1}(undef, numAllocPerRow)
            end
        else #lg == GLOBAL_INDICES
            graph.globalIndices2D = Array{Array{GID, 1}, 1}(undef, numRows)
            for row = 1:numRows
                graph.globalIndices2D[row] = Array{GID, 1}(undef, numAllocPerRow)
            end
        end
        graph.storageStatus = STORAGE_2D
    end

    graph.indicesType = lg

    if numRows > 0
        numRowEntries = zeros(LID, numRows)
        graph.numRowEntries = numRowEntries
    end
end

function allocateIndices(graph::CSRGraph{GID, <:Integer, LID},
        lg::IndexType, numAlloc, numAllocPerRow::Function) where {
        GID <: Integer, LID <: Integer}

    @assert(isLocallyIndexed(graph) == (lg == LOCAL_INDICES),
        "Graph is $(isLocallyIndexed(graph) ? "" : "not ")locally indexed, but lg=$lg")
    @assert(isGloballyIndexed(graph) == (lg == GLOBAL_INDICES),
        "Graph is $(isGloballyIndexed(graph) ? "" : "not ")globally indexed but lg=$lg")

    numRows = getLocalNumRows(graph)

    if getProfileType(graph) == STATIC_PROFILE
        rowPtrs = Array{LID, 1}(undef, numRows + 1)

        computeOffsets(rowPtrs, numAlloc)

        graph.rowOffsets = rowPtrs
        numInds = rowPtrs[numRows+1]

        if lg == LOCAL_INDICES
            graph.localIndices1D = Array{LID, 1}(undef, numInds)
        else
            graph.globalIndices1D = Array{GID, 1}(undef, numInds)
        end
        graph.storageStatus = STORAGE_1D_UNPACKED
    else
        if lg == LOCAL_INDICES
            graph.localIndices2D = Array{Array{LID, 1}, 1}(undef, numRows)
            for row = 1:numRows
                graph.localIndices2D[row] = Array{LID, 1}(undef, numAllocPerRow(row))
            end
        else #lg == GLOBAL_INDICES
            graph.globalIndices2D = Array{Array{GID, 1}, 1}(undef, numRows)
            for row = 1:numRows
                graph.globalIndices2D[row] = Array{GID, 1}(undef, numAllocPerRow(row))
            end
        end
        graph.storageStatus = STORAGE_2D
    end

    graph.indicesType = lg

    if numRows > 0
        numRowEntries = zeros(LID, numRows)
        graph.numRowEntries = numRowEntries
    end
end

function makeImportExport(graph::CSRGraph{GID, PID, LID}) where {
        GID <: Integer, PID <: Integer, LID <: Integer}
    @assert (graph.colMap != nothing) "Cannot make imports and exports without a column map"

    if graph.importer == nothing
        if !sameAs(graph.domainMap, graph.colMap)
            graph.importer = Import(graph.domainMap, graph.colMap)
        end
    end

    if graph.exporter == nothing
        if !sameAs(graph.rangeMap, graph.rowMap)
            graph.exporter = Export(graph.rowMap, graph.rangeMap)
        end
    end
end

"""
    checkInternalState(::CSRGraph)

Checks the graphs internal state for correctness.
Mainly useful for intnernal debugging and testing purposes
"""
function checkInternalState(graph::CSRGraph)
    localNumRows = getLocalNumRows(graph)

    @assert(isFillActive(graph) != isFillComplete(graph),
        "Graph must be either fill active or fill "
        * "complete$(isFillActive(graph) ? "not both" : "").")
    @assert(!isFillComplete(graph)
            || (graph.colMap != nothing
                && graph.domainMap != nothing
                && graph.rangeMap != nothing),
        "Graph is fill complete, but at least one of {column, range, domain} map is null")
    @assert((graph.storageStatus != STORAGE_1D_PACKED
                && graph.storageStatus != STORAGE_1D_UNPACKED)
        || graph.pftype != DYNAMIC_PROFILE,
        "Graph claims 1D storage, but dynamic profile")
    if graph.storageStatus == STORAGE_2D
        @assert(graph.pftype != STATIC_PROFILE ,
            "Graph claims 2D storage, but static profile")
        @assert(!isLocallyIndexed(graph)
            || length(graph.localIndices2D) == localNumRows,
            "Graph calims to be locally index and have 2D storage, "
            * "but length(graph.localIndices2D) = $(length(graph.localIndices2D)) "
            * "!= getLocalNumRows(graph) = $localNumRows")
        @assert(!isGloballyIndexed(graph)
            || length(graph.globalIndices2D) == localNumRows,
            "Graph calims to be globally index and have 2D storage, "
            * "but length(graph.globalIndices2D) = $(length(graph.globalIndices2D)) "
            * "!= getLocalNumRows(graph) = $localNumRows")
    end

    @assert(graph.haveGlobalConstants
        || (graph.globalNumEntries == 0
            && graph.globalNumDiags == 0
            && graph.globalMaxNumRowEntries == 0),
        "Graph claims to not have global constants, "
        * "but some of the global constants are not 0")

    @assert(!graph.haveGlobalConstants
        || (graph.globalNumEntries != 0
            && graph.globalMaxNumRowEntries != 0),
        "Graph claims to have global constants, but also says 0 global entries")

    @assert(!graph.haveGlobalConstants
        || (graph.globalNumEntries >= graph.nodeNumEntries
            && graph.globalNumDiags >= graph.nodeNumDiags
            && graph.globalMaxNumRowEntries >= graph.nodeMaxNumRowEntries),
        "Graph claims to have global constants, but some of the local "
        * "constants are greater than their corresponding global constants")

    @assert(!isStorageOptimized(graph)
        || graph.pftype == STATIC_PROFILE,
        "Storage is optimized, but graph is not STATIC_PROFILE")

    @assert(!isGloballyIndexed(graph)
        || length(graph.rowOffsets) == 0
        || (length(graph.rowOffsets) == localNumRows +1
            && graph.rowOffsets[localNumRows+1] == length(graph.globalIndices1D)),
        "If rowOffsets has nonzero size and the graph is globally "
        * "indexed, then rowOffsets must have N+1 rows and rowOffsets[N+1] "
        * "must equal the length of globalIndices1D")

    @assert(!isLocallyIndexed(graph)
        || length(graph.rowOffsets) == 0
        || (length(graph.rowOffsets) == localNumRows +1
            && graph.rowOffsets[localNumRows+1] == length(graph.localIndices1D)),
        "If rowOffsets has nonzero size and the graph is globally "
        * "indexed, then rowOffsets must have N+1 rows and rowOffsets[N+1] "
        * "must equal the length of localIndices1D")

    if graph.pftype == DYNAMIC_PROFILE
        @assert(localNumRows == 0
            || length(graph.localIndices2D) > 0
            || length(graph.globalIndices2D) > 0,
            "Graph has dynamic profile, the calling process has nonzero "
            * "rows, but no 2-D column index storage is present.")
        @assert(localNumRows == 0
            || length(graph.numRowEntries) != 0,
            "Graph has dynamic profiles and the calling process has "
            * "nonzero rows, but numRowEntries is not present")

        @assert(length(graph.localIndices1D) == 0
            && length(graph.globalIndices1D) == 0,
            "Graph has dynamic profile, but 1D allocations are present")

        @assert(length(graph.rowOffsets) == 0,
            "Graph has dynamic profile, but row offsets are present")

    elseif graph.pftype == STATIC_PROFILE
        @assert(length(graph.localIndices1D) != 0
            || length(graph.globalIndices1D) != 0,
            "Graph has static profile, but 1D allocations are not present")

        @assert(length(graph.localIndices2D) == 0
            && length(graph.globalIndices2D) == 0,
            "Graph has static profile, but 2D allocations are present")
    else
        error("Unknown profile type: $(graph.pftype)")
    end

    if graph.indicesType == LOCAL_INDICES
        @assert(length(graph.globalIndices1D) == 0
            && length(graph.globalIndices2D) == 0,
            "Indices are local, but global allocations are present")
        @assert(graph.nodeNumEntries == 0
            || length(graph.localIndices1D) > 0
            || length(graph.localIndices2D) > 0,
            "Indices are local and local entries exist, but there aren't local allocations present")
    elseif graph.indicesType == GLOBAL_INDICES
        @assert(length(graph.localIndices1D) == 0
            && length(graph.localIndices2D) == 0,
            "Indices are global, but local allocations are present")
        @assert(graph.nodeNumEntries == 0
            || length(graph.globalIndices1D) > 0
            || length(graph.globalIndices2D) > 0,
            "Indices are global and local entries exist, but there aren't global allocations present")
    else
        warn("Unknown indices type: $(graph.indicesType)")
    end

    #check actual allocations
    lenRowOffsets = length(graph.rowOffsets)
    if graph.pftype == STATIC_PROFILE && lenRowOffsets != 0
        @assert(lenRowOffsets == localNumRows+1,
            "Graph has static profile, rowOffsets has a nonzero length "
            * "($lenRowOffsets), but is not equal to the "
            * "local number of rows plus one ($(localNumRows+1))")
        actualNumAllocated = graph.rowOffsets[localNumRows+1]
        @assert(!isLocallyIndexed(graph)
            || length(graph.localIndices1D) == actualNumAllocated,
            "Graph has static profile, rowOffsets has a nonzero length, "
            * "but length(localIndices1D) = $(length(graph.localIndices1D)) "
            * "!= actualNumAllocated = $actualNumAllocated")
        @assert(!isGloballyIndexed(graph)
            || length(graph.globalIndices1D) == actualNumAllocated,
            "Graph has static profile, rowOffsets has a nonzero length, "
            * "but length(globalIndices1D) = $(length(graph.globalIndices1D)) "
            * "!= actualNumAllocated = $actualNumAllocated")
    end
    true #if an problem occurred, an error was already thrown
end

function setLocallyModified(graph::CSRGraph)
    graph.indicesAreSorted = false
    graph.noRedundancies = false
    graph.haveLocalConstants = false
end

function sortAndMergeAllIndices(graph::CSRGraph, sorted::Bool, merged::Bool)
    @assert(isLocallyIndexed(graph),
        "This method may only be called after makeIndicesLocal(graph)")
    @assert(merged || isStoragedOptimized(graph),
        "The graph is already storage optimized, "
        * "so we shouldn't be merging any indices.")

    if !sorted || !merged
        localNumRows = getLocalNumRows(graph)
        totalNumDups = 0
        for localRow = 1:localNumRows
            rowInfo = getRowInfo(graph, localRow)
            if !sorted
                sortRowIndices(graph, rowInfo)
            end
            if !merged
                numDups += mergeRowIndices(graph, rowInfo)
            end
            recycleRowInfo(rowInfo)
        end
        graph.nodeNumEntries -= totalNumDups
        graph.indiciesAreSorted = true
        graph.noRedunancies = true
    end
end

function sortRowIndices(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}) where {GID, PID, LID <: Integer}
    if rowInfo.numEntries > 0
        localColumnIndices = getLocalView(graph, rowInfo)
        sort!(localColumnIndices)
    end
end

function mergeRowIndices(graph::CSRGraph{GID, PID, LID}, rowInfo::RowInfo{LID}) where {GID, PID, LID <: Integer}
    localColIndices = getLocalView(graph, rowInfo)
    localColIndices[:] = unique(localColIndices)
    mergedEntries = length(localColIndices)
    graph.numRowEntries[rowInfo.localRow] = mergedEntries

    rowInfo.numEntries - mergedEntries
end

function setDomainRangeMaps(graph::CSRGraph{GID, PID, LID}, domainMap::BlockMap{GID, PID, LID}, rangeMap::BlockMap{GID, PID, LID}) where {GID, PID, LID}
    if graph.domainMap != domainMap
        graph.domainMap = domainMap
        graph.importer = nothing
    end
    if graph.rangeMap != rangeMap
        graph.rangeMap = rangeMap
        graph.exporter = nothing
    end
end

function globalAssemble(graph::CSRGraph)
    @assert isFillActive(graph) "Fill must be active before calling globalAssemble(graph)"

    comm = getComm(graph)
    myNumNonlocalRows = length(graph.nonlocals)

    maxNonlocalRows = maxAll(comm, myNumNonlocalRows)
    if maxNonlocalRows != 0
        return
    end

    #skipping: nonlocalRowMap = null

    numEntPerNonlocalRow = Array{LID, 1}(undef, myNumNonlocalRows)
    myNonlocalGlobalRows = Array{GID, 1}(undef, myNumNonlocalRows)

    for (i, (key, val)) = zip(1:length(graph.nonlocals), graph.nonlocals)
        myNonlocalGlobalRows[i] = key
        globalCols = val
        sort!(globalCols)
        globalCols[:] = unique(globalCols)
        numEntPerNonlocalRow[i] = length(globalCols)
    end

    myMinNonLocalGlobalRow = minimum(myNonLocalGlobalRows)

    globalMinNonlocalRow = minAll(comm, myMinNonlocalGlobalRow)

    nonlocalRowMap = BlockMap(-1, myNonlocalGlobalRows, comm)

    nonlocalGraph = CSRGraph(nonlocalRowMap, numEntPerNonlocalRow, STATIC_PROFILE)
    for (i, (key, val)) = zip(1:length(graph.nonlocals), graph.nonlocals)
        globalRow = key
        globalColumns = val
        numEnt = length(numEntPerNonlocalRow[i])
        insertGlobalIndices(nonLocalGraph, globalRow, numEnt, globalColumns)
    end

    origRowMap = graph.rowMap
    origRowMapIsOneToOne = isOneToOne(origRowMap)

    if origRowMapIsOneToOne
        exportToOrig = Export(nonlocalRowMap, origRowMap)
        doExport(nonLocalGraph, graph, exportToOrig, INSERT)
    else
        oneToOneRowMap = createOneToOne(origRowMap)
        exportToOneToOne = Export(nonlocalRowMap, oneToOneRowMap)

        oneToOneGraph = CSRGraph(oneToOneRowMap, 0)
        doExport(nonlocalGraph, oneToOneGraph, exportToOneToOne, INSERT)

        #keep memory highwater mark down
        #nonlocalGraph = null

        importToOrig(oneToOneRowMap, origRowMap)
        doImport(oneToOneGraph, graph, importToOrig, INSERT)
    end
    clear!(graph.nonLocals)
end

function makeIndicesLocal(graph::CSRGraph{GID, PID, LID}) where {GID, PID, LID}
    @assert hasColMap(graph) "The graph does not have a column map yet.  This method should never be called in that case"

    colMap = graph.colMap
    localNumRows = getLocalNumRows(graph)

    if isGloballyIndexed(graph) && localNumRows != 0
        numRowEntries = graph.numRowEntries

        if getProfileType(graph) == STATIC_PROFILE
            if GID == LID
                graph.localIndices1D = graph.globalIndices1D

            else
                @assert(length(graph.rowOffsets) != 0,
                    "length(graph.rowOffsets) == 0.  "
                    * "This should never happen at this point")
                numEnt = graph.rowOffsets[localNumRows+1]
                graph.localIndices1D = Array{LID, 1}(undef, numEnt)
            end

            localColumnMap = getLocalMap(colMap)

            numBad = convertColumnIndicesFromGlobalToLocal(
                        graph.localIndices1D,
                        graph.globalIndices1D,
                        graph.rowOffsets,
                        localColumnMap,
                        numRowEntries)

            if numBad != 0
                throw(InvalidArgumentError("When converting column indices from "
                        * "global to local, we encountered $numBad indices that "
                        * "do not live in the column map on this process"))
            end

            graph.globalIndices1D = Array{LID, 1}(undef, 0)
        else #graph has dynamic profile
            graph.localIndices2D = Array{Array{LID, 1}, 1}(undef, localNumRows)
            for localRow = 1:localNumRows
                if length(graph.globalIndices2D[localRow]) != 0
                    globalIndices = graph.globalIndices2D[localRow]

                    graph.localIndices2D[localRow] = [lid(colMap, gid) for gid in globalIndices]
                end
            end
            graph.globalIndices2D = Array{GID, 1}[]
        end
    end

    graph.localGraph = LocalCSRGraph(graph.localIndices1D, graph.rowOffsets)
    graph.indicesType = LOCAL_INDICES
end

function convertColumnIndicesFromGlobalToLocal(localColumnIndices::AbstractArray{LID, 1},
        globalColumnIndices::AbstractArray{GID, 1}, ptr::AbstractArray{LID, 1},
        localColumnMap::BlockMap{GID, PID, LID}, numRowEntries::AbstractArray{LID, 1}
        )::LID where {GID, PID, LID}

    localNumRows = max(length(ptr)-1, 0)
    numBad = 0
    for localRow = 1:localNumRows
        offset = ptr[localRow]

        for j = 0:numRowEntries[localRow]-1
            gid = globalColumnIndices[offset+j]
            localColumnIndices[offset+j] = lid(localColumnMap, gid)
            if localColumnIndices[offset+j] == 0
                numBad += 1
            end
        end
    end
    numBad
end

#covers the overlap between insert methods
macro insertIndicesImpl(indicesType, innards)
    indices1D = Symbol(indicesType*"Indices1D")
    indices2D = Symbol(indicesType*"Indices2D")

    esc(quote
            rowInfo = getRowInfo(graph, myRow)
            numNewIndices = length(indices)
            newNumEntries = rowInfo.numEntries + numNewIndices

            if newNumEntries > rowInfo.allocSize
                if getProfileType(graph) == STATIC_PROFILE
                    $innards
                else
                    newAllocSize = 2*rowInfo.allocSize
                    if newAllocSize < newNumEntries
                        newAllocSize = newNumEntries
                    end
                    resize!(graph.$indices2D[myRow], newAllocSize)
                end
            end

            if length(graph.$indices1D) != 0
                offset = rowInfo.offset1D + rowInfo.numEntries
                destRange = offset+1:offset+numNewIndices

                graph.$indices1D[destRange] = indices[:]
            else
                graph.$indices2D[myRow][rowInfo.numEntries+1:newNumEntries] = indices[:]
            end

            graph.numRowEntries[myRow] += numNewIndices
            graph.nodeNumEntries += numNewIndices
            setLocallyModified(graph)

            recycleRowInfo(rowInfo)
            nothing
    end)
end

function insertLocalIndicesImpl(graph::CSRGraph{GID, PID, LID},
        myRow::LID, indices::AbstractArray{LID, 1}) where {
        GID, PID, LID <: Integer}
    @insertIndicesImpl "local" begin
        throw(InvalidArgumentError("new indices exceed statically allocated graph structure"))
    end
end

#TODO figure out if this all can be moved to @insertIndicesImpl
function insertGlobalIndicesImpl(graph::CSRGraph{GID, PID, LID},
        myRow::LID, indices::AbstractArray{GID, 1}) where {
        GID <: Integer, PID, LID <: Integer}
    @insertIndicesImpl "global" begin
        @assert(rowInfo.numEntries <= rowInfo.allocSize,
            "For local row $myRow, rowInfo.numEntries = $(rowInfo.numEntries) "
            * "> rowInfo.allocSize = $(rowInfo.allocSize).")

        dupCount = 0
        if length(graph.globalIndices1D) != 0
            curOffset = rowInfo.offset1D
            @assert(length(graph.globalIndices1D) >= curOffset,
                "length(graph.globalIndices1D) = $(length(graph.globalIndices1D)) "
                * ">= curOffset = $curOffset")
            @assert(length(graph.globalIndices1D) >= curOffset + rowInfo.offset1D,
                "length(graph.globalIndices1D) = $(length(graph.globalIndices1D)) "
                * ">= curOffset+rowInfo.offset1D = $(curOffset + rowInfo.offset1D)")

            range = curOffset:curOffset+rowInfo.numEntries
            globalIndicesCur = view(graph.globalIndices1D, range)
        else
            #line 1959

            globalIndices = graph.globalIndices2D[myRow]
            @assert(rowInfo.allocSize == length(globalIndices),
                "rowInfo.allocSize = $(rowInfo.allocSize) "
                * "== length(globalIndices) = $(length(globalIndices))")
            @assert(rowInfo.numEntries <= length(globalIndices),
                "rowInfo.numEntries = $(rowInfo.numEntries) "
                * "== length(globalIndices) = $(length(globalIndices))")

            globalIndicesCur = view(globalIndices, 0, rowInfo.numEntries)
        end
        for newIndex = indices
            dupCount += count(old -> old==newIndex, globalIndicesCur)
        end

        numNewToInsert = numNewInds - dupCount
        @assert numNewToInsert >= 0 "More duplications than indices"

        if rowInfo.numEntries + numNewToInsert > rowInfo.allocSize
            throw(InvalidArgumentError("$(myPid(getComm(graph))): "
                    * "For local row $myRow, even after excluding "
                    * "$dupCount duplicate(s) in input, the new number "
                    * "of entries $(rowInfo.numEntries + numNewToInsert) "
                    * "still exceeds this row's static allocation size "
                    * "$(rowInfo.allocSize).  You must either fix the upper "
                    * "bound on number of entries in this row, or switch "
                    * "to dynamic profile."))
        end

        if length(graph.globalIndices) != 0
            curOffset = rowInfo.offset1D
            globalIndicesCur = view(graph.globalIndices1D,
                range(curOffset, step=1, length=rowInfo.numEntries))
            globalIndicesNew = view(graph.globalIndices1D,
                curOffset+rowInfo.numEntries+1 : currOffset+rowInfo.allocSize)
        else
            #line 2036

            globalIndices = graph.globalIndices2D[myRow]
            globalIndicesCur = view(globalIndices, 1:rowInfo.numEntries)
            globalIndicesNew = view(globalIndices,
                rowInfo.numEntries+1 : rowInfo.allocSize-rowInfo.numEntries)
        end

        curPos = 1
        for globalIndexToInsert = indices

            alreadyInOld = globalIndexToInsert in globalIndicesCur
            if !alreadyInOld
                @assert(curPos <= numNewToInsert,
                    "curPos = $curPos >= numNewToInsert = $newToInsert.")
                globalIndicesNew[curPos] = globalIndexToInsert
                curPos += 1
            end
        end

        graph.numRowEntries[myRow] = rowInfo.numEntries+numNewToInsert
        graph.nodeNumEntries += numNewToInsert
        setLocallyModified(graph)
        return
    end
end

#internal implementation of makeColMap, needed to handle some return and debuging stuff
#returns Tuple(error, colMap)
function __makeColMap(graph::CSRGraph{GID, PID, LID}, domMap::Union{BlockMap{GID, PID, LID}, Nothing}
        ) where {GID, PID, LID}
    error = false

    if domMap == nothing
        return nothing
    end

    if isLocallyIndexed(graph)
        colMap = graph.colMap

        if colMap == nothing
            warn("$(myPid(getComm(graph))): The graph is locally indexed, but does not have a column map")

            error = true
            myColumns = GID[]
        else
            if linearMap(colMap) #i think isContiguous(map) <=> linearMap(map)?
                numCurGIDs = numMyElements(colMap)
                myFirstGlobalIndex = minMyGIDs(colMap)

                myColumns = collect(range(myFirstGlobalIndex, step=1, length=numCurGIDs))
            else
                myColumns = copy(myGlobalElements(colMap))
            end
        end
        return (error, BlockMap(numGlobalElements(colMap), myColumns, getComm(domMap)))
    end

    #else if graph.isGloballyIndexed
    localNumRows = getLocalNumRows(graph)

    numLocalColGIDs = 0

    gidIsLocal = falses(localNumRows)
    remoteGIDSet = Set()

    #if rowMap != null
    rowMap = graph.rowMap

    for localRow = 1:localNumRows
        globalRow = gid(rowMap, localRow)
        (rowGIDPtr, numEnt) = getGlobalRowViewPtr(graph, globalRow)

        if numEnt != 0
            for k = 1:numEnt
                gid::GID = unsafe_load(rowGIDPtr, k)
                lid::LID = JuliaPetra.lid(domMap, gid)
                if lid != 0
                    @inbounds if !gidIsLocal[lid]
                        gidIsLocal[lid] = true
                        numLocalColGIDs += 1
                    end
                else
                    #don't need containment checks, set already takes care of that
                    push!(remoteGIDSet, gid)
                end
            end
        end
    end

    numRemoteColGIDs = length(remoteGIDSet)

    #line 214, abunch of explanation of serial short circuit
    if numProc(getComm(domMap)) == 1
        if numRemoteColGIDs != 0
            error = true
        end
        if numLocalColGIDs == localNumRows
            return (error, domMap)
        end
    end
    myColumns = Vector{GID}(undef, numLocalColGIDs+numRemoteColGIDs)
    localColGIDs  = view(myColumns, 1:numLocalColGIDs)
    remoteColGIDs = view(myColumns, numLocalColGIDs+1:numLocalColGIDs+numRemoteColGIDs)

    remoteColGIDs[:] = [el for el in remoteGIDSet]

    remotePIDs = Array{PID, 1}(undef, numRemoteColGIDs)

    remotePIDs = remoteIDList(domMap, remoteColGIDs)[1]
    if any(remotePIDs .== 0)
        warn("Some column indices are not in the domain Map")
        error = true
    end

    order = sortperm(remotePIDs)
    permute!(remotePIDs, order)
    permute!(remoteColGIDs, order)

    #line 333

    numDomainElts = numMyElements(domMap)
    if numLocalColGIDs == numDomainElts
        if linearMap(domMap) #I think isContiguous() <=> linearMap()
            localColGIDs[1:numLocalColGIDs] = range(minMyGID(domMap), step=1, length=numLocalColGIDs)
        else
            domElts = myGlobalElements(domMap)
            localColGIDs[1:length(domElts)] = domElts
        end
    else
        numLocalCount = 0
        if linearMap(domMap) #I think isContiguous() <=> linearMap()
            curColMapGID = minMyGID(domMap)
            for i = 1:numDomainElts
                if gidIsLocal[i]
                    numLocalCount += 1
                    localColGIDs[numLocalCount] = curColMapGID
                end
                curColMapGID += 1
            end
        else
            domainElts = myGlobalElement(domMap)
            for i = 1:numDomainElts
                if gidIsLocal[i]
                    numLocalCount += 1
                    localColGIDs[numLocalCount] = domainElts[i]
                end
                curColMapGID += 1
            end
        end

        if numLocalCount != numLocalColGIDs
            error = true
        end
    end

    return (error, BlockMap(numGlobalElements(domMap), myColumns, getComm(domMap)))
end


text:
using Test

# temporarily moved to start (for debugging)
include("testSpecialEuclidean2Mani.jl")

include("testFactorGradients.jl")
include("testpartialconstraint.jl")
include("testPartialNH.jl")

include("testSphereMani.jl")
include("testSpecialOrthogonalMani.jl")

include("testCliqSolveDbgUtils.jl")
include("TestModuleFunctions.jl")
include("testApproxConv.jl")
include("testCompareVariablesFactors.jl")
include("typeReturnMemRef.jl")
include("basicGraphsOperations.jl")
include("testMixturePrior.jl")
include("testGradientUtils.jl")
include("testPartialFactors.jl")
include("testPartialPrior.jl")
include("testSaveLoadDFG.jl")
include("testJunctionTreeConstruction.jl")
include("testBayesTreeiSAM2Example.jl")
include("testTreeFunctions.jl")

#FIXME fails on MetaBayesTree
include("testTreeSaveLoad.jl")

include("testSpecialSampler.jl")
include("saveconvertertypes.jl")
include("testgraphpackingconverters.jl")
include("testNLsolve.jl")
include("testCommonConvWrapper.jl")
include("testBasicForwardConvolve.jl")
include("testFactorMetadata.jl")
include("testStateMachine.jl")
include("testBasicCSM.jl")
include("testCliqueFactors.jl")
include("testCcolamdOrdering.jl")
include("testBasicGraphs.jl")
include("testDefaultDeconv.jl")
include("testJointEnforcement.jl")
include("testHasPriors913.jl")
include("testInitVariableOrder.jl")
include("testTreeMessageUtils.jl")
include("testCSMMonitor.jl")
include("testExpXstroke.jl")
include("testBasicRecycling.jl")
include("testSkipUpDown.jl")
include("testlocalconstraintexamples.jl")
include("testBasicTreeInit.jl")
include("testSolveOrphanedFG.jl")
include("testSolveSetPPE.jl")
include("testSolveKey.jl")
include("testEuclidDistance.jl")
include("priorusetest.jl")
include("testnullhypothesis.jl")
include("testVariousNSolveSize.jl")
include("testExplicitMultihypo.jl")
include("TestCSMMultihypo.jl")
include("testMultihypoFMD.jl")
include("testMultiHypo2Door.jl")
include("testMultimodal1D.jl")
include("testMultihypoAndChain.jl")
include("testMultithreaded.jl")
include("testmultihypothesisapi.jl")
include("fourdoortest.jl")
include("testCircular.jl")
include("testMixtureLinearConditional.jl")
include("testFluxModelsDistribution.jl")
include("testAnalysisTools.jl")
include("testDERelative.jl")

include("testBasicParametric.jl")
include("testMixtureParametric.jl")

# dont run test on ARM, as per issue #527
if Base.Sys.ARCH in [:x86_64;]
  include("testTexTreeIllustration.jl")
end

include("testMultiprocess.jl")
include("testDeadReckoningTether.jl")

#


text:
using AbstractTensors
using Test, LinearAlgebra, DirectSum

# example data
struct SpecialTensor{V} <: TensorAlgebra{V} end
a,b = SpecialTensor{ℝ}(), SpecialTensor{ℝ'}()
@test ndims(+(a)) == ndims(b)

## tensor operation (trivial test)
op(s::SpecialTensor{V},::SpecialTensor{V}) where V = s
op(a::TensorAlgebra{V},b::TensorAlgebra{W}) where {V,W} = interop(op,a,b)
(W::Signature)(s::SpecialTensor{V}) where V = SpecialTensor{W}()
@test Manifold(op(a,b)) == ℝ⊕ℝ'
@test Manifold(interop(op,a,b)) == ℝ⊕ℝ'
@test Manifold(op(a,a)) == ℝ
@test Manifold(interop(op,a,a)) == ℝ
Base.:+(s::SpecialTensor{V},::SpecialTensor{V}) where V = s
@test Manifold(+(a,b)) == ℝ⊕ℝ'
@test Manifold(interop(+,a,b)) == ℝ⊕ℝ'
@test Manifold(+(a,a)) == ℝ
@test Manifold(interop(+,a,a)) == ℝ

## tensor pseudoscalar (trivial test)
op(a::TensorAlgebra{V},b::UniformScaling) where V = op(a,V(b))
op(a::UniformScaling,b::TensorAlgebra{V}) where V = op(V(a),b)
(W::Signature)(s::UniformScaling) where V = SpecialTensor{W}()
@test Manifold(op(a,LinearAlgebra.I)) == ℝ
@test Manifold(op(LinearAlgebra.I,a)) == ℝ
@test Manifold(+(a,LinearAlgebra.I)) == ℝ
@test Manifold(+(LinearAlgebra.I,a)) == ℝ

## tensor evaluation (trivial test)
(a::SpecialTensor{V})(b::SpecialTensor{V}) where V = a
(a::SpecialTensor{W})(b::SpecialTensor{V}) where {V,W} = interform(a,b)
@test Manifold(b(a)) == ℝ⊕ℝ'
@test Manifold(interform(b,a)) == ℝ⊕ℝ'
@test Manifold(a(a)) == ℝ
@test Manifold(interform(a,a)) == ℝ

## algebraic tests
@test !I == 1
@test (2)⁻¹ == 1/2
@test (im)ǂ == -im
@test (sqrt(2))ˣ == sqrt(2)
@test (sqrt(2))₊ == sqrt(2)
@test (sqrt(2))₋ == 0


text:
module ROOTDataFrames

#include("log.jl")

#maximum size of preallocation for vector branches
const MAX_SIZE = 50

#print debugging info
const DEBUG = false

using ROOT, DataFrames, DataArrays
import Base.length, Base.getindex

import DataFrames.nrow, DataFrames.size, DataFrames.index

#holds data from TBranch
type BranchValue{T, N}
    data::Vector{T} #buffer
    size::Symbol #branch size (number or other branch name)
    branch::TBranch
end

BranchValue{T<:Real,N <: Symbol}(::Type{T}, s::N, br::TBranch=TBranch(C_NULL)) =
    BranchValue{T,N}(zeros(T, MAX_SIZE), s, br)

BranchValue{T<:AbstractVector, N <: Symbol}(::Type{T}, s::N, br::TBranch=TBranch(C_NULL)) =
    BranchValue{eltype(T),N}(zeros(eltype(T), MAX_SIZE), s, br)

#static array
BranchValue{T<:Real, N}(::Type{T}, s::Type{Val{N}}, br::TBranch=TBranch(C_NULL)) =
    BranchValue{T,N}(zeros(T, N), :nothing, br::TBranch)

#static array from vector
BranchValue{T<:AbstractVector, N}(::Type{T}, s::Type{Val{N}}, br::TBranch=TBranch(C_NULL)) =
    BranchValue{eltype(T),N}(zeros(eltype(T), N), :nothing, br::TBranch)

getsize{T, N}(bval::BranchValue{T, N}) = N
getsize{T, N <: Symbol}(bval::BranchValue{T, N}) = bval.size

import Base.eltype
# eltype{T, N <: Symbol}(bval::BranchValue{T, N}) = Vector{T}
# eltype{T, N}(bval::BranchValue{T, N}) = T
# eltype{T}(bval::BranchValue{T, 1}) = T
eltype{T, N}(bval::BranchValue{T, N}) = T

datatype{T, N <: Symbol}(bval::BranchValue{T, N}) = Vector{T}
datatype{T, N}(bval::BranchValue{T, N}) = Vector{T}
datatype{T}(bval::BranchValue{T, 1}) = T

#return value held in branch buffer
value{T, N <: Symbol}(bval::BranchValue{T, N}) = bval.data::Vector{T}
value{T, N}(bval::BranchValue{T, N}) = bval.data::Vector{T}
value{T}(bval::BranchValue{T, 1}) = bval.data[1]::T

import Base.call
call{T}(bval::BranchValue{T, 1}) = bval.data[1]::T
call{T, N}(bval::BranchValue{T, N}) = bval.data[1:N]::Vector{T}
call{T, N <: Symbol}(bval::BranchValue{T, N}) = bval.data::Vector{T}

#makes a class that represents the structure of a tree
function makeclass(
    names::Vector{Symbol},
    types::Vector{BranchValue},
    )
    typename = gensym("RowData")
    ret = :(
        type $typename

        $typename() = new()
        end
    )
    #members
    for (name, typ) in zip(names, types)
        push!(ret.args[3].args, Expr(:(::), name, typeof(typ)))
    end

    #constructor
    for (name, typ) in zip(names, types)
        leafsize = getsize(typ)
        if typeof(leafsize) <: Real #make a scalar
            push!(ret.args[3].args[2].args[2].args[2].args,
                Expr(:call, :BranchValue, eltype(typ), Val{leafsize})
            )
        else #variable length branch
            push!(ret.args[3].args[2].args[2].args[2].args,
                Expr(:call, :BranchValue, eltype(typ), QuoteNode(leafsize))
            )
        end
    end
    eval(ret)

    return eval(typename)
end

type TreeDataFrame{T} <: AbstractDataFrame
    tf::TObjectA
    tt::TTreeA
    row::T
    rowindex::Int64
    treeindex::Int64
end

import DataFrames.showcols
function showcols(io::IO, df::TreeDataFrame)
    println(io, summary(df))
    metadata = DataFrame(
        Name = names(df),
        Eltype = eltypes(df),
        sizes = [getsize(getfield(df.row, n)) for n in names(df)],
        #Missing = colmissing(df)
    )
    showall(io, metadata, true, symbol("Col #"), false)
    return
end
import DataFrames.show
function show(io::IO, df::TreeDataFrame)
    showcols(io, df)
    return
end

function update_branches!{T}(df::TreeDataFrame{T})
    for field in fieldnames(df.row)
        br = GetBranch(df.tt, string(field))
        getfield(df.row, field).branch = br
    end
end

function TreeDataFrame(fns::AbstractVector; treename="dataframe")

    #@debugprint "opening $fns:$treename"
    #DEBUG && println("[TreeDataFrame] opening $fns:$treename")
    tch = TChain(treename)

    for fn in fns
        #@debugprint "adding file $fn"
        #-1 -> load number of entries from file
        AddFile(tch, convert(ASCIIString, fn), -1)
        nentries = GetEntries(tch)
        #@debugprint "N=$nentries"
    end
    #@debugprint "TChain $tch created"
    #tf = TFile(fn)
    #tt = root_cast(TTree, Get(root_cast(ROOT.TDirectory, tf), "dataframe"))

    branches = GetListOfBranches(tch);
    branches = [root_cast(TBranch, branches[i]) for i=1:length(branches)];
    #@debugprint "TTree has $(length(branches)) branches"
    branches_d = Dict()
    branch_names = Symbol[]
    branch_types = BranchValue[]
    for b::TBranch in branches
        #println("creating branch $b")
        @assert !is_null(b)
        name = GetName(root_cast(TObject, b)) |> bytestring |> symbol;
        branches_d[name] = b
        push!(branch_names, name)
    end

    row_types = Dict{Symbol, BranchValue}()
    for name in branch_names
        #println("branch $k")
        leaves = GetListOfLeaves(branches_d[name])
        if length(leaves)!=1
            warn("branch=$k, nleaf=$(length(leaves)), skipping")
            continue
        end

        leaf = root_cast(TLeaf, leaves[1])
        #Leaf data type in ROOT types (Double_t etc)
        leaftype = GetTypeName(leaf)|>bytestring|>parse

        leaf_staticsize = ROOT.GetLenStatic(leaf)
        leafsize = 1
        if leaf_staticsize!=1
            leafsize = leaf_staticsize
            #warn("$name, nleaf size == $(leaf_staticsize)")
            #continue
        end

        #Branch is variable size, size is determined by another branch
        leafcount = ROOT.GetLeafCount(leaf)
        if leafcount != C_NULL
            leafcount = TLeaf(leafcount)
            leafname = leafcount |> GetName |> bytestring
            leafsize = symbol(leafname)
            #warn("$name, nleaf size dynamic $bname")
            #continue
        end

        if !haskey(ROOT.type_replacement, leaftype)
            warn("branch=$name with type=$leaftype does not have a julia type replacement, skipping")
            continue
        end

        #convert leaf type to julia type
        leaftype = eval(ROOT.type_replacement[leaftype])::Type

        if typeof(leafsize) <: Integer
            leafsize = Val{leafsize}
        end

        push!(
            branch_types,
            BranchValue(
                leaftype, leafsize, branches_d[name]
            )
        )
    end

    rowtype = makeclass(branch_names, branch_types)
    rowdata = rowtype()

    for field in fieldnames(rowdata)
        typ = fieldtype(rowtype, field)
        #getfield(rowdata, field).branch = branches_d[field]
        SetBranchAddress(tch, string(field), getfield(rowdata, field).data)
    end

    df = TreeDataFrame{rowtype}(
        TObject(C_NULL),
        tch,
        rowdata,
        1,
        1
    )
    update_branches!(df)
    return df
end

TreeDataFrame(fn::AbstractString) = TreeDataFrame([ascii(fn)])

function Base.length(t::TreeDataFrame)
    if t.tt != C_NULL
        @assert !is_null(t.tt)
        return GetEntries(t.tt)
    else
        return 0
    end
end
Base.size(df::TreeDataFrame) = (nrow(df), ncol(df))
Base.size(df::TreeDataFrame, n) = size(df)[n]

function load_row(df::TreeDataFrame, i::Integer)
    localentry = LoadTree(df.tt, i-1)
    localentry >= 0 || error("could not get entry $i")
    return GetEntry(df.tt, localentry)
end

function update!(df::TreeDataFrame, i::Int)
    treeindex = GetTreeNumber(df.tt) + 1
    if treeindex != df.treeindex
        update_branches!(df)
        df.treeindex = treeindex
    end
    df.rowindex = i
end

function load_row(df::TreeDataFrame, i::Integer, cols::Vector{Symbol})
    @assert !is_null(df.tt)

    localentry = LoadTree(df.tt, i-1)
    localentry >= 0 || error("could not get entry $i")
    update!(df, i)

    ntot = 0
    for col in cols
        br = getfield(df.row, col).branch
        @assert !is_null(br)
        ntot += GetEntry(br, localentry)
    end
    return ntot
end

function load_row(df::TreeDataFrame, i::Integer)
    localentry = LoadTree(df.tt, i-1)
    localentry >= 0 || error("could not get entry $i")
    update!(df, i)
    return GetEntry(df.tt, localentry)
end

#return value of scalar branch
@generated function Base.getindex{T <: Any, X <: Any, N}(
    df::TreeDataFrame{T},
    i::Int64,
    s::Symbol,
    getentry::Bool=true,
    ft::Type{BranchValue{X, N}}=fieldtype(T, s),
    )

    if N == Symbol
        ex = quote
            #get dynamic branch size
            sizebranch = getsize(getfield(df.row, s))
            if getentry
                load_row(df, i, [s, sizebranch]) > 0 || error("could not load row $i")
            end
            _size = value(getfield(df.row, sizebranch))
            ret = getfield(df.row, s)
            return value(ret)[1:_size]::Vector{X}
        end
    else
        if N == 1
            ret  = :(value(ret)::X)
        elseif N > 1
            ret = :(value(ret)[1:$N]::Vector{X})
        end

        ex = quote
            if getentry
                load_row(df, i, [s])>0 || error("could not load row $i")
            end
            ret = getfield(df.row, s)
            return $ret
        end
    end

    return ex

end
#
# #return value of scalar branch
# function Base.getindex{T <: Any, X <: Any}(
#     df::TreeDataFrame{T},
#     i::Int64,
#     s::Symbol,
#     getentry::Bool=true,
#     ::Type{BranchValue{X, N}}=fieldtype(T, s),
#     )
#     println("scalar called ", s)
#     if getentry
#         load_row(df, i, [s])>0 || error("could not load row $i")
#     end
#     ret = getfield(df.row, s)
#     @show ret
#     return value(ret)[1:_size]
# end

#return value of dynamic branch
# function Base.getindex{T <: Any, R <: Symbol, X <: Any}(
#     df::TreeDataFrame{T},
#     i::Int64,
#     s::Symbol,
#     getentry::Bool=true,
#     _size::Int64=-1,
#     ft::Type{BranchValue{X, R}}=fieldtype(T, s),
#     )
#     println(ft)
#     #get dynamic branch size
#     if _size == -1
#         sizebranch = getsize(getfield(df.row, s))
#         if getentry
#             load_row(df, i) > 0 || error("could not load row $i")
#         end
#         _size = value(getfield(df.row, sizebranch))
#     end
#     ret = getfield(df.row, s)
#     @show ret, _size
#     return value(ret)[1:_size]::Vector{X}
# end

import DataFrames.nrow, DataFrames.ncol
DataFrames.nrow(df::TreeDataFrame) = length(df)
DataFrames.ncol(df::TreeDataFrame) = length(names(df))

import Base.names
Base.names(df::TreeDataFrame) = fieldnames(df.row)

import DataFrames.eltypes
DataFrames.eltypes(df::TreeDataFrame) = Type[datatype(getfield(df.row, n)) for n in names(df)]

set_branch_status!(df, pattern, status) = SetBranchStatus(
    df.tt, pattern, status, convert(Ptr{Cuint}, 0)
)

function enable_branches(df, brs)
    #println("enabling branches $brs")
    #SetCacheSize(df.tt, 0)
    #SetCacheSize(df.tt, 256 * 1024 * 1024)
    set_branch_status!(df, "*", false)
    for b in brs
        set_branch_status!(df, "$b", true)
        #AddBranchToCache(df.tt, "$b") #cache does not work properly, causes crash
    end
end

function Base.getindex{T}(
    df::TreeDataFrame{T},
    mask::AbstractVector,
    cols::AbstractVector{Symbol}
    )
    length(mask) == nrow(df) || error("mask=$(length(mask)), nrow=$(nrow(df))")
    enable_branches(df, ["$(s)*" for s in cols])

    #prepare output array
    names_types = Dict{Symbol, DataType}([n => eltype(getfield(df.row, n)) for n in names(df)])
    sizes = Dict{Symbol, Any}([n => getsize(getfield(df.row, n)) for n in names(df)])

    #output vector length
    const n = sum(mask)

    #need to disable GC to prevent silent corruption and crash
    #still un-understood
    #gc_disable()

    arrlist = Any[]
    coltypes = Type[]
    for col in cols
        _size = sizes[col]
        _type = names_types[col]
        if typeof(_size) <: Symbol
            _size = MAX_SIZE
        end
        if _size > 1
            _type = DataVector{_type}
        end
        push!(coltypes, _type)
        a = DataFrames.DataArray(_type, n)
        push!(arrlist, a)
    end

    #construct output dataframe
    const ret = DataFrame(
        arrlist,
        DataFrames.Index(cols)
    )
    j = 1

    t0 = time()
    nloaded = 0
    for i=1:nrow(df)
        (!isna(mask[i]) && mask[i]) || continue
        nloaded += load_row(df, i, cols)
        for (icol, col) in enumerate(cols)
            ret[j, col] = df[i, col]
        end
        j += 1
    end
    t1 = time()
    enable_branches(df, ["*"])
    #gc_enable()
    return ret
end

Base.getindex(df::TreeDataFrame, s::Symbol) = df[[s]][s]
Base.getindex(df::TreeDataFrame, ss::AbstractVector{Symbol}) = df[[true for i=1:nrow(df)], ss]
Base.getindex(df::TreeDataFrame, inds::AbstractVector{Int64}, ss::AbstractVector{Symbol}) = df[[i in inds for i=1:nrow(df)], ss]
Base.getindex(df::TreeDataFrame, i::Int64, ss::AbstractVector{Symbol}) = df[[i], ss]
Base.getindex(df::TreeDataFrame, mask::AbstractVector, s::Symbol) = df[mask, [s]][s]

TreeDataFrame(
    fn::ASCIIString,
    names::AbstractVector,
    types::AbstractVector,
    colsizes::AbstractVector
    ;kwargs...
    ) = TreeDataFrame(
        ascii(fn),
        [symbol(n) for n in names],
        BranchValue[BranchValue(dt, cs) for (dt, cs) in
        zip(types, colsizes)]; kwargs...
)

function TreeDataFrame(
    fn::ASCIIString,
    names::Vector{Symbol},
    types::Vector{BranchValue}
    ;treename="dataframe"
    )

    rowdatatype = makeclass(names, types)
    rowdata = rowdatatype()

    tf = TFile(fn, "RECREATE")

    #make a TTree pointer
    const tree = TTree(
        treename, treename
    )

    for (ibranch, (name, typ)) in enumerate(
        zip(names, types)
        )
        eltyp = eltype(typ)

        sz = getsize(typ)
        rtype = SHORT_TYPEMAP[eltyp]

        brstring = ""

        #variable size branch
        if isa(sz, Symbol)
            brstring = "$name[$sz]/$rtype"
        #fixed size branch
        elseif sz > 1
            brstring = "$name[$sz]/$rtype"
        # scalar branch
        else
            brstring = "$name/$rtype"
        end

        br = Branch(
            tree, string(name),
            convert(Ptr{Void}, pointer(getfield(rowdata, name).data)),
            brstring
        )
        assert(root_pointer(br) != C_NULL)
        getfield(rowdata, name).branch = br
    end

    dtf = TreeDataFrame(
        tf,
        tree,
        rowdata,
        1,
        1
    )
end

import Base.setindex!
function Base.setindex!{T, K <: Real}(
    df::TreeDataFrame{T}, val::K, i::Integer, col::Symbol
    )
    getfield(df.row, col).data[1] = val
end

function Base.setindex!{T, K <: Real}(
    df::TreeDataFrame{T}, val::Vector{K}, i::Integer, col::Symbol
    )
    v = val[:]
    lv = length(v)
    assert(lv < MAX_SIZE)
    #fill with zeros
    getfield(df.row, col).data[:] = zero(eltype(getfield(df.row, col)))

    #fill with correct length
    getfield(df.row, col).data[1:lv] = v[1:lv]
end

function TreeDataFrame(df::AbstractDataFrame, filename; size_branches::Associative=Dict(), treename="dataframe")
    colsizes = Any[]

    _eltypes = eltypes(df)
    _names = names(df)
    _sizes = Array(Any, ncol(df))
    for i=1:ncol(df)
        _el = _eltypes[i]
        _cname = _names[i]

        size_branch = get(size_branches, _cname, symbol(:n, _cname))

        idx_size = findfirst(_names, size_branch)
        if idx_size > 0
            idx_size > i && error("found size branch $size_branch for $_cname, but must be declared before.")
            _size = size_branch

        #constant size
        elseif _el <: AbstractVector && eltype(_el) <: Real
            maxsize = maximum(map(x->length(df[x, _cname]), 1:nrow(df)))
            assert(maxsize < MAX_SIZE)
            _size = Val{maxsize}
        #scalar
        else
            _size = Val{1}
        end
        _sizes[i] = _size
    end
    dtf = TreeDataFrame(filename, _names, _eltypes, _sizes;treename=treename)

    for i=1:nrow(df)
        for j=1:ncol(df)
            dtf[i, _names[j]] = df[i, _names[j]]
        end

        Fill(dtf.tt)
    end

    Write(dtf.tt)
    Close(dtf.tf)
    return TreeDataFrame([filename]; treename=treename)
end

function writetree(
    fn::ASCIIString,
    df::AbstractDataFrame;
    size_branches=Dict(),
    treename="dataframe"
    )
    TreeDataFrame(df, fn; size_branches=size_branches, treename=treename)
end

is_selected(i, row, selector::Function) = selector(row)
is_selected(i, row, selector::AbstractArray) = selector[i]

function with{T, R <: Real}(
    df::TreeDataFrame{T}, #input data
    func::Function, #process func ev->R
    selector, #selector func ev->Bool
    branches::Vector{Symbol}, #branches to enable
    rng=1:length(df), #rows to process
    ::Type{R}=Float64 #return type
    )
    enable_branches(df, branches)
    ntot = 0
    tic()
    ret = zeros(R, length(rng))
    bitmask = BitArray(length(rng))
    for i=rng
        ntot += load_row(df, i, branches)

        @inbounds const sel = is_selected(i, df.row, selector)::Bool
        @inbounds bitmask[i] = sel
        if sel
            @inbounds const res = func(df.row)::R
            @inbounds ret[i] = res
        end
    end
    nmb = float(ntot) / 1024.0 / 1024.0
    dt = toq()
    speed = ntot/dt / 1024.0 / 1024.0
    #println("Read ", round(nmb), " Mb in ", dt, " s speed=", round(speed,2), " Mb/s")
    #re-enable branches
    enable_branches(df, ["*"])
    return ret[bitmask]
end

function loop{T}(df::TreeDataFrame{T}, f1::Function, f2::Function, branches::Vector{Symbol}, rng=1:length(df))
    enable_branches(df, branches)
    ntot = 0
    tic()
    for i=rng
        ntot += load_row(df, i, branches)
        @inbounds const sel = f2(df.row)::Bool
        if sel
            @inbounds f1(df.row)
        end
    end
    nmb = float(ntot) / 1024.0 / 1024.0
    dt = toq()
    speed = ntot/dt / 1024.0 / 1024.0
    println("Read ", round(nmb), " Mb in ", dt, " s speed=", round(speed,2), " Mb/s")
    enable_branches(df, ["*"])
end

export value
export with, loop
export BranchValue
export writetree, TreeDataFrame
export writetree_temp
export load_row
export enable_branches
end # module


text:
#!/usr/bin/env julia

using SpectralElements, OrdinaryDiffEq, Flux
using LinearAlgebra,Plots
#----------------------------------#
Ex = 2; nr1 = 8;
Ey = 2; ns1 = 8;

m1 = Mesh(nr1,ns1,Ex,Ey,[true,true])
bc = ['D','D','D','D']
M  = generateMask(bc,m1)
opM(u) = m1.Bi .* u
opB(u) = gatherScatter(mass(u,m1),m1)
#----------------------------------#
kx=1.0
ky=1.0
kt=1.0
ν =1.0

utrue(x,y,t) = @. sin(kx*pi*x)*sin.(ky*pi*y)*cos(kt*pi*t)
setIC(x,y)   = utrue(x,y,0.0)

function setForcing(x,y,t)
    ut = utrue(x,y,t)
    f  = ut*((kx^2+ky^2)*pi^2*ν)
    f -= sin(kx*pi*x)*sin(ky*pi*y)*sin(kt*pi*t)*(kt*pi)
#   f  = zero(x)
end

#----------------------------------#
function cond(u,t,integrator)
    istep = integrator.iter
    cond = (istep % 1) == 0
    return cond
end

function affect!(integrator)
    istep = integrator.iter
    u = integrator.u
    p = integrator.p
    t = integrator.t
    dt = integrator.dt

    ut = utrue(m1.x,m1.y,t)
    er = norm(ut-u,Inf)
    println("Step=$istep, Time=$t, dt=$dt, er=$er")
    plt = meshplt(u,m1)
    plt = plot!(zlims=(-1,1))
    display(plt)
    return
end

cb = DiscreteCallback(cond,affect!,save_positions=(false,false))

#----------------------------------#
function dudt!(du,u,p,t)
    f = setForcing.(m1.x,m1.y,t)

    rhs = -lapl(u,m1) + mass(f,m1)
    rhs = gatherScatter(rhs,m1)

    pcg!(du,rhs,opB;opM=opM,mult=m1.mult,ifv=false)

    return du
end

u0   = setIC(m1.x,m1.y)
dt   = 0.01
tspn = (0.0,1.0)

func = ODEFunction(dudt!)
prob = ODEProblem(func,u0,tspn)
sol  = solve(prob, Rodas5(); saveat=0.1, callback=cb)
using LinearSolve
sol  = solve(prob, Rodas5(linsolve=LUFactorization()); saveat=0.1, callback=cb)

err = sol.u - [utrue.(m1.x,m1.y,sol.t[i]) for i=axes(sol.t,1)]
ee = maximum.(abs.(err[i]) for i=axes(sol.t,1))
#----------------------------------#
return ee


text:
@testset "File formats" begin

    @testset "avr" begin
        include("avr.jl")
    end
    @testset "bsa" begin
        include("bsa.jl")
    end
    @testset "dat" begin
        include("dat.jl")
    end
    @testset "elp" begin
        include("elp.jl")
    end
    @testset "evt" begin
        include("evt.jl")
    end
    @testset "rba" begin
        include("rba.jl")
    end
    @testset "sfp" begin
        include("rba.jl")
    end
end


text:
###################################################################################
# Antony Schutz 2015, ANR - POLCA - 2016
###################################################################################
# OIFITS RELATED - Read files, create quantities of interest
###################################################################################

# path: read OIFITS data from files included in folder OIDATA.Folder
# OIDATA: structure containing all OIDATA and users informations
# indfile: restrict the number of files (as listed by dir() or ls() )
# ex: indfile = [1 3 5] will only read information from the 1st, the 3rd and the 5th file
# if empty read all the files
# indwvl: restrict the analyzed wavelengths restricted by index, ex:
# indwvl=1:20 --- indwvl= [1 5 7] --- indwvl = 1:5:20
#
# - - - - - - - - - - - - - - - - - - - restreindre les longueurs d'ondes par des valeurs en metre
#
function readoifits(OIDATA::PAINTER_Input,indfile=[],indwvl=[])

# Check for DP and T3 must be in data at least
    OIDATA.isDP = 1

# Read Directory
    Folder    = OIDATA.Folder
    CDir      = readdir(Folder)

# check files fits and oifits
    fits      = filter(r"\.fits$", CDir)
    oifits    = filter(r"\.oifits$", CDir)
    CDir      = vcat(fits,oifits)
    NbFiles   = length(CDir)

    if NbFiles == 0
      error("No fits or oifits files in folder: $Folder")
    end

# files reduction, only files indexed by "indfile" are read
    if isempty(indfile)
        indfile = 1:NbFiles
    end
    OIDATA.indwvl = indwvl

# temporary variable for closure index
    LT3           = 0
    LFil          = 0
# loop on files
# Extract Data from OIFITS or FITS Files

    for n in indfile
        name = string(Folder, "/", CDir[n])
        println(name)
        master = OIFITS.load(name)
# # # OI wavelength
        nw = 0
        for db in OIFITS.select(master, "OI_WAVELENGTH")
            nw = length(OIFITS.get_eff_wave(db))
        end
# Variables initialization
        if n == minimum(indfile)
            if isempty(OIDATA.indwvl)
                OIDATA.indwvl = 1:nw
            end
# Initialisation of "empty" array
            OIDATA.Closure_index = Array(Int, 0, 3)
            OIDATA.U = Array(Float64, 0, nw)
            OIDATA.V = Array(Float64, 0, nw)
            OIDATA.P = Array(Float64, 0, nw)
            OIDATA.W = Array(Float64, 0, nw)
            OIDATA.T3 = Array(Float64, 0, nw)
            OIDATA.T3err = Array(Float64, 0, nw)
            OIDATA.DP = Array(Float64, 0, nw)
            OIDATA.DPerr = Array(Float64, 0, nw)
            OIDATA.wvl = Array(Float64, 0, nw)
        end

# # # OI wavelength
        for db in OIFITS.select(master, "OI_WAVELENGTH")
            OIDATA.wvl = vcat(OIDATA.wvl,OIFITS.get_eff_wave(db)')
        end
# # # OI VIS
        if isempty(OIFITS.select(master, "OI_VIS"))
            println("OI_VIS field is missing, PAINTER needs differential visibilities")
            println(" --- will try to work only with closure phases ")
            OIDATA.isDP = 0
        else

            for dbvis1 in OIFITS.select(master, "OI_VIS")
                OIDATA.DP = vcat(OIDATA.DP, OIFITS.get_visphi(dbvis1)')
                OIDATA.DPerr = vcat(OIDATA.DPerr, OIFITS.get_visphierr(dbvis1)')
            end

        end

# # OI VIS 2
        if isempty(OIFITS.select(master, "OI_VIS2"))
            error("OI_VIS2 field is missing, PAINTER needs squared visibilities")

        else
            U_v21 = Array(Float64, 0)

            for dbvis2 in OIFITS.select(master,  "OI_VIS2")
                U,V  = spatialfrequencies(OIFITS.get_ucoord(dbvis2), OIFITS.get_vcoord(dbvis2), OIFITS.get_eff_wave(dbvis2))
                U_v21 = vcat(U_v21, OIFITS.get_ucoord(dbvis2)) # for closure index
                OIDATA.U = vcat(OIDATA.U, U)
                OIDATA.V = vcat(OIDATA.V, V)
                OIDATA.P = vcat(OIDATA.P, OIFITS.get_vis2data(dbvis2)')
                OIDATA.W = vcat(OIDATA.W, OIFITS.get_vis2err(dbvis2)')
            end
        end

# # OI T3
        if isempty(OIFITS.select(master, "OI_T3"))
            error("OI_T3 field is missing, PAINTER needs closure phases")

        else
            U_t31 = Array(Float64, 0)
            U_t32 = Array(Float64, 0)
            for dbvis3 in OIFITS.select(master, "OI_T3")
                OIDATA.T3 = vcat(OIDATA.T3, OIFITS.get_t3phi(dbvis3)')
                OIDATA.T3err = vcat(OIDATA.T3err, OIFITS.get_t3phierr(dbvis3)')
# Closure index
                U_t31 = OIFITS.get_u1coord(dbvis3)
                U_t32 = OIFITS.get_u2coord(dbvis3)
                (OIDATA.Closure_index, LT3, LFil) = findclosureindex(U_v21, U_t31, U_t32, OIDATA.Closure_index, LT3, LFil)
            end
        end
    end

## FINAL OUTPUT - Wavelength reduction
    unitconvert = pi / 180
    OIDATA.FilesName = CDir
    OIDATA.indfile = indfile
    OIDATA.U = OIDATA.U[:, OIDATA.indwvl]
    OIDATA.V = OIDATA.V[:, OIDATA.indwvl]
    OIDATA.P = OIDATA.P[:, OIDATA.indwvl]
    OIDATA.W = OIDATA.W[:, OIDATA.indwvl]

    OIDATA.T3 = OIDATA.T3[:, OIDATA.indwvl] * unitconvert
    OIDATA.T3err = OIDATA.T3err[:, OIDATA.indwvl] * unitconvert
    if OIDATA.isDP == 1
        OIDATA.DP = OIDATA.DP[:, OIDATA.indwvl] * unitconvert
        OIDATA.DPerr = OIDATA.DPerr[:, OIDATA.indwvl] * unitconvert
    end
    OIDATA.wvl = vec(mean(OIDATA.wvl[:, OIDATA.indwvl], 1))
# # compute Differential Phases as needed for PAINTER and adjust angle unit
    OIDATA.nb = size(OIDATA.U, 1)
    OIDATA.nw = size(OIDATA.U, 2)

# # Create cluster of independent CLosure Phases and associated
    OIDATA = independentT3data(OIDATA)

    return OIDATA
end

function independentT3data(OIDATA::PAINTER_Input)
  Closure_index = OIDATA.Closure_index
  Cluster = independentT3(Closure_index)
  baseNb = basesincluster(Cluster)
  orderedCluster = makeclusterordered(Cluster,baseNb)
  rowt3 = t3row(Closure_index,Cluster)
  Nt3indep = length(baseNb)
  Xi = Dict{}()
  K = Dict{}()
  H = Dict{}()
  for n in 1:Nt3indep
      HDP = phasetophasediff(orderedCluster[n], OIDATA.nw, length(baseNb[n]), 0, 1, OIDATA.dptype, OIDATA.dpprm)
      T3 = OIDATA.T3[rowt3[n],:]
      T3err = OIDATA.T3err[rowt3[n],:]
      if OIDATA.isDP == 1
          DP = OIDATA.DP[baseNb[n],:]
          DPerr = OIDATA.DPerr[baseNb[n],:]
          Xi[n] = vcat(vec(T3),HDP*vec(DP)  )
          K[n] = vcat(vec(T3err),abs(HDP)*vec(DPerr)  )
      else
          Xi[n] = vec(T3)
          K[n] = vec(T3err)
      end
  end
  OIDATA.Xi = Xi
  OIDATA.K = K
  OIDATA.baseNb = baseNb
  OIDATA.orderedCluster= orderedCluster
  return OIDATA
end
###################################################################################
# FUNCTIONS FOR OIFITS DATA TRANSFORMATION
###################################################################################
# Spatial Frequencies construction
# ---------------------------------------------------------------------------------
function spatialfrequencies(U_Coord::Vector,V_Coord::Vector,wvl::Vector)
# Create spatiale frequencies from bases coordinates and wavelength
# (U,V) = SpatialFrequencies(U_Coord::Vector,V_Coord::Vector,wvl::Vector)
# U_Coord:  U coordinate of V2 (bases not spatial frequencies) vector (U not divided by wavelength)
# V_Coord:  V coordinate of V2 (bases not spatial frequencies) vector (U not divided by wavelength)
# wvl    :  vector of effective used wavelength
    U_3D = repmat(U_Coord, 1, length(wvl)) ./ repmat(wvl', length(U_Coord), 1)
    V_3D = repmat(V_Coord, 1, length(wvl)) ./ repmat(wvl', length(V_Coord), 1)
    return U_3D,V_3D
end
# ---------------------------------------------------------------------------------
# PHASE CLOSURE INDEX (sample's index) creation
# ---------------------------------------------------------------------------------
function findclosureindex(V2_U_Coord::Vector,T3_U1_Coord::Vector,T3_U2_Coord::Vector,Closure_index::Matrix,LT3::Integer,LFil::Integer)
# update Closure_index array for each files
# deduce the index in sample related to phases vector, from U1 and U2 coord of phi T3 to U coord of V2
# Closure_index = CreateClosureIndex(V2_U_Coord::Vector,T3_U1_Coord::Vector,T3_U2_Coord::Vector)
# V2_U_Coord:  U coordinate of V2 (bases not spatial frequencies) vector (U not divided by wavelength)
# T3_U1_Coord: U coordinate of T3 first base (not spatial frequencies) vector (U1 not divided by wavelength)
# T3_U2_Coord: U coordinate of T3 second base (not spatial frequencies) vector (U2 not divided by wavelength)
# U3 is given by default
# LT3 and LFil are used to increment correctly information from several files
# Closure_index: index in tabs of associated visibilities
    utmp = round(1000 * V2_U_Coord )
    U1 = round(1000 * T3_U1_Coord)
    U2 = round(1000 * T3_U2_Coord)
    U3 = round(1000 * (T3_U1_Coord + T3_U2_Coord))
    for k in 1:length(U1)
        pos1 = LT3 + find(x->(x == U1[k]), utmp)
        pos2 = LT3 + find(x->(x == U2[k]), utmp)
        pos3 = LT3 + find(x->(x == U3[k]), utmp)
        Closure_index = vcat(Closure_index, hcat(pos1, pos2, pos3))
    end
    LT3 = maximum(Closure_index[:])
    LFil = LFil + length(U1)
    return Closure_index, LT3, LFil
end
# ---------------------------------------------------------------------------------
# PAINTER used Differential Visibilities, creation from Differential Visibilities
# ---------------------------------------------------------------------------------
function diffphi(DiffPhi::Array, HDP::SparseMatrixCSC)
# DiffPhi: Differential visibilities as given in OIFITS Data
# DiffPhi is matrix (Nbase*Nwvl)
# All channels reference are assumed to be equal
# defined as in equation 21 of PAINTER (1)
    DiffPhiAB = HDP * vec(DiffPhi)
end
# ---------------------------
# All channels are assumed to be independent
# variance of sum of random variables = sum of random variables variance
function diffphierr(DiffPhierr::Array, HDP::SparseMatrixCSC)
    DiffPhiABerr = abs(HDP) * vec(DiffPhierr)
end
# ---------------------------------------------------------------------------------
# INDEPENDENT Phases Closures from index -- to put in PAINTEROIFITS.JL
# ---------------------------------------------------------------------------------
function independentT3(Closure_index::Matrix)
    Cluster = Dict{}()
    Index = collect(1:size(Closure_index,1))
    m=0
    while sum(Index)>0
        m+=1
        idxinit = find(Index.>0)[1]
        Cluster[m] = Closure_index[Index[idxinit],:]
        Index[idxinit] = 0
        nn = 0
        while true
            nn+=1
            init = vec(Cluster[m][nn,:])
            Cluster[m], Index = searchfromclusterinlist(Cluster[m],init,Index,Closure_index)
            l = size(Cluster[m],1)
            if  nn==l
                break
            end

          end
      end
    return Cluster
end
function searchfromclusterinlist(Cluster::Matrix,init::Vector,Index::Vector,Closure_index::Matrix)
    for n in Index
      if n>0
        if( !isempty( find( init[1] .== Closure_index[n,:] ) )
          ||!isempty( find( init[2] .== Closure_index[n,:] ) )
          ||!isempty( find( init[3] .== Closure_index[n,:] ) ) )
            Cluster = vcat( Cluster , Closure_index[n,:] )
            Index[n] = 0
        end
      end
    end
    return Cluster, Index
end
# ---------------------------------------------------------------------------------
# Number of bases in each cluster - tools to create smaller T3 and DP Matrices
# ---------------------------------------------------------------------------------
function basesincluster(Cluster::Dict)
    Nc = length(Cluster)
    basenb = Dict{}()
    for n in 1:Nc
        vecclust = sort( vec( Cluster[n] ) )
        indpos = vcat(1, diff( vecclust ).>0)
        basenb[n] = round(Int, vecclust[ (vecclust.*indpos).>0 ])
    end
    return basenb
end
# ---------------------------------------------------------------------------------
# Sort bases in Cluster from 1 to Nb base in Cluster
# ---------------------------------------------------------------------------------
function makeclusterordered(Cluster::Dict,baseNb::Dict)
  Nc = length(Cluster)
  theorderedcluster = Dict{}()
  for n in 1:Nc
      thecluster = Cluster[n]
      baseinthecluster = baseNb[n]
      numberofindependentbaseinthecluster = length(baseinthecluster)
      nt3 = size(thecluster,1)
      orderedbaseinthecluster = collect(1:numberofindependentbaseinthecluster)
      thevectorizedcluster = vec(thecluster)
      theorderedvectorizedcluster = zeros(thevectorizedcluster)
      for m in 1:numberofindependentbaseinthecluster
          cond = find( thevectorizedcluster .==  baseinthecluster[m])
          theorderedvectorizedcluster[cond] = orderedbaseinthecluster[m]
        end
      theorderedcluster[n] = reshape(theorderedvectorizedcluster,nt3,3)
  end
  return theorderedcluster
end
# ---------------------------------------------------------------------------------
# T3 rows related to H
# ---------------------------------------------------------------------------------
function t3row(Closure_index::Matrix,Cluster::Dict)
  Nc = length(Cluster)
  rowt3 = Dict{}()
  for n in 1:Nc
      thecluster = Cluster[n]
      lc = size(thecluster,1)
      rowt3[n] = zeros(lc)
      for m in 1:lc
          rowt3[n][m] = find( prod(Closure_index .== thecluster[m,:],2))[1]
      end
      rowt3[n] = round(Int,rowt3[n])
  end
  return rowt3
end