text:
"""
    scale(f=std, x, y=f(skipmissing(x)))

Scale an array `x` by a scalar `y`.

!!! warning
    This only scales and does not center the values, unlike `scale` in R.
    See `StatsBase.zscore` for that functionality.

See also [`scale!`](@ref)
"""
scale(x) = scale(std, x)
scale(f::Function, x) = scale(x, f(skipmissing(x)))
scale(x, y) = x ./ y

"""
    scale(f=std, x, y=f(skipmissing(x)))

Scale an array `x` in place by a scalar `y`.

!!! warning
    This only scales and does not center the values, unlike `scale` in R.
    See `StatsBase.zscore` for that functionality.

See also [`scale`](@ref)
"""
scale!(x) = scale!(std, x)
scale!(f::Function, x) = scale!(x, f(skipmissing(x)))

function scale!(x, y)
    x ./= y
    return x
end

"""
    struct Scale

Represents a scaling scheme, akin to `StatsModels.AbstractContrasts`.  Pass as
value in `Dict` as hints to `schema` (or as `contrasts` kwarg for `fit`).

## Examples

Can specify scale value to use:

```
julia> schema((x=collect(1:10), ), Dict(:x => Scale(5)))
StatsModels.Schema with 1 entry:
  x => x(scaled: 5))
```

You can use a function to compute the scale value:

julia> schema((x=collect(1:10), ), Dict(:x => Scale(mad)))
StatsModels.Schema with 1 entry:
  x => x(scaled: 3.71)

Or [`scale`](@ref) will be automatically computed if left out:

```
julia> schema((x=collect(1:10), ), Dict(:x => Scale()))
StatsModels.Schema with 1 entry:
  x => x(scaled: 3.03)
```
"""
struct Scale
    scale::Any
end

Scale() = Scale(nothing)
Scale(xs::AbstractArray, s::Scale) = Scale(_standard(xs, s.scale))

"""
    struct ScaledTerm{T,S} <: AbstractTerm

A lazily scaled term.  A wrapper around an `T<:AbstractTerm` which will
produce scaled values with `modelcols` by dividing each element by `scale`.

## Fields

- `term::T`: The wrapped term.
- `scale::S`: The scale value which the resulting `modelcols` are divided by.

## Examples

Directly construct with given scale:

```
julia> d = (x=collect(1:10), );

julia> t = concrete_term(term(:x), d)
x(continuous)

julia> ts = ScaledTerm(t, 5)
x(scaled: 5)

julia> hcat(modelcols(t + ts, d)...)
10×2 Matrix{Float64}:
  1.0  0.2
  2.0  0.4
  3.0  0.6
  4.0  0.8
  5.0  1.0
  6.0  1.2
  7.0  1.4
  8.0  1.6
  9.0  1.8
 10.0  2.0
```

Construct with lazy scaling via [`Scale`](@ref)

```
julia> ts = concrete_term(term(:x), d, Scale())
x(scaled: 3.03)

julia> hcat(modelcols(t + ts, d)...)
10×2 Matrix{Float64}:
  1.0  0.330289
  2.0  0.660578
  3.0  0.990867
  4.0  1.32116
  5.0  1.65145
  6.0  1.98173
  7.0  2.31202
  8.0  2.64231
  9.0  2.9726
 10.0  3.30289
```

Or similarly via schema hints:

```
julia> sch = schema(d, Dict(:x => Scale()))
StatsModels.Schema with 1 entry:
  x => scale(x, 3.03)
```

"""
struct ScaledTerm{T,S} <: AbstractTerm
    term::T
    scale::S
end

function StatsModels.concrete_term(t::Term, xs::AbstractArray, s::Scale)
    return scale(StatsModels.concrete_term(t, xs, nothing), Scale(xs, s))
end

# run-time constructors:
scale(t::ContinuousTerm, s::Scale) = ScaledTerm(t, something(s.scale, sqrt(t.var)))
scale(t::ContinuousTerm, s) = ScaledTerm(t, s)
scale(t::ContinuousTerm) = ScaledTerm(t, sqrt(t.var))
function scale(t::AbstractTerm)
    throw(ArgumentError("can only compute scale for ContinuousTerm; must provide scale value via scale(t, s)"))
end

function scale(t::AbstractTerm, s::Scale)
    s.scale !== nothing ||
        throw(ArgumentError("can only compute scale for ContinuousTerm; must provide scale via scale(t, s)"))
    return ScaledTerm(t, s.scale)
end

StatsModels.modelcols(t::ScaledTerm, d::NamedTuple) = modelcols(t.term, d) ./ t.scale

function StatsBase.coefnames(t::ScaledTerm)
    if StatsModels.width(t.term) == 1
        return "$(coefnames(t.term))(scaled: $(_round(t.scale)))"
    elseif length(t.scale) > 1
        return string.(vec(coefnames(t.term)), "(scaled: ", _round.(vec(t.scale)), ")")
    else
        return string.(coefnames(t.term), "(scaled: ", _round(t.scale), ")")
    end
end
# coef table: "x(scaled: 5.5)"
Base.show(io::IO, t::ScaledTerm) = print(io, "$(t.term)(scaled: $(_round(t.scale)))")
# regular show: "x(scaled: 5.5)", used in displaying schema dicts
function Base.show(io::IO, ::MIME"text/plain", t::ScaledTerm)
    return print(io, "$(t.term)(scaled: $(_round(t.scale)))")
end
# long show: "x(scaled: 5.5)"

# statsmodels glue code:
StatsModels.width(t::ScaledTerm) = StatsModels.width(t.term)
# don't generate schema entries for terms which are already scaled
StatsModels.needs_schema(::ScaledTerm) = false
StatsModels.termsyms(t::ScaledTerm) = StatsModels.termsyms(t.term)
StatsModels.termvars(t::ScaledTerm) = StatsModels.termvars(t.term)


text:
struct _Hand
    hand::NTuple{4,Card}
end

_Hand(cs::Vector{Card}) = _Hand((cs...,))
_Hand(cs::Vector{<:AbstractString}) = _Hand(Card.(cs))
_Hand(cs::Card...) = _Hand((cs...,))

hand(h::_Hand) = h.hand

Base.rand(::Type{_Hand}) = _Hand(sample(DECK, 4, replace=false))

Base.sort!(h::_Hand) = sort!(hand(h))

Base.string(h::_Hand) = "(" * join(string.(hand(h)), " ") * ")"

struct Hand
    hand::_Hand
    cut::Card
end

Hand(h::Vector{<:AbstractString}, c::AbstractString) = Hand(_Hand(h), Card(c))
Hand(h::Vector{Card}, c::Card) = Hand(_Hand(h), c)
Hand(cs::Vector{Card}) = Hand(cs[1:4], cs[5])

hand(h::Hand) = hand(h.hand)
cut(h::Hand) = h.cut
cards(h::Hand) = [hand(h)..., cut(h)]

Base.rand(::Type{Hand}) = Hand(sample(DECK, 5, replace=false))

function Base.show(io::IO, h::Hand)
    print(io, "Hand([", join(string.(hand(h)), ", "), "], ", string(cut(h)), ")")
end

ext: jl


text:
struct GAPAgent
    id::Int64
    capacity::Int64
    costs::Vector{Int64}
    consumption::Vector{Int64}
end

function Base.show(io::IO, agent::GAPAgent)
    print(io, "A($(agent.id))")
end

struct GAPJob
    id::Int64
    costs::Vector{Int64}
    consumption::Vector{Int64}
end

function Base.show(io::IO, job::GAPJob)
    print(io, "J($(job.id))")
end

struct GAPData
    name::String
    agents::Vector{GAPAgent}
    jobs::Vector{GAPJob}

    lb::Int64
    ub::Int64
end

function Base.show(io::IO, data::GAPData)
    print(io, "GAP Data $(data.name)")
    print(io, " ($(length(data.agents)) agents,")
    print(io, " $(length(data.jobs)) jobs)")
    if data.lb != typemin(Int64)
        print(io, " [$(data.lb),")
    else
        print(io, " [-Inf,")
    end
    if data.ub != typemax(Int64)
        print(io, " $(data.ub)]")
    else
        print(io, " Inf]")
    end
end


text:
mutable struct MetalBufferRC{T}
    buf::MetalBuffer{T}

    freed::Bool
    refcount::Int
end

function alloc_refcount(args...; kwargs...)
    buf = alloc(args...; kwargs...)
    self = MetalBufferRC(buf, false, 0)
end

@inline function retain(a::MetalBufferRC)
  a.refcount += 1
  return
end

@inline function release(a::MetalBufferRC)
  a.refcount -= 1
  return a.refcount == 0
end

function unsafe_free!(xs::MetalBufferRC)
  # this call should only have an effect once, becuase both the user and the GC can call it
  xs.freed && return
  _unsafe_free!(xs)
  xs.freed = true
  return
end

function _unsafe_free!(xs::CuArray)
  @assert xs.refcount >= 0
  if release(xs)
    free(xs.buf)
  end
  return
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.4
#     language: julia
#     name: julia-1.6
# ---

# %%
using Distributions
using Roots
using Plots
using StatsFuns: logistic, logit

# \approx TAB → ≈
# \lessapprox TAB → ⪅
x ⪅ y = x < y || x ≈ y

# %%
# Bernoulli分布モデルでの通常の信頼区間

"""P値函数(exact版)"""
function pvalue_exact(n, p, k)
    bin = Binomial(n, p)
    p0 = pdf(bin, k)
    s = sum(pdf(bin, j) for j in support(bin) if pdf(bin, j) ⪅ p0)
    min(1, s)
end

"""P値函数(正規分布近似)"""
function pvalue_normal(n, p, k)
    bin = Binomial(n, p)
    normal = Normal(mean(bin), std(bin))
    min(1, 2cdf(normal, k), 2ccdf(normal, k))
end

"""信頼区間函数"""
function confidence_interval(pvalue_func, n, k; α = 0.05)
    f(t) = pvalue_func(n, logistic(t), k) - α
    CI = logistic.(find_zeros(f, -10.0, 10.0))
    if length(CI) < 2
        return 2k ≤ n ? [0.0, first(CI)] : [first(CI), 1.0]
    else
        return [first(CI), last(CI)]
    end
end

# %%
n = 20
x = 0:n
θ = 0:0.002:1
z_exact = pvalue_exact.(n, θ, x')
z_normal = pvalue_normal.(n, θ, x')

P = plot(; size=(500, 500), colorbar=false)
heatmap!(x, θ, z_exact; clim=(-0.1, 1.0))
plot!(; xlabel="x", ylabel="θ")
plot!(; xtick=0:n, ytick=0:0.05:1)
title!("pvalue_exact(n=$n, θ, x)"; titlefontsize=12)

Q = plot(; size=(500, 500), colorbar=false)
heatmap!(x, θ, z_normal; clim=(-0.1, 1.0))
plot!(; xlabel="x", ylabel="θ")
plot!(; xtick=0:n, ytick=0:0.05:1)
title!("pvalue_normal(n=$n, θ, x)"; titlefontsize=12)

plot(P, Q; size=(1000, 500))
plot!(; leftmargin=3Plots.mm, bottommargin=3Plots.mm)

# %%

max_stars_repo_path: 0025/heatmaps of pvalue functions.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0025/heatmaps of pvalue functions.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0025/heatmaps of pvalue functions.jl

max_forks_repo_name: genkuroki/public


text:
struct FASTA <: FileFormat end

# FASTA Parser
# ============

function _pre_readfasta(string::AbstractString)
    seqs = split(string, '>')
    N = length(seqs) - 1
    IDS  = Array{String}(undef, N)
    SEQS = Array{String}(undef, N)
    for i in 1:N
        fields = split(seqs[i+1], '\n')
        IDS[i] = fields[1]
        SEQS[i] = replace(join(fields[2:end]), r"\s+" => "")
    end
    _check_seq_len(IDS, SEQS)
    (IDS, SEQS)
end

# MethodError: no method matching seek(::TranscodingStreams.TranscodingStream, ::Int)
_pre_readfasta(io::TranscodingStreams.TranscodingStream) = _pre_readfasta(read(io, String))

function _pre_readfasta(io)
    IDS  = String[]
    SEQS = String[]
    for (name, seq) in FastaReader{String}(io) # FastaIO
        push!(IDS,  name)
        push!(SEQS, seq )
    end
    _check_seq_len(IDS, SEQS)
    (IDS, SEQS)
end

function Base.parse(io::Union{IO, AbstractString},
                    format::Type{FASTA},
                    output::Type{AnnotatedMultipleSequenceAlignment};
                    generatemapping::Bool=false,
                    useidcoordinates::Bool=false,
                    deletefullgaps::Bool=true,
                    keepinserts::Bool=false)
    IDS, SEQS = _pre_readfasta(io)
    annot = Annotations()
    _generate_annotated_msa(annot, IDS, SEQS, keepinserts, generatemapping,
        useidcoordinates, deletefullgaps)
end

function Base.parse(io::Union{IO, AbstractString},
                    format::Type{FASTA},
                    output::Type{NamedResidueMatrix{Array{Residue,2}}};
                    deletefullgaps::Bool=true)
    IDS, SEQS = _pre_readfasta(io)
    msa = _generate_named_array(SEQS, IDS)
    if deletefullgaps
        return deletefullgapcolumns(msa)
    end
    msa
end

function Base.parse(io::Union{IO, AbstractString},
                   format::Type{FASTA},
                   output::Type{MultipleSequenceAlignment};
                   deletefullgaps::Bool=true)
    msa = parse(io, format, NamedResidueMatrix{Array{Residue,2}},
                deletefullgaps=deletefullgaps)
    MultipleSequenceAlignment(msa)
end

function Base.parse(io::Union{IO, AbstractString},
                    format::Type{FASTA},
                    output::Type{Matrix{Residue}};
                    deletefullgaps::Bool=true)
    IDS, SEQS = _pre_readfasta(io)
    _strings_to_matrix_residue_unsafe(SEQS, deletefullgaps)
end

function Base.parse(io::Union{IO, AbstractString},
                    format::Type{FASTA};
                    generatemapping::Bool=false,
                    useidcoordinates::Bool=false,
                    deletefullgaps::Bool=true,
                    keepinserts::Bool=false)
    parse(io, FASTA, AnnotatedMultipleSequenceAlignment; generatemapping=generatemapping,
        useidcoordinates=useidcoordinates, deletefullgaps=deletefullgaps,
        keepinserts=keepinserts)
end

# Print FASTA
# ===========

function Base.print(io::IO, msa::AbstractMatrix{Residue}, format::Type{FASTA})
    seqnames = sequencenames(msa)
    for i in 1:nsequences(msa)
        println(io, ">", seqnames[i], "\n", stringsequence(msa, i))
    end
end

Base.print(msa::AbstractMatrix{Residue}, format::Type{FASTA}) = print(stdout, msa, FASTA)


text:
export ParAmulB!,
       AmulB!,
       AcmulB!,
       ca_add_b!

include("AmulB!.jl")
include("AcmulB!.jl")
include("a_add_b.jl")

m = 1000
A = sprand(m, m, 0.9);
x = rand(m);
y = zeros(m);
y1 = zeros(m);
nthreads = 4;
yt = zeros(m*nthreads);

@time ParAmulB!(nthreads, yt, y, A, x);
@time A_mul_B!(y1, A, x);

max_stars_repo_name: WneleiGao/ElasticWave


text:
#=
poweriter.jl
2019-06-06, Jeff Fessler, University of Michigan
=#

export poweriter

using LinearAlgebra: norm

"""
    v1,σ1 = poweriter(A; niter=?, ...)

Determine first right singular vector `v1`
and first singular value `σ1` of `A`
by applying power iteration to `A'A`

in
- `A` M × N matrix

option
- `niter` default 200
- `x0` initial guess of `v1`
- `tol` stopping tolerance for s1, default 1e-6
- `chat::Bool` verbose? default false

out
- `v1` `[N]` principal right singular vector
- `σ1` spectral norm of `A`
"""
function poweriter(
    A ;
    niter::Int=200,
    tol::Real = 1e-6,
    x0::AbstractArray{<:Number} = ones(eltype(A), size(A,2)),
    chat::Bool = true,
)

	x = copy(x0)
	ratio_old = Inf
	for iter=1:niter
		Ax = A * x
		ratio = norm(Ax) / norm(x)
		if abs(ratio - ratio_old) / ratio < tol
			chat && @info "done at iter $iter"
			break
		end
		ratio_old = ratio
		x = A' * Ax
		x /= norm(x)
	end
	return x, norm(A * x) / norm(x)
end


text:
function summarize(dd::AbstractDataFrame,Xvar::Symbol,Yvar::Symbol; Err = :MouseID , mode = :sem)
    ErrGroups = vcat(Xvar,Err)
    XaxisGroups = vcat(Xvar)
    pre_err = combine(groupby(dd, ErrGroups)) do df
        (Mean = mean(df[:,Yvar]),)
    end
    if mode == :sem
        with_err = combine(groupby(pre_err,XaxisGroups)) do df
            (Mean = mean(df.Mean), SEM = sem(df.Mean))
        end
        rename!(with_err, Xvar=>:Xaxis)
        sort!(with_err,:Xaxis)
        filter!(r -> !isnan(r.SEM), with_err)
    elseif mode == :conf
        with_err = combine(groupby(pre_err,XaxisGroups)) do df
            ci = confint(OneSampleTTest(df.Mean))
            m = mean(df.Mean)
            (Mean = m, ERRlow = m - ci[1], ERRup = ci[2] - m)
        end
        with_err[!,:ERR] = [(low,up) for (low,up) in zip(with_err.ERRlow,with_err.ERRup)]
        rename!(with_err, Xvar=>:Xaxis)
        sort!(with_err,:Xaxis)
    end
    return with_err
end

function StatsBase.ecdf(dd::AbstractDataFrame,Xvar::Symbol; Err = :MouseID, mode = :sem)
    common_xaxis = ecdf(dd[:,Xvar]).sorted_values #new
    pre_err = combine(groupby(dd, Err)) do df
        F = ecdf(df[:,Xvar])
        #(AN = F(F.sorted_values),Xaxis = F.sorted_values) #old
        (AN = F(common_xaxis),Xaxis = common_xaxis) #new
    end
    pre_err = flatten(pre_err,:AN) #results retrn in a vector within a cell
    if mode == :sem
        with_err = combine(groupby(pre_err,:Xaxis)) do df
            (Mean = mean(df.AN), ERR = sem(df.AN))
        end
        sort!(with_err,:Xaxis)
    elseif mode == :conf
        with_err = combine(groupby(pre_err,:Xaxis)) do df
            ci = confint(OneSampleTTest(df.AN))
            m = mean(df.AN)
            (Mean = m, ERRlow = m - ci[1], ERRup = ci[2] - m)
        end
        with_err[!,:ERR] = [(low,up) for (low,up) in zip(with_err.ERRlow,with_err.ERRup)]
        sort!(with_err,:Xaxis)
    end
    #dropnan!(with_err)
    return with_err
end

function effect_size(df::AbstractDataFrame,Effectvar::Symbol,Yvar::Symbol; Err = :MouseID, baseline = nothing)
    if length(union(df[:,Effectvar])) != 2
        error("Effect variable has more than 2 levels")
    end
    if baseline == nothing
        baseline = sort(union(sort(df[:,Effectvar])))[end]
    end
    gd = groupby(df,[Effectvar,Err])
    df1 = combine(Yvar => mean => :Mean,gd)
    df2 = unstack(df1,Effectvar,:Mean)
    manipulation = filter(x -> x != baseline ,union(df[:,Effectvar]))[1]
    df2[!,:Effect] = df2[:,Symbol(manipulation)] - df2[:,Symbol(baseline)]
    df2[!, :Calculation] .= string(manipulation," - ", baseline)
    deletecols!(df2,[Symbol(manipulation), Symbol(baseline)])
    return df2
end

function MVT_scatter(toplot::AbstractDataFrame; group = :Treatment)
    if !in("Average_mean",names(toplot))
        println("missing Average reward rate column")
    else
        mdl = lm(@formula(Leaving_mean ~ Average_mean),toplot)
        a,b = round.(coeftable(mdl).cols[1],digits = 3)
        label =  toplot[1, group]
        plt = @df toplot scatter(:Average_mean,:Leaving_mean,
            ylims = (0,ceil(maximum(:Leaving_mean),digits = 1)+0.1),
            xlims = (0,ceil(maximum(:Leaving_mean),digits = 1)+0.1),
            color = :grey,
            xlabel = "y = $a + $(b)x",
            title = label)
        Plots.abline!(1,0,color = :black)
        Plots.abline!(b,a,linecolor = :red, legend = false)
        return plt
    end
end

function MVT(df::AbstractDataFrame; group = :Treatment)
    gd = groupby(df,[group,:MouseID,:Day,:Trial])
    Rrate = combine([:InstRewRate, :AverageRewRate, :Reward] => (i,a,r) -> (Leaving = i[end],
        Average = a[end],
        Reward = r[end]),gd)
    filter!(r -> r.Reward, Rrate)
    gd = groupby(Rrate,[:MouseID,group])
    res = combine([:Leaving,:Average] => (l,a) ->(Leaving_mean = mean(jump_missing(l)), Average_mean = mean(jump_NaN(a))), gd)
    tp = [MVT_scatter(subdf; group = group) for subdf in groupby(res,group)]
    if group == :Treatment
        ord = [1,5,4,2,7,9,8,6,3]
        return plot(tp[ord]...)
    else
        return plot(tp...)
    end
end

function plot_wilcoxon(dd; color = nothing , sorting = nothing, showmice = true)
    if isnothing(sorting)
        sorting = Dict()
        for (v,k) in enumerate(sort(dd[:,1]))
            sorting[k] = v
        end
    end
    position = []
    original = []
    for (k,v) in sorting
        if k in dd[:,1]
            push!(position,v)
            push!(original,k)
        end
    end
    dd[!,:x] = [get(sorting,x,0) for x in dd[:,1]]
    dd2 = flatten(dd,:Vals)
    if showmice
        @df dd2 scatter(:x, :Vals,
            color = :grey,
            markeralpha = 0.4,
            markercolor = :grey)
    else
        plot()
    end
    Plots.abline!(0,0,color = :black, linestyle = :dash)
    if isnothing(color)
        c = :color in propertynames(dd) ? dd[:,:color] : :auto
    else
        c = color
    end
    @df dd scatter!(:x,:Median,
        yerror = :CI,
        linecolor = :black,
        markerstrokecolor = :black,
        markersize = 8,
        legend = false,
        tickfont = (7, :black),
        xlims = (minimum(:x) - 0.5, maximum(:x) + 0.5),
        xticks = (position,original),
        color = c,
        thickness_scaling = 2,
        guidefontsize = 8,
        ylabel = "Signed rank test - median and 95% c.i.")
end

function WebersLaw(df,x,group)
    group_vars = vcat(:MouseID,group)
    gd = groupby(df,group_vars)
    df1 = combine([x] => a -> (Mean = mean(a), STD = std(a)), gd)
    Drug_colors!(df1)
    plt = @df df1 scatter(:Mean,:STD, color = :color,
        markersize = 5, group = :Treatment, legend = :topleft,
        xlims = (0,13), ylims = (0,13))
    Plots.abline!(1,0, linestyle = :dash, xlabel = "Mean pokes after last reward", ylabel = "STD pokes after last reward", label = nothing)
    combine(groupby(df1,group)) do dd
        mdl = lm(@formula(STD ~ Mean),dd)
        a,b = round.(coeftable(mdl).cols[1],digits = 3)
        Plots.abline!(b,a, linewidth = 2, linecolor = dd[1,:color], label = "y = $(b)x + $a")
    end
    return plt
end

function plot_wilcoxon_odc(data::AbstractDataFrame,Phase::AbstractString,Allignment::Symbol; limit = 6, color = nothing)

    df0 = dropmissing(data,Allignment)
    filter!(Allignment => a -> -limit < a < limit, df0)

    df1 =  filter(r -> !ismissing(r.ODC) &&
        !r.Reward &&
        r.Phase in [Phase],
        df0)

    df2 = unstack(df1,:Treatment,:ODC)

    gd = groupby(df2,[Allignment,:MouseID])
    df3 = combine([:Control,Symbol(Phase)] => (c,d) ->
        (Control = mean(skipmissing(c)), Treatment = mean(skipmissing(d))),gd)

    df4 = combine(groupby(df3,[Allignment])) do dd
        wilcoxon(dd,:Treatment, :Control; f = x -> mean(skipmissing(x)))
    end
    plot_wilcoxon(df4; color = color)
end

function plot_wilcoxon_odc(data::AbstractDataFrame,Phases::AbstractVector{<:AbstractString},Allignment::Symbol; limit = 6, color = nothing)
    ps = []
    for phase in Phases
        c = get(drug_colors,phase,:grey)
        p = plot_wilcoxon_odc(data,phase,Allignment; limit = limit, color = c)
        title!(phase)
        ylabel!("")
        push!(ps,p)
    end
    hight_layout = Int64(round(sqrt(length(ps))))
    width_layout = Int64(ceil(length(ps) / hight_layout))

    for i in (1:width_layout:length(ps) )
        ylabel!(ps[i], "Signed rank test ODC")
    end
    plot(ps..., xlabel = string(Allignment))
end

function plot_odc(data::AbstractDataFrame,Phase::AbstractString,Allignment::Symbol; limit = 6)

    df0 = dropmissing(data,Allignment)
    filter!(Allignment => a -> -limit < a < limit, df0)

    df1 =  filter(r -> !ismissing(r.ODC) &&
        !r.Reward &&
        r.Phase in [Phase],
        df0)

    df2 = combine(groupby(df1,:Treatment)) do dd
        summarize(dd,Allignment,:ODC)
    end

    Drug_colors!(df2)
    @df df2 scatter(:Xaxis, :Mean,
        yerror = :SEM,
        group = :Treatment,
        color = :color,
        legend = false)
end

function plot_odc(data::AbstractDataFrame,Phases::AbstractVector{<:AbstractString},Allignment::Symbol; limit = 6)
    ps = []
    for phase in Phases
        p = plot_odc(data, phase,Allignment; limit = limit)
        title!(phase)
        push!(ps,p)
    end
    hight_layout = Int64(round(sqrt(length(ps))))
    width_layout = Int64(ceil(length(ps) / hight_layout))
    for i in (1:width_layout:length(ps) )
        ylabel!(ps[i], "ODC")
    end
    plot(ps..., xlabel = string(Allignment))
end

function odc_regression(data::AbstractDataFrame,phase::AbstractVector{<:AbstractString})
    df1 = filter(r -> r.Phase in phase,data)
    df2 = combine(groupby(df1,:Manipulation)) do dd
        summarize(dd, :PokeFromLastRew,:ODC)
    end
    df2[!, :color] = [get(drug_colors,r,:grey) for r in df2.Manipulation]
    @df df2 plot(:Xaxis, :Mean, ribbon = :SEM,
        group = :Manipulation,
        xlabel = "Poke from Last Reward",
        ylabel = "Omission duty cycle",
        xticks = 1:30,
        fillalpha = 0.2,
        linewidth = 3,
        color = :color,
        linecolor = :color)
end

function plot_deltafromrew(df, phase; limit = 15)
    gd = groupby(df,[:PokeFromLastRew])
    df1 = combine(gd,:DeltaCitalopram => skipmean => :DeltaCitalopram,
        :DeltaOptogenetic => skipmean => :DeltaOptogenetic,
        :DeltaAltanserin=> skipmean => :DeltaAltanserin,
        :DeltaSB242084 => skipmean => :DeltaSB242084)
    c_name = "Delta" * phase
    c_symbol = Symbol(c_name)
     @df filter(r -> r.PokeFromLastRew < limit, df) scatter(:PokeFromLastRew, cols(c_symbol), markercolor = :grey)
    @df filter(r -> r.PokeFromLastRew < limit, df1) scatter!(:PokeFromLastRew, cols(c_symbol), markersize = 8, markercolor = get(drug_colors,phase,:red))
    Plots.abline!(0,0, legend = false,
        ylims = (-0.1,0.1),
        xticks = 1:20, xlabel = "Poke from last reward", ylabel = "ODC delta " * phase)
end

function plot_QODC(df, allignment; xflip = false, legend = false, ylims = (3,12))
    gd1 = groupby(df,[:Phase,:MouseID,allignment,:Treatment])
    df1 = combine(gd1, :QODC => mean)
    gd2 = groupby(df1,[:Phase,allignment,:Treatment])
    df2 = combine(gd2, :QODC_mean => mean => :QODC_mean, :QODC_mean => sem => :QODC_sem)
    plts = []
    Drug_colors!(df2)
    combine(groupby(df2,:Phase)) do dd
        p = @df dd scatter(cols(allignment), :QODC_mean, yerror = :QODC_sem,
        group = :Treatment, color = :color, xflip = xflip, legend = legend,
        xticks = 0:20, yticks = 0:20, ylims = ylims,
        markersize = 4,
        thickness_scaling = 2,
        guidefontsize = 8,
        ylabel = "Omissions duty cycle quantile")
        push!(plts,p)
    end
    plot(plts...)
end


text:
#import GeometryBasics: Point, Point3, Triangle
import GeometryBasics: Point, Point3, Point2, Vec, Vec2, Vec3, TriangleFace, Rectangle
using LinearAlgebra
import StaticArrays
import GeneralPolygonClipper as gpc
using VoronoiCells

struct SProject2d end

function boundingbox2(pts)
    i = firstindex(pts)
    xmin = pts[i][1]
    ymin = pts[i][2]

    xmax = xmin
    ymax = ymin

    for p in pts
        x,y = p
        if x < xmin
            xmin = x
        elseif x > xmax
            xmax = x
        end
        if y < ymin
            ymin = y
        elseif y > ymax
            ymax = y
        end
    end
    return (xmin, ymin), (xmax, ymax)

end

function boundingbox3(pts)

    i = firstindex(pts)
    xmin = pts[i][1]
    ymin = pts[i][2]
    zmin = pts[i][3]

    xmax = xmin
    ymax = ymin
    zmax = zmin

    for p in pts
        x,y,z = p
        if x < xmin
            xmin = x
        elseif x > xmax
            xmax = x
        end
        if y < ymin
            ymin = y
        elseif y > ymax
            ymax = y
        end
        if z < zmin
            zmin = z
        elseif z > zmax
            zmax = z
        end

    end
    return (xmin, ymin, zmin), (xmax, ymax, zmax)

end

function get_basis(n::Vec{3,T}) where {T}
    # We want to project the polygons on a surface plane
    # but we still need two two independent vector to form
    # an orthonormal basis. It would be better if the
    # axis are partially aligned with the original
    # 3D reference frame.

    ex = Vec{3,T}(1, 0, 0)
    ey = Vec{3,T}(0, 1, 0)
    ez = Vec{3,T}(0, 0, 1)

    e = (ex, ey, ez)
    dx = ex ⋅ n
    dy = ey ⋅ n
    dz = ez ⋅ n

    # If the surface is almost horizontal idxmax=3
    idxmax = argmax( abs.((dx, dy, dz)) )
    if idxmax==3
        idx = 2 # Y axis will be vertical
    else
        # Figure out which axis is more perpendicular to the normal:
        idx = argmin( abs.((dx, dy, dz)) )
    end
    # idx corresponds to the axis more nearly perpendicular.
    # Since in most buildings most faces are vertical or
    # approximately so, the upward direction will be the y axis
    # for the projection.

    vy = e[idx] - (e[idx] ⋅ n) * n  # Vector normal to n and e[idx]
    vx = vy × n

    # return the new *orthonormal* basis
    return vx ./ hypot(vx...), vy ./ hypot(vy...), n ./ hypot(n...)

end

"""
`mean_mesh_normal(msh)`

Compute the mean normal of a mesh
"""
function mean_mesh_normal(msh)

    np = length(msh) # Number of polygons (or triangles...)

    nrm <- sum(normal(p) for p in msh)
    n⃗ = nrm ./ hypot(nrm...) # Unit vector

    return get_basis(n⃗)

end

normal1(x) = (x[2]-x[1]) × (x[3]-x[1]) / 2

"""
`discrsurface(type, faces, pext, fext, pint, fint)`

Building surface discretization.

This function discretizes a face of a building around
 internal and external pressure taps. A face is made up of subfaces
which are made up of triangles.

 * `type` A type specifying the algorithm used
 * `faces` a collection of faces
 * `pext` External pressure taps
 * `fext` Index of faces where the pressure taps are located
 * `pint` Internal pressure taps
 * `fint` Index of faces where internal taps are located. Should be a subset of `fext`
"""
function discrsurface(::Type{SProject2d}, faces, pext, fext,  pint, fint)
    # Check if `fint` is an subset of `fext`
    for f in fint
        if f ∉ fext
            error("External surfaces should include every internal surface!")
        end
    end

    # Calculate the mean normal of all external faces
    nrm = Vec3{Float64}(0,0,0)
    for f in faces[fext]
        for tri in f
            n = normal1(tri)
            nrm += n
        end
    end

    nrm = nrm / norm(n)

    # Project points and faces on a plane normal to `nrm`
    e = get_basis(nrm) # Basis

    P₀ = Point3(0,0,0)

    # Discretize external surface:
    tri_ext = partition_faces(pext, faces, fext, e, P₀)

    # Discretize internal surfaces, if there are any
    if length(pint) > 0
        # Get the mesh
        tri_int = partition_faces(pint, faces, fint, e, P₀)

        # Now we should intersec
    end

end

"""
`partition_faces(pts, facelist, fidx, e, p0)`

Partition a collection of faces in regions corresponding to a collection of points.

The faces are collections of triangles and they are projected into a plane defined
by point `p0` and normal `e[3]`. This plane has two direction vectors:

 * `e[1]` that correspnds to the 2d axis direction vector î
 * `e[2]` that correspnds to the 2d axis direction vector ĵ

Once the points and faces are projected into the plane, a Voronoi tesselation
is carried out with respect to the points and each Voronoi cell is intersected
with each triangle that compose the faces. This intersection is represented
as a set of triangles.

## Arguments

 * `pts` Vector of points
 * `facelist` A vector containing every face in the model
 * `fidx` Vector with indices of faces that should be processed
 * `e` Basis defining the projection plane
 * `p0` Origin of the 2d coordinate system on the projection plane

## Return value

This function returns a vector with a tuple containing a triangle, the index of
the point, the index of the face that this triangle belongs to and the index of the
model triangle to which the triangle belongs.

"""
function partition_faces(pts, facelist, fidx, e, p0)

    # Let's work with the necessary faces only!
    faces = facelist[fidx]

    # A plane is specified by a point and a normal. Just use the origin as the point
    P₀ = p0
    # Project the triangles
    faces2d = [[project_triangle(e, P₀, tri3) for tri3 in ff] for ff in faces]

    # Project the points on thge same plane
    pts2 = [Point2( (p - P₀) ⋅ e[1], (p-P₀)⋅e[2] ) for p in pts]

    # Get a bounding box of all the projected surfaces to create the Voronoi tesselation
    xmin, ymin = faces2d[1][1][1]
    xmax, ymax = xmin, ymin
    for f in faces2d
        for tri in f
            for i in 1:3
                x,y = tri[i]
                if x < xmin
                    xmin = x
                elseif x > xmax
                    xmax = x
                end
                if y < ymin
                    ymin = y
                elseif y > ymax
                    ymax = y
                end
            end
        end
    end
    Δx = xmax - xmin; Δy = ymax - ymin

    # Bounding box. Let's make it bigger to ensure everything fits...
    rect = Rectangle(Point2(xmin-Δx, ymin-Δy), Point2(xmax+Δx, ymax+Δy))

    # Voronoi tesselation of external points
    tessex = voronoicells(pts2, rect)

    # Array to store the 2d triangles, the face index and triangle index.
    tri2d = Tuple{TriangleFace{Point2{Float64}}, Int, Int, Int}[]

    # Let's go through every Voronoi cell and intersect it with every triangle
    # in the faces
    for i in eachindex(pts)
        vpts = tessex.Cells[i] # Voronoi cell
        (xmin1,ymin1), (xmax1,ymax1) = boundingbox2(vpts)
        # Create a GenericPolygonClipper.jl polygon
        pgpc = gpc.GPCPolygon([false], [[gpc.Vertex(v[1], v[2]) for v in vpts]])
        for (k, face)   in enumerate(faces2d)
            for (j, tri) in enumerate(face)
                (xmin2, ymin2), (xmax2, ymax2) = boundingbox2(tri)
                if !(xmax2 < xmin1 || xmax1 < xmin2 ||
                    ymax2 < ymin1 || ymax1 < ymin2)
                    # There *might* be an intersection!
                    # If there is an intersection, we want it decomposed in triangle
                    # strips.
                    trigpc = gpc.GPCPolygon([false],
                                            [[gpc.Vertex(tri[i][1], tri[i][2]) for i in 1:3]])
                    # Now we carry out the intersection between the Voronoi cell
                    # and the face triangle
                    trilst = gpc.intersect_strip(pgpc, trigpc)
                    for s in trilst
                        for l in 1:length(s)
                            vv = s[l]
                            tri_new = TriangleFace(Point2{Float64}(vv[1].x, vv[1].y),
                                                   Point2{Float64}(vv[2].x, vv[2].y),
                                                   Point2{Float64}(vv[3].x, vv[3].y))
                            push!(tri2d, (tri_new, i, fidx[k], j))
                        end
                    end

                end

            end
        end
    end
    # Project back to 3d.
    # And remember that the orientation might have been inverted
    return [ (project_back(e, P₀, x[1], facelist[x[3]][x[4]]),
              x[2], x[3], x[4]) for x in tri2d ]

end

"""
`project_triangle(e, p, tri3)`

Project a 3d triangle into a plane defined by basis `e` with origin `p`.

## Arguments

 * `e` Tuple containing the basis. `e[1]` and `e[2]` are the vectors definind the plane and `e[3]` is the vector normal to the plane
 * `p` 3d point which is used as origin in the plane
 * `tri3` 3d Triangle that should be projected.

## Return value

A `2d TriangleFace` object containing the projected triangle.
"""
function project_triangle(e, p, tri3)
    v1 = Point2( (tri3[1]-p) ⋅ e[1], (tri3[1]-p) ⋅ e[2])
    v2 = Point2( (tri3[2]-p) ⋅ e[1], (tri3[2]-p) ⋅ e[2])
    v3 = Point2( (tri3[3]-p) ⋅ e[1], (tri3[3]-p) ⋅ e[2])
    return TriangleFace(v1, v2, v3)

end

"""
`project_back(e, p, tri2d, tri0)`

Project back a 2d triangle into the 3d surface of the model.

This function projects a 2d triangle back into 3d space where the triangle
belongs to the surface defined by 3d triangle `tri0`.

The 2d plane where `tri2d` is defined has as origin the point `p` (3d) and
unit direction vectors `î = e[1]` and `ĵ = e[2]`.

The function will reorder the output vertices so that it has the same orientation
as the 3d triangle `tri0`.

## Arguments
 * `e` Tuple with direction vectors defining the projection plane
 * `p` 3d point which is used as origin in the 2d projection plane
 * `tri2d` The 2d triangle that should be projected back
 * `tri0` triangle that defines the 3d plane where the output triangle should be located.

## Return value

A `TriangleFace` corresponding to `tri2d` in the plane defined by `tri0`.
"""
function project_back(e, p, tri2d, tri0)
    # Any point P in the same plane of a triangle
    # can be computed as P - P₁ = α⋅u + β⋅v
    # Where:
    #  * P₁ is the first vertex of the triangle
    #  * u is the vector from P₁ to P₂
    #  * v is the vector from P₁ to P₃
    # In the projected plane (basis `e`), this results
    # in a system of equations M = [ux vx; uy vy]
    # M⋅[α; β] = [x-x₁, y-y₁]
    # Inverting M we can calculate α and β and now
    # we can use the original definintion above to compute P.

    # Project first vertex on the plane
    vx₁ = (tri0[1]-p) ⋅ e[1]
    vy₁ = (tri0[1]-p) ⋅ e[2]

    u = tri0[2] - tri0[1]
    v = tri0[3] - tri0[1]
    ux = u ⋅ e[1]
    uy = u ⋅ e[2]

    vx = v ⋅ e[1]
    vy = v ⋅ e[2]

    Δ = ux * vy - uy * vx

    A11 =  vy / Δ; A12 = -vx / Δ
    A21 = -uy / Δ; A22 =  ux / Δ

    δx = tri2d[1][1] - vx₁
    δy = tri2d[1][2] - vy₁
    α = A11*δx + A12*δy
    β = A21*δx + A22*δy
    p1 = tri0[1] + α * u + β * v

    δx = tri2d[2][1] - vx₁
    δy = tri2d[2][2] - vy₁
    α = A11*δx + A12*δy
    β = A21*δx + A22*δy
    p2 = tri0[1] + α * u + β * v

    δx = tri2d[3][1] - vx₁
    δy = tri2d[3][2] - vy₁
    α = A11*δx + A12*δy
    β = A21*δx + A22*δy
    p3 = tri0[1] + α * u + β * v

    # Check orientation of the nodes.
    # If the normals are in the same direction, keep the order p1, p2, p3.
    # Otherwise invert: p1, p3, p2.
    if (u × v) ⋅ ( (p2-p1)×(p3-p1) ) > 0
        return TriangleFace(p1, p2, p3)
    else
        return TriangleFace(p1, p3, p2)
    end

end


text:
using Clang

GAGE_INCLUDE = raw"C:\Program Files (x86)\Gage\CompuScope\include"

clang_includes = String[]
push!(clang_includes, GAGE_INCLUDE)
push!(clang_includes, raw"C:\Program Files\LLVM\include\clang-c",raw"C:\Program Files\LLVM\include\llvm-c" )
clang_extraargs = ["-v"]
clang_extraargs = ["-D", "__STDC_CONSTANT_MACROS","-D", "__STDC_LIMIT_MACROS"]
function wrap_header(top_hdr, cursor_hdr)
    return true#startswith(dirname(cursor_hdr), GAGE_INCLUDE)
end

function lib_file(hdr)
     if startswith(hdr, "CsAppSupport")
         return :CsSsm
     else
         return :CsSsm
     end
 end

output_file(hdr) = "Gage.jl"

const wc = wrap_c.init(;
headers = [joinpath(GAGE_INCLUDE,"CsAppSupport.h")],
output_file= "gage_headers.jl",
common_file = "cs_common.jl",
clang_includes = clang_includes,
clang_args = clang_extraargs,
header_wrapped = wrap_header,
header_library = lib_file,
header_outputfile = output_file)
run(wc)

# @show Clang.wrap_c.debug_cursors


text:
using JuMP,LightGraphs,DataStructures,DelimitedFiles,LinearAlgebra
using Gurobi

include("data_MMM.jl")
include("graph.jl")
include("combi_k_resist.jl")

#-----------------------------------------------------------------------------------

function create_model()
  if Optimizer=="CPLEX"
    M = Model(with_optimizer(CPLEX.Optimizer))
    MOI.set(M, MOI.RawParameter("CPX_PARAM_THREADS"), 1)
    MOI.set(M, MOI.RawParameter("CPX_PARAM_SCRIND"), 0)
    MOI.set(M, MOI.RawParameter("CPX_PARAM_MIPDISPLAY"), 1)
    #MOI.set(M, MOI.RawParameter("CPX_PARAM_TILIM"), max(0,TIME_LIM))
  elseif Optimizer=="Gurobi"
    M = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(GUROBI_ENV)))
    set_optimizer_attribute(M, "OutputFlag", 0)
    #set_optimizer_attribute(M, "TimeLimit", TIME_LIM)
    set_optimizer_attribute(M, "Threads", 1)
    #M = Model(with_optimizer(Gurobi.Optimizer,GUROBI_ENV,OutputFlag= 0,TimeLimit=max(0,TIME_LIM), Threads=1))
  end
  return M
end

#-----------------------------------------------------------------------------------

const Optimizer = "Gurobi"
const GUROBI_ENV = Gurobi.Env()

TIME_LIMIT = 7200
problem = ("SP")
separation = "binary"
count1, count2, count3 = 0, 0, 0 # used to assess the quality of the upper bounds on the objective function
minmaxmincost = 0
LB = 0
maxmintime, inter_computation_time, test_time = 0, 0, 0

appfolder = dirname(@__FILE__)

heuristic = Dict()
for K in [2,3,4], Γ in [3,6], n in 20:5:30
  file = readdlm("$(dirname(@__FILE__))/heuristic_results/mip_vs_heuristic_n_$(n)_B_$(Γ)_K_$(K).txt")
  for i in 1:100
      heuristic[K,n,Γ,i] = file[i,2]
  end
end

K = parse(Int, ARGS[1])
Γ = parse(Int, ARGS[2])
discrete = ARGS[3] == "discrete"
convex = ARGS[3] == "convex"
#convex = true
#outputfilename2 = "res/SP_resistLP_$(ARGS[3])_B_$(Γ)_K_$(K).csv"
outputfilename2 = "res/SP_resistLP_convex_B_$(Γ)_K_$(K).csv"
#outputfilename = "res/SP_$(ARGS[3])_B_$(Γ)_K_$(K).csv"
n1,n2 = 20, 30
i1, i2 = 1, 100
if length(ARGS) > 3
    n1 = parse(Int, ARGS[4])
    n2 = parse(Int, ARGS[5])
    i1 = parse(Int, ARGS[6])
    i2 = parse(Int, ARGS[7])
end
granularity = 1 # If greater than 1, means fractional deviations are considered in X_ω, e.g. = 2 => deviations of 0.5 are considered. Keep this equal to 1 for the discrete model, otherwise not valid
if convex
    granularity = 10
else
    granularity = 1
end
#@info "Solving $(ARGS[3]) problems with K=$K and Γ=$Γ, granularity=$granularity"
@info "Solving convex problems with K=$K and Γ=$Γ, granularity=$granularity"
σ = Vector{Int}(undef,K)
ω = Vector{Int}(undef,K)
σ_first = Vector{Int}(undef,K)
ω_first = Vector{Int}(undef,K)
σ_last = Vector{Int}(undef,K)
ω_last = Vector{Int}(undef,K)
@warn "warming up"
STARTTIME = time()
data = read_dataMMM_SP("/../data/SP/1_20.txt")
algo_combi_k_resist(heuristic[2,20,3,1],data)
@warn "main loop"
for n in n1:5:n2, i in i1:i2
    GC.gc()
    instance = "/../data/SP/$(i)_$(n).txt"
    global data = read_dataMMM_SP(instance)
    global starttime = time()
    global STARTTIME = time()

    (runtime, opt, r, X_time, prep_time, count1, count2, count3, maxmintime, inter_computation_time) = algo_combi_k_resist(heuristic[K,n,Γ,i],data)
    solved = runtime < TIME_LIMIT
    if !solved
        runtime = TIME_LIMIT
    end
    output = "$solved $n $i $opt $r $count1 $count2 $count3 $runtime $(round(100*maxmintime/runtime)) $(round(100*inter_computation_time/runtime)) $(round(100*X_time/runtime)) $(round(100*prep_time/runtime))"
    @info output
    outputfile = open(outputfilename2, "a")
    println(outputfile, output)
    close(outputfile)
end

max_stars_repo_name: mjposs/min-max-min

max_issues_repo_name: mjposs/min-max-min

max_forks_repo_name: mjposs/min-max-min


text:
using vOptSpecific

length(ARGS) != 1 && error("run julia parseLAP.jl filename")
!isfile(ARGS[1]) && error("can't find $(ARGS[1])")

id = load2LAP(ARGS[1])
z1, z2, solutions = vSolve(id)

for i = 1:length(z1)
    println("(" , z1[i], " | ", z2[i], ") : ", solutions[i,:])
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
using StructuresKit

#Ix Iy Ixy J Cw
section_properties = [(3.230E6,449530,-865760, 397.09, 3.4104E9)]

#E  ν
material_properties = [(200,0.30)]

#kx kϕ
spring_stiffness = [(0.0,300/1000)]

#ay_kx
spring_location = [(101.6)]

#qx qy start end
loads = [(0.0001, 0.00002, 0.0, 7620.0)]

#ax ay
load_locations = [(27.826,101.6)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
end_boundary_conditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 7620]

#member information
#L(1) dL(2) section_properties(3) material_properties(4) spring_stiffness(5) spring_location(6) load(7) load_location(8)
member_definitions = [(7620, 7620/12, 1, 1, 1, 1, 1, 1)]

#Discretize beam.
z, m = ThinWalledBeam.discretize(member_definitions)

#Apply loads.
num_nodes = length(z)
qx = 0.001 * ones(num_nodes)
qy = 0.001 * ones(num_nodes)
ax = 27.8 * ones(num_nodes)
ay = 101.6 * ones(num_nodes)

#Define springs.
kx = 0.0 * ones(num_nodes)
kϕ = 0.3 * ones(num_nodes)
ay_kx = 101.6 * ones(num_nodes)

#Define model.
model = ThinWalledBeam.define(z, m, member_definitions, section_properties, material_properties, kx, kϕ, ay_kx, qx, qy, ax, ay, end_boundary_conditions, supports)

#Solve model.
model = ThinWalledBeam.solve(model)

max_stars_repo_name: runtosolve/PurlinLine

max_issues_repo_name: runtosolve/PurlinLine

max_forks_repo_name: runtosolve/PurlinLine


text:
# Defines structus and functions relating to display.
# Defines `render` for all game objects.

struct ScreenPixelPos  # 0,0 == top-left
    x::Int
    y::Int
end
ScreenPixelPos(x::Number, y::Number) = ScreenPixelPos(convert.(Int, floor.((x,y)))...)
struct UIPixelPos  # 0,0 == top-left (Same as ScreenPixelPos but not scaled.)
    x::Int
    y::Int
end
UIPixelPos(x::Number, y::Number) = UIPixelPos(convert.(Int, floor.((x,y)))...)

+(a::UIPixelPos, b::UIPixelPos) = UIPixelPos(a.x+b.x, a.y+b.y)
# TODO: 1.0?
+(a::UIPixelPos, x::Number) = UIPixelPos(a.x+x, a.y+x)

mutable struct Camera
    pos::WorldPos
    w::Threads.Atomic{Int32}   # Note: These are Atomics, since they can be modified by the
    h::Threads.Atomic{Int32}   # windowEventWatcher callback, which can run in another thread!
end
Camera() = Camera(WorldPos(0,0),100,100)

screenCenter() = UIPixelPos(winWidth[]/2, winHeight[]/2)
screenCenterX() = winWidth[]/2
screenCenterY() = winHeight[]/2
screenOffsetFromCenter(x::Int,y::Int) = UIPixelPos(screenCenterX()+x,screenCenterY()+y)

worldScale(c::Camera) = cam.w[] / winWidth[];
function toScreenPos(p::WorldPos, c::Camera)
    scale = worldScale(c)
    ScreenPixelPos(
        floor(c.w[]/2.0 + scale*p.x), floor(c.h[]/2.0 - scale*p.y))
end
function toScreenPos(p::UIPixelPos, c::Camera)
    scale = worldScale(c)
    ScreenPixelPos(floor(scale*p.x), floor(scale*p.y))
end
function screenToUI(p::ScreenPixelPos, c::Camera)
    scale = worldScale(c)
    ScreenPixelPos(floor(p.x/scale), floor(p.y/scale))
end
function screenScaleDims(w,h,c::Camera)
    scale = worldScale(c)
    scale*w, scale*h
end

function render(o::Ball, cam::Camera, renderer)
    ballW = ballWidth; ballH = ballWidth;
    topLeft = WorldPos(o.pos.x - ballW/2.0, o.pos.y + ballH/2.0)
    screenPos = toScreenPos(topLeft, cam)
    rect = SDL2.Rect(screenPos.x, screenPos.y, screenScaleDims(ballW, ballH, cam)...)
    color = kBallColor
    SDL2.SetRenderDrawColor(renderer, Int64(color.r), Int64(color.g), Int64(color.b), Int64(color.a))
    SDL2.RenderFillRect(renderer, Ref(rect) )
end
function render(o::Paddle, cam::Camera, renderer)
    paddleW = o.length; paddleH = kPaddleRenderH;
    # Move up/down so the edge matches the "center" of paddle.
    edgeShift = if o.pos.y > 0; paddleH/2.0; else -paddleH/2.0; end
    topLeft = WorldPos(o.pos.x - paddleW/2.0, o.pos.y + paddleH/2.0 + edgeShift)
    screenPos = toScreenPos(topLeft, cam)
    rect = SDL2.Rect(screenPos.x, screenPos.y, screenScaleDims(paddleW, paddleH, cam)...)
    color = kPaddleColor
    SDL2.SetRenderDrawColor(renderer, Int64(color.r), Int64(color.g), Int64(color.b), Int64(color.a))
    SDL2.RenderFillRect(renderer, Ref(rect) )
end

abstract type AbstractButton end
mutable struct MenuButton <: AbstractButton
    enabled::Bool
    pos::UIPixelPos
    w::Int
    h::Int
    text::String
    callBack
end
mutable struct KeyButton <: AbstractButton
    enabled::Bool
    pos::UIPixelPos
    w::Int
    h::Int
    text::String
    callBack
end

mutable struct CheckboxButton
    toggled::Bool
    button::MenuButton
end

import Base.run
run(b::AbstractButton) = b.callBack()
function run(b::CheckboxButton)
    b.toggled = !b.toggled
    b.button.callBack(b.toggled)
end

# pointwise subtraction with bounds checking (floors to 0)
-(a::SDL2.Color, b::Int) = SDL2.Color(a.r-min(b,a.r), a.g-min(b,a.g), a.b-min(b,a.b), a.a-min(b,a.a))
SDL2.Color(1,5,1,1) - 2 == SDL2.Color(0,3,0,0)

function render(b::AbstractButton, cam::Camera, renderer, color, fontSize)
    if (!b.enabled)
         return
    end
    topLeft = UIPixelPos(b.pos.x - b.w/2.0, b.pos.y - b.h/2.0)
    screenPos = toScreenPos(topLeft, cam)
    rect = SDL2.Rect(screenPos.x, screenPos.y, screenScaleDims(b.w, b.h, cam)...)
    x,y = Int[0], Int[0]
    SDL2.GetMouseState(pointer(x), pointer(y))
    if clickedButton == b
        if mouseOnButton(UIPixelPos(x[],y[]),b,cam)
            color = color - 50
        else
            color = color - 30
        end
    else
        if mouseOnButton(UIPixelPos(x[],y[]),b,cam)
            color = color - 10
        end
    end
    SDL2.SetRenderDrawColor(renderer, Int64(color.r), Int64(color.g), Int64(color.b), 255)
    SDL2.RenderFillRect(renderer, Ref(rect) )
    renderText(renderer, cam, b.text, b.pos; fontSize = fontSize)
end

function render(b::MenuButton, cam::Camera, renderer)
    render(b, cam, renderer, kMenuButtonColor, kMenuButtonFontSize)
end
function render(b::KeyButton, cam::Camera, renderer)
    render(b, cam, renderer, kKeySettingButtonColor, kKeyButtonFontSize)
end

function render(b::CheckboxButton, cam::Camera, renderer)
    # Hack: move button text offcenter before rendering to accomodate checkbox
    offsetText = " "
    text_backup = b.button.text
    b.button.text = offsetText * b.button.text
    render(b.button, cam, renderer)
    b.button.text = text_backup

    # Render checkbox
    render_checkbox_square(b.button, 6, SDL2.Color(200,200,200, 255), cam, renderer)

    if b.toggled
        # Inside checkbox "fill"
        render_checkbox_square(b.button, 8, SDL2.Color(100,100,100, 255), cam, renderer)
    end
end

function render_checkbox_square(b::AbstractButton, border, color, cam, renderer)
    checkbox_radius = b.h/2.0 - border  # (checkbox is a square)
    topLeft = UIPixelPos(b.pos.x - b.w/2.0, b.pos.y - b.h/2.0)
    topLeft = topLeft + border
    screenPos = toScreenPos(topLeft, cam)
    rect = SDL2.Rect(screenPos.x, screenPos.y, screenScaleDims(checkbox_radius*2, checkbox_radius*2, cam)...)
    SDL2.SetRenderDrawColor(renderer, Int64(color.r), Int64(color.g), Int64(color.b), 255)
    SDL2.RenderFillRect(renderer, Ref(rect) )
end

# ---- Text Rendering ----

fonts_cache = Dict()
txt_cache = Dict()
function sizeText(cam, txt, fontName, fontSize)
    scale = worldScale(cam)
    sizeText(scale, txt, loadFont(scale, fontName, fontSize))
end
function sizeText(scale, txt, font::Ptr{SDL2.TTF_Font})
   fw,fh = Cint[1], Cint[1]
   SDL2.TTF_SizeText(font, txt, pointer(fw), pointer(fh))
   return fw[1]/scale,fh[1]/scale
end
function loadFont(scale, fontName, fontSize)
   fontSize = scale*fontSize
   fontKey = (fontName, fontSize)
   if haskey(fonts_cache, fontKey)
       font = fonts_cache[fontKey]
   else
       font = SDL2.TTF_OpenFont(fontKey...)
       font == C_NULL && throw(ErrorException("Failed to load font '$fontKey'"))
       fonts_cache[fontKey] = font
   end
   return font
end
function createText(renderer, cam, txt, fontName, fontSize)
   font = loadFont(worldScale(cam), fontName, fontSize)
   txtKey = (font, txt)
   if haskey(txt_cache, txtKey)
       tex = txt_cache[txtKey]
   else
       text = SDL2.TTF_RenderText_Blended(font, txt, SDL2.Color(20,20,20,255))
       tex = SDL2.CreateTextureFromSurface(renderer,text)
       #SDL2.FreeSurface(text)
       txt_cache[txtKey] = tex
   end

   fw,fh = Cint[1], Cint[1]
   SDL2.TTF_SizeText(font, txt, pointer(fw), pointer(fh))
   fw,fh = fw[1],fh[1]

   return tex, fw, fh
end
@enum TextAlign centered leftJustified rightJustified
function renderText(renderer, cam::Camera, txt::String, pos::UIPixelPos
                    ; fontName = defaultFontName,
                     fontSize=defaultFontSize, align::TextAlign = centered)
   tex, fw, fh = createText(renderer, cam, txt, fontName, fontSize)
   renderTextSurface(renderer, cam, pos, tex, fw, fh, align)
end

function renderTextSurface(renderer, cam::Camera, pos::UIPixelPos,
                           tex::Ptr{SDL2.Texture}, fw::Integer, fh::Integer, align::TextAlign)
   screenPos = toScreenPos(pos, cam)
   renderPos = SDL2.Rect(0,0,0,0)
   if align == centered
       renderPos = SDL2.Rect(Int(floor(screenPos.x-fw/2.0)), Int(floor(screenPos.y-fh/2.0)), fw,fh)
   elseif align == leftJustified
       renderPos = SDL2.Rect(Int(floor(screenPos.x)), Int(floor(screenPos.y-fh/2.0)), fw,fh)
   else # align == rightJustified
       renderPos = SDL2.Rect(Int(floor(screenPos.x-fw)), Int(floor(screenPos.y-fh/2.0)), fw,fh)
   end
   SDL2.RenderCopy(renderer, tex, C_NULL, pointer_from_objref(renderPos))
   #SDL2.DestroyTexture(tex)
end

function hcat_render_text(lines, renderer, cam, gap, pos::UIPixelPos;
         fixedWidth=nothing, fontName=defaultFontName, fontSize=defaultFontSize)
    numLines = size(lines)[1]
    if fixedWidth != nothing
        widths = fill(fixedWidth,size(lines))
    else
        widths = [sizeText(cam, line, fontName, fontSize)[1] for line in lines]
    end
    totalWidth = sum(widths) + gap*(numLines-1)
    runningWidth = 0
    leftMostPos = pos.x - totalWidth/2.0
    text_centers = []
    for i in 1:numLines
        linePos = leftMostPos + runningWidth
        leftPos = UIPixelPos(linePos, pos.y)
        renderText(renderer, cam, lines[i], leftPos; fontName=fontName, fontSize=fontSize, align=leftJustified)
        runningWidth += widths[i] + gap
        push!(text_centers, UIPixelPos(leftPos.x + widths[i]÷2, leftPos.y))
    end
    return text_centers
end

#  ------- Image rendering ---------

function render(t::Ptr{SDL2.Texture}, pos, cam::Camera, renderer; size = nothing)
    if (t == C_NULL) return end
    pos = toScreenPos(pos, cam)
    if size != nothing
        size = toScreenPos(size, cam)
        w = size.x
        h = size.y
    else
        w,h,access = Cint[1], Cint[1], Cint[1]
        format = Cuint[1]
        SDL2.QueryTexture( t, format, access, w, h );
        w,h = w[], h[]
    end
    rect = SDL2.Rect(pos.x - w÷2,pos.y - h÷2,w,h)
    SDL2.RenderCopy(renderer, t, C_NULL, pointer_from_objref(rect))
end

max_stars_repo_name: NHDaly/pongClone

max_issues_repo_name: NHDaly/pongClone

max_forks_repo_name: NHDaly/pongClone


text:
# build factor graph from sandshark data
addprocs(7)
using Caesar
using RoME
using IncrementalInference, KernelDensityEstimate
using MAT
using JLD, HDF5
using Gadfly

vars = matread("/home/dehann/data/sandshark/cbf.mat");

len = length(vars["imu_gyr_x"]);

accs = zeros(len,3);
gyrs = zeros(len,3);

lbl = ["x";"y";"z"]
for i in 1:3 accs[:,i] = vars["imu_acc_$(lbl[i])"][:]; end
for i in 1:3 gyrs[:,i] = vars["imu_gyr_$(lbl[i])"][:]; end
imutime = vars["imu_time"];

nativeOdo = zeros(length(vars["nav_time"][:]), 5);
nativeOdo[:,1] = vars["nav_time"][:];
nativeOdo[:,2] = 0.68*vars["nav_speed"][:];
nativeOdo[:,3] = vars["nav_roll"][:]*pi/180.0;
nativeOdo[:,4] = vars["nav_pitch"][:]*pi/180.0;
nativeOdo[:,5] = (vars["nav_yaw"][:]-14.72)*pi/180.0; # assumed to be heading

nativeOdo[:,1] -= vars["nav_time"][1]
nativeOdo[:,1] *= 1e-9

for i in 1:size(nativeOdo,1)
  nativeOdo[i,5] = wrapRad( (2pi-nativeOdo[i,5])-3.0*pi/2.0 )
end

nativeOdo[2:end,5] = diff(nativeOdo[:,5])
nativeOdo[1,5] = 0
nativeOdo[nativeOdo[:,5] .> 2pi-0.1, 5] -= 2pi
nativeOdo[nativeOdo[:,5] .< -2pi+0.1, 5] += 2pi

nativeMeas = zeros(length(vars["CBF_MF_time"]),5);
nativeMeas[:,1] = (vars["CBF_MF_time"][:]-(vars["nav_time"][1]*1e-9));
# quality estimate of acoustic measurement
nativeMeas[:,2] = vars["MF_std"]
# elevation' and azimuth -- TODO remove roll, pitch from body relative azimuth
nativeMeas[:,3:4] = vars["CBF_max"][:,1:2]
# argmax range
for i in 1:length(vars["CBF_MF_time"])
  nativeMeas[i,5] =findmax(sum(abs(vars["MF_prob"][:,:,i]),2))[2]
end
# nativeMeas[:,5] = findmax(sum(abs(vars["MF_prob"]),2),1)[2][:]
nativeMeas[:,5] -= 400.0
nativeMeas[:,5] *= 1500.0/37500

@parallel for i in  1:905
  pls = []
  for j in 1:4
    push!(pls, plot(y=(abs(vars["MF_prob"][:,j,i])),Geom.line,Coord.Cartesian(ymax=0.2)));
  end
  push!(pls, plot(y=sum(abs(vars["MF_prob"][:,:,i]),2),Geom.line,Coord.Cartesian(ymax=0.2)));
  pl = vstack(pls...)
  draw(PNG("/home/dehann/Desktop/here/MF_prob_$(i).png",10cm,35cm),pl);
end

# integrate native measurements for dead reckoning
# Dxi_k1 is previous k-1 position in pose i frame, x_{k-1}=xk1
function sandsharkDeadReckoning2D(dt, x, speed, pitch, dyaw)
  Dxi_k1 = SE2(x)
  v = speed*cos(pitch)
  D = SE2([0.0;0.0;dyaw])*SE2(dt*[v;0.0;0.0])
  Dxi_k = Dxi_k1*D
  return vec(se2vee(Dxi_k))
end

function deltaOdom(xprev::Array{Float64,1}, xnow::Array{Float64,1})
  xk1 = SE2(xprev)
  xk = SE2(xnow)
  return se2vee((xk1 \ xk))
end

# odo has form odo[1:n,:] = [t, x, y, theta], with n rows of data in time
function advOdoByRules(odo::Array{Float64,2}, meas::Array{Float64,2};
                        distrule=10.0, timerule=30.0, yawrule=pi/6.0, trkfeats=false)
  # DBG = Dict{Int,Dict{Int64,Array{BallTreeDensity,1}}}()
  # DBG = Dict{Int,Dict{Int64,Feature}}()
  DBG = nothing
  cntr = 1

  # bfts0 = Array{Float64,2}() #lsrFeats[1].feats
  instSys = makeGenericInSituSys(zeros(3))
  instSys.T0 = odo[1,1] # in nanoseconds
  fdict = Dict{Int64,Array{Float64,1}}()
  # for f in instSys.trackers
  #   # if length(f[2].lastz) < 3
  #   #   error("HEY!! $(f[1])")
  #   # end
  #   fdict[f[2].id] = f[2].lastz
  # end
  # instSys.FeatAssc[instSys.poseid] = fdict
  for i in 1:size(odo,1)
    dt = (odo[i,1]-instSys.T0)
    bTbm = SE2(instSys.x)
    # instSys.x is the real time odometry estimate
    instSys.x = sandsharkDeadReckoning2D(dt, instSys.x, odo[i,2], odo[i,4],odo[i,5])
    # bTbp = SE2(instSys.x)

    trigmeas = false
    if cntr <= size(meas,1)
      if odo[i,1] >= meas[cntr,1]
        # TODO -- relax to include at least range on large azimuth
        if meas[cntr,2] < 5.0 &&  (meas[cntr,4] < 95.0 || 265.0 < meas[cntr,4])
          trigmeas = true
        end
        cntr+=1
      end
    end
  #   # subsample motion tracks for constructing the factor graph SLAM solution
    if poseTrigAndAdd!(instSys, (odo[i,1]), distrule, timerule, yawrule, auxtrig=trigmeas)
      fdict = Dict{Int64,Array{Float64,1}}()
      if trigmeas
        bearing = wrapRad(meas[cntr-1,4]*pi/180.0)
        fdict[1] = [meas[cntr-1,5];bearing;0.0]
      end
      instSys.FeatAssc[instSys.poseid] = fdict
    end
  #   # advance previous time point
    instSys.T0 = odo[i,1]
  end
  return instSys.dOdo, instSys.FeatAssc, DBG
end

XX = advOdoByRules(nativeOdo, nativeMeas);
XX[2][1] = Dict{Int64,Array{Float64,1}}()

progressExamplePlot(XX[1],XX[2])

MM = Dict{Int, Int}()

fg = emptyFactorGraph();
idx = appendFactorGraph!(fg, XX[1], XX[2], lcmode=:unimodal, MM=MM);

tree = prepBatchTree!(fg, ordering=:qr,drawpdf=true);
inferOverTree!(fg,tree, N=100);

isamdict, fgu = doISAMSolve(XX[1],XX[2], savejld=false, MM=MM, retfg=true)
isamdict[1] = zeros(3)
# @save "odoSandsharkISAM1.jld" isamdict fgu
@load "odoSandsharkISAM1.jld"

@save "sandsharkSolve01.jld" fg nativeOdo nativeMeas

pl = drawCompPosesLandm(fg, isamdict, fgu, lbls=false)
using Colors, Cairo
draw(PDF("ss01.pdf",20cm,40cm),pl)

function testSandsharkOdo(;iters::Int=10)
  odo = zeros(iters, 5)
  dt = 0.1
  odo[:,1] = collect(1:iters)*dt
  odo[:,2] = 1.0
  odo[4:end,5] = pi/4
  odo[8:end,5] = pi/2.0

  prevpose = zeros(3)
  prevx = zeros(3)
  X = Array{Float64,2}(iters,3)
  for i in 1:iters
    x = sandsharkDeadReckoning2D(dt, prevx, odo[i,2], odo[i,4], odo[i,5])
    X[i,:] = deepcopy(x)
    prevx = x
    # do pose trigger test
      # new pose
  end

  return X
end

# X = testSandsharkOdoReal(nativeOdo,iters=10000);
# plot(x=X[:,1],y=X[:,2], Geom.path())
function testSandsharkOdoReal(odo;iters::Int=100)

  Tprev = odo[1,1]

  prevx = zeros(3)
  X = Array{Float64,2}(iters,3)
  for i in 1:iters
    dt = (odo[i,1] - Tprev)
    x = sandsharkDeadReckoning2D(dt, prevx, odo[i,2], odo[i,4], odo[i,5])
    X[i,:] = deepcopy(x)
    prevx = x
    # do pose trigger test
      # new pose
    Tprev = odo[i,1]
  end

  return X
end


text:
abstract type FiniteLengthAnimation{T} end

Base.Broadcast.broadcastable(f::FiniteLengthAnimation) = Ref(f)

"""
`Animation{T}`

An Animation that contains a `Vector` of `Keyframe`s and a `Vector` of `Easing`s
"""
struct Animation{T} <: FiniteLengthAnimation{T}
    frames::Vector{Keyframe{T}}
    easings::Vector{<:Easing}

    function Animation(kfs::Vector{Keyframe{T}}, easings::Vector{<:Easing}) where T

        validate_keyframe_times(kfs)

        if length(kfs) - length(easings) != 1
            error(
                """
                There must be one more keyframe than easings.
                There are $(length(kfs)) keyframes but only $(length(easings)) easings.
                """)
        end

        new{T}(kfs, easings)
    end
end

function Animation(differenttype_kfs::Vector{Keyframe}, easings::Vector{Easing})
    types = unique(map(x -> typeof(x), differenttype_kfs))
    error("""All keyframes need the same parametric type. Types are:\n$types""")
end

"""
`Animation(kfs::Vector{Keyframe{T}}, easing::Easing) where T`

Create an `Animation` from a `Vector` of `Keyframe`s and one `Easing` that is repeated
for every pair of keyframes.
"""
function Animation(kfs::Vector{Keyframe{T}}, easing::Easing) where T
    if length(kfs) <= 1
        error("There must be at least two keyframes.")
    end
    Animation(kfs, Easing[easing for _ in 1:length(kfs) - 1])
end

function Animation(timestamps::AbstractVector{<:Real}, values::AbstractVector{T}, easings::AbstractVector{<:Easing}) where T
    keyframes = Keyframe{T}.(timestamps, values)
    Animation(keyframes, easings)
end

function Animation(timestamps::AbstractVector{<:Real}, values::AbstractVector{T}, easing::Easing=linear()) where T
    Animation(timestamps, values, Easing[easing for _ in 1:(length(timestamps) - 1)])
end

function Animation(args...; defaulteasing=Easing())

    timestamps = Float64[]
    valtype = typeof(args[2])
    values = Vector{valtype}()
    easings = Easing[]

    push!(timestamps, args[1])
    last = :time
    i = 2
    while i <= length(args)
        if last == :time
            if typeof(args[i]) <: valtype
                push!(values, args[i])
            else
                error("Value with type <: $valtype expected after timestamp at index $i. Got $(typeof(args[i])) instead.")
            end
            last = :value
        elseif last == :value
            if typeof(args[i]) <: Easing
                push!(easings, args[i])
                last = :easing
            elseif typeof(args[i]) <: Real
                push!(easings, defaulteasing)
                push!(timestamps, args[i])
                last = :time
            elseif typeof(args[i]) <: Relative
                push!(easings, defaulteasing)
                push!(timestamps, timestamps[end] + args[i].t)
                last = :time
            else
                error("Timestamp with type <: Real or Relative, or easing with type <: Easing expected after value at index $i. Got $(typeof(args[i])) instead.")
            end
        elseif last == :easing
            if typeof(args[i]) <: Real
                push!(timestamps, args[i])
            elseif typeof(args[i]) <: Relative
                push!(timestamps, timestamps[end] + args[i].t)
            else
                error("Timestamp with type <: Real or Relative expected after easing at index $i. Got $(typeof(args[i])) instead.")
            end
            last = :time
        end

        i += 1
    end

    if last == :time
        error("Last value can't be a timestamp.")
    elseif last == :easing
        error("Last value can't be an easing")
    end

    Animation(timestamps, values, easings)
end

Base.:+(a::Animation, t::Real) = Animation(a.frames .+ t, a.easings)
Base.:-(a::Animation, t::Real) = Animation(a.frames .- t, a.easings)
Base.:*(a::Animation, stretch::Real) = Animation(a.frames .* stretch, a.easings)
Base.:/(a::Animation, compress::Real) = Animation(a.frames ./ compress, a.easings)

easings(a::Animation) = a.easings
timestamps(a::Animation) = [kf.t for kf in a.frames]
keyframes(a::Animation) = a.frames
keyvalues(a::Animation) = [kf.value for kf in a.frames]

function at(a::Animation, t::Real)
    # the first keyframe with a higher time is the second one of the two with
    # t in between (except when t is before the first or after the last keyframe)
    i_first_after_t = findfirst(kf -> kf.t >= t, a.frames)

    if isnothing(i_first_after_t)
        # t lies after the last keyframe
        # return a.frames[end].value
        return interpolate(a.easings[end], t, a.frames[end-1], a.frames[end])
    elseif i_first_after_t == 1
        # t lies before the first keyframe
        # return a.frames[1].value
        return interpolate(a.easings[1], t, a.frames[1], a.frames[2])
    else
        # t lies between two keyframes
        i_from = i_first_after_t - 1
        i_to = i_first_after_t
        return interpolate(a.easings[i_from], t, a.frames[i_from], a.frames[i_to])
    end
end

(a::Animation)(t::Real) = at(a, t)

struct Relative
    t::Float64

    function Relative(t::Real)
        if t <= 0
            error("a relative has to be positive, not $t")
        end
        new(convert(Float64, t))
    end
end

rel(t::Real) = Relative(t)

start(a::Animation) = a.frames[1].t
stop(a::Animation) = a.frames[end].t
duration(a::Animation) = stop(a) - start(a)


text:
""" Finite presentations of a model of a generalized algebraic theory (GAT).

We support two methods for defining models of a GAT: as Julia objects using the
`@instance` macro and as syntactic objects using the `@present` macro.
Instances are useful for casting generic data structures, such as matrices,
abstract tensor systems, and wiring diagrams, in categorical language.
Presentations define small categories by generators and relations and are useful
in applications like knowledge representation.
"""
module Present
export @present, Presentation, generator, generators, has_generator, equations,
  add_generator!, add_generators!, add_definition!, add_equation!, generator_index

using Base.Meta: ParseError
using MLStyle: @match

using ..Meta, ..Syntax
import ..GAT
import ..Syntax: parse_json_sexpr, to_json_sexpr

# Data types
############

struct Presentation{Theory,Name}
  syntax::Module
  generators::NamedTuple
  generator_name_index::Dict{Name,Pair{Symbol,Int}}
  equations::Vector{Pair}
end

function Presentation{Name}(syntax::Module) where Name
  Theory = syntax.theory()
  theory = GAT.theory(Theory)
  names = Tuple(type.name for type in theory.types)
  vectors = ((getfield(syntax, name){:generator})[] for name in names)
  Presentation{Theory,Name}(syntax, NamedTuple{names}(vectors),
                            Dict{Name,Pair{Symbol,Int}}(), Pair[])
end
Presentation(syntax::Module) = Presentation{Symbol}(syntax)

function Base.:(==)(pres1::Presentation, pres2::Presentation)
  pres1.syntax == pres2.syntax && pres1.generators == pres2.generators &&
    pres1.equations == pres2.equations
end

function Base.copy(pres::Presentation{T,Name}) where {T,Name}
  Presentation{T,Name}(pres.syntax, map(copy, pres.generators),
                       copy(pres.generator_name_index), copy(pres.equations))
end

# Presentation
##############

""" Get all generators of a presentation.
"""
generators(pres::Presentation) = collect(Iterators.flatten(pres.generators))
generators(pres::Presentation, type::Symbol) = pres.generators[type]
generators(pres::Presentation, type::Type) = generators(pres, nameof(type))

""" Retrieve generators by name.

Generators can also be retrieved using indexing notation, so that
`generator(pres, name)` and `pres[name]` are equivalent.
"""
function generator(pres::Presentation, name)
  type, index = pres.generator_name_index[name]
  pres.generators[type][index]
end
Base.getindex(pres::Presentation, name) = generator.(Ref(pres), name)

""" Does the presentation contain a generator with the given name?
"""
function has_generator(pres::Presentation, name)
  haskey(pres.generator_name_index, name)
end

""" Add a generator to a presentation.
"""
function add_generator!(pres::Presentation, expr)
  name, type = first(expr), gat_typeof(expr)
  generators = pres.generators[type]
  if !isnothing(name)
    if haskey(pres.generator_name_index, name)
      error("Name $name already defined in presentation")
    end
    pres.generator_name_index[name] = type => length(generators)+1
  end
  push!(generators, expr)
  expr
end

""" Add iterable of generators to a presentation.
"""
function add_generators!(pres::Presentation, exprs)
  for expr in exprs
    add_generator!(pres, expr)
  end
end

""" Get all equations of a presentation.
"""
equations(pres::Presentation) = pres.equations

""" Add an equation between terms to a presentation.
"""
function add_equation!(pres::Presentation, lhs::GATExpr, rhs::GATExpr)
  push!(pres.equations, lhs => rhs)
end

""" Add a generator defined by an equation.
"""
function add_definition!(pres::Presentation, name::Symbol, rhs::GATExpr)
  generator = Syntax.generator_like(rhs, name)
  add_generator!(pres, generator)
  add_equation!(pres, generator, rhs)
  generator
end

""" Get the index of a generator
"""
generator_index(pres::Presentation, x::Symbol) = pres.generator_name_index[x].second

# Presentation macro
####################

""" Define a presentation using a convenient syntax.
"""
macro present(head, body)
  name, pres = @match head begin
    Expr(:call, name::Symbol, syntax_name::Symbol) =>
      (name, :($(module_ref(:Presentation))($syntax_name)))
    Expr(:(<:), name::Symbol, parent::Symbol) => (name, :(copy($parent)))
    _ => throw(ParseError("Ill-formed presentation header $head"))
  end
  let_expr = Expr(:let, Expr(:block), translate_presentation(pres, body))
  esc(Expr(:(=), name, let_expr))
end

""" Translate a presentation in the DSL to a block of Julia code.
"""
function translate_presentation(pres::Expr0, body::Expr)::Expr
  @assert body.head == :block
  code = Expr(:block)
  append_expr!(code, :(_presentation = $pres))
  for expr in strip_lines(body).args
    append_expr!(code, translate_statement(expr))
  end
  append_expr!(code, :(_presentation))
  code
end

""" Translate a single statement in the presentation DSL to Julia code.
"""
function translate_statement(expr::Expr)::Expr
  @match expr begin
    Expr(:(::), name::Symbol, type_expr) =>
      translate_generator(name, type_expr)
    Expr(:(::), Expr(:tuple, names...), type_expr) =>
      Expr(:block, (translate_generator(name, type_expr) for name in names)...)
    Expr(:(::), type_expr) =>
      translate_generator(nothing, type_expr)
    Expr(:(:=), name::Symbol, def_expr) =>
      translate_definition(name, def_expr)
    Expr(:call, :(==), lhs, rhs) => translate_equation(lhs, rhs)
    _ => throw(ParseError("Ill-formed presentation statement $expr"))
  end
end

""" Translate declaration of a generator.
"""
function translate_generator(name::Union{Symbol,Nothing}, type_expr)::Expr
  type_name, args = @match type_expr begin
    name::Symbol => (name, [])
    Expr(:call, name::Symbol, args...) => (name, args)
    _ => throw(ParseError("Ill-formed type expression $type_expr"))
  end
  gen_expr = Expr(:call, type_name,
                  isnothing(name) ? nothing : QuoteNode(name), args...)
  Expr(:call, module_ref(:add_generator!), :_presentation,
       translate_expr(gen_expr))
end

""" Translate definition of generator in terms of other generators.
"""
function translate_definition(name::Symbol, def_expr)::Expr
  Expr(:call, module_ref(:add_definition!), :_presentation,
       QuoteNode(name), translate_expr(def_expr))
end

""" Translate declaration of equality between terms.
"""
function translate_equation(lhs, rhs)::Expr
  Expr(:call, module_ref(:add_equation!), :_presentation,
       translate_expr(lhs), translate_expr(rhs))
end

""" Translate GAT expression appearing in presentation.
"""
function translate_expr(expr)
  @match expr begin
    Expr(:call, name::Symbol, args...) =>
      Expr(:call, GlobalRef(Syntax, :invoke_term), :(_presentation.syntax),
           QuoteNode(name), map(translate_expr, args)...)
    name::Symbol =>
      Expr(:call, module_ref(:generator), :_presentation, QuoteNode(name))
    _ => expr
  end
end

module_ref(sym::Symbol) = GlobalRef(Present, sym)

# Serialization
###############

function to_json_sexpr(pres::Presentation, expr::GATExpr)
  to_json_sexpr(expr;
    by_reference = name -> has_generator(pres, name))
end

function parse_json_sexpr(pres::Presentation{Theory,Name},
                          syntax::Module, sexpr) where {Theory,Name}
  parse_json_sexpr(syntax, sexpr;
    symbols = Name == Symbol,
    parse_reference = name -> generator(pres, name))
end

end


text:
# This file was generated, do not modify it. # hide
ŷ = predict_mean(cmach, rows=test)
ŷ[1:3]


text:
function Base.getindex(d::XML, name::AbstractString)
    name = lowercase(name)
    for elem in d
        if isa(elem, XML) && lowercase(elem.name) == name
            return elem
        end
    end
end
function Base.getindex(d::XML, name::AbstractString, ::Colon)
    name = lowercase(name)
    found = []
    for elem in d
        if isa(elem, XML) && lowercase(elem.name) == name
            push!(found, elem)
        end
    end
    return found
end

Base.start(d::XML) = reverse(d.elems)

Base.done(::XML, s) = isempty(s)

function Base.next(d::XML, s)
    v = pop!(s)
    if isa(v, XML)
        append!(s, reverse(v.elems))
    end
    return v, s
end


text:
module LotMassing

using NPFinancial, JuMP, BlackBoxOptim, Polyhedra, Cbc #,GAMS

mutable struct RotInfo
    mat
    cr
    theta
end

mutable struct SubPoly
    points
    ladoComun
end

mutable struct PolyData
    V
    A
    b
    proyeccion
end

mutable struct ConstraintData
    V
    k
    A
    b
    norm
end

mutable struct PolyShape
    Vertices
    NumRegions
end

mutable struct FlagPlotEdif3D
    predio
    volTeor
    restSombra
    edif
    sombraVolTeor_p
    sombraVolTeor_o
    sombraVolTeor_s
    sombraEdif_p
    sombraEdif_o
    sombraEdif_s
end

mutable struct ResultadoCabida
    SalidaNormativa
    SalidaArquitectonica
    SalidaIndicadores
    SalidaTerreno
    SalidaMonetaria
    SalidaFlujoCaja
    Xopt
end

mutable struct datosCabidaPredio
    x
    y
    ladosConCalle
    ANCHOESPACIOPUBLICO
    FUSIONTERRENOS
    DISTANCIAMETRO
end

mutable struct datosCabidaNormativa
    SEPMIN
    ANTEJARDIN
    RASANTE
    ALTURAMAX
    MAXPISOS
    COEFOCUPACION
    SUBPREDIALMIN
    DENSIDADMAX
    FLAGDENSIDADBRUTA
    COEFCONSTRUCTIBILIDAD
    ESTACIONAMIENTOSPORVIV
    PORCADICESTACVISITAS
    SUPPORESTACIONAMIENTO
    ESTBICICLETAPOREST
    BICICLETASPOREST
    FLAGCAMBIOESTPORBICICLETA
    MAXSUBTE
    COEFOCUPACIONEST
    SEPESTMIN
    REDUCCIONESTPORDISTMETRO
end

mutable struct datosCabidaArquitectura
    ALTURAPISO
    PORCSUPCOMUN
    PORCTERRAZA
    ANCHOMAX
    MATCELDASCONF
    MATCONFHOR
end

mutable struct datosCabidaComercial
    SUPDEPTOUTIL
    MAXPORCTIPODEPTO
    PRECIOVENTA
    PRECIOVENTAEST
end

mutable struct datosCabidaUnit
    duracionProyecto
    CostoTerreno
    ComisionCorredor
    Demolicion
    LosaSNT
    LosaBNT
    EstacVisita
    ExtrasPostVenta
    Arquitectura
    GestionAdministracion
    ITO # por mes
    Contabilidad # por mes
    Legales
    Calculo
    MecanicaSuelo
    Topografia
    ProyectoElectrico
    ProyectoSanitario
    ProyectoSeguridad
    ProyectoPavimentacion
    ProyectoBasura
    ProyectoCalefaccion
    Marketing
    PilotSalaVentas
    NotarialesLegales
    CBR
    EmpalmesAportes
    DerechosPermisosMunicipales
    InspeccionCertificacionTecnica # por mes
    Contribuciones
    CopiaPlanosArquitctura
    CopiaPlanosObra
    HallAcceso
    JardinesPaisajismo
    HabilitacionSalaUsosMultiples
    HabilitacionPiscinaEnfermeria
    CuentasServiciosDeptos
    GastosComunes
    ContribucionesDeptos
    Seguros
    ResponsabilidadCivil
    PostVentaInmobiliaria
    SeguroVentaEnVerde
    Imprevistos
end

# flujos corresponden a meses: 0 6 12 18 24 30
mutable struct datosCabidaFlujo
    IngresosVentas
    CostoTerreno
    CostoConstruccion
    CostoHonorariosProyectos
    CostoVenta
    CostoInmobiliarioObra
    CostosHabilitacion
    CostosPuestaEnMarcha
    CostosAtencionCliente
    Imprevistos
    ingresoLineaCredito
    pagoLineaCredito
    TasaImpuestoRenta
    TasaInteresLineaCredito
    DuracionLineaCredito
    Tasaciones
end

# Rentabilidad exigida
mutable struct datosCabidaRentabilidad
    RetornoExigido
end

struct salidaNormativa
    maxNumDeptos
    maxOcupacion
    maxConstructibilidad
    maxPisos
    maxAltura
    minEstacionamientosVendibles
    minEstacionamientosVisita
    minEstacionamientosDiscapacitados
end

struct salidaArquitectonica
    numDeptosTipo
    numDeptos
    ocupacion
    constructibilidad
    numPisos
    altura
    superficieUtilSNT
    superficieComunSNT
    superficieEdificadaSNT
    superficiePorPiso
    estacionamientosVendibles
    estacionamientosVisita
    numEstacionamientos
    numBicicleteros
end

struct salidaTerreno
    superficieTerreno
    superficieBruta
    superficieEdificacion
    costoTerrenoAntesCorr
    costoUnitTerrenoAntesCorr
    CostoCorredor
    costoTerreno
    costoUnitTerreno
end

struct salidaIndicadores
    IngresosVentas
    CostoTotal
    MargenAntesImpuesto
    TirAntesImpuestos
    ImpuestoRenta
    UtilidadDespuesImpuesto
    RentabilidadTotalBruta
    RentabilidadTotalNeta
    IncidenciaTerreno
end

struct salidaMonetaria
    IngresosVentas
    CostoTotal
    CostoTerreno
    CostoUnitarioTerreno
    CostoConstruccion
    CostoHonorariosProyectos
    CostoVenta
    CostoInmobiliarioObra
    CostosHabilitacion
    CostosPuestaEnMarcha
    CostosAtencionCliente
    Imprevistos
    CostoPagoIVA
end

struct salidaFlujoCaja
    IngresosVentas
    CostoTerreno
    CostoConstruccion
    CostoHonorariosProyectos
    CostoVenta
    CostoInmobiliarioObra
    CostosHabilitacion
    CostosPuestaEnMarcha
    CostosAtencionCliente
    Imprevistos
    CostoPagoIVA
    flujoCajaNetoAntesImpuestos
    TirAntesImpuestos
    ImpuestoRenta
    NetoDespuesImpuestos
    TirDespuesImpuestos
    LineaCredito
    MontoUtilizadoLineaCredito
    InteresLineaCredito
    CostoOperacionalCredito
    ApalancadoNetoAntesImpuesto
    TirApalancadoAntesImpuestos
end

struct salidaGeometria
    config
    celdasConfig
    p1
    p2
    p3
    p4
    pos_x
    pos_y
    largo
end

export datosCabidaPredio, datosCabidaNormativa, datosCabidaArquitectura, datosCabidaComercial, datosCabidaUnit,
         datosCabidaFlujo, datosCabidaRentabilidad, RotInfo, SubPoly, PolyData, salidaArquitectonica, salidaIndicadores, salidaMonetaria,
         salidaTerreno, salidaNormativa, salidaFlujoCaja, salidaGeometria, PolyShape, FlagPlotEdif3D, ConstraintData,
         ResultadoCabida

struct particulaPso
    pos_x
    pos_y
    largo_x
    largo_y
    theta
    alt
end

struct polyPso
    poly
    Rmat
    cr
    p_1
    p_2
    p_3
    p_4
end

export particulaPso, polyPso

MATCELDASCONF_ = [[1 3 5];
                  [2 3 4]; #
                  [1 3 NaN];
                  [3 5 NaN]; #
                  [2 3 NaN]; #
                  [3 4 NaN]; #
                  [1 NaN NaN];
                  [2 NaN NaN]]; #
MATCONFHOR_ = [[1 0 1];
               [0 1 0];
               [1 0 NaN];
               [0 1 NaN];
               [0 1 NaN];
               [0 1 NaN];
               [1 NaN NaN];
               [0 NaN NaN]];

"""
MATCELDASCONF_ = [[1 3 5];
                  [1 3 NaN];
                  [1 NaN NaN]];
MATCONFHOR_ = [[1 0 1];
               [1 0 NaN];
               [1 NaN NaN]];
"""

export MATCELDASCONF_, MATCONFHOR_

include("generaMatCeldasConf.jl")
include("infoPredio.jl")
include("calculaAnguloRotacion.jl")
include("nonconv2sumofconv_v2.jl")
include("ismembern.jl")
include("pso.jl")
include("plotCabidaOptima.jl")
include("generaCalles.jl")
include("generaSombraEdificio.jl")
include("optiEdificio.jl")
include("displayResults.jl")
include("plotBaseEdificio3d_v2.jl")
include("poly2D.jl")
include("polyShape.jl")
include("evol.jl")
include("resultConverter_v2.jl")
include("executaCalculoCabidas.jl")
include("generaVol3d_v4.jl")
include("generaSombraTeor_v3.jl")
include("randomPointsNearVertices.jl")
include("ajusteArea.jl")
include("generaSupBruta.jl")
include("combinaSoluciones_v2.jl")

export generaMatCeldasConf, infoPredio, calculaAnguloRotacion, ismembern, plotCabidaOptima,
       pso, generaCalles, generaSombraEdificio,
       optiEdificio,  displayResults, evol, poly2D, polyShape,
       nonconv2sumofconv_v2, resultConverter_v2, plotBaseEdificio3d_v2,
       executaCalculoCabidas, generaVol3d_v4, generaSombraTeor_v3, randomPointsNearVertices, ajusteArea, generaSupBruta,
       combinaSoluciones_v2
end


text:
function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DP8ConstantCache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<7 || calcVal
    @unpack c7,c8,c9,c10,c11,c6,c5,c4,c3,c2,b1,b6,b7,b8,b9,b10,b11,b12,bhh1,bhh2,bhh3,er1,er6,er7,er8,er9,er10,er11,er12,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache
    @unpack c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615 = cache
    @unpack d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 = cache
    k1 = f(t,uprev)
    k2 = f(t+c2*dt,uprev+dt*(a0201*k1))
    k3 = f(t+c3*dt,uprev+dt*(a0301*k1+a0302*k2))
    k4 = f(t+c4*dt,uprev+dt*(a0401*k1       +a0403*k3))
    k5 = f(t+c5*dt,uprev+dt*(a0501*k1       +a0503*k3+a0504*k4))
    k6 = f(t+c6*dt,uprev+dt*(a0601*k1                +a0604*k4+a0605*k5))
    k7 = f(t+c7*dt,uprev+dt*(a0701*k1                +a0704*k4+a0705*k5+a0706*k6))
    k8 = f(t+c8*dt,uprev+dt*(a0801*k1                +a0804*k4+a0805*k5+a0806*k6+a0807*k7))
    k9 = f(t+c9*dt,uprev+dt*(a0901*k1                +a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8))
    k10= f(t+c10*dt,uprev+dt*(a1001*k1                +a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9))
    k11= f(t+c11*dt,uprev+dt*(a1101*k1                +a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10))
    k12= f(t+dt,uprev+dt*(a1201*k1                +a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11))
    kupdate= b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12
    update = dt*kupdate
    utmp = uprev + update
    k13 = f(t+dt,utmp)
    k14 = f(t+c14*dt,uprev+dt*(a1401*k1         +a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13))
    k15 = f(t+c15*dt,uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8                   +a1511*k11+a1512*k12+a1513*k13+a1514*k14))
    k16 = f(t+c16*dt,uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9                              +a1613*k13+a1614*k14+a1615*k15))
    udiff = kupdate
    copyat_or_push!(k,1,udiff)
    bspl = k1 - udiff
    copyat_or_push!(k,2,bspl)
    copyat_or_push!(k,3,udiff - k13 - bspl)
    copyat_or_push!(k,4,(d401*k1+d406*k6+d407*k7+d408*k8+d409*k9+d410*k10+d411*k11+d412*k12+d413*k13+d414*k14+d415*k15+d416*k16))
    copyat_or_push!(k,5,(d501*k1+d506*k6+d507*k7+d508*k8+d509*k9+d510*k10+d511*k11+d512*k12+d513*k13+d514*k14+d515*k15+d516*k16))
    copyat_or_push!(k,6,(d601*k1+d606*k6+d607*k7+d608*k8+d609*k9+d610*k10+d611*k11+d612*k12+d613*k13+d614*k14+d615*k15+d616*k16))
    copyat_or_push!(k,7,(d701*k1+d706*k6+d707*k7+d708*k8+d709*k9+d710*k10+d711*k11+d712*k12+d713*k13+d714*k14+d715*k15+d716*k16))
  end
end

function ode_addsteps!{calcVal,calcVal2,calcVal3}(k,t,uprev,u,dt,f,cache::DP8Cache,always_calc_begin::Type{Val{calcVal}} = Val{false},allow_calc_end::Type{Val{calcVal2}} = Val{true},force_calc_end::Type{Val{calcVal3}} = Val{false})
  if length(k)<7 || calcVal
    @unpack c7,c8,c9,c10,c11,c6,c5,c4,c3,c2,b1,b6,b7,b8,b9,b10,b11,b12,bhh1,bhh2,bhh3,er1,er6,er7,er8,er9,er10,er11,er12,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab
    @unpack c14,c15,c16,a1401,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1613,a1614,a1615 = cache.tab
    @unpack d401,d406,d407,d408,d409,d410,d411,d412,d413,d414,d415,d416,d501,d506,d507,d508,d509,d510,d511,d512,d513,d514,d515,d516,d601,d606,d607,d608,d609,d610,d611,d612,d613,d614,d615,d616,d701,d706,d707,d708,d709,d710,d711,d712,d713,d714,d715,d716 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,update,udiff,bspl,dense_tmp3,dense_tmp4,dense_tmp5,dense_tmp6,dense_tmp7,kupdate,utilde,tmp,atmp,atmp2 = cache
    utmp = utilde
    k = [cache.udiff,cache.bspl,cache.dense_tmp3,cache.dense_tmp4,cache.dense_tmp5,cache.dense_tmp6,cache.dense_tmp7]
    uidx = eachindex(u)
    f(t,uprev,k1)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0201*k1[i])
    end
    f(t+c2*dt,tmp,k2)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0301*k1[i]+a0302*k2[i])
    end
    f(t+c3*dt,tmp,k3)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0401*k1[i]+a0403*k3[i])
    end
    f(t+c4*dt,tmp,k4)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0501*k1[i]+a0503*k3[i]+a0504*k4[i])
    end
    f(t+c5*dt,tmp,k5)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0601*k1[i]+a0604*k4[i]+a0605*k5[i])
    end
    f(t+c6*dt,tmp,k6)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0701*k1[i]+a0704*k4[i]+a0705*k5[i]+a0706*k6[i])
    end
    f(t+c7*dt,tmp,k7)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0801*k1[i]+a0804*k4[i]+a0805*k5[i]+a0806*k6[i]+a0807*k7[i])
    end
    f(t+c8*dt,tmp,k8)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a0901*k1[i]+a0904*k4[i]+a0905*k5[i]+a0906*k6[i]+a0907*k7[i]+a0908*k8[i])
    end
    f(t+c9*dt,tmp,k9)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a1001*k1[i]+a1004*k4[i]+a1005*k5[i]+a1006*k6[i]+a1007*k7[i]+a1008*k8[i]+a1009*k9[i])
    end
    f(t+c10*dt,tmp,k10)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a1101*k1[i]+a1104*k4[i]+a1105*k5[i]+a1106*k6[i]+a1107*k7[i]+a1108*k8[i]+a1109*k9[i]+a1110*k10[i])
    end
    f(t+c11*dt,tmp,k11)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a1201*k1[i]+a1204*k4[i]+a1205*k5[i]+a1206*k6[i]+a1207*k7[i]+a1208*k8[i]+a1209*k9[i]+a1210*k10[i]+a1211*k11[i])
    end
    f(t+dt,tmp,k12)
    for i in uidx
      kupdate[i] = b1*k1[i]+b6*k6[i]+b7*k7[i]+b8*k8[i]+b9*k9[i]+b10*k10[i]+b11*k11[i]+b12*k12[i]
      update[i] = dt*kupdate[i]
      utmp[i] = uprev[i] + update[i]
    end
    f(t+dt,utmp,k13)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a1401*k1[i]+a1407*k7[i]+a1408*k8[i]+a1409*k9[i]+a1410*k10[i]+a1411*k11[i]+a1412*k12[i]+a1413*k13[i])
    end
    f(t+c14*dt,tmp,k14)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a1501*k1[i]+a1506*k6[i]+a1507*k7[i]+a1508*k8[i]+a1511*k11[i]+a1512*k12[i]+a1513*k13[i]+a1514*k14[i])
    end
    f(t+c15*dt,tmp,k15)
    for i in uidx
      tmp[i] = uprev[i]+dt*(a1601*k1[i]+a1606*k6[i]+a1607*k7[i]+a1608*k8[i]+a1609*k9[i]+a1613*k13[i]+a1614*k14[i]+a1615*k15[i])
    end
    f(t+c16*dt,tmp,k16)
    for i in uidx
      udiff[i]= kupdate[i]
      bspl[i] = k1[i] - udiff[i]
      k[3][i] = udiff[i] - k13[i] - bspl[i]
      k[4][i] = (d401*k1[i]+d406*k6[i]+d407*k7[i]+d408*k8[i]+d409*k9[i]+d410*k10[i]+d411*k11[i]+d412*k12[i]+d413*k13[i]+d414*k14[i]+d415*k15[i]+d416*k16[i])
      k[5][i] = (d501*k1[i]+d506*k6[i]+d507*k7[i]+d508*k8[i]+d509*k9[i]+d510*k10[i]+d511*k11[i]+d512*k12[i]+d513*k13[i]+d514*k14[i]+d515*k15[i]+d516*k16[i])
      k[6][i] = (d601*k1[i]+d606*k6[i]+d607*k7[i]+d608*k8[i]+d609*k9[i]+d610*k10[i]+d611*k11[i]+d612*k12[i]+d613*k13[i]+d614*k14[i]+d615*k15[i]+d616*k16[i])
      k[7][i] = (d701*k1[i]+d706*k6[i]+d707*k7[i]+d708*k8[i]+d709*k9[i]+d710*k10[i]+d711*k11[i]+d712*k12[i]+d713*k13[i]+d714*k14[i]+d715*k15[i]+d716*k16[i])
    end
  end
end


text:
using ResumableFunctions
using SimJulia
using Distributions
using DataFrames
using Random
using StatsPlots
using BenchmarkTools

function increment!(a::Array{Int64})
    push!(a,a[length(a)]+1)
end

function decrement!(a::Array{Int64})
    push!(a,a[length(a)]-1)
end

function carryover!(a::Array{Int64})
    push!(a,a[length(a)])
end;

mutable struct SIRPerson
    id::Int64 # numeric ID
    status::Symbol # :S,I,R
end;

mutable struct SIRModel
    sim::Simulation
    β::Float64
    c::Float64
    γ::Float64
    ta::Array{Float64}
    Sa::Array{Int64}
    Ia::Array{Int64}
    Ra::Array{Int64}
    allIndividuals::Array{SIRPerson}
end;

function infection_update!(sim::Simulation,m::SIRModel)
    push!(m.ta,now(sim))
    decrement!(m.Sa)
    increment!(m.Ia)
    carryover!(m.Ra)
end;

function recovery_update!(sim::Simulation,m::SIRModel)
    push!(m.ta,now(sim))
    carryover!(m.Sa)
    decrement!(m.Ia)
    increment!(m.Ra)
end;

@resumable function live(sim::Simulation, individual::SIRPerson, m::SIRModel)
  while individual.status==:S
      # Wait until next contact
      @yield timeout(sim,rand(Distributions.Exponential(1/m.c)))
      # Choose random alter
      alter=individual
      while alter==individual
          N=length(m.allIndividuals)
          index=rand(Distributions.DiscreteUniform(1,N))
          alter=m.allIndividuals[index]
      end
      # If alter is infected
      if alter.status==:I
          infect = rand(Distributions.Uniform(0,1))
          if infect < m.β
              individual.status=:I
              infection_update!(sim,m)
          end
      end
  end
  if individual.status==:I
      # Wait until recovery
      @yield timeout(sim,rand(Distributions.Exponential(1/m.γ)))
      individual.status=:R
      recovery_update!(sim,m)
  end
end;

function MakeSIRModel(u0,p)
    (S,I,R) = u0
    N = S+I+R
    (β,c,γ) = p
    sim = Simulation()
    allIndividuals=Array{SIRPerson,1}(undef,N)
    for i in 1:S
        p=SIRPerson(i,:S)
        allIndividuals[i]=p
    end
    for i in (S+1):(S+I)
        p=SIRPerson(i,:I)
        allIndividuals[i]=p
    end
    for i  in (S+I+1):N
        p=SIRPerson(i,:R)
        allIndividuals[i]=p
    end
    ta=Array{Float64,1}(undef,0)
    push!(ta,0.0)
    Sa=Array{Int64,1}(undef,0)
    push!(Sa,S)
    Ia=Array{Int64,1}(undef,0)
    push!(Ia,I)
    Ra=Array{Int64,1}(undef,0)
    push!(Ra,R)
    SIRModel(sim,β,c,γ,ta,Sa,Ia,Ra,allIndividuals)
end;

function activate(m::SIRModel)
     [@process live(m.sim,individual,m) for individual in m.allIndividuals]
end;

function sir_run(m::SIRModel,tf::Float64)
    SimJulia.run(m.sim,tf)
end;

function out(m::SIRModel)
    result = DataFrame()
    result[!,:t] = m.ta
    result[!,:S] = m.Sa
    result[!,:I] = m.Ia
    result[!,:R] = m.Ra
    result
end;

tmax = 40.0;

u0 = [990,10,0];

p = [0.05,10.0,0.25];

Random.seed!(1234);

des_model = MakeSIRModel(u0,p)
activate(des_model)
sir_run(des_model,tmax)

data_des=out(des_model);

@df data_des plot(:t,
    [:S :I :R],
    labels = ["S" "I" "R"],
    xlab="Time",
    ylab="Number")

@resumable function live_delay(sim::Simulation, individual::SIRPerson, m::SIRModel)
  while individual.status==:S
      # Wait until next contact
      @yield timeout(sim,rand(Distributions.Exponential(1/m.c)))
      # Choose random alter
      alter=individual
      while alter==individual
          N=length(m.allIndividuals)
          index=rand(Distributions.DiscreteUniform(1,N))
          alter=m.allIndividuals[index]
      end
      # If alter is infected
      if alter.status==:I
          infect = rand(Distributions.Uniform(0,1))
          if infect < m.β
              individual.status=:I
              infection_update!(sim,m)
          end
      end
  end
  # The following is modified to capture a fixed infection delay
  if individual.status==:I
      # Wait until recovery, here fixed
      @yield timeout(sim,1/m.γ)
      individual.status=:R
      recovery_update!(sim,m)
  end
end;

function activate_delay(m::SIRModel)
     [@process live_delay(m.sim,individual,m) for individual in m.allIndividuals]
end;

des_delay_model = MakeSIRModel(u0,p)
activate_delay(des_delay_model)
sir_run(des_delay_model,tmax)

data_delay_des=out(des_delay_model);

@df data_delay_des plot(:t,
    [:S :I :R],
    labels = ["S" "I" "R"],
    xlab="Time",
    ylab="Number")

@benchmark begin
    des_model = MakeSIRModel(u0,p)
    activate(des_model)
    sir_run(des_model,tmax)
end