text:
export MobileNet

import Knet
using Knet.Layers21: Conv, BatchNorm, Linear, Block, Op, Add, Mul
using Knet.Ops21: relu, hardswish, hardsigmoid, pool, mean, reshape2d
using Artifacts

"""
--DRAFT--

    MobileNet(; kwargs...)
    MobileNet(name::String; pretrained=true)

Return a MobileNet model. The first call above returns a randomly initialized model, the
second loads a pretrained model. The models satisfy the indexing interface, e.g. model[1] is
the preprocessing layer and model[1:end-1] omits the imagenet classification top.

Pretrained models:

    name                           top1   ref1   size  time  settings
    ----                           -----  -----  ----  ----  --------
    mobilenet_v1_100_224_tf        70.59  70.6   17M
    mobilenet_v2_100_224_tf        71.88  71.8   14M
    mobilenet_v2_100_224_pt        70.76  71.88  14M
    mobilenet_v3_small_100_224_pt  66.78  67.67  10M
    mobilenet_v3_large_100_224_pt  73.48  74.04  22M

Keyword arguments:

* `width = 1`
* `resolution = 224`
* `input = 16`
* `output = (960,1280)`
* `classes = 1000`
* `activation = hardswish`
* `tfpadding = false`
* `bnupdate = 0.01`
* `bnepsilon = 0.001`
* `dropout = 0.2`
* `block = mobilenet_v3_block`
* `layout = mobilenet_v3_large_layout`
* `normalize = mobilenet_normalize_pt`

References:
* https://arxiv.org/abs/1704.04861
* https://arxiv.org/abs/1801.04381
* https://arxiv.org/abs/1905.02244
* https://pytorch.org/vision/master/models.html
* https://keras.io/api/applications

"""
function MobileNet(
    ;
    width = 1,
    resolution = 224,
    input = 16,
    output = (960,1280),
    classes = 1000,
    activation = hardswish,      # TODO: document kwargs
    tfpadding = false,           # torch:1 keras:((0,1),(0,1)); torch:2, keras:((1,2),(1,2))
    bnupdate = 0.01,             # torchV2:0.1, torchV3:0.01, kerasV1:0.01, kerasV2:0.001
    bnepsilon = 0.001,           # torchV2:1e-5, torchV3:0.001, keras:0.001
    dropout = 0.2,               # kerasV1:0.001, others:0.2
    block = mobilenet_v3_block,
    layout = mobilenet_v3_large_layout,
    normalize = mobilenet_normalize_pt,
)
    global _mobilenet_config = (; tfpadding, bnupdate, bnepsilon)
    α(x) = round(Int, width*x)
    s = Block()
    push!(s, Op(imagenet_preprocess; normalize, resolution))
    push!(s, mobilenet_conv_bn(3, 3, 3, α(input); stride=2, activation))
    channels = input
    for l in layout
        for r in 1:l.repeat
            push!(s, block(channels, l.output; layout=l, repeat=r, width))
            channels = l.output
        end
    end
    top = Block()
    if isempty(output)
        push!(top, Op(pool; op=mean, window=typemax(Int)))
        push!(top, reshape2d)
        push!(top, Linear(α(channels), classes; binit=zeros, dropout))
    else
        push!(top, mobilenet_conv_bn(1, 1, α(channels), α(output[1]); activation))
        push!(top, Op(pool; op=mean, window=typemax(Int)))
        push!(top, reshape2d)
        for o in 2:length(output)
            push!(top, Linear(α(output[o-1]), α(output[o]); binit=zeros, activation))
        end
        push!(top, Linear(α(output[end]), classes; binit=zeros, dropout))
    end
    push!(s, top)
    return s
end

mobilenet_normalize_pt(x) = (x .- [0.485, 0.456, 0.406]) ./ [0.229, 0.224, 0.225]
mobilenet_normalize_tf(x) = ((a,b)=extrema(x); x .* (2/(b-a)) .- ((b+a)/(b-a)))

function mobilenet_conv_bn(w,h,x,y; groups=1, stride=1, activation=nothing)
    c = _mobilenet_config
    p = (w-1)÷2
    padding = (w > 1 && stride > 1 && c.tfpadding ? ((p-1,p),(p-1,p)) : p)
    normalization=BatchNorm(; update=c.bnupdate, epsilon=c.bnepsilon)
    Conv(w,h,x,y; normalization, groups, stride, padding, activation)
end

function relu6(x)
    relu(x; max_value=6)
end

function mobilenet_v1_block(x, y; layout, repeat, width)
    stride = (repeat == 1 ? layout.stride : 1)
    activation = relu6
    x, y = round(Int,width*x), round(Int,width*y)
    Block(
        mobilenet_conv_bn(3, 3, 1, x; groups=x, stride, activation),
        mobilenet_conv_bn(1, 1, x, y; activation),
    )
end

const mobilenet_v1_layout = (
    (repeat=1, output=64, stride=1),
    (repeat=2, output=128, stride=2),
    (repeat=2, output=256, stride=2),
    (repeat=6, output=512, stride=2),
    (repeat=2, output=1024, stride=2),
)

function mobilenet_v2_block(x, y; layout, repeat, width)
    stride = (repeat == 1 ? layout.stride : 1)
    activation = relu6
    x, y = round(Int,width*x), round(Int,width*y)
    b = layout.expansion * x
    s = Block()
    b != x && push!(s, mobilenet_conv_bn(1, 1, x, b; activation))
    push!(s, mobilenet_conv_bn(3, 3, 1, b; groups=b, stride, activation))
    push!(s, mobilenet_conv_bn(1, 1, b, y))
    x == y ? Add(s, identity) : s
end

const mobilenet_v2_layout = (
    (repeat=1, output=16, stride=1, expansion=1),
    (repeat=2, output=24, stride=2, expansion=6),
    (repeat=3, output=32, stride=2, expansion=6),
    (repeat=4, output=64, stride=2, expansion=6),
    (repeat=3, output=96, stride=1, expansion=6),
    (repeat=3, output=160, stride=2, expansion=6),
    (repeat=1, output=320, stride=1, expansion=6),
)

function mobilenet_v3_block(input, output; layout, repeat, width)
    stride = (repeat == 1 ? layout.stride : 1)
    kernel, activation = layout.kernel, layout.activation
    input, output, expand, squeeze = (x->round(Int,width*x)).((input, output, layout.expand, layout.squeeze))
    s = Block()
    channels = input
    if expand > 0
        push!(s, mobilenet_conv_bn(1, 1, channels, expand; activation))
        channels = expand
    end
    push!(s, mobilenet_conv_bn(kernel, kernel, 1, channels; activation, groups=channels, stride))
    if squeeze > 0
        push!(s, Mul(Block(
            Op(pool; op=mean, window=typemax(Int)),
            Conv(1, 1, channels, squeeze; binit=zeros, activation=relu),
            Conv(1, 1, squeeze, channels; binit=zeros, activation=hardsigmoid),
        ), identity))
    end
    push!(s, mobilenet_conv_bn(1, 1, channels, output))
    return input == output && stride == 1 ? Add(s, identity) : s
end

mobilenet_v3_large_layout = (
    # first input=16, last output through 160=>960 hardswish before pool, then linear 960=>1280=>1000
    (repeat=1, input=16, expand=0, squeeze=0, output=16, kernel=3, stride=1, activation=relu),
    (repeat=1, input=16, expand=64, squeeze=0, output=24, kernel=3, stride=2, activation=relu),
    (repeat=1, input=24, expand=72, squeeze=0, output=24, kernel=3, stride=1, activation=relu),
    (repeat=1, input=24, expand=72, squeeze=24, output=40, kernel=5, stride=2, activation=relu),
    (repeat=2, input=40, expand=120, squeeze=32, output=40, kernel=5, stride=1, activation=relu),
    (repeat=1, input=40, expand=240, squeeze=0, output=80, kernel=3, stride=2, activation=hardswish),
    (repeat=1, input=80, expand=200, squeeze=0, output=80, kernel=3, stride=1, activation=hardswish),
    (repeat=2, input=80, expand=184, squeeze=0, output=80, kernel=3, stride=1, activation=hardswish),
    (repeat=1, input=80, expand=480, squeeze=120, output=112, kernel=3, stride=1, activation=hardswish),
    (repeat=1, input=112, expand=672, squeeze=168, output=112, kernel=3, stride=1, activation=hardswish),
    (repeat=1, input=112, expand=672, squeeze=168, output=160, kernel=5, stride=2, activation=hardswish),
    (repeat=2, input=160, expand=960, squeeze=240, output=160, kernel=5, stride=1, activation=hardswish),
)

mobilenet_v3_small_layout = (
    # first input=16, last output through 96=>576 hardswish before pool, then linear 576=>1024=>1000
    (repeat=1, input=16, expand=0, squeeze=8, output=16, kernel=3, stride=2, activation=relu),
    (repeat=1, input=16, expand=72, squeeze=0, output=24, kernel=3, stride=2, activation=relu),
    (repeat=1, input=24, expand=88, squeeze=0, output=24, kernel=3, stride=1, activation=relu),
    (repeat=1, input=24, expand=96, squeeze=24, output=40, kernel=5, stride=2, activation=hardswish),
    (repeat=2, input=40, expand=240, squeeze=64, output=40, kernel=5, stride=1, activation=hardswish),
    (repeat=1, input=40, expand=120, squeeze=32, output=48, kernel=5, stride=1, activation=hardswish),
    (repeat=1, input=48, expand=144, squeeze=40, output=48, kernel=5, stride=1, activation=hardswish),
    (repeat=1, input=48, expand=288, squeeze=72, output=96, kernel=5, stride=2, activation=hardswish),
    (repeat=2, input=96, expand=576, squeeze=144, output=96, kernel=5, stride=1, activation=hardswish),
)

## Pretrained models

mobilenet_models = Dict{String,NamedTuple}(
    "mobilenet_v1_100_224_tf" => (width=1, resolution=224, input=32, output=(), classes=1000, activation=relu6, tfpadding=true, bnupdate=0.01, bnepsilon=0.001, dropout=0.001, block=mobilenet_v1_block, layout=mobilenet_v1_layout, normalize=mobilenet_normalize_tf),
    "mobilenet_v2_100_224_tf" => (width=1, resolution=224, input=32, output=1280, classes=1000, activation=relu6, tfpadding=true, bnupdate=0.001, bnepsilon=0.001, dropout=0.2, block=mobilenet_v2_block, layout=mobilenet_v2_layout, normalize=mobilenet_normalize_tf),
    "mobilenet_v2_100_224_pt" => (width=1, resolution=224, input=32, output=1280, classes=1000, activation=relu6, tfpadding=false, bnupdate=0.1, bnepsilon=1e-5, dropout=0.2, block=mobilenet_v2_block, layout=mobilenet_v2_layout, normalize=mobilenet_normalize_pt),
    "mobilenet_v3_large_100_224_pt" => (width=1, resolution=224, input=16, output=(960,1280), classes=1000, activation=hardswish, tfpadding=false, bnupdate=0.01, bnepsilon=0.001, dropout=0.2, block=mobilenet_v3_block, layout=mobilenet_v3_large_layout, normalize=mobilenet_normalize_pt),
    "mobilenet_v3_small_100_224_pt" => (width=1, resolution=224, input=16, output=(576,1024), classes=1000, activation=hardswish, tfpadding=false, bnupdate=0.01, bnepsilon=0.001, dropout=0.2, block=mobilenet_v3_block, layout=mobilenet_v3_small_layout, normalize=mobilenet_normalize_pt),
)

function MobileNet(s::String; pretrained=true)
    @assert haskey(mobilenet_models, s)  "Please choose from known MobileNet models:\n$(collect(keys(mobilenet_models)))"
    kwargs = mobilenet_models[s]
    model = MobileNet(; kwargs...)
    res = kwargs.resolution
    model(Knet.atype(zeros(Float32,res,res,3,1)))
    pretrained && setweights!(model, joinpath(@artifact_str(s), "$s.jld2"))
    return model
end


text:
mutable struct Parameters
    n1fctrs::Tuple{Real,Real}
    n2ton1fctrs::Tuple{Real,Real}
end

struct Problem
    grids
    n1fctrs
    n1values
    n2mincontinueabs
    n2ton1fctrs
    maxmultipleonestage
    nmax
    curtail_stage_one_buffer
    x1min
    partial
    unimodal
    type
    objective
    power
    toer

    function Problem(
        objective,
        toer,
        power;
        type                  = :TwoStage,
        unimodal              = false,
        partial::Tuple{TI,TI} = (0, 0),
        maxmultipleonestage   = 2.,
        α                     = toer.α,
        β                     = power.β,
        pnull                 = mean(toer.score.prior),
        palt                  = mean(power.score.prior),
        nmax                  = min(
                150,
                maxmultipleonestage*one_stage_sample_size(pnull, α, palt, β)
            ) |> ceil |> Int,
        n1fctrs::Tuple{Real,Real} = (.2/maxmultipleonestage, .51),
        n1min                 = max(partial[2], n1fctrs[1] * nmax) |> ceil |> Int,
        n1max                 = n1fctrs[2] * nmax |> ceil |> Int,
        n1values              = collect(n1min:n1max),
        n2mincontinueabs      = 5,
        n2ton1fctrs::Tuple{Real,Real} = (1.1, 4.0),
        curtail_stage_one_buffer = 2
    ) where {TI<:Integer}

        @assert 0 <= partial[1] <= partial[2]

        x1min = partial[1]

        # prebuild sparse grid space

        function n2(n1, x1)

            if type == :OneStage; return [0] end
            n2min = max(n2mincontinueabs, n1*(n2ton1fctrs[1] - 1)) |> ceil |> Int
            n2max = min(nmax - n1, n1*(n2ton1fctrs[2] - 1)) |> floor |> Int
            # curtail if error rates are too low in stage one
            if (1 - cdf(x1 - curtail_stage_one_buffer, n1, toer.score.prior; partial = partial) <= α) |
                (cdf(x1 + curtail_stage_one_buffer, n1, power.score.prior; partial = partial) <= β)
                return [0]
            end
            # add 0 for early stopping
            return n2min > 0 ? vcat( [0], collect(n2min:n2max) ) : collect(n2min:n2max)
        end

        function c2(n1, x1, n2)

            if n2 == 0; return [-Inf, Inf] end
            candidates = collect(0:(n2 - 1))
            # filter Pr[X1 == x1] * toer[x1] > alpha
            prior      = toer.score.prior
            candidates = candidates[ pmf(x1, n1, prior; partial = partial) .* (1 .- cdf.(candidates, n2, update(prior, x1, n1))) .<= α ]
            # filter Pr[X1 == x1] * tter[x1]) > beta
            prior      = power.score.prior
            candidates = candidates[ pmf(x1, n1, prior; partial = partial) .* cdf.(candidates, n2, update(prior, x1, n1)) .<= β ]
            # check conditional power and type one error rate
            valid = trues(length(candidates))
            for i in 1:length(candidates)
                for cnstr in (power, toer)
                    min, max = cnstr.conditional
                    valid[i] = !(min <= (1 - cdf(candidates[i], n2, update(cnstr.score.prior, x1, n1))) <= max) ? false : valid[i]
                end
            end
            return vcat([-Inf], candidates[valid], [Inf]) # required to be safe!
        end

        function n1_grid(n1)
            [ (x1, n2, c2)
                for x1 in x1min:n1
                for n2 in n2(n1, x1)
                for c2 in c2(n1, x1, n2)
            ]
        end

        grids = Dict(n1 => n1_grid(n1) for n1 in n1values)

        return new(
            grids,
            n1fctrs,
            sort(n1values),
            n2mincontinueabs,
            n2ton1fctrs,
            maxmultipleonestage,
            nmax,
            curtail_stage_one_buffer,
            x1min,
            partial,
            unimodal,
            type,
            objective,
            power,
            toer
        )
    end

end

Base.iterate(problem::Problem, state = 0) = state > 0 ? nothing : (problem, state + 1)
Base.length(problem::Problem) = 1

Base.show(io::IO, problem::Problem) = print(io, string(problem))
Base.show(io::IO, ::MIME"application/prs.juno.inline", problem::Problem) = print(io, string(problem))

Base.string(problem::Problem) = @sprintf "Problem<x1>=%i,%i<=n1<=%i,n<=%i>" problem.x1min problem.n1values[1] problem.n1values[end] problem.nmax

# query precomputed grid values
grid(problem::Problem, n1)  = problem.grids[n1]
grid(problem::Problem)      = [[ (n1, x1, n2, c2) for (x1, n2, c2) in grid] for (n1, grid) in problem.grids] |> x -> vcat(x...)

Base.size(problem::Problem, n1) = size(problem.grids[n1],1)
Base.size(problem::Problem) = [size(problem, n1) for n1 = problem.n1values] |> sum

function build_model(problem::Problem; verbose::Bool = true)

    n1values = problem.n1values
    if verbose
        prog = ProgressMeter.Progress(
            1 + length(n1values) + 3,
            desc      = "Building ILP Problem: ",
            dt        = 0.1,
            barglyphs = ProgressMeter.BarGlyphs("[=> ]"),
            barlen    = 20,
            color     = :gray
        )
    end
    update_progress(str) = if verbose; ProgressMeter.next!(prog; showvalues = [(Symbol(str), "...")], valuecolor = :gray) end

    update_progress("initialising JuMP variables")
    m = Model()
    # ind[x1, n1, n2, c2] == 1 iff n1 = n1, n2(x1) = n2, c2(x1) = c2
    @variable(m, ind[(n1, x1, n2, c2) in grid(problem)], Bin)
    # need to make sure that exactly one n1 value is selected, IP or trick via
    # auxiliary variables: n1_selected[n1] == 1 iff n1 = n1
    @variable(m, n1_selected[n1 in n1values], Bin)
    @constraint(m, sum(n1_selected[n1] for n1 in n1values) == 1)
    if (problem.type == :TwoStage) & problem.unimodal
        @variable(m, is_mode[n1 in n1values, x1 in collect(problem.x1min:n1)], Bin)
        for n1_ in n1values
            # at least one mode (can be on boundary as well!)
            @constraint(m,
              sum( is_mode[n1_, x1] for x1 in collect(problem.x1min:n1_) ) >= 1
            )
        end
    end
    # implement all constraints conditional on n1
    for n1_ in n1values
        update_progress(@sprintf "building constraints for n1=%i" n1_)
        # make sure that n1_selected[n1] is 1 iff any of the other values is assigned
        @constraint(m,
            sum( ind[(n1_, x1, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) )
            <= 5*problem.nmax*n1_selected[n1_]
        )
        for x1_ in collect(problem.x1min:n1_)
            # make it a function in x1
            @constraint(m,
                sum( ind[(n1_, x1_, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) if x1 == x1_ )
                == n1_selected[n1_]
            )
            # make sure that n1(x1) is constant
            if problem.x1min < x1_
                @constraint(m,
                    sum( ind[(n1_, x1_, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) if (x1 == x1_) )
                    == sum( ind[(n1_, x1_ - 1, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) if (x1 == x1_ - 1) )
                )
                if problem.type == :GroupSequential
                    for n2_ in unique([n2 for (x1, n2, c2) in grid(problem, n1_) if n2 > 0])
                        @constraint(m,
                            sum( ind[(n1_, x1_, n2_, c2)] # same n2
                                for (x1, n2, c2) in grid(problem, n1_)
                                if (x1 == x1_) & (n2 == n2_) & isfinite(c2)
                            ) +
                            ind[(n1_, x1_, 0, -Inf)] # or early efficacy
                            >= sum( ind[(n1_, x1_ - 1, n2_, c2)]
                                for (x1, n2, c2) in grid(problem, n1_)
                                if (x1 == x1_ - 1) & (n2 == n2_) & isfinite(c2)
                            )
                        )
                    end
                end
            end
            # make sure we have contiguous stopping for futility ...
            if x1_ > problem.x1min
                @constraint(m, ind[(n1_, x1_ - 1, 0, Inf)] >= ind[(n1_, x1_, 0, Inf)] )
            end
            # ... and for efficacy
            if x1_ < n1_
                @constraint(m, ind[(n1_, x1_ + 1, 0, -Inf)] >= ind[(n1_, x1_, 0, -Inf)] )
            end
            # optional unimodality
            if (problem.type == :TwoStage) & problem.unimodal
                for xx1_ in (problem.x1min + 1):x1_
                    # for a mode at x1_, n must be increasing before ...
                    @constraint(m,
                           sum( n2 * ind[(n1_, xx1_, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) if (x1 == xx1_) )
                         - sum( n2 * ind[(n1_, xx1_ - 1, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) if (x1 == xx1_ - 1) )
                        >= 10*problem.nmax*(is_mode[n1_, x1_] - 1)
                    )
                end
                for xx1_ in (x1_ + 1):n1_
                    # ... and decreasing after
                    @constraint(m,
                           sum( n2 * ind[(n1_, xx1_, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) if (x1 == xx1_) )
                         - sum( n2 * ind[(n1_, xx1_ - 1, n2, c2)] for (x1, n2, c2) in grid(problem, n1_) if (x1 == xx1_ - 1) )
                        <= -10*problem.nmax*(is_mode[n1_, x1_] - 1)
                    )
                end
            end
        end
    end
    update_progress("adding type one error rate constraint")
    @constraint(m,
        sum( integrand_x1(problem.toer, x1, n1, n2, c2; partial_stage_one = problem.partial) * ind[(n1, x1, n2, c2)]
            for (n1, x1, n2, c2) in grid(problem)
        ) <= problem.toer.α
    )
    update_progress("adding power constraint")
    @constraint(m,
        sum( integrand_x1(problem.power, x1, n1, n2, c2; partial_stage_one = problem.partial) * ind[(n1, x1, n2, c2)]
            for (n1, x1, n2, c2) in grid(problem)
        ) >= 1 - problem.power.β
    )
    update_progress("adding objective")
    add!((m, ind), problem.objective, problem; partial = problem.partial)
    update_progress("ready to start ILP solver ")
    if verbose; ProgressMeter.next!(prog) end
    return m, ind, n1_selected
end

function optimise!(m, verbosity::Integer, timelimit::Integer)
    set_optimizer(
        m,
        optimizer_with_attributes(GLPK.Optimizer,
            "msg_lev" => verbosity,
             "tm_lim" => 1000*timelimit
        )
    )
    optimize!(m)
    # termination_status(model.jump_model)
end

function extract_solution(problem::Problem, ind, n1_selected)

    vals = value.(ind)
    # find n1
    x1obs = problem.partial[1]
    n1_ = problem.n1values[findfirst(value.(n1_selected).data .== 1)]
    c2_res = repeat([Inf], n1_ + 1 - x1obs)
    n2_res = zeros(n1_ + 1 - x1obs)
    for (x1, n2, c2) in grid(problem, n1_)
        if vals[(n1_, x1, n2, c2)] == 1
            c2_res[x1 + 1 - x1obs] = c2
            n2_res[x1 + 1 - x1obs] = n2
        end
    end
    return collect(x1obs:n1_), n2_res, c2_res
end

mutable struct OptimalDesign{TI<:Integer,TR<:Real} <: AbstractDesign
    n2::Vector{TI}
    c2::Vector{TR}
    problem::Problem
    info::Dict{String,Any}
end
function OptimalDesign(design::Design, problem::Problem, info::Dict{String,Any})
    return OptimalDesign{eltype(design.n2),eltype(design.c2)}(design.n2, design.c2, problem, info)
end

function optimise(problem::Problem; verbosity = 3, timelimit = 180)

    info = Dict{String,Any}()
    _, info["model build time [s]"], _, _, _ = @timed begin
        m, ind, n1_selected = build_model(problem; verbose = verbosity > 0)
    end
    _, info["model ILP solution time [s]"], _, _, _ = @timed begin
        optimise!(m, verbosity, timelimit)
    end
    _, info["solution extraction time [s]"], _, _, _ = @timed begin
        xx1, nn2, cc2 = extract_solution(problem, ind, n1_selected)
    end
    info["total time [s]"] = sum([val for (key, val) in info])
    info["number of variables"] = size(problem)
    if problem.partial[2] > 0
        return xx1, nn2, cc2, info
    else
        design  = Design(nn2, cc2)
        pnull   = problem.toer.score.bounds[2]
        α       = problem.toer.α
        mlecomp = mlecompatible(design, pnull, α)
        if !mlecomp["compatible"] & (verbosity > 0)
            @warn @sprintf "design is not compatible with MLE-ordering, incompatibility degree is %i/%i" mlecomp["incompatibility degree"] size(mlecomp["details"], 1)
        end
        info["MLE-compatible"] = mlecomp
        return OptimalDesign(design, problem, info)
    end
end

function adapt(design::OptimalDesign, prior::TP, partial::Tuple{TI,TI};
        α         = design.problem.toer.score(design, partial_stage_one = partial),
        β         = 1 - design.problem.power.score(design, partial_stage_one = partial),
        verbosity = 0,
        timelimit = 300
    ) where {TP<:Prior,TI<:Integer}

        objective   = update(design.problem.objective, prior)
        toer        = design.problem.toer
        toer.α      = α
        power       = update(design.problem.power, prior)
        power.β     = β
        pnull, palt = mean(toer.score.prior), mean(power.score.prior)

        n1fctrs     = design.problem.n1fctrs
        n2mincontinueabs = design.problem.n2mincontinueabs
        n2ton1fctrs = design.problem.n2ton1fctrs
        maxmultipleonestage = design.problem.maxmultipleonestage
        curtail_stage_one_buffer = design.problem.curtail_stage_one_buffer
        type = design.problem.type
        unimodal = design.problem.unimodal

        nmax = max(
            min(
                150,
                maxmultipleonestage * one_stage_sample_size(pnull, α, palt, β)
            ) |> ceil |> Int,
            design.problem.nmax
        )
        n1min = max(partial[2], n1fctrs[1] * nmax) |> floor |> Int
        n1max = n1fctrs[2] * nmax |> ceil |> Int

        adaptation_problem = Problem(
            objective,
            toer,
            power;
            partial               = partial,
            maxmultipleonestage   = maxmultipleonestage,
            α                     = α,
            β                     = β,
            pnull                 = pnull,
            palt                  = palt,
            nmax                  = nmax,
            n1fctrs               = n1fctrs,
            n1min                 = n1min,
            n1max                 = n1max,
            n1values              = collect(n1min:n1max),
            n2mincontinueabs      = n2mincontinueabs,
            n2ton1fctrs           = n2ton1fctrs,
            curtail_stage_one_buffer = curtail_stage_one_buffer,
            type                  = design.problem.type,
            unimodal              = unimodal
        )

        m, ind, n1_selected = build_model(adaptation_problem; verbose = verbosity > 0)
        optimise!(m, verbosity, timelimit)
        xx1, nn2, cc2 = extract_solution(adaptation_problem, ind, n1_selected)
        return xx1, nn2, cc2
end

ext: jl


text:
import Base.+;

const MarketDataSet = Dict{String, AbstractMonteCarloProcess};
export MarketDataSet;

#Link an underlying to a single name model
function →(x::String, y::AbstractMonteCarloProcess)
    out = MarketDataSet(x => y)
    return out
end

#Link an multiname underlying to a multi variate model
function →(x::String, y::NDimensionalMonteCarloProcess)
    sep = findfirst("_", x)
    if (isnothing(sep))
        error("Nvariate process must follow the format: INDEX1_INDEX2")
    end
    len_ = length(y.models)
    sep = split(x, "_")
    @assert len_ == length(sep)
    out = MarketDataSet(x => y)

    return out
end

#Joins two different MarketDataSet s
function +(x::MarketDataSet, y::MarketDataSet)
    out = copy(x)
    y_keys = keys(y)
    for y_key in y_keys
        if haskey(out, y_key) || any(out_el -> any(out_el_sub -> out_el_sub == y_key, split(out_el, "_")), keys(out)) || any(y_key_el -> haskey(out, y_key_el), split(y_key, "_"))
            error("check keys")
        else
            out[y_key] = y[y_key]
        end
    end
    return out
end


text:
####
#### Defines list of theory documents
####

theory_docs = Any[
    "Common" => Any["SurfaceFluxes" => "Theory/Common/SurfaceFluxes.md",],
    "Atmos" => Any[
        "AtmosModel" => "Theory/Atmos/AtmosModel.md",
        "Microphysics" => "Theory/Atmos/Microphysics.md",
        "EDMF equations" => "Theory/Atmos/EDMFEquations.md",
        "Turbulence" => "Theory/Atmos/Model/turbulence.md",
        "Tracers" => "Theory/Atmos/Model/tracers.md",
    ],
]


text:
#---------------------------------------------------------------
function temperatureVI(posterior, σ2, T; maxiter = 1, S = 10)
#---------------------------------------------------------------

    D = round(Int, length(posterior)/2)

    gT = observedunitgalaxyvector(T)

    function unpack(param)

        @assert(length(param) == 2D+2)

        local a     = param[1:D]
        local b     = param[1+D:2*D]
        local t     = param[end-1]
        local scale = exp(param[end])

        return a, b, t, scale

    end

    function logp(a, b, t, scale)

        local x = scale * gT

        logpdf(MvNormal(a*t + b, Diagonal(σ2)), x) + logpdf(posterior,[a;b])

    end

    objective(param) = logp(unpack(param)...)

    optimoptions = Optim.Options(show_trace = false, iterations = maxiter, g_tol=1e-6)

    VI(objective, [5*randn(2D+2) for i=1:100]; maxiter = maxiter, S = S, optimiser = LBFGS(), optimoptions = optimoptions)

end

#---------------------------------------------------------------
function temperatureVI2(posterior, σ2, T; maxiter = 1, S = 10)
#---------------------------------------------------------------

    D = round(Int, length(posterior)/2)

    gT = observedunitgalaxyvector(T)

    function unpack(param)

        @assert(length(param) == 2D+1)

        local a     = param[1:D]
        local b     = param[1+D:2*D]
        local scale = exp(param[end])

        return a, b, scale

    end

    function logp(a, b, scale)

        local x = scale * gT

        logpdf(MvNormal(a*(dot(a,x-b)) + b, Diagonal(σ2)), x) + logpdf(posterior,[a;b])

    end

    objective(param) = logp(unpack(param)...)

    optimoptions = Optim.Options(show_trace = false, iterations = maxiter, g_tol=1e-6)

    VI(objective, [5*randn(2D+1) for i=1:100]; maxiter = maxiter, S = S, optimiser = LBFGS(), optimoptions = optimoptions)

end

#---------------------------------------------------------------
function temperatureVI3(posterior, σ2, T; S = S, maxiter = maxiter)
#---------------------------------------------------------------

    # dimension of problem is D

    D = round(Int, length(posterior)/2)

    # create galaxy vector

    g = observedunitgalaxyvector(T)

    # marginals, conditionals

    linemeansigma = marginaline(posterior, σ2)

    function logp(p)

        local s, t = exp(p[1]), p[2]

        local mu, Sigma = linemeansigma(t)

        # local logprior = logpdf(Normal(0.0, 1000.0), s) + logpdf(Normal(0.0, 1000.0), p[2])

        logpdf(MvNormal(mu, Sigma), s*g)# + logprior

    end

    initoptions  = Optim.Options(iterations = 5, g_tol = 1e-6)

    optimoptions = Optim.Options(show_trace = false, show_every = 1, iterations = maxiter, g_tol = 1e-6)

    initsol      = estimatetemperature4(posterior, σ2, T)[1]

    VI(logp, push!([randn(2) for i=1:5], initsol); S = S, optimiser = LBFGS(), optimoptions = optimoptions, initoptions = initoptions)

end


text:
mutable struct StateAM{F <: AbstractFloat, V <: AbstractVector{<:Real}}
  gammaC::F
  gammaT::F
  etaC::V  # K
  etaT::V  # K
  mu::V  # K. mixture loc.
  sigma::V  # K. mixture scale.
  nu::V  # K. mixture df.
  phi::V  # K. mixture skew.
end

function StateAM(data::Data, prior::PriorAM)
  gammaC = rand(prior.gamma)
  gammaT = rand(prior.gamma)
  etaC = rand(prior.eta)
  etaT = rand(prior.eta)
  mu = rand(prior.mu, prior.K)
  sigma = rand(prior.sigma, prior.K)
  nu = rand(prior.nu, prior.K)
  phi = rand(prior.phi, prior.K)

  return StateAM(gammaC, gammaT, etaC, etaT, mu, sigma, nu, phi)
end

# TODO
function tovec(s::StateAM)
  return [s.mu; log.(s.sigma); log.(s.nu); s.phi;
          logit(s.gammaC); logit(s.gammaT);
          MCMC.simplex_transform(s.etaC);
          MCMC.simplex_transform(s.etaT)]
end

mufromvec(K::Integer, v::AbstractVector{<:Real}) = v[1:K]
sigmafromvec(K::Integer, v::AbstractVector{<:Real}) = v[K+1:2K]
nufromvec(K::Integer, v::AbstractVector{<:Real}) = v[2K+1:3K]
phifromvec(K::Integer, v::AbstractVector{<:Real}) = v[3K+1:4K]
gammaCfromvec(K::Integer, v::AbstractVector{<:Real}) = v[4K+1]
gammaTfromvec(K::Integer, v::AbstractVector{<:Real}) = v[4K+2]
etaCfromvec(K::Integer, v::AbstractVector{<:Real}) = v[4K+3:5K+1]
etaTfromvec(K::Integer, v::AbstractVector{<:Real}) = v[5K+2:6K]

function fromvec!(s::StateAM, v::AbstractVector{<:Real})
  K = length(s.mu)
  s.gammaC = logistic(gammaCfromvec(K, v))
  s.gammaT = logistic(gammaTfromvec(K, v))
  s.etaC = MCMC.simplex_invtransform(etaCfromvec(K, v))
  s.etaT = MCMC.simplex_invtransform(etaTfromvec(K, v))
  s.mu = mufromvec(K, v)
  s.sigma = exp.(sigmafromvec(K, v))
  s.nu = exp.(nufromvec(K, v))
  s.phi = phifromvec(K, v)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
#=
Takes an array of Janelia whiskers and finds the angle and curvature using
their method
woi is array of whiskers
tracked specifies if a frame should be included or not
=#
function get_curv_and_angle(woi,tracked=trues(length(woi)),follicle=(400.0f0,50.0f0);face_axis='x')
    curv=zeros(Float64,length(woi))
    aa=zeros(Float64,length(woi))

    #Get angle and curvature from Janelia
    f1=round(Int32,follicle[1])
    f2=round(Int32,follicle[2])
    for i=1:length(woi)

        if (length(woi[i].x)>3)&(tracked[i])
            mymeas=JT_measure(woi[i],f1,f2,face_axis)
            curv[i]=unsafe_wrap(Array,mymeas.data,8)[4]
            aa[i]=unsafe_wrap(Array,mymeas.data,8)[3]

            if isnan(curv[i])|isnan(aa[i])
                tracked[i]=false
            end
        end
    end

    #Interpolate missing data points
    A_x = find(tracked)
    knots = (A_x,)

    itp_a = interpolate(knots, aa[tracked], Gridded(Linear()))
    itp_c = interpolate(knots, curv[tracked], Gridded(Linear()))

    for i=1:length(woi)

        if (!tracked[i])&((i>A_x[1])&(A_x[end]>i))
            curv[i]=itp_c(i)
            aa[i]=itp_a(i)
        end
    end

    for i=1:length(woi)
        if isnan(curv[i])
            curv[i]=itp_c(i)
        end
        if isnan(aa[i])
            aa[i]=itp_a(i)
        end
    end

    (curv,aa)
end

function culm_dist(x,y,thres)
    tot=0.0
    outind=1
    for i=(length(x)-1):-1:1

        tot += sqrt((x[i]-x[i+1])^2+(y[i]-y[i+1])^2)

        if tot>thres
            outind=i
            break
        end
    end
    outind
end

#=
Get curvature of just particular segment of whisker (almost certainly a faster
not janelia way to do this)
-xx is array of x positions for each whisker
-yy is array of y positions for each whisker
-woi is the janelia array of whiskers
-tracked specifies of the whisker should be included or not.
  this variable will be modified to false if the curvature of the high SNR segment cannot be determined
=#
function get_ip_curv_array(xx,yy,woi,tracked,p1,p2)

    c_ip=zeros(Float64,length(xx))
    for i=1:length(xx)
        if tracked[i]
            c_ip[i]=get_ip_curv(xx[i],yy[i],woi[i],p1,p2)
            if isnan(c_ip[i])
                tracked[i]=false
                c_ip[i]=0.0
            end
        end
    end
    c_ip
end

function get_ip_curv(xx,yy,woi,p1,p2)

    ip_1=WhiskerTracking.culm_dist(xx,yy,p1)
    ip_2=WhiskerTracking.culm_dist(xx,yy,p2)

    new_wx=xx[ip_2:ip_1]
    new_wy=yy[ip_2:ip_1]

    new_mywhiskers=deepcopy(woi)
    new_mywhiskers.len=length(new_wx)
    new_mywhiskers.x=convert(Array{Float32,1},new_wx)
    new_mywhiskers.y=convert(Array{Float32,1},new_wy)
    new_mywhiskers.thick=ones(Float32,length(new_wx))
    new_mywhiskers.scores=ones(Float32,length(new_wx))

    mymeas=WhiskerTracking.JT_measure(new_mywhiskers,400.0f0,50.0f0)
    unsafe_wrap(Array,mymeas.data,8)[4]

end


text:
module Ops21

using AutoGrad: AutoGrad, @primitive
using SpecialFunctions: erf

include("activation.jl"); export elu, gelu, hardsigmoid, hardswish, relu, selu, sigm, swish
include("linear.jl"); export linear
include("batchnorm.jl"); export batchnorm
include("conv.jl"); export conv
include("pool.jl"); export pool, mean
include("zeropad.jl"); export zeropad
include("reshape2d.jl"); export reshape2d

import Knet.Ops20
const dropout = Ops20.dropout; export dropout
const bmm = Ops20.bmm; export bmm

end


text:
#!/usr/bin/env julia

using EzXMLPatched

# Parse an XML string
# (use `readxml(<filename>)` to read a document from a file).
doc = parsexml("""
<primates>
    <genus name="Homo">
        <species name="sapiens">Human</species>
    </genus>
    <genus name="Pan">
        <species name="paniscus">Bonobo</species>
        <species name="troglodytes">Chimpanzee</species>
    </genus>
</primates>
""")

# Get the root element from `doc`.
primates = root(doc)

# Iterate over child elements.
for genus in eachelement(primates)
    # Get an attribute value by name.
    genus_name = genus["name"]
    println("- ", genus_name)
    for species in eachelement(genus)
        # Get the content within an element.
        species_name = nodecontent(species)
        println("  └ ", species["name"], " (", species_name, ")")
    end
end
println()

# Find texts using XPath query.
for species_name in nodecontent.(findall("//species/text()", primates))
    println("- ", species_name)
end


text:
"""
Returns a symbolic name for the color interpretation.

This is derived from the enumerated item name with the `GCI_` prefix removed,
but there are some variations. So `GCI_GrayIndex` returns "Gray" and
`GCI_RedBand` returns "Red". The returned strings are static strings and should
not be modified or freed by the application.
"""
_getcolorinterpretationname(colorinterp::GDALColorInterp) =
    GDALGetColorInterpretationName(colorinterp)::Ptr{UInt8}

colorinterpname(colorinterp::GDALColorInterp) =
    unsafe_string(_getcolorinterpretationname(colorinterp))

"Get color interpretation corresponding to the given symbolic name."
_getcolorinterpretationbyname(pszName::Ptr{UInt8}) =
    GDALGetColorInterpretationByName(pszName)::GDALColorInterp

colorinterp(name::String) = _getcolorinterpretationbyname(pointer(name))

"Color Interpretation value for band"
_getrastercolorinterpretation(band::GDALRasterBandH) =
    GDALGetRasterColorInterpretation(band)::GDALColorInterp

getcolorinterp(band::GDALRasterBandH) = _getrastercolorinterpretation(band)

"Set color interpretation of a band."
_setrastercolorinterpretation(band::GDALRasterBandH,
                              colorinterp::GDALColorInterp) =
    GDALSetRasterColorInterpretation(band, colorinterp)

function setrastercolorinterp!(band::GDALRasterBandH, colorinterp::GDALColorInterp)
    result = _setrastercolorinterpretation(band, colorinterp)
    if result == CE_Failure
        error("$colorinterpname(colorinterp) is unsupported by raster.")
    end
end

"""
Fetch the color table associated with band.

If there is no associated color table, the return result is `NULL`. The
returned color table remains owned by the `GDALRasterBand`, and can't be
depended on for long, nor should it ever be modified by the caller.
"""
_getrastercolortable(band::GDALRasterBandH) =
    GDALGetRasterColorTable(band)::GDALColorTableH

function getrastercolortable(band::GDALRasterBandH)
    result = _getrastercolortable(band)
    (result == C_NULL) && error("Colortable not available")
    result
end

"""
Set the raster color table.

The driver will make a copy of all desired data in the colortable. It remains
owned by the caller after the call.

### Parameters
* `poCT`    the color table to apply. This may be NULL to clear the color table
(where supported).

### Returns
`CE_None` on success, or `CE_Failure` on failure. If the action is unsupported
by the driver, a value of `CE_Failure` is returned, but no error is issued.
"""
_setrastercolortable(band::GDALRasterBandH, poCT::GDALColorTableH) =
    GDALSetRasterColorTable(band, poCT)::CPLErr

function setrastercolortable!(band::GDALRasterBandH, poCT::GDALColorTableH)
    result = _setrastercolortable(band, poCT)
    (result == CE_Failure) && error("failed to set raster colortable")
end

"Construct a new color table."
_createcolortable(eInterp::GDALPaletteInterp = GPI_RGB) =
    GDALCreateColorTable(eInterp)::GDALColorTableH

"Destroys a color table."
_destroycolortable(hTable::GDALColorTableH) =
    GDALDestroyColorTable(hTable)

"Make a copy of a color table."
_clonecolortable(hTable::GDALColorTableH) =
    GDALCloneColorTable(hTable)::GDALColorTableH

"""
Fetch palette interpretation.

The returned value is used to interpret the values in the `GDALColorEntry`.

### Returns
palette interpretation enumeration value, usually `GPI_RGB`.
"""
_getpaletteinterpretation(hTable::GDALColorTableH) =
    GDALGetPaletteInterpretation(hTable)::GDALPaletteInterp

"Get number of color entries in table."
_getcolorentrycount(hTable::GDALColorTableH) =
    GDALGetColorEntryCount(hTable)::Cint

"""
Fetch a color entry from table.

### Parameters
* `i`   entry offset from `0` to `GetColorEntryCount()-1`.

### Returns
pointer to internal color entry, or NULL if index is out of range.
"""
_getcolorentry(htable::GDALColorTableH, i::Integer) =
    GDALGetColorEntry(hTable, i)::Ptr{GDALColorEntry}

getcolorentry(htable::GDALColorTableH, i::Integer) =
    _getcolorentry(hTable, i-1)

"""
Fetch a table entry in RGB format.

In theory this method should support translation of color palettes in non-RGB
color spaces into RGB on the fly, but currently only works on RGB color tables.

### Parameters
* `i`           entry offset from `0` to `GetColorEntryCount()-1`.
* `poEntry`     existing GDALColorEntry to be overrwritten with the RGB values.

### Returns
`TRUE` on success, or `FALSE` if the conversion isn't supported.
"""
_getcolorentryasrgb(hTable::GDALColorTableH, i::Integer,
                    poEntry::Ptr{GDALColorEntry}) =
    GDALGetColorEntryAsRGB(hTable, i, poEntry)::Cint

function getcolorentryasrgb!(hTable::GDALColorTableH, i::Integer,
                             poEntry::Ptr{GDALColorEntry})
    result = Bool(_getcolorentryasrgb(hTable, i-1, poEntry))
    result || error("conversion of color entry to RGB not supported")
end

"""
Set entry in color table.

Note that the passed in color entry is copied, and no internal reference to it
is maintained. Also, the passed in entry must match the color interpretation of
the table to which it is being assigned.

The table is grown as needed to hold the supplied offset.

### Parameters
* `i`           entry offset from `0` to `GetColorEntryCount()-1`.
* `poEntry`     value to assign to table.
"""
_setcolorentry(hTable::GDALColorTableH, i::Integer,
               poEntry::Ptr{GDALColorEntry}) =
    GDALSetColorEntry(hTable, i, poEntry)

setcolorentry!(hTable::GDALColorTableH, i::Integer,
               poEntry::Ptr{GDALColorEntry}) =
    _setcolorentry(hTable, i-1, poEntry)


text:
# Test the "Actualize" actually does what it's supposed to.
@testset "Testing utils" begin
    backend = AutoTM.Backend("CPU")

    #####
    ##### Actualize, @closure and friends
    #####

    # Inference
    f = AutoTM.@closure AutoTM.Zoo.vgg19_inference(16)
    @test isa(f(), AutoTM.Utils.Actualizer)
    F = AutoTM.actualize(backend, f)

    @test isa(F(), nGraph.TensorView)

    # Training
    f = AutoTM.@closure AutoTM.Zoo.vgg_training(AutoTM.Zoo.Vgg19(), 16)
    @test isa(f(), AutoTM.Utils.Actualizer)
    F = AutoTM.actualize(backend, f)

    @test isa(F(), nGraph.TensorView)

    #####
    ##### findonly
    #####

    x = [2, 3, 4]
    # Examples where just one object if found
    @test AutoTM.Utils.findonly(isequal(2), x) == 1
    @test AutoTM.Utils.findonly(x -> iszero(mod(x, 3)), x) == 2
    # No items found
    @test_throws ArgumentError AutoTM.Utils.findonly(x -> iszero(x), x)
    # More than 1 item found.
    @test_throws ArgumentError AutoTM.Utils.findonly(x -> x > 0, x)

    #####
    ##### dict_push!
    #####

    d = Dict{Int,Vector{Int}}()
    AutoTM.Utils.dict_push!(d, 1, 2)
    @test d == Dict(1 => [2])

    AutoTM.Utils.dict_push!(d, 2, 2)
    @test d == Dict(1 => [2], 2 => [2])

    AutoTM.Utils.dict_push!(d, 1, 10)
    @test d == Dict(1 => [2, 10], 2 => [2])

    #####
    ##### vflatten
    #####

    x = [1,2]
    y = [3,4]
    @test collect(AutoTM.Utils.vflatten(x, y)) == [1,2,3,4]
end

@testset "Testion IOConfig" begin
    DRAM = AutoTM.Utils.DRAM
    PMEM = AutoTM.Utils.PMEM

    io = AutoTM.Utils.IOConfig((DRAM, PMEM), (PMEM, DRAM))

    @test length(io) == 4
    @test collect(io) == [DRAM, PMEM, PMEM, DRAM]
    @test io[1] == DRAM
    @test io[3] == PMEM
    @test_throws BoundsError io[0]
    @test_throws BoundsError io[5]

    # Test `isless` functionality.
    @test io < AutoTM.Utils.IOConfig((DRAM, PMEM), (PMEM, PMEM))
    @test io > AutoTM.Utils.IOConfig((DRAM, DRAM), (PMEM, PMEM))

    # Test `setindex`
    @test Base.setindex(io, 2, DRAM) == AutoTM.Utils.IOConfig((DRAM, DRAM), (PMEM, DRAM))
    @test Base.setindex(io, 4, PMEM) == AutoTM.Utils.IOConfig((DRAM, PMEM), (PMEM, PMEM))
end

@testset "Testing metagraph" begin
    # Create a dummy graph like this:
    #
    # 1 --> 2 --> 3
    #       ^     |
    #       |     |
    #       *-----*
    #
    # Associate metadata like this
    # Vertex:
    #   1 -- 'a'
    #   2 -- 'b'
    #   3 -- 'c'
    #
    # Edge:
    #   1 => 2 -- 1
    #   2 => 3 -- 2
    #   3 => 2 -- 3

    LG = AutoTM.Utils.LightGraphs
    AU = AutoTM.Utils

    mg = AU.MetaGraph(LG.SimpleDiGraph(),Int,Char)
    # Add the vertices with their metadata
    LG.add_vertex!(mg, 'a')
    LG.add_vertex!(mg, 'b')
    LG.add_vertex!(mg, 'c')

    @test AU.getmeta(mg, 1) == 'a'
    @test AU.getmeta(mg, 2) == 'b'
    @test AU.getmeta(mg, 3) == 'c'

    LG.add_edge!(mg, 1, 2, 1)
    LG.add_edge!(mg, 2, 3, 2)
    LG.add_edge!(mg, 3, 2, 3)

    @test AU.getmeta(mg, LG.edgetype(mg)(1, 2)) == 1
    @test AU.getmeta(mg, LG.edgetype(mg)(2, 3)) == 2
    @test AU.getmeta(mg, LG.edgetype(mg)(3, 2)) == 3

    # Test some of the standard lightgraphs functions to make sure things are forwarding
    # properly.
    @test LG.ne(mg) == 3
    @test LG.nv(mg) == 3
    @test LG.outneighbors(mg, 1) == [2]
    @test LG.inneighbors(mg, 2) == [1, 3]

    et = LG.edgetype(mg)
    @test Set(AU.inedges(mg, 2)) == Set((et(1, 2), et(3, 2)))
    @test Set(AU.outedges(mg, 2)) == Set((et(2, 3),))

    # Check some of the find functions.
    v = AU.find_vertex((g,v) -> AU.getmeta(g, v) == 'a', mg)
    @test v == 1
    @test AU.getmeta(mg, v) == 'a'

    # Check when number of matches is greater than 1
    @test_throws ArgumentError AU.find_vertex((g,v) -> true, mg)

    # Check when number of matches is 0
    @test_throws ArgumentError AU.find_vertex((g,v) -> false, mg)

    # Do the same thing for edges.
    e = AU.find_edge((g,v) -> AU.getmeta(g, v) == 3, mg)
    @test e == et(3, 2)
    @test_throws ArgumentError AU.find_edge((g,v) -> true, mg)
    @test_throws ArgumentError AU.find_edge((g,v) -> false, mg)
end

@testset "Testing Allocator" begin
    # Get to the AllocatorModel ... Module ...
    AllocatorModel = AutoTM.Utils.AllocatorModel

    # Instantiate a memory allocator with a memory limit of 100 and alignment of 10
    M = AllocatorModel.MemoryAllocator(100, 10)
    node_list = M.node_list
    @test length(node_list) == 1

    # Allocate a block of size 1. We expect it to be located at the beginning of the
    # block.
    offset = AllocatorModel.allocate(M, 1)
    @test offset == 0
    @test length(node_list) == 2

    # The first block should have size 10 due to alignment
    @test AllocatorModel.isfree(node_list[1]) == false
    @test sizeof(node_list[1]) == 10
    @test AllocatorModel.isfree(node_list[2]) == true
    @test sizeof(node_list[2]) == 90

    # Try freeing this block - make sure we get just a single block back
    AllocatorModel.free(M, offset)

    @test length(node_list) == 1
    @test AllocatorModel.isfree(node_list[1]) == true
    @test sizeof(node_list[1]) == 100

    # Allocate three blocks - free the middle one, then free the last
    a = AllocatorModel.allocate(M, 11)
    b = AllocatorModel.allocate(M, 9)
    c = AllocatorModel.allocate(M, 20)

    # Consistency Checks
    @test a == 0
    @test b == 20
    @test c == 30
    @test length(node_list) == 4
    @test AllocatorModel.isfree(node_list[1]) == false
    @test AllocatorModel.isfree(node_list[2]) == false
    @test AllocatorModel.isfree(node_list[3]) == false
    @test AllocatorModel.isfree(node_list[4]) == true
    @test sizeof(node_list[1]) == 20
    @test sizeof(node_list[2]) == 10
    @test sizeof(node_list[3]) == 20
    @test sizeof(node_list[4]) == 50

    AllocatorModel.free(M, b)
    @test length(node_list) == 4
    @test AllocatorModel.isfree(node_list[1]) == false
    @test AllocatorModel.isfree(node_list[2]) == true
    @test AllocatorModel.isfree(node_list[3]) == false
    @test AllocatorModel.isfree(node_list[4]) == true
    @test sizeof(node_list[1]) == 20
    @test sizeof(node_list[2]) == 10
    @test sizeof(node_list[3]) == 20
    @test sizeof(node_list[4]) == 50

    AllocatorModel.free(M, c)
    @test length(node_list) == 2
    @test AllocatorModel.isfree(node_list[1]) == false
    @test AllocatorModel.isfree(node_list[2]) == true
    @test sizeof(node_list[1]) == 20
    @test sizeof(node_list[2]) == 80

    # Test that an oversized allocation fails
    @test isnothing(AllocatorModel.allocate(M, 200))
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using SymbolicWedderburn
using PermutationGroups
using Cyclotomics

using SparseArrays

using DynamicPolynomials
using SumOfSquares
using SCS

const OPTIMIZER = optimizer_with_attributes(
    SCS.Optimizer,
    "acceleration_lookback" => 10,
    "max_iters" => 10_000,
    "eps" => 1e-6,
    "linear_solver" => SCS.DirectSolver,
)

include(joinpath(@__DIR__, "action_polynomials.jl"))

N = 4

@polyvar x[1:N]

f = sum(x) +
    sum(x .^ 2) +
    (sum((x .+ 1) .^ 2 .* (x .+ x') .^ 2))^2 * (1 + sum(x .^ 2))

no_symmetry = let f = f, basis = monomials(x, 0:(DynamicPolynomials.maxdegree(f)÷2))
    stats = @timed let m = SOSModel(OPTIMIZER)
        @variable m t
        @objective m Max t
        @variable m sos SOSPoly(basis)
        @constraint m f - t == sos
        m
    end
    m = stats.value
    optimize!(m)
    @info (m,) termination_status(m) objective_value(m) creation_time=(stats.time) solve_time(m)
    (
        status=termination_status(m),
        objective=objective_value(m),
        creation_t=stats.time,
        solve_t=solve_time(m)
    )
end

semisimple_dec = let f = f, G = PermGroup([perm"(1,2)", Perm([2:N; 1])])
    basis = monomials(x, 0:(DynamicPolynomials.maxdegree(f)÷2))

    stats = @timed begin
        sa_basis, symmetry_adaptation_time, = @timed let
            ssimple_basis = SymbolicWedderburn.symmetry_adapted_basis(
                Rational{Int},
                G,
                VariablePermutation(),
                basis,
                semisimple=true
            )
            sp = sparse.(ssimple_basis)
            droptol!.(sp, 1e-12)
        end

        m = SOSModel(OPTIMIZER)

        let basis = basis, R = sa_basis
            @variable m t
            @objective m Max t

            soses = @variable m [r in R] SOSPoly(FixedPolynomialBasis(r * basis))
            @constraint m f - t == sum(soses)
        end
        m
    end
    m = stats.value
    optimize!(m)
    @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time=(stats.time) symmetry_adaptation_time
    (
        status=termination_status(m),
        objective=objective_value(m),
        creation_t=stats.time,
        solve_t=solve_time(m)
    )
end

@assert no_symmetry.status == MOI.OPTIMAL
@assert semisimple_dec.status == MOI.OPTIMAL
@assert isapprox(no_symmetry.objective, semisimple_dec.objective, atol=1e-4)
@assert semisimple_dec.solve_t < no_symmetry.solve_t

# This is faster, but not quite the speedup we'd like to see.
# Let us try now to use the decomposition into simple summands
include(joinpath(@__DIR__, "util.jl"))

wedderburn_dec = let f=f, G = PermGroup([perm"(1,2)", Perm([2:N; 1])]), T = Float64

    basis_half = monomials(x, 0:(DynamicPolynomials.maxdegree(f)÷2))

    wedderburn, symmetry_adaptation_time, = @timed WedderburnDecomposition(
        T,
        G,
        VariablePermutation(),
        monomials(x, 0:DynamicPolynomials.maxdegree(f)),
        basis_half
    )

    M = [SymbolicWedderburn.basis(wedderburn)[x*y] for x in basis_half, y in basis_half]

    stats = @timed begin
        m = let m = JuMP.Model(OPTIMIZER)
            JuMP.@variable m t
            JuMP.@objective m Max t
            psds = [
                JuMP.@variable(m, [1:d, 1:d] in PSDCone())
                for d in size.(direct_summands(wedderburn), 1)
            ]

            # preallocating
            Mπs = zeros.(T, size.(psds))
            M_orb = similar(M, T)
            tmps = SymbolicWedderburn._tmps(wedderburn)

            C = DynamicPolynomials.coefficients(f-t, SymbolicWedderburn.basis(wedderburn))

            for iv in invariant_vectors(wedderburn)
                c = dot(C, iv)
                # C needs to be invariant under G,
                # which is eqivalent to being constant on the orbits
                let cc = C[findfirst(!iszero, iv)]
                    @assert all(C[m] == cc for m in findall(!iszero, iv))
                end

                # invariant_constraint! is defined locally in util.jl
                # it essentially amounts to averaging constraint matrices defined by M
                # with weights provided by iv and storing the result in (dense) matrix M_orb
                M_orb = invariant_constraint!(M_orb, M, iv)
                Mπs = SymbolicWedderburn.diagonalize!(Mπs, M_orb, wedderburn, tmps)

                JuMP.@constraint m sum(
                    dot(Mπ, Pπ) for (Mπ, Pπ) in zip(Mπs, psds) if !iszero(Mπ)
                ) == c
            end
            m
        end
    end

    m = stats.value
    optimize!(m)
    @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time=(stats.time) symmetry_adaptation_time
    (
        status=termination_status(m),
        objective=objective_value(m),
        creation_t=stats.time,
        solve_t=solve_time(m)
    )
end

@assert wedderburn_dec.status == MOI.OPTIMAL
@assert isapprox(wedderburn_dec.objective, semisimple_dec.objective, atol=1e-4)
@assert semisimple_dec.solve_t > wedderburn_dec.solve_t
@assert wedderburn_dec.solve_t < 2.0


text:
type DataSet{m}
    X::Matrix{Float64}        # n x m matrix of n input variables by m samples
    Y::Matrix{Float64}        # k x m matrix of k output variables by m samples
    metric::PreMetric         # error measurement method (from Distances package)
    weights::Vector{Float64}  # samples weighting

    function DataSet(X::Matrix{Float64},
        Y::Matrix{Float64},
        metric::PreMetric,
        weights::Vector{Float64})

        @assert size(X,2) == size(Y,2) "incompatible number of samples between input and output"
        @assert size(X,2) == length(weights) "incompatible number of samples and weights"

        new(X, Y, metric, weights)
    end
end

DataSet(X::Matrix{Float64},
        Y::Matrix{Float64},
        metric::PreMetric=SqEuclidean(),
        weights::Vector{Float64}=ones(size(X,2))) =
    DataSet{size(X,2)}(X, Y, metric, weights)

# if vector, assume the unique dimension is along samples
function DataSet(X::Union(Vector{Float64}, Float64),
                 Y::Union(Vector{Float64}, Float64),
                 metric::PreMetric=SqEuclidean(),
                 weights::Vector{Float64}=length(X))

    DataSet{1}(reshape(X, 1, length(X)),
               reshape(Y, 1, length(Y)),
               metric, weights)
end

compatible(nn::NN, ts::DataSet) = (inputsize(nn)==size(ts.X,1)) && (size(nn)[end]==size(ts.Y,1))
nsamples{m}(ts::DataSet{m}) = m

function shuffle!(ts::DataSet; seed::Integer=0)
    (seed != 0) && srand(seed)
    idx = shuffle(collect(1:nsamples(ts)))

    ts.X       = ts.X[:,idx]
    ts.Y       = ts.Y[:,idx]
    ts.weights = ts.weights[idx]
end

function split(ts::DataSet, fraction=0.8)
    pos = round(Int, nsamples(ts) * fraction)

    ( DataSet(ts.X[:,1:pos-1],
              ts.Y[:,1:pos-1],
              ts.metric,
              ts.weights[1:pos-1]),
      DataSet(ts.X[:,pos:end],
              ts.Y[:,pos:end],
              ts.metric,
              ts.weights[pos:end]) )
end

function show{m}(io::IO, ts::DataSet{m})
    print(io, "training data : ")
    print(io, "[$m] samples, ")
    print(io, "input size [$(size(ts.X,1))], ")
    print(io, "output size [$(size(ts.Y,1))], ")
    print(io, "metric [$(typeof(ts.metric))]")
    any(ts.weights .!= 1.) && print(io, ", with weights")
    println(io,"")
end

max_stars_repo_name: fredo-dedup/SimpleNeuralNets

max_issues_repo_name: fredo-dedup/SimpleNeuralNets

max_forks_repo_name: fredo-dedup/SimpleNeuralNets


text:
"""
    T2partSEcorr([io = stderr,] T2distributions::Array{T,4}; <keyword arguments>)
    T2partSEcorr([io = stderr,] T2distributions::Array{T,4}, opts::T2partOptions{T})

Analyzes T2 distributions produced by [`T2mapSEcorr`](@ref) to produce data maps of a series of parameters.

# Arguments
- `T2distributions`: 4D array with data as `(row, column, slice, T2 amplitude)`
- A series of optional keyword argument settings which will be used to construct a [`T2partOptions`](@ref) struct internally, or a [`T2partOptions`](@ref) struct directly

# Ouputs
- `maps`: a dictionary containing the following 3D data maps as fields:
    - `"sfr"`: small pool fraction, e.g. myelin water fraction (`MatrixSize` 3D array)
    - `"sgm"`: small pool geometric mean T2 (`MatrixSize` 3D array)
    - `"mfr"`: medium pool fraction, e.g. intra/extracellular water fraction (`MatrixSize` 3D array)
    - `"mgm"`: medium pool geometric mean T2 (`MatrixSize` 3D array)

# Examples
```julia-repl
julia> dist = DECAES.mock_T2_dist(MatrixSize = (100,100,1), nT2 = 40); # mock distribution with size 100x100x1x40

julia> maps = T2partSEcorr(dist; T2Range = (10e-3, 2.0), SPWin = (10e-3, 25e-3), MPWin = (25e-3, 200e-3), Silent = true); # compute T2-parts maps

julia> maps # MWF is contained in maps["sfr"]
Dict{String, Any} with 4 entries:
  "sgm" => [0.0105638 0.0128673 … 0.0120685 0.0117591; 0.0108218 0.0139536 … 0.0118711 0.011…
  "mfr" => [0.879499 0.889248 … 0.886412 0.897167; 0.927041 0.887472 … 0.872513 0.89026; … ;…
  "sfr" => [0.120501 0.110752 … 0.113588 0.102833; 0.0729592 0.112528 … 0.127487 0.10974; … …
  "mgm" => [0.0621055 0.0617419 … 0.0612094 0.0615634; 0.0607907 0.0625588 … 0.0629157 0.062…
```

See also:
* [`T2mapSEcorr`](@ref)
"""
T2partSEcorr(T2distributions::Array{T,4}; kwargs...) where {T} = T2partSEcorr(stderr, T2distributions; kwargs...)
T2partSEcorr(io::IO, T2distributions::Array{T,4}; kwargs...) where {T} = T2partSEcorr(io, T2distributions, T2partOptions(T2distributions; kwargs...))
T2partSEcorr(T2distributions::Array{T,4}, opts::T2partOptions{T}) where {T} = T2partSEcorr(stderr, T2distributions, opts)

function T2partSEcorr(io::IO, T2distributions::Array{T,4}, opts::T2partOptions{T}) where {T}
    @assert size(T2distributions) == (opts.MatrixSize..., opts.nT2)

    # Print settings to terminal
    !opts.Silent && printbody(io, _show_string(opts))

    # Initial output
    maps = init_output_t2parts(T2distributions, opts)
    thread_buffers = [thread_buffer_maker(opts) for _ in 1:Threads.nthreads()]

    # Run T2-Part analysis
    LinearAlgebra.BLAS.set_num_threads(1) # Prevent BLAS from stealing julia threads
    @inbounds Threads.@threads for I in CartesianIndices(opts.MatrixSize)
        thread_buffer = thread_buffers[Threads.threadid()]
        voxelwise_T2_parts!(thread_buffer, maps, T2distributions, opts, I)
    end
    LinearAlgebra.BLAS.set_num_threads(Threads.nthreads()) # Reset BLAS threads

    return maps
end

function init_output_t2parts(T2distributions::Array{T,4}, opts::T2partOptions{T}) where {T}
    maps = Dict{String,Any}()
    maps["sfr"] = fill(T(NaN), opts.MatrixSize...)
    maps["sgm"] = fill(T(NaN), opts.MatrixSize...)
    maps["mfr"] = fill(T(NaN), opts.MatrixSize...)
    maps["mgm"] = fill(T(NaN), opts.MatrixSize...)
    return maps
end

# =========================================================
# Save thread local results to output maps
# =========================================================
function voxelwise_T2_parts!(thread_buffer, maps, T2distributions, o::T2partOptions{T}, I) where {T}
    @unpack dist, T2_times, sp_range, mp_range, logT2_times_sp, logT2_times_mp, weights = thread_buffer

    # Return nothing if distribution contains NaN entries
    @inbounds for j in 1:o.nT2
        isnan(T2distributions[I, j]) && return nothing
    end

    # Load in voxel T2 distribution
    @inbounds @simd for j in 1:o.nT2
        dist[j] = T2distributions[I, j]
    end

    # Precompute sums and dot products over small pool, medium pool, and entire ranges
    Σ_dist, Σ_dist_sp, Σ_dist_mp = zero(T), zero(T), zero(T)
    dot_sp, dot_mp = zero(T), zero(T)
    @inbounds @simd for j in 1:length(dist)
        Σ_dist += dist[j]
    end
    @inbounds @simd for j in 1:length(sp_range)
        sp = sp_range[j]
        dot_sp    += dist[sp] * logT2_times_sp[j]
        Σ_dist_sp += dist[sp]
    end
    @inbounds @simd for j in 1:length(mp_range)
        mp = mp_range[j]
        dot_mp    += dist[mp] * logT2_times_mp[j]
        Σ_dist_mp += dist[mp]
    end

    # Compute T2 distribution parts
    if Σ_dist > 0
        @inbounds maps["sfr"][I] = o.Sigmoid !== nothing ? dot(dist, weights) / Σ_dist : Σ_dist_sp / Σ_dist
        @inbounds maps["mfr"][I] = Σ_dist_mp / Σ_dist
    end
    if Σ_dist_sp > 0
        @inbounds maps["sgm"][I] = exp(dot_sp / Σ_dist_sp)
    end
    if Σ_dist_mp > 0
        @inbounds maps["mgm"][I] = exp(dot_mp / Σ_dist_mp)
    end

    return nothing
end

# =========================================================
# Utility functions
# =========================================================
function thread_buffer_maker(o::T2partOptions{T}) where {T}
    dist           = zeros(T, o.nT2)
    T2_times       = logrange(o.T2Range..., o.nT2)
    sp_range       = findfirst(>=(o.SPWin[1]), T2_times) : findlast(<=(o.SPWin[2]), T2_times)
    mp_range       = findfirst(>=(o.MPWin[1]), T2_times) : findlast(<=(o.MPWin[2]), T2_times)
    logT2_times    = log.(T2_times)
    logT2_times_sp = logT2_times[sp_range]
    logT2_times_mp = logT2_times[mp_range]
    weights        = sigmoid_weights(o)
    return (; dist, T2_times, sp_range, mp_range, logT2_times_sp, logT2_times_mp, weights)
end

function sigmoid_weights(o::T2partOptions{T}) where {T}
    if o.Sigmoid !== nothing
        # Curve reaches 50% at T2_50perc and is (k and 1-k)*100 percent at T2_50perc +/- T2_kperc
        k, T2_kperc, T2_50perc = T(0.1), o.Sigmoid, o.SPWin[2]
        sigma = abs(T2_kperc / (sqrt(T(2)) * erfinv(2*k-1)))
        (x -> x <= eps(T) ? zero(T) : x).(normccdf.((logrange(o.T2Range..., o.nT2) .- T2_50perc) ./ sigma))
    else
        nothing
    end
end


text:
@diffusion_process FitzHughDiffusionConjug begin
    :dimensions
    process --> 2
    wiener --> 1

    :parameters
    (ϵ, s, γ, β, σ) --> Float64

    :conjugate
    phi(t, x) -->(
        (-x[2],),
        (x[1]-x[1]^3+(1-3*x[1]^2)*x[2],),
        (one(x[1]),),
        (-x[1],),
        (zero(x[1]),),
        (zero(x[1]),),
    )
    nonhypo(x) --> x[2:2]
    num_non_hypo --> 1

    :additional
    constdiff --> true
end

function b(t, x, P::FitzHughNagumoConjug)
    ℝ{2}(
        x[2],
        (
            (P.ϵ - P.γ)*x[1]
            - P.ϵ*(x[1]^3 + (3.0*x[1]^2 - 1.0)*x[2])
            + P.s
            - P.β
            - x[2]
        )
    )
end

σ(t, x, P::FitzHughNagumoConjug) = ℝ{2}(0.0, P.σ)


text:
using ModelingToolkit
using NonlinearSolve
using Test

@variables x y z u
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]

par = [
    σ => 1,
    ρ => 0.1+σ,
    β => ρ*1.1
]
u0 = Pair{Num, Any}[
    x => u,
    y => σ, # default u0 from default p
    z => u-0.1,
]
ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β], name=:ns, defaults=[par; u0])
ns.y = u*1.1
resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns), defaults=ModelingToolkit.defaults(ns))
@test resolved == [1, 0.1+1, (0.1+1)*1.1]

prob = NonlinearProblem(ns, [u=>1.0], Pair[])
@test prob.u0 == [1.0, 1.1, 0.9]
@show sol = solve(prob,NewtonRaphson())

@variables a
@parameters b
top = NonlinearSystem([0 ~ -a + ns.x+b], [a], [b], systems=[ns], name=:top)
top.b = ns.σ*0.5
top.ns.x = u*0.5

res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top), defaults=ModelingToolkit.defaults(top))
@test res == [0.5, 1, 0.1+1, (0.1+1)*1.1]

prob = NonlinearProblem(top, [states(ns, u)=>1.0, a=>1.0], Pair[])
@test prob.u0 == [1.0, 0.5, 1.1, 0.9]
@show sol = solve(prob,NewtonRaphson())

# test NullParameters+defaults
prob = NonlinearProblem(top, [states(ns, u)=>1.0, a=>1.0])
@test prob.u0 == [1.0, 0.5, 1.1, 0.9]
@show sol = solve(prob,NewtonRaphson())


text:
type BarplotGraphics{R<:Real} <: GraphicsArea
  bars::Vector{R}
  color::Symbol
  charWidth::Int
  maxFreq::R
  maxFreqLen::R
  symb::(@compat AbstractString)

  function BarplotGraphics(bars::Vector{R},
                           charWidth::Int,
                           color::Symbol = :blue,
                           symb = "▪")
    charWidth = max(charWidth, 5)
    maxFreq = maximum(bars)
    maxFreqLen = length(string(maxFreq))
    new(bars, color, charWidth, maxFreq, maxFreqLen, symb)
  end
end

function BarplotGraphics{R<:Real}(bars::Vector{R},
                                  charWidth::Int;
                                  color::Symbol = :blue,
                                  symb = "▪")
  BarplotGraphics{R}(bars, charWidth, color, symb)
end

function addRow!{R<:Real}(c::BarplotGraphics{R}, bars::Vector{R})
  append!(c.bars, bars)
  c.maxFreq = maximum(c.bars)
  c.maxFreqLen = length(string(c.maxFreq))
  c
end

function addRow!{R<:Real}(c::BarplotGraphics{R}, bar::R)
  push!(c.bars, bar)
  c.maxFreq = max(c.maxFreq, bar)
  c.maxFreqLen = length(string(c.maxFreq))
  c
end

nrows(c::BarplotGraphics) = length(c.bars)
ncols(c::BarplotGraphics) = c.charWidth

function printRow(io::IO, c::BarplotGraphics, row::Int)
  numrows = nrows(c)
  0 < row <= numrows || throw(ArgumentError("Argument row out of bounds: $row"))
  bar = c.bars[row]
  maxBarWidth = max(c.charWidth - 2 - c.maxFreqLen, 1)
  barLen = c.maxFreq > 0 ? safeRound(bar / c.maxFreq * maxBarWidth): 0
  barStr = c.maxFreq > 0 ? repeat(c.symb, barLen): ""
  barLbl = string(bar)
  print_with_color(c.color, io, barStr)
  print_with_color(:white, io, spceStr, barLbl)
  panLen = max(safeRound(maxBarWidth + 1 + c.maxFreqLen - barLen - length(barLbl)), 0)
  pad = repeat(spceStr, panLen)
  print(io, pad)
end

function barplot{T<:(@compat AbstractString),N<:Real}(text::Vector{T}, heights::Vector{N};
                                    border = :solid,
                                    title::(@compat AbstractString) = "",
                                    margin::Int = 3,
                                    padding::Int = 1,
                                    color::Symbol = :blue,
                                    width::Int = 40,
                                    labels::Bool = true,
                                    symb = "▪")
  margin >= 0 || throw(ArgumentError("Margin must be greater than or equal to 0"))
  length(text) == length(heights) || throw(DimensionMismatch("The given vectors must be of the same length"))
  minimum(heights) >= 0 || throw(ArgumentError("All values have to be positive. Negative bars are not supported."))
  width = max(width, 5)

  area = BarplotGraphics(heights, width, color = color, symb = symb)
  newPlot = Plot(area, title=title, margin=margin,
                 padding=padding, border=border, showLabels=labels)
  for i in 1:length(text)
    annotate!(newPlot, :l, i, text[i])
  end
  newPlot
end

function barplot!{C<:BarplotGraphics,T<:(@compat AbstractString),N<:Real}(plot::Plot{C},
                                                        text::Vector{T},
                                                        heights::Vector{N};
                                                        args...)
  length(text) == length(heights) || throw(DimensionMismatch("The given vectors must be of the same length"))
  !isempty(text)|| throw(ArgumentError("Can't append empty array to barplot"))
  curIdx = nrows(plot.graphics)
  addRow!(plot.graphics, heights)
  for i = 1:length(heights)
    annotate!(plot, :l, curIdx + i, text[i])
  end
  plot
end

function barplot{T,N<:Real}(dict::Dict{T,N}; args...)
  barplot(collect(keys(dict)), collect(values(dict)); args...)
end

function barplot{T<:Real,N<:Real}(labels::Vector{T}, heights::Vector{N}; args...)
  labelsStr = map(string, labels)
  barplot(labelsStr, heights; args...)
end

function barplot{T<:Symbol,N<:Real}(labels::Vector{T}, heights::Vector{N}; args...)
  labelsStr = map(string, labels)
  barplot(labelsStr, heights; args...)
end


text:
# polyhelper function guesses what type of poly you want
export polyhelper

function polyhelper(y)
    LGRPoly(y)
end

function polyhelper(x,y)
    if length(x) < 6
        vandpoly(x,y)
    else
        LagrangePoly(x,y) # test speed of this against fitting function of Polynomials package
    end
end

#function polyhelper(x,y,b::Bound)
#    if length(x) < 6
#        vandpoly(x,y,b)
#    else
#        LagrangePoly(x,y,b)
#    end
#end

# add automatically gen GlobalPoly of local polys


text:
using CuBatchedRoutines
using Test

include("blas.jl")


text:
using BinDeps

@BinDeps.setup

mpg123 = library_dependency("libmpg123", aliases=["libmpg123", "libmpg123-0"])
lame = library_dependency("libmp3lame")

provides(AptGet, "libmpg123-dev", mpg123)
provides(AptGet, "libmp3lame-dev", lame)

provides(Pacman, "mpg123", mpg123)
provides(Pacman, "lame", lame)

@static if is_apple()
    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    provides(Homebrew.HB, "mpg123", mpg123)
    provides(Homebrew.HB, "lame", lame)
end

@static if is_windows()
    if Sys.WORD_SIZE == 32
        provides(Binaries, URI("http://www.rarewares.org/files/mp3/libmp3lame-3.99.5x86.zip"), lame,
                 SHA="7bf2a33de715d968e5ef5049b7b18e34b9ae98224842293adcdc6f6f23ad9b76",
                 unpacked_dir=".")
        provides(Binaries, URI("http://www.mpg123.org/download/win32/1.24.0/mpg123-1.24.0-x86.zip"), mpg123,
                 SHA="331a48db899f6405b7728a4fcd9ee3ef0282f1d27d983bdc8b0738314bb780cb",
                 unpacked_dir="mpg123-1.24.0-x86")
    elseif Sys.WORD_SIZE == 64
        provides(Binaries, URI("http://www.rarewares.org/files/mp3/libmp3lame-3.99.5x64.zip"), lame,
                 SHA="97959980490caa1669b1c4af386c3e23fe4ceed0835c04e328c9a5966e4d2521",
                 unpacked_dir=".")
        provides(Binaries, URI("http://www.mpg123.org/download/win64/1.24.0/mpg123-1.24.0-x86-64.zip"), mpg123,
                 SHA="d3c21a6f62204a8c93c3f48288f6f7b4d7609394daa966f2d1e6c75b30476f85",
                 unpacked_dir="mpg123-1.24.0-x86-64")
    else
        error("System not 32-bit or 64-bit. Something is wrong")
    end
end

@BinDeps.install Dict("libmpg123" => :libmpg123, "libmp3lame" => :libmp3lame)