text:
"""
    emb_sz_rule(n_cat)

Compute an embedding size corresponding to the number of classes for a
categorical variable using the rule of thumb present in python fastai.
(see https://github.com/fastai/fastai/blob/2742fe844573d06e700f869839fb9ec5f3a9bca9/fastai/tabular/model.py#L12)
"""
emb_sz_rule(n_cat) = min(600, round(Int, 1.6 * n_cat^0.56))

"""
    get_emb_sz(cardinalities::AbstractVector, [size_overrides::AbstractVector])

Given a vector of `cardinalities` of each categorical column
(i.e. each element of `cardinalities` is the number of classes in that categorical column),
compute the output embedding size according to [`emb_sz_rule`](#).
Return a vector of tuples where each element is `(in_size, out_size)` for an embedding layer.

## Keyword arguments

- `size_overrides`: A collection of integers (or `nothing` to skip override) where the value present at any index
    will be used to as the output embedding size for that column.
"""
get_emb_sz(cardinalities::AbstractVector{<:Integer}, size_overrides=fill(nothing, length(cardinalities))) =
    map(zip(cardinalities, size_overrides)) do (cardinality, override)
        emb_dim = isnothing(override) ? emb_sz_rule(cardinality + 1) : Int64(override)
        return (cardinality + 1, emb_dim)
    end

sigmoidrange(x, low, high) = @. Flux.sigmoid(x) * (high - low) + low

function tabular_embedding_backbone(embedding_sizes, dropout_rate=0.)
    embedslist = [Flux.Embedding(ni, nf) for (ni, nf) in embedding_sizes]
    emb_drop = iszero(dropout_rate) ? identity : Dropout(dropout_rate)
    Chain(
        x -> tuple(eachrow(x)...),
        Parallel(vcat, embedslist),
        emb_drop
    )
end

tabular_continuous_backbone(n_cont) = BatchNorm(n_cont)

"""
    TabularModel(catbackbone, contbackbone, [finalclassifier]; kwargs...)

Create a tabular model which operates on a tuple of categorical values
(label or one-hot encoded) and continuous values.
The categorical backbones (`catbackbone`) and continuous backbone (`contbackbone`) operate on each element of the input tuple.
The output from these backbones is then passed through a series of linear-batch norm-dropout layers before a `finalclassifier` block.

## Keyword arguments

- `outsize`: The output size of the final classifier block. For single classification tasks,
    this would be the number of classes, and for regression tasks, this would be the
    number of target continuous variables.
- `layersizes`: A vector of sizes for each hidden layer in the sequence of linear layers.
- `dropout_rates`: Dropout probabilities for the linear-batch norm-dropout layers.
    This could either be a single number which would be used for for all the layers,
    or a collection of numbers which are cycled through for each layer.
- `batchnorm`: Set to `false` to skip each batch norm in the linear-batch norm-dropout sequence.
- `activation`: The activation function to use in the classifier layers.
- `linear_first`: Controls if the linear layer comes before or after batch norm and dropout.
"""
function TabularModel(
        catbackbone,
        contbackbone;
        outsize,
        layersizes=(200, 100),
        kwargs...)
    TabularModel(catbackbone, contbackbone, Dense(layersizes[end], outsize); layersizes=layersizes, kwargs...)
end

function TabularModel(
        catbackbone,
        contbackbone,
        finalclassifier;
        layersizes=(200, 100),
        dropout_rates=0.,
        batchnorm=true,
        activation=Flux.relu,
        linear_first=true)

    tabularbackbone = Parallel(vcat, catbackbone, contbackbone)

    classifierin = mapreduce(layer -> size(layer.weight)[1], +, catbackbone[2].layers;
                             init = contbackbone.chs)
    dropout_rates = Iterators.cycle(dropout_rates)
    classifiers = []

    first_ps, dropout_rates = Iterators.peel(dropout_rates)
    push!(classifiers, linbndrop(classifierin, first(layersizes);
                                 use_bn=batchnorm, p=first_ps, lin_first=linear_first, act=activation))

    for (isize, osize, p) in zip(layersizes[1:(end-1)], layersizes[2:end], dropout_rates)
        layer = linbndrop(isize, osize; use_bn=batchnorm, p=p, act=activation, lin_first=linear_first)
        push!(classifiers, layer)
    end

    Chain(
        tabularbackbone,
        classifiers...,
        finalclassifier
    )
end

"""
    TabularModel(n_cont, outsize, [layersizes; kwargs...])

Create a tabular model which operates on a tuple of categorical values
(label or one-hot encoded) and continuous values. The default categorical backbone (`catbackbone`) is
a [`Flux.Parallel`](https://fluxml.ai/Flux.jl/stable/models/layers/#Flux.Parallel) set of `Flux.Embedding` layers corresponding to each categorical variable.
The default continuous backbone (`contbackbone`) is a single [`Flux.BatchNorm`](https://fluxml.ai/Flux.jl/stable/models/layers/#Flux.BatchNorm).
The output from these backbones is concatenated then passed through a series of linear-batch norm-dropout layers before a `finalclassifier` block.

## Arguments

- `n_cont`: The number of continuous columns.
- `outsize`: The output size of the model.
- `layersizes`: A vector of sizes for each hidden layer in the sequence of linear layers.

## Keyword arguments

- `cardinalities`: A collection of sizes (number of classes) for each categorical column.
- `size_overrides`: An optional argument which corresponds to a collection containing
    embedding sizes to override the value returned by the "rule of thumb" for a particular index
    corresponding to `cardinalities`, or `nothing`.
"""
function TabularModel(
        n_cont::Number,
        outsize::Number,
        layersizes=(200, 100);
        cardinalities,
        size_overrides=fill(nothing, length(cardinalities)))
    embedszs = get_emb_sz(cardinalities, size_overrides)
    catback = tabular_embedding_backbone(embedszs)
    contback = tabular_continuous_backbone(n_cont)

    TabularModel(catback, contback; layersizes=layersizes, outsize=outsize)
end


text:
include("../common.jl")
build_julia(ARGS, v"1.3.1")

ext: jl

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# module containing patterns used for Metamaterial waveguide structure
module MetamaterialWaveguide
    using Devices, Clipper
    import Devices: Meta, GDSMeta
    export cross_polygon, mwg_unit_cell!, mwg_last_cell!, mwg_input_cap!, mwg_cell_with_claw!

    """
        cross_polygon(hline_width, hline_height, vline_width,
        vline_height; center=true)
    Create a cross-shaped polygon of specified widths and heights of horizontal
    line and vertical line. The keyword argument `center=true` determines whether
    the returned polygon is centered at the origin or at distance hline_width/2
    along the x-axis away from the origin.
    """
    function cross_polygon(hline_width, hline_height, vline_width,
        vline_height; center=true)

        if hline_width < vline_width
            error("The width of the horizontal line must be larger than that of
                the vertical line")
        elseif hline_height > vline_height
            error("The height of the vertical line must be larger than that of
                the horizontal line")
        end

        offset =  (~center) * Point(hline_width/2, typeof(hline_width)(0))
        hline = centered(Rectangle(hline_width, hline_height)) + offset
        vline = centered(Rectangle(vline_width, vline_height)) + offset

        poly = clip(Clipper.ClipTypeUnion, hline, vline)[1]
        return poly
    end

    """
        mwg_unit_cell!{T,S}(c::Cell{T,S}, Lr_width, Lr_height, Cr_overlap,
        res_ground_length, res_conductor_width, res_conductor_gap,
        Ck_terminal_width, Ck_overlap, wg_meander_length, wg_meander_nlegs::Int,
        wg_ground_gap, wg_conductor_width, wg_conductor_gap, centerline_width,
        centerline_length, meta=GDSMeta(0,0))
    Create a unit cell of the metamaterial waveguide.
    - `c` : a cell to render metamaterial waveguide in
    - `Lr_width` : width of the rectangular coil Lr
    - `Lr_height` : height of the rectangular coil Lr
    - `Cr_overlap` : overlap length of the interdigital capacitor Cr forming
        the resonator
    - `res_ground_length` : length of the line from the top of the resonator
        to the ground. This also sets the gap between the resonator
        and the ground.
    - `res_conductor_width` : width of the conductor forming a resonator
    - `res_conductor_gap` : gap between conductor traces forming a resonator
    - `Ck_terminal_width` : terminal width of the coupling capacitor Ck between
        the resonator and the centerline
    - `Ck_overlap` : overlap length of the coupling capacitor Ck
    - `wg_meander_length` : the length of meandering of waveguide. The trace of
        the meander goes back and forth vertically with this distance.
    - `wg_meander_nlegs` : number of legs in the meander forming the waveguide
    - `wg_ground_gap` : the distance between the meander trace of the waveguide
        and the ground
    - `wg_conductor_width` : width of the conductor forming the meander
    - `wg_conductor_gap` : gap between conductor traces forming the meander
    - `centerline_width` : width of the centerline that connects the resonator
        with waveguide.
    - `centerline_length` : length of the centerline. This has to be equal to or
        greater than `Lr_width`
    """
    function mwg_unit_cell!(c::Cell{T,S}, Lr_width, Lr_height, Cr_overlap,
        res_ground_length, res_conductor_width, res_conductor_gap, Ck_terminal_width,
        Ck_overlap, wg_meander_length, wg_meander_nlegs::Int, wg_ground_gap,
        wg_conductor_width, wg_conductor_gap, centerline_width, centerline_length,
        meta=GDSMeta(0,0)) where {T,S}
        if Lr_width > centerline_length
            error("The length of the centerline must be larger than the width of the coil")
        end
        wg_length = (wg_meander_nlegs - 1) * (wg_conductor_width + wg_conductor_gap) + wg_conductor_width / 2
        res_height = 2 * (Lr_height + res_conductor_gap) + Cr_overlap

        # define a mask for unit cell
        wg_mask_width = 2 * wg_length + centerline_length
        wg_mask_height = wg_meander_length + 2 * (wg_conductor_width + wg_ground_gap)

        res_mask_width = Lr_width + 2 * res_ground_length
        res_mask_height = centerline_width + 2 * (res_ground_length + res_height + Ck_terminal_width + Ck_overlap + 3 * res_conductor_gap)

        # global mask = res_mask ∪ wg_mask
        # mask = Cell{T,S}(uniquename("mask"))
        # mask_poly = cross_polygon(wg_mask_width, wg_mask_height, res_mask_width,
            # res_mask_height; center=false)
        # render!(mask, mask_poly)

        # lumped-element resonator with coupling capacitor
        res_with_Ck = Cell{T,S}(uniquename("lumped_elem_resonator_with_Ck"))
        lumped_elem_resonator_with_Ck!(res_with_Ck, Lr_width, Lr_height, Cr_overlap, Ck_overlap,
            Ck_terminal_width, res_conductor_width, res_conductor_gap, meta,
            ground_length=res_ground_length, Ck_terminal_width2=zero(T))

        res = Cell{T,S}(uniquename("resonator"))
        push!(res.refs, CellReference(res_with_Ck, Point(wg_mask_width/2, centerline_width/2)))
        push!(res.refs, CellReference(res_with_Ck,
            Point(wg_mask_width/2, -centerline_width/2), xrefl=true))

        # meander forming waveguide part of MWG
        meander = Cell{T,S}(uniquename("meander"))
        rectangular_meander!(meander, wg_meander_length, wg_meander_nlegs,
            wg_conductor_width, wg_conductor_gap, meta;
            first_length=wg_meander_length + wg_conductor_width,
            last_length=(wg_meander_length + centerline_width)/2,
            start_upward=true)

        wg = Cell{T,S}(uniquename("waveguide"))
        push!(wg.refs, CellReference(meander, Point(-wg_conductor_width/2,
            -wg_conductor_width-wg_meander_length/2)))
        push!(wg.refs, CellReference(meander,
            Point(wg_conductor_width/2+ wg_mask_width,
            wg_conductor_width + wg_meander_length/2), rot=π))

        # centerline
        centerline = Cell{T,S}(uniquename("centerline"))
        render!(centerline, centered(Rectangle(centerline_length, centerline_width)) +
            Point(wg_mask_width/2, zero(T)), meta)

        # # take the negative trace of each element
        # res_neg = negative(res, mask)
        # wg_neg = negative(wg, mask)
        # centerline_neg = negative(centerline, mask)

        # take the intersection res_neg ∩ wg_neg ∩ centerline_neg
        # unit_cell_poly = clip(Clipper.ClipTypeIntersection,
        #     polygon.(elements(res_neg)), polygon.(elements(wg_neg)))
        # unit_cell_poly = clip(Clipper.ClipTypeIntersection,
        #     unit_cell_poly, polygon.(elements(centerline_neg)))
        #
        # for p in unit_cell_poly
        #     render!(c, p, meta)
        # end
        for part in [res, wg, centerline]
            push!(c.refs, CellReference(part, Point(zero(T), zero(T))))
        end
        return c
    end

    """
    Metamaterial waveguide cell with claw coming out from the resonator node
    """
    function mwg_cell_with_claw!(c::Cell{T,S}, Lr_width, Lr_height, Cr_overlap,
        res_ground_length, res_conductor_width, res_conductor_gap, Ck_terminal_width,
        Ck_overlap, wg_meander_length, wg_meander_nlegs::Int, wg_ground_gap,
        wg_conductor_width, wg_conductor_gap, centerline_width, centerline_length,
        res_claw_gap, res_claw_conductor_width, claw_width, claw_length,
        claw_ground_gap, claw_dist, meta::Meta=GDSMeta(0,0); rot=0) where {T,S}

        wg_length = ((wg_meander_nlegs - 1) * (wg_conductor_width + wg_conductor_gap)
            + wg_conductor_width/2)
        res_height = 2 * (Lr_height + res_conductor_gap) + Cr_overlap

        # define a mask for unit cell
        wg_mask_width = 2 * wg_length + centerline_length
        wg_mask_height = wg_meander_length + 2 * (wg_conductor_width + wg_ground_gap)

        res_mask_width = Lr_width + 2 * res_ground_length
        res_mask_height = centerline_width + 2 * (res_ground_length + res_height
            + Ck_terminal_width + Ck_overlap + 3 * res_conductor_gap)

        claw_mask_left = wg_mask_width/2 + Lr_width/2
        claw_mask_right = (claw_mask_left + res_claw_gap + res_claw_conductor_width
            + claw_length + claw_ground_gap)
        claw_mask_vdist = claw_dist/2 - claw_ground_gap
        claw_mask_hdist = res_claw_gap + res_claw_conductor_width + claw_ground_gap

        wg_mask = Rectangle(Point(zero(T), -wg_mask_height/2),
            Point(wg_mask_width/2, wg_mask_height/2))
        res_mask = centered(Rectangle(res_mask_width, res_mask_height)) + Point(wg_mask_width/2, zero(T))

        claw_mask_bulk = Rectangle(Point(claw_mask_left, -res_mask_height/2),
            Point(claw_mask_right, res_mask_height/2))
        claw_mask_rm = Rectangle(Point(claw_mask_left + claw_mask_hdist,
            -claw_mask_vdist), Point(claw_mask_right, claw_mask_vdist))
        claw_mask = clip(Clipper.ClipTypeDifference, claw_mask_bulk, claw_mask_rm)[1]

        # global mask = res_mask ∪ wg_mask ∪ claw_mask
        # mask = Cell{T,S}(uniquename("mask"))
        # mask_poly = clip(Clipper.ClipTypeUnion, wg_mask, res_mask)[1]
        # mask_poly = clip(Clipper.ClipTypeUnion, mask_poly, claw_mask)[1]
        # render!(mask, mask_poly, meta_mask)

        # lumped-element resonator with coupling capacitor
        res_with_Ck = Cell{T,S}(uniquename("lumped_elem_resonator_with_Ck"))
        lumped_elem_resonator_with_Ck!(res_with_Ck, Lr_width, Lr_height, Cr_overlap, Ck_overlap,
        Ck_terminal_width, res_conductor_width, res_conductor_gap, meta,
        ground_length=res_ground_length, Ck_terminal_width2=zero(T))

        res = Cell{T,S}(uniquename("resonator"))
        push!(res.refs, CellReference(res_with_Ck,
            Point(wg_mask_width/2, centerline_width/2)))
        push!(res.refs, CellReference(res_with_Ck,
            Point(wg_mask_width/2, -centerline_width/2), xrefl=true))

            # meander forming waveguide part of MWG
        meander = Cell{T,S}(uniquename("meander"))
        rectangular_meander!(meander, wg_meander_length, wg_meander_nlegs,
            wg_conductor_width, wg_conductor_gap, meta;
            first_length=wg_meander_length + wg_conductor_width,
            last_length=(wg_meander_length + centerline_width)/2,
            start_upward=true)

        wg = Cell{T,S}(uniquename("waveguide"))
        push!(wg.refs, CellReference(meander, Point(-wg_conductor_width/2,
            -wg_conductor_width-wg_meander_length/2)))

        # centerline
        centerline = Cell{T,S}(uniquename("centerline"))
        render!(centerline, centered(Rectangle(centerline_length, centerline_width)) +
            Point(wg_mask_width/2, zero(T)), meta)

        # claw
        claw = Cell{T,S}(uniquename("claw"))
        res_claw_conn_vdist = (res_mask_height/2 - Lr_height - res_ground_length
            - 2 * res_conductor_gap - Cr_overlap)
        render!(claw, Rectangle(Point(claw_mask_left, res_claw_conn_vdist - 2*res_conductor_width),
            Point(claw_mask_left + res_claw_gap, res_claw_conn_vdist)), meta)
        render!(claw, Rectangle(Point(claw_mask_left, -res_claw_conn_vdist + 2*res_conductor_width),
            Point(claw_mask_left + res_claw_gap, -res_claw_conn_vdist)), meta)

            claw_left = claw_mask_left + res_claw_gap + res_claw_conductor_width
            claw_right = claw_left + claw_length

        render!(claw, Rectangle(Point(claw_left - res_claw_conductor_width,
            -res_claw_conn_vdist), Point(claw_left, res_claw_conn_vdist)), meta)

        render!(claw, Rectangle(Point(claw_left, claw_mask_vdist + claw_ground_gap),
            Point(claw_right, claw_mask_vdist + claw_width + claw_ground_gap)), meta)
        render!(claw, Rectangle(Point(claw_left, -claw_mask_vdist - claw_ground_gap),
            Point(claw_right, -claw_mask_vdist - claw_width - claw_ground_gap)), meta)

        # This part discarded due to instability of Clipper
        # take the negative trace of each element
        # res_neg = negative(res, mask)
        # wg_neg = negative(wg, mask)
        # centerline_neg = negative(centerline, mask)
        # claw_neg = negative(claw, mask)
        # take the intersection res_neg ∩ wg_neg ∩ centerline_neg
        # cell_poly = clip(Clipper.ClipTypeIntersection,
        #     polygon.(elements(res_neg)), polygon.(elements(wg_neg)))
        # cell_poly = clip(Clipper.ClipTypeIntersection,
        #     cell_poly, polygon.(elements(centerline_neg)))
        # cell_poly = clip(Clipper.ClipTypeIntersection,
        #     cell_poly, polygon.(elements(claw_neg)))
        # for p in cell_poly
            # render!(c, p, meta)
        # end

        for part in [res, wg, centerline, claw]# ,mask]
            push!(c.refs, CellReference(part, Point(zero(T), zero(T)), rot=rot))
        end

        # render!(c, mask, meta_mask)
        return c
    end

    #
    # """
    # First cell for the metamaterial waveguide including the input capacitor
    # """
    # function mwg_first_cell!(c::Cell{T,S}, Lr_width, Lr_height, Cr_overlap_first,
    #     res_ground_length, res_conductor_width, res_conductor_gap, Ck_terminal_width,
    #     Ck_overlap, wg_meander_length, wg_meander_nlegs::Int, wg_ground_gap,
    #     wg_conductor_width, wg_conductor_gap, centerline_width, centerline_length,
    #     res_claw_gap, res_claw_conductor_width_first, claw_width_first,
    #     claw_length_first, claw_ground_gap_first, input_cap_width, input_cap_gap, ground_width_first,
    #     meta::Meta=GDSMeta(0,0)) where {T,S}
    #     tmp = Cell{T,S}(uniquename("last_cell"))
    #     mwg_last_cell!(tmp, Lr_width, Lr_height, Cr_overlap_first,
    #         res_ground_length, res_conductor_width, res_conductor_gap,
    #         Ck_terminal_width, Ck_overlap, wg_meander_length, wg_meander_nlegs,
    #         wg_ground_gap, wg_conductor_width, wg_conductor_gap,
    #         centerline_width, centerline_length, res_claw_gap,
    #         res_claw_conductor_width_first, claw_width_first, claw_length_first,
    #         claw_ground_gap_first, input_cap_width, input_cap_gap,
    #         ground_width_first, meta)
    #     for p in polygon.(tmp.elements)
    #         render!(c, Rotation(π)(p)+Point(width(bounds(tmp)), zero(T)), meta)
    #     end
    #     return c
    # end

    # """
    # Input capacitor for the metamaterial waveguide
    # """
    # function mwg_input_cap!(c::Cell{T,S}, res_mask_width, res_mask_height,
    #     wg_mask_width, wg_mask_height, wg_meander_length, wg_meander_nlegs::Int,
    #     wg_ground_gap, wg_conductor_width, wg_conductor_gap, centerline_width,
    #     input_cap_terminal_width, input_cap_terminal_width2, input_cap_overlap,
    #     input_cap_finger_width, input_cap_finger_spacing, input_cap_end_gap,
    #     cpw_gap, cpw_trace, meta=GDSMeta(0,0)) where {T,S}
    #
    #     wg_length = (wg_meander_nlegs - 1) * (wg_conductor_width + wg_conductor_gap) + wg_conductor_width/2
    #
    #     input_cap_terminal_length = 2 * wg_conductor_width + wg_meander_length
    #
    #     # define a mask for input capacitor
    #     wg_mask = (centered(Rectangle(wg_mask_width/2, wg_mask_height))
    #         + Point((wg_mask_width/2+res_mask_width)/2, zero(T)))
    #     res_mask = (centered(Rectangle(res_mask_width, res_mask_height))
    #         + Point(res_mask_width/2, zero(T)))
    #
    #     # global mask = res_mask ∪ wg_mask
    #     mask = clip(Clipper.ClipTypeUnion, wg_mask, res_mask)[1]
    #
    #     # meander forming waveguide part of MWG
    #     meander = Cell{T,S}(uniquename("meander"))
    #     rectangular_meander!(meander, wg_meander_length, wg_meander_nlegs,
    #         wg_conductor_width, wg_conductor_gap;
    #         first_length=wg_meander_length + wg_conductor_width,
    #         last_length=(wg_meander_length + centerline_width)/2,
    #         start_upward=true)
    #     wg = Cell{T,S}(uniquename("waveguide"))
    #     push!(wg.refs, CellReference(meander,
    #         Point(wg_conductor_width/2 + wg_mask_width/2 + res_mask_width/2,
    #         wg_conductor_width + wg_meander_length/2), rot=π))
    #
    #     # input capacitor
    #     cap = Cell{T,S}("capacitor")
    #     interdigit = Cell{T,S}("interdigit")
    #     interdigit_capacitor!(interdigit, input_cap_terminal_width, input_cap_terminal_length,
    #         input_cap_overlap, input_cap_finger_width, input_cap_finger_spacing, input_cap_end_gap,
    #         terminal_width2=input_cap_terminal_width2)
    #     push!(cap.refs, CellReference(interdigit, Point(zero(T), -input_cap_terminal_length/2)))
    #
    #     # centerline
    #     conn_line = Cell{T,S}(uniquename("connecting_line"))
    #     render!(conn_line, Rectangle(Point(width(bounds(interdigit)), -centerline_width/2),
    #         Point(bounds(wg).ll.x, centerline_width/2)), meta)
    #     push!(cap.refs, CellReference(conn_line, Point(zero(T), zero(T))))
    #
    #     input_cap_neg = negative(cap, mask)
    #     input_wg_neg = negative(wg, mask)
    #
    #     # take the intersection res_neg ∩ wg_neg ∩ centerline_neg
    #     input_cap_poly = clip(Clipper.ClipTypeIntersection,
    #         polygon.(elements(input_cap_neg)), polygon.(elements(input_wg_neg)))
    #
    #     input_path = Path(Point(typeof(input_cap_terminal_length)(0),
    #         typeof(input_cap_terminal_length)(0)), α0 = 0)
    #     input_taper_len = (wg_mask_width - res_mask_width)/2
    #     input_trace = input_cap_terminal_length
    #     input_gap = input_trace * cpw_gap/cpw_trace
    #
    #     straight!(input_path, input_taper_len, Paths.CPW(
    #         t->(cpw_trace * (input_taper_len - t)/input_taper_len + t/input_taper_len * input_trace),
    #         t->(cpw_gap * (input_taper_len - t)/input_taper_len + t/input_taper_len * input_gap)))
    #
    #     render!(c, input_path, meta)
    #     trans = Translation(Point(input_taper_len, zero(T)))
    #     for p in input_cap_poly
    #         render!(c, trans(p), meta)
    #     end
    #     return c
    # end
end

# package code goes here


text:
"""Type constructor

frame - stores ArrayNArrays{T,N}
bls - Int, size of ordinary block
bln - Int, number of blocks
datasize - Int, size of data stored (in each direction the same)
sqr - Bool, is the last block size a same as ordinary's block size
"""
mutable struct SymmetricTensor{T <: AbstractFloat, N}
    frame::ArrayNArrays{T,N}
    bls::Int
    bln::Int
    dats::Int
    sqr::Bool
    function (::Type{SymmetricTensor})(frame::ArrayNArrays{T,N};
        testdatstruct::Bool = true) where {T, N}
        bls = size(frame[fill(1, N)...,], 1)
        bln = size(frame, 1)
        last_block = size(frame[end], 1)
        dats = bls * (bln - 1) + last_block
        if testdatstruct
          frtest(frame)
        end
        new{T, N}(frame, bls, bln, dats, bls == last_block)
    end
end

"""
  unfold(ar::Array{N}, mode::Int)

Returns an matrix being an unfold of N dims array in given mode.

```jldoctest
julia> A = reshape(collect(1.:8.), 2, 2, 2);

julia> unfold(A, 1)
2×4 Array{Float64,2}:
 1.0  3.0  5.0  7.0
 2.0  4.0  6.0  8.0

 julia> unfold(A, 2)
 2×4 Array{Float64,2}:
  1.0  2.0  5.0  6.0
  3.0  4.0  7.0  8.0

  julia> unfold(A, 3)
  2×4 Array{Float64,2}:
   1.0  2.0  3.0  4.0
   5.0  6.0  7.0  8.0
```
"""
function unfold(ar::Array{T,N}, mode::Int) where {T <: Real, N}
    C = [1:mode-1; mode+1:N]
    i = size(ar)
    k = prod(i[C])
    return reshape(permutedims(ar,[mode; C]), i[mode], k)
end

"""
  issymetric(ar::Array{N}, atol::Float64)

Returns: Assertion Error if not symmetric given a tolerance.

```jldoctest
julia> A = reshape(collect(1.:8.), 2, 2, 2);

julia> julia> issymetric(A)
ERROR: AssertionError: not symmetric
```
"""
function issymetric(ar::Array{T, N}, atol::Float64 = 1e-7) where {T <: AbstractFloat, N}
  for i=2:N
     maximum(abs.(unfold(ar, 1)-unfold(ar, i))) < atol ||throw(AssertionError("not symmetric"))
  end
end

"""
  frtest(data::ArrayNArrays{T,N})

Returns assertion error if: all sizes of nullable array not equal, at least
  some undergiagonal block not null, some blocks (not last) not squared,
  some diagonal blocks not symmetric.
"""
function frtest(data::ArrayNArrays{T,N}) where {T <: AbstractFloat, N}
  bln = size(data, 1)
  bls = size(data[fill(1, N)...,], 1)
  all(collect(size(data)) .== bln) || throw(AssertionError("frame not has non-equal dimensions"))

  for i in CartesianIndices(size(data))
      if data[i]!=nothing && !issorted(i.I)
          throw(AssertionError("underdiagonal block not null"))
      end
  end

  for i in pyramidindices(N, bln-1)
    @inbounds all(collect(size(data[i...])) .== bls)||
        throw(AssertionError("$i block not square"))
  end
  for i=1:bln
    @inbounds issymetric(data[fill(i, N)...,])
  end
end

"""
  pyramidindices(dims::Int, tensize::Int)

```jldoctest
julia> pyramidindices(2,3)
6-element Array{Tuple{Int64,Int64},1}:
 (1,1)
 (1,2)
 (1,3)
 (2,2)
 (2,3)
 (3,3)
```
"""
function pyramidindices(dims::Int, tensize::Int)
    multinds = Tuple{fill(Int,dims)...,}[]
    @eval begin
        @nloops $dims i x -> (x==$dims) ? (1:$tensize) : (i_{x+1}:$tensize) begin
            @inbounds multind = @ntuple $dims x -> i_{$dims-x+1}
            push!($multinds, multind)
        end
    end
    multinds
end

"""

    sizetest(dats::Int, bls::Int)

Returns: DimensionMismatch if blocks size is grater than data size.

```jldoctest
julia> SymmetricTensors.sizetest(2,3)
ERROR: DimensionMismatch("bad block size 3 > 2")
```
"""
sizetest(dats::Int, bls::Int) =
  (dats >= bls > 0)|| throw(DimensionMismatch("bad block size $bls > $dats"))

"""
  getblockunsafe(st::SymmetricTensor, i::Tuple)

Returns a block from Symmetric Tensor, unsafe works only if multi-index is sorted
"""
getblockunsafe(st::SymmetricTensor, mulind::Tuple) = st.frame[mulind...]

"""
    getblock(st::SymmetricTensor, i::Tuple)

Returns a block from Symmetric Tensor, works for all multi-indices also not sorted
"""
function getblock(st::SymmetricTensor, mulind::Tuple)
  ind = sortperm([mulind...])
  permutedims(getblockunsafe(st, mulind[ind]), invperm(ind))
end

"""
    getindex(st::SymmetricTensor, mulind::Tuple)

Returns a Symmetric Tensor element for a given multi-index
"""
function getindex(st::SymmetricTensor, mulind::Int...)
  b = st.bls
  j = map(k -> div((k-1), b)+1, mulind)
  i = map(k -> ((k-1)%b)+1, mulind)
  getblock(st, j)[i...]
end

"""
    setindexunsafe!(st::SymmetricTensor{T,N}, x::T,  mulind::Int...)

Unsafe change a SymmetricTensors value at the given multi-index
"""
function setindexunsafe!(st::SymmetricTensor{T,N}, x::T,  mulind::Int...) where {T <: AbstractFloat, N}
    b = st.bls
    j = map(k -> div((k-1), b)+1, mulind)
    i = map(k -> ((k-1)%b)+1, mulind)
    if j == (unique(j)...,)
        getblockunsafe(st, j)[i...] = x
    else
        for p in fixpointperms(j)
            getblockunsafe(st, j)[i[p]...] = x
        end
    end
end

"""
    setindex!(st::SymmetricTensor{T,N}, x::T,  mulind::Int...)

Change a SymmetricTensors value at the given multi-index
"""
function setindex!(st::SymmetricTensor{T,N}, x::T,  mulind::Int...) where {T <: AbstractFloat, N}
    setindexunsafe!(st, x, sort([mulind...])...)
end
"""

    ind2range(i::Int, bls::Int, dats::Int)

Returns a range given index i, size of a block and size of data

```jldoctest
julia> ind2range(2,3,5)
4:5
```
"""
ind2range(i::Int, bls::Int, dats::Int) = (i-1)*bls+1: ((i*bls <= dats) ? i*bls : dats)

"""
  SymmetricTensor(data::Array{N}, bls::Int)

Returns: data in SymmetricTensor form.
```jldoctest
julia> a = reshape(collect(1.:16.), 4, 4);

julia> SymmetricTensor(a*a')
SymmetricTensors.SymmetricTensor{Float64,2}(Union{Array{Float64,2}, Void}[[276.0 304.0; 304.0 336.0][332.0 360.0; 368.0 400.0]; nothing [404.0 440.0; 440.0 480.0]], 2, 2, 4, true)
```
"""
function SymmetricTensor(data::Array{T, N}, bls::Int = 2) where {T <: AbstractFloat, N}
  issymetric(data)
  dats = size(data,1)
  sizetest(dats, bls)
  bln = mod(dats,bls)==0 ?  dats÷bls : dats÷bls + 1
  symten = arraynarrays(T, fill(bln, N)...,)
  for writeind in pyramidindices(N, bln)
      readind = map(k::Int -> ind2range(k, bls, dats), writeind)
      @inbounds symten[writeind...] = data[readind...]
  end
  SymmetricTensor(symten)
end

"""
  Array(st::SymmetricTensor{N})

Return N dims array converted from SymmetricTensor type

"""
function Array(st::SymmetricTensor{T,N}) where {T<:AbstractFloat, N}
  array = zeros(T, fill(st.dats, N)...,)
  for i = 1:(st.bln^N)
    dims = (fill(st.bln, N)...,)
    readind= Tuple(CartesianIndices(dims)[i])
    writeind = map(k -> ind2range(readind[k], st.bls, st.dats), 1:N)
    @inbounds array[writeind...] = getblock(st, readind)
  end
  array
end
# ---- basic operations on Symmetric Tensors ----

"""
  diag(st::SymmetricTensor{N})

Return vector of floats, the super-diag of st

"""
diag(st::SymmetricTensor{T,N}) where {T<: AbstractFloat, N} = map(i->st[fill(i, N)...,], 1:st.dats)

# implements simple operations on bs structure

for f = (:+, :-)
  @eval function ($f)(st::SymmetricTensor{T,N}...) where {T <: AbstractFloat, N}
    for s in st[2:end]
      s.dats == st[1].dats || throw(DimensionMismatch("dimensions must match"))
    end
    narg = size(st, 1)
    stret = similar(st[1].frame)
    for i in pyramidindices(N, st[1].bln)
      @inbounds stret[i...] = broadcast($f, map(t -> getblockunsafe(t, i), st)...)
    end
    SymmetricTensor(stret; testdatstruct = false)
  end
end

for f = (:+, :-, :*, :/)
  @eval function ($f)(st::SymmetricTensor{T, N}, numb::S) where {T <: AbstractFloat, S <: Real, N}
      stret = similar(st.frame)
      for i in pyramidindices(N, st.bln)
        @inbounds stret[i...] = broadcast($f, getblockunsafe(st, i), numb)
      end
      SymmetricTensor(stret; testdatstruct = false)
  end
end

for f = (:+, :*)
  @eval ($f)(numb::S, st::SymmetricTensor{T}) where {T <: AbstractFloat, S <: Real} =
  ($f)(st::SymmetricTensor{T}, numb::S)
end


text:
export DiagBlockSparse, DiagBlockSparseTensor

# DiagBlockSparse can have either Vector storage, in which case
# it is a general DiagBlockSparse tensor, or scalar storage,
# in which case the diagonal has a uniform value
struct DiagBlockSparse{ElT,VecT,N} <: TensorStorage{ElT}
  data::VecT
  diagblockoffsets::BlockOffsets{N}  # Block number-offset pairs

  # Nonuniform case
  function DiagBlockSparse(
    data::VecT, blockoffsets::BlockOffsets{N}
  ) where {VecT<:AbstractVector{ElT},N} where {ElT}
    return new{ElT,VecT,N}(data, blockoffsets)
  end

  # Uniform case
  function DiagBlockSparse(data::VecT, blockoffsets::BlockOffsets{N}) where {VecT<:Number,N}
    return new{VecT,VecT,N}(data, blockoffsets)
  end
end

function DiagBlockSparse(
  ::Type{ElT}, boffs::BlockOffsets, diaglength::Integer
) where {ElT<:Number}
  return DiagBlockSparse(zeros(ElT, diaglength), boffs)
end

function DiagBlockSparse(boffs::BlockOffsets, diaglength::Integer)
  return DiagBlockSparse(Float64, boffs, diaglength)
end

function DiagBlockSparse(
  ::Type{ElT}, ::UndefInitializer, boffs::BlockOffsets, diaglength::Integer
) where {ElT<:Number}
  return DiagBlockSparse(Vector{ElT}(undef, diaglength), boffs)
end

function DiagBlockSparse(::UndefInitializer, boffs::BlockOffsets, diaglength::Integer)
  return DiagBlockSparse(Float64, undef, boffs, diaglength)
end

diagblockoffsets(D::DiagBlockSparse) = D.diagblockoffsets

blockoffsets(D::DiagBlockSparse) = D.diagblockoffsets

function findblock(
  D::DiagBlockSparse{<:Number,<:Union{Number,AbstractVector},N}, block::Block{N}; vargs...
) where {N}
  return findblock(diagblockoffsets(D), block; vargs...)
end

const NonuniformDiagBlockSparse{ElT,VecT} =
  DiagBlockSparse{ElT,VecT} where {VecT<:AbstractVector}
const UniformDiagBlockSparse{ElT,VecT} = DiagBlockSparse{ElT,VecT} where {VecT<:Number}

@propagate_inbounds function getindex(D::NonuniformDiagBlockSparse, i::Int)
  return data(D)[i]
end

getindex(D::UniformDiagBlockSparse, i::Int) = data(D)

@propagate_inbounds function setindex!(D::DiagBlockSparse, val, i::Int)
  data(D)[i] = val
  return D
end

function setindex!(D::UniformDiagBlockSparse, val, i::Int)
  return error("Cannot set elements of a uniform DiagBlockSparse storage")
end

#fill!(D::DiagBlockSparse,v) = fill!(data(D),v)

copy(D::DiagBlockSparse) = DiagBlockSparse(copy(data(D)), copy(diagblockoffsets(D)))

setdata(D::DiagBlockSparse, ndata) = DiagBlockSparse(ndata, diagblockoffsets(D))

## convert to complex
## TODO: this could be a generic TensorStorage function
#complex(D::DiagBlockSparse) = DiagBlockSparse(complex(data(D)), diagblockoffsets(D))

#conj(D::DiagBlockSparse{<:Real}) = D
#conj(D::DiagBlockSparse{<:Complex}) = DiagBlockSparse(conj(data(D)), diagblockoffsets(D))

# TODO: make this generic for all storage types
eltype(::DiagBlockSparse{ElT}) where {ElT} = ElT
eltype(::Type{<:DiagBlockSparse{ElT}}) where {ElT} = ElT

# Deal with uniform DiagBlockSparse conversion
#convert(::Type{<:DiagBlockSparse{ElT,VecT}},D::DiagBlockSparse) where {ElT,VecT} = DiagBlockSparse(convert(VecT,data(D)))

size(D::DiagBlockSparse) = size(data(D))

# TODO: write in terms of ::Int, not inds
similar(D::NonuniformDiagBlockSparse) = setdata(D, similar(data(D)))

similar(D::NonuniformDiagBlockSparse, ::Type{S}) where {S} = setdata(D, similar(data(D), S))
#similar(D::NonuniformDiagBlockSparse,inds) = DiagBlockSparse(similar(data(D),minimum(dims(inds))), diagblockoffsets(D))
#function similar(D::Type{<:NonuniformDiagBlockSparse{ElT,VecT}},inds) where {ElT,VecT}
#  return DiagBlockSparse(similar(VecT,diaglength(inds)), diagblockoffsets(D))
#end

similar(D::UniformDiagBlockSparse) = setdata(D, zero(eltype(D)))
similar(D::UniformDiagBlockSparse, inds) = similar(D)
function similar(::Type{<:UniformDiagBlockSparse{ElT}}, inds) where {ElT}
  return DiagBlockSparse(zero(ElT), diagblockoffsets(D))
end

similar(D::DiagBlockSparse, n::Int) = setdata(D, similar(data(D), n))

function similar(D::DiagBlockSparse, ::Type{ElR}, n::Int) where {ElR}
  return setdata(D, similar(data(D), ElR, n))
end

# TODO: make this work for other storage besides Vector
function zeros(::Type{<:NonuniformDiagBlockSparse{ElT}}, dim::Int64) where {ElT}
  return DiagBlockSparse(zeros(ElT, dim))
end
function zeros(::Type{<:UniformDiagBlockSparse{ElT}}, dim::Int64) where {ElT}
  return DiagBlockSparse(zero(ElT))
end

#
# Type promotions involving DiagBlockSparse
# Useful for knowing how conversions should work when adding and contracting
#

function promote_rule(
  ::Type{<:UniformDiagBlockSparse{ElT1}}, ::Type{<:UniformDiagBlockSparse{ElT2}}
) where {ElT1,ElT2}
  ElR = promote_type(ElT1, ElT2)
  return DiagBlockSparse{ElR,ElR}
end

function promote_rule(
  ::Type{<:NonuniformDiagBlockSparse{ElT1,VecT1}},
  ::Type{<:NonuniformDiagBlockSparse{ElT2,VecT2}},
) where {ElT1,VecT1<:AbstractVector,ElT2,VecT2<:AbstractVector}
  ElR = promote_type(ElT1, ElT2)
  VecR = promote_type(VecT1, VecT2)
  return DiagBlockSparse{ElR,VecR}
end

# This is an internal definition, is there a more general way?
#promote_type(::Type{Vector{ElT1}},
#                  ::Type{ElT2}) where {ElT1<:Number,
#                                       ElT2<:Number} = Vector{promote_type(ElT1,ElT2)}
#
#promote_type(::Type{ElT1},
#                  ::Type{Vector{ElT2}}) where {ElT1<:Number,
#                                               ElT2<:Number} = promote_type(Vector{ElT2},ElT1)

# TODO: how do we make this work more generally for T2<:AbstractVector{S2}?
# Make a similartype(AbstractVector{S2},T1) -> AbstractVector{T1} function?
function promote_rule(
  ::Type{<:UniformDiagBlockSparse{ElT1,VecT1}},
  ::Type{<:NonuniformDiagBlockSparse{ElT2,Vector{ElT2}}},
) where {ElT1,VecT1<:Number,ElT2}
  ElR = promote_type(ElT1, ElT2)
  VecR = Vector{ElR}
  return DiagBlockSparse{ElR,VecR}
end

function promote_rule(
  ::Type{BlockSparseT1}, ::Type{<:NonuniformDiagBlockSparse{ElT2,VecT2,N2}}
) where {BlockSparseT1<:BlockSparse,ElT2<:Number,VecT2<:AbstractVector,N2}
  return promote_type(BlockSparseT1, BlockSparse{ElT2,VecT2,N2})
end

function promote_rule(
  ::Type{BlockSparseT1}, ::Type{<:UniformDiagBlockSparse{ElT2,ElT2}}
) where {BlockSparseT1<:BlockSparse,ElT2<:Number}
  return promote_type(BlockSparseT1, ElT2)
end

# Convert a DiagBlockSparse storage type to the closest Dense storage type
dense(::Type{<:NonuniformDiagBlockSparse{ElT,VecT}}) where {ElT,VecT} = Dense{ElT,VecT}
dense(::Type{<:UniformDiagBlockSparse{ElT}}) where {ElT} = Dense{ElT,Vector{ElT}}

const DiagBlockSparseTensor{ElT,N,StoreT,IndsT} =
  Tensor{ElT,N,StoreT,IndsT} where {StoreT<:DiagBlockSparse}
const NonuniformDiagBlockSparseTensor{ElT,N,StoreT,IndsT} =
  Tensor{ElT,N,StoreT,IndsT} where {StoreT<:NonuniformDiagBlockSparse}
const UniformDiagBlockSparseTensor{ElT,N,StoreT,IndsT} =
  Tensor{ElT,N,StoreT,IndsT} where {StoreT<:UniformDiagBlockSparse}

function DiagBlockSparseTensor(
  ::Type{ElT}, ::UndefInitializer, blocks::Vector, inds
) where {ElT}
  blockoffsets, nnz = diagblockoffsets(blocks, inds)
  storage = DiagBlockSparse(ElT, undef, blockoffsets, nnz)
  return tensor(storage, inds)
end

function DiagBlockSparseTensor(::UndefInitializer, blocks::Vector, inds)
  return DiagBlockSparseTensor(Float64, undef, blocks, inds)
end

function DiagBlockSparseTensor(::Type{ElT}, blocks::Vector, inds) where {ElT}
  blockoffsets, nnz = diagblockoffsets(blocks, inds)
  storage = DiagBlockSparse(ElT, blockoffsets, nnz)
  return tensor(storage, inds)
end

DiagBlockSparseTensor(blocks::Vector, inds) = DiagBlockSparseTensor(Float64, blocks, inds)

# Uniform case
function DiagBlockSparseTensor(x::Number, blocks::Vector, inds)
  blockoffsets, nnz = diagblockoffsets(blocks, inds)
  storage = DiagBlockSparse(x, blockoffsets)
  return tensor(storage, inds)
end

diagblockoffsets(T::DiagBlockSparseTensor) = diagblockoffsets(storage(T))

"""
blockview(T::DiagBlockSparseTensor, block::Block)

Given a block in the block-offset list, return a Diag Tensor
that is a view to the data in that block (to avoid block lookup if the position
is known already).
"""
function blockview(T::DiagBlockSparseTensor, blockT::Block)
  return blockview(T, BlockOffset(blockT, offset(T, blockT)))
end

getindex(T::DiagBlockSparseTensor, block::Block) = blockview(T, block)

function blockview(T::DiagBlockSparseTensor, bof::BlockOffset)
  blockT, offsetT = bof
  blockdimsT = blockdims(T, blockT)
  blockdiaglengthT = minimum(blockdimsT)
  dataTslice = @view data(storage(T))[(offsetT + 1):(offsetT + blockdiaglengthT)]
  return tensor(Diag(dataTslice), blockdimsT)
end

function blockview(T::UniformDiagBlockSparseTensor, bof::BlockOffset)
  blockT, offsetT = bof
  blockdimsT = blockdims(T, blockT)
  return tensor(Diag(getdiagindex(T, 1)), blockdimsT)
end

IndexStyle(::Type{<:DiagBlockSparseTensor}) = IndexCartesian()

# TODO: this needs to be better (promote element type, check order compatibility,
# etc.
function convert(
  ::Type{<:DenseTensor{ElT,N}}, T::DiagBlockSparseTensor{ElT,N}
) where {ElT<:Number,N}
  return dense(T)
end

# These are rules for determining the output of a pairwise contraction of NDTensors
# (given the indices of the output tensors)
function contraction_output_type(
  TensorT1::Type{<:DiagBlockSparseTensor}, TensorT2::Type{<:BlockSparseTensor}, IndsR::Type
)
  return similartype(promote_type(TensorT1, TensorT2), IndsR)
end

function contraction_output_type(
  TensorT1::Type{<:BlockSparseTensor}, TensorT2::Type{<:DiagBlockSparseTensor}, IndsR::Type
)
  return contraction_output_type(TensorT2, TensorT1, IndsR)
end

# This performs the logic that DiagBlockSparseTensor*DiagBlockSparseTensor -> DiagBlockSparseTensor if it is not an outer
# product but -> DenseTensor if it is
# TODO: if the tensors are both order 2 (or less), or if there is an Index replacement,
# then they remain diagonal. Should we limit DiagBlockSparseTensor*DiagBlockSparseTensor to cases that
# result in a DiagBlockSparseTensor, for efficiency and type stability? What about a general
# SparseTensor result?
function contraction_output_type(
  TensorT1::Type{<:DiagBlockSparseTensor{<:Number,N1}},
  TensorT2::Type{<:DiagBlockSparseTensor{<:Number,N2}},
  IndsR::Type,
) where {N1,N2}
  if ValLength(IndsR) === Val{N1 + N2}
    # Turn into is_outer(inds1,inds2,indsR) function?
    # How does type inference work with arithmatic of compile time values?
    return similartype(dense(promote_type(TensorT1, TensorT2)), IndsR)
  end
  return similartype(promote_type(TensorT1, TensorT2), IndsR)
end

# The output must be initialized as zero since it is sparse, cannot be undefined
function contraction_output(T1::DiagBlockSparseTensor, T2::Tensor, indsR)
  return zero_contraction_output(T1, T2, indsR)
end
function contraction_output(T1::Tensor, T2::DiagBlockSparseTensor, indsR)
  return contraction_output(T2, T1, indsR)
end

function contraction_output(T1::DiagBlockSparseTensor, T2::DiagBlockSparseTensor, indsR)
  return zero_contraction_output(T1, T2, indsR)
end

function array(T::DiagBlockSparseTensor{ElT,N}) where {ElT,N}
  return array(dense(T))
end
matrix(T::DiagBlockSparseTensor{<:Number,2}) = array(T)
vector(T::DiagBlockSparseTensor{<:Number,1}) = array(T)

function Array{ElT,N}(T::DiagBlockSparseTensor{ElT,N}) where {ElT,N}
  return array(T)
end

function Array(T::DiagBlockSparseTensor{ElT,N}) where {ElT,N}
  return Array{ElT,N}(T)
end

# Needed to get slice of DiagBlockSparseTensor like T[1:3,1:3]
function similar(
  T::DiagBlockSparseTensor{<:Number,N}, ::Type{ElR}, inds::Dims{N}
) where {ElR<:Number,N}
  return tensor(similar(storage(T), ElR, minimum(inds)), inds)
end

getdiagindex(T::DiagBlockSparseTensor{<:Number}, ind::Int) = storage(T)[ind]

# XXX: handle case of missing diagonal blocks
function setdiagindex!(T::DiagBlockSparseTensor{<:Number}, val, ind::Int)
  storage(T)[ind] = val
  return T
end

function setdiag(T::DiagBlockSparseTensor, val, ind::Int)
  return tensor(DiagBlockSparse(val), inds(T))
end

@propagate_inbounds function getindex(
  T::DiagBlockSparseTensor{ElT,N}, inds::Vararg{Int,N}
) where {ElT,N}
  if all(==(inds[1]), inds)
    return storage(T)[inds[1]]
  else
    return zero(eltype(ElT))
  end
end

@propagate_inbounds function getindex(T::DiagBlockSparseTensor{<:Number,1}, ind::Int)
  return storage(T)[ind]
end

@propagate_inbounds function getindex(T::DiagBlockSparseTensor{<:Number,0})
  return storage(T)[1]
end

# Set diagonal elements
# Throw error for off-diagonal
@propagate_inbounds function setindex!(
  T::DiagBlockSparseTensor{<:Number,N}, val, inds::Vararg{Int,N}
) where {N}
  all(==(inds[1]), inds) ||
    error("Cannot set off-diagonal element of DiagBlockSparse storage")
  storage(T)[inds[1]] = val
  return T
end

@propagate_inbounds function setindex!(T::DiagBlockSparseTensor{<:Number,1}, val, ind::Int)
  storage(T)[ind] = val
  return T
end

@propagate_inbounds function setindex!(T::DiagBlockSparseTensor{<:Number,0}, val)
  storage(T)[1] = val
  return T
end

function setindex!(
  T::UniformDiagBlockSparseTensor{<:Number,N}, val, inds::Vararg{Int,N}
) where {N}
  return error("Cannot set elements of a uniform DiagBlockSparse storage")
end

# TODO: make a fill!! that works for uniform and non-uniform
#fill!(T::DiagBlockSparseTensor,v) = fill!(storage(T),v)

function dense(
  ::Type{<:Tensor{ElT,N,StoreT,IndsT}}
) where {ElT,N,StoreT<:DiagBlockSparse,IndsT}
  return Tensor{ElT,N,dense(StoreT),IndsT}
end

# convert to Dense
function dense(T::TensorT) where {TensorT<:DiagBlockSparseTensor}
  R = zeros(dense(TensorT), inds(T))
  for i in 1:diaglength(T)
    setdiagindex!(R, getdiagindex(T, i), i)
  end
  return R
end

# convert to BlockSparse
function denseblocks(D::Tensor)
  nzblocksD = nzblocks(D)
  T = BlockSparseTensor(eltype(D), nzblocksD, inds(D))
  for b in nzblocksD
    for n in 1:diaglength(D)
      setdiagindex!(T, getdiagindex(D, n), n)
    end
  end
  return T
end

function outer!(
  R::DenseTensor{<:Number,NR},
  T1::DiagBlockSparseTensor{<:Number,N1},
  T2::DiagBlockSparseTensor{<:Number,N2},
) where {NR,N1,N2}
  for i1 in 1:diaglength(T1), i2 in 1:diaglength(T2)
    indsR = CartesianIndex{NR}(ntuple(r -> r ≤ N1 ? i1 : i2, Val(NR)))
    R[indsR] = getdiagindex(T1, i1) * getdiagindex(T2, i2)
  end
  return R
end

# TODO: write an optimized version of this?
function outer!(R::DenseTensor{ElR}, T1::DenseTensor, T2::DiagBlockSparseTensor) where {ElR}
  R .= zero(ElR)
  outer!(R, T1, dense(T2))
  return R
end

function outer!(R::DenseTensor{ElR}, T1::DiagBlockSparseTensor, T2::DenseTensor) where {ElR}
  R .= zero(ElR)
  outer!(R, dense(T1), T2)
  return R
end

# Right an in-place version
function outer(
  T1::DiagBlockSparseTensor{ElT1,N1}, T2::DiagBlockSparseTensor{ElT2,N2}
) where {ElT1,ElT2,N1,N2}
  indsR = unioninds(inds(T1), inds(T2))
  R = tensor(Dense(zeros(promote_type(ElT1, ElT2), dim(indsR))), indsR)
  outer!(R, T1, T2)
  return R
end

function permutedims!(
  R::DiagBlockSparseTensor{<:Number,N},
  T::DiagBlockSparseTensor{<:Number,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {N}
  # TODO: check that inds(R)==permute(inds(T),perm)?
  for i in 1:diaglength(R)
    @inbounds setdiagindex!(R, f(getdiagindex(R, i), getdiagindex(T, i)), i)
  end
  return R
end

function permutedims(
  T::DiagBlockSparseTensor{<:Number,N}, perm::NTuple{N,Int}, f::Function=identity
) where {N}
  R = similar(T, permute(inds(T), perm))
  permutedims!(R, T, perm, f)
  return R
end

function permutedims(
  T::UniformDiagBlockSparseTensor{ElT,N}, perm::NTuple{N,Int}, f::Function=identity
) where {ElR,ElT,N}
  R = tensor(DiagBlockSparse(f(getdiagindex(T, 1))), permute(inds(T), perm))
  return R
end

# Version that may overwrite in-place or may return the result
function permutedims!!(
  R::NonuniformDiagBlockSparseTensor{<:Number,N},
  T::NonuniformDiagBlockSparseTensor{<:Number,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {N}
  RR = convert(promote_type(typeof(R), typeof(T)), R)
  permutedims!(RR, T, perm, f)
  return RR
end

function permutedims!!(
  R::UniformDiagBlockSparseTensor{ElR,N},
  T::UniformDiagBlockSparseTensor{ElT,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {ElR,ElT,N}
  RR = convert(promote_type(typeof(R), typeof(T)), R)
  RR = tensor(DiagBlockSparse(f(getdiagindex(RR, 1), getdiagindex(T, 1))), inds(RR))
  return RR
end

function permutedims!(
  R::DenseTensor{ElR,N},
  T::DiagBlockSparseTensor{ElT,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {ElR,ElT,N}
  for i in 1:diaglength(T)
    @inbounds setdiagindex!(R, f(getdiagindex(R, i), getdiagindex(T, i)), i)
  end
  return R
end

function permutedims!!(
  R::DenseTensor{ElR,N},
  T::DiagBlockSparseTensor{ElT,N},
  perm::NTuple{N,Int},
  f::Function=(r, t) -> t,
) where {ElR,ElT,N}
  permutedims!(R, T, perm, f)
  return R
end

function _contract!!(
  R::UniformDiagBlockSparseTensor{ElR,NR},
  labelsR,
  T1::UniformDiagBlockSparseTensor{<:Number,N1},
  labelsT1,
  T2::UniformDiagBlockSparseTensor{<:Number,N2},
  labelsT2,
) where {ElR,NR,N1,N2}
  if NR == 0  # If all indices of A and B are contracted
    # all indices are summed over, just add the product of the diagonal
    # elements of A and B
    R = setdiag(R, diaglength(T1) * getdiagindex(T1, 1) * getdiagindex(T2, 1), 1)
  else
    # not all indices are summed over, set the diagonals of the result
    # to the product of the diagonals of A and B
    R = setdiag(R, getdiagindex(T1, 1) * getdiagindex(T2, 1), 1)
  end
  return R
end

function contraction_output(
  T1::TensorT1, labelsT1, T2::TensorT2, labelsT2, labelsR
) where {TensorT1<:BlockSparseTensor,TensorT2<:DiagBlockSparseTensor}
  indsR = contract_inds(inds(T1), labelsT1, inds(T2), labelsT2, labelsR)
  TensorR = contraction_output_type(TensorT1, TensorT2, typeof(indsR))
  blockoffsetsR, contraction_plan = contract_blockoffsets(
    blockoffsets(T1),
    inds(T1),
    labelsT1,
    blockoffsets(T2),
    inds(T2),
    labelsT2,
    indsR,
    labelsR,
  )
  R = zeros(TensorR, blockoffsetsR, indsR)
  return R, contraction_plan
end

function contract(
  T1::BlockSparseTensor,
  labelsT1,
  T2::DiagBlockSparseTensor,
  labelsT2,
  labelsR=contract_labels(labelsT1, labelsT2),
)
  R, contraction_plan = contraction_output(T1, labelsT1, T2, labelsT2, labelsR)
  R = contract!(R, labelsR, T1, labelsT1, T2, labelsT2, contraction_plan)
  return R
end

function contract(
  T1::DiagBlockSparseTensor,
  labelsT1,
  T2::BlockSparseTensor,
  labelsT2,
  labelsR=contract_labels(labelsT2, labelsT1),
)
  return contract(T2, labelsT2, T1, labelsT1, labelsR)
end

function contract!(
  R::BlockSparseTensor{ElR,NR},
  labelsR,
  T1::BlockSparseTensor,
  labelsT1,
  T2::DiagBlockSparseTensor,
  labelsT2,
  contraction_plan,
) where {ElR<:Number,NR}
  already_written_to = Dict{Block{NR},Bool}()
  # In R .= α .* (T1 * T2) .+ β .* R
  α = one(ElR)
  for (block1, block2, blockR) in contraction_plan
    T1block = T1[block1]
    T2block = T2[block2]
    Rblock = R[blockR]

    # <fermions>
    α = compute_alpha(
      ElR, labelsR, blockR, inds(R), labelsT1, block1, inds(T1), labelsT2, block2, inds(T2)
    )

    β = one(ElR)
    if !haskey(already_written_to, blockR)
      already_written_to[blockR] = true
      # Overwrite the block of R
      β = zero(ElR)
    end
    contract!(Rblock, labelsR, T1block, labelsT1, T2block, labelsT2, α, β)
  end
  return R
end

function contract!(
  C::BlockSparseTensor,
  Clabels,
  A::BlockSparseTensor,
  Alabels,
  B::DiagBlockSparseTensor,
  Blabels,
)
  return contract!(C, Clabels, B, Blabels, A, Alabels)
end

function show(io::IO, mime::MIME"text/plain", T::DiagBlockSparseTensor)
  summary(io, T)
  for (n, block) in enumerate(keys(diagblockoffsets(T)))
    blockdimsT = blockdims(T, block)
    println(io, block)
    println(io, " [", _range2string(blockstart(T, block), blockend(T, block)), "]")
    print_tensor(io, blockview(T, block))
    n < nnzblocks(T) && print(io, "\n\n")
  end
end

show(io::IO, T::DiagBlockSparseTensor) = show(io, MIME("text/plain"), T)


text:
####
#### generic API
####

export MinimizationProblem, local_minimization, multistart_minimization

struct MinimizationProblem{F,T<:AbstractVector}
    "The function to be minimized."
    objective::F
    "Lower bounds (a vector of real numbers)."
    lower_bounds::T
    "Upper bounds (a vector of real numbers)."
    upper_bounds::T
    @doc """
    $(SIGNATURES)

    Define a minimization problem.

    - `objective` is a function that maps `ℝᴺ` vectors to real numbers. It should accept all
       `AbstractVector`s of length `N`.
    - `lower_bounds` and `upper_bounds` define the bounds, and their lengths implicitly
      define the dimension `N`.

    The fields of the result should be accessible with the names above.
    """
    function MinimizationProblem(objective::F, lower_bounds::T, upper_bounds::T) where {F,T}
        @argcheck all(lower_bounds .< upper_bounds)
        new{F,T}(objective, lower_bounds, upper_bounds)
    end
end

"""
`$(FUNCTIONNAME)(minimization_problem, local_method, x)`

Perform a local minimization using `local_method`, starting from `x`.

The objective and the bounds are provided in `minimization_problem`, see
[`MinimizationProblem`](@ref). `x` should be an `AbstractVector` of conforming dimension.

`local_method` can be a type for which a method is defined (recommended). However, it can
also be a closure, in which case it will be called as `local_method(minimization_problem, x)`.

In both cases, it should return `nothing`, or a value which has the following properties:

- `location`, an `AbstractVector` for the minimizer.

- `value`, for the value of the objective at `location`. `Inf` (or equivalent) should be
  used for infeasible regions.

The returned value may have other properties too, these are useful for convergence
diagnostics, debugging information, etc, these depend on the `local_method`.

Returning `nothing` is equivalent to `value == Inf`, but in some cases can work better for
type inference as the method won't have to construct a counterfactual type.
"""
function local_minimization(local_method, minimization_problem::MinimizationProblem, x)
    local_method(minimization_problem, x)
end

"""
`$(FUNCTIONNAME)(multistart_method, local_method, minimization_problem)`

Multistart minimization using the given multistart and local methods.
"""
function multistart_minimization end


text:
module TestExtra
using Base.Test

export runtests

function testfiles(dir=".")
    filter(f -> ismatch(r"^[^\.].*\.jl$",f) && !(f == "runtests.jl"), readdir(dir))
end

function testnworkers()
    if haskey(ENV,"JULIA_TEST_NWORKERS")
        parse(Int, ENV["JULIA_TEST_NWORKERS"])
    else
        Sys.CPU_CORES
    end
end

function runtests(tests::AbstractVector, nworkers::Int=testnworkers(); verbose=true)
    verbose && println("Testing with $nworkers workers")
    @testset "TOTAL" begin
        mts = Base.Test.get_testset()
        @sync for p in 1:nworkers
            @async while !isempty(tests)
                Base.Test.push_testset(mts)
                test = shift!(tests)
                verbose && println("Running $test")
                pid = addprocs(1)[1]
                try
                    remotecall_fetch(Main.eval, pid, :(using Base.Test))
                    rts = remotecall_fetch(Main.eval, pid, quote
                                     ts = Base.Test.DefaultTestSet($test)
                                     Base.Test.push_testset(ts)
                                     Base.Test.get_testset()
                                     try
                                         include($test)
                                     catch err
                                         Base.Test.record(ts, Error(:nontest_error, :(), err, Base.Test.catch_backtrace()))
                                     end
                                     Base.Test.pop_testset()
                                     ts
                    end)
                    Base.Test.finish(rts)
                finally
                    rmprocs(pid)
                end
            end
        end
    end
end

runtests(nworkers::Int=testnworkers(); kwargs...) = runtests(testfiles(), nworkers; kwargs...)

end # module


text:
"""
    qzero(N, dims=(N,))

Generate a zero operator for a Hilbert space of size `N`. It is possible to specify the subspace dimensions with the `dims` argument. Returns a sparse matrix.

# Example
```jldoctest
julia> qzero(4,(2,2))
4×4 Operator{SparseMatrixCSC{Float64,Int64},2} with dimensions 2⊗2
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
```
"""
function qzero(::Type{T}, N::Integer, dims::Dims=(N,)) where {T<:Number}
    return Operator(SparseMatrixCSC{T,Int}(N,N,ones(Int,N+1),Int[],T[]),dims)
end
qzero(::Type{T}, dims::Dims) where {T<:Number} = qzero(T,prod(dims),dims)
qzero(N::Integer, dims::Dims=(N,)) = qzero(Float64,N,dims)
qzero(dims::Dims) = qzero(Float64,dims)

"""
    qeye(N, dims=(N,))

Generate an identity operator for a Hilbert space of size `N`. It is possible to specify the subspace dimensions with the `dims` argument. Returns a sparse matrix.

# Example
```jldoctest
julia> qeye(4,(2,2))
4×4 Operator{SparseMatrixCSC{Float64,Int64},2} with dimensions 2⊗2
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
```
"""
function qeye(::Type{T}, N::Integer, dims::Dims=(N,)) where {T<:Number}
    return Operator(SparseMatrixCSC{T,Int}(N,N,[1:N+1;],[1:N;],ones(T,N)),dims)
end
qeye(::Type{T}, dims::Dims) where {T<:Number} = qeye(T,prod(dims),dims)
qeye(N::Integer, dims::Dims=(N,)) = qeye(Float64,N,dims)
qeye(dims::Dims) = qeye(Float64,prod(dims),dims)

"""
    destroy(N)

Generate a quantum harmonic oscillator lowering (annihilation) operator ``\\hat{a}`` in a truncated Hilbert space of size `N`. Returns a sparse matrix.

# Example
```jldoctest
julia> destroy(4)
4×4 Operator{SparseMatrixCSC{Float64,Int64},1} with dimensions 4
 0.0  1.0  0.0      0.0
 0.0  0.0  1.41421  0.0
 0.0  0.0  0.0      1.73205
 0.0  0.0  0.0      0.0
```
"""
function destroy(::Type{T}, N::Integer) where {T<:Number}
    colptr = fill(1,N+1); @inbounds for i in 1:N; colptr[i+1]=i; end
    nzval = T[√(T(i)) for i in 1:N-1]
    return Operator(SparseMatrixCSC{T,Int}(N,N,colptr,[1:N-1;],nzval),(N,))
end
destroy(N::Integer) = destroy(Float64,N)

"""
    create(N)

Generate a quantum harmonic oscillator raising (creation) operator ``\\hat{a}^†`` in a truncated Hilbert space of size `N`. Returns a sparse matrix.

# Example
```jldoctest
julia> create(4)
4×4 Operator{SparseMatrixCSC{Float64,Int64},1} with dimensions 4
 0.0  0.0      0.0      0.0
 1.0  0.0      0.0      0.0
 0.0  1.41421  0.0      0.0
 0.0  0.0      1.73205  0.0
```
"""
function create(::Type{T}, N::Integer) where {T<:Number}
    colptr = fill(N,N+1); @inbounds for i in 1:N; colptr[i]=i; end
    nzval = T[√(T(i)) for i in 1:N-1]
    return Operator(SparseMatrixCSC{T,Int}(N,N,colptr,[2:N;],nzval),(N,))
end
create(N::Integer) = create(Float64,N)

"""
    numberop(N)

Generate a number operator ``\\hat{n}`` in a Hilbert space of size `N`. Returns a sparse matrix.

# Example
```jldoctest
julia> numberop(4)
4×4 Operator{SparseMatrixCSC{Float64,Int64},1} with dimensions 4
 0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  2.0  0.0
 0.0  0.0  0.0  3.0
```
"""
function numberop(::Type{T}, N::Integer) where {T<:Number}
    # nzval includes a structural 0 for the [1,1] entry
    nzval = T[0:N;]
    return Operator(SparseMatrixCSC{T,Int}(N,N,[1:N+1;],[1:N;],nzval),(N,))
end
numberop(N::Integer) = numberop(Float64,N)

"""
    displacementop(N, α)

Generate a quantum harmonic oscillator displacement operator ``\\hat{D}(α)`` in a truncated Hilbert space of size `N`. Returns a dense matrix.

```math
\\hat{D}(α) = \\exp\\left(α\\hat{a}^† - α^*\\hat{a}\\right)
```

# Example
```jldoctest
julia> displacementop(3,0.5im)
3×3 Operator{Array{Complex{Float64},2},1} with dimensions 3
   0.88262+0.0im            0.0+0.439802im  -0.166001+0.0im
       0.0+0.439802im  0.647859+0.0im             0.0+0.621974im
 -0.166001+0.0im            0.0+0.621974im    0.76524+0.0im
```
"""
function displacementop(N::Integer, α::Number)
    a = Array(destroy(N))
    return Operator(LinearAlgebra.exp!(α.*a' .- α'.*a),(N,))
end

"""
    squeezeop(N, z)

Generate a quantum harmonic oscillator squeeze operator ``\\hat{S}(z)`` in a truncated Hilbert space of size `N`. Returns a dense matrix.

```math
\\hat{S}(z) = \\exp\\left(\\frac{1}{2}\\left(z^*\\hat{a}^2 - z\\hat{a}^{†2}\\right)\\right)
```

# Example
```jldoctest
julia> squeezeop(3,0.5im)
3×3 Operator{Array{Complex{Float64},2},1} with dimensions 3
 0.938148-0.0im       0.0-0.0im       0.0-0.346234im
      0.0+0.0im       1.0+0.0im       0.0+0.0im
      0.0-0.346234im  0.0-0.0im  0.938148-0.0im
```
"""
function squeezeop(N::Integer, z::Number)
    a = Array(destroy(N))
    return Operator(LinearAlgebra.exp!((z'.*a^2 .- z.*a'^2)./2),(N,))
end

"""
    projectorop(N,S)

Generate a projector on the subspaces defined by an integer or a vector/range of integers `S`:

```math
P = ∑_{i∈S} |i⟩⟨i|.
```

# Example
```jldoctest
julia> projectorop(5,[1,3])
5×5 Operator{SparseMatrixCSC{Float64,Int64},1} with dimensions 5
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0
```
"""
function projectorop(N::Integer,S::AbstractVector{<:Integer})
    maximum(S)<N || throw(ArgumentError("a $N-d space cannot be projected on level $(maximum(S))"))
    I = S.+1
    V = ones(length(S))
    return Operator(sparse(I,I,V,N,N),(N,))
end
projectorop(N::Integer,S::Integer) = projectorop(N,S:S)

"""
    sylvesterop(N,k,l)

Generate the ``(k,j)^\\text{th}`` Sylvester generalized Pauli matrix in N-d.
https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices
"""
function sylvesterop(N::Integer,k::Integer,j::Integer)
    ωʲ = Complex(cospi(2j/N),sinpi(2j/N))
    rowval = [mod1(i+k,N) for i in 1:N]
    nzval  = [ωʲ^m for m = 0:N-1]
    return Operator(SparseMatrixCSC(N,N,[1:N+1;],rowval,nzval),(N,))
end

"""
    Sigma1(N)

Generate the N-d shift matrix ``Σ₁``.
"""
function Sigma1(N)
    rowval = circshift(1:N,-1)
    nzval  = ones(N)
    return Operator(SparseMatrixCSC(N,N,[1:N+1;],rowval,nzval),(N,))
end

"""
    Sigma3(N)

Generate the N-d clock matrix ``Σ₃``.
"""
function Sigma3(N)
    ω = Complex(cospi(2/N),sinpi(2/N))
    rowval = [1:N;]
    nzval  = [ω^m for m = 0:N-1]
    return Operator(SparseMatrixCSC(N,N,[1:N+1;],rowval,nzval),(N,))
end


text:
if !isdefined(Main,:IsReduced); IsReduced = FastRationals.IsReduced; end;
if !isdefined(Main,:MayReduce); MayReduce = FastRationals.MayReduce; end;
if !isdefined(Main,:basistype); basistype = FastRationals.basistype; end;

function sums_toward_half(::Type{T}, n; details::Bool=false) where {T}
    one_third = onethird(T)
    one_ninth = one_third * one_third
    s = one_third + one_ninth
    t = one_third * one_ninth

    n = abs(n)
    while n>1
     n -= 1
     s += t
     t *= one_third
    end

    return !details ? s : (t, s) # debugging info
end

function onethird(::Type{Q}) where Q
    T = basistype(Q)
    return Q(one(T), one(T)+one(T)+one(T))
end

slow17 = sums_toward_half(Rational{Int64},17)
# (1743392200//3486784401, 1//3486784401)

fast17 = sums_toward_half(FastRational{Int64, IsReduced},17)
# (1743392200//3486784401, 1//3486784401)


text:
using Gadfly, DataFrames, DataArrays
df = readtable("C:/Users/Peter.Mancini/Documents/Datafiles/playersessiongap_post.csv")

light_panel = Theme(
    panel_fill="white",
    default_color="darkblue",
    panel_stroke="darkblue",
    background_color="old lace")

hmdata3 = by(df, [:hour], x -> [mean(x[:similarity]) std(x[:similarity])])
Ns = by(df, [:hour], x-> size(x)[1])
ymins = hmdata3[:x1] .- (1.96 * Array(hmdata3[:x2]) ./ Array(sqrt(Ns[:x1])))
ymaxs = hmdata3[:x1] .+ (1.96 * Array(hmdata3[:x2]) ./ Array(sqrt(Ns[:x1])))

#l1 = layer(hmdata[:x1], background)
p3 = plot(x=0:23, y=hmdata3[:x1], ymin=ymins, ymax = ymaxs,
  Guide.title("Mean likelihood of stickiness by hour"), Guide.xlabel("Hour"),
  Guide.ylabel("Percent"),
  Geom.point, Geom.smooth(method=:loess, smoothing=0.9), Geom.errorbar, light_panel)

draw(SVG("Mean likelihood of stickiness by hour with error.svg", 10inch, 5inch), p3)
draw(PNG("Mean likelihood of stickiness by hour with error.png", 10inch, 5inch), p3)

hmdata4 = by(df, [:hour], x -> [mean(x[:gap]) std(x[:gap])])
Ns = by(df, [:hour], x-> size(x)[1])
ymins = hmdata4[:x1] .- (1.96 * Array(hmdata4[:x2]) ./ Array(sqrt(Ns[:x1])))
ymaxs = hmdata4[:x1] .+ (1.96 * Array(hmdata4[:x2]) ./ Array(sqrt(Ns[:x1])))
p4 = plot(x=0:23, y=hmdata4[:x1], ymin=ymins, ymax = ymaxs,
  Guide.title("Mean gap between microsessions by hour"), Guide.xlabel("Hour"),
  Guide.ylabel("Seconds"),
  Geom.point, Geom.smooth(method=:loess,smoothing=0.45), Geom.errorbar, light_panel)

draw(SVG("Mean gap between microsessions by hour with error.svg", 10inch, 5inch), p4)
draw(PNG("Mean gap between microsessions by hour with error.png", 10inch, 5inch), p4)


text:
using StaticArrays
include("PeriodicNumber.jl")

#=
Definition of norm of vectors whose elements are PeriodicNumbers.
Note that in this way, the same function can be applied to usual 'Vector' arrays, but also to 'StaticVector' types as will be used afterwards
=#
import LinearAlgebra.norm
"norm function for a vector of `PeriodicNumber` elements"
norm(v::AbstractVector{<:PeriodicNumber}) = norm(value.(v))
norm(PeriodicNumber.(3.0*rand(5), 2.0))

#########################################################################################################
#########################################################################################################
### Functions to create static vectors whose components are of type 'PeriodicNumber'
#########################################################################################################
#########################################################################################################

"""
    PeriodicVector(vec::Vector{T}, L::T=1.0)

Convert 'vec' to a Static Vector of the same size, but with elements of 'PeriodicNumber{T}' type.

The periodicity of the numbers, 'L', is given as second argument and defaults to 1.0.
"""
function PeriodicVector(vec::Vector{T}, L::T=1.0) where {T<:AbstractFloat}
    dim=length(vec)
    if typeof(L)==T
        SVector{dim}(PeriodicNumber.(vec, L))
    else
        SVector{dim}(PeriodicNumber.(vec, convert(T, L)))
    end
end
PeriodicVector(rand(4), 2.0)

"""
    PeriodicVectors(mat::Matrix{T}, L::T=1.0) where {T<:AbstractFloat}

Convert 'mat' to a Vector of elements of type SVector{d, PeriodicNumber{T}}.

mat should be of size d x N. The output is a 1-dimensional array of N elements, each of which
consists of StaticVector's of size d, and PeriodicNumber{T} as data.
The periodicity of the numbers, 'L', is given as second argument and defaults to 1.0.
"""
function PeriodicVectors(mat::Matrix{T}, L::T=1.0) where {T<:AbstractFloat}
    dim, N = size(mat)
    @assert dim<N "d>N. N should be larger than the dimensionality of vectors"
    [PeriodicVector(mat[:,i], L) for i in 1:N]
end
PeriodicVectors(rand(4, 10), 2.0)

#### Format the output of periodic vectors
function Base.show(io::IO, PV::SVector{d, <:PeriodicNumber}) where {d}
    L = PV[1].L
    vals = value.(PV)
    if L==1.0
        print(io, "PV(", vals, ")")
    else
        print(io, "PV(", vals, "; L=", L, ")" )
    end
end

function value(Xs::Vector{SVector{d, PeriodicNumber{T}}}) where {d, T<:AbstractFloat}
    N = length(Xs)
    value.(reshape(reinterpret(PeriodicNumber{T}, Xs), (d, N)))
end
#########################################################################################################
#########################################################################################################
### Functions to compute the distance between two vectors with 'PeriodicNumber' components
#########################################################################################################
#########################################################################################################

"""
    trit_to_shift(x::Char, L::Real=1.0)

Associate a base-3 digit (as a character) to a shift.

The idea is to use a ternary representation (trit is for base 3, what bit is for base 2) to
identify a system's images by shifts.
The association rule is the following:
    0: for no shift
    1: for a shift of size L; so the image to the right
    2: for a shift of size -L; so the image to the left

In Julia 0.7 and later transforming a number to a different base outputs a string, that's why
it is useful to assume that x is type 'Char'
"""
function trit_to_shift(x::Char, L::Float64=1.0)
    if x=='0'
        return 0.0
    elseif x=='1'
        return L
    elseif x=='2'
        return -L
    else
        return error("Wrong digit as input! Only '0', '1', '2' are possible.")
    end
end

"""
    trit_to_shift_vector(s::String, L::T=1.0)

Split 's' into characters and apply 'trit_to_shift' to each of them.

The output is thus a vector of the same number of elements as characters in 's'.
The value of the shift, 'L', defaults to 1.0.

See also: [`trit_to_shift`](@ref)
"""
function trit_to_shift_vector(s::String, L::Float64=1.0)
    trit_to_shift.(collect(s), L)
end

"""
    generate_system_images(d::Int64, L::T=1.0)

Generate the 3^d-1 (virtual) images of a d-dimensional system of periodicity 'L' in each side.

The output is a Vector whose elements are of type SVector{d,T}.
The value of 'L' defaults to 1.0.0
For a d-dimensional system, there are 3^d-1 virtual images. Each of them is identified with a
single number to which a shift vector is associated.

See also: [`trit_to_shift`](@ref), [`trit_to_shift_vector`](@ref)
"""
function generate_system_images(d::Int64, L::Float64=1.0)
    images = Vector{SVector{d, Float64}}(undef, 3^d-1)
    # Note that i does not begin in '0'; so we're excluding the case of no shift (the most common one for big systems)
    for i in 1:3^d-1::Int
        str = string(i, base=3, pad=d)
        images[i] = SVector{d}(trit_to_shift_vector(str, L) )
    end
    images
end
generate_system_images(2, 5.0)

"""
    MIC_distance(V1::SVector{d, PeriodicNumber{T}}, V2::SVector{d, PeriodicNumber{T}}, images::Vector{SVector{d, T}})

Compute the distance between 'V1' and 'V2' considering periodic boundary conditions and using
the so called Minimum Image Convention (MIC).

To avoid unnecessary calls in other functions, the system's set of virtual images should be
given as third argument, as an array of SVector's.
It is assumed that both vectors are contained in a d-dimensional box of size L.
The MIC guarantees that all the possible periodic shifts or 'virtual images' of the relative vectors
are considered when calculating their distance.

The output is the index of the image of minimum distance (0 if no shift is needed) and the value of such distance.

See also: [`MIC_vector`](@ref).
"""
function MIC_distance(V1::SVector{d, PeriodicNumber{T}}, V2::SVector{d, PeriodicNumber{T}}, images::Vector{SVector{d, T}}) where {d, T<:AbstractFloat}
    L1 = V1[1].L; L2 = V2[1].L # periodicities of each vector
    #check that both vectors have the same periodicity
    (L1 != L2) && error("Comparing periodic vectors of different periodicity")

    X1 = value.(V1); X2 = value.(V2) # values of the components of the first and second vector
    diffs =  abs.(X1 .- X2) # absolute value of the difference of each of the vectors components

    if all(x -> x<=L1/2, diffs) # if all the components differ by less than half the box size, there is no need to consider the system images
        return 0, norm(diffs) # the output is the 0-th image (so the system itself), and the usual norm of the vector difference

    else # if at least one of the components is larger than L/2, then check all the distances with all the possible images
        #TODO: Try something smarter, so that not all the 3^d-1 shifts are evaluated, but only the relevant ones
        min_dist, close_image = findmin( [norm(X1 .- (X2.+shift) ) for shift in images])
        return close_image, min_dist # the output is the index of the image that produced the smallest distance, and the value of such distance
    end
end
MIC_distance(PeriodicVector(rand(4), 0.5), PeriodicVector(rand(4), 0.5), generate_system_images(4, 0.5))

"""
    MIC_vector(V1::SVector{d, PeriodicNumber{T}}, V2::SVector{d, PeriodicNumber{T}}, images::Vector{SVector{d, T}})

Compute the vector joining V2 to V1 (i.e. V1-V2) considering periodic boundary conditions
and using the so called Minimum Image Convention (MIC).

To avoid unnecessary calls in other functions, the system's set of virtual images should be
given as third argument, as an array of SVector's.
It is assumed that both vectors are contained in a d-dimensional box of size L.
The MIC guarantees that all the possible periodic shifts or 'virtual images' of the vectors
are considered when calculating their distance.

The output is the MIC vector (as a `SVector` type) and the index of the image of minimum
 distance (0 if no such image is needed).

See also: [`MIC_distance`](@ref).
"""
function MIC_vector(V1::SVector{d, PeriodicNumber{T}}, V2::SVector{d, PeriodicNumber{T}}, images::Vector{SVector{d, T}}) where {d, T<:AbstractFloat}
    L1 = V1[1].L; L2 = V2[1].L # periodicities of each vector
    #check that both vectors have the same periodicity
    (L1 != L2) && error("Comparing periodic vectors of different periodicity")

    X1 = value.(V1); X2 = value.(V2) # values of the components of the first and second vector
    diffs =  abs.(X1 .- X2) # absolute value of the difference of each of the vectors components

    if all(x -> x<=L1/2, diffs) # if all the components differ by less than half the box size, there is no need to consider the system images
        return SVector{d}(X1 .- X2), 0   # the output is the usual vector difference and the 0-th image (so the system itself), and
    else # if at least one of the components is larger than L/2, then check all the distances with all the possible images
        #TODO: Try something smarter, so that not all the 3^d-1 shifts are evaluated, but only the relevant ones
        min_dist, close_image = findmin( [norm(X1 .- (X2.+shift) ) for shift in images])

        # the output is the vector difference using a virtual image that produces the smallest distance, and the index of such image
        return SVector{d}(X1 .- (X2 .+ images[close_image])), close_image
    end
end
MIC_vector(PeriodicVector(rand(4), 0.5), PeriodicVector(rand(4), 0.5), generate_system_images(4, 0.5))

#########################################################################################################
#########################################################################################################
### Functions to check distances between set of vectors
#########################################################################################################
#########################################################################################################

"""
    distances_between_centers(Xs::Vector{SVector{d, PeriodicNumber{T}}})

Compute the distance between all pair of SVector's.

If 'Xs' is an array of length N, the output is a NxN symmetric matrix.
When this function is called each element of 'Xs' plays the role of the center of a particle.

We make use of the 'Symmetric' constructor, so only N*(N-1)/2 distances are calculated.

See also [`distances_between_particles`](@ref)
"""
function distances_between_centers(Xs::Vector{SVector{d, PeriodicNumber{T}}}) where {d, T<:AbstractFloat}
    N = length(Xs)
    distances=zeros(N,N)
    @inbounds for i in 1:N, j in i+1:N
        distances[i,j]=norm(Xs[i] .- Xs[j])
    end
    return Symmetric(distances)
end

function distances_between_centers!(distances::Symmetric{T, Matrix{T}}, Xs::Vector{SVector{d, PeriodicNumber{T}}}, is::Vector{Int64}) where {d, T<:AbstractFloat}
    N = length(Xs)
    dists=distances.data
    @inbounds for i in is, j in i+1:N
        dists[i,j]=norm(Xs[i] .- Xs[j])
    end
    copy!(distances, Symmetric(dists))
end

"""
    distances_between_centers(Xs::Vector{SVector{d, PeriodicNumber{T}}},  images::Vector{SVector{d, T}})

Compute the distance between all pair of SVector's and the image of MIC distance between them.

If 'Xs' is an array of length N, the output is a NxN symmetric matrix containing the distances
and a NxN symmetric matrix containing the MIC images.
When this function is called each element of 'Xs' plays the role of the center of a particle.

We make use of the 'Symmetric' constructor, so only N*(N-1)/2 distances are calculated.

See also: [`MIC_distance`](@ref).
"""
function distances_between_centers(Xs::Vector{SVector{d, PeriodicNumber{T}}}, images::Vector{SVector{d, T}}) where {d, T<:AbstractFloat}
    N = length(Xs)
    distances=zeros(N,N)
    image_indices = zeros(Int64, N, N)
    for i in 1:N, j in i+1:N
        img, dist = MIC_distance(Xs[i], Xs[j], images)
        distances[i,j] = dist
        image_indices[i,j] = img

        if img != 0
            min_val, mirror_img = findmin(norm.( images .+ [images[img]] ))
            image_indices[j, i] = mirror_img
        else
            image_indices[j, i] = 0
        end
    end
    return Symmetric(distances), image_indices
end

"""
    check_for_overlaps(Xs::Vector{SVector{d, PeriodicNumber{T}}}, R::T, tolerance::T)

Check whether there is an overlap between all pairs of particles centred at 'Xs' and of radius
'R'.

Each element of 'Xs' plays the role of the position of a particle's center, and all of them
are assumed to be of the same size, i.e. radius 'R'.
A tolerance to determine whether there is an overlap or not should  be passed as third argument.

# Output

- 'overlap': a boolean that is 'true' only when an overlap is present
- 'message': a string  that contains some information about the overlapping particles.
- 'particles': a tuple containing the indices of the overlapping particles; '(0, 0)' if no overlap
"""
function check_for_overlaps(Xs::Vector{SVector{d, PeriodicNumber{T}}}, R::T, tolerance::T) where {d, T<:AbstractFloat}
    σ = 2*R
    N = length(Xs)
    distances = distances_between_centers(Xs)
    overlap=false; message = "No overlap is present"; particles = (0, 0)
    @inbounds for i in 1:N, j in i+1:N
        if distances[i,j]-σ < -tolerance
            gap = σ - distances[i,j]
            Xi = value.(Xs[i])
            Xj = value.(Xs[j])
            message = "Overlap between particles $i and $j;\n Centers at: $Xi; \t $Xj. \n\t\tradius:  $R; \t Overlap = $gap ."

            overlap = true
            particles = (i, j)
            break
        end
    end
    return overlap, message, particles
end

"""
    check_for_overlaps(Xs::Vector{SVector{d, PeriodicNumber{T}}}, Rs::Vector{T}, tolerance::T)

Check whether there is an overlap between all pairs of particles centred at 'Xs' and of radii
'Rs'.

Each element of 'Xs' plays the role of the position of a particle's center, while its associated
radius is stored in corresponding element of 'Rs'
A tolerance to determine whether there is an overlap or not should  be passed as third argument.

# Output

- 'overlap': a boolean that is 'true' only when an overlap is present
- 'message': a string  that contains some information about the overlapping particles.
- 'particles': a tuple containing the indices of the overlapping particles; '(0, 0)' if no overlap
"""
function check_for_overlaps(Xs::Vector{SVector{d, PeriodicNumber{T}}}, Rs::Vector{T}, tolerance::T) where {d, T<:AbstractFloat}
    N = length(Xs)
    distances = distances_between_centers(Xs)
    overlap=false; message = "No overlap is present"; particles = (0, 0)
    @inbounds for i in 1:N, j in i+1:N
        σ = Rs[i] + Rs[j]
        if distances[i,j]-σ < -tolerance
            gap = σ - distances[i,j]
            Xi = value.(Xs[i]); Ri = Rs[i]
            Xj = value.(Xs[j]); Rj = Rs[j]
            message = "Overlap between particles $i and $j;\n Centers at: $Xi; \t $Xj. \n\t\tradii:  ($Ri, $Rj); \t Overlap = $gap ."

            overlap = true
            particles = (i, j)
            break
        end
    end
    return overlap, message, particles
end

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
function setup_square_corrector_test()
    F = SPSystem(equations(katsura(6)))
    x = rand(Complex{Float64}, 7)
    xnext = copy(x)
    t = rand()
    H = HomotopyWithCache(StraightLineHomotopy(F, F), x, t)
    Jac = HC.Jacobian(HC.jacobian(H, x, t))
    H, x, xnext, t, Jac
end

function setup_overdetermined_corrector_test()
    @polyvar x y
    F = SPSystem([x^2+y, y^2-3x*y, y+x+3, y+2x-5])
    x = rand(Complex{Float64}, 2)
    xnext = copy(x)
    t = rand()
    H = HomotopyWithCache(StraightLineHomotopy(F, F), x, t)
    Jac = HC.Jacobian(HC.jacobian(H, x, t))
    H, x, xnext, t, Jac
end

@testset "Correctors" begin
    @testset "Newton - square" begin
        H, x, xnext, t, Jac = setup_square_corrector_test()

        corrector = NewtonCorrector()
        @test corrector isa NewtonCorrector
        corrector_cache = cache(corrector, H, x, t)
        @test corrector_cache isa HomotopyContinuation.NewtonCorrectorCache

        # check that this doesn't throw
        out = correct!(xnext, corrector, corrector_cache, H, x, t, euclidean_norm, Jac, 1e-7, 3,
        HC.StepSizeModel())
        @test out isa CorrectorResult

        out = correct!(xnext, corrector, corrector_cache, H, x, t,
                    euclidean_norm, Jac, 1e-7, 3, HC.StepSizeModel(); update_jacobian_infos=true)
        @test out isa CorrectorResult
    end

    @testset "Newton - overdetermined" begin
        H, x, xnext, t, Jac = setup_overdetermined_corrector_test()

        corrector = NewtonCorrector()
        @test corrector isa NewtonCorrector
        corrector_cache = cache(corrector, H, x, t)
        @test corrector_cache isa HomotopyContinuation.NewtonCorrectorCache

        # check that this doesn't throw
        out = correct!(xnext, corrector, corrector_cache, H, x, t, euclidean_norm, Jac, 1e-7, 3, HC.StepSizeModel())
        @test out isa CorrectorResult

        out = correct!(xnext, corrector, corrector_cache, H, x, t,
                    euclidean_norm, Jac, 1e-7, 3, HC.StepSizeModel(); update_jacobian_infos=true)
        @test out isa CorrectorResult
    end
end


text:
#**************************************************************************************
# Input_Interface.jl
# =============== part of the GeoEfficiency.jl package.
#
# all the input either from the console or from the csv files to the package is handled by some function here.
#
#**************************************************************************************

#------------------ consts - globals - imports ----------------------------

using Compat
using Compat.MathConstants
using Compat.DelimitedFiles
import Compat: @info, @warn, @error

@compat isconst(@__MODULE__, :dataFolder ) 	||	const dataFolder = string(@__MODULE__)
@compat isconst(@__MODULE__, :dataDir )		||	const  dataDir    = joinpath(homedir(), dataFolder);
mkpath(dataDir)

const Detectors  = "Detectors.csv";
const srcHeights = "srcHeights.csv";
const srcRhos    = "srcRhos.csv";
const srcRadii   = "srcRadii.csv";
const srcLengths = "srcLengths.csv";

@enum SrcType srcUnknown=-1 srcPoint=0 srcLine=1 srcDisk=2 srcVolume=3 srcNotPoint =4
global srcType = srcUnknown

#------------------ typeofSrc --------------------------------------

"""
    typeofSrc()::SrcType

return the current value of the global `GeoEfficiency.srcType`.

"""
typeofSrc()::SrcType = srcType  # srcType !== SrcType, but

"""
    typeofSrc(x::Int)::SrcType

set and return the value of the global `GeoEfficiency.srcType` corresponding to `x`.

  *  srcUnknown = -1 also any negative integer treated as so,
  *  srcPoint   = 0,
  *  srcLine    = 1,
  *  srcDisk    = 2,
  *  srcVolume  = 3,
  *  srcNotPoint = 4 also any greater than 4 integer treated as so.

"""
function typeofSrc(x::Int)::SrcType
	global srcType = if x < 0
					SrcType(-1)
				elseif x > 4
					SrcType(4)
				else
					SrcType(x)
				end
end #function

#------------------------ setSrcToPoint ---------------------------

"""
    setSrcToPoint()::Bool

return whether the source type is a point or not.
"""
setSrcToPoint()::Bool = srcType === srcPoint

"""

    setSrcToPoint(yes::Bool)::Bool

return whether the source type is a point or not after setting `srcType` to ``srcPoint`` if
`yes` = ``true`` else if `yes` = ``false`` setting it to ``srcNotPoint`` if it was not already
set to other non-point type (``srcDisk``, ``srcLine``, ``srcVolume``).

!!! note
    *  The user can use this function to change the source type any time.
    *  The source type is set the fist time asked for source.

**see also:** [`typeofSrc(::Int)`](@ref).

"""
function setSrcToPoint(yes::Bool)::Bool
	global srcType = if yes
						srcPoint
					elseif srcType in [srcUnknown, srcPoint]
						srcNotPoint
					else
						srcType
					end
	return srcType === srcPoint
end

"""
	setSrcToPoint(prompt::AbstractString)::Bool

return whether the source type is a point or not. only prompt the user to set the source
type if it were not already set before.

**see also:** [`typeofSrc(::Int)`](@ref), [`setSrcToPoint(::Bool)`](@ref).

"""
setSrcToPoint(prompt::AbstractString)::Bool = srcType != srcUnknown ?	setSrcToPoint() :
											setSrcToPoint(input(prompt) |> lowercase != "n" )

#---------------------------- input ---------------------------------

""" # UnExported

    input(prompt::AbstractString = "? ", incolor::Symbol = :green)

return a string delimited by new line excluding the new line. prompt the user with the massage `prompt` defaults to `? `.
`incolor` specify the prompt text color, default to ``green``.

"""
function input(prompt::AbstractString = "? ", incolor::Symbol = :green)
    printstyled(prompt, color=incolor); chomp(readline())
end # function

#---------------------------- getfloat -----------------------------------

"""# UnExported

	getfloat(prompt::AbstractString = "? ", from::Real = 0.0, to::Real = Inf; value::AbstractString="nothing")::Float64

prompts the user with the massage `prompt` defaults to `? ` to input a numerical expression
evaluate to a numerical value and asserts that the value is in the semi open interval [`from`, `to`[
before returning it as a `Float64`.

!!! note
    *  a blank input (i.e just a return) is considered as being ``0.0``.
    *  input from the `console` can be numerical expression not just a number.
    *  All ``5/2``, ``5//2``, ``exp(2)``, ``pi``, ``1E-2``, ``5.2/3``, ``sin(1)``, ``pi/2/3``
       are valid mathematical expressions.
    *  the key word  argument `value` , if provided the function will not ask for input from the
       `console` and take it as if it where inputted from the  `console` [``for test propose mainly``].

# Examples
```
julia> getfloat("input a number:", value="3")
3.0

julia> getfloat("input a number:", value="")
0.0

julia> getfloat("input a number:", value="5/2")
2.5

julia> getfloat("input a number:", value="5//2")
2.5

julia> getfloat("input a number:", value="pi")
3.141592653589793
```

"""
function getfloat(prompt::AbstractString = "? ", from::Real = 0.0, to::Real = Inf; value::AbstractString="nothing") ::Float64
	try
		"nothing" == value ? value = input(prompt) : nothing
		"" == value && return 0.0		# just pressing return is interapted as <0.0>
		val::Float64 =  Meta.parse(value) |> eval |> float
		@assert from <= val < to
		return val

    catch err
        if isa(err, AssertionError)
			@warn("""input `$value` evaluated to be outside the semi open interval [$from, $to[,
			\n Please: provide an adequate value""")
        else
			@warn("""input `$value` cannot be parsed to a valid numerical value!,
			\n Please: provide a valid expression""")
        end #if

        return getfloat(prompt, from, to)

    end #try
end

#---------------------------- detector_info_from_csvFile ------------------------------

"""# UnExported

	 detector_info_from_csvFile(detectors::AbstractString = Detectors,
                                      datadir::AbstractString = dataDir)
return a vector{Detector} based on information in the file of name `detectors` found in the
directory `datadir`.

!!! note
    *  if no path is given the second argument `datadir` is default to ``$(dataDir)`` as set by
       the constant ``dataDir``.
    *  if no file name is specified the name of the predefined file ``$Detectors`` as set by
       the constant ``Detectors``.
    *  the no argument method is the most useful; other methods are mainly for ``test propose``.

"""
function detector_info_from_csvFile(detectors::AbstractString = Detectors,
                                      datadir::AbstractString = dataDir)
    detector_info_array::Matrix{Float64} = Matrix{Float64}(undef, 0, 0)
    @info("opening '$(detectors)'......")
    try
        detector_info_array = readdlm(joinpath(datadir, detectors), ',', header=true)[1];
        return getDetectors(detector_info_array)

    catch err
        if isa(err, SystemError)
		    @error("detector_info_from_csvFile: Some thing went wrong, may be the file '$(joinpath( datadir, detectors))' can't be found")
		end
        rethrow()

    end #try

end #function

#---------------------------- read_from_csvFile --------------------------------

"""# UnExported

	read_from_csvFile(csv_data::AbstractString,
                       datadir::AbstractString = dataDir)::Vector{Float64}

return Vector{Float64} based on data in csv file named `csv_data`. directory `datadir` point to
where the file is located default to ``$(dataDir)`` as set by the constant `dataDir`.

"""
function read_from_csvFile(csv_data::AbstractString, datadir::AbstractString = datadir)::Vector{Float64}
	@info("Opening `$(csv_data)`......")
	try
		indata = readdlm(joinpath(datadir, csv_data), ',',  header=true)[1][:,1]
		return float(indata ) |> sort;

	catch err
	    if isa(err, SystemError)
		    @error("Some thing went wrong, may be `$(csv_data)` can't be found in `$(datadir)`")

		else
		    #println(err)

		end
		return Float64[0.0]

	end #try
end #function

#--------------------------- read_batch_info ------------------------------------

"""# UnExported

	read_batch_info()

read `detectors` and `sources` parameters from the predefined csv files.

Return a tuple
	   (detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		GeoEfficiency_isPoint)

"""
read_batch_info() = read_batch_info(datadir,
                                  detectors,
								 srcHeights,
								    srcRhos,
								   srcRadii,
								 srcLengths)

"""# UnExported

	read_batch_info(datadir::AbstractString,
                  detectors::AbstractString,
                 srcHeights::AbstractString,
                    srcRhos::AbstractString,
                   srcRadii::AbstractString,
                 srcLengths::AbstractString)

read `detectors` and `sources` parameters from the location given in the argument list.

Return a tuple

	   (detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		isPoint)

"""
function read_batch_info(datadir::AbstractString,
                       detectors::AbstractString,
					  srcHeights::AbstractString,
					     srcRhos::AbstractString,
					    srcRadii::AbstractString,
					  srcLengths::AbstractString)

	@info("Starting the Batch Mode ....")
	isPoint = setSrcToPoint("\n Is it a point source {Y|n} ?")

	@info("Reading data from `CSV files` at $datadir .....")
	detectors_array ::Vector{RadiationDetector} = try  detector_info_from_csvFile(detectors, datadir); catch err; getDetectors(); end
	srcHeights_array::Vector{Float64} = read_from_csvFile(srcHeights, datadir)
	srcRhos_array   ::Vector{Float64} = [0.0]
	srcRadii_array  ::Vector{Float64} = [0.0]
	srcLengths_array::Vector{Float64} = [0.0]

	function batchfailure(err::AbstractString)
		@warn(err, ", transfer to direct data input via the `console`......")
		sleep(3); src = source()
		srcHeights_array, srcRhos_array, srcRadii_array  , srcLengths_array   =
		[src[1].Height] , [src[1].Rho] , [src[2]]        , [src[3]]
		return nothing
	end #function

	if srcHeights_array == [0.0]
		batchfailure("`$(srcHeights)` is not found in `$(datadir)`)")

	elseif isPoint
		srcRhos_array =	read_from_csvFile(srcRhos, datadir)

	else
		srcRadii_array = read_from_csvFile(srcRadii, datadir)
		if srcRadii_array == [0.0]
			batchfailure("`$(srcRadii)` is not found in `$(datadir)`)")

		else
			srcLengths_array = read_from_csvFile(srcLengths, datadir)

		end #if
	end #if
	#println("\n Results log\n=============")
	return (
		detectors_array,
		srcHeights_array,
		srcRhos_array,
		srcRadii_array,
		srcLengths_array,
		isPoint,
		)
end #function

#------------------------- getDetectors -------------------------------------

"""

    getDetectors(detectors_array::Vector{<:Detector} = Detector[])

return the  `detectors_array` extended by the entered detectors and sorted according to the
detector volume.
prompt the user to input detector parameters from the `console`.

!!! note
    If no array received in the input an empty array will be created to receive the converted detectors.

"""
function getDetectors(detectors_array::Vector{<:RadiationDetector} = RadiationDetector[])
	Vector{RadiationDetector}(detectors_array); @info("Please, input the detector information via the console")
	while(true)
		try push!(detectors_array, RadiationDetector()); catch err	println(err); @warn("Please: Enter a New Detector"); continue; end
		lowercase(input(
			"""\n
    	                - To add a new detector press return\n
    	                - To quit press 'q'|'Q' then return\n
			\n\t your Choice: """, :blue))  == "q" && break
	end #while
	return detectors_array |> sort
end #function

"""

	getDetectors(detector_info_array::Matrix{<:Real}, detectors_array::Vector{<:Detector} = Detector[] ;console_FB=true)

return `detectors_array`, after extending it with the successfully converted detectors. while,
attempt to convert detectors from the information in `detector_info_array`.

!!! note
    if `console_FB` argument is set to true , the function will call `getDetectors()` to take input
    from the `console` if the `detector_info_array` is empty or contain no numerical element.

"""
function getDetectors(detector_info_array::Matrix{<:Real}, detectors_array::Vector{<:Detector} = Detector[] ; console_FB=true)

	if isempty(detector_info_array)
		if console_FB
			@info("The new detectors information may entred via the console")
			return getDetectors(detectors_array)
		else
		 	error("getDetectors: Empty `detector_info_array`")
		end

	else
		Vector{RadiationDetector}(detectors_array)
		for i_th_line = 1:size(detector_info_array)[1]
			try push!(detectors_array, RadiationDetector((detector_info_array[i_th_line,:])...)) end #try
		end #for

		return detectors_array |> sort
	end   #if
end #function


text:
function correctorstep!(A::Vector{Float64},Q::Vector{Float64},Aprev::Vector{Float64},
    Qprev::Vector{Float64},Fbarforward::Vector{Float64},Fbarbackward::Vector{Float64},
    Qforward::Vector{Float64},Aforward::Vector{Float64},Qbackward::Vector{Float64},Abackward::Vector{Float64},
    colsint::Vector{Int8},acolscor::Vector{Int8},qcolscor::Vector{Int8},
    h::Float64,k::Float64,diff::Float64,μ::Float64,rho::Float64)
    A[:] .= Aprev[colsint] .- h./k.*(Fbarforward[acolscor] .- Fbarbackward[acolscor]);
        # system.branches.A[i][system.solverparams.colsint,n+2] .= (
        #     system.branches.A[i][system.solverparams.colsint,n+1] .-
        #     system.solverparams.h./system.branches.k[i].*
        #     (system.branches.Fbarforward[i][system.solverparams.acolscor] .-
        #     system.branches.Fbarbackward[i][system.solverparams.acolscor]));
    Q[:] .= Qprev[colsint] .- h./k.*(Fbarforward[qcolscor] .- Fbarbackward[qcolscor]) .-
        h.*diff.*π.*μ./rho.*(0.5.*(Qforward./Aforward .+ Qbackward./Abackward));
        # system.branches.Q[i][system.solverparams.colsint,n+2] .= (
        #     system.branches.Q[i][system.solverparams.colsint,n+1] .-
        #     system.solverparams.h./system.branches.k[i].*
        #     (system.branches.Fbarforward[i][system.solverparams.qcolscor] .-
        #     system.branches.Fbarbackward[i][system.solverparams.qcolscor]) .-
        #     system.solverparams.h.*system.solverparams.diffusioncoeff.*π.*
        #     system.solverparams.mu./system.solverparams.rho.*(0.5.*(
        #     system.branches.Qforward[i]./system.branches.Aforward[i] .+
        #     system.branches.Qbackward[i]./system.branches.Abackward[i])));
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling