text:
export ConditionBlock, condition

struct ConditionBlock{N, BTT<:AbstractBlock{N}, BTF<:AbstractBlock{N}} <: CompositeBlock{N}
    m::Measure
    block_true::BTT
    block_false::BTF
    function ConditionBlock(m::Measure, block_true::BTT, block_false::BTF) where {N, BTT<:AbstractBlock{N}, BTF<:AbstractBlock{N}}
        new{N, BTT, BTF}(m, block_true, block_false)
    end
end

subblocks(c::ConditionBlock) = (c.m, c.block_true, c.block_false)
chsubblocks(c::ConditionBlock, blocks) = ConditionBlock(blocks...)

function apply!(reg::AbstractRegister{B}, c::ConditionBlock) where B
    if !isdefined(c.m, :results)
        println("Conditioned on a measurement that has not been performed.")
        throw(UndefRefError())
    end
    for i = 1:B
        viewbatch(reg, i) |> (c.m.results[i] == 0 ? c.block_false : c.block_true)
    end
    reg
end

condition(m, a::AbstractBlock{N}, b::Nothing) where N = ConditionBlock(m, a, eyeblock(N))
condition(m, a::Nothing, b::AbstractBlock{N}) where N = ConditionBlock(m, eyeblock(N), b)
mat(c::ConditionBlock) = throw(ArgumentError("ConditionBlock does not has matrix representation, try `mat(c.block_true)` or `mat(c.block_false)`"))

function print_block(io::IO, c::ConditionBlock)
    print(io, "if result(id = $(objectid(c.m)))")
end

function print_subblocks(io::IO, tree::ConditionBlock, depth, charset, active_levels)
    child_active_levels = active_levels

    print_prefix(io, depth+1, charset, active_levels)
    print_tree(
        io, tree.block_true;
        depth=depth+1,
        active_levels=child_active_levels,
        charset=charset,
        roottree=tree,
    )
    print_prefix(io, depth, charset, active_levels)
    println(io, "else")
    print_prefix(io, depth+1, charset, active_levels)
    print_tree(
        io, tree.block_false;
        depth=depth+1,
        active_levels=child_active_levels,
        charset=charset,
        roottree=tree,
    )
    print_prefix(io, depth, charset, active_levels)
    println(io, "end")
end


text:
const tutorialpath = normpath(joinpath(@__DIR__, "..", "examples", "tutorial"))
const multiregionpath = joinpath(tutorialpath, "02-multi-region-model")
const oneregionpath = joinpath(tutorialpath, "01-one-region-model")

function run_oneregion()
    include(joinpath(oneregionpath, "one-region-model.jl"))
end

function run_multiregion()
    include(joinpath(multiregionpath, "main.jl"))
end

ext: jl


text:
using Parameters

function compute_parent_face_categorical_weights(parent_contact_face, mixture_prob)
    @assert 0 < mixture_prob < 1
    if isnothing(parent_contact_face)
        probs = ones(length(S.BOX_SURFACE_IDS))
        probs /= sum(probs)
        return S.BOX_SURFACE_IDS, probs
    else
        faces = [:top :bottom;:left :right;:front :back]
        i,j = let
            idxs = findall(faces .== parent_contact_face)
            @assert length(idxs) == 1
            Tuple(idxs[1])
        end
        j = (j % 2) + 1
        weights = zeros(size(faces))
        weights[i,j] = 1.0

        uniform_weight = ones(size(faces))
        uniform_weight /= sum(uniform_weight)

        weights = weights * mixture_prob + uniform_weight * (1.0 - mixture_prob)

        return reshape(faces, (6,)), reshape(weights, (6,))
    end
end

@gen function sample_contact(boxes, parent_idx, parent_contact_face, hyperparams)
    # For both the parent object and child object, sample which of the 6 faces
    # of the bounding box are in contact.

    faces, probs = compute_parent_face_categorical_weights(parent_contact_face,
            hyperparams.parent_face_mixture_prob)
    parent_face ~ labeled_categorical(faces, probs)

    probs = ones(length(S.BOX_SURFACE_IDS))
    probs /= sum(probs)
    child_face ~ labeled_categorical(S.BOX_SURFACE_IDS, probs)

    parent_shape = boxes[parent_idx]
    if parent_face in [:top, :bottom]
        x ~ uniform(-parent_shape.sizeX/2, parent_shape.sizeX/2)
        y ~ uniform(-parent_shape.sizeY/2, parent_shape.sizeY/2)
    elseif parent_face in [:front, :back]
        x ~ uniform(-parent_shape.sizeX/2, parent_shape.sizeX/2)
        y ~ uniform(-parent_shape.sizeZ/2, parent_shape.sizeZ/2)
    elseif parent_face in [:left, :right]
        x ~ uniform(-parent_shape.sizeZ/2, parent_shape.sizeZ/2)
        y ~ uniform(-parent_shape.sizeY/2, parent_shape.sizeY/2)
    else
        @error "Invalid Parent Face"
    end

    # Sample the (x,y) coordinate along the parent's surface plane at which the contact occurs.
    # if length(shapes) == parent_idx
    # x ~ uniform(-parent_size/2, parent_size/2)
    # y ~ uniform(-parent_size/2, parent_size/2)
    # else
    #     x ~ normal(0.0, 0.25 * parent_size)
    #     y ~ normal(0.0, 0.25 * parent_size)
    # end

    # Sample the angle of rotation of the child object about the contact plane's normal vector.
    angle = R.rotation_angle({:angle} ~ GenDirectionalStats.uniform_rot2())
    # Since the contact may not be perfectly flush, we additionally sample a unit vector
    # representing the child object's outward contact normal (in the parent object's contact frame).
    # (This distribution is concentratedat [0,0,-1], a vector pointing exactly downward, which would
    # mean the contact planes are exactly aligned.
    slack_dir ~ GenDirectionalStats.vmf_3d_direction(GenDirectionalStats.UnitVector3(0.0, 0.0, -1.0), hyperparams.slack_dir_conc)
    # Sample a distance to shift the contact point along the sampled slack_dir
    slack_offset ~ normal(0, hyperparams.slack_offset_var)
    return S.ShapeContact(parent_face, Float64[], child_face, Float64[],
                          make_hopf_planar_contact(x, y, angle, slack_offset, slack_dir))
end

@gen function sample_floating_pose(hyperparams)
    # Sample a floating pose with position uniformly over a 3D region and orientation uniformly.
    pose ~ uniformPose(hyperparams.floating_position_bounds...)
end

@gen function sample_scene_graph_params(structure, boxes, hyperparams)
    scene_graph = S.SceneGraph()

    for (idx, bbox) in enumerate(boxes)
        # Each Shape object has an associated bounding box.
        S.addObject!(scene_graph, obj_name_from_idx(idx), bbox)
    end

    roots = S.rootsOfForest(structure)
    dists = LG.dijkstra_shortest_paths(structure, roots).dists
    unique_dists = sort(unique(dists))

    # Given the structure (a graph) sample parameters of the scene graph.
    for d in unique_dists
        for idx in findall(dists .== d)
            if isFloating(structure, idx)
                # For floating (root nodes in the graph), sample a floating pose and set the pose in the scene graph.
                pose = {idx => :floating} ~ sample_floating_pose(hyperparams)
                S.setPose!(scene_graph, obj_name_from_idx(idx), pose)
            else
                # For contacting objects, sample parameters of the contact and set the contact relationship in the scene graph.
                parent_idx = parent(structure, idx)
                parent_contact_face = let
                    parents_parent_idx = parent(structure, parent_idx)
                    face=nothing
                    if !isnothing(parents_parent_idx)
                        shape_contact = S.getContact(scene_graph,
                            obj_name_from_idx(parents_parent_idx),obj_name_from_idx(parent_idx))
                        face = shape_contact.childFamilyId
                    end
                    face
                end
                contact = {idx => :contact} ~ sample_contact(boxes, parent_idx, parent_contact_face, hyperparams)
                S.setContact!(scene_graph, obj_name_from_idx(parent_idx), obj_name_from_idx(idx), contact)
            end
        end
    end
    return scene_graph
end

@gen (static) function scene_graph_prior(model_params)
    structure ~ uniform_diforest(length(model_params.boxes))
    scene_graph = {:params} ~ sample_scene_graph_params(structure, model_params.boxes, model_params.hyperparams)
    return scene_graph
end

@with_kw struct Hyperparams
    slack_dir_conc::Float32
    slack_offset_var::Float32
    p_outlier::Real
    noise::Real
    resolution::Real
    parent_face_mixture_prob::Real
    floating_position_bounds::Tuple
end

@with_kw struct SceneGraphPriorModelParameters
    boxes::Vector{S.Box}
    hyperparams::Hyperparams
end

@with_kw struct SceneModelParameters
    boxes::Vector{S.Box}
    hyperparams::Hyperparams
    get_cloud_from_poses_and_idx::Function
    camera_pose::Pose
    N::Int
end

export Hyperparams, SceneGraphPriorModelParameters, SceneModelParameters

function render_clouds(model_params::SceneModelParameters, scene_graph)
    # Since each type of object has an associated shape distribution, we approximately integrate out the shape variable
    # by sampling shapes from each object type's shape distribution and then summing the weighted scores for these samples.
    # This sampling does not occur in this generative function, instead the samples are stored in the `meshes` field of each
    # Shape object and we expect there to be `N` of those samples in the `meshes` list. Below, we render `N` point clouds of the
    # scene, where the i-th point cloud is generated by taking the i-th sample from each object's `meshes` list and rendering a
    # scene using those shapes (at the poses sampled in the above code).

    # Floating poses of all objects as sampled above.
    poses = floatingPosesOf(scene_graph)
    # N rendered point clouds
    rendered_clouds = [
        let
            rendered_point_cloud_in_camera_frame = model_params.get_cloud_from_poses_and_idx(
                poses, i, model_params.camera_pose
            )
            rendered_point_cloud_in_world_frame = move_points_to_frame_b(rendered_point_cloud_in_camera_frame,model_params.camera_pose)
            voxelize(rendered_point_cloud_in_world_frame, model_params.hyperparams.resolution)
        end
        for i in 1:model_params.N
    ]
    rendered_clouds
end

"""
Sample a scene
    shapes : Vector of `Shape` objects
    get_cloud_from_objects : Function to render Objects (consisting of a pose and a mesh) from a specifried camera viewpoint
    camera_pose : Pose of camera
    likelihood_params : Hyperparameters of the likelihood model.
    N : Number of shape samples to approximately integrate over
"""
@gen (static) function scene(model_params::SceneModelParameters)
    scene_graph ~ scene_graph_prior(model_params)
    rendered_clouds = render_clouds(model_params, scene_graph)
    p_outlier ~ exponential(1/model_params.hyperparams.p_outlier)
    noise ~ exponential(1/model_params.hyperparams.noise)
    obs_cloud = {:obs} ~ uniform_mixture_from_template_multi_cloud(
                                                rendered_clouds,
                                                p_outlier,
                                                noise,
                                                (-100.0,100.0,-100.0,100.0,-100.0,300.0))
    return (scene_graph=scene_graph,
            rendered_clouds=rendered_clouds,
            obs_cloud=obs_cloud)
end

@gen (static) function scene_graph_prior_model(model_params)
    scene_graph ~ scene_graph_prior(model_params)
    return (scene_graph=scene_graph)
end

Gen.@load_generated_functions

export scene, scene_graph_prior_model


text:
module Trebuchet

using WebIO, JSExpr
using OrdinaryDiffEq

export TrebuchetState, simulate, visualise, run, shoot

include("defs.jl")
include("utils.jl")
include("derive.jl")
include("transition.jl")
include("simulate.jl")
include("visualise.jl")

end


text:
# Copyright (c) 2020 California Institute of Technology (“Caltech”). U.S.
# Government sponsorship acknowledged.
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#   • Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   • Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   • Neither the name of Caltech nor its operating division, the Jet
#   Propulsion Laboratory, nor the names of its contributors may be used to
#   endorse or promote products derived from this software without specific
#   prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
# OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
# WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
# OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
# ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

module SphericalFourierBesselDecompositions

include("HealPy.jl")
include("healpix_helpers.jl")
include("Splines.jl")
include("gnl.jl")  # SphericalBesselGnls
include("modes.jl")  # AnlmModes, ClnnModes, ClnnBinnedModes
include("SeparableArrays.jl")
include("windows.jl")  # window function related
include("theory.jl")  # mostly for testing the package, may be split at some point
include("NDIterators.jl")
include("wigner_chains.jl")  # Wigner symbol related
include("window_chains.jl")  # window function related
include("covariance.jl")  # mostly theory, may be split at some point
include("cat2anlm.jl")

using Statistics
#using Roots
using Healpix
using .HealPy
using .HealpixHelpers
using .Splines
using .GNL
using .Modes
using .Windows
using .SeparableArrays
using .Theory
using .Covariance
using .WignerChains
using .WindowChains
using .NDIterators
using .Cat2Anlm

#using QuadGK
#using FastGaussQuadrature
using Distributed
#using Base.Threads

## for testing only
#using MeasureAngularPowerSpectra
#using Profile
#using PyPlot

###################### utility functions #########################

# gen_withcache()
function gen_withcache(nnn, fname, gen_fn)
    dims = length(nnn)
    @assert dims == 2
    nnn_cache = fill(Int64(0), dims)
    arr = fill(Float64(NaN), nnn_cache...)
    if isfile(fname)  # check what's in the cache
        read!(fname, nnn_cache)
    end
    if all(nnn_cache .>= nnn)  # cache is sufficient
        open(fname, "r") do f
            read!(f, nnn_cache)
            arr = read!(f, Array{Float64}(undef, nnn_cache...))
            arr = collect(arr[1:nnn[1],1:nnn[2]])  # can this be generalized to N-dims
        end
    else  # cache is not sufficient
        arr = gen_fn(nnn)
        @assert all(size(arr) .== nnn)
        open(fname, "w") do f
            nnn_cache = [nnn[i] for i=1:dims]  # 'write()' cannot do tuples
            write(f, Int64.(nnn_cache))
            write(f, Float64.(arr))
        end
    end
    return arr
end

####################### a_nlm -> c_lnn #####################################

function alm2cl(alm1, alm2, lmax, lmin=0)
    # Note: our alm are in a different order than Healpix. In our array the
    # (l,m) pairs occur in the order
    # [(0,0),(1,0),(1,1),(2,0),(2,1),(2,2),(3,0),...].
    cl = fill(NaN, lmax+1 - lmin)
    b = 1
    for l=lmin:lmax
        c = real(alm1[b] * alm2[b])
        for m=1:l
            c += 2 * real(alm1[b+m] * conj(alm2[b+m]))
        end
        cl[l+1-lmin] = c / (2 * l + 1)
        b += l + 1
    end
    return cl
end

function amln2clnn(anlm, cmodes::ClnnModes)
    clnn = fill(NaN, getlnnsize(cmodes))
    nmax = cmodes.amodes.nmax
    for n̄=1:nmax, Δn=0:cmodes.Δnmax
        n1 = n̄ + Δn
        n2 = n̄
        if n1 > nmax || n2 > nmax
            continue
        end
        lmax = minimum(cmodes.amodes.lmax_n[[n1,n2]])
        lmsize = getlmsize(lmax)
        n1_idxs = getnlmsize(cmodes.amodes, n1 - 1) .+ (1:lmsize)
        n2_idxs = getnlmsize(cmodes.amodes, n2 - 1) .+ (1:lmsize)
        lnn_idxs = getidx.(cmodes, 0:lmax, n1, n2)
        clnn[lnn_idxs] .= alm2cl(anlm[n1_idxs], anlm[n2_idxs], lmax)
    end
    return clnn
end

@doc raw"""
    pixwin(cmodes)

    Return the angular pixel window for application to a Clnn object.
"""
function pixwin(cmodes)
    hppwin = Healpix.pixwin(cmodes.amodes.nside, pol=false)
    lnnsize = getlnnsize(cmodes)
    sfbpwin = fill(NaN, lnnsize)
    for i=1:lnnsize
        l, = getlnn(cmodes, i)
        sfbpwin[i] = hppwin[l+1]
    end
    return sfbpwin
end

############################ convenience functions for plotting #####################

function get_ck(knl, clkk, l, dn; allmodes=false)
    #@show l, dn, lmin
    lmax, nmax, nmax = size(clkk) .+ (-1,0,0)
    if nmax - dn < 0
        return Float64[], Float64[]
    end
    nmin = (!allmodes && l == 0) ? 2 : 1
    k = fill(NaN, nmax-dn-nmin+1)
    ck = fill(NaN, nmax-dn-nmin+1)
    for n1=nmin:nmax-dn
        n2 = n1 + dn
        k[n1-nmin+1] = (knl[n1,l+1] + knl[n2,l+1]) / 2
        ck[n1-nmin+1] = (clkk[l+1,n1,n2] + clkk[l+1,n2,n1]) / 2
    end
    return k, ck
end

function get_ck_err(knl, clkk, l, dn; allmodes=false)
    k0, cl0 = get_ck(knl, clkk, l, 0; allmodes=allmodes)
    k, cl = get_ck(knl, clkk, l, dn; allmodes=allmodes)

    cl0spline = Spline1D(k0, cl0)
    Δcl² = @. (cl0spline(k)^2 + cl^2) / (2*l + 1)

    return @. √Δcl²
end

function predict_clnn(knl, pk, rmax, lmin=0)
    kernel(k, k_n, l) = begin
        #@show typeof(k) typeof(k_n) typeof(l)
        k^2 * sphericalbesselj(l-1, k * rmax) / (k_n^2 - k^2)
    end
    nmax, lmax = size(knl) .+ (0,lmin-1)
    clnn = fill(NaN, lmax+1-lmin, nmax, nmax)
    for l=lmin:lmin, n1=1:nmax, n2=max(n1-1,1):min(n1+1,nmax)
        n2 = n1
        k_n1 = knl[n1,l+1-lmin]
        k_n2 = knl[n2,l+1-lmin]
        integrand(k) = pk(k) * kernel(k, k_n1, l) * kernel(k, k_n2, l)
        I,E = quadgk(integrand, 0.0, 10 * maximum(knl))  # TODO: infinity?
        clnn[l+1-lmin,n1,n2] = I
    end
    return (4 * rmax / π) * clnn
end

############################ Window function stuff #####################

function gen_mask(nside, fsky)
    npix = nside2npix(nside)
    mask = fill(0.0, npix)
    θ, ϕ = pix2angRing(nside, 1:npix)
    @show extrema(θ) extrema(ϕ)
    θₘₐₓ = acos(1 - 2*fsky)
    @show θₘₐₓ
    for i=1:length(θ)
        if θ[i] <= θₘₐₓ
            mask[i] = 1.0
        end
    end
    return mask
end

function gen_win_insep_cossin(mask, rmin, rmax, r, l, m)
    nr = length(r)
    win = fill(0.0, nr, length(mask))
    nside = npix2nside(length(mask))
    pix = (1:length(mask))[mask .!= 0]
    θ, ϕ = pix2angRing(nside, pix)
    @show extrema(ϕ)
    θmin, θmax = extrema(θ)
    dr = rmax - rmin
    dθ = θmax - θmin
    rmid = (rmax + rmin) / 2
    #@show l,m dr dθ
    r0fn(θ, ϕ) = begin
        rmid - (dr/2) * cos(θ * l*π/dθ) * cos(m*ϕ)
    end
    phifn(r,θ,ϕ) = exp(- ((r - r0fn(θ,ϕ)) / dr)^2)
    for i=1:length(pix)
        #win[:,pix[i]] .= r0fn(θ[i], ϕ[i])
        win[:,pix[i]] .= phifn.(r, θ[i], ϕ[i])
    end
    win ./= maximum(win)
    return win
end

function make_window(wmodes::ConfigurationSpaceModes, features...)
    rmin = wmodes.rmin
    rmax = wmodes.rmax
    r = wmodes.r
    Δr = wmodes.Δr
    nr = wmodes.nr
    nside = wmodes.nside
    @show rmin rmax nr Δr nside

    #features = features..., :separable

    phi = fill(1.0, nr)
    mask = fill(1.0, wmodes.npix)
    features = filter(i -> i != :fullsky, features)

    if :ang_75 in features
        mask = gen_mask(nside, 0.75)
        features = filter(i -> i != :ang_75, features)
    end

    if :ang_half in features
        mask = gen_mask(nside, 1/2)
        features = filter(i -> i != :ang_half, features)
    end

    if :ang_quarter in features
        mask = gen_mask(nside, 1/4)
        features = filter(i -> i != :ang_quarter, features)
    end

    if :ang_eighth in features
        mask = gen_mask(nside, 1/8)
        features = filter(i -> i != :ang_eighth, features)
    end

    if :ang_sixteenth in features
        mask = gen_mask(nside, 1/16)
        features = filter(i -> i != :ang_sixteenth, features)
    end

    if :rotate in features
        r = hp.Rotator(coord=["E", "G"])
        mask = r.rotate_map_pixel(mask)
        features = filter(i -> i != :rotate, features)
    end

    if :linear_mask in features
        pix = @. mask .!= 0
        features = filter(i -> i != :linear_mask, features)
    end

    if :radial in features
        r0 = rmax * 0.55
        for i=1:nr
            phi[i] = exp(- (r[i] / r0)^2)
        end
        features = filter(i -> i != :radial, features)
    end

    if :step_rmin in features
        r0 = rmin + 100.0
        for i=1:nr
            phi[i] = (r[i] < r0) ? 0.0 : phi[i]
        end
        features = filter(i -> i != :step_rmin, features)
    end

    if :step_rmax in features
        r0 = rmax - 100.0
        for i=1:nr
            phi[i] = (r[i] > r0) ? 0.0 : phi[i]
        end
        features = filter(i -> i != :step_rmax, features)
    end

    win = phi * mask'
    maxwin = maximum(win)

    if :separable in features
        phi ./= maxwin
        win = @SeparableArray phi mask
        features = filter(i -> i != :separable, features)
    else
        win ./= maxwin
    end

    # :radial_cossin_l05_m0, :radial_cossin_l1_m1, ...
    lmax = mmax = 5  # independent l and m, should probably call them something else
    for l=0:0.5:lmax, m=0:0.5:mmax
        lshort = filter(c -> c != '.', string(l))
        mshort = filter(c -> c != '.', string(m))
        sym = Symbol("radial_cossin_l$(lshort)_m$(mshort)")
        #@show sym
        if sym in features
            win = gen_win_insep_cossin(mask, rmin, rmax, r, l, m)
            features = filter(i -> i != sym, features)
        end
    end

    @assert maximum(win) == 1

    if length(features) != 0
        @warn "Did not recognize all features" features
    end

    return win
end

@doc raw"""
    calc_fsky(win, wmodes)

This functions returns a measure of the sky fraction covered by the survey with
window `win`. The exact implementation is considered an implementation detail
and can change in the future.
"""
function calc_fsky(win, wmodes)
    meanmask = mean(win, dims=1)[:]
    fsky = sum(meanmask .> 0) ./ length(meanmask)
end

end

# vim: set sw=4 et sts=4 :


text:
# Lazy concatenation of AbstractVector's.
# Similar to Iterators.Flatten and some code has been reused from julia/base/iterators.jl

# This fixes type inference issues in v1.0
for op in  (:hcat, :vcat)
    @eval call(::ApplyLayout{typeof($op)}, _) = $op
end

const Vcat{T,N,I<:Tuple} = ApplyArray{T,N,typeof(vcat),I}

Vcat(A...) = ApplyArray(vcat, A...)
Vcat{T}(A...) where T = ApplyArray{T}(vcat, A...)
Vcat() = Vcat{Any}()

Vcat(A::AbstractVector...) = ApplyVector(vcat, A...)
Vcat{T}(A::AbstractVector...) where T = ApplyVector{T}(vcat, A...)

function instantiate(A::Applied{DefaultApplyStyle,typeof(vcat)})
    isempty(A.args) && return A
    m = size(A.args[1],2)
    for k=2:length(A.args)
        size(A.args[k],2) == m || throw(ArgumentError("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
    end
    Applied{DefaultApplyStyle}(A.f,map(instantiate,A.args))
end

@inline eltype(A::Applied{<:Any,typeof(vcat)}) = promote_type(map(eltype,A.args)...)
@inline eltype(A::Applied{<:Any,typeof(vcat),Tuple{}}) = Any
@inline ndims(A::Applied{<:Any,typeof(vcat),I}) where I = max(1,maximum(map(ndims,A.args)))
@inline ndims(A::Applied{<:Any,typeof(vcat),Tuple{}}) = 1
@inline axes(f::Vcat{<:Any,1,Tuple{}}) = (OneTo(0),)
@inline axes(f::Vcat{<:Any,1}) = tuple(oneto(+(map(length,f.args)...)))
@inline axes(f::Vcat{<:Any,2}) = (oneto(+(map(a -> size(a,1), f.args)...)), axes(f.args[1],2))
Base.IndexStyle(::Type{<:Vcat{T,1}}) where T = Base.IndexLinear()
Base.IndexStyle(::Type{<:Vcat{T,2}}) where T = Base.IndexCartesian()

function ==(a::Vcat{T,N}, b::Vcat{T,N}) where {N,T}
    a_args = arguments(vcat, a)
    b_args = arguments(vcat, b)
    if length(a_args) ≠ length(b_args) || any(map(size,a_args) .≠ map(size,b_args))
        return Base.invoke(==, NTuple{2,AbstractArray}, a, b)
    end
    all(a_args .== b_args)
end

@propagate_inbounds @inline vcat_getindex(f, idx...) =
    vcat_getindex_recursive(f, idx, f.args...)

@propagate_inbounds @inline function vcat_getindex_recursive(
        f, idx::Tuple{Integer}, A, args...)
    k, = idx
    T = eltype(f)
    n = length(A)
    k ≤ n && return convert(T, A[k])::T
    vcat_getindex_recursive(f, (k - n, ), args...)
end

@propagate_inbounds @inline function vcat_getindex_recursive(
        f, idx::Tuple{Integer,Integer}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 1)
    k ≤ n && return convert(T, A[k, j])::T
    vcat_getindex_recursive(f, (k - n, j), args...)
end

@propagate_inbounds @inline function vcat_getindex_recursive(
        f, idx::Tuple{Integer,Union{Colon,AbstractVector}}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 1)
    k ≤ n && return convert(AbstractVector{T}, A[k, j])
    vcat_getindex_recursive(f, (k - n, j), args...)
end

@inline vcat_getindex_recursive(f, idx) = throw(BoundsError(f, idx))

@propagate_inbounds @inline getindex(f::Vcat{<:Any,1}, k::Integer) = vcat_getindex(f, k)
@propagate_inbounds @inline getindex(f::Vcat{<:Any,2}, k::Integer, j::Integer) = vcat_getindex(f, k, j)
@propagate_inbounds @inline getindex(f::Vcat{<:Any,2}, k::Integer, j::AbstractVector) = vcat_getindex(f, k, j)
@propagate_inbounds @inline getindex(f::Vcat{<:Any,2}, k::Integer, j::Colon) = vcat_getindex(f, k, j)

getindex(f::Applied{DefaultArrayApplyStyle,typeof(vcat)}, k::Integer)= vcat_getindex(f, k)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(vcat)}, k::Integer, j::Integer)= vcat_getindex(f, k, j)
getindex(f::Applied{<:Any,typeof(vcat)}, k::Integer)= vcat_getindex(f, k)
getindex(f::Applied{<:Any,typeof(vcat)}, k::Integer, j::Integer)= vcat_getindex(f, k, j)

# since its mutable we need to make a copy
copy(f::Vcat) = Vcat(map(copy, f.args)...)
map(::typeof(copy), f::Vcat) = Vcat(map.(copy, f.args)...)

@propagate_inbounds @inline function vcat_setindex_recursive!(
        f::Vcat{T,1} where T, v, idx::NTuple{1}, A, args...)
    k, = idx
    n = length(A)
    k ≤ n && return setindex!(A, v, idx...)
    vcat_setindex_recursive!(f, v, (k - n, ), args...)
end

@propagate_inbounds @inline function vcat_setindex_recursive!(
        f::Vcat{T,2} where T, v, idx::NTuple{2}, A, args...)
    k, j = idx
    n = size(A, 1)
    k ≤ n && return setindex!(A, v, idx...)
    vcat_setindex_recursive!(f, v, (k - n, j), args...)
end

@inline vcat_setindex_recursive!(f, v, idx) = throw(BoundsError(f, idx))

@propagate_inbounds @inline function setindex!(
        f::Vcat{T,N}, v, idx::Vararg{Integer,N}) where {T,N}
    vcat_setindex_recursive!(f, v, idx, f.args...)
end

reverse(f::Vcat{<:Any,1}) = Vcat((reverse(itr) for itr in reverse(f.args))...)

####
# Hcat
####

const Hcat{T,I<:Tuple} = ApplyArray{T,2,typeof(hcat),I}

Hcat(A...) = ApplyArray(hcat, A...)
Hcat{T}(A...) where T = ApplyArray{T}(hcat, A...)

function instantiate(A::Applied{DefaultApplyStyle,typeof(hcat)})
    isempty(A.args) && return A
    m = size(A.args[1],1)
    for k=2:length(A.args)
        size(A.args[k],1) == m || throw(ArgumentError("number of rows of each array must match (got $(map(x->size(x,1), A)))"))
    end
    Applied{DefaultApplyStyle}(A.f,map(instantiate,A.args))
end

@inline eltype(A::Applied{<:Any,typeof(hcat)}) = promote_type(map(eltype,A.args)...)
ndims(::Applied{<:Any,typeof(hcat)}) = 2
size(f::Applied{<:Any,typeof(hcat)}) = (size(f.args[1],1), +(map(a -> size(a,2), f.args)...))
Base.IndexStyle(::Type{<:Hcat}) where T = Base.IndexCartesian()

@inline hcat_getindex(f, k, j::Integer) = hcat_getindex_recursive(f, (k, j), f.args...)

@inline function hcat_getindex_recursive(f, idx::Tuple{Integer,Integer}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 2)
    j ≤ n && return convert(T, A[k, j])::T
    hcat_getindex_recursive(f, (k, j - n), args...)
end

@inline function hcat_getindex_recursive(f, idx::Tuple{Union{Colon,AbstractVector},Integer}, A, args...)
    kr, j = idx
    T = eltype(f)
    n = size(A, 2)
    j ≤ n && return convert(AbstractVector{T}, A[kr, j])
    hcat_getindex_recursive(f, (kr, j - n), args...)
end

@inline hcat_getindex_recursive(f, idx) = throw(BoundsError(f, idx))

getindex(f::Hcat, k::Integer, j::Integer) = hcat_getindex(f, k, j)
getindex(f::Hcat, k::AbstractVector, j::Integer) = hcat_getindex(f, k, j)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(hcat)}, k::Integer, j::Integer)= hcat_getindex(f, k, j)
getindex(f::Applied{<:Any,typeof(hcat)}, k::Integer, j::Integer)= hcat_getindex(f, k, j)

# since its mutable we need to make a copy
copy(f::Hcat) = Hcat(map(copy, f.args)...)

@inline function hcat_setindex_recursive!(f, v, idx::NTuple{2}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 2)
    j ≤ n && return setindex!(A, v, k, j)
    hcat_setindex_recursive!(f, v, (k, j - n), args...)
end

@inline hcat_setindex_recursive!(f, v, idx) = throw(BoundsError(f, idx))

function setindex!(f::Hcat{T}, v, k::Integer, j::Integer) where T
    hcat_setindex_recursive!(f, v, (k, j), f.args...)
end

####
# Hvcat
####

@inline eltype(A::Applied{<:Any,typeof(hvcat)}) = promote_type(map(eltype,tail(A.args))...)
ndims(::Applied{<:Any,typeof(hvcat)}) = 2
function size(f::Applied{<:Any,typeof(hvcat),<:Tuple{Int,Vararg{Any}}})
    n = f.args[1]
    sum(size.(f.args[2:n:end],1)),sum(size.(f.args[2:1+n],2))
end

function size(f::Applied{<:Any,typeof(hvcat),<:Tuple{NTuple{N,Int},Vararg{Any}}}) where N
    n = f.args[1]

    as = tuple(2, (2 .+ cumsum(most(n)))...)
    sum(size.(getindex.(Ref(f.args), as),1)),sum(size.(f.args[2:1+n[1]],2))
end

@inline hvcat_getindex(f, k, j::Integer) = hvcat_getindex_recursive(f, (k, j), f.args...)

_hvcat_size(A) = size(A)
_hvcat_size(A::Number) = (1,1)
_hvcat_size(A::AbstractVector) = (size(A,1),1)

@inline function hvcat_getindex_recursive(f, (k,j)::Tuple{Integer,Integer}, N::Int, A, args...)
    T = eltype(f)
    m,n = _hvcat_size(A)
    N ≤ 0 && throw(BoundsError(f, (k,j))) # ran out of arrays
    k ≤ m && j ≤ n && return convert(T, A[k, j])::T
    k ≤ m && return hvcat_getindex_recursive(f, (k, j - n), N-1, args...)
    hvcat_getindex_recursive(f, (k - m, j), N, args[N:end]...)
end

@inline function hvcat_getindex_recursive(f, (k,j)::Tuple{Integer,Integer}, N::NTuple{M,Int}, A, args...) where M
    T = eltype(f)
    m,n = _hvcat_size(A)
    k ≤ m && return hvcat_getindex_recursive(f, (k, j), N[1], A, args...)
    hvcat_getindex_recursive(f, (k - m, j), tail(N), args[N[1]:end]...)
end

@inline hvcat_getindex_recursive(f, idx, N) = throw(BoundsError(f, idx))

getindex(f::ApplyMatrix{<:Any,typeof(hvcat)}, k::Integer, j::Integer) = hvcat_getindex(f, k, j)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(hvcat)}, k::Integer, j::Integer)= hvcat_getindex(f, k, j)
getindex(f::Applied{<:Any,typeof(hvcat)}, k::Integer, j::Integer)= hvcat_getindex(f, k, j)

#####
# copyto!
####
# based on Base/array.jl, Base/abstractarray.jl
_copyto!(_, LAY::ApplyLayout{typeof(vcat)}, dest::AbstractArray{<:Any,N}, V::AbstractArray{<:Any,N}) where N =
    vcat_copyto!(dest, arguments(LAY, V)...)
function vcat_copyto!(dest::AbstractMatrix, arrays...)
    nargs = length(arrays)
    nrows = size(dest,1)
    nrows == sum(a->size(a, 1), arrays) || throw(DimensionMismatch("sum of rows each matrix must equal $nrows"))
    ncols = size(dest, 2)
    for a in arrays
        if size(a, 2) != ncols
            throw(DimensionMismatch("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
        end
    end
    pos = 1
    for a in arrays
        p1 = pos+size(a,1)-1
        copyto!(view(dest,pos:p1, :), a)
        pos = p1+1
    end
    return dest
end

function vcat_copyto!(arr::AbstractVector, arrays...)
    n = 0
    for a in arrays
        n += length(a)
    end
    n == length(arr) || throw(DimensionMismatch("destination must have length equal to sums of concatenated vectors"))

    i = 0
    @inbounds for a in arrays
        m = length(a)
        copyto!(view(arr,i+1:i+m), a)
        i += m
    end
    arr
end

function vcat_copyto!(arr::Vector{T}, arrays::Vector{T}...) where T
    n = 0
    for a in arrays
        n += length(a)
    end
    n == length(arr) || throw(DimensionMismatch("destination must have length equal to sums of concatenated vectors"))

    ptr = pointer(arr)
    if isbitstype(T)
        elsz = Core.sizeof(T)
    elseif isbitsunion(T)
        elsz = bitsunionsize(T)
        selptr = convert(Ptr{UInt8}, ptr) + n * elsz
    else
        elsz = Core.sizeof(Ptr{Cvoid})
    end
    t = @_gc_preserve_begin arr
    for a in arrays
        na = length(a)
        nba = na * elsz
        if isbitstype(T)
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr, a, nba)
        elseif isbitsunion(T)
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr, a, nba)
            # copy selector bytes
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  selptr, convert(Ptr{UInt8}, pointer(a)) + nba, na)
            selptr += na
        else
            ccall(:jl_array_ptr_copy, Cvoid, (Any, Ptr{Cvoid}, Any, Ptr{Cvoid}, Int),
                  arr, ptr, a, pointer(a), na)
        end
        ptr += nba
    end
    @_gc_preserve_end t
    return arr
end

# special case for adjoints of hcat. This is useful for catching fast paths
# for vector case, e.g., _fast_blockbradcast_copyto! in BlockArrays.jl
function vcat_copyto!(dest::AbstractMatrix, arrays::Adjoint{<:Any,<:AbstractVector}...)
    hcat_copyto!(dest', map(adjoint, arrays)...)
    dest
end

_copyto!(_, LAY::ApplyLayout{typeof(hcat)}, dest::AbstractMatrix, H::AbstractMatrix) =
    hcat_copyto!(dest, arguments(LAY,H)...)
function hcat_copyto!(dest::AbstractMatrix, arrays...)
    nargs = length(arrays)
    nrows = size(dest, 1)
    ncols = 0
    dense = true
    for a in arrays
        dense &= isa(a,Array)
        nd = ndims(a)
        ncols += (nd==2 ? size(a,2) : 1)
    end

    nrows == size(first(arrays),1) || throw(DimensionMismatch("Destination rows must match"))
    ncols == size(dest,2) || throw(DimensionMismatch("Destination columns must match"))

    pos = 1
    if dense
        for a in arrays
            n = length(a)
            copyto!(dest, pos, a, 1, n)
            pos += n
        end
    else
        for a in arrays
            p1 = pos+(isa(a,AbstractMatrix) ? size(a, 2) : 1)-1
            copyto!(view(dest,:, pos:p1), a)
            pos = p1+1
        end
    end
    return dest
end

function hcat_copyto!(dest::AbstractMatrix, arrays::AbstractVector...)
    height = size(dest, 1)
    for j = 1:length(arrays)
        if length(arrays[j]) != height
            throw(DimensionMismatch("vectors must have same lengths"))
        end
    end
    for j=1:length(arrays)
        copyto!(view(dest,:,j), arrays[j])
    end

    dest
end

_copyto!(_, lay::ApplyLayout{typeof(hvcat)}, dest::AbstractMatrix, src::AbstractMatrix) = hvcat_copyto!(dest, arguments(lay, src)...)

function hvcat_copyto!(out::AbstractMatrix{T}, nbc::Integer, as...) where T
    # nbc = # of block columns
    n = length(as)
    mod(n,nbc) != 0 &&
        throw(ArgumentError("number of arrays $n is not a multiple of the requested number of block columns $nbc"))
    nbr = div(n,nbc)
    hvcat_copyto!(out, ntuple(i->nbc, nbr), as...)
end

function hvcat_copyto!(a::AbstractMatrix{T}, rows::Tuple{Vararg{Int}}, xs::T...) where T<:Number
    nr = length(rows)
    nc = rows[1]

    size(a) == (nc,nr) || throw(DimensionMismatch())

    if length(a) != length(xs)
        throw(ArgumentError("argument count does not match specified shape (expected $(length(a)), got $(length(xs)))"))
    end
    k = 1
    @inbounds for i=1:nr
        if nc != rows[i]
            throw(ArgumentError("row $(i) has mismatched number of columns (expected $nc, got $(rows[i]))"))
        end
        for j=1:nc
            a[i,j] = xs[k]
            k += 1
        end
    end
    a
end

function hvcat_copyto!(out::AbstractMatrix{T}, rows::Tuple{Vararg{Int}}, as::AbstractVecOrMat...) where T
    nbr = length(rows)  # number of block rows

    nc = 0
    for i=1:rows[1]
        nc += size(as[i],2)
    end

    nr = 0
    a = 1
    for i = 1:nbr
        nr += size(as[a],1)
        a += rows[i]
    end

    size(out) == (nr,nc) || throw(DimensionMismatch())

    a = 1
    r = 1
    for i = 1:nbr
        c = 1
        szi = size(as[a],1)
        for j = 1:rows[i]
            Aj = as[a+j-1]
            szj = size(Aj,2)
            if size(Aj,1) != szi
                throw(ArgumentError("mismatched height in block row $(i) (expected $szi, got $(size(Aj,1)))"))
            end
            if c-1+szj > nc
                throw(ArgumentError("block row $(i) has mismatched number of columns (expected $nc, got $(c-1+szj))"))
            end
            out[r:r-1+szi, c:c-1+szj] = Aj
            c += szj
        end
        if c != nc+1
            throw(ArgumentError("block row $(i) has mismatched number of columns (expected $nc, got $(c-1))"))
        end
        r += szi
        a += rows[i]
    end
    out
end

#####
# adjoint/transpose
#####

for adj in (:adjoint, :transpose)
    @eval begin
        $adj(A::Hcat{T}) where T = Vcat{T}(map($adj,A.args)...)
        $adj(A::Vcat{T,2}) where T = Hcat{T}(map($adj,A.args)...)
    end
end

_vec(a) = a
_vec(a::AbstractArray) = vec(a)
_vec(a::Adjoint{<:Number,<:AbstractVector}) = _vec(parent(a))
vec(A::Hcat) = Vcat(map(_vec,A.args)...)

copy(f::Adjoint{<:Any,<:Union{Vcat,Hcat}}) = copy(parent(f))'
copy(f::Transpose{<:Any,<:Union{Vcat,Hcat}}) = transpose(copy(parent(f)))

_permutedims(a) = a
_permutedims(a::AbstractArray) = permutedims(a)

permutedims(A::Hcat{T}) where T = Vcat{T}(map(_permutedims,A.args)...)
permutedims(A::Vcat{T}) where T = Hcat{T}(map(_permutedims,A.args)...)

transposelayout(::ApplyLayout{typeof(vcat)}) = ApplyLayout{typeof(hcat)}()
transposelayout(::ApplyLayout{typeof(hcat)}) = ApplyLayout{typeof(vcat)}()

arguments(::ApplyLayout{typeof(vcat)}, A::Adjoint) = map(adjoint, arguments(ApplyLayout{typeof(hcat)}(), parent(A)))
arguments(::ApplyLayout{typeof(hcat)}, A::Adjoint) = map(adjoint, arguments(ApplyLayout{typeof(vcat)}(), parent(A)))
arguments(::ApplyLayout{typeof(vcat)}, A::Transpose) = map(transpose, arguments(ApplyLayout{typeof(hcat)}(), parent(A)))
arguments(::ApplyLayout{typeof(hcat)}, A::Transpose) = map(transpose, arguments(ApplyLayout{typeof(vcat)}(), parent(A)))

function arguments(::ApplyLayout{typeof(vcat)}, C::CachedVector)
    data = cacheddata(C)
    Vcat(data, C.array[length(data)+1:end])
end

#####
# broadcasting
#
# We want broadcasting for numbers with concaenations to pass through
# to take advantage of special implementations of the sub-components
######

BroadcastStyle(::Type{<:Vcat{<:Any,N}}) where N = LazyArrayStyle{N}()
BroadcastStyle(::Type{<:Hcat{<:Any}}) where N = LazyArrayStyle{2}()

# This is if we broadcast a function on a mixed concat f.([1; [2,3]])
# such that f returns a vector, e.g., f(1) == [1,2], we don't want
# to have the concat return [f(1); [f(2),f(3)]] but rather [[f(1)]; [f(2),f(3)]]

_flatten_nums(args::Tuple{}, bc::Tuple{}) = ()
_flatten_nums(args::Tuple, bc::Tuple) = (bc[1], _flatten_nums(tail(args), tail(bc))...)
_flatten_nums(args::Tuple{Number, Vararg{Any}}, bc::Tuple{AbstractArray, Vararg{Any}}) = (Fill(bc[1],1), _flatten_nums(tail(args), tail(bc))...)

broadcasted(::LazyArrayStyle, op, A::Vcat) = Vcat(_flatten_nums(A.args, broadcast(x -> broadcast(op, x), A.args))...)

for Cat in (:Vcat, :Hcat)
     @eval begin
        broadcasted(::LazyArrayStyle, op, A::$Cat, c::Number) = $Cat(_flatten_nums(A.args, broadcast((x,y) -> broadcast(op, x, y), A.args, c))...)
        broadcasted(::LazyArrayStyle, op, c::Number, A::$Cat) = $Cat(_flatten_nums(A.args, broadcast((x,y) -> broadcast(op, x, y), c, A.args))...)
        broadcasted(::LazyArrayStyle, op, A::$Cat, c::Ref) = $Cat(_flatten_nums(A.args, broadcast((x,y) -> broadcast(op, x, Ref(y)), A.args, c))...)
        broadcasted(::LazyArrayStyle, op, c::Ref, A::$Cat) = $Cat(_flatten_nums(A.args, broadcast((x,y) -> broadcast(op, Ref(x), y), c, A.args))...)
     end
 end

# determine indices of components of a vcat
_vcat_axes(::Tuple{}) = (1,)
_vcat_axes(a::Tuple{<:AbstractUnitRange}) = (first(a),)
_vcat_axes(::Tuple{}, b, c...) = tuple(1, broadcast(x -> broadcast(+, 1, x), _vcat_axes(b, c...))...)
_vcat_axes(a::Tuple{<:AbstractUnitRange}, b, c...) = tuple(first(a), broadcast((α,x) -> broadcast(+, α, x), last(first(a)),
                                                            _vcat_axes(b, c...))...)

_vcat_getindex_eval(y) = ()
_vcat_getindex_eval(y, a, b...) = tuple(y[a], _vcat_getindex_eval(y, b...)...)

# let it stay lazy
layout_broadcasted(::ApplyLayout{typeof(vcat)}, ::AbstractLazyLayout, op, A::AbstractVector, B::AbstractVector) =
    Broadcasted{LazyArrayStyle{1}}(op, (A, B))
layout_broadcasted(::AbstractLazyLayout, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) =
    Broadcasted{LazyArrayStyle{1}}(op, (A, B))

layout_broadcasted(::ApplyLayout{typeof(vcat)}, lay::CachedLayout, op, A::AbstractVector, B::AbstractVector) = layout_broadcasted(UnknownLayout(), lay, op, A, B)
layout_broadcasted(lay::CachedLayout, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) = layout_broadcasted(lay, UnknownLayout(), op, A, B)

function layout_broadcasted(::ApplyLayout{typeof(vcat)}, _, op, A::AbstractVector, B::AbstractVector)
    kr = _vcat_axes(map(axes,A.args)...)  # determine how to break up B
    B_arrays = _vcat_getindex_eval(B,kr...)    # evaluate B at same chunks as A
    ApplyVector(vcat, broadcast((a,b) -> broadcast(op,a,b), A.args, B_arrays)...)
end

function layout_broadcasted(_, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector)
    kr = _vcat_axes(axes.(B.args)...)
    A_arrays = _vcat_getindex_eval(A,kr...)
    Vcat(broadcast((a,b) -> broadcast(op,a,b), A_arrays, B.args)...)
end

######
# Special Vcat broadcasts
#
# We use Vcat for infinite padded vectors, so we need to special case
# two arrays. This may be generalisable in the future
######

layout_broadcasted(lay::ApplyLayout{typeof(vcat)}, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) =
    _vcat_layout_broadcasted(arguments(lay, A), arguments(lay, B), op, A, B)

_vcat_layout_broadcasted(Aargs, Bargs, op, A, B) = Broadcasted{LazyArrayStyle{1}}(op, (A,B))

function _vcat_layout_broadcasted((Ahead,Atail)::Tuple{AbstractVector,Any}, (Bhead,Btail)::Tuple{AbstractVector,Any}, op, A, B)
    T = Broadcast.combine_eltypes(op, (eltype(A), eltype(B)))

    if length(Ahead) ≥ length(Bhead)
        M,m = length(Ahead), length(Bhead)
        Chead = Vector{T}(undef,M)
        view(Chead,1:m) .= op.(view(Ahead,1:m), Bhead)
        view(Chead,m+1:M) .= op.(view(Ahead,m+1:M),Btail[1:M-m])

        Ctail = op.(Atail, Btail[M-m+1:end])
    else
        m,M = length(Ahead), length(Bhead)
        Chead = Vector{T}(undef,M)
        view(Chead,1:m) .= op.(Ahead, view(Bhead,1:m))
        view(Chead,m+1:M) .= op.(Atail[1:M-m],view(Bhead,m+1:M))

        Ctail = op.(Atail[M-m+1:end], Btail)
    end

    Vcat(Chead, Ctail)
end

function _vcat_layout_broadcasted((Ahead,Atail)::Tuple{AbstractVector,Any}, (Bhead,Btail)::Tuple{Number,Any}, op, A, B)
    T = Broadcast.combine_eltypes(op, (eltype(A), eltype(B)))
    M = length(Ahead)
    Chead = Vector{T}(undef,max(1,M))
    Chead[1] = op(A[1], Bhead)
    view(Chead,2:M) .= op.(view(Ahead,2:M),Btail[1:M-1])
    Ctail = op.(Atail, Btail[M:end])
    Vcat(Chead, Ctail)
end

function _vcat_layout_broadcasted((Ahead,Atail)::Tuple{Number,Any}, (Bhead,Btail)::Tuple{AbstractVector,Any}, op, A, B)
    T = Broadcast.combine_eltypes(op, (eltype(A), eltype(B)))
    M = length(Bhead)
    Chead = Vector{T}(undef,max(1,M))
    Chead[1] = op(Ahead, B[1])
    view(Chead,2:M) .= op.(Atail[1:M-1],view(Bhead,2:M))
    Ctail = op.(Atail[M:end],Btail)
    Vcat(Chead, Ctail)
end

_vcat_layout_broadcasted((Ahead,Atail)::Tuple{Number,Any}, (Bhead,Btail)::Tuple{Number,Any}, op, A, B) = Vcat(op.(Ahead,Bhead), op.(Atail,Btail))
_vcat_layout_broadcasted((Ahead,Atail)::Tuple{SVector{M},Any}, (Bhead,Btail)::Tuple{SVector{M},Any}, op, A, B) where M = Vcat(op.(Ahead,Bhead), op.(Atail,Btail))

broadcasted(::LazyArrayStyle, op, a::Vcat, b::AbstractArray) = layout_broadcasted(op, a, b)
broadcasted(::LazyArrayStyle, op, a::AbstractArray, b::Vcat) = layout_broadcasted(op, a, b)
broadcasted(::LazyArrayStyle{1}, op, a::Vcat{<:Any,1}, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), op, a, b)
broadcasted(::LazyArrayStyle{1}, op, a::Zeros{<:Any,1}, b::Vcat{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), op, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(\), a::Vcat{<:Any,1}, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), \, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(/), a::Zeros{<:Any,1}, b::Vcat{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), /, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::Vcat{<:Any,1}, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), *, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::Zeros{<:Any,1}, b::Vcat{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), *, a, b)

# Cannot broadcast Vcat's in a lazy way so stick to BroadcastArray
broadcasted(::LazyArrayStyle, op, A::Vcat, B::Vcat) = layout_broadcasted(op, A, B)

# ambiguities
broadcasted(::LazyArrayStyle, op, A::Vcat{<:Any,1}, B::CachedVector) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::CachedVector, B::Vcat{<:Any,1}) = layout_broadcasted(op, A, B)

function +(A::Vcat, B::Vcat)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end
function +(A::Vcat, B::AbstractArray)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end
function +(A::AbstractArray, B::Vcat)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end

####
# Cumsum
####

sum(V::Vcat) = mapreduce(sum, +, V.args)
all(V::Vcat) = all(all.(V.args))
any(V::Vcat) = any(any.(V.args))
all(f::Function, V::Vcat) = all(all.(f, V.args))
any(f::Function, V::Vcat) = any(any.(f, V.args))

_dotplus(a,b) = broadcast(+, a, b)

@inline _cumsum(x::Number) = x
@inline _cumsum(x) = cumsum(x)
_cumsum_last(x::AbstractVector{T}) where T = isempty(x) ? zero(T) : last(x)
_cumsum_last(x) = last(x)

_tuple_cumsum() = ()
_tuple_cumsum(a) = (a,)
_tuple_cumsum(a, b...) = (a, broadcast(+,a,_tuple_cumsum(b...))...)
function _vcat_cumsum(x...)
    cs = map(_cumsum,x)
    cslasts = tuple(0,_tuple_cumsum(map(_cumsum_last,most(cs))...)...)
    map((a,b) -> broadcast(+,a,b), cslasts, cs)
end

@inline cumsum(V::Vcat{<:Any,1}) = ApplyVector(vcat,_vcat_cumsum(V.args...)...)

_vcat_diff(x::Number) = ()
_vcat_diff(x) = (diff(x),)

_vcat_diff(a::Number, b, c...) = (first(b)-a, _vcat_diff(b,c...)...)
_vcat_diff(a, b, c...) = (diff(a), first(b)-last(a), _vcat_diff(b,c...)...)
@inline diff(V::Vcat{T,1}) where T = ApplyVector{T}(vcat,_vcat_diff(V.args...)...)

####
# maximum/minimum
####

sum(v::Vcat{<:Any,1}) = sum(map(sum,v.args))

for op in (:maximum, :minimum)
    @eval $op(V::Vcat) = $op($op.(V.args))
end

function in(x, V::Vcat)
    for a in V.args
        in(x, a) && return true
    end
    false
end

_fill!(a, x) = fill!(a,x)
function _fill!(a::Number, x)
    a == x || throw(ArgumentError("Cannot set $a to $x"))
    a
end

function fill!(V::Union{Vcat,Hcat}, x)
    for a in V.args
        _fill!(a, x)
    end
    V
end

###
# *
###

copy(M::Mul{ApplyLayout{typeof(vcat)},<:AbstractLazyLayout}) = vcat((arguments(vcat, M.A) .* Ref(M.B))...)

_all_prods(a::Tuple{}, ::Tuple) = ()
_all_prods(a::Tuple, b::Tuple) = tuple((Ref(first(a)) .* b)..., _all_prods(tail(a), b)...)
function copy(M::Mul{ApplyLayout{typeof(vcat)},ApplyLayout{typeof(hcat)}})
    b = arguments(hcat,M.B)
    ApplyArray(hvcat, length(b), _all_prods(arguments(vcat, M.A), b)...)
end

function materialize!(M::MatMulVecAdd{ApplyLayout{typeof(hcat)},ApplyLayout{typeof(vcat)}})
    α,A,B,β,C =  M.α,M.A,M.B,M.β,M.C
    T = eltype(C)
    _fill_lmul!(β,C) # this is temporary until strong β = false is supported
    for (a,b) in zip(A.args,B.args)
        materialize!(MulAdd(α,a,b,one(T),C))
    end
    C
 end

 function materialize!(M::MatMulMatAdd{ApplyLayout{typeof(hcat)},ApplyLayout{typeof(vcat)}})
    α,A,B,β,C =  M.α,M.A,M.B,M.β,M.C
    T = eltype(C)
    _fill_lmul!(β,C) # this is temporary until strong β = false is supported
    for (a,b) in zip(A.args,B.args)
        materialize!(MulAdd(α,a,b,one(T),C))
    end
    C
 end

 ####
 # col/rowsupport
 ####

most(a) = reverse(tail(reverse(a)))
function colsupport(lay::ApplyLayout{typeof(vcat)}, M::AbstractArray, j)
    args = arguments(lay, M)
    first(colsupport(first(args),j)):(size(Vcat(most(args)...),1)+last(colsupport(last(args),j)))
end

function rowsupport(lay::ApplyLayout{typeof(vcat)}, V::AbstractArray, k::Integer)
    ξ = k
    for A in arguments(lay, V)
        n = size(A,1)
        ξ ≤ n && return rowsupport(A, ξ)
        ξ -= n
    end
    return 1:0
end

function colsupport(lay::ApplyLayout{typeof(hcat)}, H::AbstractArray, j::Integer)
    ξ = j
    for A in arguments(lay,H)
        n = size(A,2)
        ξ ≤ n && return colsupport(A, ξ)
        ξ -= n
    end
    return 1:0
end

function rowsupport(lay::ApplyLayout{typeof(hcat)}, M::AbstractArray, k)
    args = arguments(lay, M)
    first(rowsupport(first(args),k)):(size(Hcat(most(args)...),2)+last(rowsupport(last(args),k)))
end

include("padded.jl")

###
# norm
###

for Cat in (:Vcat, :Hcat)
    for (op,p) in ((:norm1,1), (:norm2,2), (:normInf,Inf))
        @eval $op(a::$Cat) = $op(norm.(a.args,$p))
    end
    @eval normp(a::$Cat, p) = norm(norm.(a.args, p), p)
end

###
# subarrays
###

sublayout(::ApplyLayout{typeof(vcat)}, _) = ApplyLayout{typeof(vcat)}()
sublayout(::ApplyLayout{typeof(hcat)}, _) = ApplyLayout{typeof(hcat)}()
# a row-slice of an Hcat is equivalent to a Vcat
sublayout(::ApplyLayout{typeof(hcat)}, ::Type{<:Tuple{Number,AbstractVector}}) = ApplyLayout{typeof(vcat)}()

_vcat_lastinds(sz) = _vcat_cumsum(sz...)
_vcat_firstinds(sz) = (1, (1 .+ most(_vcat_lastinds(sz)))...)

_argsindices(sz) = broadcast(:, _vcat_firstinds(sz), _vcat_lastinds(sz))

_view_vcat(a::Number, kr) = Fill(a,length(kr))
_view_vcat(a::Number, kr, jr) = Fill(a,length(kr), length(jr))
_view_vcat(a, kr...) = _viewifmutable(a, kr...)

function _vcat_sub_arguments(::ApplyLayout{typeof(vcat)}, A, V, kr)
    sz = size.(arguments(A),1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    map(_view_vcat, arguments(A), skr2)
end

function _vcat_sub_arguments(::ApplyLayout{typeof(vcat)}, A, V, kr, jr)
    sz = size.(arguments(A),1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    _view_vcat.(arguments(A), skr2, Ref(jr))
end

_vcat_sub_arguments(LAY::ApplyLayout{typeof(vcat)}, A, V) = _vcat_sub_arguments(LAY, A, V, parentindices(V)...)

function _vcat_sub_arguments(L::ApplyLayout{typeof(hcat)}, A, V)
    A = parent(V)
    args = arguments(L, A)
    k,jr = parentindices(V)
    sz = size.(args,2)
    sjr = intersect.(_argsindices(sz), Ref(jr))
    sjr2 = broadcast((a,b) -> a .- b .+ 1, sjr, _vcat_firstinds(sz))
    _view_hcat.(args, k, sjr2)
end

_vcat_sub_arguments(::DualLayout{ML}, A, V) where ML = _vcat_sub_arguments(ML(), A, V)
_vcat_sub_arguments(A, V) = _vcat_sub_arguments(MemoryLayout(typeof(A)), A, V)
arguments(::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,1}) = _vcat_sub_arguments(parent(V), V)

function arguments(L::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2})
    A = parent(V)
    args = arguments(L, A)
    kr,jr = parentindices(V)
    sz = size.(args,1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    _view_vcat.(args, skr2, Ref(jr))
end

@inline _view_hcat(a::Number, kr, jr) = Fill(a,length(kr),length(jr))
@inline _view_hcat(a::Number, kr::Number, jr) = Fill(a,length(jr))
@inline _view_hcat(a, kr, jr) = _viewifmutable(a, kr, jr)
@inline _view_hcat(a::AbstractVector, kr, jr::Colon) = _viewifmutable(a, kr)

# equivalent to broadcast but written to maintain type stability
__view_hcat(::Tuple{}, _, ::Tuple{}) = ()
__view_hcat(::Tuple{}, _, ::Colon) = ()
@inline __view_hcat(args::Tuple, kr, jrs::Tuple) = (_view_hcat(args[1], kr, jrs[1]), __view_hcat(tail(args), kr, tail(jrs))...)
@inline __view_hcat(args::Tuple, kr, ::Colon) = (_view_hcat(args[1], kr, :), __view_hcat(tail(args), kr, :)...)

function arguments(L::ApplyLayout{typeof(hcat)}, V::SubArray)
    A = parent(V)
    args = arguments(L, A)
    kr,jr = parentindices(V)
    sz = size.(args,2)
    sjr = intersect.(_argsindices(sz), Ref(jr))
    sjr2 = broadcast((a,b) -> a .- b .+ 1, sjr, _vcat_firstinds(sz))
    __view_hcat(args, kr, sjr2)
end

arguments(::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{<:Slice,<:Any}}) =
    _viewifmutable.(arguments(parent(V)), Ref(:), Ref(parentindices(V)[2]))
arguments(::ApplyLayout{typeof(hcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{<:Any,<:Slice}}) =
    __view_hcat(arguments(parent(V)), parentindices(V)[1], :)

function sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, _)
    ret = similar(V)
    n = 0
    _,jr = parentindices(V)
    for a in arguments(V)
        m = size(a,1)
        copyto!(view(ret,n+1:n+m,:), a)
        n += m
    end
    ret
end

sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractVector, _) = ApplyVector(V)

function sub_materialize(::ApplyLayout{typeof(hcat)}, V, _)
    ret = similar(V)
    n = 0
    kr,_ = parentindices(V)
    for a in arguments(V)
        m = size(a,2)
        copyto!(view(ret,:,n+1:n+m), a)
        n += m
    end
    ret
end

# temporarily allocate. In the future, we add a loop over arguments
materialize!(M::MatMulMatAdd{<:AbstractColumnMajor,<:ApplyLayout{typeof(vcat)}}) =
    materialize!(MulAdd(M.α,M.A,Array(M.B),M.β,M.C))
materialize!(M::MatMulVecAdd{<:AbstractColumnMajor,<:ApplyLayout{typeof(vcat)}}) =
    materialize!(MulAdd(M.α,M.A,Array(M.B),M.β,M.C))

## print

_replace_in_print_matrix(A::AbstractArray, k, j, s) = replace_in_print_matrix(A, k, j, s)
_replace_in_print_matrix(_, k, j, s) = s

function layout_replace_in_print_matrix(LAY::ApplyLayout{typeof(vcat)}, f::AbstractVecOrMat, k, j, s)
    κ = k
    for A in arguments(LAY, f)
        n = size(A,1)
        κ ≤ n && return _replace_in_print_matrix(A, κ, j, s)
        κ -= n
    end
    throw(BoundsError(f, (k,j)))
end

# searchsorted

_searchsortedfirst(a, x) = searchsortedfirst(a, x)
_searchsortedfirst(a::Number, x) = 1 + (x > a)
_searchsortedlast(a, x) = searchsortedlast(a, x)
_searchsortedlast(a::Number, x) = 0 + (x ≥ a)

searchsortedfirst(f::Vcat{<:Any,1}, x) =
    searchsortedfirst_recursive(0, x, arguments(f)...)

searchsortedlast(f::Vcat{<:Any,1}, x) =
    searchsortedlast_recursive(length(f), x, reverse(arguments(f))...)

@inline searchsortedfirst_recursive(n, x) = n + 1

@inline function searchsortedfirst_recursive(n, x, a, args...)
    m = length(a)
    r = _searchsortedfirst(a, x)
    r ≤ m && return n + r
    return searchsortedfirst_recursive(n + m, x, args...)
end

@inline searchsortedlast_recursive(n, x) = n

@inline function searchsortedlast_recursive(n, x, a, args...)
    n -= length(a)
    r = _searchsortedlast(a, x)
    r > 0 && return n + r
    return searchsortedlast_recursive(n, x, args...)
end

searchsorted(f::Vcat{<:Any,1}, x) = searchsortedfirst(f, x):searchsortedlast(f,x)


text: qsort(L) = isempty(L) ? L : vcat(qsort(filter(x -> x < L[1], L[2:end])), L[1:1], qsort(filter(x -> x >= L[1], L[2:end])))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Clang.cindex
using Base.Test

top = cindex.parse_header("cxx/cbasic.h",
                          flags = TranslationUnit_Flags.DetailedPreprocessingRecord |
                          TranslationUnit_Flags.SkipFunctionBodies)

# Tokenizer

CONST = tokenize(cindex.search(
                    top,
                    x-> (isa(x, cindex.MacroDefinition) && name(x) == "CONST") )[1])
CONSTADD = tokenize(cindex.search(
                    top,
                    x-> (isa(x, cindex.MacroDefinition) && name(x) == "CONSTADD"))[1])

@test (typeof(CONST[1]) == cindex.Identifier && typeof(CONST[2]) == cindex.Literal)
@test (CONSTADD[1].text == "CONSTADD" && CONSTADD[2].text == "CONST" &&
       CONSTADD[3].text == "+" && CONSTADD[4].text == "2")

# Function arguments

func1 = cindex.search(top, "func1")[1]
@test cindex.getNumArgTypes(cu_type(func1)) == 4
func1_args = cindex.function_args(func1)
@test map(spelling, func1_args) == ASCIIString["a","b","c","d"]


text:
# s^n = sum(a.^n) + error
struct Solution{T<:Integer,N}
    s :: T
    n :: Int
    a :: NTuple{N,T}
    function Solution(s::T,n::Int,a,AZsPC::Bool=true) where {T<:Integer}
        N = length(a)
        @assert N<s
        if AZsPC
            @assert N>1
            @assert n>15
            @assert n<41
        end
        for x in a
            @assert x<s
        end
        sorted_a=sort([a...])
        @assert allunique(sorted_a)
        new{T,N}(s,n,Tuple(sorted_a))
    end
end
Base.show(io::IO, x::Solution) = print(io,string(x))
function Base.string(sol::Solution,witherror::Bool= false)
    to_string((sol.s,sol.n,sol.a),witherror)
end
function to_string(x,witherror::Bool=false)
    s,n,a = x
    out = string(s,"^",n,"=>{",a[1])
    for value in a[2:end]
        out = out*string(",",value)
    end
    out = out*"}"
    if witherror
        out = out*",e="*string(err(s,n,a))
    end
    out
end
function err(sol::Solution, acc=BigInt(0))
    err((sol.s,sol.n,sol.a),acc)
end
function err(x, acc=BigInt(0))
    s,n,a = x
    acc = acc + BigInt(s)^n
    for ak in a
        acc = acc - BigInt(ak)^n
    end
    acc
end


text:
include("DM1_2.jl")

function GRASP(
    cost, #Une array représentant les couts, de taille n
    matrix, #une matrice de taille m*n representant les contrantes
    n, #nombre de variables
    m, #nombbre de contraintes
    alpha # le alpha que nous allons utiliser
     )

     (desactive_condition, stop1, variables_actives, stop2, util, SOL) = initial(m,n)

    #Création de la solution
    while desactive_condition!=stop1 && variables_actives!=stop2
        #Calcul des utilitées
        util = Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2)
        #Choix du candidat
        PosCandidat::Int64 = choixcandidat(util, alpha)
        #On ajoute le candidat à la solution
        SOL[PosCandidat] = true
        #Desactive! le candidat selectionné et les variables des contraintes où il apparait
        Desactive!(PosCandidat, matrix, desactive_condition, m, variables_actives,n)
    end
    z = calcul_z(SOL,cost,n)
#    println("Valeur Solution GRASP: " , z)
    return(SOL, z, desactive_condition)
end

#===========================================================================#
function initial(m,n)
    #Creating a set of lines that will be avaluated
    desactive_condition= zeros(Bool, m) #si desactive_condition[i]=0 la ligne sera évalué
    stop1= ones(Bool,m) #quand atteint, toutes les conditions ont été évaluées
    #Creation the an array of activated variables
    variables_actives = ones(Bool, n) #si variables_actives[j]=1 la variable xi sera considéré
    stop2= zeros(Bool , n) #quand atteint, toutes les variables ont été traité
    #creating my utility array for my m variables
    util = zeros(Float64, n)
    #initialisation de la Solution
    SOL = zeros(Bool, n)
    return(desactive_condition, stop1, variables_actives, stop2, util, SOL)
end

# =========================================================================== #
#Détermine l'utilitée des élément via leur contribution divisé par leur nombre d'apparition dans les contraintes
function Utilite(cost, matrix, desactive_condition, n, m , variables_actives, stop1, stop2)
    util = zeros(Float64, n) #réinitialisation du vecteur

    for j=1:n # Pour chaque variable
        if variables_actives[j] == 1 #on vérifie si elle est active
                for i=1:m #Pour chaque contrainte
                    if desactive_condition[i]==0 #on vérifie si celle-ci est évaluée
                        K=matrix[i,j]
                        util[j]=util[j]+K #On compte le nombre d'apprationt de xj
                    end #fin if
                end #fin for
            #On calcul la contribution moyenne
            if util[j]!=false #On évide de divisé par 0
                util[j] =  cost[j] / util[j]
            else #Si la variable ne contribue pas, on la désactive
                variables_actives[j] = 0
            end #fin if
        end #fin if
    end #fin for
    return util
end #fin Utilite

# =========================================================================== #
function Desactive!(PosCandidat, matrix, desactive_condition, m , variables_actives , n) #On désactive les conditions où apparaissent le candidat
    for i=1:m #pour chaque contraintes
        if matrix[i,PosCandidat] == 1
            desactive_condition[i] = 1
            for j=1:n #On désactive les variables de la même contrainte
                if matrix[i,j] ==1
                    variables_actives[j] = 0
                end #fin if
            end #fin for des variables
        end #fin if
    end #fin for de contrainte
end #fin Desactive!

#=========================================================================#
function calcul_z(SOL,costs,n)
    z = 0
    for i in 1:n
        z += SOL[i]*costs[i]
    end
    return z
end #fin calcul_z

#==========================================#
function choixcandidat(util, alpha) #On dératemine la liste de candidati restreinte et en choisit un au hasard
    max = maximum(util)
    min = minimum(util)
    liste_candidat_restreint = Float64[]
    cpt::Int = 1
    for i in 1:length(util)
        if util[i] >= trunc(min+(alpha*(max-min))*1000)/1000 #évite les erreur d'apporximation quand alpha = 1, pourrait biaisé le résultat, mais ce n'est arrivé dans aucune instances
            push!(liste_candidat_restreint,i)
        end
    end
    prob = rand(1:length(liste_candidat_restreint)) #tirage aléatoire
    return(liste_candidat_restreint[prob])
end


text:
using Sockets
using SHA
username = "revox" # Replace this with your username

socket = Sockets.connect("51.15.127.80", 2811)
println("Connected to Duino-Coin server")

server_ver = String(read(socket, 3))
println("Server is on version: ", server_ver)

while true
	write(socket, string("JOB,", String(username), ",MEDIUM"))
	job = String(read(socket, 87))
	#println("Job received: ", job)

	job = split(job, ",")
	lastBlockHash = job[1]
	result = job[2]
	difficulty = parse(Int32, job[3]) * 100

	for i = 0:difficulty
		stringToHash = string(lastBlockHash, string.(i))
		ducos1 = bytes2hex(sha1(stringToHash))

		if ducos1 == result
			write(socket, string(i, ",,Julia Miner"))
			feedback = String(read(socket, 4))
			if feedback == "GOOD"
				println("Accepted share ", i, "\tDifficulty ", difficulty)
				break
			else
				println("Rejected share ", i, "\tDifficulty ", difficulty)
				break
			end
		end
		i += 1
	end
end

max_stars_repo_path: Unofficial miners/Julia_Miner.jl

max_stars_repo_name: BastelPichi/duino-coin

max_issues_repo_path: Unofficial miners/Julia_Miner.jl

max_issues_repo_name: BastelPichi/duino-coin

max_forks_repo_path: Unofficial miners/Julia_Miner.jl

max_forks_repo_name: BastelPichi/duino-coin


text:
module IOcdf

"""
    module EmpiricalCDFS.IOcdf

Binary IO functions for `EmpiricalCDFs`.

Types: `CDFfile`

Functions: `save`, `save`, `readcdf`, `readcdfinfo`, `header`, `version`, `getcdf`, `data`

`CDFfile` contains an `AbstractEmpiricalCDF` as a field. Many functions on the latter are forwarded to
the former.
"""
IOcdf

export CDFfile, save, readcdf, readcdfinfo, header, version, getcdf, data

import Statistics
using ..EmpiricalCDFs: AbstractEmpiricalCDF, EmpiricalCDF, EmpiricalCDFHi

import ..EmpiricalCDFs: data, linprint, logprint, getindex, getcdfindex, counts

include("readstring.jl")

#### CDFfile type

@enum CDFTYPE EmpiricalCDFtype=1 EmpiricalCDFHitype=2

const CDFfileVersion = v"0.0.2"

"""
    CDFfile(cdf::AbstractEmpiricalCDF, header="")

    struct CDFfile{T <: AbstractEmpiricalCDF}
        cdf::T
        header::String
    end

Binary data file for `AbstractEmpiricalCDF`

The file format is

- Identifying string
- `n::Int64` number of bytes in the header string
- `s::String` The header string
- `t::Int64` Type of `AbstractEmpiricalCDF`, 1 or 2. 1 for `EmpiricalCDF`, 2 for `EmpiricalCDFHi`.
- `lowreject::Float64` the lower cutoff, only for `EmpiricalCDFHi`.
- `npts::Int64` number of data points in the CDF
- `npts` data points of type `Float64`
"""
struct CDFfile{T <: AbstractEmpiricalCDF}
    cdf::T
    header::String
    vn::VersionNumber
end

CDFfile(cdf) = CDFfile(cdf,"",CDFfileVersion)

data(cdff::CDFfile) = data(cdff.cdf)

function Base.show(io::IO, cdff::CDFfile)
    print(io, string(typeof(cdff)), "(")
    show(io,cdff.cdf)
    print(io, ",")
    print(io, "v", cdff.vn)
    print(io, ")")
    nothing
end

"""
    header::String = header(cdff::CDFfile)

Return the header from `cdff`.
"""
header(cdff::CDFfile) = cdff.header

"""
    cdf::AbstractEmpiricalCDF = getcdf(cdff::CDFfile)

Return the CDF from `cdff`.
"""
getcdf(cdff::CDFfile) = cdff.cdf

"""
    version(cdff::CDFfile)

Return the version number of the file format.
"""
version(cdff::CDFfile) = cdff.vn

for f in (:sort!, :push!, :append!, :getindex, :length, :size, :firstindex, :lastindex, :rand, :minimum, :maximum)
    @eval Base.$(f)(cdff::CDFfile, args...) = $(f)(cdff.cdf, args...)
end

for f in (:mean, :std, :quantile )
    @eval Statistics.$(f)(cdff::CDFfile, args...) = $(f)(cdff.cdf, args...)
end

for f in (:getcdfindex, :counts)
#           :mle, :KSstatistic, :mleKS, :scanmle)
    @eval $(f)(cdff::CDFfile, args...) = $(f)(cdff.cdf, args...)
end

for f in (:linprint, :logprint )
    @eval $(f)(ioorfile, cdff::CDFfile, args...) = $(f)(ioorfile, cdff.cdf, args...)
end

(cdff::CDFfile)(x::Real) = getcdf(cdff)(x)
(cdff::CDFfile)(v::AbstractArray) = getcdf(cdff)(v)

####

function make_CDFfile_version_string(v)
    nchars = 100
    s = "CDFfile " * string(v)
    s = s * " "^(nchars-length(s))
end

make_CDFfile_version_string() = make_CDFfile_version_string(CDFfileVersion)

function read_CDFfile_version_string(io::IO)
    nchars = 100
    local s
    try
        s = readstring(io,nchars)
    catch
        error("Failed reading identifying string and version number from CDFfile")
    end
    (filetype,versionstring) = split(s)
    if filetype != "CDFfile"
        error("Failed reading version number from CDFfile")
    end
    local vn
    try
        vn = VersionNumber(versionstring)
    catch
        error("Unable to parse version number string '$versionstring'")
    end
    vn
end

#### Writing CDFfile and CDF

function save(io::IO,cdff::CDFfile)
    write(io, make_CDFfile_version_string())
    write(io,sizeof(cdff.header))
    write(io,cdff.header)
    save(io,cdff.cdf)
end

function save(io::IO,cdf::EmpiricalCDF)
    write(io, Int64(EmpiricalCDFtype))
    _writedata(io,cdf)
end

function save(io::IO,cdf::EmpiricalCDFHi)
    write(io, Int64(EmpiricalCDFHitype))
    write(io, convert(Float64, cdf.lowreject))
    _writedata(io, cdf)
end

function _writedata(io::IO,cdf::AbstractEmpiricalCDF)
    write(io, convert(Int64, length(cdf)))
    for x in cdf.xdata
        write(io,x)
    end
end

"""
    save(fn::String, cdf::AbstractEmpiricalCDF, header::String="")

write `cdf` to file `fn` in a fast binary format.
"""
function save(fn::String, cdf::AbstractEmpiricalCDF, header::String="")
    io = open(fn,"w")
    save(io, CDFfile(cdf,header,CDFfileVersion))
    close(io)
end

#save(fn::String, cdf::AbstractEmpiricalCDF) = save(fn,cdf,"")

#### Reading CDFfile and CDF

function readcdfdata(io::IO, cdf::AbstractEmpiricalCDF)
    npts = read(io,Int64)
    resize!(cdf.xdata,npts)
@inbounds for i in 1:npts
        x = read(io,Float64)
        cdf.xdata[i] = x
    end
end

function readcdf(io::IO)
    info = readcdfinfo(io)
    (vn,header) = (info.vn,info.header)
    cdf = typeof(info) == CDFInfo ?  EmpiricalCDF() : EmpiricalCDFHi(info.lowreject)
    readcdfdata(io,cdf)
    CDFfile(cdf,header,vn)
end

struct CDFInfo
    vn::VersionNumber
    header::String
    npts::Int
end

function Base.show(io::IO, i::CDFInfo)
    print(io,typeof(i),"(")
    print(io,"length(header)=", length(i.header))
    print(io,",vn=v", i.vn)
    print(io,",npts=",i.npts,")")
end

struct CDFHiInfo
    vn::VersionNumber
    header::String
    npts::Int
    lowreject::Float64
end

function Base.show(io::IO, i::CDFHiInfo)
    print(io,typeof(i),"(")
    print(io,"length(header)=", length(i.header))
    print(io,",vn=v", i.vn)
    print(io,",npts=",i.npts)
    print(io,",lowreject=",i.lowreject,")")
end

function readcdfinfo(io::IO)
    vn = read_CDFfile_version_string(io)
    header = readlengthandstring(io)
    cdftype = CDFTYPE(read(io,Int64))
    local lowreject
    if cdftype == EmpiricalCDFtype
        npts = peektype(io,Int64)
        CDFInfo(vn,header,npts)
    elseif cdftype == EmpiricalCDFHitype
        lowreject = read(io,Float64)
        npts = peektype(io,Int64)
        CDFHiInfo(vn,header,npts,lowreject)
    else
        error("Uknown cdf type ", cdftype)
    end
end

"""
    readcdfinfo(fn::String)

Return an object containing information about the cdf saved in the binary
file `fn`. The data itself is not read.
"""
function readcdfinfo(fn::String)
    io = open(fn,"r")
    info = readcdfinfo(io)
    close(io)
    info
end

"""
    readcdf(fn::String)

Read an empirical CDF from file `fn`. Return an object
`cdff` of type `CDFfile`. The header is in field `header`.
The cdf is in in field `cdf`.
"""
function readcdf(fn::String)
    io = open(fn,"r")
    cdffile = readcdf(io)
    close(io)
    cdffile
end

end # module IO


text:
using Random, Statistics, DataFrames, Distributed

# Run multiple instances of a learner
# (type and parameter values defined in 'settings')
# on Categorical environments (Section 2.3.2)
# (defined in 'specs')
function comparison(settings, specs::SpecsMultinomial)

    df = @distributed vcat for setting in settings
        parameterset = getparamset(setting)
        # @show parameterset
        callback = run!(specs.envtype, setting[1][1].type, parameterset,
                        ns = specs.ns,
                        stochasticity = setting[1][2],
                        changeprobability = setting[1][3],
                        nsteps = specs.nsteps,
                        iseedenvpolicy = specs.seedenvpolicy[setting[2]],
                        iseedlearner = specs.seedlearner[setting[2]]
                        )
        terror = computesquareddiff(callback.true_history,
                                    callback.estimated_history)
        DataFrame(learner = setting[1][1].name,
                type = setting[1][1].type,
                param = parameterset,
                stochasticity = setting[1][2],
                changeprobability = setting[1][3],
                seed = specs.seedenvpolicy[setting[2]],
                seedlearner = specs.seedlearner[setting[2]],
                terrors = [terror],
                states = [callback.state_history],
                switches = [callback.switchflag])#,
                #nposteriors = [callback.nposteriors[2:end]])
        end
    df
end


text:
abstract AbstractNoiseModel

noise(s::Symbol,args...) = _noise(Val{s},args...)
noisemodelname(::AbstractNoiseModel) = "AbstractNoiseModel"

function show(io::IO,n::AbstractNoiseModel)
  println(io,noisemodelname(n))
  for f in fieldnames(n)
    print(io," ",f,": ",typeof(getfield(n,f)))
    println(io)
  end
  nothing
end


text:
"""
This package provides some tools for working with ColorSchemes:

    colorscheme_to_image(), colorscheme_to_text(), colorscheme_weighted(),
    compare_colors(), extract(), extract_weighted_colors(), convert_to_scheme(),
    image_to_swatch(), sortcolorscheme(), get_linear_segment_color(),
    make_colorscheme(), get_indexed_list_color()
"""
module ColorSchemeTools

using Images, ColorSchemes, Colors, Clustering, FileIO, Dates

import Base.get

export
    colorscheme_to_image,
    colorscheme_to_text,
    colorscheme_weighted,
    compare_colors,
    extract,
    extract_weighted_colors,
    convert_to_scheme,
    image_to_swatch,
    sortcolorscheme,
    make_colorscheme,
    get_linear_segment_color,
    get_indexed_list_color

"""
    extract(imfile, n=10, i=10, tolerance=0.01; shrink=n)

`extract()` extracts the most common colors from an image from the image file
`imfile` by finding `n` dominant colors, using `i` iterations. You can (and
probably should) shrink larger images before running this function.

Returns a ColorScheme.
"""
function extract(imfile, n=10, i=10, tolerance=0.01; kwargs...)
    ewc = extract_weighted_colors(imfile, n, i, tolerance; kwargs...)[1] # throw away the weights
    return ewc
end

"""
    extract_weighted_colors(imfile, n=10, i=10, tolerance=0.01; shrink = 2)

Extract colors and weights of the clusters of colors in an image file. Returns a
ColorScheme and weights.

Example:

```
pal, wts = extract_weighted_colors(imfile, n, i, tolerance; shrink = 2)
```
"""
function extract_weighted_colors(imfile, n=10, i=10, tolerance=0.01; shrink = 2.0)
    img = load(imfile)
    # TODO this is the wrong to do errors
    (!@isdefined img) && error("Can't load the image file \"$imfile\"")
    w, h = size(img)
    neww = round(Int, w/shrink)
    newh = round(Int, h/shrink)
    smaller_image = Images.imresize(img, (neww, newh))
    w, h = size(smaller_image)
    imdata = convert(Array{Float64}, channelview(smaller_image))

    !any(n -> n == 3, size(imdata)) && error("Image file \"$imfile\" doesn't have three color channels; perhaps it has an alpha channel as well?")

    d = reshape(imdata, 3, :)
    R = kmeans(d, n, maxiter=i, tol=tolerance)
    cols = RGB{Float64}[]
    for i in 1:3:length(R.centers)
        push!(cols, RGB(R.centers[i], R.centers[i+1], R.centers[i+2]))
    end
    return ColorScheme(cols), R.cweights/sum(R.cweights)
end

"""
    colorscheme_weighted(colorscheme, weights, length)

Returns a new ColorScheme of length `length` (default 50) where the proportion
of each color in `colorscheme` is represented by the associated weight of each
entry.

Examples:

```
colorscheme_weighted(extract_weighted_colors("hokusai.jpg")...)
colorscheme_weighted(extract_weighted_colors("filename00000001.jpg")..., 500)
```
"""
function colorscheme_weighted(cscheme::ColorScheme, weights, l = 50)
    iweights = map(n -> convert(Integer, round(n * l)), weights)
    #   adjust highest or lowest so that length of result is exact
    while sum(iweights) < l
        val, ix = findmin(iweights)
        iweights[ix]=val+1
    end
    while sum(iweights) > l
        val,ix = findmax(iweights)
        iweights[ix]=val-1
    end
    a = Array{RGB{Float64}}(undef, 0)
    for n in 1:length(cscheme)
        a = vcat(a, repeat([cscheme[n]], iweights[n]))
    end
    return ColorScheme(a)
end

"""
    compare_colors(color_a, color_b, field = :l)

Compare two colors, using the Luv colorspace. `field` defaults to luminance `:l`
but could be `:u` or `:v`. Return true if the specified field of `color_a` is
less than `color_b`.
"""
function compare_colors(color_a, color_b, field = :l)
    if 1 < color_a.r < 255
        fac = 255
    else
        fac = 1
    end
    luv1 = convert(Luv, RGB(color_a.r/fac, color_a.g/fac, color_a.b/fac))
    luv2 = convert(Luv, RGB(color_b.r/fac, color_b.g/fac, color_b.b/fac))
    return getfield(luv1, field) < getfield(luv2, field)
end

"""
    sortcolorscheme(colorscheme::ColorScheme, field; kwargs...)

Sort (non-destructively) a colorscheme using a field of the LUV colorspace.

The less than function is `lt = (x,y) -> compare_colors(x, y, field)`.

The default is to sort by the luminance field `:l` but could be by `:u` or `:v`.

Returns a new ColorScheme.
"""
function sortcolorscheme(colorscheme::ColorScheme, field = :l; kwargs...)
    cols = sort(colorscheme.colors, lt = (x,y) -> compare_colors(x, y, field); kwargs...)
    return ColorScheme(cols)
end

"""
    convert_to_scheme(cscheme, img)

Converts `img` from its current color values to use only the colors defined in the ColorScheme `cscheme`.

```
image = nonTransparentImg
convert_to_scheme(ColorSchemes.leonardo, image)
convert_to_scheme(ColorSchemes.Paired_12, image)
```
"""
convert_to_scheme(cscheme::ColorScheme,img) =
    map(c->get(cscheme, getinverse(cscheme, c)), img)

"""
    colorscheme_to_image(cs, nrows=50, tilewidth=5)

Make an image from a ColorScheme by repeating the colors in `nrows` rows, repeating each pixel `tilewidth` times.

Returns the image as an array.

Examples:

```
using FileIO

img = colorscheme_to_image(ColorSchemes.leonardo, 50, 200)
save("/tmp/cs_image.png", img)

save("/tmp/blackbody.png", colorscheme_to_image(ColorSchemes.blackbody, 10, 100))
```
"""
function colorscheme_to_image(cs::ColorScheme, nrows=50, tilewidth=5)
    ncols = tilewidth * length(cs)
    a = Array{RGB{Float64}}(undef, nrows, ncols)
    for row in 1:nrows
        for col in 1:ncols
            a[row, col] = cs.colors[div(col-1, tilewidth) + 1]
        end
    end
    return a
end

"""
    image_to_swatch(imagefilepath, samples, destinationpath;
        nrows=50,
        tilewidth=5)

Extract a ColorsSheme from the image in `imagefilepath` to a swatch image PNG in
`destinationpath`. This just runs `sortcolorscheme()`, `colorscheme_to_image()`,
and `save()` in sequence.

Specify the number of colors. You can also specify the number of rows, and how
many times each color is repeated.

```
image_to_swatch("monalisa.jpg", 10, "/tmp/monalisaswatch.png")
```
"""
function image_to_swatch(imagefilepath, n::Int64, destinationpath;
    nrows=50,
    tilewidth=5)
    tempcs = sortcolorscheme(extract(imagefilepath, n))
    img = colorscheme_to_image(tempcs, nrows, tilewidth)
    save(destinationpath, img)
end

"""
    colorscheme_to_text(cscheme::ColorScheme, schemename, filename;
        category="dutch painters",   # category
        notes="it's not really lost" # notes
    )

Write a ColorScheme to a Julia text file.

## Example

```
colorscheme_to_text(ColorSchemes.vermeer,
    "the_lost_vermeer",          # name
    "/tmp/the_lost_vermeer.jl",  # file
    category="dutch painters",   # category
    notes="it's not really lost" # notes
    )
```

and read it back in with:

```
include("/tmp/the_lost_vermeer.jl")
```
"""
function colorscheme_to_text(cs::ColorScheme, schemename::String, file::String;
        category="",
        notes="")
    fhandle = open(file, "w")
    write(fhandle, string("loadcolorscheme(:$(schemename), [\n"))
    for c in cs.colors
        write(fhandle, string("\tColors.$(c), \n"))
    end
    write(fhandle, string("], \"$(category)\", \"$(notes)\")"))
    close(fhandle)
end

"""
    get_linear_segment_color(dict, n)

Get the RGB color for value `n` from a dictionary of linear color segments.

This following is a dictionary where red increases from 0 to 1 over the bottom
half, green does the same over the middle half, and blue over the top half:

```
cdict = Dict(:red  => ((0.0,  0.0,  0.0),
                       (0.5,  1.0,  1.0),
                       (1.0,  1.0,  1.0)),
            :green => ((0.0,  0.0,  0.0),
                       (0.25, 0.0,  0.0),
                       (0.75, 1.0,  1.0),
                       (1.0,  1.0,  1.0)),
            :blue =>  ((0.0,  0.0,  0.0),
                       (0.5,  0.0,  0.0),
                       (1.0,  1.0,  1.0)))
```

The value of RGB component at every value of `n` is defined by a set of tuples.
In each tuple, the first number is `x`. Colors are linearly interpolated in
bands between consecutive values of `x`; if the first tuple is given by `(Z, A,
B)` and the second tuple by `(X, C, D)`, the color of a point `n` between Z and
X will be given by `(n - Z) / (X - Z) * (C - B) + B`.

For example, given an entry like this:

```
:red  => ((0.0, 0.0, 0.0),
          (0.5, 1.0, 1.0),
          (1.0, 1.0, 1.0))
```

and if `n` = 0.75, we return 1.0; 0.75 is between the second and third segments,
but we'd already reached 1.0 (segment 2) when `n` was 0.5.

"""
function get_linear_segment_color(dict, n)
    result = Float64[]
    for c in [:red, :green, :blue]
        listoftuples = dict[c]
        n = clamp(n, 0.0, 1.0)
        upper = max(2, findfirst(f -> n <= first(f), listoftuples))
        lower = max(1, upper - 1)
        lowersegment = listoftuples[lower]
        uppersegment = listoftuples[upper]
        Z, A, B = lowersegment
        X, C, D = uppersegment
        if X > Z
            color_at_n = (n - Z)/(X - Z) * (C - B) + B
        else
            color_at_n = 0.0
        end
        push!(result, color_at_n)
    end
    return result
end

"""
    lerp((x, from_min, from_max, to_min=0.0, to_max=1.0)

Linear interpolation of `x` between `from_min` and `from_max`.

Example

```
ColorSchemeTools.lerp(128, 0, 256)
0.5
```
"""
function lerp(x, from_min, from_max, to_min=0.0, to_max=1.0)
   if !isapprox(from_max, from_min)
       return ((x - from_min) / (from_max - from_min)) * (to_max - to_min) + to_min
   else
       return from_max
   end
end

"""
    get_indexed_list_color(indexedlist, n)

Get the color at a point `n` given an indexed list of triples like this:

```
gist_rainbow = (
       (0.000,   (1.00, 0.00, 0.16)),
       (0.030,   (1.00, 0.00, 0.00)),
       (0.215,   (1.00, 1.00, 0.00)),
       (0.400,   (0.00, 1.00, 0.00)),
       (0.586,   (0.00, 1.00, 1.00)),
       (0.770,   (0.00, 0.00, 1.00)),
       (0.954,   (1.00, 0.00, 1.00)),
       (1.000,   (1.00, 0.00, 0.75))
    )
```

To make a ColorScheme from this type of list, use:

```
make_colorscheme(gist_rainbow)
```
"""
function get_indexed_list_color(indexedlist, n)
   m = clamp(n, 0.0, 1.0)
   upper = max(2, findfirst(f -> m <= first(first(f)), indexedlist))
   lower = max(1, upper - 1)
   lowercolorvalues = last(indexedlist[lower])
   uppercolorvalues = last(indexedlist[upper])
   lowerv = first(indexedlist[lower])
   upperv = first(indexedlist[upper])
   lr, lg, lb = lowercolorvalues
   ur, ug, ub = uppercolorvalues
   r = lerp(m, lowerv, upperv, lr, ur)
   g = lerp(m, lowerv, upperv, lg, ug)
   b = lerp(m, lowerv, upperv, lb, ub)
   return round.((r, g, b), digits=6)
end

"""
    make_colorscheme(dict;
        length=100,
        category="",
        notes="")

Make a new ColorScheme from a dictionary of linear-segment information. Calls
`get_linear_segment_color(dict, n)` with `n` for every `length` value between 0 and 1.
"""
function make_colorscheme(dict::Dict;
        length=100,
        category="",
        notes="")
    cs = ColorScheme([RGB(get_linear_segment_color(dict, i)...)
        for i in range(0, stop=1, length=length)],
        category, notes)
    return cs
end

"""
    make_colorscheme(indexedlist;
        length=100,
        category="",
        notes="")

Make a ColorScheme using an 'indexed list' like this:

```
gist_rainbow = (
       (0.000, (1.00, 0.00, 0.16)),
       (0.030, (1.00, 0.00, 0.00)),
       (0.215, (1.00, 1.00, 0.00)),
       (0.400, (0.00, 1.00, 0.00)),
       (0.586, (0.00, 1.00, 1.00)),
       (0.770, (0.00, 0.00, 1.00)),
       (0.954, (1.00, 0.00, 1.00)),
       (1.000, (1.00, 0.00, 0.75))
)

make_colorscheme(gist_rainbow)
```

The first element of each item is the point on the color scheme.
"""
function make_colorscheme(indexedlist::Tuple;
        length=100,
        category="",
        notes="indexed list")
    cs = ColorScheme([RGB(get_indexed_list_color(indexedlist, i)...)
        for i in range(0, stop=1, length=length)],
        category, notes)
    return cs
end

"""
    make_colorscheme_new(f1::Function, f2::Function, f3::Function;
        model    = :RGB,
        length   = 100,
        category = "",
        notes    = "functional ColorScheme")

Make a ColorScheme using functions. Each function should take a value between 0
and 1 and return for that color component at each point on the ColorScheme,
depending on the color model.

The default color model is `:RGB`, and the functions should return values in the
appropriate range:

- f1 - [0.0 - 1.0]   - red
- f2 - [0.0 - 1.0]   - green
- f3 - [0.0 - 1.0]   - blue

For the `:HSV` color model:

- f1 - [0.0 - 360.0] - hue
- f2 - [0.0 - 1.0]   - saturataion
- f3 - [0.0 - 1.0]   - value (brightness)

For the `:LCHab` color model:

- f1 - [0.0 - 100.0] - luminance
- f2 - [0.0 - 100.0] - chroma
- f3 - [0.0 - 360.0] - hue
"""
function make_colorscheme(f1::Function, f2::Function, f3::Function;
        model    = :RGB,
        length   = 100,
        category = "",
        notes    = "functional ColorScheme")
    # output is always RGB for the moment
    cs = RGB[]
    if model == :LCHab
        clamp1 = (0.0, 100.0) #
        clamp2 = (0.0, 100.0) #
        clamp3 = (0.0, 360.0) # Hue is 0-360
    elseif model == :HSV
        clamp1 = (0.0, 360.0) # Hue is 0-360
        clamp2 = (0.0, 1.0) #
        clamp3 = (0.0, 1.0) #
    elseif model == :RGB
        clamp1 = (0.0, 1.0) #
        clamp2 = (0.0, 1.0) #
        clamp3 = (0.0, 1.0) #
    end
    counter = 0
    for i in range(0.0, stop=1.0, length=length)
        raw1, raw2, raw3 = f1(i), f2(i), f3(i)
        final1           = clamp(raw1, clamp1...)
        final2           = clamp(raw2, clamp2...)
        final3           = clamp(raw3, clamp3...)
        if model == :LCHab
            push!(cs, convert(RGB, LCHab(final1, final2, final3)))
        elseif model == :RGB
            push!(cs, RGB(final1, final2, final3))
        elseif model == :HSV
            push!(cs, convert(RGB, HSV(final1, final2, final3)))
        end
    end
    return ColorScheme(cs, category, notes)
end

end


text:
include("../server/simulation/stimulus/stimulus.jl")

function test_generating()
    stream = PoissonStream(UInt64(13163))
    println(next(stream))
    println(next(stream))
    println(next(stream))
    println(next(stream))
end

const TOTAL_SPIKES = 25

function test_stream_reset()
    println("Testing PoissonStream with reset")

    # Create a stream that generates spikes
    stream = PoissonStream(UInt64(13163), 0.2)
    # Capture isi value for testing after a reset
    isi = stream.isi
    # println("isi: ", stream.isi)

    # Create array to capture spike stream for testing after reset
    spikes = []

    # Fill array with output of stream
    for i in 1:TOTAL_SPIKES
        step!(stream)
        push!(spikes, output(stream))
        # println(output(stream))
    end

    # Now reset stream. This should reset the stream back to its begining.
    reset!(stream)

    # They should be the same
    @assert isi == stream.isi "ISI should be equal"

    # println("isi: ", stream.isi)
    # And the spikes should be too.
    for i in 1:TOTAL_SPIKES
        step!(stream)
        @assert output(stream) == spikes[i] "Invalid spike at index: $i"
        # println(output(stream))
    end

    println("Passed")
end

function test_stream()
    println("Testing PoissonStream with no reset")
    stream = PoissonStream(UInt64(13163), 0.2)
    isi = stream.isi
    # println("isi: ", stream.isi)
    spikes = []
    for i in 1:TOTAL_SPIKES
        step!(stream)
        push!(spikes, output(stream))
        # println(output(stream))
    end

    nisi = stream.isi
    @assert isi != nisi "ISI should NOT be equal"
    # println("isi: ", stream.isi)
    atleast_one_ne = false
    for i in 1:TOTAL_SPIKES
        step!(stream)

        if output(stream) != spikes[i]
            atleast_one_ne = true
            break;
        end
    end
    @assert atleast_one_ne "Atleast one spike should be different"
    println("Passed")
end

function test_stream_two()
    println("Testing two different poisson streams")

    stream = PoissonStream(UInt64(13163), 0.2)
    isi = stream.isi

    spikes = []
    for i in 1:TOTAL_SPIKES
        step!(stream)
        # println(output(stream))
        push!(spikes, output(stream))
    end

    stream2 = PoissonStream(UInt64(11862), 0.2)
    atleast_one_ne = false
    # println("----------")

    for i in 1:TOTAL_SPIKES
        step!(stream2)
        # println(output(stream))
        if output(stream2) != spikes[i]
            atleast_one_ne = true
            break;
        end
    end
    @assert atleast_one_ne "Atleast one spike should be different between streams"

    println("Passed")
end

test_stream_two()
test_stream()
test_stream_reset()


text:
module SchurPolynomials

export schur, semi_standard_young_tableaux

using AbstractAlgebra, Transducers, Dictionaries
using YAPP: Monomial, Polynomial
using BangBang
using AbstractAlgebra: Partition
using Transducers: next, complete, Foldable, Map

const RowsType = Vector{SubArray{Int, 1, Vector{Int}, Tuple{UnitRange{Int}}, true}}

"""
Iterates over all semi-standard Young tableaux of shape `λ`.
"""
struct SemiStandardYoungTableaux <: Foldable
    n::Int
    T::Generic.YoungTableau{Int}
    rows::RowsType
    function SemiStandardYoungTableaux(n, λ)
        fill = Vector{Int}(undef, sum(λ))
        rows = map(λ, Iterators.accumulate(+, λ)) do len, stop
            view(fill, stop-len+1:stop)
        end
        return new(n, YoungTableau(λ, fill), rows)
    end
end

"""
Iterates over all configurations of row `i` that are allowed by the previous row.
Iterator returns `nothing` but mutates `row`.
"""
struct RowConfigs <: Foldable
    rows::RowsType
    i::Int
    n::Int
end
Base.eltype(::RowConfigs) = Nothing
function Transducers.__foldl__(rf::RF, acc, (; rows, i, n)::RowConfigs) where {RF}
    row = @inbounds rows[i]
    acc = _row!(i, 1, row, rows, n, rf, acc)
    return complete(rf, acc)
end
# row i, column j
function _row!(i, j, row, rows, n, rf::RF, acc::A) where {RF, A}
    j > length(row) && return next(rf, acc, nothing)
    above = i == 1 ? 0 : rows[i-1][j]
    left = get(row, j-1, 1)
    @inbounds for row[j] in max(above + 1, left):n
        acc = _row!(i, j + 1, row, rows, n, rf, acc)
        acc isa Reduced && return acc
    end
    acc
end

function _worker!(i, rf::RF, acc, ssyt) where {RF}
    (; n, T, rows) = ssyt
    i > length(rows) && return next(rf, acc, nothing)
    foldl(RowConfigs(rows, i, n); init=acc) do acc, _
        _worker!(i + 1, rf, acc, ssyt)
    end
end

function Transducers.__foldl__(rf::RF, acc, ssyt::SemiStandardYoungTableaux) where {RF}
    acc = _worker!(1, rf, acc, ssyt)
    return complete(rf, acc)
end

function Base.length((; n, T)::SemiStandardYoungTableaux)
    return Int(prod(
        (n-i+j) // hooklength(T, i, j)
        for (i, j_max) in enumerate(T.part)
        for j in 1:j_max
    ))
end
Base.eltype(::SemiStandardYoungTableaux) = Generic.YoungTableau

semi_standard_young_tableaux(n, λ) = SemiStandardYoungTableaux(n, λ) |> Map() do T
    YoungTableau(T.part, copy(T.fill))
end

function schur(n, λ)
    SSYT = SemiStandardYoungTableaux(n, λ)
    (; T) = SSYT
    foldl(SSYT; init=Polynomial{Int}()) do p, _
        exponents = map(1:n) do i
            count(==(i), T.fill)
        end
        add!!(p, Monomial(exponents))
    end
end

end


text:
# Harris corner features from image files

using SensorFeatureTracking
using TestImages

# using Images, ImageView, ImageDraw, ImageFeatures, Gtk.ShortNames, VideoIO
#
#
# @show projdir = joinpath(dirname(@__FILE__), "..")
# srcdir = joinpath(projdir,"src")
# include(joinpath(srcdir, "SensorFeatureTracking.jl"))
# include(joinpath(srcdir, "Common.jl"))

# trying out different stuff
# img = testimage("walkbridge")
# img = testimage("livingroom")

# io = VideoIO.open("sometest.mp4")
# f = VideoIO.openvideo(io)
# f = VideoIO.testvideo("annie_oakley")
f = VideoIO.opencamera()
img = read(f);
close(f)
# size(img)

imgg = Gray.(img);
imshow(img)

# Just-in-time compiling, see usage below
@time Features(fastcorners(imgg, 12, 0.2)); # compile first
@time getapproxbestharris(imgg, 500, k=0.04); # compile first

# Fast corner features
@time begin
  corners = fastcorners(imgg, 12, 0.15)
  feats1 = Features(corners)
end

# selected best harris features
@time feats2 = getapproxbestharris(imgg, 500)

# separate drawing image
imcf = deepcopy(imgg);
imch = deepcopy(imgg);
for ft in feats1
  drawfeaturecircle2d!(imcf, ft, radius=3)
end
for ft in feats2
  drawfeaturecircle2d!(imch, ft, radius=3)
end

imshow(imcf)
imshow(imch)

# cs = getharriscorners(imgg,fracmax=0.6, k=0.04)
## corners = BitArray{2}(cs)
# fts1 = Features(cs)

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking