text:
"""
    plot(net::HybridNetwork, method::Symbol)

Plot a network using R graphics.
`method` should be `:R` (actually, any symbol would do, for now!).

optional arguments, shared with the Gadfly-based plot function:
- `useEdgeLength = false` : if true, the tree edges and major hybrid edges are
  drawn proportionally to their length. Minor hybrid edges are not, however.
  Note that edge lengths in coalescent units may scale very poorly with time.
- `showTipLabel = true` : if true, taxon labels are shown. You may need to zoom out to see them.
- `showNodeNumber = false` : if true, nodes are labelled with the number used internally.
- `showEdgeLength = false` : if true, edges are labelled with their length (above)
- `showGamma = false` : if true, hybrid edges are labelled with their heritability (below)
- `edgeColor = "black"` : color for tree edges.
- `majorHybridEdgeColor = "deepskyblue4"` : color for major hybrid edges
- `minorHybridEdgeColor = "deepskyblue"` : color for minor hybrid edges
- `showEdgeNumber = false` : if true, edges are labelled with the number used internally.
- `showIntNodeLabel = false` : if true, internal nodes are labelled with their names.
  Useful for hybrid nodes, which do have tags like 'H1'.
- `edgeLabel = DataFrame()` : dataframe with two columns: the first with edge numbers, the second with labels
  (like bootstrap values) to annotate edges. empty by default.
- `nodeLabel = DataFrame()` : dataframe with two columns: the first with node numbers, the second with labels
  (like bootstrap values for hybrid relationships) to annotate nodes. empty by default.
- `style = :fulltree` : symbol indicating the style of the diagram
  * `:majortree` will simply draw minor edges onto the major tree.
  * `:fulltree` will draw minor edges as their own branches in the tree (like in icytree.org),
    useful for overlapping or confusing networks.
- `arrowlen` : the length of the arrow tips in the full tree style. if `style = :fulltree`, then
  `arrowlen = 0.2`. otherwise, `arrowlen = 0`, which makes the arrows appear as segments.

optional arguments specific to this function:
- `xlim`, `ylim` : array of 2 values
- `tipOffset = 0.0` : to offset tip labels

plot() returns the following tuple:
`(xmin, xmax, ymin, ymax, node_x, node_y, node_yB, node_yE,
edge_xB, edge_xE, edge_yB, edge_yE, ndf, edf)`

1. `xmin` : the minimum x value of the plot
2. `xmax` : the maximum x value of the plot
3. `ymin` : the minimum y value of the plot
4. `ymax` : the maximum y value of the plot
5. `node_x` : the x values of the nodes in net.node in their respective order
6. `node_y` : the y values of the nodes
7. `node_yB` : the y value of the beginning of the verticle bar
8. `node_yE` : the y value of the end of the verticle bar
9. `node_yE` : the y value of the end of the verticle bar
10. `edge_xB` : the x value of the beginning of the edges in net.edge in their respective order
11. `edge_xE` : the x value of the end of the edges
12. `edge_yB` : the y value of the beginning of the edges
13. `edge_yE` : the y value of the end of the edges
14. `ndf` : the node data frame: see section [Adding labels](@ref) for more
15. `edf` : the edge data frame

Note that `plot` actually modifies some (minor) attributes of the network,
as it calls `directEdges!` and `preorder!`.

If hybrid edges cross tree and major edges, you may choose to rotate some tree
edges to eliminate crossing edges, using `rotate!`
(in [`PhyloNetworks`](http://crsl4.github.io/PhyloNetworks.jl/latest/lib/public/#PhyloNetworks.rotate!)).

**Alternative**: a tree or network can be exported with [`sexp`](@ref)
and then displayed with R's "plot" and all its options.
"""
function plot(net::HybridNetwork, ::Symbol; useEdgeLength=false::Bool,
    mainTree=false::Bool, showTipLabel=true::Bool, showNodeNumber=false::Bool,
    showEdgeLength=false::Bool, showGamma=false::Bool,
    edgeColor="black"::String,
    majorHybridEdgeColor="deepskyblue4"::String,
    minorHybridEdgeColor="deepskyblue"::String,
    showEdgeNumber=false::Bool, showIntNodeLabel=false::Bool,
    edgeLabel=DataFrame()::DataFrame, nodeLabel=DataFrame()::DataFrame,
    xlim=Float64[]::Array{Float64,1}, ylim=Float64[]::Array{Float64,1},
    tipOffset=0.0::Float64, tipcex=1.0::Float64,
    style=:fulltree::Symbol, arrowlen=(style==:majortree ? 0 : 0.1)::Real)

    (edge_xB, edge_xE, edge_yB, edge_yE, node_x, node_y, node_yB, node_yE,
     hybridedge_xB, hybridedge_xE, hybridedge_yB, hybridedge_yE,
     xmin, xmax, ymin, ymax) = getEdgeNodeCoordinates(net, useEdgeLength, style==:majortree)
    labelnodes, nodeLabel = checkNodeDataFrame(net, nodeLabel)
    ndf = prepareNodeDataFrame(net, nodeLabel, showNodeNumber,
            showIntNodeLabel, labelnodes, node_x, node_y)
    if (showTipLabel || showNodeNumber || showIntNodeLabel || labelnodes)
        expfac = 0.1  # force 10% more space to show tip/node/root name
        expfacy = 0.5 # additive expansion for y axis
        xmin -= (xmax-xmin)*expfac
        xmax += (xmax-xmin)*expfac
        ymin -= expfacy
        ymax += expfacy
    end
    xmax += tipOffset
    if length(xlim)==2
        xmin=xlim[1]; xmax=xlim[2]
    end
    if length(ylim)==2
        ymin=ylim[1]; ymax=ylim[2]
    end
    leaves = [n.leaf for n in net.node]
    eCol = fill(edgeColor, length(net.edge))
    eCol[ [ e.hybrid  for e in net.edge] ] .= majorHybridEdgeColor
    eCol[ [!e.isMajor for e in net.edge] ] .= minorHybridEdgeColor

     # this makes the arrows dashed if :fulltree is used
     arrowstyle = style==:majortree ? "solid" : "longdash"

    if !(style in [:fulltree, :majortree])
      @warn "Style $style is unknown. Defaulted to :fulltree."
      style = :fulltree
    end

    R"""
    plot($(node_x[leaves]), $(node_y[leaves]), type='n',
         xlim=c($xmin,$xmax), ylim=c($ymin,$ymax),
         axes=FALSE, xlab='', ylab='')
    segments($edge_xB, $edge_yB, $edge_xE, $edge_yE, col=$eCol)
    arrows($hybridedge_xB, $hybridedge_yB, $hybridedge_xE, $hybridedge_yE, length=$arrowlen, angle = 20, col=$minorHybridEdgeColor, lty=$arrowstyle)
    segments($node_x, $node_yB, $node_x, $node_yE, col=$edgeColor,)
    """
    if showTipLabel
      R"text"(node_x[leaves] .+ tipOffset, node_y[leaves],
              tipLabels(net), adj=0, font=3, cex=tipcex)
    end
    if showIntNodeLabel
      R"text"(ndf[.!ndf[!,:lea],:x], ndf[.!ndf[!,:lea],:y],
              ndf[.!ndf[!,:lea],:name], font=3, cex=tipcex, adj=[.5,0])
    end
    if showNodeNumber
      R"text"(ndf[!,:x], ndf[!,:y], ndf[!,:num], adj=1)
    end
    if labelnodes
      R"text"(ndf[!,:x], ndf[!,:y], ndf[!,:lab], adj=1)
    end
    labeledges, edf = prepareEdgeDataFrame(net, edgeLabel, mainTree,
                        edge_xB, edge_xE, edge_yB, edge_yE,
                        hybridedge_xB, hybridedge_xE, hybridedge_yB, hybridedge_yE)
    if labeledges
      R"text"(edf[!,:x], edf[!,:y], edf[!,:lab], adj=[.5,0])
    end
    if showEdgeLength
      R"text"(edf[!,:x], edf[!,:y], edf[!,:len], adj=[.5,1.])
    end
    if (showGamma && net.numHybrids>0)
      im = edf[!,:hyb] .& edf[!,:min]
      iM = edf[!,:hyb] .& .!edf[!,:min]
      R"text"(edf[im,:x], edf[im,:y], edf[im,:gam],
              adj=[.5,1], col=minorHybridEdgeColor)
      R"text"(edf[iM,:x], edf[iM,:y], edf[iM,:gam],
              adj=[.5,1], col=majorHybridEdgeColor)
    end
    if showEdgeNumber
      R"text"(edf[!,:x], edf[!,:y], edf[!,:num], adj=[.5,0])
    end
    return (xmin, xmax, ymin, ymax, node_x, node_y, node_yB, node_yE,
      edge_xB, edge_xE, edge_yB, edge_yE, ndf, edf)
end

max_stars_repo_name: cecileane/PhyloPlots

max_issues_repo_name: cecileane/PhyloPlots

max_forks_repo_name: cecileane/PhyloPlots


text:
using NiLang, NiLang.AD
using Test, Random
using ForwardDiff
using FixedPointNumbers
using NiLangCore: default_constructor

@testset "ULogarithmic" begin
	@test check_grad(PlusEq(gaussian_log), (1.0, 2.0); iloss=1)
	function muleq(f, x, y, z)
		x = default_constructor(ULogarithmic, x)
		y = default_constructor(ULogarithmic, y)
		z = default_constructor(ULogarithmic, z)
		x *= f(y, z)
		x.log
	end
	g1 = ForwardDiff.gradient(arr->muleq(+, arr...), [7.0, 5.0, 3.0])
	x, y, z = default_constructor(ULogarithmic, 7.0),
		default_constructor(ULogarithmic, 5.0),
		default_constructor(ULogarithmic, 3.0)
	@instr (MulEq(+))(x, y, z)
	@instr GVar(x)
	@instr GVar(y)
	@instr GVar(z)
	@instr x.log.g += 1
	@instr (~MulEq(+))(x, y, z)
	@test grad(x.log) ≈ g1[1]
	@test grad(y.log) ≈ g1[2]
	@test grad(z.log) ≈ g1[3]

	g2 = ForwardDiff.gradient(arr->muleq(-, arr...), [7.0, 5.0, 3.0])
	x, y, z = default_constructor(ULogarithmic, 2.0),
		default_constructor(ULogarithmic, 5.0),
		default_constructor(ULogarithmic, 3.0)
	@instr (MulEq(-))(x, y, z)
	@instr GVar(x)
	@instr GVar(y)
	@instr GVar(z)
	@instr x.log.g += 1
	@instr (~MulEq(-))(x, y, z)
	@test grad(x.log) ≈ g2[1]
	@test grad(y.log) ≈ g2[2]
	@test grad(z.log) ≈ g2[3]
end

@testset "iexp" begin
	@i function i_exp(y!::T, x::T) where T<:Union{Fixed, GVar{<:Fixed}}
	    @invcheckoff begin
	        @routine begin
	            s ← one(ULogarithmic{T})
	            lx ← one(ULogarithmic{T})
	            k ← 0
	        end
	        lx *= convert(x)
	        y! += convert(s)
	        while (s.log > -20, k != 0)
	            k += 1
	            s *= lx / k
	            y! += convert(s)
	        end
	        ~(while (s.log > -20, k != 0)
	            k += 1
	            s *= x / k
	        end)
	        lx /= convert(x)
	        ~@routine
	    end
	end

	x = Fixed43(3.5)
	res = i_exp(Fixed43(0.0), x)[1]
	gx = grad(Grad(i_exp)(Val(1), Fixed43(0.0), x)[3])
	@test res ≈ exp(3.5)
	@test gx ≈ exp(3.5)
end

ext: jl


text:
module KitFort

export flux_kfvs!,
       flux_kcu!,
       flux_ugks!,
       gfortran,
       ifort

include("kfvs.jl")
include("kcu.jl")
include("ugks.jl")

COMPILER = :gfortran

gfortran(file::AbstractString) = "__kinetic_MOD_" * file

ifort(file::AbstractString) = "kinetic_mp_" * file * "_"

end # module

max_stars_repo_name: vavrines/Kifort

max_issues_repo_name: vavrines/Kifort

max_forks_repo_name: vavrines/Kifort


text:
using Documenter, MajorizationExtrema

makedocs(;
    modules=[MajorizationExtrema],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        # "Subentropy counterexample" => "subentropy.md",
    ],
    repo="https://github.com/ericphanson/MajorizationExtrema.jl/blob/{commit}{path}#L{line}",
    sitename="MajorizationExtrema.jl",
    authors="Eric",
    assets=String[],
)

deploydocs(;
    repo="github.com/ericphanson/MajorizationExtrema.jl",
)


text:
cd("C:\\Users\\victo\\source\\repos\\Silly_Julia_stuff\\Julia Packages\\LassoPath\\usage_examples\\crime_data")
using DelimitedFiles, CSV
using DataFrames

D = readdlm("communities.csv", ',')
Header = readdlm("header_info.txt", ' ')
Df = DataFrame(D, :auto)
rename!(Df, Symbol.(reshape(Header[:, 1], :)))

# convert the types
function FilterTheDataFrame(Df::DataFrame)
    GoodColumns = []
    for (II, Col) in enumerate(eachcol(Df))
        if !("?" ∈ Df[:, II])
            push!(GoodColumns, II)
        end
        convertedType = Header[II, 2] == "numeric" ? Float64 : String
        convertFunc = entry -> begin
            if entry == "?"
                return NaN64
            else
                return convert(convertedType, entry)
            end

        end
        Df[!, II] = convertFunc.(Df[:, II])  # To include type, `!` will be needed.
    end
    return Df[!, GoodColumns]
end

Df = FilterTheDataFrame(Df)
Df = Df[!, Not(2)]
describe(Df)

# Save the processed DataFrame
CSV.write("process_data.csv", Df)


text:
"""
JuMP constraints
"""
sense_to_set(::Function, ::Val{:!=}) = CPE.DifferentFrom(0.0)
sense_to_set(::Function, ::Val{:<}) = CPE.Strictly(MOI.LessThan(0.0))
sense_to_set(::Function, ::Val{:>}) = CPE.Strictly(MOI.GreaterThan(0.0))

include("indicator.jl")
include("reified.jl")
include("bool.jl")
include("complement.jl")

"""
MOI constraints
"""

"""
Linear constraints
"""
MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{MOI.EqualTo{T}},
) where {T<:Real} = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{MOI.LessThan{T}},
) where {T<:Real} = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{CPE.AllEqual},
) = true
MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{CPE.AllDifferent},
) = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{VAF{T}},
    ::Type{CPE.AllDifferent},
) where {T <: Real} = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{TableSetInternal},
) = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{VAF{T}},
    ::Type{TableSetInternal},
) where {T <: Real} = true

MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{CPE.DifferentFrom{T}},
) where {T<:Real} = true

# Don't directly support StrictlyGreaterThan => use a bridge
MOI.supports_constraint(
    ::Optimizer,
    ::Type{SAF{T}},
    ::Type{CPE.Strictly{MOI.LessThan{T}, T}},
) where {T<:Real} = true

function MOI.supports_constraint(
    optimizer::Optimizer,
    fct::Type{<:MathOptInterface.AbstractFunction},
    ::Type{<:ComplementSet{F,S}},
) where {F,S}
    return MOI.supports_constraint(optimizer, F, S)
end

function MOI.supports_constraint(
    ::Optimizer,
    func::Type{VAF{T}},
    set::Type{OS},
) where {A,T<:Real,IS<:MOI.AbstractScalarSet,OS<:MOI.IndicatorSet{A,IS}}
    if IS <: MOI.GreaterThan || IS <: CPE.Strictly{MOI.GreaterThan{T}}
        return false
    end
    return A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO
end

supports_inner_constraint(optimizer::Optimizer, func, set) = MOI.supports_constraint(optimizer, func, set)
supports_inner_constraint(optimizer::Optimizer, func::Type{VAF{T}}, ::Type{CPE.AllDifferent}) where T = false
supports_inner_constraint(optimizer::Optimizer, func::Type{VAF{T}}, ::Type{TableSetInternal}) where T = false

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{<:MOI.AbstractFunction},
    set::Type{OS},
) where {A,F,IS,OS<:CS.IndicatorSet{A,F,IS}}
    !(A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO) && return false
    return supports_inner_constraint(optimizer, F, IS)
end

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{<:MOI.AbstractFunction},
    set::Type{OS}
) where {A,F,IS,OS<:CS.ReifiedSet{A,F,IS}}
    !(A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO) && return false
    return supports_inner_constraint(optimizer, F, IS)
end

function MOI.supports_constraint(
    optimizer::Optimizer,
    func::Type{<:Union{VAF{T},MOI.VectorOfVariables}},
    set::Type{BS},
) where {T,F1,F2,F1dim,F2dim,S1,S2,BS<:CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}}
    return is_boolset_supported(optimizer, set)
end

MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{GeqSetInternal},
) = true

"""
    Return whether the two constraint inside the `AbstractBoolSet` are supported directly by the solver
"""
function is_boolset_supported(optimizer::Optimizer, ::Type{<:CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}}) where {F1, F2, F1dim, F2dim, S1, S2}
    is_supported_left = MOI.supports_constraint(optimizer, F1, S1)
    is_supported_right = MOI.supports_constraint(optimizer, F2, S2)
    return is_supported_left && is_supported_right
end

function check_inbounds(model::Optimizer, aff::SAF{T}) where {T<:Real}
    for term in aff.terms
        check_inbounds(model, term.variable_index)
    end
    return
end

function check_inbounds(model::Optimizer, vov::MOI.VectorOfVariables)
    for var in vov.variables
        check_inbounds(model, var)
    end
    return
end

"""
    add_constraint!(model::Optimizer, constraint::Constraint)

Calls `add_constraint!` on the `CoM` without setting pvals.
The variable bounds might not be ready at this stage.
"""
function add_constraint!(model::Optimizer, constraint::Constraint)
    add_constraint!(model.inner, constraint; set_pvals = false)
end

"""
    create_internals(com::CoM, vars, set)

Create ConstraintInternals for a vector of variables constraint
"""
function create_internals(com::CoM, vars, set)
    internals = create_internals(vars, set)
    internals.idx = length(com.constraints) + 1
    return internals
end

"""
    create_internals(vars::MOI.VectorOfVariables, set)

Create ConstraintInternals for a vector of variables constraint
"""
function create_internals(vars::MOI.VectorOfVariables, set)
    ConstraintInternals(
        0,
        vars,
        set,
        Int[v.value for v in vars.variables],
    )
end

"""
    create_internals(func::VAF{T}, set) where {T}

Create ConstraintInternals for a vector of variables constraint
"""
function create_internals(func::VAF{T}, set) where {T}
    ConstraintInternals(
        0,
        func,
        set,
        get_indices(func)
    )
end

"""
    get_complement_constraint(com, constraint)

Return the complement constraint if it exists and `nothing` otherwise.
The complement constraint is the constraint that expresses the opposite i.e
input: 2x + 7 <= 5 => 2x + 7 > 5
 - it will actually output only less than constraints not great than as it's not supported
input 5x == 2 => 5x != 2

Currently it's only implemented for linear constraints
"""
function get_complement_constraint(com, constraint::Constraint)
    return nothing
end

function get_complement_constraint(com, constraint::LinearConstraint{T}) where T
    set = constraint.set
    complement_fct = nothing
    complement_set = nothing
    if constraint.set isa MOI.LessThan
        complement_fct = MOIU.operate(-, T, constraint.fct)
        complement_set = CPE.Strictly(MOI.LessThan(-set.upper))
    elseif constraint.set isa CPE.Strictly{MOI.LessThan{T}, T}
        complement_fct = MOIU.operate(-, T, constraint.fct)
        complement_set = MOI.LessThan(-set.set.upper)
    elseif constraint.set isa MOI.EqualTo
        complement_fct = copy(constraint.fct)
        complement_set = CPE.DifferentFrom(set.value)
    elseif constraint.set isa CPE.DifferentFrom
        complement_fct = copy(constraint.fct)
        complement_set = MOI.EqualTo(set.value)
    end
    complement_lc = get_linear_constraint(complement_fct, complement_set)
    return complement_lc
end

function get_complement_constraint(com, constraint::BoolConstraint)
    return complement_bool_constraint(com, typeof(constraint.set), constraint.fct, constraint.lhs, constraint.rhs)
end

"""
    complement_bool_constraint(com, bst::Type{<:AbstractBoolSet}, fct, lhs_constraint::Constraint, rhs_constraint::Constraint)

Return the complement constraint
"""
function complement_bool_constraint(com, bset::Type{<:AbstractBoolSet}, fct, lhs_constraint::Constraint, rhs_constraint::Constraint)
    cs = complement_set(bset)
    # if there is no complement set check for demorgan complement sets
    if cs === nothing
        dcs = demorgan_complement_set(bset)
        if dcs === nothing
            return nothing # if no complement set option exists
        end
        # demorgan rule -> complement the inners and then apply the demorgan set
        lhs_complement_constraint = get_complement_constraint(com, lhs_constraint)
        rhs_complement_constraint = get_complement_constraint(com, rhs_constraint)
        if lhs_complement_constraint === nothing || rhs_complement_constraint === nothing
            return nothing
        end
        T = parametric_type(com)
        fct = MOIU.operate(vcat, T, lhs_complement_constraint.fct, rhs_complement_constraint.fct)
        set = dcs{typeof(lhs_complement_constraint.fct), typeof(rhs_complement_constraint.fct)}(lhs_complement_constraint.set, rhs_complement_constraint.set)
        internals = ConstraintInternals(0,fct,set,get_indices(fct))
        return demorgan_complement_constraint_type(bset)(com, internals, lhs_complement_constraint, rhs_complement_constraint)
    else
        set = cs{typeof(lhs_constraint.fct), typeof(rhs_constraint.fct)}(lhs_constraint.set, rhs_constraint.set)
        internals = ConstraintInternals(0,fct,set,get_indices(fct))
        return complement_constraint_type(bset)(com, internals, lhs_constraint, rhs_constraint)
    end
end

"""
    MOI.add_constraint(
        model::Optimizer,
        vars::MOI.VectorOfVariables,
        set::MOI.AbstractVectorSet,
    )

Add all kinds of vector of variables constraints like:
TableConstraint and AllDifferentConstraint
"""
function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::MOI.AbstractVectorSet,
)
    check_inbounds(model, vars)
    com = model.inner

    internals = create_internals(com, vars, set)

    constraint = init_constraint_struct(com, set, internals)

    add_constraint!(model, constraint)
    if set isa CPE.AllDifferent
        com.info.n_constraint_types.alldifferent += 1
    elseif set isa TableSetInternal
        com.info.n_constraint_types.table += 1
    elseif set isa CPE.AllEqual
        com.info.n_constraint_types.equality += 1
    end

    return MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(length(com.constraints))
end

"""
    function MOI.add_constraint(
        model::Optimizer,
        vaf::VAF{T},
        set::Union{CPE.AllDifferent, TableSetInternal},
    ) where T

`VectorAffineFunction` constraint for `AllDifferentSet` or `TableSet` to support for things like
`[a+1, b+2, c+3] in CS.AllDifferent()` or
`[a, b, 4] in TableSet(...)`
"""
function MOI.add_constraint(
    model::Optimizer,
    vaf::VAF{T},
    set::Union{CPE.AllDifferent, TableSetInternal},
) where T
    fs = MOIU.eachscalar(vaf)
    variables = Vector{MOI.VariableIndex}(undef, length(fs))
    for (i,f) in enumerate(fs)
        # we need to create a new variable and SAF constraint when it's not a SVF
        if !is_svf(f)
            discrete, non_continuous_value = is_discrete_saf(f)
            if !discrete
                throw(DomainError(non_continuous_value, "The constant and all coefficients need to be discrete"))
            end
            vidx = MOI.add_variable(model)
            variables[i] = vidx
            min_val, max_val = get_extrema(model, f)
            svf = MOI.SingleVariable(vidx)
            MOI.add_constraint(model, svf, MOI.Integer())
            MOI.add_constraint(model, svf, MOI.Interval(min_val, max_val))
            new_constraint_fct = MOIU.operate(-, T, f, svf)
            MOI.add_constraint(model, new_constraint_fct, MOI.EqualTo(0.0))
        else
            variables[i] = f.terms[1].variable_index
        end
    end
    ci = MOI.add_constraint(model, MOI.VectorOfVariables(variables), set)
    return MOI.ConstraintIndex{VAF{T},typeof(set)}(ci.value)
end

"""
    MOI.add_constraint(
        model::Optimizer,
        vars::MOI.AbstractFunction,
        set::ComplementSet,
    )

Add a complement constraint
"""
function MOI.add_constraint(
    model::Optimizer,
    fct::F,
    set::ComplementSet{CF},
) where {F<:MOI.AbstractFunction,CF}
    com = model.inner
    inner_set = set.set
    if CF <: SAF
        fct = get_saf(fct)
    end
    constraint = get_constraint(com, fct, inner_set)
    complement_constraint = get_complement_constraint(com, constraint)
    complement_constraint.idx = length(com.constraints)+1

    add_constraint!(model, complement_constraint)
    return MOI.ConstraintIndex{F,typeof(set)}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::MOI.EqualTo{T},
) where {T<:Real}
    check_inbounds(model, func)

    if length(func.terms) == 1
        vidx = func.terms[1].variable_index.value
        val = convert(Int, (set.value - func.constant) / func.terms[1].coefficient)
        push!(model.inner.init_fixes, (vidx, val))
        return MOI.ConstraintIndex{SAF{T},MOI.EqualTo{T}}(0)
    elseif length(func.terms) == 2 && set.value == zero(T) && func.constant == zero(T)
        if func.terms[1].coefficient == -func.terms[2].coefficient
            # we have the form a == b
            vecOfvar = MOI.VectorOfVariables([
                func.terms[1].variable_index,
                func.terms[2].variable_index,
            ])
            com = model.inner
            internals = ConstraintInternals(
                length(com.constraints) + 1,
                vecOfvar,
                CS.CPE.AllEqual(2),
                Int[v.value for v in vecOfvar.variables],
            )
            constraint = EqualConstraint(internals, ones(Int, 2))

            add_constraint!(model, constraint)
            com.info.n_constraint_types.equality += 1
            return MOI.ConstraintIndex{SAF{T},MOI.EqualTo{T}}(length(model.inner.constraints))
        end
    end

    lc = new_linear_constraint(model.inner, func, set)

    add_constraint!(model, lc)
    model.inner.info.n_constraint_types.equality += 1

    return MOI.ConstraintIndex{SAF{T},MOI.EqualTo{T}}(length(model.inner.constraints))
end

function add_variable_less_than_variable_constraint(
    model::Optimizer,
    func::SAF{T},
    set::MOI.LessThan{T},
) where {T<:Real}
    reverse_order = false
    if func.terms[1].coefficient != 1.0 || func.terms[2].coefficient != -1.0
        if func.terms[1].coefficient == -1.0 && func.terms[2].coefficient == 1.0
            # rhs is lhs and other way around
            reverse_order = true
        end
    end

    com = model.inner

    if reverse_order
        lhs = func.terms[2].variable_index.value
        rhs = func.terms[1].variable_index.value
    else
        lhs = func.terms[1].variable_index.value
        rhs = func.terms[2].variable_index.value
    end

    internals = ConstraintInternals(
        length(model.inner.constraints) + 1, # constraint idx
        func,
        set,
        [lhs, rhs],
    )
    svc = SingleVariableConstraint(internals, lhs, rhs)

    add_constraint!(model, svc)
    com.info.n_constraint_types.inequality += 1

    return MOI.ConstraintIndex{SAF{T},typeof(set)}(length(model.inner.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::MOI.LessThan{T},
) where {T<:Real}
    check_inbounds(model, func)

    # support for a <= b which is written as a-b <= 0
    # currently only supports if the coefficients are 1 and -1
    if set.upper == 0.0 &&
       length(func.terms) == 2 &&
       abs(func.terms[1].coefficient) == 1.0 &&
       abs(func.terms[2].coefficient) == 1.0 &&
       func.terms[1].coefficient == -func.terms[2].coefficient
        return add_variable_less_than_variable_constraint(model, func, set)
    end

    lc = new_linear_constraint(model.inner, func, set)

    add_constraint!(model, lc)
    model.inner.info.n_constraint_types.inequality += 1

    return MOI.ConstraintIndex{SAF{T},typeof(set)}(length(model.inner.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::CPE.DifferentFrom{T},
) where {T<:Real}
    com = model.inner
    com.info.n_constraint_types.notequal += 1

    # support for cx != a where a and c are constants
    if length(func.terms) == 1
        # reformulate to x != a where x is variable and a a constant
        rhs = set.value / func.terms[1].coefficient
        if isapprox_discrete(com, rhs)
            rhs = get_approx_discrete(rhs)
            rm!(
                com,
                com.search_space[func.terms[1].variable_index.value],
                rhs;
                changes = false,
            )
        end
        return MOI.ConstraintIndex{SAF{T},CPE.DifferentFrom{T}}(0)
    end

    lc = new_linear_constraint(model.inner, func, set)

    add_constraint!(model, lc)

    return MOI.ConstraintIndex{SAF{T},CPE.DifferentFrom{T}}(length(com.constraints))
end
function MOI.add_constraint(
    model::Optimizer,
    func::SAF{T},
    set::CPE.Strictly{MOI.LessThan{T}},
) where {T<:Real}
    check_inbounds(model, func)

    lc = new_linear_constraint(model.inner, func, set)

    add_constraint!(model, lc)
    model.inner.info.n_constraint_types.inequality += 1

    return MOI.ConstraintIndex{SAF{T},typeof(set)}(length(model.inner.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::IS,
) where {A,T<:Real,ASS<:MOI.AbstractScalarSet,IS<:MOI.IndicatorSet{A,ASS}}
    com = model.inner
    com.info.n_constraint_types.indicator += 1

    indices = get_indices(func)

    internals = ConstraintInternals(
        length(com.constraints)+1,
        func,
        MOI.IndicatorSet{A}(set.set),
        indices,
    )

    activator_internals = get_activator_internals(A, indices)

    lc = get_inner_constraint(com, func, set, set.set)

    constraint = IndicatorConstraint(internals, activator_internals, lc)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},MOI.IndicatorSet{A,ASS}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::IS,
) where {A,F,S,IS<:CS.IndicatorSet{A,F,S}}
    com = model.inner
    com.info.n_constraint_types.indicator += 1

    internals = create_internals(com, vars, set)

    inner_constraint = get_inner_constraint(com, vars, set, set.set)

    indices = internals.indices
    activator_internals = get_activator_internals(A, indices)
    constraint =
        IndicatorConstraint(internals, activator_internals, inner_constraint)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{MOI.VectorOfVariables,CS.IndicatorSet{A,F,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::IS,
) where {T,A,F,S<:MOI.AbstractVectorSet,IS<:CS.IndicatorSet{A,F,S}}
    com = model.inner
    com.info.n_constraint_types.reified += 1

    internals = create_internals(com, func, set)

    inner_constraint = get_inner_constraint(com, func, set, set.set)
    indices = internals.indices

    activator_internals = get_activator_internals(A, indices)

    constraint =
        IndicatorConstraint(internals, activator_internals, inner_constraint)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},CS.IndicatorSet{A,F,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::RS,
) where {A,F,S,T<:Real,RS<:ReifiedSet{A,F,S}}
    com = model.inner
    com.info.n_constraint_types.reified += 1

    indices = get_indices(func)

    internals = ConstraintInternals(
        length(com.constraints)+1,
        func,
        typeof(set)(set.set, set.dimension),
        indices,
    )

    inner = get_inner_constraint(com, func, set, set.set)
    complement_inner = get_complement_constraint(com, inner)

    activator_internals = get_activator_internals(A, indices)
    constraint = ReifiedConstraint(internals, activator_internals, inner, complement_inner)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},CS.ReifiedSet{A,F,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::RS,
) where {A,F,S<:MOI.AbstractVectorSet,RS<:CS.ReifiedSet{A,F,S}}
    com = model.inner
    com.info.n_constraint_types.reified += 1
    internals = create_internals(com, vars, set)

    inner_constraint = get_inner_constraint(com, vars, set, set.set)
    complement_constraint = get_complement_constraint(model, inner_constraint)
    indices = internals.indices
    activator_internals = get_activator_internals(A, indices)
    constraint =
        ReifiedConstraint(internals, activator_internals, inner_constraint, complement_constraint)

    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{MOI.VectorOfVariables,CS.ReifiedSet{A,F,S}}(length(com.constraints))
end

function MOI.add_constraint(
    model::Optimizer,
    func::VAF{T},
    set::BS,
) where {T,BS<:AbstractBoolSet}
    com = model.inner
    internals = create_internals(com, func, set)
    constraint = init_constraint_struct(com, set, internals)
    add_constraint!(model, constraint)

    return MOI.ConstraintIndex{VAF{T},typeof(set)}(length(com.constraints))
end


text:
### Streaming method
#   only implemented for quadratic objectives
#   TODO: add quadratic regularization

export StreamingParams, streaming_fit!, streaming_impute!

mutable struct StreamingParams<:AbstractParams
    T0::Int # number of rows to use to initialize Y before streaming begins
    stepsize::Float64 # stepsize (inverse of memory)
    Y_update_interval::Int # how often to prox Y
end
function StreamingParams(
    T0::Int=1000; # number of rows to use to initialize Y before streaming begins
    stepsize::Number=1/T0, # (inverse of memory)
    Y_update_interval::Int=10 # how often to prox Y

)
    return StreamingParams(T0, convert(Float64, stepsize), Y_update_interval)
end

### FITTING
function streaming_fit!(glrm::GLRM, params::StreamingParams=StreamingParams();
    ch::ConvergenceHistory=ConvergenceHistory("StreamingGLRM"),
    verbose=true)

  # make sure everything is quadratic
  @assert all(map(l->isa(l, QuadLoss), glrm.losses))
  @assert all(map(l->isa(l, QuadReg), glrm.rx))
  @assert all(map(l->isa(l, QuadReg), glrm.ry))

  # initialize Y and first T0 rows of X
  init_glrm = keep_rows(glrm, params.T0)
  init_svd!(init_glrm)
  copy!(glrm.Y, init_glrm.Y)
  copy!(view(glrm.X, :, 1:params.T0), init_glrm.X)

  ### initialization
  A = glrm.A # rename these for easier local access
  rx = glrm.rx
  ry = glrm.ry
  X = glrm.X; Y = glrm.Y
  k = glrm.k
  m,n = size(A)

  # yscales = map(r->r.scale, ry)

  for i=params.T0+1:m
    # update x_i
    obs = glrm.observed_features[i]
    Yobs = Y[:, obs]
    Aobs = A[i, obs]
    xi = view(X, :, i)

    copy!(xi, (Yobs * Yobs' + 2 * rx[i].scale * I) \ (Yobs * Aobs))

    # update objective
    r = Yobs'*xi - Aobs
    push!(ch.objective, norm(r) ^ 2)

    # # update Y
    # TODO verify this is stochastic proximal gradient (with constant stepsize) for the problem
    # TODO don't prox Y at every iteration
    # TODO don't assume scales on all the rys are equal
    # gY[:, jj] = xi * r' == r[jj] * xi # gradient of ith row objective wrt Y
    for jj in 1:length(obs)
      Y[:,obs[jj]] -= params.stepsize * r[jj] * xi
    end
    if i%params.Y_update_interval == 0
      # prox!(ry, Y, params.stepsize * params.Y_update_interval)
      Y ./= (1 + 2 * params.stepsize * params.Y_update_interval * ry[1].scale)
    end
  end

  return X, Y, ch
end

### FITTING
function streaming_impute!(glrm::GLRM, params::StreamingParams=StreamingParams();
    ch::ConvergenceHistory=ConvergenceHistory("StreamingGLRM"),
    verbose=true)

  # make sure everything is quadratic
  @assert all(map(l->isa(l, QuadLoss), glrm.losses))
  @assert all(map(l->isa(l, QuadReg), glrm.rx))
  @assert all(map(l->isa(l, QuadReg), glrm.ry))

  # initialize Y and first T0 rows of X
  init_glrm = keep_rows(glrm, params.T0)
  init_svd!(init_glrm)
  copy!(glrm.Y, init_glrm.Y)
  copy!(view(glrm.X, :, 1:params.T0), init_glrm.X)

  ### initialization
  A = glrm.A # rename these for easier local access
  Ahat = copy(glrm.A)
  rx = glrm.rx
  ry = glrm.ry
  X = glrm.X; Y = glrm.Y
  k = glrm.k
  m,n = size(A)

  # yscales = map(r->r.scale, ry)

  for i=params.T0+1:m
    # update x_i
    obs = glrm.observed_features[i]
    Yobs = Y[:, obs]
    Aobs = A[i, obs]
    xi = view(X, :, i)

    copy!(xi, (Yobs * Yobs' + 2 * rx[i].scale * I) \ (Yobs * Aobs))

    # impute
    not_obs = setdiff(Set(1:n), Set(obs))
    if length(not_obs)>0
      ahat = xi'*Y
      Ahat[i, not_obs] = ahat[not_obs]
    end

    # update objective
    r = Yobs'*xi - Aobs
    push!(ch.objective, norm(r) ^ 2)

    # # update Y
    # TODO verify this is stochastic proximal gradient (with constant stepsize) for the problem
    # TODO don't prox Y at every iteration
    # TODO don't assume scales on all the rys are equal
    # gY[:, jj] = xi * r' == r[jj] * xi # gradient of ith row objective wrt Y
    for jj in 1:length(obs)
      Y[:,obs[jj]] -= params.stepsize * r[jj] * xi
    end
    if i%params.Y_update_interval == 0
      # prox!(ry, Y, params.stepsize * params.Y_update_interval)
      Y ./= (1 + 2 * params.stepsize * params.Y_update_interval * ry[1].scale)
    end
  end

  return Ahat
end

""" Constructs new GLRM on subset of rows of the data from input glrm """
function keep_rows(glrm, r::UnitRange{Int})
  @assert maximum(r) <= size(glrm.A, 1)
  obs = flatten_observations(glrm.observed_features)
  first_row = minimum(r)
  if first_row > 1
    new_obs = map( t -> (t[1]-first_row+1, t[2]), filter( t -> (t[1] in r), obs))
  else
    new_obs = filter( t -> (t[1] in r), obs)
  end
  of, oe = sort_observations(new_obs, length(r), size(glrm.A, 2))
  new_glrm = GLRM(glrm.A[r,:], glrm.losses, glrm.rx[r], glrm.ry, glrm.k,
                  observed_features = of, observed_examples = oe)
  return new_glrm
end
keep_rows(glrm, T::Int) = keep_rows(glrm, 1:T)


text:
module SeisDetective

using SeisIO
using Statistics

include("trigger.jl")

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule hsa_rocr_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("hsa_rocr")
JLLWrappers.@generate_main_file("hsa_rocr", UUID("dd59ff1a-a01a-568d-8b29-0669330f116a"))
end  # module hsa_rocr_jll


text:
@testset "circle" begin
    p = Point([-75.343, 39.984])

    @test circle(center=p, radius=5).coordinates ≈ [
					[[
							-75.343,
							40.028966
						],
						[
							-75.348756,
							40.028749
						],
						[
							-75.354456,
							40.028101
						],
						[
							-75.360046,
							40.027029
						],
						[
							-75.365472,
							40.025541
						],
						[
							-75.37068,
							40.023653
						],
						[
							-75.375622,
							40.021383
						],
						[
							-75.380248,
							40.018753
						],
						[
							-75.384516,
							40.015788
						],
						[
							-75.388383,
							40.012517
						],
						[
							-75.391813,
							40.008972
						],
						[
							-75.394772,
							40.005185
						],
						[
							-75.397232,
							40.001195
						],
						[
							-75.399169,
							39.997039
						],
						[
							-75.400565,
							39.992758
						],
						[
							-75.401406,
							39.988393
						],
						[
							-75.401685,
							39.983985
						],
						[
							-75.401399,
							39.979578
						],
						[
							-75.40055,
							39.975213
						],
						[
							-75.399148,
							39.970934
						],
						[
							-75.397204,
							39.96678
						],
						[
							-75.39474,
							39.962792
						],
						[
							-75.391777,
							39.959008
						],
						[
							-75.388345,
							39.955465
						],
						[
							-75.384477,
							39.952197
						],
						[
							-75.380211,
							39.949235
						],
						[
							-75.375586,
							39.946608
						],
						[
							-75.370648,
							39.94434
						],
						[
							-75.365444,
							39.942455
						],
						[
							-75.360025,
							39.940969
						],
						[
							-75.354442,
							39.939897
						],
						[
							-75.348748,
							39.93925
						],
						[
							-75.343,
							39.939034
						],
						[
							-75.337252,
							39.93925
						],
						[
							-75.331558,
							39.939897
						],
						[
							-75.325975,
							39.940969
						],
						[
							-75.320556,
							39.942455
						],
						[
							-75.315352,
							39.94434
						],
						[
							-75.310414,
							39.946608
						],
						[
							-75.305789,
							39.949235
						],
						[
							-75.301523,
							39.952197
						],
						[
							-75.297655,
							39.955465
						],
						[
							-75.294223,
							39.959008
						],
						[
							-75.29126,
							39.962792
						],
						[
							-75.288796,
							39.96678
						],
						[
							-75.286852,
							39.970934
						],
						[
							-75.28545,
							39.975213
						],
						[
							-75.284601,
							39.979578
						],
						[
							-75.284315,
							39.983985
						],
						[
							-75.284594,
							39.988393
						],
						[
							-75.285435,
							39.992758
						],
						[
							-75.286831,
							39.997039
						],
						[
							-75.288768,
							40.001195
						],
						[
							-75.291228,
							40.005185
						],
						[
							-75.294187,
							40.008972
						],
						[
							-75.297617,
							40.012517
						],
						[
							-75.301484,
							40.015788
						],
						[
							-75.305752,
							40.018753
						],
						[
							-75.310378,
							40.021383
						],
						[
							-75.31532,
							40.023653
						],
						[
							-75.320528,
							40.025541
						],
						[
							-75.325954,
							40.027029
						],
						[
							-75.331544,
							40.028101
						],
						[
							-75.337244,
							40.028749
						],
						[
							-75.343,
							40.028966
						]
					]] atol = 0.05

	p2 = Point([11.343, 44.495])
	@test sector(p2, 5., 20., 60.).coordinates ≈ [[[
							11.343,
							44.495
						],
						[
							11.364576,
							44.537252
						],
						[
							11.370282,
							44.53554
						],
						[
							11.375724,
							44.533437
						],
						[
							11.380851,
							44.530964
						],
						[
							11.385612,
							44.528145
						],
						[
							11.389963,
							44.525006
						],
						[
							11.39386,
							44.521578
						],
						[
							11.397267,
							44.517894
						],
						[
							11.397614,
							44.51747
						],
						[
							11.343,
							44.495]]]
end


text:
using JuMP
using Gurobi
using Plasmo
using PlasmoAlgorithms

##Place MP and SP into ModelGraph
mp = Model(solver = GurobiSolver())
sp = Model(solver = GurobiSolver())

@variable(mp,y>=1)
@objective(mp,Min,2y)

@variable(sp,x[1:2]>=0)
@variable(sp,y>=0)
@constraint(sp,x[1:2].<=1)
@constraint(sp,x[1]+2x[2]+y>=5)
@objective(sp,Min,2x[1]+3x[2])

## Plasmo Graph
g = ModelGraph()
setsolver(g, GurobiSolver())
n1 = add_node(g)
setmodel(n1,mp)
n2 = add_node(g)
setmodel(n2,sp)

##Set n2 as a child node of n1
edge = Plasmo.add_edge(g,n1,n2)

## Linking constraints between MP and SP
@linkconstraint(g, n1[:y] == n2[:y])

bendersolve(g, max_iterations=20)


text:
module RTMidi

using rtmidi_jll
using CEnum
using MIDI

include("api.jl")

export MidiOut, MidiIn, port_count, open_virtual_port!, listen!, stop_listening!

mutable struct MidiOut
    name::String
    ptr::Ptr{RtMidiWrapper}
    function MidiOut(name, api::RtMidiApi = RTMIDI_API_UNSPECIFIED)
        name = convert(String, name)
        ptr = rtmidi_out_create(api, name)
        midiout = new(name, ptr)
        finalizer(midiout) do m
            rtmidi_out_free(m.ptr)
        end
        midiout
    end
end

mutable struct MidiIn
    name::String
    ptr::Ptr{RtMidiWrapper}
    timer::Union{Nothing, Timer}
    function MidiIn(name, api::RtMidiApi = RTMIDI_API_UNSPECIFIED; queue_size_limit = 1024)
        name = convert(String, name)
        ptr = rtmidi_in_create(api, name, queue_size_limit)
        midiin = new(name, ptr, nothing)
        finalizer(midiin) do m
            rtmidi_out_free(m.ptr)
        end
        midiin
    end
end

function _check_status(m)
    wrapper = unsafe_load(m.ptr)
    if wrapper.ok == 0
        msg = unsafe_string(wrapper.msg)
        error(msg)
    end
    nothing
end

function port_count(m)::Int
    count = rtmidi_get_port_count(m.ptr)
    _check_status(m)
    count
end

function open_virtual_port!(m, name)
    rtmidi_open_virtual_port(m.ptr, name)
    _check_status(m)
    nothing
end

function listen!(callback, mi::MidiIn; interval_seconds = 0.001)
    if mi.timer !== nothing
        error("MidiIn is already listening.")
    end

    bufferlength = 1024
    vector = zeros(UInt8, bufferlength) # 1024 is maximal SYSEX size
    sz = Ref{Cint}()
    laststatus = 0x00 # is this correct?

    mi.timer = Timer(0, interval = interval_seconds) do timer
        sz[] = bufferlength
        deltatime = RTMidi.rtmidi_in_get_message(mi.ptr, vector, sz)
        sz[] == 0 && return
        event = MIDI.readMIDIevent(0, IOBuffer(vector), laststatus)
        laststatus = event.status
        callback(event)
    end
    return
end

function stop_listening!(mi::MidiIn)
    mi.timer === nothing && return
    close(mi.timer)
    mi.timer = nothing
    return
end

####################

end


text:
using PackageCompiler

builddir = ARGS[1]
mkpath(builddir)

build_shared_lib(
    joinpath(@__DIR__, "staticrbd.jl"), "libstaticrbdjl",
    snoopfile = joinpath(@__DIR__, "snoopfile.jl"),
    startup_file = "no",
    builddir = builddir,
    Release = true,
    check_bounds = "no",
    optimize = 3
)

max_stars_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_issues_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_forks_repo_name: CobbledSteel/rbd-benchmarks-riscv


text:
#
# This file just illustrates the nuts-and-bolts of the approximate inference API provided
# by Stheno. It's worth mentally contrasting this with what you would have written to
# perform exact inference (hint: it's remarkably similar).
#

# Set up the environment to run this example. Make sure you're within the folder that this
# file lives in.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

# The time-to-first-plot issue means that this might take a little while.
using AbstractGPs, BenchmarkTools, LinearAlgebra, Plots, Random, Stheno
rng = MersenneTwister(123456);

# Construct a Matern-5/2 kernel with lengthscale 0.5 and variance 1.2.
k = 1.2 * transform(Matern52Kernel(), 2.0);

# Construct a zero-mean GP with a simple kernel. Don't worry about the GPC object.
f = Stheno.wrap(GP(k), GPC());

# Specify some locations at which to consider the GP.
N = 5_000
x = vcat(rand(rng, div(N, 2)) * 3, rand(rng, div(N, 2)) * 4 .+ 6)

# Specify the variance of the noise under which we'll make observations of the GP.
# We could also have made this a `Real` to specify isotropic noise.
Σ = Diagonal(rand(rng, N) .+ 0.1)

# Construct marginal distribution over `f` at `x` added to some independent zero-mean
# Gaussian noise with covariance matrix `Σ`.
fx = f(x, Σ)

# Generate a sample from the prior.
y = rand(rng, fx)

# Now pretend that we can't perform inference because `N` is too large.

# Specify some pseudo-point locations.
z = range(-3.0, 13.0; length=150)

# Specify inducing-points. Add a small amount of jitter for numerical stability.
u = f(z, 1e-6)

# Compute the Evidence Lower BOund.
println("The elbo is reasonable tight.")
@show elbo(fx, y, u), logpdf(fx, y)

println("Benchmark logpdf")
display(@benchmark logpdf($fx, $y))
println()

println("Benchmark elbo")
display(@benchmark elbo($fx, $y, $u))
println()

# Compute the approximate posterior process.
f_post = approx_posterior(VFE(), fx, y, u);

# Specify some points at which to plot the approximate posterior.
Npr = 1000
xpr = range(-3.0, 13.0; length=Npr)

# It's possible to work efficently with the approximate posterior marginals.
# Unfortunately it's not possible to work efficiently with the entire posterior process --
# this is limitation of pseudo-point approximations generally, as opposed to a limitation
# of Stheno.
posterior_marginals = marginals(f_post(xpr))

# Visualise the posterior. At the time of writing, it remains important to do this manually
# for the sake of efficiency. If you would like to have a high-level interface similar to
# the one available in the exact inference setting, please feel free to raise an issue, or
# implement it yourself and open a PR!

# Specify our plotting backend.
gr();

# Construct a new plot object.
posterior_plot = plot(legend=nothing);

# Generate several samples from the posterior predictive distribution.
plot!(posterior_plot, f_post(xpr); color=:blue);

# Plot the observations.
scatter!(posterior_plot, x, y;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=1,
    markeralpha=0.2,
);

# Show the plot.
display(posterior_plot);


text:
# -----------------------------------------------------------------------------
# JuLIP, ACE, etc : just use import throughout, this avoids bugs

import Base: ==, length

import JuLIP

import JuLIP: SitePotential,
              cutoff,
              z2i, i2z, numz,
              AbstractCalculator,
              Atoms,
              chemical_symbol, AtomicNumber,
              JVec, JMat,
              energy, forces, virial

import JuLIP.MLIPs: IPBasis

import JuLIP.Potentials: ZList, SZList, zlist

import ACEbase: evaluate, evaluate_d,
                evaluate!, evaluate_d!,
                read_dict, write_dict,
                valtype,
                combine,
                ScalarACEBasis, _allfieldsequal

import ACE
import ACE: AbstractState,
            AbstractDiscreteState,
            AbstractContinuousState,
            AbstractConfiguration,
            scaling, PolyTransform,
            transformed_jacobi


text:
@testset "constant" begin
    @testset "ZeroKernel" begin
        k = ZeroKernel()
        @test eltype(k) == Any
        @test kappa(k,2.0) == 0.0
        @test KernelFunctions.metric(ZeroKernel()) == KernelFunctions.Delta()
        @test repr(k) == "Zero Kernel"
        test_ADs(ZeroKernel)
    end
    @testset "WhiteKernel" begin
        k = WhiteKernel()
        @test eltype(k) == Any
        @test kappa(k,1.0) == 1.0
        @test kappa(k,0.0) == 0.0
        @test EyeKernel == WhiteKernel
        @test metric(WhiteKernel()) == KernelFunctions.Delta()
        @test repr(k) == "White Kernel"
        test_ADs(WhiteKernel)
    end
    @testset "ConstantKernel" begin
        c = 2.0
        k = ConstantKernel(c=c)
        @test eltype(k) == Any
        @test kappa(k,1.0) == c
        @test kappa(k,0.5) == c
        @test metric(ConstantKernel()) == KernelFunctions.Delta()
        @test metric(ConstantKernel(c=2.0)) == KernelFunctions.Delta()
        @test repr(k) == "Constant Kernel (c = $(c))"
        test_ADs(c->ConstantKernel(c=first(c)), [c])
    end
end


text:
using HMCExamples

HMCExamples.main_FVGQ_1_joint(ARGS)

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
module DiodeClipper

using Interpolations
using NLsolve: nlsolve
using NumericalIntegration: cumul_integrate

export runsim, runsimaaout, runsimaa

# component values
const R = 1e3
const C = 33e-9

# nonlinear function mapping p = x+u to y (depending on used sampling rate)
function nlfunction(fS, p)
    Is = 1e-15
    Vt = 25e-3
    function f!(F, x)
        vD=x[1]
        F[1] = p - (2*R*C*fS + 1)*vD - 2*R*Is*sinh(vD/Vt)
    end
    return nlsolve(f!, [0.0]).zero[1]
end

# precompute antiderivate table for given sampling rate
function makeantiderivtable(fS)
    pmax = 1.5 * (4 + 8R*C*fS) # guesstimate of relevant range
    ps = range(-pmax, pmax, length=1024)
    zs = [nlfunction(fS, p) for p in ps]
    zsint = cumul_integrate(ps, zs)
    return scale(interpolate(zsint, BSpline(Cubic(Line(OnGrid())))), ps)
end

# simulation without any antialiasing
function runsim(fS, u)
    y = similar(u) # preallocate output
    x = 0.0 # initial state
    for n in eachindex(u)
        p = x + u[n]
        y[n] = nlfunction(fS, p)
        x = -x + 4*C*R*fS*y[n]
    end
    return y
end

# simulation with antialiasing of the output only
function runsimaaout(fS, u)
    y = similar(u)
    adtab = makeantiderivtable(fS)
    x = 0.0
    last_p = 0.0
    last_ad = adtab(last_p)
    for n in eachindex(u)
        p = R*x + u[n]
        ad = adtab(p)
        if abs(p-last_p) > 1e-15
            y[n] = (ad - last_ad) / (p - last_p)
        else
            y[n] = nlfunction(fS, 0.5 * (p + last_p))
        end
        last_p = p
        last_ad = ad
        x = -x + 4*C*fS*nlfunction(fS, p)
    end
    return y
end

# simulation with full antialiasing
function runsimaa(fS, u)
    fS = fS/1.5
    y = similar(u)
    adtab = makeantiderivtable(fS)
    x = 0.0
    last_p = 0.0
    last_ad = adtab(last_p)
    prev_x = 0.0
    for n in eachindex(u)
        p = R*x + u[n]
        ad = adtab(p)
        if abs(p-last_p) > 1e-15
            y[n] = (ad - last_ad) / (p - last_p)
        else
            y[n] = nlfunction(fS, 0.5 * (p + last_p))
        end
        last_p = p
        last_ad = ad
        prev_x, x = x, -0.5*(x + prev_x) + 4fS*C*y[n]
    end
    return y
end

end


text:
using NLExample
using Test

@testset "NLExample.jl" begin
    # Write your own tests here.
end


text:
type Node
    parent # parent Id of this node, -1 means root.
    child # child node Ids.
    uList # neighbour list (including itself), only leaf has
    vList # far field list
    wList # near field, only leaf has
    xList # far field list

    nLevel # level of this node
    nUList # number of neighbours
    nVList # number of V far field
    nWList # number of W near field
    nXList # number of X far field

    nodeIndex # numbering of this node in its parent's child

    center # center coordinate of this node
    radius # radius of this node

    nSource # source points number
    nTarget # target points  number

    sourceIndex # index for source for locating source
    targetIndex # index for target for locating target

    source # source list
    # normal # normal list
    target # target list

    scaledCnode # cnodes
    cnodeCharge # charge on cnodes
    cnodePotential # potential on cnodes

    R # transfer matrix from cnode to cluster
    L # transfer matrix from cnode to targets

    charge # associated with source
    potential # associated with target

    isLeaf
    isEmpty
    chargeComputed

    Node(nl, ni) = new(-1, Int[], Int[], Int[], Int[], Int[], nl, 0,0,0,0, ni, [Inf, Inf], [Inf, Inf], 0, 0, Int[], Int[], [],[], [],[], Float64[], Float64[], Float64[], Float64[], Float64[], false, false, false)
end


text:
module SimpleGraphs

import Base:
    eltype, show, ==, Pair, Tuple, copy, length, start, next, done, issubset, zero

import LightGraphs:
    _NI, _insert_and_dedup!, AbstractGraph, AbstractEdge, AbstractEdgeIter,
    src, dst, edgetype, nv, ne, vertices, edges, is_directed,
    add_vertex!, add_edge!, rem_vertex!, rem_edge!,
    has_vertex, has_edge, in_neighbors, out_neighbors,

    indegree, outdegree, degree, has_self_loops, num_self_loops

export AbstractSimpleGraph, AbstractSimpleDiGraph, AbstractSimpleEdge,
    SimpleEdge, SimpleGraph, SimpleGraphEdge,
    SimpleDiGraph, SimpleDiGraphEdge

"""
    AbstractSimpleGraph

An abstract type representing a simple graph structure.
AbstractSimpleGraphs must have the following elements:
- vertices::UnitRange{Integer}
- fadjlist::Vector{Vector{Integer}}
- ne::Integer
"""
abstract type AbstractSimpleGraph <: AbstractGraph end

function show(io::IO, g::AbstractSimpleGraph)
    if is_directed(g)
        dir = "directed"
    else
        dir = "undirected"
    end
    if nv(g) == 0
        print(io, "empty $dir simple $(eltype(g)) graph")
    else
        print(io, "{$(nv(g)), $(ne(g))} $dir simple $(eltype(g)) graph")
    end
end

nv(g::AbstractSimpleGraph) = eltype(g)(length(fadj(g)))
vertices(g::AbstractSimpleGraph) = one(eltype(g)):nv(g)

edges(g::AbstractSimpleGraph) = SimpleEdgeIter(g)

fadj(g::AbstractSimpleGraph) = g.fadjlist
fadj(g::AbstractSimpleGraph, v::Integer) = g.fadjlist[v]

badj(x...) = _NI("badj")

# handles single-argument edge constructors such as pairs and tuples
has_edge(g::AbstractSimpleGraph, x) = has_edge(g, edgetype(g)(x))
add_edge!(g::AbstractSimpleGraph, x) = add_edge!(g, edgetype(g)(x))

# handles two-argument edge constructors like src,dst
has_edge(g::AbstractSimpleGraph, x, y) = has_edge(g, edgetype(g)(x, y))
add_edge!(g::AbstractSimpleGraph, x, y) = add_edge!(g, edgetype(g)(x, y))

in_neighbors(g::AbstractSimpleGraph, v::Integer) = badj(g, v)
out_neighbors(g::AbstractSimpleGraph, v::Integer) = fadj(g, v)

function issubset(g::T, h::T) where T<:AbstractSimpleGraph
    (gmin, gmax) = extrema(vertices(g))
    (hmin, hmax) = extrema(vertices(h))
    return (hmin <= gmin <= gmax <= hmax) && issubset(edges(g), edges(h))
end

has_vertex(g::AbstractSimpleGraph, v::Integer) = v in vertices(g)

ne(g::AbstractSimpleGraph) = g.ne

function rem_edge!(g::AbstractSimpleGraph, u::Integer, v::Integer)
    T = eltype(g)
    rem_edge!(g, edgetype(g)(T(u), T(v)))
end

@doc_str """
    rem_vertex!(g, v)

Remove the vertex `v` from graph `g`. Return false if removal fails
(e.g., if vertex is not in the graph); true otherwise.

### Performance
Time complexity is ``\\mathcal{O}(k^2)``, where ``k`` is the max of the degrees
of vertex ``v`` and vertex ``|V|``.

### Implementation Notes
This operation has to be performed carefully if one keeps external
data structures indexed by edges or vertices in the graph, since
internally the removal is performed swapping the vertices `v`  and ``|V|``,
and removing the last vertex ``|V|`` from the graph. After removal the
vertices in `g` will be indexed by ``1:|V|-1``.
"""
function rem_vertex!(g::AbstractSimpleGraph, v::Integer)
    v in vertices(g) || return false
    n = nv(g)

    # remove the in_edges from v
    srcs = copy(in_neighbors(g, v))
    for s in srcs
        rem_edge!(g, edgetype(g)(s, v))
    end
    # remove the in_edges from the last vertex
    neigs = copy(in_neighbors(g, n))
    for s in neigs
        rem_edge!(g, edgetype(g)(s, n))
    end
    if v != n
        # add the edges from n back to v
        for s in neigs
            add_edge!(g, edgetype(g)(s, v))
        end
    end

    if is_directed(g)
        # remove the out_edges from v
        dsts = copy(out_neighbors(g, v))
        for d in dsts
            rem_edge!(g, edgetype(g)(v, d))
        end
        # remove the out_edges from the last vertex
        neigs = copy(out_neighbors(g, n))
        for d in neigs
            rem_edge!(g, edgetype(g)(n, d))
        end
        if v != n
            # add the out_edges back to v
            for d in neigs
                add_edge!(g, edgetype(g)(v, d))
            end
        end
    end

    pop!(g.fadjlist)
    if is_directed(g)
        pop!(g.badjlist)
    end
    return true
end

zero(g::T) where T<:AbstractSimpleGraph = T()

include("simpleedge.jl")
include("simpledigraph.jl")
include("simplegraph.jl")
include("simpleedgeiter.jl")

end # module


text:
using LinearAlgebra: norm

using Colors: HSV

"""
    colorflow(flo; maxflow)

Visualize the given flow field `flo` by creating an image encoding
direction as color and length as saturation (scaled using `maxflow`).
Missing values are encoded as black.
If `maxflow` is omitted, it will be computed from the input.
"""
function colorflow(flo; maxflow = _maxflow(flo))
    CT = HSV{Float32}
    color(x1, x2) = ismissing(x1) || ismissing(x2) ?
        CT(0, 1, 0) :
        CT(180f0 / pi * atan(x1, x2), norm((x1, x2)) / maxflow, 1)
    x1 = selectdim(flo, 1, 1)
    x2 = selectdim(flo, 1, 2)
    return color.(x1, x2)
end

function _maxflow(flo)
    return maximum(CartesianIndices(axes(flo)[2:end])) do I
        v = view(flo, :, I)
        any(ismissing, v) && return zero(eltype(flo))
        return norm(v)
    end
end


text:
"""
    Parareal <: AbstractTimeParallelSolver

A composite type for the parareal algorithm.

# Constructors
```julia
Parareal(finesolver, coarsolver, P, K, control)
Parareal(finesolver, coarsolver; P=10, K=P, control=ErrorControl())
```

# Arguments
- `finesolver :: AbstractInitialValueSolver` : fine solver.
- `coarsolver :: AbstractInitialValueSolver` : coarse solver.
- `P :: Integer` : number of time chunks.
- `K :: Integer` : maximum number of iterations.
- `control :: ErrorControl` : tolerance and error mechanism.

# Functions
- [`show`](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
struct Parareal{finesolver_T, coarsolver_T, P_T, K_T, control_T} <: AbstractTimeParallelSolver
    finesolver::finesolver_T
    coarsolver::coarsolver_T
    P::P_T
    K::K_T
    control::control_T
end

function Parareal(finesolver, coarsolver; P=10, K=P, control=ErrorControl())
    return Parareal(finesolver, coarsolver, P, K, control)
end


text:
# edyna Julia module begins -->
module edyna

# ellipsoids data
ellcount = 0
ellcenter = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
ellradii = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
ellorient = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
             Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
             Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
elldefo = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
           Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
           Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
elldensity = Array{Float32}(undef, ellcount);
ellyoung = Array{Float32}(undef, ellcount);
ellpoisson = Array{Float32}(undef, ellcount);
ellmass = Array{Float32}(undef, ellcount);
elleuler = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
            Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
            Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
ellinveuler = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
               Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
               Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
ellinertia = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
              Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
              Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
ellinvinertia = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
                 Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount),
                 Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
elllinear = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
ellangular = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]

# insert ellipsoid
function ellipsoid(x, a, b, c, R, density, young = 1E9, poisson = 0.25, linear = [0.0 0.0 0.0], angular = [0.0 0.0 0.0])
  push!(ellcenter[1],x[1])
  push!(ellcenter[2],x[2])
  push!(ellcenter[3],x[3])
  push!(ellradii[1],a)
  push!(ellradii[2],b)
  push!(ellradii[3],c)
  push!(ellorient[1],R[1,1])
  push!(ellorient[2],R[2,1])
  push!(ellorient[3],R[3,1])
  push!(ellorient[4],R[1,2])
  push!(ellorient[5],R[2,2])
  push!(ellorient[6],R[3,2])
  push!(ellorient[7],R[1,3])
  push!(ellorient[8],R[2,3])
  push!(ellorient[9],R[3,3])
  push!(elldefo[1],1.0)
  push!(elldefo[2],0.0)
  push!(elldefo[3],0.0)
  push!(elldefo[4],0.0)
  push!(elldefo[5],1.0)
  push!(elldefo[6],0.0)
  push!(elldefo[7],0.0)
  push!(elldefo[8],0.0)
  push!(elldefo[9],1.0)

  (mass, volume, E, J) = ellipsoid_properties(a, b, c, R, density)

  push!(elldensity, density)
  push!(ellyoung, young)
  push!(ellpoisson, poisson)
  push!(ellmass, mass)
  push!(elleuler[1], E[1,1])
  push!(elleuler[2], E[2,1])
  push!(elleuler[3], E[3,1])
  push!(elleuler[4], E[1,2])
  push!(elleuler[5], E[2,2])
  push!(elleuler[6], E[3,2])
  push!(elleuler[7], E[1,3])
  push!(elleuler[8], E[2,3])
  push!(elleuler[9], E[3,3])
  invE = inv(E)
  push!(ellinveuler[1], invE[1,1])
  push!(ellinveuler[2], invE[2,1])
  push!(ellinveuler[3], invE[3,1])
  push!(ellinveuler[4], invE[1,2])
  push!(ellinveuler[5], invE[2,2])
  push!(ellinveuler[6], invE[3,2])
  push!(ellinveuler[7], invE[1,3])
  push!(ellinveuler[8], invE[2,3])
  push!(ellinveuler[9], invE[3,3])
  push!(ellinertia[1], J[1,1])
  push!(ellinertia[2], J[2,1])
  push!(ellinertia[3], J[3,1])
  push!(ellinertia[4], J[1,2])
  push!(ellinertia[5], J[2,2])
  push!(ellinertia[6], J[3,2])
  push!(ellinertia[7], J[1,3])
  push!(ellinertia[8], J[2,3])
  push!(ellinertia[9], J[3,3])
  invJ = inv(J)
  push!(ellinvinertia[1], invJ[1,1])
  push!(ellinvinertia[2], invJ[2,1])
  push!(ellinvinertia[3], invJ[3,1])
  push!(ellinvinertia[4], invJ[1,2])
  push!(ellinvinertia[5], invJ[2,2])
  push!(ellinvinertia[6], invJ[3,2])
  push!(ellinvinertia[7], invJ[1,3])
  push!(ellinvinertia[8], invJ[2,3])
  push!(ellinvinertia[9], invJ[3,3])

  push!(elllinear[1],linear[1])
  push!(elllinear[2],linear[2])
  push!(elllinear[3],linear[3])
  push!(ellangular[1],angular[1])
  push!(ellangular[2],angular[2])
  push!(ellangular[3],angular[3])

  global ellcount
  ellcount += 1
end

# triangles data
tricount = 0
trip = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
triq = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
trir = [Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount), Array{Float32}(undef, ellcount)]
trilinear = Array{Function}(undef, ellcount) # [vx, vy, vz] = linear[i](t)
triangular = Array{Function}(undef, ellcount) # [ox, oy, oz] = angular[i](t)

# insert an obstacle
# triangles format: [[[p1, q1, r1], [p2, q2, r2], [p3, q3, r3]], [[p1, q1, r1], [p2, q2, r2], [p3, q3, r3]], ...]
zerofunc(t) = Float32[0.0; 0.0; 0.0]
function obstacle(triangles, linear = zerofunc, angular = zerofunc)
  for i = 1:size(triangles)[1]
    push!(trip[1], triangles[i][1][1])
    push!(trip[2], triangles[i][1][2])
    push!(trip[3], triangles[i][1][3])
    push!(triq[1], triangles[i][2][1])
    push!(triq[2], triangles[i][2][2])
    push!(triq[3], triangles[i][2][3])
    push!(trir[1], triangles[i][3][1])
    push!(trir[2], triangles[i][3][2])
    push!(trir[3], triangles[i][3][3])
    push!(trilinear, linear)
    push!(triangular, angular)

    global tricount
    tricount += 1
  end
end

# set gravity
grav = [0, 0, 0]
function gravity(g1, g2, g3)
 grav[1] = g1
 grav[2] = g2
 grav[3] = g3
end

# run rigid simulation
function run_rigid(duration, step, outstep, outdir)

  output(0.0, outdir)
end

# run pseudo-rigid simulation
function run_pseudorigid(duration, step, outstep, outdir)
end

export ellipsoid, obstacle, gravity, run_rigid, run_pseudorigid

# --- utility functions ----

# calculate boundingx box
function ellbox()

  amax = maximum(ellradii[1])
  bmax = maximum(ellradii[2])
  cmax = maximum(ellradii[3])

  xlo = minimum(ellcenter[1]) - amax
  ylo = minimum(ellcenter[2]) - bmax
  zlo = minimum(ellcenter[3]) - cmax
  xhi = maximum(ellcenter[1]) + amax
  yhi = maximum(ellcenter[2]) + bmax
  zhi = maximum(ellcenter[3]) + cmax

  return (xlo, ylo, zlo, xhi, yhi, zhi)
end

# output current time step
function output(curtime, outdir)

  mkpath(outdir)

  path = joinpath(outdir,"ellip.dump")

  open(path, "w") do f
    write(f, "ITEM: TIMESTEP\n")
    write(f, curtime)
    write(f, "ITEM: NUMBER OF ATOMS\n")
    write(f, ellcount)
    (x0, y0, z0, x1, y1, z1) = ellbox()
    write(f, "ITEM: BOX BOUNDS\n");
    write(f, "$x0 $x1\n")
    write(f, "$y0 $y1\n")
    write(f, "$z0 $z1\n")
    write(f, "ITEM: ATOMS id x y z radius vx vy vz omegax omegay omegaz\n")
    for i = 1:ellcount
      println(f, i, " ", ellcenter[1][i], " ", ellcenter[2][i], " ", ellcenter[3][i], " ",
                (ellradii[1][i]+ellradii[2][i]+ellradii[3][i])/3.0, " ",
		 elllinear[1][i], " ", elllinear[2][i], " ", elllinear[3][i], " ",
		 ellangular[1][i], " ", ellangular[2][i], " ", ellangular[3][i], " ")
    end
  end
end

# ellipsoid properties
include("ellpro.jl")

end
# <-- edyna Julia module ends

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
@testset "91.decode-ways.jl" begin
    @test num_decodings("12") == 2
    @test num_decodings("226") == 3
    @test num_decodings("0") == 0
    @test num_decodings("06") == 0
end


text:
Base.load_juliarc()

using MATLAB

const libmex = dlopen(MATLAB.matlab_library("libmex"), Libdl.RTLD_GLOBAL | Libdl.RTLD_LAZY)
mexfn(fun::Symbol) = dlsym(libmex::Ptr{Void}, fun)

function call_matlab(fn::ASCIIString, nout::Integer, args...)
  mxs = map(mxarray, args)
  ins = Ptr{Void}[mx.ptr for mx in mxs]
  nin = length(ins)

  outs = Vector{Ptr{Void}}(nout)

  @assert 0 == ccall(mexfn(:mexCallMATLAB), Int32,
    (Int32, Ptr{Ptr{Void}}, Int32, Ptr{Ptr{Void}}, Ptr{Uint8}),
     nout,  outs,           nin,   ins,            fn        )

  [jvariable(MxArray(o)) for o in outs]
end

mex_write(fid, bytes) = call_matlab("fwrite", 0, Float64(fid), bytes)

# redirect STDOUT, STDERR -- wish I knew how to do that!
#const pout = redirect_stdout()[1]
#const perr = redirect_stderr()[1]
# pout.readcb = p -> mex_write(1, readavailable(p))
# perr.readcb = p -> mex_write(2, readavailable(p))
# Base.start_reading(pout)
# Base.start_reading(perr)

function mex_args(ins)
  [ jvariable(MxArray(mx, false)) for mx in ins ]
end

function mex_return(outs, vs...)
  nouts = length(outs)
  @assert nouts == length(vs)
  for i in 1:nouts
    mx = mxarray(vs[i])
    mx.own = false
    outs[i] = mx.ptr
  end
end

function mex_showerror(e)
  buf = IOBuffer()
  showerror(buf, e)
  seek(buf, 0)
  ccall(mexfn(:mexErrMsgTxt), Void, (Ptr{Uint8},), readall(buf))
end

# a fancier eval
function mex_eval(outs::Vector{Ptr{Void}}, ins::Vector{Ptr{Void}})
  try
    @assert length(outs) == length(ins)
    mex_return(outs, [ eval(parse(e)) for e in mex_args(ins) ]...)
  catch e
    mex_showerror(e)
  end
end

# call an arbitrary julia function (or other callable)
function mex_call(outs::Vector{Ptr{Void}}, ins::Vector{Ptr{Void}})
  try
    @assert length(outs) == 1 && length(ins) >= 1
    args = mex_args(ins)
    mex_return(outs, eval(parse(args[1]))(args[2:end]...))
  catch e
    mex_showerror(e)
  end
end

max_stars_repo_name: bermanmaxim/jlcall

max_issues_repo_name: bermanmaxim/jlcall

max_forks_repo_name: bermanmaxim/jlcall


text:
#=

    This is a port of the SPPL model student-interviews.pynb

    Here is the output from the SPPL model of gpa_prior and gpa_posterior for gpa_student[1]:

        gpa                   gpa_prior (cumsum < gpa)  gpa_posterior (cumsum < gpa)
        -----------------------------------------------------------------------------
        0.0                   0.0                    0.0
        0.21052631578947367   2.9321229365490116e-08 3.938526167438616e-08
        0.42105263157894735   3.399923850680484e-06  4.566892092513705e-06
        0.631578947368421     5.235987400911553e-05  7.03315441989632e-05
        0.8421052631578947    0.0003515966799506029  0.0004722764885157954
        1.0526315789473684    0.0014935106089171875  0.0020061336928425886
        1.263157894736842     0.004734542792041406   0.0063595971522462345
        1.4736842105263157    0.012227795632612991   0.01642478644698973
        1.6842105263157894    0.027098290512422067   0.03639933542542442
        1.894736842105263     0.05325350529777253    0.07153189980839067
        2.1052631578947367    0.09485332168120815    0.12741017262720986
        2.3157894736842106    0.15542018995700252    0.20876562761532186
        2.526315789473684     0.23665218641802727    0.31787917797401755
        2.7368421052631575    0.3371086961787463     0.45281574132555596
        2.9473684210526314    0.45107070187375614    0.6058933411522036
        3.1578947368421053    0.5680350963535515     0.7630038505974432
        3.3684210526315788    0.6734850651989842     0.9046477961757059
        3.578947368421052     0.7517864026032522     0.983997710206192
        3.789473684210526     0.7922926319251767     0.9970552758034202
        4.0                   0.8                    0.9995398186972839

        The mean value of the prior is about: 3.050505050505051
        The mean value of the posterior is about: 2.8080808080808084
        (I am not sure how to get an exact value of the mean in SPPL.
         For continous distributions, SPPL requires that one work with
         intervals.)

    This Turing.jl model:

    * Prior()
    Summary Statistics
             parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
                 Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

     student_perfect[1]    0.1982    0.3987     0.0040    0.0036    9428.5188    0.9999    38018.2211
     student_perfect[2]    0.1940    0.3954     0.0040    0.0039   10299.3911    0.9999    41529.8029
         student_gpa[1]    3.0403    0.6872     0.0069    0.0069    8852.1375    0.9999    35694.1029
         student_gpa[2]    3.0270    0.6943     0.0069    0.0070   10029.1311    1.0001    40440.0447
         num_recruiters   24.9399    4.9591     0.0496    0.0502    9860.3856    1.0002    39759.6194
  student_interviews[1]   14.6656    5.4594     0.0546    0.0522    9365.7231    1.0001    37765.0125
  student_interviews[2]   14.5799    5.4370     0.0544    0.0518   10197.5782    1.0000    41119.2671
      student_offers[1]    5.8633    2.8757     0.0288    0.0262    9664.9849    1.0001    38971.7134
      student_offers[2]    5.8544    2.8591     0.0286    0.0278   10359.6432    0.9999    41772.7549

    * With the observation/condition: student_offers[1] == 1 && num_recruiters > 30
      using PG(5)

      Note: I am not happy at all with some of the ess and rhat of this model.

    Summary Statistics
                parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
                    Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

        student_perfect[1]    0.0008    0.0283     0.0003    0.0008   2745.1436    1.0007      114.0626
        student_perfect[2]    0.5617    0.4962     0.0050    0.0493     21.6873    1.0500        0.9011
            student_gpa[1]    2.8155    0.3181     0.0032    0.0312     23.9498    1.0122        0.9951
            student_gpa[2]    3.4779    0.8543     0.0085    0.0847     21.8331    1.0939        0.9072
            num_recruiters   32.0304    1.1016     0.0110    0.1015     24.3081    1.3990        1.0100
     student_interviews[1]   12.7146    0.6928     0.0069    0.0576     34.6829    1.0038        1.4411
     student_interviews[2]   22.6132    3.1529     0.0315    0.3070     23.8147    1.0007        0.9895
         student_offers[1]    1.0187    0.3429     0.0034    0.0187    312.0347    1.0029       12.9653
         student_offers[2]    8.4442    1.9881     0.0199    0.1962     21.7784    1.2301        0.9049
    student_perfect1_prior    0.2032    0.4024     0.0040    0.0051   6213.4045    1.0000      258.1711
        student1_gpa_prior    3.0350    0.6952     0.0070    0.0086   6705.7080    1.0006      278.6267

    Quantiles
                parameters      2.5%     25.0%     50.0%     75.0%     97.5%
                    Symbol   Float64   Float64   Float64   Float64   Float64

        student_perfect[1]    0.0000    0.0000    0.0000    0.0000    0.0000
        student_perfect[2]    0.0000    0.0000    1.0000    1.0000    1.0000
            student_gpa[1]    1.9660    2.6852    2.6852    3.1288    3.1889
            student_gpa[2]    1.6329    3.6500    4.0000    4.0000    4.0000
            num_recruiters   31.0000   31.0000   32.0000   32.0000   34.0000
     student_interviews[1]   11.0000   13.0000   13.0000   13.0000   13.0000
     student_interviews[2]   16.0000   22.0000   25.0000   25.0000   25.0000
         student_offers[1]    1.0000    1.0000    1.0000    1.0000    1.0000
         student_offers[2]    4.0000    6.0000   10.0000   10.0000   10.0000
    student_perfect1_prior    0.0000    0.0000    0.0000    0.0000    1.0000
        student1_gpa_prior    1.6481    2.5444    3.0229    3.5678    4.0000

      Given that student 1 only get one offer (and from quite many interviews) makes it probable
      that he/she get a worse GPA (2.8155) than the prior (about 3.04).

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function student_interviews(num_students=2)
    # Students.
    student_perfect    = tzeros(num_students)
    student_gpa        = tzeros(num_students)
    student_interviews = tzeros(num_students)
    student_offers     = tzeros(num_students)

    # For each student, sample GPA.
    for s in 1:num_students
        student_perfect[s] ~ Bernoulli(0.2)
        if student_perfect[s] == false
            # Note: The SPPL model has (a=7,b=3,scale=4), but this form is not supported in Turing.jl
            #       or Distributions.jl
            # This works and is translated to LocationScale automatically.
            student_gpa[s] ~ Beta(7, 3)*4
        else
            student_gpa[s] ~ Dirac(4.0)
        end
    end

    # Sample number of recruiters
    num_recruiters ~ truncated(Poisson(25),10,40)
    # num_recruiters ~ truncated(Poisson(25),30.001,40)  # Alternative version

    # For each student, sample interviews and offers.
    for j in 1:num_students
        student_interviews[j] ~ (student_gpa[j] == 4)      ? Binomial(num_recruiters, 0.9) :
                                (3.5 < student_gpa[j] < 4) ? Binomial(num_recruiters, 0.6) :
                                                             Binomial(num_recruiters, 0.5)

        student_offers[j] ~ Binomial(student_interviews[j], 0.4)
    end

    # Observation/condition
    true ~ Dirac(student_offers[1] == 1 && num_recruiters > 30)
    # true ~ Dirac(student_offers[1] == 1)
    # true ~ Dirac(num_recruiters > 30)

    # Priors
    student_perfect1_prior ~ Bernoulli(0.2)
    if student_perfect1_prior == false
        student1_gpa_prior ~ Beta(7, 3)*4
    else
        student1_gpa_prior ~ Dirac(4.0)
    end

end

model = student_interviews()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
chns = sample(model, PG(5),  10_000)
# chns = sample(model, PG(5),  MCMCThreads(), 10_000, 4)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)
# display(plot(chns))

chains_params = Turing.MCMCChains.get_sections(chns, :parameters)

# show_var_dist_pct(chns,Symbol("student_offers[1]"))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Zygote
using Zygote: @adjoint, literal_getproperty
import Zygote: accum

import Distances: pairwise, colwise
const dtol = 1e-12 # threshold value for precise recalculation of distances

@nograd MersenneTwister, propertynames, broadcast_shape

function accum(D::Diagonal{T}, B::AbstractMatrix) where {T}
    A = Matrix{Union{T, Nothing}}(undef, size(D))
    A[diagind(A)] .= D.diag
    return accum(A, B)
end
accum(A::AbstractMatrix, D::Diagonal) = accum(D, A)
accum(A::Diagonal, B::Diagonal) = Diagonal(accum(diag(A), diag(B)))

@adjoint function ZygoteRules.literal_getproperty(C::Cholesky, ::Val{:factors})
    error("@adjoint not implemented for :factors as is unsafe.")
    return literal_getproperty(C, Val(:factors)), function(Δ)
        error("@adjoint not implemented for :factors. (I couldn't make it work...)")
    end
end

import LinearAlgebra: HermOrSym, diag, Diagonal

diag(S::Symmetric{T, <:Diagonal{T}} where T) = S.data.diag
Zygote._symmetric_back(Δ::Diagonal) = Δ

# Diagonal matrices are always symmetric...
cholesky(A::HermOrSym{T, <:Diagonal{T}} where T) = cholesky(Diagonal(diag(A)))

#
# Some very specific broadcasting hacks while Zygote has crappy broadcasting.
#

import Base.Broadcast: broadcasted

function rrule(::typeof(broadcasted), ::typeof(-), x::AbstractArray)
    function broadcasted_minus_pullback(Δ)
        return (NO_FIELDS, DoesNotExist(), .-Δ)
    end
    return .-x, broadcasted_minus_pullback
end

function rrule(::typeof(broadcasted), ::typeof(exp), x::AbstractArray)
    y = exp.(x)
    return y, Δ->(NO_FIELDS, DoesNotExist(), Δ .* y)
end


text:
@testset "Test ForecastCache" begin
    sys = IS.SystemData()
    name = "Component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(sys, component)
    file = joinpath(FORECASTS_DIR, "DateTimeAsColumnDeterministic.csv")
    forecast = IS.Deterministic("test", file, component, Dates.Hour(1))
    initial_timestamp = IS.get_initial_timestamp(forecast)
    initial_times = collect(IS.get_initial_times(forecast))
    interval = IS.get_interval(forecast)
    IS.add_time_series!(sys, component, forecast)

    cache = IS.ForecastCache(IS.Deterministic, component, "test")
    @test cache.in_memory_count == 168
    @test IS.get_next_time(cache) == initial_timestamp
    @test length(cache) == cache.common.num_iterations == 168

    # Iterate over all initial times with default cache size.
    cache = IS.ForecastCache(IS.Deterministic, component, "test")
    for (i, ta) in enumerate(cache)
        it = initial_times[i]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    IS.reset!(cache)
    for it in initial_times
        ta = IS.get_next_time_series_array!(cache)
        @test first(TimeSeries.timestamp(ta)) == it
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end
    @test IS.get_next_time(cache) === nothing

    # Iterate over all initial times with custom cache size.
    cache = IS.ForecastCache(IS.Deterministic, component, "test"; cache_size_bytes = 1024)
    @test length(cache) == cache.common.num_iterations == 168
    for (i, ta) in enumerate(cache)
        it = initial_times[i]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    IS.reset!(cache)
    for it in initial_times
        ta = IS.get_next_time_series_array!(cache)
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    # Start at an offset.
    cache = IS.ForecastCache(
        IS.Deterministic,
        component,
        "test";
        start_time = Dates.DateTime("2020-01-02T00:00:00"),
    )
    for (i, ta) in enumerate(cache)
        it = initial_times[i + 24]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    # Test caching internals.
    cache = IS.ForecastCache(IS.Deterministic, component, "test"; cache_size_bytes = 1024)
    @test cache.in_memory_count == 5
    @test IS.get_next_time(cache) == initial_timestamp
    for it in initial_times[1:(cache.in_memory_count)]
        ta = IS.get_next_time_series_array!(cache)
        @test IS._get_last_cached_time(cache) == initial_times[5]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    # The next access should trigger a read.
    ta = IS.get_next_time_series_array!(cache)
    @test IS._get_last_cached_time(cache) == initial_times[10]
    @test TimeSeries.timestamp(ta) ==
          IS.get_time_series_timestamps(component, forecast, initial_times[6])
    @test TimeSeries.values(ta) ==
          IS.get_time_series_values(component, forecast, initial_times[6])
    @test IS.get_next_time(cache) == initial_times[7]
end

@testset "Test StaticTimeSeriesCache" begin
    sys = IS.SystemData()
    name = "Component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(sys, component)

    initial_timestamp = Dates.DateTime("2020-09-01")
    resolution = Dates.Hour(1)

    data = TimeSeries.TimeArray(
        range(initial_timestamp; length = 365, step = resolution),
        rand(365),
    )
    ts = IS.SingleTimeSeries(data = data, name = "test")
    IS.add_time_series!(sys, component, ts)

    cache = IS.StaticTimeSeriesCache(IS.SingleTimeSeries, component, "test")
    @test cache.in_memory_rows == 365
    @test IS.get_next_time(cache) == initial_timestamp
    @test length(cache) == cache.common.num_iterations == 1

    # Iterate over all initial times with default cache size.
    cache = IS.StaticTimeSeriesCache(IS.SingleTimeSeries, component, "test")
    for (i, ta) in enumerate(cache)
        it = initial_timestamp + (i - 1) * resolution
        @test TimeSeries.timestamp(ta) == IS.get_time_series_timestamps(component, ts, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, ts, it)
    end

    ta = IS.get_next_time_series_array!(cache)
    @test first(TimeSeries.timestamp(ta)) == initial_timestamp
    @test TimeSeries.timestamp(ta) ==
          IS.get_time_series_timestamps(component, ts, initial_timestamp)
    @test TimeSeries.values(ta) ==
          IS.get_time_series_values(component, ts, initial_timestamp)

    # Iterate over all initial times with custom cache size.
    cache = IS.StaticTimeSeriesCache(
        IS.SingleTimeSeries,
        component,
        "test";
        cache_size_bytes = 1024,
    )
    @test cache.in_memory_rows == 128
    @test length(cache) == cache.common.num_iterations == 3
    ta = IS.get_next_time_series_array!(cache)
    @test first(TimeSeries.timestamp(ta)) == initial_timestamp
    @test length(ta) == 128
    @test TimeSeries.values(ta) == TimeSeries.values(data)[1:128]

    IS.reset!(cache)
    for (i, ta) in enumerate(cache)
        it = initial_timestamp + (i - 1) * cache.in_memory_rows * resolution
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, ts, it; len = length(ta))
        @test TimeSeries.values(ta) ==
              IS.get_time_series_values(component, ts, it; len = length(ta))
    end

    IS.reset!(cache)
    for i in 1:3
        ta = IS.get_next_time_series_array!(cache)
        it = initial_timestamp + (i - 1) * cache.in_memory_rows * resolution
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, ts, it; len = length(ta))
        @test TimeSeries.values(ta) ==
              IS.get_time_series_values(component, ts, it; len = length(ta))
    end

    cache = IS.StaticTimeSeriesCache(
        IS.SingleTimeSeries,
        component,
        "test";
        start_time = Dates.DateTime("2020-09-06T08:00:00"),
        cache_size_bytes = 1024,
    )
    @test cache.in_memory_rows == 128
    @test cache.common.num_iterations == 2
    @test IS._get_length_remaining(cache) == 365 - 128
    for i in 1:2
        ta = IS.get_next_time_series_array!(cache)
        it = initial_timestamp + i * cache.in_memory_rows * resolution
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, ts, it; len = length(ta))
        @test TimeSeries.values(ta) ==
              IS.get_time_series_values(component, ts, it; len = length(ta))
    end
end

@testset "Test DeterministicSingleTimeSeries with ForecastCache" begin
    sys = IS.SystemData()
    component = IS.TestComponent("Component1", 5)
    IS.add_component!(sys, component)

    resolution = Dates.Minute(5)
    dates = create_dates("2020-01-01T00:00:00", resolution, "2020-01-31T23:00:00")
    data = collect(1:length(dates))
    ta = TimeSeries.TimeArray(dates, data, [IS.get_name(component)])
    name = "val"
    ts = IS.SingleTimeSeries(name, ta)
    IS.add_time_series!(sys, component, ts)
    horizon = 24
    interval = Dates.Hour(1)
    IS.transform_single_time_series!(
        sys,
        IS.DeterministicSingleTimeSeries,
        horizon,
        interval,
    )

    forecast = IS.get_time_series(IS.AbstractDeterministic, component, name)
    initial_times = collect(IS.get_initial_times(forecast))
    cache =
        IS.ForecastCache(IS.AbstractDeterministic, component, name; cache_size_bytes = 1024)

    for (i, ta) in enumerate(cache)
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, initial_times[i])
        @test TimeSeries.values(ta) ==
              IS.get_time_series_values(component, forecast, initial_times[i])
    end
end

@testset "Test Probabilistic with ForecastCache" begin
    initial_time = Dates.DateTime("2020-09-01")
    interval = Dates.Hour(1)
    resolution = Dates.Hour(1)
    name = "test"
    horizon = 24
    data = SortedDict{Dates.DateTime, Matrix{Float64}}()
    for (i, it) in enumerate(range(initial_time, step = interval, length = 100))
        data[it] = ones(horizon, 99) * i
    end
    sys = IS.SystemData()
    component_name = "Component1"
    component = IS.TestComponent(component_name, 5)
    IS.add_component!(sys, component)
    forecast = IS.Probabilistic(name, data, ones(99), resolution)
    IS.add_time_series!(sys, component, forecast)

    # Iterate over all initial times with custom cache size.
    sz = 1024 * 1024
    cache = IS.ForecastCache(IS.Probabilistic, component, "test"; cache_size_bytes = sz)
    initial_times = collect(keys(data))
    @test cache.in_memory_count == trunc(Int, sz / (99 * 8 * 24))
    for (i, ta) in enumerate(cache)
        it = initial_times[i]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end
end


text:
# mask operator for matrix completion
struct MaskOP{T1<:Int64, T2<:Vector{Int64}, T3<:SparseMatrixCSC{Float64, Int64}} <: AbstractOperator
    m::T1
    n::T1
    nnz::T1
    I::T2
    J::T2
    M::T3
    function MaskOP(M::SparseMatrixCSC{Float64, Int64})
        m, n = size(M)
        I, J, ~ = findnz(M)
        nnz = length(I)
        new{Int64, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}(m, n, nnz, I, J, M)
    end
end

size(Mop::MaskOP) = Mop.nnz, Mop.m * Mop.n

struct TMaskOP{T1<:Int64, T2<:Vector{Int64}, T3<:SparseMatrixCSC{Float64, Int64}} <: AbstractOperator
    m::T1
    n::T1
    nnz::T1
    colptr::T2
    rowval::T2
    M::T3
    function TMaskOP(M::SparseMatrixCSC{Float64, Int64})
        m, n = size(M)
        colptr = M.colptr
        rowval = M.rowval
        nnz = length(rowval)
        new{Int64, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}(m, n, nnz, colptr, rowval, M)
    end
end

size(TMop::TMaskOP) = TMop.m * TMop.n, TMop.nnz

Base.adjoint(Mop::MaskOP) = TMaskOP(Mop.M)

function Base.:(*)(Mop::MaskOP, a::NucBallAtom)
    u, v = a.u, a.v
    I, J = Mop.I, Mop.J
    V = []
    for k=1:length(I)
        i, j = I[k], J[k]
        append!(V, u[i]*v[j])
    end
    A = sparse(I, J, V)
    return A.nzval
end

function LinearAlgebra.mul!(Ma::Vector{Float64}, Mop::MaskOP, a::NucBallAtom)
    u, v = a.u, a.v
    I, J = Mop.I, Mop.J
    for k=1:length(I)
        i, j = I[k], J[k]
        Ma[k] = u[i]*v[j]
    end
end

function Base.:(*)(Mop::MaskOP, X::Tuple{Matrix{Float64}, Adjoint{Float64, Matrix{Float64}}})
    U, V = X[1], X[2]
    I, J = Mop.I, Mop.J
    Val = []
    for k=1:length(I)
        i, j = I[k], J[k]
        append!(Val, U[i,:]'*V[j,:])
    end
    A = sparse(I, J, Val)
    return A.nzval
end

function Base.:(*)(TMop::TMaskOP, y::Vector{Float64})
    m, n = TMop.m, TMop.n
    colptr, rowval = TMop.colptr, TMop.rowval
    return SparseMatrixCSC{Float64, Int64}(m, n, colptr, rowval, y)
end

function LinearAlgebra.mul!(z::SparseMatrixCSC{Float64, Int64}, TMop::TMaskOP, y::Vector{Float64})
    z.nzval .= y
end

function Base.:(*)(M::MaskOP, F::NucBallFace)
    r = F.r
    U, V = F.U, F.V
    f = p->M * (U*reshape(p, r, r), V)
    fc = q->vec(F.U'*(M'*q)*V)
    return LinearMap(f, fc, nnz(M.M), r*r)
end

################################################################################
# phaselift operator for phase retrieval
struct PhaseLiftOP{T1<:Int64, T2<:Matrix{Float64}} <: AbstractOperator
    n::T1
    p::T1
    M::T2
    Mt::T2
    function PhaseLiftOP(M::Matrix{Float64})
        n, p = size(M)
        Mt = copy(M')
        new{Int64, Matrix{Float64}}(n, p, M, Mt)
    end
end

size(PLop::PhaseLiftOP) = PLop.p, Plop.n*PLop.n

struct TPhaseLiftOP{T1<:Int64, T2<:Matrix{Float64}} <: AbstractOperator
    n::T1
    p::T1
    M::T2
    Mt::T2
    function TPhaseLiftOP(M::Matrix{Float64})
        n, p = size(M)
        Mt = copy(M')
        new{Int64, Matrix{Float64}}(n, p, M, Mt)
    end
end

size(TPLop::TPhaseLiftOP) = TPlop.n*TPLop.n, TPLop.p

Base.adjoint(PLop::PhaseLiftOP) = TPhaseLiftOP(PLop.M)

function Base.:(*)(PLop::PhaseLiftOP, a::TraceBallAtom)
    y = PLop.Mt * a.v
    y .^= 2
    return y
end

function LinearAlgebra.mul!(Ma::Vector{Float64}, PLop::PhaseLiftOP, a::TraceBallAtom)
    Ma .= PLop.Mt * a.v
    Ma .^= 2
    return nothing
end

function Base.:(*)(TPLop::TPhaseLiftOP, y::Vector{Float64})
    L = LinearMap(v->TPLop.M*Diagonal(y)*TPLop.Mt*v, TPLop.n; issymmetric=true, isposdef=true)
    return L
end

function LinearAlgebra.mul!(z::LinearMap, TPLop::TPhaseLiftOP, y::Vector{Float64})
    z = TPLop*y
    return nothing
end

function Base.:(*)(PLop::PhaseLiftOP, F::TraceBallFace)
    n, p = size(PLop.M)
    r = size(F.U, 2)
    MtU = PLop.Mt * F.U
    function f(w::Vector{Float64})
        W = reshape(w, r, r)
        MtUW = MtU*W
        v = zeros(p)
        for i in 1:p
            v[i] = MtUW[i,:]'*MtU[i,:]
        end
        return v
    end
    f = p->M * (U*reshape(p, r, r), V)
    fc = q->vec(F.U'*(M'*q)*V)
    return LinearMap(f, fc, nnz(M.M), r*r)
end


text:
##Reineke SDI chart
#using Gadfly
using Plots; gr()

function sdi_chart(tpa,qmd;max_sdi=450)
sdi=max_sdi
diarng=collect(1:20)
tparng=[1,1200]
maxline=sdi./(((diarng)./10.0).^1.605)
comp_mortline=maxline*0.55
crown_closeline=maxline*0.35

myxticks=[1,5,10,50,100,500,1000]
myyticks=[1,2,5,10,20]

plot(crown_closeline,diarng,color=:green)
plot!(comp_mortline,diarng,color=:orange)
plot!(maxline,diarng,color=:red,
xscale=:log10,
yscale=:log10,
xlabel="Trees Per Acre",
ylabel="Quadratic Mean Diameter",
legend=false)
annotate!([(80, 18, text("35%",10)),
           (160, 18, text("55%",10)),
           (320, 18, text("100%",10))
           ])
scatter!((tpa,qmd))
end #end function


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

using DifferentialEquations
using LinearAlgebra
using Plots
using RecursiveArrayTools

function main()
# Set the required model parameters for the SIRS model with two levels of transmission - Within and between households
N = 10; # Household size - Change to 10 for final analysis
betaHH = 6; # Within household transmission parameter
betaG = 1; # Population wide transmission
gamma = 1; # Rate of recovery from infection
tau = 1; # Rate of loss of protection
params = [betaHH,gamma,tau,betaG,N]; # Put all the parameters together
time = (0.0, 30.0) # Simulation time - note it defined as a float
dim = dim = 0.5*(N+1)*(N+2); # Number of possible configurations - works for three epidemiological classes
y0 = zeros(Float64, 1, Int(dim)); # Initial condition vector
y0[end-1] = 0.00000001;
y0[end] = 0.99999999;

function hhTransitions(N,dim)
    # Function to generate transition matrices for household model
    # Input: N is the household size

    # Initialize things
    Qinf = zeros(Float64, Int(dim), Int(dim));
    Qrec = zeros(Float64, Int(dim), Int(dim));
    Qext = zeros(Float64, Int(dim), Int(dim));
    Qwane = zeros(Float64, Int(dim), Int(dim));
    dataI = zeros(Int64, Int(dim), 3)
    m = 0;
    I = zeros(Int64, Int(N+1), Int(N+1))

    # To help remember where to store the variables
    for ss = 0:N
        for ii = 0:(N-ss)
            m = m + 1;
            I[ss+1,ii+1] = m
        end
    end

    # Describe the epidemiological transitions

    # Counter for susceptibles
    for ss = 0:N
        # Counter for infecteds
        for ii = 0:(N-ss)
            # If susceptibles and infecteds are more than 1, then infection within the household can occur
            if (ss > 0 && ii > 0)
                Qinf[I[ss+1,ii+1],I[ss,ii+2]] = ii*ss/(N-1);
            end

            # If infecteds are more than 1, recovery can occur
            if ii > 0
                # Rate of recovery
                Qrec[I[ss+1,ii+1],I[ss+1,ii]] = ii;
            end

            # For external infection - just keep track of susceptibles
            if ss > 0
                # Rate of within household infection
                Qext[I[ss+1,ii+1],I[ss,ii+2]] = ss;
            end

            # Loss of protection hence becoming susceptible again. Possible if N-ss-ii = rr > 0
            if (N-ss-ii) > 0
                # Rate of loss of protection
                Qwane[I[ss+1,ii+1],I[ss+2,ii+1]] = N-ss-ii;
            end

            # Store the relevant indices to help identify the household configurations
            dataI[I[ss+1,ii+1],:] = [ss, ii, N-ss-ii];
        end
    end

    Qinf = Qinf - diagm(0=>vec(sum(Qinf,dims= 2)));
    Qrec = Qrec - diagm(0=>vec(sum(Qrec,dims = 2)));
    Qext = Qext - diagm(0=>vec(sum(Qext,dims=2)));
    Qwane = Qwane - diagm(0=>vec(sum(Qwane,dims=2)));

    # Return
    return Qinf, Qrec, Qext, Qwane, dataI
end
Qinf, Qrec, Qext, Qwane, dataI = hhTransitions(N,dim);

# +
function rateSIRS(dy_dt,y0,params,time)

    # Extract the parameters
    betaHH = params[1];
    gamma = params[2];
    tau = params[3];
    betaG = params[4];
    N = params[5];

    # Generate the transition matrices
    Qinf, Qrec, Qext, Qwane, HHconfig = hhTransitions(N,dim);

    # Combine within and external transitions
    Q = betaHH*Qinf + gamma*Qrec + tau*Qwane + (betaG*((HHconfig[:,2]*y0)/N)*Qext);

end

# Define the ODE problem
prob = ODEProblem(rateSIRS,y0,time,params);

# Solve
sol = solve(prob);
# -

# +
Iconfig = dataI[:,2];
infProp = zeros(Float64,length(sol),1);

# Prepare the plots
for i = 1:length(sol.t)
    infProp[i,1] = *Iconfig/N;
end
# -

# Total infectious in the population
plot(sol.t,infProp,color="blue",xlabel="Time",ylabel="Proportion infectious",label=["Mean infection"],ylims=[0, 1])

# Household profile at endemic prevalence
# Prepare the plots
using PyPlot
yprop = zeros(N+1,length(sol.t))
for j = 1:length(sol.t)
    for i = 1:N+1
        index = find(Iconfig.==i-1);
        yprop[i,j] = sum(sol[index,j])
    end
end
step(-0.5:1:10.5,[yprop[:,length(sol.t)];0])

# Household profile at peak prevalence
# Prepare the plots
x = find(infProp.==maximum(infProp))
step(-0.5:1:10.5,[yprop[:,x];0]);

end