text:
using BlockArrays, FillArrays, OffsetArrays, Test, Base64
import BlockArrays: BlockIndex, BlockIndexRange, BlockSlice

@testset "Blocks" begin
    @test Int(Block(2)) === Integer(Block(2)) === Number(Block(2)) === 2
    @test Block((Block(3), Block(4))) === Block(3,4)

    @testset "Block arithmetic" begin
        @test +(Block(1)) == Block(1)
        @test -(Block(1)) == Block(-1)
        @test Block(2) + Block(1) == Block(3)
        @test Block(2) + 1 == Block(3)
        @test 2 + Block(1) == Block(3)
        @test Block(2) - Block(1) == Block(1)
        @test Block(2) - 1 == Block(1)
        @test 2 - Block(1) == Block(1)
        @test 2*Block(1) == Block(2)
        @test Block(1)*2 == Block(2)

        @test isless(Block(1), Block(2))
        @test !isless(Block(1), Block(1))
        @test !isless(Block(2), Block(1))
        @test Block(1) < Block(2)
        @test Block(1) ≤ Block(1)
        @test Block(2) > Block(1)
        @test Block(1) ≥ Block(1)
        @test min(Block(1), Block(2)) == Block(1)
        @test max(Block(1), Block(2)) == Block(2)

        @test +(Block(1,2)) == Block(1,2)
        @test -(Block(1,2)) == Block(-1,-2)
        @test Block(1,2) + Block(2,3) == Block(3,5)
        @test Block(1,2) + 1 == Block(2,3)
        @test 1 + Block(1,2) == Block(2,3)
        @test Block(2,3) - Block(1,2) == Block(1,1)
        @test Block(1,2) - 1 == Block(0,1)
        @test 1 - Block(1,2) == Block(0,-1)
        @test 2*Block(1,2) == Block(2,4)
        @test Block(1,2)*2 == Block(2,4)

        @test isless(Block(1,1), Block(2,2))
        @test isless(Block(1,1), Block(2,1))
        @test !isless(Block(1,1), Block(1,1))
        @test !isless(Block(2,1), Block(1,1))
        @test Block(1,1) < Block(2,1)
        @test Block(1,1) ≤ Block(1,1)
        @test Block(2,1) > Block(1,1)
        @test Block(1,1) ≥ Block(1,1)
        @test min(Block(1,2), Block(2,2)) == Block(1,2)
        @test max(Block(1,2), Block(2,2)) == Block(2,2)

        @test convert(Int, Block(2)) == 2
        @test convert(Float64, Block(2)) == 2.0

        @test_throws MethodError convert(Int, Block(2,1))
        @test convert(Tuple{Int,Int}, Block(2,1)) == (2,1)
        @test convert(Tuple{Float64,Int}, Block(2,1)) == (2.0,1)

        @test Block(1)[:] ≡ Block(1)[Base.Slice(1:2)] ≡ Block(1)
    end

    @testset "BlockIndex" begin
        @test Block(1)[1] == BlockIndex((1,),(1,))
        @test Block(1)[1:2] == BlockIndexRange(Block(1),(1:2,))
        @test Block(1,1)[1,1] == BlockIndex((1,1),(1,1))
        @test Block(1,1)[1:2,1:2] == BlockIndexRange(Block(1,1),(1:2,1:2))
        @test Block(1)[1:3][1:2] == BlockIndexRange(Block(1),1:2)
    end

    @testset "BlockRange" begin
        @test Block.(2:5) isa BlockRange
        @test Block.(Base.OneTo(5)) isa BlockRange
        @test Block.(2:5) == [Block(2),Block(3),Block(4),Block(5)]
        b = Block.(2:5)
        @test Int.(b) === 2:5
        @test Base.OneTo.(1:5) isa Vector{Base.OneTo{Int}} #98
    end
end

@testset "BlockedUnitRange" begin
    @testset "Block indexing" begin
        b = blockedrange([1,2,3])
        @test axes(b) == (b,)
        @test blockaxes(b,1) isa BlockRange

        @test @inferred(b[Block(1)]) == 1:1
        @test b[Block(2)] == 2:3
        @test b[Block(3)] == 4:6
        @test_throws BlockBoundsError b[Block(0)]
        @test_throws BlockBoundsError b[Block(4)]

        o = OffsetArray([2,2,3],-1:1)
        b = blockedrange(o)
        @test axes(b) == (b,)
        @test @inferred(b[Block(-1)]) == 1:2
        @test b[Block(0)] == 3:4
        @test b[Block(1)] == 5:7
        @test_throws BlockBoundsError b[Block(-2)]
        @test_throws BlockBoundsError b[Block(2)]

        b = BlockArrays._BlockedUnitRange(-1,[-1,1,4])
        @test axes(b,1) == blockedrange([1,2,3])
        @test b[Block(1)] == -1:-1
        @test b[Block(2)] == 0:1
        @test b[Block(3)] == 2:4
        @test_throws BlockBoundsError b[Block(0)]
        @test_throws BlockBoundsError b[Block(4)]

        o = OffsetArray([2,2,3],-1:1)
        b = BlockArrays._BlockedUnitRange(-3, cumsum(o) .- 4)
        @test axes(b,1) == blockedrange([2,2,3])
        @test b[Block(-1)] == -3:-2
        @test b[Block(0)] == -1:0
        @test b[Block(1)] == 1:3
        @test_throws BlockBoundsError b[Block(-2)]
        @test_throws BlockBoundsError b[Block(2)]

        b = blockedrange(Fill(3,1_000_000))
        @test b isa BlockedUnitRange{StepRange{Int,Int}}
        @test b[Block(100_000)] == 299_998:300_000
        @test_throws BlockBoundsError b[Block(0)]
        @test_throws BlockBoundsError b[Block(1_000_001)]
    end

    @testset "firsts/lasts/lengths" begin
        b = blockedrange([1,2,3])
        @test blockfirsts(b) == [1,2,4]
        @test blocklasts(b) == [1,3,6]
        @test blocklengths(b) == [1,2,3]
    end

    @testset "convert" begin
        b = blockedrange(Fill(2,3))
        c = blockedrange([2,2,2])
        @test convert(BlockedUnitRange, b) === b
        @test blockisequal(convert(BlockedUnitRange, Base.OneTo(5)), blockedrange([5]))
        @test blockisequal(convert(BlockedUnitRange, Base.Slice(Base.OneTo(5))), blockedrange([5]))
        @test blockisequal(convert(BlockedUnitRange, Base.IdentityUnitRange(-2:2)), BlockArrays._BlockedUnitRange(-2,[2]))
        @test convert(BlockedUnitRange{Vector{Int}}, c) === c
        @test blockisequal(convert(BlockedUnitRange{Vector{Int}}, b),b)
        @test blockisequal(convert(BlockedUnitRange{Vector{Int}}, Base.OneTo(5)), blockedrange([5]))
    end

    @testset "findblock" begin
        b = blockedrange([1,2,3])
        @test @inferred(findblock(b,1)) == Block(1)
        @test @inferred(findblockindex(b,1)) == Block(1)[1]
        @test findblock.(Ref(b),1:6) == Block.([1,2,2,3,3,3])
        @test findblockindex.(Ref(b),1:6) == BlockIndex.([1,2,2,3,3,3], [1,1,2,1,2,3])
        @test_throws BoundsError findblock(b,0)
        @test_throws BoundsError findblock(b,7)
        @test_throws BoundsError findblockindex(b,0)
        @test_throws BoundsError findblockindex(b,7)

        o = OffsetArray([2,2,3],-1:1)
        b = blockedrange(o)
        @test @inferred(findblock(b,1)) == Block(-1)
        @test @inferred(findblockindex(b,1)) == Block(-1)[1]
        @test findblock.(Ref(b),1:7) == Block.([-1,-1,0,0,1,1,1])
        @test findblockindex.(Ref(b),1:7) == BlockIndex.([-1,-1,0,0,1,1,1], [1,2,1,2,1,2,3])
        @test_throws BoundsError findblock(b,0)
        @test_throws BoundsError findblock(b,8)
        @test_throws BoundsError findblockindex(b,0)
        @test_throws BoundsError findblockindex(b,8)

        b = BlockArrays._BlockedUnitRange(-1,[-1,1,4])
        @test @inferred(findblock(b,-1)) == Block(1)
        @test @inferred(findblockindex(b,-1)) == Block(1)[1]
        @test findblock.(Ref(b),-1:4) == Block.([1,2,2,3,3,3])
        @test findblockindex.(Ref(b),-1:4) == BlockIndex.([1,2,2,3,3,3],[1,1,2,1,2,3])
        @test_throws BoundsError findblock(b,-2)
        @test_throws BoundsError findblock(b,5)
        @test_throws BoundsError findblockindex(b,-2)
        @test_throws BoundsError findblockindex(b,5)

        o = OffsetArray([2,2,3],-1:1)
        b = BlockArrays._BlockedUnitRange(-3, cumsum(o) .- 4)
        @test @inferred(findblock(b,-3)) == Block(-1)
        @test @inferred(findblockindex(b,-3)) == Block(-1)[1]
        @test findblock.(Ref(b),-3:3) == Block.([-1,-1,0,0,1,1,1])
        @test findblockindex.(Ref(b),-3:3) == BlockIndex.([-1,-1,0,0,1,1,1], [1,2,1,2,1,2,3])
        @test_throws BoundsError findblock(b,-4)
        @test_throws BoundsError findblock(b,5)
        @test_throws BoundsError findblockindex(b,-4)
        @test_throws BoundsError findblockindex(b,5)

        b = blockedrange(Fill(3,1_000_000))
        @test @inferred(findblock(b, 1)) == Block(1)
        @test @inferred(findblockindex(b, 1)) == Block(1)[1]
        @test findblock.(Ref(b),299_997:300_001) == Block.([99_999,100_000,100_000,100_000,100_001])
        @test findblockindex.(Ref(b),299_997:300_001) == BlockIndex.([99_999,100_000,100_000,100_000,100_001],[3,1,2,3,1])
        @test_throws BoundsError findblock(b,0)
        @test_throws BoundsError findblock(b,3_000_001)
        @test_throws BoundsError findblockindex(b,0)
        @test_throws BoundsError findblockindex(b,3_000_001)
    end

    @testset "BlockIndex indexing" begin
       b = blockedrange([1,2,3])
       @test b[Block(3)[2]] == b[Block(3)][2] == 5
       @test b[Block(3)[2:3]] == b[Block(3)][2:3] == 5:6
    end

    @testset "BlockRange indexing" begin
       b = blockedrange([1,2,3])
       @test b[Block.(1:2)] == blockedrange([1,2])
       @test b[Block.(1:3)] == b
       @test_throws BlockBoundsError b[Block.(0:2)]
       @test_throws BlockBoundsError b[Block.(1:4)]
    end

    @testset "misc" begin
        b = blockedrange([1,2,3])
        @test axes(b) == Base.unsafe_indices(b) == (b,)
        @test Base.dataids(b) == Base.dataids(blocklasts(b))
        @test_throws ArgumentError BlockedUnitRange(b)
    end

    @testset "OneTo interface" begin
        b = Base.OneTo(5)
        @test blockaxes(b) == (Block.(1:1),)
        @test blocksize(b) == (1,)
        @test b[Block(1)] == b
        @test b[Block(1)[2]] == 2
        @test b[Block(1)[2:3]] == 2:3
        @test blockfirsts(b) == [1]
        @test blocklasts(b) == [5]
        @test blocklengths(b) == [5]
        @test_throws BlockBoundsError b[Block(0)]
        @test_throws BlockBoundsError b[Block(2)]
        @test findblock(b,1) == Block(1)
        @test_throws BoundsError findblock(b,0)
        @test_throws BoundsError findblock(b,6)
        @test stringmime("text/plain",blockedrange([1,2,2])) == "3-blocked 5-element BlockedUnitRange{Array{Int64,1}}:\n 1\n ─\n 2\n 3\n ─\n 4\n 5"
    end
end

@testset "BlockSlice" begin
    b = BlockSlice(Block(5),1:3)
    @test b[Base.Slice(1:3)] ≡ b
    @test b[1:2] ≡ b[1:2][1:2] ≡ BlockSlice(Block(5)[1:2],1:2)
end

#=
[1,1  1,2] | [1,3  1,4  1,5]
--------------------------
[2,1  2,2] | [2,3  2,4  2,5]
[3,1  3,2] | [3,3  3,4  3,5]
----------------------------
[4,1  4,2] | [4,3  4,4  4,5]
[5,1  5,2] | [5,3  5,4  5,5]
[6,1  6,2] | [6,3  6,4  6,5]
=#

# @testset " BlockIndices" begin
#     A = BlockVector([1,2,3],[1,2])
#     @test A[Block(2)[2]] == 3
#     @test A[Block(2)[1:2]] == [2,3]
#     @test A[getindex.(Block.(1:2), 1)] == [1,2]

#     @test_throws BlockBoundsError A[Block(3)]
#     @test_throws BlockBoundsError A[Block(3)[1]]
#     @test_throws BoundsError A[Block(3)[1:1]] # this is likely an error
#     @test_throws BoundsError A[Block(2)[3]]
#     @test_throws BoundsError A[Block(2)[3:3]]
# end

@testset "sortedin" begin
    v = [1,3,4]
    @test BlockArrays.sortedin(1,v)
    @test !BlockArrays.sortedin(2,v)
    @test !BlockArrays.sortedin(0,v)
    @test !BlockArrays.sortedin(5,v)
end

@testset "eachblock" begin
    v = Array(reshape(1:6, (2, 3)))
    A = BlockArray(v, [1,1], [2,1])
    B = PseudoBlockArray(v, [1,1], [2,1])

    # test that contents match
    @test collect(eachblock(A)) == collect(eachblock(B)) == A.blocks

    # test that eachblock returns views
    first(eachblock(A))[1,2] = 0
    @test A[1,2] == 0
    first(eachblock(B))[1,2] = 0
    @test B[1,2] == 0
end


text:
import FirstOrderLp, Plots, CSV
using LaTeXStrings, DataFrames
include("../src/PrimalDualMethods.jl")

@assert length(ARGS) == 2

output_dir = ARGS[1]
method = RecoverPrimalDualMethodFromString(ARGS[2])

mkpath(output_dir)
println("Creating 1D bilinear problem\n")
""" Returns a 1D bilinear problem.
"""
simple_lp = FirstOrderLp.linear_programming_problem(
  [-Inf],  # variable_lower_bound
  [Inf],  # variable_upper_bound
  [0.0],  # objective_vector
  0.0,                 # objective_constant
  ones(1, 1),           # constraint_matrix
  [0.0],      # right_hand_side
  1,                     # num_equalities
)
restart_length = 25

# Current iterate with no restarts
params = PrimalDualOptimizerParameters(
  method, # method
  0.2,  # step_size
  1.0,  # primal_weight
  1, # record every
  20,   # printing frequency
  true, # verbose
  51,  # iteration limit
  NoRestarts(), # restart scheme
  TwoDimensionalSubspace([1.0, 0.0], [0.0, 1.0]),
  1e-10,
)

# No restarts
print("About to start optimizing\n")
solver_output = optimize(params, simple_lp, [1.0], [1.0])

no_restarts_csv_path = joinpath(output_dir, "no_restarts.csv")
CSV.write(no_restarts_csv_path, solver_output.iteration_stats)

# Fixed frequency restarts
params.restart_scheme = FixedFrequencyRestarts(restart_length, true)
print("About to start optimizing\n")
solver_output = optimize(params, simple_lp, [1.0], [1.0])

fixed_frequency_csv_path = joinpath(output_dir, "fixed_frequency_results.csv")
CSV.write(fixed_frequency_csv_path, solver_output.iteration_stats)

# Plots
fixed_frequency_df = CSV.read(fixed_frequency_csv_path, DataFrame)
no_restarts_df = CSV.read(no_restarts_csv_path, DataFrame)

Plots.plot(
  [0.0],
  [0.0],
  label = "Solution",
  markershape = :star,
  linealpha = 0.0,
  markerstrokewidth = 0.0,
  markersize = 6,
  color = "black",
  legend = :bottomleft,
)

iterates_plt = Plots.plot!(
  no_restarts_df.current_subspace_coordinate1,
  no_restarts_df.current_subspace_coordinate2,
  label = "No restarts (current iterate)",
  markershape = :diamond,
  markersize = 8,
  markerstrokewidth = 0.0,
  linealpha = 0.0,
  markerstrokealpha = 0.0,
  color = "red",
  xlabel = L"x_1",
  ylabel = L"y_1",
  markerstrokecolor = "red",
)

Plots.plot!(
  no_restarts_df.average_subspace_coordinate1[2:end],
  no_restarts_df.average_subspace_coordinate2[2:end],
  label = "No restarts (average iterate)",
  markershape = :cross,
  linealpha = 0.0,
  markerstrokewidth = 2.0,
  color = "red",
)

Plots.plot!(
  fixed_frequency_df.average_subspace_coordinate1[2:end],
  fixed_frequency_df.average_subspace_coordinate2[2:end],
  label = "Fixed frequency restarts (average iterate)",
  markershape = :xcross,
  linealpha = 0.0,
  markerstrokealpha = 0.0,
  markerstrokewidth = 2.0,
  color = "blue",
)

Plots.plot!(
  fixed_frequency_df.current_subspace_coordinate1,
  fixed_frequency_df.current_subspace_coordinate2,
  label = "Fixed frequency restarts (current iterate)",
  markershape = :circle,
  markersize = 4,
  linealpha = 0.0,
  markerstrokewidth = 0.0,
  markerstrokecolor = "blue",
  color = "blue",
)

Plots.savefig(iterates_plt, joinpath(output_dir, "bilinear-iterates.pdf"))

ext: jl

max_stars_repo_name: pedersor/google-research

max_issues_repo_name: pedersor/google-research

max_forks_repo_name: pedersor/google-research


text:
module DiferentialEvolution

using Random

export DE_classic
export DE_jither
export DE_dither
export DE_JADE

include("SupportFunctions.jl")
include("PopulationFunctions.jl")
include("MutationFunctions.jl")
include("CrossoverFunctions.jl")
include("SelectionFunctions.jl")

### Classic Diferential Evolution ###
"""
    classic_DE(fob::Function,generations::Integer,individues::Integer, lower_limits::AbstractArray, upper_limits::AbstractArray,F::Real,Cr::Real)
Also known as **DE/rand/1/bin** is the simplest differential evolution algorithm. It has a constant value of **F** and **Cr**.
"""
function DE_classic(fob::Function,generations::Integer,individues::Integer, lower_limits::AbstractArray, upper_limits::AbstractArray,F::Real,Cr::Real)
    #Error checks

    if length(lower_limits) != length(upper_limits)
        throw(ErrorException("The limits must have the same length: lower and upper limits of length $(length(lower_limits)) and $(length(upper_limits))."))
    end

    if false in [lower_limits .< upper_limits]
        throw(ErrorException("Lower limits must be less than upper limits. This happens at indexes $(collect(1:length(lower_limits))[lower_limits.>upper_limits]...)"))
    end

    if F == 0
        throw(ErrorException("F cannot be zero"))
    end

    if (Cr < 0) | (Cr > 1)
        throw(ErrorException("Cr must be in the interval [0,1]"))
    end

    #DE Algorithm

    #First generation
    init_population = initial_population(individues,lower_limits,upper_limits)
    mutated = mutation_classic(init_population,lower_limits,upper_limits,F)
    crossed = crossover_classic(init_population,mutated,Cr)
    selected = selection_classic(fob,crossed,init_population)
    #Futher generations
    for generation in 2:generations
        init_population = selected
        mutated = mutation_classic(init_population,lower_limits,upper_limits,F)
        crossed = crossover_classic(init_population,mutated,Cr)
        selected = selection_classic(fob,crossed,init_population)
    end

    return best_individue(fob,selected)
end

### Jither Diferential Evolution ###
"""
DE_jither(fob::Function,generations::Integer,individues::Integer, lower_limits::AbstractArray, upper_limits::AbstractArray,F::Real,Cr::Real)
"""
function DE_jither(fob::Function,generations::Integer,individues::Integer, lower_limits::AbstractArray, upper_limits::AbstractArray,F::Real,Cr::Real)
    #Error checks

    if length(lower_limits) != length(upper_limits)
        throw(ErrorException("The limits must have the same length: lower and upper limits of length $(length(lower_limits)) and $(length(upper_limits))."))
    end

    if false in [lower_limits .< upper_limits]
        throw(ErrorException("Lower limits must be less than upper limits. This happens at indexes $(collect(1:length(lower_limits))[lower_limits.>upper_limits]...)"))
    end

    if F == 0
        throw(ErrorException("F cannot be zero"))
    end

    if (Cr < 0) | (Cr > 1)
        throw(ErrorException("Cr must be in the interval [0,1]"))
    end

    #DE Algorithm

    #First generation
    init_population = initial_population(individues,lower_limits,upper_limits)
    mutated = mutation_jither(init_population,lower_limits,upper_limits,F)
    crossed = crossover_classic(init_population,mutated,Cr)
    selected = selection_classic(fob,crossed,init_population)
    #Futher generations
    for generation in 2:generations
        init_population = selected
        mutated = mutation_jither(init_population,lower_limits,upper_limits,F)
        crossed = crossover_classic(init_population,mutated,Cr)
        selected = selection_classic(fob,crossed,init_population)
    end

    return best_individue(fob,selected)
end

"""
DE_dither(fob::Function,generations::Integer,individues::Integer, lower_limits::AbstractArray, upper_limits::AbstractArray,F::Real,Cr::Real)
"""
function DE_dither(fob::Function,generations::Integer,individues::Integer, lower_limits::AbstractArray, upper_limits::AbstractArray,F::Real,Cr::Real)
    #Error checks

    if length(lower_limits) != length(upper_limits)
        throw(ErrorException("The limits must have the same length: lower and upper limits of length $(length(lower_limits)) and $(length(upper_limits))."))
    end

    if false in [lower_limits .< upper_limits]
        throw(ErrorException("Lower limits must be less than upper limits. This happens at indexes $(collect(1:length(lower_limits))[lower_limits.>upper_limits]...)"))
    end

    if F == 0
        throw(ErrorException("F cannot be zero"))
    end

    if (Cr < 0) | (Cr > 1)
        throw(ErrorException("Cr must be in the interval [0,1]"))
    end

    #DE Algorithm

    #First generation
    init_population = initial_population(individues,lower_limits,upper_limits)
    mutated = mutation_dither(init_population,lower_limits,upper_limits,F)
    crossed = crossover_classic(init_population,mutated,Cr)
    selected = selection_classic(fob,crossed,init_population)
    #Futher generations
    for generation in 2:generations
        init_population = selected
        mutated = mutation_dither(init_population,lower_limits,upper_limits,F)
        crossed = crossover_classic(init_population,mutated,Cr)
        selected = selection_classic(fob,crossed,init_population)
    end

    return best_individue(fob,selected)
end

function DE_JADE(fob::Function, generations::Integer, individues::Integer, lower_limits::AbstractArray, upper_limits::AbstractArray, p::Real, c::Real, μ_F=0.9, μ_Cr=0.5)

    init_population = initial_population(individues,lower_limits,upper_limits)
    number_of_parameters = length(lower_limits)

    #Initializating archive
    A = Array{Real}(undef, 0, number_of_parameters)

    for generation in 1:generations
        S_F = Vector{Real}(undef, 0)
        S_Cr = Vector{Real}(undef, 0)
        p_best = pbest(fob, init_population, p)

        for individue in 1:individues

            original_individue = init_population[individue,:]

            # Cr and F
            Cr = rand_normal_trunc(μ_Cr,0.1)
            F = rand_cauchy_trunc(μ_F,0.1)

            # Index of random individue (rnd_best)
            rnd_best = rand(1:size(p_best,1))

            # Union of population and archive
            PuA = vcat(init_population, A)
            PuA_size = size(PuA,1)

            # Random individues indexes of init_population
            rnd1 = rand_exclusive(deleteat!(collect(1:individues),individue),1)[1]
            rnd2 = rand_exclusive(deleteat!(collect(1:PuA_size),sort([individue,rnd1])),1)[1]

            # MUTATION
            mut_individue = mutation_current2pbest(original_individue, p_best[rnd_best,:], init_population[rnd1,:], PuA[rnd2,:], F)
            # CROSSOVER
            crossed_individue = crossover_JADE(original_individue, mut_individue, Cr)
            # SELECTION
            if fob(mut_individue...) <= fob(original_individue...)
                init_population[individue,:] = mut_individue
                A = vcat(A,reshape(original_individue,1,:))
                push!(S_F, F)
                push!(S_Cr, Cr)
            end
        end

        # Randomly removing solutions from A so that size(A,1) ≤ individues
        A_old_size = size(A,1)
        if A_old_size > individues
            A_new_size = rand(1:individues)
            A = A[rand_exclusive(collect(1:A_old_size), A_new_size ),:]
        end
        # Updating μ_F and μ_Cr
        μ_Cr = (1-c)*μ_Cr + c*(sum(S_Cr)/length(S_Cr))

        if isempty(S_F)
            continue
        end

        μ_F = (1-c)*μ_F + c*Lehmer_mean(S_F,2)
    end

    return best_individue(fob,init_population)
end

end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using DataStructures
using PyCall
p = PyCall.pyimport("pyplasm")
import Base.view

"""
	cuboidGrid(shape::Array{Int64,1}[, full=false])::Union{LAR,LARmodel}

Compute a *cellular complex* (mesh) with *cuboidal cells* of either `LARmodel`
or `LAR` type, depending on optional `full` parameter.

The default is
for returning a `LAR` value, i.e. a pair `(Points, Cells)`.
The *dimension* of `Cells` is the one of the number `M` of rows of
cell `Points`. The dimensions of `Array{Cells,1}` in `LARmodel` run
from ``1`` to ``M``.

# Example
```
julia> V, CV = Plasm.cuboidGrid([2,2,1])
([0.0 0.0 … 2.0 2.0; 0.0 0.0 … 2.0 2.0; 0.0 1.0 … 0.0 1.0],
Array{Int64,1}[[1,2,3,4,7,8,9,10], [3,4,5,6,9,10,11,12], [7,8,9,10,13,14,15,16],
[9,10,11,12,15,16,17,18]])

julia> V, (VV,EV,FV,CV) = Plasm.cuboidGrid([2,2,1],true);

julia> V
3×18 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0  2.0  2.0  2.0  2.0  2.0
 0.0  0.0  1.0  1.0  2.0  2.0  0.0  0.0  1.0  1.0  2.0  2.0  0.0  0.0  1.0  1.0  2.0  2.0
 0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0

julia> VV
18-element Array{Array{Int64,1},1}:
[1],[2],[3],[4],...,[12],[13],[14],[15],[16],[17],[18]

julia> EV
33-element Array{Array{Int64,1},1}:
[1,2],[3,4],[5,6],[7,8],[9,10],[11,12],,[7,13],[8,14],[9,15],[10,16],[11,17],[12,18]

julia> FV
20-element Array{Array{Int64,1},1}:
[1,2,3,4],[3,4,5,6],[7,8,9,10],[9,10,11,12],[13,14,15,16],[15,16,17,18],[1,2,7,8],[3,4,9,10],⋮,[2,4,8,10],[3,5,9,11],[4,6,10,12],[7,9,13,15],[8,10,14,16],[9,11,15,17],[10,12,16,18]

julia> CV
4-element Array{Array{Int64,1},1}:
[1,2,3,4,7,8,9,10],[3,4,5,6,9,10,11,12],[7,8,9,10,13,14,15,16],[9,10,11,12,15,16,17,18]
```
"""
#cuboidGrid = Lar.larCuboids

"""
    centroid( V::Points )::Array{Float64,1}

*Geometric center* (or *barycenter*) of a `Points` 2-array of `size` ``(M,N)``.

Each of the
``M`` coordinates of *barycenter* of the dense array of ``N`` `points` is the *mean*
of the corresponding `Points` coordinates.
"""
function centroid( V::Points )::Array{Float64,2}
	return sum(V, dims=2)/size(V,2)
end

"""
  centroid(V::Array{Float64,2})::Array{Float64,1}

*Geometric center* of a `Points` 2-array of `size` ``(M,N)``. Each of the
``M`` coordinates of *barycenter* of the dense array of ``N`` `points` is the *mean*
of the corresponding `Array` coordinates.
"""
function centroid(V::Array{Float64,2})::Array{Float64,2}
	return sum(V,2)/size(V,2)
end

"""
	cells2py(cells::Lar.Cells)::PyObject

Return a `Cells` object in a *Python* source text format. The returned `PyObject` is
 a list of lists of integers.

# Example
``` julia
julia> FV = Lar.cuboid([1,1,1],true)[2][3]
6-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
 [5, 6, 7, 8]
 [1, 2, 5, 6]
 [3, 4, 7, 8]
 [1, 3, 5, 7]
 [2, 4, 6, 8]

julia> Plasm.cells2py(FV)
PyObject [[1, 2, 3, 4], [5, 6, 7, 8], [1, 2, 5, 6], [3, 4, 7, 8],
[1, 3, 5, 7], [2, 4, 6, 8]]
```
"""
function cells2py(cells::Lar.Cells)::PyObject
	return PyObject([Any[cell[h] for h=1:length(cell)] for cell in cells])
end

"""
	points2py(V::Points)::PyObject

Return a `Points` object in a *Python* source text format. The returned `PyObject` is
 a `list of lists of float`.

# Example
``` julia
julia> V = Lar.cuboid([1,1,1])[1]
3×8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  0.0  0.0  1.0  1.0
 0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0

julia> Plasm.points2py(V)
PyObject [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 1.0, 0.0], [0.0, 1.0, 1.0],
[1.0, 0.0, 0.0], [1.0, 0.0, 1.0], [1.0, 1.0, 0.0], [1.0, 1.0, 1.0]]
```
"""
function points2py(V::Points)::PyObject
	return PyObject([Any[V[h,k] for h=1:size(V,1)] for k=1:size(V,2)])
end

"""
	mkpol(verts::Points, cells::Cells)::Hpc

Return an `Hpc` object starting from a `Points` and a `Cells` object. *HPC =
Hierarchical Polyhedral Complex* is the geometric data structure
used by `PLaSM` (Programming LAnguage for Solid Modeling). See the Wiley's book
[*Geometric Programming for Computer-Aided Design*]
(https://onlinelibrary.wiley.com/doi/book/10.1002/0470013885) and its
current `Python` library [*https://github.com/plasm-language/pyplasm*]
(https://github.com/plasm-language/pyplasm).

```julia
julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

julia> hpc = Plasm.mkpol(V,EV)
PyObject <pyplasm.xgepy.Hpc; proxy of <Swig Object of type
'std::shared_ptr< Hpc > *' at 0x12cf45d50> >
	```
"""
function mkpol( verts::Plasm.Points, cells::Plasm.Cells )::Plasm.Hpc
	p = PyCall.pyimport("pyplasm")
	Verts = Plasm.points2py(verts)
	Cells = Plasm.cells2py(cells)
	hpc = p["MKPOL"]([Verts,Cells,PyObject(true)])
	return hpc
end
function mkpol( hpcs::Array{Plasm.Hpc,1} )::Plasm.Hpc
	p = PyCall.pyimport("pyplasm")
	return p["MKPOL"](PyCall.PyVector(hpcs))
end

"""
	view(hpc::PyCall.PyObject)

Base.view extension.
Display a *Python*  `HPC` (Hierarchical Polyhedral Complex) `object` using
the *`PyPlasm` viewer*, written in C++ with `OpenGL` and acceleration algorithms
for *big geometric data* structures.

# Example
``` julia
julia> m = Lar.cuboidGrid([2,2],true)
([0.0 0.0 … 2.0 2.0; 0.0 1.0 … 1.0 2.0], Array{Array{Int64,1},1}[Array{Int64,1}[[1],
[2], [3], [4], [5], [6], [7], [8], [9]], Array{Int64,1}[[1, 2], [2, 3], [4, 5], [5,
6], [7, 8], [8, 9], [1, 4], [2, 5], [3, 6], [4, 7], [5, 8], [6, 9]],
Array{Int64,1}[[1, 2, 4, 5], [2, 3, 5, 6], [4, 5, 7, 8], [5, 6, 8, 9]]])

julia> hpc = Plasm.mkpol(m[1],m[2][2])
PyObject <pyplasm.xgepy.Hpc; proxy of <Swig Object of type 'std::shared_ptr< Hpc > *'
at 0x140d6c780> >

julia> Plasm.view(hpc)
```
"""
function view(hpc)
	p = PyCall.pyimport("pyplasm")
	p["VIEW"](hpc)
end

"""
	view(hpcs::Array{PyCall.PyObject})

Base.view extension.
Display a *Python*  `HPC` (Hierarchical Polyhedral Complex) `object`, starting from
a *Julia* array of `PyCall.PyObject`s.

# Example

```
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

EV = [[1,2],[2,3],[3,4],[4,5],[5,6],[1,6],
[6,7],[4,7],[2,7],[5,8],[1,8],[3,8]]
V = hcat([[2.,1],[2,2],[1,2],[0,3],[0,0],[3,0],[3,3],[1,1]]...)

spanningtree, fronds = Lar.depth_first_search(EV)

hpc1 = Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)], spanningtree]))
hpc2 = Plasm.numbering(1.25)((V,[[[k] for k=1:size(V,2)], fronds]))
Plasm.view([ Plasm.color("cyan")(hpc1) , Plasm.color("magenta")(hpc2) ])
```
"""
function view(hpcs::Array{PyCall.PyObject})
	p = PyCall.pyimport("pyplasm")
	hpc = p["STRUCT"](hpcs)
	p["VIEW"](hpc)
end

"""
	view(V::Points, CV::Cells)

Base.view extension.
Display a *Python*  `HPC` (Hierarchica Polyhedral Complex) `object` using
the *`PyPlasm` viewer*, written in C++ with `OpenGL` and acceleration algorithms
for *big geometric data* structures. Input parameters are of `Points` and `Cells`
type.

# Example

```julia
julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

julia> Plasm.mkpol(V,CV)
PyObject <pyplasm.xgepy.Hpc; proxy of <Swig Object of type
'std::shared_ptr< Hpc > *' at 0x12cf45d50> >

julia>
Plasm.view(Plasm.mkpol(V,CV))
```
"""
function view(V::Points, CV::Cells)
	p = PyCall.pyimport("pyplasm")
	hpc = lar2hpc(V::Points, CV::Cells)
	p["VIEW"](hpc)
end

"""
  build_K(FV::Lar.Cells)::ChainOp

The *characteristic matrix* of type `ChainOp` from 1-cells (edges) to 0-cells (vertices)
"""
function build_K(FV::Lar.Cells)
	   I = Int64[]; J = Int64[]; V = Int64[]
	   for (i,face) in enumerate(FV)
			   for v in face
					   push!(I,i)
					   push!(J,v)
					   push!(V,1)
			   end
	   end
	   kEV = SparseArrays.sparse(I,J,V)
end

"""
	view(V::Lar.Points, CV::Lar.ChainOp)

# Example

```julia
julia> V, (VV,EV,FV,CV) = Plasm.cuboidGrid([10,10,1],true);

julia> copCV = convert(Plasm.ChainOp,Plasm.build_K(FV));

julia> view(V,copCV)

```
"""
function view(V::Points, copCV::ChainOp)
	CV = [findnz(copCV[k,:])[1] for k=1:size(copCV,1)]
	view(V,CV)
end

"""
	view(model::LARmodel)

Base.view extension.
Display a *Python*  `HPC` (Hierarchica Polyhedral Complex) `object` using
the *`PyPlasm` viewer*, written in C++ with `OpenGL` and acceleration algorithms
for *big geometric data* structures. The input is a `LARmodel` object.

# Example

```julia
julia> typeof( Lar.cuboid([1,1,1], true) )
Tuple{Array{Float64,2},Array{Array{Int64,1},1}}

julia> V,(VV,EV,FV,CV) = Lar.cuboid([.5,.5,.5], true);

julia> Plasm.view( (V,[VV,EV,FV,CV]) )
```
"""
function view(model::LARmodel)
	p = PyCall.pyimport("pyplasm")
	hpc = hpc_exploded(model::LARmodel)(1.2,1.2,1.2)
	p["VIEW"](hpc)
end

"""
	view(pair::Tuple{Points,Cells})

Base.view extension.
Display a *Python*  `HPC` (Hierarchica Polyhedral Complex) `object` using
the *`PyPlasm` viewer*, written in C++ with `OpenGL` and acceleration algorithms
for *big geometric data* structures. The input is a `pair` of type
`Tuple{Points,Cells}`.

# Example
```julia
julia> typeof(Lar.cuboid([1,1,1])::LAR)
Tuple{Array{Float64,2},Array{Array{Int64,1},1}}

julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1], true);

julia> Plasm.view( (V,FV) );
```
"""
function view(pair::Tuple{Points,Cells})
	V,CV = pair
	p = PyCall.pyimport("pyplasm")
	hpc = lar2hpc(V::Points, CV::Cells)
	p["VIEW"](hpc)
end

"""
	view(obj::Lar.Struct)

Display a geometric value of `Struct` type, via conversion to `LAR`
and then to `Hpc` values.

# Example
```julia
cube = Lar.apply( Lar.t(-.5,-.5,0), Lar.cuboid([1,1,1]));
tableTop = Lar.Struct([ Lar.t(0,0,.85), Lar.s(1,1,.05), cube ]);
tableLeg = Lar.Struct([ Lar.t(-.475,-.475,0), Lar.s(.1,.1,.89), cube ]);
tablelegs = Lar.Struct( repeat([ tableLeg, Lar.r(0,0,pi/2) ],outer=4) );
table = Lar.Struct([ tableTop, tablelegs ]);

Plasm.view(table)
```
"""
function view(obj::Lar.Struct)
	lar = Lar.struct2lar(obj)
	view(lar)
end

"""
	view(scene::Array{Any,1})

Display a geometric `scene`, defined as `Array{Any,1}` of geometric objects
defined in the *same* coordinate system, i.e. in *World Coordinates*.

# Example

A hierarchical `scene` defined in *Local Coordinates* as value of `Struct` type,
must be converted to `Array{Any,1}` by the expression

	`evalStruct(scene::Struct)::Array{Any,1}`

```julia
cube = Lar.apply( Lar.t(-.5,-.5,0), Lar.cuboid([1,1,1]));
tableTop = Lar.Struct([ Lar.t(0,0,.85), Lar.s(1,1,.05), cube ]);
tableLeg = Lar.Struct([ Lar.t(-.475,-.475,0), Lar.s(.1,.1,.89), cube ]);
tablelegs = Lar.Struct( repeat([ tableLeg, Lar.r(0,0,pi/2) ],outer=4) );
table = Lar.Struct([ tableTop, tablelegs ]);

scene = Lar.evalStruct(table);
# output
# 5-element Array{Any,1}

Plasm.view(scene)
```
"""
function view(scene::Array{Any,1})
	if prod([isa(item[1:2],Lar.LAR) for item in scene])
		p = PyCall.pyimport("pyplasm")
		p["VIEW"](p["STRUCT"]([Plasm.lar2hpc(item[1],item[2]) for item in scene]))
	end
end

"""
	hpc_exploded( model::LARmodel )( sx=1.2, sy=1.2, sz=1.2 )::Hpc

Convert a `LARmodel` into a `Hpc` object, after exploding all-dimensional cells with
scale `sx,sy,sz` parameters. Every cell is *translated* by the vector difference
between its *scaled centroid* and its *centroid*. Every cell is transformed in a
single `LAR` object before explosion.

# Example
```julia
julia> hpc = Plasm.hpc_exploded(
	Lar.cuboidGrid([3,2,1], true))(1.5,1.5,1.5)

julia> view(hpc)
```
"""
function hpc_exploded( model )
		function hpc_exploded0( sx=1.2, sy=1.2, sz=1.2 )
		p = PyCall.pyimport("pyplasm")
		verts,cells = model
		out = []
		for skeleton in cells
			for cell in skeleton
				vcell = hcat([[verts[h,k] for h=1:size(verts,1)] for k in cell]...)

				center = sum([verts[:,v] for v in cell])/length(cell)
				scaled_center = length(center)==2 ? center.*[sx,sy] :
													center.*[sx,sy,sz]
				translation_vector = scaled_center-center
				vcell = vcell .+ translation_vector

				py_verts = Plasm.points2py(vcell)
				py_cells = Plasm.cells2py( [collect(1:size(vcell,2))] )

				hpc = p["MKPOL"]([ py_verts, py_cells, [] ])
				push!(out, hpc)
			end
		end
		hpc = p["STRUCT"](out)
		return hpc
	end
	return hpc_exploded0
end

"""
	lar_exploded(model)

Generate an exploded `Hpc` value starting from a `LARmodel` instance.
Useful to show the explosion of a 2-complex made by possibly non-convex polygons.

# Example

```
V, copEV, copFE = Lar.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)
EVs = FV2EVs(copEV, copFE) # polygonal face fragments

Plasm.viewcolor(V::Lar.Points, EVs::Array{Lar.Cells})
model = V,EVs
Plasm.view(lar_exploded(model)(1.2,1.2,1.2))
```
"""
function lar_exploded(model)
	function lar_exploded0( sx=1.2, sy=1.2, sz=1.2 )
		p = PyCall.pyimport("pyplasm")
		verts,arrayofcells = model
		out = []
		for cell in arrayofcells
			vcell = convert(Lar.Cells,[collect(Set(cat(cell)))])

			center = sum([verts[:,v] for v in vcell[1]])/length(vcell[1])
			scaled_center = size(center,1)==2 ? center .* [sx,sy] :
												center .* [sx,sy,sz]
			translation_vector = scaled_center - center
			vertcell = [verts[:,k]+translation_vector for k in vcell[1]]
			cellverts = hcat(vertcell...)

			py_verts = Plasm.points2py( cellverts )
			vdict = DataStructures.OrderedDict( zip( vcell[1],
						[k for k=1:length(vcell[1])] ))
			edges = [[vdict[e[1]], vdict[e[2]]] for e in cell]
			py_cells = Plasm.cells2py( edges )

			hpc = p["MKPOL"]([ py_verts, py_cells, [] ])
			push!(out, hpc)
		end
		hpc = p["STRUCT"](out)
		return hpc
	end
	return lar_exploded0
end

"""
	lar2hpc(V::Points, CV::Cells)::Hpc

Return an `Hpc` object starting from a `Points` and a `Cells` object. *HPC =
Hierarchical Polyhedral Complex* is the geometric data structure
used by `PLaSM` (Programming LAnguage for Solid Modeling). See the Wiley's book
[*Geometric Programming for Computer-Aided Design*]
(https://onlinelibrary.wiley.com/doi/book/10.1002/0470013885) and its
current `Python` library [*https://github.com/plasm-language/pyplasm*]
(https://github.com/plasm-language/pyplasm).
# Example

```julia
julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

julia> hpc = lar2hpc( (V, CV)::Plasm.LAR ... )::Hpc
PyObject <pyplasm.xgepy.Hpc; proxy of <Swig Object of type
'std::shared_ptr< Hpc > *' at 0x12cf45d50> >

julia> view(hpc)
```
"""
function lar2hpc(V::Points, CV::Cells)::Hpc
	hpc = Plasm.mkpol(V,CV)
	return hpc
end

"""
	lar2hpc(model::LARmodel)::Hpc

Return an `Hpc` object starting from a `LARmodel` object. *HPC =
Hierarchical Polyhedral Complex* is the geometric data structure
used by `PLaSM` (Programming LAnguage for Solid Modeling). See the Wiley's book
[*Geometric Programming for Computer-Aided Design*]
(https://onlinelibrary.wiley.com/doi/book/10.1002/0470013885) and its
current `Python` library [*https://github.com/plasm-language/pyplasm*]
(https://github.com/plasm-language/pyplasm).

# Example
```julia
julia> model = Lar.cuboid([1,1,1],true);

julia> view( Plasm.lar2hpc(model) )

julia> view( Plasm.hpc_exploded(model)(1.5,1.5,1.5) )
```
"""
function lar2hpc(model::LARmodel)::Hpc
	verts = model[1]
	cells = Array{Int,1}[]
	for item in model[2]
		append!(cells, item)
	end
	hpc = mkpol(verts,cells)
end

"""
	lar2hpc(scene::Array{Any,1})::Hpc

# Example
```julia
square = ([[0; 0] [0; 1] [1; 0] [1; 1]], [[1, 2, 3, 4]],
[[1,2], [1,3], [2,4], [3,4]])
V,FV,EV  =  square
model  =  V,([[1],[2],[3],[4]],EV,FV)
table = Lar.apply(Lar.t(-0.5,-0.5), square)
chair = Lar.Struct([Lar.t(0.75,0),Lar.s(0.35,0.35),table])
structo = Lar.Struct([Lar.t(2,1),table,repeat([Lar.r(pi/2),chair],
		outer = 4)...])
structo1 = Lar.Struct(repeat([structo,Lar.t(0,2.5)],outer = 10));
structo2 = Lar.Struct(repeat([structo1,Lar.t(3,0)],outer = 10));
scene = Lar.evalStruct(structo2);

Plasm.view(Plasm.lar2hpc(scene))
```

"""
function lar2hpc(scene::Array{Any,1})::Hpc
	p = PyCall.pyimport("pyplasm")
	hpc = p["STRUCT"]([ mkpol(item[1],item[2]) for item in scene ])
end

"""
	lar2exploded_hpc(V::Lar.Points,CV::Lar.Cells)::Hpc

Input `V::Points` and `CV::Cells`. Output an *exploded* `Hpc`
object,  exploding cells in `CV` with scale `sx,sy,sz` parameters.

Every cell is *translated* by the vector difference
between its *scaled centroid* and its *centroid*.

*HPC = Hierarchical Polyhedral Complex* is the geometric data structure
used by `PLaSM` (Programming LAnguage for Solid Modeling).

#	Example

```
julia> V,cells = Lar.cuboidGrid([3,3,1], true)

julia> hpc = Plasm.lar2exploded_hpc(V::Lar.Points, cells[4]::Lar.Cells)()

julia> Plasm.view(hpc)
```
"""
function lar2exploded_hpc(V::Lar.Points, cells::Lar.Cells)
	function lar2exploded_hpc0(sx=1.2, sy=1.2, sz=1.2)
		hpc = Plasm.hpc_exploded( (V,[cells]) )(sx,sy,sz)
	end
	return lar2exploded_hpc0
end

function viewexploded(V::Lar.Points, cells::Lar.Cells)
	function lar2exploded_hpc0(sx=1.2, sy=1.2, sz=1.2)
		hpc = Plasm.hpc_exploded( (V,[cells]) )(sx,sy,sz)
		Plasm.view(hpc)
	end
	return lar2exploded_hpc0
end

function viewexploded(W::Lar.Points, cells::Lar.ChainOp) # W by rows ...
	V = convert(Lar.Points, transpose(W))
	EV = [findnz(cells[k,:])[1] for k=1:size(cells,1)]
	function lar2exploded_hpc0(sx=1.2, sy=1.2, sz=1.2)
		Plasm.viewexploded(V::Lar.Points, EV::Lar.Cells)(sx,sy,sz)
	end
	return lar2exploded_hpc0
end

"""
	viewcolor(V::Lar.Points, CVs::Array{Lar.Cells})

Display a colored view of the function parameters, drawing each `Cells` array in function parameters with a random color.

# Example

```

```
"""
function viewcolor(V::Lar.Points, CVs::Array{Lar.Cells})
	hpcs = [ Plasm.lar2hpc(V,CVs[i]) for i=1:length(CVs) ]
	Plasm.view([ Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k])
		for k=1:(length(hpcs)) ])
end

function viewcolorexploded(V::Lar.Points, CVs::Array{Lar.Cells})
	n = length(CVs)
	Vs = Array{Array{Float64,2},1}(undef, n)

	for (k,CV) in enumerate(CVs)
		vs = sort(collect(Set(cat(CV))))
		vdict = OrderedDict(zip(vs,1:length(vs)))
		CV = [[vdict[h] for h in cell] for cell in CV]
		Vs[k] = hcat([V[:,i] for i in vs]...)
	end

	hpcs = [ Plasm.lar2hpc(Vs[i],CVs[i]) for i=1:n ]
	Plasm.view([ Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k])
		for k=1:(length(hpcs)) ])
end

"""
	cubicbezier2D(curvePts::Array{Array{Float,1},1})

Produce the two ``coordinate functions`` for a ``cubic`` *Bézier curve* in ``2D``.
The input is the array `curvePts` of 4 control points with 2 coordinates.

# Example

```
julia> curvePts = [[1232.24, 1340.84],[1259.53, 1119.86],
	[1417.91, 1015.16],[1133.47, 1010.63]]

julia> bx,by = cubicbezier2D(curvePts);

julia> (bx(0), by(0))
(1232.24, 1340.84)

julia> (bx(1), by(1))
(1133.47, 1010.63)
```
"""
function cubicbezier2D( curvePts::Array{Array{Float64,1},1} )
	b1(u) = (1 - u)^3
	b2(u) = 3*u*(1 - u)^2
	b3(u) = 3*u^2*(1 - u)
	b4(u) = u^3
	cntrlverts = hcat(curvePts...)
	x = cntrlverts[1,:]
	y = cntrlverts[2,:]
	Bx = u -> x[1]*b1(u) + x[2]*b2(u) + x[3]*b3(u) + x[4]*b4(u)
	By = u -> y[1]*b1(u) + y[2]*b2(u) + y[3]*b3(u) + y[4]*b4(u)
	return Bx,By
end

"""
	const cmdsplit

Regex for splitting the `d` element of a `<path` element into
a sequence of graphics commnds.
"""
const cmdsplit = r"\s*([mMzZlLhHvVcCsSqQtTaA])\s*"

"""
	const digitRegEx

Regex for extracting numeric params from from a sequence of graphics
commnds of a `<path` element.
"""
const digitRegEx = r"\s*(-?[0-9]*\.?\d+)\s*"

"""
	pathparse(elements)

Parse graphics commands in SVG `<path` tagged element.

# Example

```julia

```
"""
function pathparse(data, cnrtlpolygon=false)
	tokens = split(data,'\"')
	pathstring = tokens[2]

	iterator = eachmatch(cmdsplit, pathstring)
	offsets = [m.offset for m in iterator]
	commands = [m.match for m in iterator]

	substrings = [pathstring[offsets[k]+1:offsets[k+1]-1]
		for k=1:length(offsets)-1]
	push!(substrings, pathstring[offsets[end]+1:end])

	params = Array{Float64,1}[]
	for substring in substrings
		numbers = []
		for match in eachmatch(digitRegEx, substring)
			push!(numbers,String(match.match))
		end
		numbers = map(x->parse(Float64,x),numbers)
		push!(params,numbers)
	end

	lines = Array{Float64,1}[]
	global startpoint, endpoint = 0.0, 0.0
	for (command, args) in zip(commands, params)
		if command == "M"
			global startpoint = args
		elseif command == "L"
			endpoint = args
			line = vcat([startpoint,endpoint]...)
			startpoint = endpoint
			push!(lines, line)
		elseif command == "C"
			contrlpt1 = args[1:2]
			contrlpt2 = args[3:4]
			endpoint = args[5:6]
			curvePts = [startpoint,contrlpt1,contrlpt2,endpoint]
			println(curvePts)

			Bx,By = cubicbezier2D(curvePts)
			pts = [[Bx(u),By(u)] for u=0:.1:1]
			curvelines = [ vcat([pts[k],pts[k+1]]...) for k=1:length(pts)-1 ]
			for line in curvelines
				push!(lines, line)
			end

 			if cnrtlpolygon
				push!(lines,vcat(curvePts[1:2]...),vcat(curvePts[2:3]...),
					vcat(curvePts[3:4]...))
			end
			startpoint = endpoint
		end

	end
	return lines
end

"""
	lines2lar(lines)

LAR model construction from array of float quadruples.
Each `line` in input array stands for `x1,y1,x2,y2`.
"""
function lines2lar(lines)
	vertdict = OrderedDict{Array{Float64,1}, Int64}()
	EV = Array{Int64,1}[]
	idx = 0
	for h=1:size(lines,2)
		x1,y1,x2,y2 = lines[:,h]

		if ! haskey(vertdict, [x1,y1])
			idx += 1
			vertdict[[x1,y1]] = idx
		end
		if ! haskey(vertdict, [x2,y2])
			idx += 1
			vertdict[[x2,y2]] = idx
		end
		v1,v2 = vertdict[[x1,y1]],vertdict[[x2,y2]]
		push!(EV, [v1,v2])
	end
	V = hcat(collect(keys(vertdict))...)
	return V,EV
end

"""
	normalize(V::Lar.Points; flag=true::Bool)::Lar.Points

2D normalization transformation (isomorphic by defaults) of model
vertices to normalized coordinates ``[0,1]^2``. Used with SVG importing.
"""
function normalize(V::Lar.Points; flag=true)
	m,n = size(V)
	if m > n # V by rows
		V = convert(Lar.Points, V')
	end

	xmin = minimum(V[1,:]); ymin = minimum(V[2,:]);
	xmax = maximum(V[1,:]); ymax = maximum(V[2,:]);
	box = [[xmin; ymin] [xmax; ymax]]	# containment box
	aspectratio = (xmax-xmin)/(ymax-ymin)
	if flag
		if aspectratio > 1
			umin = 0; umax = 1
			vmin = 0; vmax = 1/aspectratio; ty = vmax
		elseif aspectratio < 1
			umin = 0; umax = aspectratio
			vmin = 0; vmax = 1; ty = vmax
		end
		T = Lar.t(0,ty) * Lar.s(1,-1) * Lar.s((umax-umin), (vmax-vmin)) *
			Lar.s(1/(xmax-xmin),1/(ymax-ymin)) * Lar.t(-xmin,-ymin)
	else
		T = Lar.t(0, ymax-ymin) * Lar.s(1,-1)
	end
	dim = size(V,1)
	W = T[1:dim,:] * [V;ones(1,size(V,2))]
	#V = map( x->round(x,digits=8), W )
	V = map(Lar.approxVal(8), W)

	if m > n # V by rows
		V = convert(Lar.Points, V')
	end
	return V
end

function normalize3D(V::Lar.Points; flag=true)
	Vxy = Plasm.normalize(V[1:2,:], flag=flag)
	zmin = minimum(V[3,:]); zmax = maximum(V[3,:]);
	Vz = (V[3,:] .- zmin) / (zmax - zmin)
	return [Vxy; Vz']
end

"""
	svg2lar(filename::String; flag=true)::Lar.LAR

Parse a SVG file to a `LAR` model `(V,EV)`.
Only  `<line >` and `<rect >` and `<path >` SVG primitives are currently translated.
TODO:  interpretation of transformations.
"""
function svg2lar(filename::String; flag=true)::Lar.LAR
	outlines = Array{Float64,1}[]
	for line in eachline(filename)
		parts = split(line, ' ')
		elements = [part for part in parts if part≠""]
		tag = elements[1]

		# SVG <line > primitives
		if tag == "<line"
			regex = r"""(<line )(.+)(" x1=")(.+)(" y1=")(.+)(" x2=")(.+)(" y2=")(.+)("/>)"""
			coords = collect(match( regex , line)[k] for k in (4,6,8,10))
			outline = [ parse(Float64, string) for string in coords ]
			push!(outlines, outline)
		# SVG <rect > primitives
		elseif tag == "<rect"
			regex = r"""(<rect x=")(.+?)(" y=")(.+?)(" )(.*?)( width=")(.+?)(" height=")(.+?)("/>)"""
			coords = collect(match( regex , line)[k] for k in (2,4,8,10))
			x, y, width, height = [ parse(Float64, string) for string in coords ]
			line1 = [ x, y, x+width, y ]
			line2 = [ x, y, x, y+height ]
			line3 = [ x+width, y, x+width, y+height ]
			line4 = [ x, y+height, x+width, y+height ]
			push!(outlines, line1, line2, line3, line4)
		# SVG <path  > primitives
		# see https://github.com/chebfun/chebfun/issues/1617
		elseif tag == "<path"
			dataregex = r"""d=(".*?")(.*?)"""
			data = string(match( dataregex , line)[1])
			polyline = pathparse(data)
			append!(outlines, polyline)
		end
	end
	lines = hcat(outlines...)
	lines = map( x->round(x,sigdigits=8), lines )
	# LAR model construction
	V,EV = lines2lar(lines)
	# normalization
	V = normalize(V,flag=flag)
	return V,EV
end


text:
"""
    DiscreteNonParametricSampler(xs, ps)

Data structure for efficiently sampling from an arbitrary probability mass
function defined by support `xs` and probabilities `ps`.
"""
struct DiscreteNonParametricSampler{T<:Real, S<:AbstractVector{T}, A<:AliasTable} <: Sampleable{Univariate,Discrete}
    support::S
    aliastable::A

    function DiscreteNonParametricSampler{T,S}(support::S, probs::AbstractVector{<:Real}
    ) where {T<:Real,S<:AbstractVector{T}}
        aliastable = AliasTable(probs)
        new{T,S,typeof(aliastable)}(support, aliastable)
    end
end

DiscreteNonParametricSampler(support::S, probs::AbstractVector{<:Real}
    ) where {T<:Real,S<:AbstractVector{T}} =
    DiscreteNonParametricSampler{T,S}(support, probs)

rand(rng::AbstractRNG, s::DiscreteNonParametricSampler) =
    (@inbounds v = s.support[rand(rng, s.aliastable)]; v)


text:
"""
`iscolumn`
"""
function iscolumn(df::DataFrames.AbstractDataFrame,what::Symbol)
    if in(what,names(df))
        return true
    else
        return false
    end
end

"""
`count_sequence`
count when along an array the value of i+1 is different from i
"""
function count_sequence(series::AbstractArray)
    v_count = [1]
    count = 1
    for i = 2:size(series,1)
        if series[i] == series[i-1]
            count = count
        else
            count = count+1
        end
        push!(v_count,count)
    end
    return v_count
end

"""
"count_series"
useful to count an ordered series of categorical values, wheter they are the same or
they change
"""
function count_trues(count::T, same) where {T <: Number}
    if same
        count+1
    else
        1
    end
end

function check_series(x)
    Vector = [false]
    for i = 2:size(x,1)
        push!(Vector,x[i] == x[i-1] ? true : false)
    end
    return Vector
end

function count_series(vector)
    x = check_series(vector)
    res = accumulate(count_trues,x;init=0)
    return res
end

"""
"checktype"
check if a column supposed to be boolean is a string and convert it
"""

function checktype(v::AbstractArray)
    eltype(v) == Bool ? v : occursin(r"ue","True")
end

function checktype(t::IndexedTables.AbstractIndexedTable,booleans::AbstractArray)
    for x in booleans
        println(x)
        t = setcol(t,x,checktype((column(t,x))))
    end
    return t
end

"""
"by_summary"
retrieve the first value of a column from a table split by by
"""
function by_summary(df::IndexedTables.IndexedTable,by,x::Symbol)
    # JuliaDB.groupby(NamedTuple{(x,)}(t-> getindex(t,1),),df,by,select = x)
    # JuliaDBMeta.@groupby df by {MouseID = cols(x)[1]}
    t = JuliaDB.groupby(t-> getindex(t,1),df,by,select = x)
    JuliaDB.rename(t,(colnames(t)[end]=>x))
end

function by_summary(df::IndexedTables.IndexedTable,by,x::AbstractArray{Symbol})
    summary_table = table((a=[1],b=[2]))
    for (i,c) in enumerate(x)
        if c in colnames(df)
            if i == 1
                summary_table = by_summary(df,by,c)
            else
                ongoing = by_summary(df,by,c)
                summary_table = join(summary_table, ongoing, lkey = by, rkey = by)
            end
        end
    end
    return summary_table
end


text:
"""
"""
abstract type SparseMatrixAssembler <: Assembler end

"""
"""
function get_matrix_builder(a::SparseMatrixAssembler)
  @abstractmethod
end

"""
"""
function get_vector_builder(a::SparseMatrixAssembler)
  @abstractmethod
end

"""
"""
function symbolic_loop_matrix!(A,a::SparseMatrixAssembler,matdata)
  @abstractmethod
end

"""
"""
function symbolic_loop_vector!(b,a::SparseMatrixAssembler,vecdata)
  @abstractmethod
end

"""
"""
function symbolic_loop_matrix_and_vector!(A,b,a::SparseMatrixAssembler,data)
  @abstractmethod
end

"""
"""
function numeric_loop_matrix!(A,a::SparseMatrixAssembler,matdata)
  @abstractmethod
end

"""
"""
function numeric_loop_vector!(b,a::SparseMatrixAssembler,vecdata)
  @abstractmethod
end

"""
"""
function numeric_loop_matrix_and_vector!(A,b,a::SparseMatrixAssembler,data)
  @abstractmethod
end

get_matrix_type(a::SparseMatrixAssembler) = get_array_type(get_matrix_builder(a))
get_vector_type(a::SparseMatrixAssembler) = get_array_type(get_vector_builder(a))

function allocate_vector(a::SparseMatrixAssembler,vecdata)
  v1 = nz_counter(get_vector_builder(a),(get_rows(a),))
  symbolic_loop_vector!(v1,a,vecdata)
  v2 = nz_allocation(v1)
  symbolic_loop_vector!(v2,a,vecdata)
  v3 = create_from_nz(v2)
  v3
end

function assemble_vector!(b,a::SparseMatrixAssembler,vecdata)
  fill_entries!(b,zero(eltype(b)))
  assemble_vector_add!(b,a,vecdata)
end

function assemble_vector_add!(b,a::SparseMatrixAssembler,vecdata)
  numeric_loop_vector!(b,a,vecdata)
  create_from_nz(b)
end

function assemble_vector(a::SparseMatrixAssembler,vecdata)
  v1 = nz_counter(get_vector_builder(a),(get_rows(a),))
  symbolic_loop_vector!(v1,a,vecdata)
  v2 = nz_allocation(v1)
  numeric_loop_vector!(v2,a,vecdata)
  v3 = create_from_nz(v2)
  v3
end

function allocate_matrix(a::SparseMatrixAssembler,matdata)
  m1 = nz_counter(get_matrix_builder(a),(get_rows(a),get_cols(a)))
  symbolic_loop_matrix!(m1,a,matdata)
  m2 = nz_allocation(m1)
  symbolic_loop_matrix!(m2,a,matdata)
  m3 = create_from_nz(m2)
  m3
end

function assemble_matrix!(mat,a::SparseMatrixAssembler,matdata)
  fill_entries!(mat,zero(eltype(mat)))
  assemble_matrix_add!(mat,a,matdata)
end

function assemble_matrix_add!(mat,a::SparseMatrixAssembler,matdata)
  numeric_loop_matrix!(mat,a,matdata)
  create_from_nz(mat)
end

function assemble_matrix(a::SparseMatrixAssembler,matdata)
  m1 = nz_counter(get_matrix_builder(a),(get_rows(a),get_cols(a)))
  symbolic_loop_matrix!(m1,a,matdata)
  m2 = nz_allocation(m1)
  numeric_loop_matrix!(m2,a,matdata)
  m3 = create_from_nz(m2)
  m3
end

function allocate_matrix_and_vector(a::SparseMatrixAssembler,data)
  m1 = nz_counter(get_matrix_builder(a),(get_rows(a),get_cols(a)))
  v1 = nz_counter(get_vector_builder(a),(get_rows(a),))
  symbolic_loop_matrix_and_vector!(m1,v1,a,data)
  m2 = nz_allocation(m1)
  v2 = nz_allocation(v1)
  symbolic_loop_matrix_and_vector!(m2,v2,a,data)
  m3 = create_from_nz(m2)
  v3 = create_from_nz(v2)
  m3,v3
end

function assemble_matrix_and_vector!(A,b,a::SparseMatrixAssembler, data)
  fill_entries!(A,zero(eltype(A)))
  fill_entries!(b,zero(eltype(b)))
  assemble_matrix_and_vector_add!(A,b,a,data)
end

function assemble_matrix_and_vector_add!(A,b,a::SparseMatrixAssembler,data)
  numeric_loop_matrix_and_vector!(A,b,a,data)
  create_from_nz(A), create_from_nz(b)
end

function assemble_matrix_and_vector(a::SparseMatrixAssembler, data)
  m1 = nz_counter(get_matrix_builder(a),(get_rows(a),get_cols(a)))
  v1 = nz_counter(get_vector_builder(a),(get_rows(a),))
  symbolic_loop_matrix_and_vector!(m1,v1,a,data)
  m2 = nz_allocation(m1)
  v2 = nz_allocation(v1)
  numeric_loop_matrix_and_vector!(m2,v2,a,data)
  m3 = create_from_nz(m2)
  v3 = create_from_nz(v2)
  m3,v3
end

function test_sparse_matrix_assembler(a::SparseMatrixAssembler,matdata,vecdata,data)
  test_assembler(a,matdata,vecdata,data)
  _ = get_matrix_builder(a)
  _ = get_vector_builder(a)
end

struct GenericSparseMatrixAssembler <: SparseMatrixAssembler
  matrix_builder
  vector_builder
  rows::AbstractUnitRange
  cols::AbstractUnitRange
  strategy::AssemblyStrategy
end

function SparseMatrixAssembler(
  mat,
  vec,
  trial::FESpace,
  test::FESpace,
  strategy::AssemblyStrategy=DefaultAssemblyStrategy())

  rows = get_free_dof_ids(test)
  cols = get_free_dof_ids(trial)
  GenericSparseMatrixAssembler(
    SparseMatrixBuilder(mat),
    ArrayBuilder(vec),
    rows,
    cols,
    strategy)
end

function SparseMatrixAssembler(mat,trial::FESpace,test::FESpace)
  mat_builder = SparseMatrixBuilder(mat)
  T = eltype(get_array_type(mat_builder))
  SparseMatrixAssembler(mat_builder,Vector{T},trial,test)
end

"""
"""
function SparseMatrixAssembler(trial::FESpace,test::FESpace)
  T = get_dof_value_type(trial)
  matrix_type = SparseMatrixCSC{T,Int}
  vector_type = Vector{T}
  SparseMatrixAssembler(matrix_type,vector_type,trial,test)
end

get_rows(a::GenericSparseMatrixAssembler) = a.rows

get_cols(a::GenericSparseMatrixAssembler) = a.cols

get_matrix_builder(a::GenericSparseMatrixAssembler) = a.matrix_builder

get_vector_builder(a::GenericSparseMatrixAssembler) = a.vector_builder

get_assembly_strategy(a::GenericSparseMatrixAssembler) = a.strategy

function symbolic_loop_matrix!(A,a::GenericSparseMatrixAssembler,matdata)
  get_mat(a::Tuple) = a[1]
  get_mat(a) = a
  if LoopStyle(A) == DoNotLoop()
    return A
  end
  for (cellmat,_cellidsrows,_cellidscols) in zip(matdata...)
    cellidsrows = map_cell_rows(a.strategy,_cellidsrows)
    cellidscols = map_cell_cols(a.strategy,_cellidscols)
    rows_cache = array_cache(cellidsrows)
    cols_cache = array_cache(cellidscols)
    @assert length(cellidscols) == length(cellidsrows)
    if length(cellidscols) > 0
      mat1 = get_mat(first(cellmat))
      rows1 = getindex!(rows_cache,cellidsrows,1)
      cols1 = getindex!(cols_cache,cellidscols,1)
      touch! = TouchEntriesMap()
      touch_cache = return_cache(touch!,A,mat1,rows1,cols1)
      caches = touch_cache, rows_cache, cols_cache
      _symbolic_loop_matrix!(A,caches,cellidsrows,cellidscols,mat1)
    end
  end
  A
end

@noinline function _symbolic_loop_matrix!(A,caches,cell_rows,cell_cols,mat1)
  touch_cache, rows_cache, cols_cache = caches
  touch! = TouchEntriesMap()
  for cell in 1:length(cell_cols)
    rows = getindex!(rows_cache,cell_rows,cell)
    cols = getindex!(cols_cache,cell_cols,cell)
    evaluate!(touch_cache,touch!,A,mat1,rows,cols)
  end
end

function numeric_loop_matrix!(A,a::GenericSparseMatrixAssembler,matdata)
  for (cellmat,_cellidsrows,_cellidscols) in zip(matdata...)
    cellidsrows = map_cell_rows(a.strategy,_cellidsrows)
    cellidscols = map_cell_cols(a.strategy,_cellidscols)
    rows_cache = array_cache(cellidsrows)
    cols_cache = array_cache(cellidscols)
    vals_cache = array_cache(cellmat)
    @assert length(cellidscols) == length(cellidsrows)
    @assert length(cellmat) == length(cellidsrows)
    if length(cellmat) > 0
      mat1 = getindex!(vals_cache,cellmat,1)
      rows1 = getindex!(rows_cache,cellidsrows,1)
      cols1 = getindex!(cols_cache,cellidscols,1)
      add! = AddEntriesMap(+)
      add_cache = return_cache(add!,A,mat1,rows1,cols1)
      caches = add_cache, vals_cache, rows_cache, cols_cache
      _numeric_loop_matrix!(A,caches,cellmat,cellidsrows,cellidscols)
    end
  end
  A
end

@noinline function _numeric_loop_matrix!(mat,caches,cell_vals,cell_rows,cell_cols)
  add_cache, vals_cache, rows_cache, cols_cache = caches
  add! = AddEntriesMap(+)
  for cell in 1:length(cell_cols)
    rows = getindex!(rows_cache,cell_rows,cell)
    cols = getindex!(cols_cache,cell_cols,cell)
    vals = getindex!(vals_cache,cell_vals,cell)
    evaluate!(add_cache,add!,mat,vals,rows,cols)
  end
end

function symbolic_loop_vector!(b,a::GenericSparseMatrixAssembler,vecdata)
  @notimplementedif LoopStyle(b) == Loop()
  b
end

function numeric_loop_vector!(b,a::GenericSparseMatrixAssembler,vecdata)
  for (cellvec, _cellids) in zip(vecdata...)
    cellids = map_cell_rows(a.strategy,_cellids)
    rows_cache = array_cache(cellids)
    vals_cache = array_cache(cellvec)
    if length(cellvec) > 0
      vals1 = getindex!(vals_cache,cellvec,1)
      rows1 = getindex!(rows_cache,cellids,1)
      add! = AddEntriesMap(+)
      add_cache = return_cache(add!,b,vals1,rows1)
      caches = add_cache, vals_cache, rows_cache
      _numeric_loop_vector!(b,caches,cellvec,cellids)
    end
  end
  b
end

@noinline function _numeric_loop_vector!(vec,caches,cell_vals,cell_rows)
  add_cache, vals_cache, rows_cache = caches
  @assert length(cell_vals) == length(cell_rows)
  add! = AddEntriesMap(+)
  for cell in 1:length(cell_rows)
    rows = getindex!(rows_cache,cell_rows,cell)
    vals = getindex!(vals_cache,cell_vals,cell)
    evaluate!(add_cache,add!,vec,vals,rows)
  end
end

function symbolic_loop_matrix_and_vector!(A,b,a::GenericSparseMatrixAssembler,data)
  matvecdata, matdata, vecdata = data
  symbolic_loop_matrix!(A,a,matvecdata)
  symbolic_loop_matrix!(A,a,matdata)
  symbolic_loop_vector!(b,a,vecdata)
  A, b
end

function numeric_loop_matrix_and_vector!(A,b,a::GenericSparseMatrixAssembler,data)
  matvecdata, matdata, vecdata = data
  for (cellmatvec,_cellidsrows,_cellidscols) in zip(matvecdata...)
    cellidsrows = map_cell_rows(a.strategy,_cellidsrows)
    cellidscols = map_cell_cols(a.strategy,_cellidscols)
    rows_cache = array_cache(cellidsrows)
    cols_cache = array_cache(cellidscols)
    vals_cache = array_cache(cellmatvec)
    @assert length(cellidscols) == length(cellidsrows)
    @assert length(cellmatvec) == length(cellidsrows)
    if length(cellmatvec) > 0
      mat1, vec1 = getindex!(vals_cache,cellmatvec,1)
      rows1 = getindex!(rows_cache,cellidsrows,1)
      cols1 = getindex!(cols_cache,cellidscols,1)
      add! = AddEntriesMap(+)
      add_mat_cache = return_cache(add!,A,mat1,rows1,cols1)
      add_vec_cache = return_cache(add!,b,vec1,rows1)
      caches = add_mat_cache, add_vec_cache, vals_cache, rows_cache, cols_cache
      _numeric_loop_matvec!(A,b,caches,cellmatvec,cellidsrows,cellidscols)
    end
  end
  numeric_loop_matrix!(A,a,matdata)
  numeric_loop_vector!(b,a,vecdata)
  A, b
end

@noinline function _numeric_loop_matvec!(A,b,caches,cell_vals,cell_rows,cell_cols)
  add_mat_cache, add_vec_cache, vals_cache, rows_cache, cols_cache = caches
  add! = AddEntriesMap(+)
  for cell in 1:length(cell_cols)
    rows = getindex!(rows_cache,cell_rows,cell)
    cols = getindex!(cols_cache,cell_cols,cell)
    vals = getindex!(vals_cache,cell_vals,cell)
    matvals, vecvals = vals
    evaluate!(add_mat_cache,add!,A,matvals,rows,cols)
    evaluate!(add_vec_cache,add!,b,vecvals,rows)
  end
end


text:
funcs = ( dof,
          dof_residual,
          confint,
          deviance,
          nulldeviance,
          loglikelihood,
          nullloglikelihood,
          dispersion,
          nobs,
          stderror,
          vcov,
          residuals,
          predict,
          response,
          modelmatrix,
          weights,
          leverage
        )

@testset "Quantile regression: high-level function" begin
    τs = range(0.1, 0.9, step=0.1)
    βs = hcat(map(τ->RobustModels.interiormethod(X, y, τ)[1], τs)...)
    println("Coefficients: $(vcat(τs', βs))")
    @test size(βs) == (size(X, 2), length(τs))
end

@testset "Quantile regression: fit method" begin
    τ = 0.5
    # Formula, dense and sparse entry  and methods :cg and :chol
    @testset "Argument type: $(typeof(A))" for (A, b) in ((form, data), (X, y), (sX, y))
        m1 = fit(QuantileRegression, A, b; quantile=τ, verbose=false)
        m2 = quantreg(A, b; quantile=τ, verbose=false)
        @test_nowarn println(m2)
        @test all(coef(m1) .== coef(m2))

        # refit
        β = copy(coef(m2))
        refit!(m2, y; quantile=τ, verbose=false)
        @test all(coef(m2) .== β)

        # interface
        @testset "interface method: $(f)" for f in funcs
            # make sure the method is defined
            robvar = f(m1)
        end

        # later fit!
        m3 = fit(QuantileRegression, A, b; quantile=τ, dofit=false)
        @test all(0 .== coef(m3))
        fit!(m3; verbose=false)
        @test all(β .== coef(m3))

        # leverage weights
        @test_nowarn refit!(m3; correct_leverage=true)
    end
end

@testset "Quantile regression: different quantiles" begin
    τs = range(0.1, 0.9, step=0.1)
    m2 = fit(QuantileRegression, form, data; quantile=0.5, verbose=false)

    @testset "$(τ) quantile" for τ in τs
        m1 = fit(QuantileRegression, form, data; quantile=τ, verbose=false)
        println(m1)
        β = coef(m1)
        res = residuals(m1)
        ## The quantile regression line exactly passes through p points, with p number of columns of X.
        @test count(x->isapprox(x, 0; atol=1e-7), res) == length(β)
#        @test count(iszero, res) == length(β)

        # refit with new quantile
        refit!(m2; quantile=τ)
        @test all(coef(m1) .== coef(m2))
    end
end

@testset "Quantile regression: sparsity estimation" begin
    m2 = fit(QuantileRegression, form, data; quantile=0.25, verbose=false)
    s = RobustModels.location_variance(m2.model, false)

    τs = range(0.25, 0.75, step=0.25)
    @testset "(q, method, kernel): $(τ), $(method), $(kernel)" for τ in τs, method in (:jones, :bofinger, :hall_sheather), kernel in (:epanechnikov, :triangle, :window)
        if τ != m2.model.τ
            refit!(m2; quantile=τ)
            s = RobustModels.location_variance(m2.model, false)
        end

        si = RobustModels.location_variance(m2.model, false; bw_method=method, α=0.05, kernel=kernel)
        if method==:jones && kernel==:epanechnikov
            @test isapprox(s, si; rtol=1e-4)
        else
            @test !isapprox(s, si; rtol=1e-4)
        end
    end
end


text:
# Hotelling T² test

export OneSampleHotellingT2, EqualCovHotellingT2, UnequalCovHotellingT2

abstract type HotellingT2Test <: HypothesisTest end

default_tail(::HotellingT2Test) = :right
pvalue(T::HotellingT2Test; tail=:right) = pvalue(FDist(dof(T)...), T.F, tail=tail)

function show_params(io::IO, T::HotellingT2Test, indent="")
    println(io, indent, "number of observations: ", nobs(T))
    println(io, indent, "number of variables:    ", T.p)
    println(io, indent, "T² statistic:           ", T.T²)
    println(io, indent, "transformed statistic:  ", T.F)
    println(io, indent, "degrees of freedom:     ", dof(T))
    println(io, indent, "covariance matrix:")
    Base.print_matrix(io, T.S, indent^2)
end

## Utility functions

function checkdims(X::AbstractMatrix, Y::AbstractMatrix)
    nx, p = size(X)
    ny, q = size(Y)
    p == q || throw(DimensionMismatch("Inconsistent number of variables: $p, $q"))
    (nx > 0 && ny > 0) || throw(ArgumentError("Inputs must be non-empty"))
    return (p, nx, ny)
end

# Helper function for computing A'B⁻¹A when B is a covariance matrix
At_Binv_A(A::AbstractArray, B::AbstractArray) = A'*(B \ A)

## One sample

struct OneSampleHotellingT2 <: HotellingT2Test
    T²::Real
    F::Real
    n::Int
    p::Int
    μ₀::Vector
    x̄::Vector
    S::Matrix
end

StatsBase.nobs(T::OneSampleHotellingT2) = T.n
StatsBase.dof(T::OneSampleHotellingT2) = (T.p, T.n - T.p)

"""
    OneSampleHotellingT2(X::AbstractMatrix, μ₀=<zero vector>)

Perform a one sample Hotelling's ``T^2`` test of the hypothesis that the vector of
column means of `X` is equal to `μ₀`.
"""
function OneSampleHotellingT2(X::AbstractMatrix{T},
                              μ₀::AbstractVector=fill(middle(zero(T)), size(X, 2))) where T
    n, p = size(X)
    p == length(μ₀) ||
        throw(DimensionMismatch("Number of variables does not match number of means"))
    n > 0 || throw(ArgumentError("The input must be non-empty"))
    x̄ = vec(mean(X, dims=1))
    S = cov(X)
    T² = n * At_Binv_A(x̄ .- μ₀, S)
    F = (n - p) * T² / (p * (n - 1))
    return OneSampleHotellingT2(T², F, n, p, μ₀, x̄, S)
end

"""
    OneSampleHotellingT2(X::AbstractMatrix, Y::AbstractMatrix, μ₀=<zero vector>)

Perform a paired Hotelling's ``T^2`` test of the hypothesis that the vector of mean
column differences between `X` and `Y` is equal to `μ₀`.
"""
function OneSampleHotellingT2(X::AbstractMatrix{T}, Y::AbstractMatrix{S},
                              μ₀::AbstractVector=fill(middle(zero(T), zero(S)), size(X, 2))) where {T,S}
    p, nx, ny = checkdims(X, Y)
    nx == ny || throw(DimensionMismatch("Inconsistent number of observations: $nx, $ny"))
    p == length(μ₀) ||
        throw(DimensionMismatch("Number of variables does not match number of means"))
    return OneSampleHotellingT2(X .- Y, μ₀)
end

testname(::OneSampleHotellingT2) = "One sample Hotelling's T² test"
population_param_of_interest(T::OneSampleHotellingT2) = ("Mean vector", T.μ₀, T.x̄)

## Two sample

## Two sample: equal covariance

struct EqualCovHotellingT2 <: HotellingT2Test
    T²::Real
    F::Real
    nx::Int
    ny::Int
    p::Int
    Δ::Vector
    S::Matrix
end

"""
    EqualCovHotellingT2(X::AbstractMatrix, Y::AbstractMatrix)

Perform a two sample Hotelling's ``T^2`` test of the hypothesis that the difference in
the mean vectors of `X` and `Y` is zero, assuming that `X` and `Y` have equal covariance
matrices.
"""
function EqualCovHotellingT2(X::AbstractMatrix, Y::AbstractMatrix)
    p, nx, ny = checkdims(X, Y)
    Δ = vec(mean(X, dims=1) .- mean(Y, dims=1))
    S = poolcov!(cov(X), nx - 1, cov(Y), ny - 1) .* (inv(nx) .+ inv(ny))
    T² = At_Binv_A(Δ, S)
    F = T² * (nx + ny - p - 1) / (p * (nx + ny - 2))
    return EqualCovHotellingT2(T², F, nx, ny, p, Δ, S)
end

StatsBase.nobs(T::EqualCovHotellingT2) = (T.nx, T.ny)
StatsBase.dof(T::EqualCovHotellingT2) = (T.p, T.nx + T.ny - T.p - 1)

testname(::EqualCovHotellingT2) =
    "Two sample Hotelling's T² test (equal covariance matrices)"
population_param_of_interest(T::EqualCovHotellingT2) =
    ("Difference in mean vectors", zeros(eltype(T.Δ), T.p), T.Δ)

## Two sample: unequal covariance

# Store the denominator degrees of freedom in the type, since the computation
# is expensive and we don't want to redo it every time the user calls dof
struct UnequalCovHotellingT2 <: HotellingT2Test
    T²::Real
    F::Real
    nx::Int
    ny::Int
    p::Int
    ν::Int
    Δ::Vector
    S::Matrix
end

"""
    UnequalCovHotellingT2(X::AbstractMatrix, Y::AbstractMatrix)

Perform a two sample Hotelling's ``T^2`` test of the hypothesis that the difference in
the mean vectors of `X` and `Y` is zero, without assuming that `X` and `Y` have equal
covariance matrices.
"""
function UnequalCovHotellingT2(X::AbstractMatrix, Y::AbstractMatrix)
    p, nx, ny = checkdims(X, Y)
    Δ = vec(mean(X, dims=1) .- mean(Y, dims=1))
    Sx = cov(X) ./ nx
    Sy = cov(Y) ./ ny
    ST = Sx .+ Sy
    T² = At_Binv_A(Δ, ST)
    F = (nx + ny - p - 1) * T² / (p * (nx + ny - 2))
    tmp = Symmetric(ST) \ Δ
    iν = (dot(tmp, Sx * tmp) / T²)^2 / (nx - 1) + (dot(tmp, Sy * tmp) / T²)^2 / (ny - 1)
    ν = trunc(Int, inv(iν))
    return UnequalCovHotellingT2(T², F, nx, ny, p, ν, Δ, ST)
end

StatsBase.nobs(T::UnequalCovHotellingT2) = (T.nx, T.ny)
StatsBase.dof(T::UnequalCovHotellingT2) = (T.p, T.ν)

testname(::UnequalCovHotellingT2) =
    "Two sample Hotelling's T² test (unequal covariance matrices)"
population_param_of_interest(T::UnequalCovHotellingT2) =
    ("Difference in mean vectors", zeros(eltype(T.Δ), T.p), T.Δ)


text:
# This routine ABSOLUTELY should not allocate any memory, called inside time loop.
@inbounds @fastmath function record!(it::Int64, issp::Int64, iss::Int64, pac, pap)

    for ipw in pac.activepw
        rinterpolatew = pap[ipw].ss[issp].rinterpolatew
        rindices = pap[ipw].ss[issp].rindices

        for rfield in pac.rfields
            recs = pap[ipw].ss[issp].records[rfield]

            # copyto! before scalar indexing if GPU
            if (_fd.use_gpu)
                pw1 = pap[ipw].w1[:t][rfield]
                pw = pap[ipw].wr[rfield]
                copyto!(pw, pw1)
            else
                pw = pap[ipw].w1[:t][rfield]
            end
            for ir = 1:pac.ageom[ipw][iss].nr
                recs[it, ir] = 0.0
                for (i, ri) in enumerate(rindices[rfield][ir])
                    recs[it, ir] += pw[ri] * rinterpolatew[rfield][ir][i]
                end
            end
        end
    end
end

# receiver(pw, weights, rindices[ir] spray, pac::P_common{FdtdAcou}, si, ::p) =
# source_term * spray * pac.mod[:K][si]


text:
#!/usr/bin/env julia

using SimradEK60
using SimradEK60TestData
using Test

# Basic tests

@test octet2deg(-128) == -180

@test power2db(25.6f0) == log10(2f0)

# SONAR equation for EK60

α = 0.009841439f0
Pt = 1000.0f0
c = 1448.2969f0
G = 25.92f0

λ = 0.038113076f0

Ψ = -20.7f0
τ = 0.001024f0
Pr = -111.22823f0
_R = [18.1677435f0]
Sa = -0.49f0

_Sv = Sv(Pr,  λ, G, Ψ, c, α, Pt, τ, Sa,
       _R)

@test typeof(_Sv[1]) == Float32
@test _Sv[1] ≈ -94.1832

# All pings

ps =collect(pings(EK60_SAMPLE));
@test length(ps) == 1716

# 38 kHz pings

ps38 = [p for p in ps if p.frequency == 38000];

@test length(ps38) == 572

# 70 kHz pings

ps70 = [p for p in ps if p.frequency == 70000];

@test length(ps70) == 0

# 120 kHz pings

ps120 = [p for p in ps if p.frequency == 120000];

@test length(ps120) == 572

# 200 kHz pings

ps200 = [p for p in ps if p.frequency == 200000];

@test length(ps200) == 572

# 38 kHz volume backscatter

Sv38 = Sv(ps38);

m, n = size(Sv38)
@test m == 3782
@test n == 572

@test typeof(Sv38[1,1]) == Float32

# Power

p38 = power(ps38)
m, n = size(p38)
@test m == 3782
@test n == 572
@test typeof(p38[1,1]) == Int16

pdb38 = powerdb(ps38)
m, n = size(pdb38)
@test m == 3782
@test n == 572

@test typeof(pdb38[1,1]) == Float32

pdb120 = powerdb(ps120)
pdb200 = powerdb(ps200)

# Volume backscatter

Sv120 = Sv(ps120)
Sv200 = Sv(ps200)

m, n = size(Sv200)
@test m == 3782
@test n == 572

# Athwartships

at38 = athwartshipangle(ps38);
at120 = athwartshipangle(ps120);
at200 = athwartshipangle(ps200);

m, n = size(at38)
@test m == 3782
@test n == 572

# Alongships

al38 = alongshipangle(ps38);
al120 = alongshipangle(ps120);
al200 = alongshipangle(ps200);

m, n = size(al38)
@test m == 3782
@test n == 572

# Range

_R = R(ps38[1])
@test typeof(_R[1]) == Float32
@test _R[1] == 0
@test _R[end] == 699.4475f0
_R = R(ps38)
@test _R[1] == 0
@test _R[end] == 699.4475f0

function rmse(A, B)
    a = A .- B
    a = a.^2
    return sqrt(maximum(a))
end

function compare(A, B)
    x = rmse(A, B)
    x < 1e-4
end

function load_echoview_matrix(filename)

    # This is a crappy but fast CSV reader. I tried CSV.jl but it
    # segfaults. I tried CSVFiles but it takes (literally) minutes to
    # run.

    a = []
    open(filename) do file
        for ln in eachline(file)
            s = split(ln, ",")
            s = s[14:end]
            s = map(x->parse(Float64,x), s)
            push!(a,s)
        end
    end

    hcat(a[2:end]...)

end

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917.raw")

ps = collect(pings(filename))
ps38 = [p for p in ps if p.frequency == 38000]
ps120 = [p for p in ps if p.frequency == 120000]
ps200 = [p for p in ps if p.frequency == 200000]

# N.B. Echoview drops the first sample, so we'll do the same for
# comparison purposes.

Sv38 = Sv(ps38)[2:end,:]
Sv120 = Sv(ps120)[2:end,:]
Sv200 = Sv(ps200)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-38-uncal.sv.csv")
A = load_echoview_matrix(filename)

@test compare(A, Sv38)

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-120-uncal.sv.csv")
A = load_echoview_matrix(filename)
@test compare(A, Sv120)

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-200-uncal.sv.csv")
A = load_echoview_matrix(filename)
@test compare(A, Sv200)

TS38 = TS(ps38)[2:end,:]
TS120 = TS(ps120)[2:end,:]
TS200 = TS(ps200)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-38-uncal.ts.csv")
A = load_echoview_matrix(filename)
@test compare(A, TS38)

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-120-uncal.ts.csv")
A = load_echoview_matrix(filename)
@test compare(A, TS120)

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-200-uncal.ts.csv")
A = load_echoview_matrix(filename)
@test compare(A, TS200)

# Calibration values from JR230.ecs Beware that Echoview
# measures pulse length in mS not seconds!

Sv38cal = Sv(ps38;
             frequency=38000,
             gain=25.9400,
             equivalentbeamangle=-20.700001,
             soundvelocity=1462.0,
             absorptioncoefficient=0.001014,
             transmitpower=2000.00000,
             pulselength=0.001024,
             sacorrection=0.12)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-38-cal.sv.csv")
A = load_echoview_matrix(filename)
@test compare(A, Sv38cal)

Sv120cal = Sv(ps120;
              frequency=120000,
              gain=21.950,
              equivalentbeamangle=-20.700001,
              soundvelocity=1462.0,
              absorptioncoefficient=0.02683,
              transmitpower=500.00000,
              pulselength=0.001024,
              sacorrection=-0.05)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-120-cal.sv.csv")
A = load_echoview_matrix(filename)
@test compare(A, Sv120cal)

Sv200cal = Sv(ps200;
              frequency=200000,
              gain=23.9900,
              equivalentbeamangle=-19.600000,
              soundvelocity=1462.0,
              absorptioncoefficient=0.04023,
              transmitpower=300.00000,
              pulselength=0.001024,
              sacorrection=0.080)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-200-cal.sv.csv")
A = load_echoview_matrix(filename)
@test compare(A, Sv200cal)

TS38cal = TS(ps38;
             frequency=38000,
             gain=25.9400,
             soundvelocity=1462.0,
             absorptioncoefficient=0.001014,
             transmitpower=2000.000002)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-38-cal.ts.csv")
A = load_echoview_matrix(filename)
@test compare(A, TS38cal)

TS120cal = TS(ps120;
              frequency=120000,
              gain=21.950,
              soundvelocity=1462.0,
              absorptioncoefficient=0.02683,
              transmitpower=500.00000)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-120-cal.ts.csv")
A = load_echoview_matrix(filename)
@test compare(A, TS120cal)

TS200cal = TS(ps200;
              frequency=200000,
              gain=23.9900,
              soundvelocity=1462.0,
              absorptioncoefficient=0.04023,
              transmitpower=300.00000)[2:end,:]

filename = joinpath(EK60_DATA, "JR230-D20091215-T121917-200-cal.ts.csv")
A = load_echoview_matrix(filename)
@test compare(A, TS200cal)

# Test range for unequal ping lengths

filename = joinpath(EK60_DATA, "JR245-D20110116-T182142.raw")

ps = pings(filename) # Get the pings
ps38 = [p for p in ps if p.frequency == 38000] # Just 38 kHz pings

m, n = size(R(ps38))
@test m == 5355 # not 5100 !
@test n == 412

#


text:
# ---
# Test
# ---
using Test, StatsModels
using DataFrames

using Unfold
include("test_utilities.jl")

data,evts = loadtestdata("test_case_3a") #
f  = @formula 0~1+conditionA+continuousA # 1

# prepare data
data_r = reshape(data,(1,:))
data_r = vcat(data_r,data_r)#add second channel

#--------------------------#
## Mass Univariate Linear ##
#--------------------------#
data_e,times = Unfold.epoch(data=data_r,tbl=evts,τ=(-1.,1.9),sfreq=20) # cut the data into epochs
m_mul,m_mul_results = fit(UnfoldLinearModel,f,evts,data_e,times)
@test m_mul_results[(m_mul_results.channel.==1).&(m_mul_results.colname_basis .==0.1),:estimate] ≈ [2,3,4]

data_e_noreshape,times = Unfold.epoch(data=data,tbl=evts,τ=(-1.,1.9),sfreq=20) # cut the data into epochs
m_mul_noreshape,m_mul_results_noreshape = fit(UnfoldLinearModel,f,evts,data_e_noreshape,times)
@test m_mul_results_noreshape[(m_mul_results_noreshape.channel.==1).&(m_mul_results_noreshape.colname_basis .==0.1),:estimate] ≈ [2,3,4]
@test size(m_mul_results_noreshape)[1] ==size(m_mul_results)[1]/2

# Add Missing in Data
data_e_missing = data_e
data_e_missing[1,25,end-5:end] .= missing
m_mul_missing,m_mul_missing_results = Unfold.fit(UnfoldLinearModel,f,evts,data_e_missing,times)
@test m_mul_missing_results.estimate ≈ m_mul_results.estimate

# Special solver solver_lsmr_se with Standard Error
se_solver = solver=(x,y)->Unfold.solver_default(x,y,stderror=true)
m_mul_se,m_mul_results_se = Unfold.fit(UnfoldLinearModel,f,evts,data_e,times,solver=se_solver)
@test all(m_mul_results_se.estimate .≈ m_mul_results.estimate)
@test !all(isnothing.(m_mul_results_se.stderror ))

#---------------------------------#
## Timexpanded Univariate Linear ##
#---------------------------------#
basisfunction = firbasis(τ=(-1,1),sfreq=20,name="basisA")
m_tul,m_tul_results = fit(UnfoldLinearModel,f,evts,data_r,basisfunction)
@test isapprox(m_tul_results[(m_tul_results.channel.==1).&(m_tul_results.colname_basis .==0.1),:estimate],[2,3,4],atol=0.0001)

# test without reshape, i.e. 1 channel vector e.g. size(data) = (1200,)
m_tul_noreshape,m_tul_results_noreshape = fit(UnfoldLinearModel,f,evts,data,basisfunction);
@test size(m_tul_results_noreshape)[1] ==size(m_tul_results)[1]/2

# Test under missing data
data_missing = Array{Union{Missing,Number}}(undef,size(data_r))
data_missing .= data_r
data_missing[4500:4600] .= missing

m_tul_missing,m_tul_missing_results = fit(UnfoldLinearModel,f,evts,data_missing,basisfunction)
@test  isapprox(m_tul_missing_results.estimate , m_tul_results.estimate,atol=1e-4)  # higher tol because we remove stuff

## Test multiple basisfunctions
b1 = firbasis(τ=(-1,1),sfreq=20,name="basisA")
b2 = firbasis(τ=(-1,1),sfreq=20,name="basisB")

f1  = @formula 0~1+continuousA # 1
f2  = @formula 0~1+continuousA # 1

# Fast-lane new implementation
m,res = fit(UnfoldLinearModel,Dict(0=>(f1,b1),1=>(f1,b1)),evts,data,eventcolumn="conditionA")

# slow manual
X1 = designmatrix(UnfoldLinearModel,f1,filter(x->(x.conditionA==0),evts),b1)
X2 = designmatrix(UnfoldLinearModel,f2,filter(x->(x.conditionA==1),evts),b2)

@time m = unfoldfit(UnfoldLinearModel,X1+X2,data')
tmp = condense_long(m)

# test fast way & slow way to be identical
@test all(tmp.estimate .== res.estimate)

# runntime tests - does something explode?
for k in 1:4
    local f
    if k == 1
        f  = @formula 0~1
    elseif k == 2
        f  = @formula 0~1+conditionA
    elseif k == 3
        f  = @formula 0~0+conditionA
    elseif k == 4
        f  = @formula 0~1+continuousA
    end
    println("Testing Runtime $k with Formula:$f")
    Xs = designmatrix(UnfoldLinearModel,f,evts,basisfunction)

    # Fit the model
    df = unfoldfit(UnfoldLinearModel,Xs,data_e)
    c = condense_long(df,times)
    fit(UnfoldLinearModel,f,evts,data_e,times)
    fit(UnfoldLinearModel,f,evts,data,basisfunction)
end

# Special solver solver_lsmr_se with Standard Error
se_solver = solver=(x,y)->Unfold.solver_default(x,y,stderror=true)
m_tul_se,m_tul_results_se = fit(UnfoldLinearModel,f,evts,data_r,basisfunction,solver=se_solver)
@test all(m_tul_results_se.estimate .== m_tul_results.estimate)
@test !all(isnothing.(m_tul_results_se.stderror ))

#m_mul_se,m_mul_results_se = fit(UnfoldLinearModel,f,evts,data_e.+randn(size(data_e)).*5,times,solver=se_solver)
#plot_results(m_mul_results_se[m_mul_results_se.channel.==1,:],se=true)
#m_tul_se,m_tul_results_se = fit(UnfoldLinearModel,f,evts,data_r.+randn(size(data_r)).*5,basisfunction,solver=se_solver)
#plot_results(m_tul_results_se[m_tul_results_se.channel.==1,:],se=true)

##
data_long,evts_long = loadtestdata("test_case_1c") #
data_long = reshape(data_long,(1,:))

data_long = vcat(data_long,data_long)
f_long  = @formula 0~1
basisfunction_long = firbasis(τ=(-1,1),sfreq=1000,name="basisA")

@time designmatrix(UnfoldLinearModel,f_long,evts_long,basisfunction_long)
@time a,m_tul_long = fit(UnfoldLinearModel,f_long,evts_long,data_long,basisfunction_long);

@test isapprox(m_tul_long[(m_tul_long.channel.==1).&(m_tul_long.colname_basis .==0.1),:estimate],[2],atol=0.0001)
# ~21s, test_case_1c, sfreq = 1000, 6000 events
#
## Older numbers for "designmatrix" only. But see "benchmark/benchmarkjl" for better benchmarks
# new version 7s-10s, dataset4, sfreq=1000, 1200stim,
# ~13s, test_case_1c, sfreq = 1000, 6000 events

###############################
##  Mixed Model tests
###############################
data,evts = loadtestdata("testCase3",dataPath=(@__DIR__)*"/data") #
append!(data,zeros(1000))
data = reshape(data,1,:)
data = vcat(data,data)
data = data.+ 1*randn(size(data)) # we have to add minimal noise, else mixed models crashes.
data_missing = Array{Union{Missing,Number}}(undef,size(data))
data_missing .= data
data_missing[4500:4600] .= missing

categorical!(evts,:subject)
f  = @formula 0~1+condA+condB + (1+condA+condB|subject)
#f  = @formula 0~1 + (1|subject)

# cut the data into epochs
# TODO This ignores subject bounds
data_e,times = Unfold.epoch(data=data,tbl=evts,τ=(-1.,1.9),sfreq=10)
data_missing_e,times = Unfold.epoch(data=data_missing,tbl=evts,τ=(-1.,1.9),sfreq=10)
evts_e,data_e = Unfold.dropMissingEpochs(copy(evts),data_e)
evts_missing_e,data_missing_e = Unfold.dropMissingEpochs(copy(evts),data_missing_e)

######################
##  Mass Univariate Mixed
@time m_mum = fit(UnfoldLinearMixedModel,f,evts_e,data_e    ,times,contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()))
df = m_mum[2]
@test isapprox(df[(df.channel .== 1).&(df.term.=="condA: 1").&(df.colname_basis.==0.0),:estimate], [5.618,9.175],rtol=0.1)
@time m_mum = fit(UnfoldLinearMixedModel,f,evts_missing_e,data_missing_e    ,times,contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()))
df = m_mum[2]
@test isapprox(df[(df.channel .== 1).&(df.term.=="condA: 1").&(df.colname_basis.==0.0),:estimate], [5.618,9.175],rtol=0.1)

# Timexpanded Univariate Mixed
f  = @formula 0~1+condA+condB + (1+condA+condB|subject)
basisfunction = firbasis(τ=(-0.2,0.3),sfreq=10,name="ABC")
@time m_tum = fit(UnfoldLinearMixedModel,f,evts,data,basisfunction, contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()) )
df = m_tum[2]
@test isapprox(df[(df.channel .== 1).&(df.term.=="condA: 1").&(df.colname_basis.==0.0),:estimate], [5.618,9.175],rtol=0.1)

# missing data in LMMs
# not yet implemented
Test.@test_broken  m_tum = fit(UnfoldLinearMixedModel,f,evts,data_missing,basisfunction, contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()) )

evts.subjectB = evts.subject;
evts1 = evts[evts.condA.==0,:]
evts2 = evts[evts.condA.==1,:]

f0_lmm  = @formula 0~1+condB+(1|subject) + (1|subjectB)
@time m_tum = fit(UnfoldLinearMixedModel,f0_lmm,evts,data,basisfunction)

f1_lmm  = @formula 0~1+condB+(1|subject)
f2_lmm  = @formula 0~1+condB+(1|subjectB)

b1 = firbasis(τ=(-0.2,0.3),sfreq=10,name="A")
b2 = firbasis(τ=(-0.1,0.3),sfreq=10,name="B")

X1_lmm  = designmatrix(UnfoldLinearMixedModel,f1_lmm,evts1,b1)
X2_lmm  = designmatrix(UnfoldLinearMixedModel,f2_lmm,evts2,b2)

r = Unfold.unfoldfit(UnfoldLinearMixedModel,X1_lmm+X2_lmm,data);
df = condense_long(r)

@test isapprox(df[(df.channel .== 1).&(df.term.=="condB").&(df.colname_basis.==0.0),:estimate],[18.18,10.4],rtol=0.1)

# Fast-lane new implementation
m,res = fit(UnfoldLinearMixedModel,Dict(0=>(f1_lmm,b1),1=>(f2_lmm,b2)),evts,data,eventcolumn="condA")

if 1 == 0
    using WGLMakie,AlgebraOfGraphics
    m = mapping(:colname_basis,:estimate,color=:term,layout_x=:term,layout_y=:basisname)

    AlgebraOfGraphics.data(df[df.channel.==2,:]) * visual(Lines) * m  |> draw

end

##------

if 1 == 0
    # Fit mass-univariate 1st level for all subjects
    basisfunction = firbasis(τ=(-.1,.4),sfreq=10,name="A")
    resAll = DataFrame()
    f  = @formula 0~1+condA+condB
    for s in unique(evts.subject)
        ##
        from = minimum(evts[evts.subject.==s,:latency])-10
        to = maximum(evts[evts.subject.==s,:latency])+40
        to = min(to,size(data,1))
        evts_s = evts[evts.subject.==s,:]
        evts_s.latency .-= from
        m = fit(UnfoldLinearModel,f,evts_s,data[from:to],basisfunction,contrasts=Dict(:condA => EffectsCoding(), :condB => EffectsCoding()))
        m.results.subject = s
        append!(resAll,m.results)
    end

    results = resAll[resAll.term.=="(Intercept)",:]

    results[results.time .== .1,:]
end


text:
using MCMCChains
using Test

# Tests for missing values.

@testset "utils" begin
    x = [1, missing, 3, 2]
    @test MCMCChains.cummean(x) == [1., 1., 2., 2.]
end

@testset "stats" begin
    chn = Chains(randn(1000, 2, 2))

    # Call describe without missing values.
    describe(devnull, chn; showall=true)
    s1, s2 = MCMCChains.summarystats(chn),
      MCMCChains.quantile(chn)

    # Add missing values.
    chn_m = Chains(cat(chn.value, ones(1, 2, 2) .* missing, dims = 1))

    # Call describe with missing values.
    describe(devnull, chn_m; showall=true)
    m1, m2 = MCMCChains.summarystats(chn_m),
      MCMCChains.quantile(chn_m)

    @test s1[:,2:4] == m1[:,2:4]
end

@testset "diagnostic functions" begin
    nchains = 2
    chn = Chains(randn(5000, 2, nchains))

    # Add missing values.
    chn_m = Chains(cat(chn.value, ones(1, 2, nchains) .* missing, dims = 1))

    # Currently we only check if diag. functions throw an assertion if a value is missing.
    @test_throws AssertionError gelmandiag(chn_m)

    gw_1 = gewekediag(chn)
    gw_2 = gewekediag(chn_m)
    hd_1 = heideldiag(chn)
    hd_2 = heideldiag(chn_m)
    rf_1 = rafterydiag(chn)
    rf_2 = rafterydiag(chn_m)

    @testset "diagnostics missing tests" for i in 1:nchains
        @test all(gw_1[i].df == gw_2[i].df)
        @test all(hd_1[i].df == hd_2[i].df)
        @test all(rf_1[i].df == rf_2[i].df)
    end

    @test_throws AssertionError discretediag(chn_m)
end


text:
using Documenter, AugmentedGaussianProcesses

makedocs(modules=[AugmentedGaussianProcesses],
         format = Documenter.Writers.HTMLWriter.HTML(
         assets = ["assets/icon.ico"],
         analytics="UA-129106538-2"),
         sitename= "AugmentedGaussianProcesses",
         authors="Theo Galy-Fajou",
         pages = [
         "Home"=>"index.md",
         "Background"=>"background.md",
         "User Guide"=>"userguide.md",
         "Kernels"=>"kernel.md",
         "Examples"=>"examples.md",
         "Julia GP Packages"=>"comparison.md",
         "API"=>"api.md"
         ]
         )

deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/theogf/AugmentedGaussianProcesses.jl.git",
    target = "build"
)


text:
"""
    tenxpca(;tenxfile::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="sqrt", rowmeanlist::AbstractString="", rowvarlist::AbstractString="", colsumlist::AbstractString="", dim::Number=3, noversamples::Number=5, niter::Number=3, chunksize::Number=5000, group::AbstractString, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)

A randomized SVD.

Input Arguments
---------
- `tenxfile` : Julia Binary file generated by `OnlinePCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {sqrt,log,raw}-scaling of the value.
- `rowmeanlist` : The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `rowvarlist` : The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `colsumlist` : The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `dim` : The number of dimension of PCA.
- `noversamples` : The number of over-sampling.
- `niter` : The number of power interation.
- `chunksize` is the number of rows reading at once (e.g. 5000).
- `group` : The group name of 10XHDF5 (e.g. mm10).
- `initW` : The CSV file saving the initial values of eigenvectors.
- `initV` : The CSV file saving the initial values of loadings.
- `logdir` : The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration.
- `perm` : Whether the data matrix is shuffled at random.
- `cper` : Count per X (e.g. CPM: Count per million <1e+6>)

Output Arguments
---------
- `V` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `U` : Loading vectors of covariance matrix (No. rows of the data matrix × dim)
- `Scores` : Principal component scores
- `ExpVar` : Explained variance by the eigenvectors
- `TotalVar` : Total variance of the data matrix
"""

# Total Variance
function tv(TotalVar::Number, X::AbstractArray)
    l = size(X)[2]
    progress = Progress(l)
    for i in 1:l
        TotalVar = TotalVar + X[:,i]'*X[:,i]
        # Progress Bar
        next!(progress)
    end
    TotalVar
end

# Normalization of X
function tenxnormalizex(X, scale, cper, colsumvec)
    if scale == "sqrt"
        X = sqrt.(cper .* X ./ colsumvec')
    end
    if scale == "log"
        X = log10.(cper .* X ./ colsumvec' .+ 1)
    end
    return X
end

# Initialization (only TENXPCA)
function tenxinit(tenxfile::AbstractString, dim::Number, chunksize::Number, group::AbstractString, rowmeanlist::AbstractString, rowvarlist::AbstractString, colsumlist::AbstractString, initW::Union{Nothing,AbstractString}, initV::Union{Nothing,AbstractString}, logdir::Union{Nothing,AbstractString}, pca::TENXPCA, cper::Number, scale::AbstractString="sqrt", perm::Bool=false)
    N, M = tenxnm(tenxfile, group)
    # Eigen vectors
    if initW == nothing
        W = zeros(Float32, M, dim)
        for i=1:dim
            W[i, i] = 1
        end
    end
    if typeof(initW) == String
        if initV == nothing
            W = readcsv(initW, Float32)
        else
            error("initW and initV are not specified at once. You only have one choice.")
        end
    end
    if typeof(initV) == String
            V = readcsv(initV, Float32)
            V = V[:,1:dim]
    end
    D = Diagonal(reverse(1:dim)) # Diagonal Matrix
    rowmeanvec = zeros(Float32, N, 1)
    rowvarvec = zeros(Float32, N, 1)
    colsumvec = zeros(Float32, M, 1)
    if rowmeanlist != ""
        rowmeanvec = readcsv(rowmeanlist, Float32)
    end
    if rowvarlist != ""
        rowvarvec = readcsv(rowvarlist, Float32)
    end
    if colsumlist != ""
        colsumvec = readcsv(colsumlist, Float32)
    end
    # N, M, All Variance
    TotalVar = 0.0
    # Index Pointer
    idp = indptr(tenxfile, group)
    # Each chunk
    if N > chunksize
        lasti = fld(N, chunksize)+1
    else
        lasti = 1
    end
    for i in 1:lasti
        startp = Int64((i-1)*chunksize+1)
        endp = Int64(i*chunksize)
        if N - endp + chunksize < chunksize
            endp = N
        end
        println("Loading a chunk of sparse matrix from 10XHDF5...")
        X = loadchromium(tenxfile, group, idp, startp, endp, M, perm)
        X = tenxnormalizex(X, scale, cper, colsumvec)
        println("Calculating the total variance...")
        TotalVar = tv(TotalVar, X)
        # Creating W from V
        if typeof(initV) == String
            W = W .+ (V[n,:]*X')'
        end
    end
    TotalVar = TotalVar / M
    # directory for log file
    if logdir isa String
        if(!isdir(logdir))
            mkdir(logdir)
        end
    end
    return W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, idp
end

function tenxpca(;tenxfile::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="sqrt", rowmeanlist::AbstractString="", rowvarlist::AbstractString="", colsumlist::AbstractString="", dim::Number=3, noversamples::Number=5, niter::Number=3, chunksize::Number=5000, group::AbstractString, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)
    # Initial Setting
    # Input
    if !(scale in ["sqrt", "log", "raw"])
        error("scale must be specified as sqrt, log, or raw")
    end
    pca = TENXPCA()
    println("Initial Setting...")
    W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, idp = tenxinit(tenxfile, dim, chunksize, group, rowmeanlist, rowvarlist, colsumlist, initW, initV, logdir, pca, cper, scale, perm)
    # Perform PCA
    out = tenxpca(tenxfile, outdir, scale, rowmeanlist, rowvarlist, colsumlist, dim, noversamples, niter, chunksize, logdir, pca, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, perm, idp, group, cper)
    # Output
    if outdir isa String
        writecsv(joinpath(outdir, "Eigen_vectors.csv"), out[1])
        writecsv(joinpath(outdir, "Eigen_values.csv"), out[2])
        writecsv(joinpath(outdir, "Loadings.csv"), out[3])
        writecsv(joinpath(outdir, "Scores.csv"), out[4])
        writecsv(joinpath(outdir, "ExpVar.csv"), out[5])
        writecsv(joinpath(outdir, "TotalVar.csv"), out[6])
    end
    return out
end

function tenxpca(tenxfile, outdir, scale, rowmeanlist, rowvarlist, colsumlist, dim, noversamples, niter, chunksize, logdir, pca, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, perm, idp, group, cper)
    # Setting
    N, M = tenxnm(tenxfile, group)
    l = dim + noversamples
    @assert 0 < dim ≤ l ≤ min(N, M)
    Ω = rand(Float32, M, l)
    XΩ = zeros(Float32, N, l)
    Y = zeros(Float32, N, l)
    L = zeros(Float32, N, l)
    Q = zeros(Float32, N, l)
    B = zeros(Float32, l, M)
    QtX = zeros(Float32, l, M)
    Scores = zeros(Float32, M, dim)
    lasti = 0
    if N > chunksize
        lasti = fld(N, chunksize)+1
    else
        lasti = 1
    end

    # If not 0 the calculation is converged
    # Each epoch s
    println("Random Projection : Y = A Ω")
    for i in 1:lasti
        startp = Int64((i-1)*chunksize+1)
        endp = Int64(i*chunksize)
        if N - endp + chunksize < chunksize
            endp = N
        end
        X = loadchromium(tenxfile, group, idp, startp, endp, M, perm)
        X = tenxnormalizex(X, scale, cper, colsumvec)
        if rowvarlist != ""
            XΩ[startp:endp,:] .= (X ./ rowvarvec[startp:endp])*Ω
        else
            XΩ[startp:endp,:] .= X*Ω
        end
    end
    if rowvarlist != ""
        Y .= XΩ .- (rowmeanvec ./ rowvarvec) .* sum(Ω, dims=1) # N * l
    else
        Y .= XΩ .- rowmeanvec .* sum(Ω, dims=1) # N * l
    end

    # LU factorization
    println("LU factorization : L = lu(Y)")
    L .= lu!(Y).L # N * l

    for i in 1:niter
        println("##### "*string(i)*" / "*string(niter)*" niter #####")
        # Temporal matrix
        XL = zeros(Float32, M, l)
        AtL = zeros(Float32, M, l)
        XAtL = zeros(Float32, N, l)

        println("Normalized power iterations (1/3) : A' L")
        for j in 1:lasti
            startp = Int64((j-1)*chunksize+1)
            endp = Int64(j*chunksize)
            if N - endp + chunksize < chunksize
                endp = N
            end
            X = loadchromium(tenxfile, group, idp, startp, endp, M, perm)
            X = tenxnormalizex(X, scale, cper, colsumvec)
            if rowvarlist != ""
                XL .+= (X ./ rowvarvec[startp:endp])'*L[startp:endp,:]
            else
                XL .+= X'*L[startp:endp,:]
            end
        end
        if rowvarlist != ""
            AtL .= XL .- (rowmeanvec ./ rowvarvec)' * L
        else
            AtL .= XL .- rowmeanvec' * L
        end

        println("Normalized power iterations (2/3) : A A' L")
        for j in 1:lasti
            startp = Int64((j-1)*chunksize+1)
            endp = Int64(j*chunksize)
            if N - endp + chunksize < chunksize
                endp = N
            end
            X = loadchromium(tenxfile, group, idp, startp, endp, M, perm)
            X = tenxnormalizex(X, scale, cper, colsumvec)
            if rowvarlist != ""
                XAtL[startp:endp,:] .= (X ./ rowvarvec[startp:endp]) * AtL
            else
                XAtL[startp:endp,:] .= X * AtL
            end
        end
        if rowvarlist != ""
            Y .= XAtL .- (rowmeanvec ./ rowvarvec) .* sum(AtL, dims=1)
        else
            Y .= XAtL .- rowmeanvec .* sum(AtL, dims=1)
        end

        if i < niter
            println("Normalized power iterations (3/3) : L = lu(A A' L)")
            # Renormalize with LU factorization
            L .= lu!(Y).L
        else
            println("QR factorization  (3/3) : Q = qr(A A' L)")
            # Renormalize with QR factorization
            Q .= Array(qr!(Y).Q)
        end
    end

    println("Calculation of small matrix : B = Q' A")
    for j in 1:lasti
        startp = Int64((j-1)*chunksize+1)
        endp = Int64(j*chunksize)
        if N - endp + chunksize < chunksize
            endp = N
        end
        X = loadchromium(tenxfile, group, idp, startp, endp, M, perm)
        X = tenxnormalizex(X, scale, cper, colsumvec)
        if rowvarlist != ""
            QtX .+= ((X ./ rowvarvec[startp:endp])' * Q[startp:endp,:])'
        else
            QtX .+= (X' * Q[startp:endp,:])'
        end
    end
    if rowvarlist != ""
        B = QtX .- Q'*(rowmeanvec ./ rowvarvec)
    else
        B = QtX .- Q'*rowmeanvec
    end

    # SVD with small matrix
    println("SVD with small matrix : svd(B)")
    W, σ, V = svd(B)
    U = Q*W
    λ = σ .* σ ./ M
    # PC scores, Explained Variance
    for n = 1:dim
        Scores[:, n] .= λ[n] .* V[:, n]
    end
    ExpVar = sum(λ) / TotalVar
    # Return
    return (V[:,1:dim], λ[1:dim], U[:,1:dim], Scores[:,1:dim], ExpVar, TotalVar)
end