text:
# -*- coding: utf-8 -*-
using LinearAlgebra

using Plots
using Distributions

using GeometryObjects

Base.@kwdef struct CollisionSystem
    xmin::Float64 = -4
    xmax::Float64 = 4
    ymin::Float64 = -4
    ymax::Float64 = 4
    maxT::Int = 500
    num_particles::Int = 20
end

function setup_particle(sys::CollisionSystem)
    particles = Particle2D[]
    num_particles = sys.num_particles
    while length(particles) != num_particles
        θ = 2π * rand()
        v_abs = 0.1
        v = Velocity2D(v_abs * cos(θ), v_abs * sin(θ))
        m = r = rand(Uniform(0.1, 0.3))
        x = rand(Uniform(sys.xmin + r, sys.xmax - r))
        y = rand(Uniform(sys.ymin + r, sys.ymax - r))

        c = Particle2D(
            Circle2D(Point2D(x, y), r),
            v,
            m,
        )

        c_plus_r = Particle2D(
            Circle2D(Point2D(x, y), r + v_abs),
            v,
            m,
        )

        any(has_contact.(Ref(c_plus_r), particles)) && continue
        push!(particles, c)
    end
    return particles
end

function move!(sys::CollisionSystem, pt::Point2D, v::Velocity2D, r::Real=0, Δt=1.0)
    pt .+= Δt .* v # do not use `pt += v`
    if pt.y ≥ sys.ymax - r
        pt.y = 2(sys.ymax - r) - pt.y
        v.y = -v.y
    end
    if pt.y ≤ (sys.ymin + r)
        pt.y = 2(sys.ymin + r) - pt.y
        v.y = -v.y
    end

    if pt.x ≥ (sys.xmax - r)
        pt.x = 2(sys.xmax - r) - pt.x
        v.x = -v.x
    end
    if pt.x ≤ (sys.xmin + r)
        pt.x = 2(sys.xmin + r) - pt.x
        v.x = -v.x
    end
    return pt, v
end

move!(sys::CollisionSystem, p::Particle2D; Δt=1.0) = move!(sys, p.c.pt, p.v, p.c.r, Δt)

dist(pt1, pt2) = sqrt((pt1.x - pt2.x)^2 + (pt1.y - pt2.y)^2)
has_contact(c1::Circle2D, c2::Circle2D) = c1.r + c2.r ≥ dist(c1.pt, c2.pt)
has_contact(p1::Particle2D, p2::Particle2D) = has_contact(p1.c, p2.c)

isijulia() = isdefined(Main, :IJulia)

function create_animation(sys::CollisionSystem; disp::Bool=true)
    particles = setup_particle(sys)
    Δt = 0.1
    accumΔt = 0.0
    cnt = 0
    plt = [
        plot(
            size=(400, 300),
            xlim=(sys.xmin, sys.xmax), ylim=(sys.ymin, sys.ymax),
            aspect_ratio=:equal, legend=false,
            grid=false, ticks=false,
            framestyle=:box,
        ),
    ]
    while cnt < sys.maxT
        accumΔt += Δt
        disp && isijulia() && IJulia.clear_output(true)
        move!.(Ref(sys), particles; Δt)
        # naive collision algorithm
        for i in 1:length(particles)
            for j in (i+1):length(particles)
                particle1 = particles[i]
                particle2 = particles[j]
                c1 = particle1.c
                c2 = particle2.c
                # update c1.v and c2.v
                if has_contact(c1, c2)
                    pt1 = c1.pt
                    pt2 = c2.pt

                    v1 = particle1.v
                    v2 = particle2.v
                    m1 = particle1.m
                    m2 = particle2.m

                    # correct position so that c1 and c2 does not contact
                    while has_contact(c1, c2)
                        c1.pt .-= 0.015 .* v1
                        c2.pt .-= 0.015 .* v2
                    end

                    # update velocity
                    v1_next =
                        v1 -
                        2m2 / (m1 + m2) * dot(v1 - v2, pt1 - pt2) * (pt1 - pt2) /
                        dot(pt1 - pt2, pt1 - pt2)
                    v2_next =
                        v2 -
                        2m1 / (m2 + m1) * dot(v2 - v1, pt2 - pt1) * (pt2 - pt1) /
                        dot(pt2 - pt1, pt2 - pt1)
                    particle1.v .= v1_next
                    particle2.v .= v2_next
                end
            end
        end

        if accumΔt ≥ 1.0
            cnt += 1
            s = plot(
                size=(400, 300),
                xlim=(sys.xmin, sys.xmax), ylim=(sys.ymin, sys.ymax),
                aspect_ratio=:equal, legend=false,
                grid=false, ticks=false,
                framestyle=:box,
            )
            plot!(s, particles)
            plot!(s, title="$cnt")
            accumΔt = zero(accumΔt)
            disp && isijulia() && display(s)
            push!(plt, s)
        end
    end

    anim = @animate for p in plt[2:end]
        plot(p)
    end
    return anim
end

anim = create_animation(CollisionSystem(; maxT=500, xmin=-8, xmax=8), disp=isijulia())
gif(anim, "result.gif")


text:
export normcols, normrows, rank, qrfactors, isorth
export tolabsdefault, tolreldefault, adjustsize
export QRFactorization, QRWrapper

const QRFactorization = Union{QRPivoted,
                              LinearAlgebra.QRCompactWY,
                              SuiteSparse.SPQR.QRSparse}

struct QRWrapper{F<:Union{QRPivoted,LinearAlgebra.QRCompactWY,SuiteSparse.SPQR.QRSparse}}
    parent::F
end

"""
    calculate 2-norms all columns of a sparse matrix
"""
function normcols(A::SparseMatrixCSC{T}) where T<:Number
    i = 1
    TT = typeof(sqrt(one(eltype(A))))
    z = zero(T)
    sum = z
    r = Vector{TT}(undef, A.n)
    for j = 1:A.n
        k = A.colptr[j+1]
        while i < k
            sum += abs2(A.nzval[i])
            i += 1
        end
        r[j] = sum
        sum = z
        end
        sqrt.(real(r))
end

function normcols(A::AbstractMatrix{T}) where T<:Number
    [norm(view(A, :, i)) for i in 1:size(A, 2)]
end

"""
    calculate 2-norms of all rows of a sparse matrix
"""
function normrows(A::SparseMatrixCSC{T}) where T<:Number
    TT = typeof(sqrt(one(eltype(A))))
    r = zeros(TT, size(A, 1))
    for i = 1:nnz(A)
        j = A.rowval[i]
        r[j] += abs2(A.nzval[i])
    end
    sqrt.(r)
end

function normrows(A::AbstractMatrix{T}) where T<:Number
    [norm(view(A, i, :)) for i in 1:size(A, 1)]
end

"""
    default tolerance absolute for qr-calculations
"""
function tolabsdefault(A::AbstractMatrix{T}) where T<:Number
    20maximum(normcols(A)) * eps(real(T)) * max(size(A)...)
end

"""
    default tolerance relative for rank-calculations
"""
function tolreldefault(A)
    T = eltype(A)
    sqrt(size(A,1)) * eps(real(T)) * 20
end

"""
    absolute tolerance as a combination of tolrel and tolabs
"""
function tolerance(A::AbstractMatrix, pivot::Bool,
                   tolrel::AbstractFloat, tolabs::AbstractFloat)

    tolabs >= 0 || throw(ArgumentError("absolute tolerance must not be negative"))
    tolrel >= 0 || throw(ArgumentError("relative tolerance must not be negative"))
    if tolrel == 0 == tolabs && pivot
        tolrel = tolreldefault(A)
    end
    tolrel == 0 ? tolabs : max(tolrel * maximum(normcols(A)), tolabs)
end

"""
    rank of a QR-factorization of a matrix
"""
function LinearAlgebra.rank(F::Union{QRFactorization,QRWrapper};
                            tolrel::AbstractFloat=0.0, tolabs::AbstractFloat=0.0)

    FR = F.R
    tol = tolerance(FR, true, tolrel, tolabs)
    count(abs2.(diag(FR)) .> tol^2)
end

"""
    `adjustsize(A::Matrix, m:Integer, n::Integer)`

Return a Matrix with `m` rows and `n` columns. If `m > size(A, 1)` new zero rows are added,
if `n > size(A, 2)` zero columns are added. The left upper corner of the resulting
matrix is filled with the corresponding elements of `A`.
"""
function adjustsize(A::AbstractMatrix, m::Integer=size(A,1), n::Integer=size(A,2))

    m0, n0 = size(A)
    if n == n0 && m == m0
        A
    elseif A isa SparseMatrixCSC &&
        ( n0 <= n && ( m > m0 || maximum(view(A.rowval, 1:nnz(A))) <= m < m0 ))
        if n == n0
            acolptr = A.colptr
        else
            acolptr = resize!(copy(A.colptr), n+1)
            fill!(view(acolptr, n0+2:n+1), acolptr[n0+1])
        end
        SparseMatrixCSC(m, n, acolptr, A.rowval, A.nzval)
    elseif m <= m0 && n <= n0
        A[1:m,1:n]
    else
        B = similar(A, m, n)
        fill!(B, 0)
        m0 = min(m, m0)
        n0 = min(n, n0)
        B[1:m0,1:n0] = A[1:m0,1:n0]
        B
    end
end

Base.show(io::IO, qr::QRWrapper) = show(io, qr.parent)
Base.propertynames(qr::QRWrapper, private::Bool) = (:R, :Q, :pcol, :prow)
import Base.getproperty

function getproperty(qr::QRWrapper{<:SuiteSparse.SPQR.QRSparse}, d::Symbol)
    if d == :R || d == :Q || d == :pcol || d == :prow
        getproperty(qr.parent, d)
    elseif d == :p
        qr.parent.pcol
    else
        getfield(qr, d)
    end
end
function getproperty(qr::QRWrapper{<:QRPivoted}, d::Symbol)
    if d == :R || d == :Q || d == :p || d == :P
        getproperty(qr.parent, d)
    elseif d == :pcol
        getproperty(qr.parent, :p)
    elseif d == :prow
        collect(1:size(qr.parent.Q, 1))
    else
        getfield(qr, d)
    end
end
function getproperty(qr::QRWrapper{<:LinearAlgebra.QRCompactWY}, d::Symbol)
    if d == :R || d == :Q
        getproperty(qr.parent, d)
    elseif d == :pcol
        collect(1:size(qr.parent.R, 2))
    elseif d == :prow
        collect(1:size(qr.parent.Q, 1))
    else
        getfield(qr, d)
    end
end

Base.size(qrw::QRWrapper, arg...) = size(qrw.parent, arg...)

"""
    qrfactors(A::AbstractMatrix; tolrel, tolabs, pivot::Bool=true)

Perform QR-factorization of matrix `A` with various matrix types (dense and sparse).
Allow to specify if pivoting is done or not and a relative - and absolute tolerance
bounds on the diagonal elements of R.
Return a common abstract result type `q::QRWrapper` in all cases, which allows to
access the elements of the factorization `q.Q q.R q.pcol q.prow` in all cases.
In the case of dense `A` with `pivot=true` and `tolerance != 0`, the tiny elements
in `R` are set to zero and the corresponding Householder-factors to unity.
"""
function qrfactors(A::SparseMatrixCSC{<:Union{ComplexF64,Float64}};
                 tolrel::AbstractFloat=0.0, tolabs::AbstractFloat=0.0, pivot=true)

    tol = tolerance(A, pivot, tolrel, tolabs)
    QRWrapper(qrfact(A, tol=tol))
end

function qrfactors(A::AbstractSparseMatrix;
                 tolrel::AbstractFloat=0.0, tolabs::AbstractFloat=0.0, pivot=true)

    throw(ArgumentError("qrfact does not support $(typeof(A))"))
end

function qrfactors(A::AbstractMatrix;
                 tolrel::AbstractFloat=0.0, tolabs::AbstractFloat=0.0, pivot=true)

    tol = pivot ? tolerance(A, pivot, tolrel, tolabs) : 0.0
    if tol > 0 && !pivot
        throw(ArgumentError("tol > 0 requires pivot=true"))
    end
    qr = qrfact(A, Val(pivot))
    if tol > 0
        m = size(qr.factors, 1)
        k = rank(qr, tolrel=tolrel, tolabs=tolabs)
        for i = k+1:length(qr.τ)
            qr.τ[i] = 0
            qr.factors[i:m,i] = 0
        end
    end
    QRWrapper(qr)
end

function qrfactors(A::Union{Adjoint, Transpose};
                 tolrel::AbstractFloat=0.0, tolabs::AbstractFloat=0.0, pivot=true)

    qrfactors(copy(A), tolrel=tolrel, tolabs=tolabs)
end

"""
    isorth(A) -> Bool

Check if matrix `A` is orthogonal/unitary.
If `A` is not square, check if all columns/rows are mutually orthogonal.
"""
function isorth(A::AbstractMatrix)
    m, n = size(A)
    tol = eps(eltype(A)) * max(size(A)...) * maximum(abs.(collect(extrema(A))))
    if m >= n
        maximum(abs.(A'A - I)) < tol
    else
        maximum(abs.(A*A' - I)) < tol
    end
end

ext: jl


text:
using DimensionalData, Test, Unitful
using DimensionalData: identify

struct Unsortable
    val::Int
end

@testset "identify IndexMode" begin
   @testset "identify Categorical from AutoMode" begin
        @test identify(AutoMode(), X, [:a, :b]) == Categorical(Unordered())
        @test identify(AutoMode(), X, ["a", "b"]) == Categorical(Unordered())
        @test identify(AutoMode(), X, ['b', 'a']) == Categorical(Unordered())
        @test identify(AutoMode(Ordered(index=ReverseIndex())), X, ['b', 'a']) ==
            Categorical(Ordered(index=ReverseIndex()))
        # Mixed types are categorical
        @test identify(AutoMode(Ordered(index=ReverseIndex())), X, ['b', 2]) ==
            Categorical(Unordered())
        @test identify(AutoMode(Ordered(index=ReverseIndex())), X, ['b', 2]) ==
            Categorical(Unordered())
    end

    @testset "identify Categorical order" begin
        @test identify(Categorical(Ordered()), X, [1, 2]) == Categorical(Ordered())
        @test identify(Categorical(), X, [1, 2]) == Categorical(Unordered())
    end

    @testset "identify Sampled Order, Span and Sampling from AutoMode" begin
        @testset "identify vectors" begin
            @test identify(AutoMode(), X, [1, 2, 3, 4, 5]) ==
                Sampled(Ordered(), Irregular(), Points())
            @test identify(AutoMode(), X, [5, 4, 3, 2, 1]) ==
                Sampled(Ordered(index=ReverseIndex()), Irregular(), Points())
            @test identify(AutoMode(), X, [500, 3, 7, 99, 1]) ==
                Sampled(Unordered(), Irregular(), Points())
        end
        @testset "identify range" begin
            @test identify(AutoMode(), X, 1:2:10) == Sampled(Ordered(), Regular(2), Points())
            @test identify(AutoMode(), X, 10:-2:1) ==
                Sampled(Ordered(ReverseIndex(), ForwardArray(), ForwardRelation()),
                        Regular(-2), Points())
        end
    end

    @testset "identify Sampled" begin
        @testset "identify Locus" begin
            @test identify(Sampled(sampling=Intervals()), X, 1:2:9) ==
                Sampled(Ordered(), Regular(2), Intervals(Center()))
            @test identify(Sampled(span=Regular(), sampling=Intervals()), Ti, 1:2:9) ==
                Sampled(Ordered(), Regular(2), Intervals(Start()))
        end
        @testset "identify Regular span step" begin
            @test identify(Regular(), X, 1:2:9) == Regular(2)
            @test identify(Sampled(span=Regular()), X, 1:2:9) ==
                Sampled(Ordered(), Regular(2), Points())
            @test identify(Sampled(span=Regular()), X, 9:-2:1) ==
                Sampled(Ordered(index=ReverseIndex()), Regular(-2), Points())
            @test identify(Sampled(span=Regular(1)), X, [1, 2, 3]) ==
                Sampled(Ordered(), Regular(1), Points())
            @test identify(Sampled(span=Irregular()), X, [1, 2, 3]) ==
                Sampled(Ordered(), Irregular(nothing, nothing), Points())
            @test_throws ArgumentError identify(Regular(), X, Val((1, 2, 3)))
            @test_throws ArgumentError identify(Regular(2.0), X, 1:3)
            @test_throws ArgumentError identify(Sampled(span=Regular()), X, [1, 2, 3])
            @test_throws ArgumentError identify(Sampled(span=Regular(2)), X, 1:4)
        end
        @testset "identify Irregular span step" begin
            # TODO clarify this. For `Points` the bounds in
            # Irregular aren't used so `nothing` is a reasonable value.
            # for `Intervals` they are used, and will return `(nothing, nothing)`
            # from `bounds`. After slicing the `bounds` will be correct, so
            # it may be fine to leave this behaviour up to the user.
            @test identify(Sampled(span=Irregular()), X, [2, 4, 8]) ==
                Sampled(Ordered(), Irregular((nothing, nothing)), Points())
        end
    end

    @testset "An unsortable index is Unordered" begin
        @test identify(AutoMode(), X, [Unsortable(1), Unsortable(2)]) ==
            Sampled(Unordered(), Irregular(), Points())
    end

end


text:
## read_nctiles function
#
#examples:
#  fileName="nctiles_grid/GRID"
#  Depth=read_nctiles(fileName,"Depth")
#  hFacC=read_nctiles(fileName,"hFacC")

"""
    read_nctiles(fileName,fldName)

Read model output from Netcdf / NCTiles file and convert to gcmfaces structure.
"""
function read_nctiles(fileName,fldName)

if ~(MeshArrays.grTopo=="llc");
  error("non-llc cases not implemented yet");
end;

fileIn=@sprintf("%s.%04d.nc",fileName,1);
x = ncread(fileIn,fldName);
ndims=length(size(x));

#initialize f
if ndims==2;
  f0=Array{Float64}(undef,90,0);
  f00=Array{Float64}(undef,0,90);
else;
  nr=size(x,3)
  f0=Array{Float64}(undef,90,0,nr);
  f00=Array{Float64}(undef,0,90,nr);
end;
f=[f0,f0,NaN,f00,f00];

#fill in f
for ff=1:13;
  #read one tile
  fileIn=@sprintf("%s.%04d.nc",fileName,ff);
  x = ncread(fileIn,fldName);
  #combine tiles
  if ff<=3;
    f[1]=cat(f[1],x;dims=2);
  elseif ff<=6;
    f[2]=cat(f[2],x;dims=2);
  elseif ff==7;
    f[3]=x;
  elseif ff<=10;
    f[4]=cat(f[4],x;dims=1);
  elseif ff<=13;
    f[5]=cat(f[5],x;dims=1);
  end;

end;

#return f gcmfaces object fld
fld=gcmfaces(5,"llc",f);

end


text:
using Mustache
using DataFrames
tpl = read("grid.tpl", String)
dic = Dict()
dic["offset"] = 1
dic["nx"] = 9
dic["ny"] = 80
dic["nz"] = 5
dic["dxv"] = [9.518768370414794 2.384362757059747 0.5972606471782251 0.14960822534723633 0.1 0.14960822534723633 0.5972606471782251 2.384362757059747 9.518768370414794]
dic["dyv"] = [40.0 60.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 40.0 60.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 40.0 60.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0]
dic["dzv"] = string.([50, 50, 100, 50, 50])
dic["tops"] = 9640.0

#view = dic
vview = Dict()
vview["grid"] = dic

open("grid.ixf", "w") do io
    render(io, tpl, vview)
end

#tpl = "{{#:vec}}{{.}}{{^.[end]}}, {{/.[end]}}{{/:vec}}"
#Ans = render(tpl, vec = ["A1", "B2", "C3"])  # "A1, B2, C3"

I1 = [1, 1, 1]
I2 = [9, 9, 9]
J1 = [80, 22, 51]
df = DataFrame(I1 = Int64[], I2=Int64[])
push!(df, [1, 2])
df = DataFrame([Int64, Int64, Int64, Int64, Int64, Int64, String, String, Int64],
["I1", "I2", "J1", "J2", "K1", "K2", "Grid", "Property", "Expression"])
push!(df, [1, 9, 1, 80, 1, 5, "CoarseGrid", "SAT_TAB", 1])
push!(df, [1, 9, 1, 22, 1, 5, "CoarseGrid", "SAT_TAB", 2])

dic = Dict();
dic["rock_region"] = df
tpl = read("rock.tpl", String)
open("rock.ixf", "w") do io
    render(io, tpl, dic)
end


text:
# """
#   pt_sol = get_knitro_sol(filename)
#
# Compute a local optima of the filename complex QCQP.
# """
# function get_knitro_sol(filename::String)
#   open(joinpath("src_ampl", "minlp.dat"), "w") do f
#     println(f, "param: KEYS: LEFT RIGHT := include \"$(joinpath("..", "$filename"))\";")
#   end
#
#   cd("src_ampl")
#
#   instancename = split(filename, "\\")[end]
#   instancename = split(instancename, ".")[1]
#   out_log = "Knitro_$(instancename).log"
#   run(`cmd /c ampl minlp.run '>' $out_log`)
#   cd("..")
#
#   sol_cplx = read_Knitro_output("Knitro_sol.txt")
#   return sol_cplx
# end

export run_knitro

"""
sol = read_Knitro_output(filepath)

Return the point corresponding to Knitro output stored in the filepath text
file.
"""
function read_Knitro_output(filepath::String)
  filepath = joinpath("src_ampl", "Knitro_sol.txt")
  data = readdlm(filepath)

  sol_cplx = Point()

  for i=1:size(data, 1)
    sol_cplx[Variable(data[i, 1], Complex)] = data[i, 2] + im*data[i, 3]
  end

  return sol_cplx
end

"""
  run_knitro(pb_path::String, src_ampl_path::String)

Run knitro on the `real_minlp.run`, `real_minlp.run` template scripts from the
`src_ampl_path` folder and `pb_path` files ("real_minlp_instance.dat" and
"real_minlp_precond_cstrs.dat").
"""
function run_knitro(pb_path::String; src_ampl_path::String=joinpath(Pkg.dir("MathProgComplex"), "src", "export_dat"))
    root = pwd()
    date = Dates.format(now(), "yy_u_dd_HH_MM_SS")
    outlog = "Knitro_$(date).log"

    cd(pb_path)

    open("real_minlp.run", "w") do f
      println(f, "include $(joinpath(src_ampl_path, "real_minlp.run"));")
    end

    open("real_minlp.mod", "w") do f
      println(f, "include $(joinpath(src_ampl_path, "real_minlp.mod"));")
    end

    try
        # run(`cmd /c ampl real_minlp.run '>' $(outlog)`)
        run(`cmd /c ampl real_minlp.run `)
    catch
        warn(LOGGER, "AMPL/Knitro failed, returning.")
    end

    cd(root)
end

"""
pt_knitro, pt_GOC = read_Knitro_output(pb_path, pb)

Read knitro output files at `pb_path` and return points according to variables
types in pb.
"""
function read_Knitro_output(pb_path::String, pb::Problem)
  function build_pt(vararray, valarray)
    pt = Point()
    for i=1:length(vararray)
      varname, vartype = vararray[i], get_variabletype(pb, String(vararray[i]))
      val = valarray[i]
      if (vartype == Bool)
        val = round(val)
      end
      pt[Variable(varname, vartype)] = val
    end
    return pt
  end

  root = pwd()
  cd(pb_path)

  ## Read Knitro solution
  files = filter(x->ismatch(r"solution\S+csv", x), readdir(pwd()))
  length(files) == 1 || warn(LOGGER, "get_knitro_solutions(): $(length(files)) .csv files found in $(pb_path).\nExtracting solution from $(files[1]).")

  sol_data = readdlm(files[1], ';')

  pt_knitro = build_pt(sol_data[:, 1], sol_data[:, 2])
  pt_GOC = build_pt(sol_data[:, 1], sol_data[:, 3])

  cd(root)
  return pt_knitro, pt_GOC
end

function read_knitro_info_csvfile(filepath::String)
  lines = readdlm(joinpath(filepath, "knitro_info.csv"), ';')
  solve_result_1 = lines[2,2]
  opterror1 = lines[3,2]
  solve_result_2 = lines[5,2]
  opterror2 = lines[6,2]
  solve_result_3 = lines[8,2]
  opterror3 = lines[9,2]

  return solve_result_1, opterror1, solve_result_2, opterror2, solve_result_3, opterror3
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function arrow_test(fname)
  pagewidth, pageheight = 2000, 2000
  Drawing(pagewidth, pageheight, fname)
  origin() # move 0/0 to center
  background("ivory")
  setopacity(0.8)
  setline(2)
  pagetiles = Tiler(pagewidth, pageheight, 40, 8, margin=50)
  for (pos, n) in pagetiles
    setlinecap()
    setlinejoin()
    randomhue()
    gsave()
    arrow(pos, Point(pos.x + 200, pos.y),
      decorate = () -> begin
        setline(0.2)
        circle(O, 5, :stroke)
      end,
      decoration = range(0, stop=1, length=rand(2:10)),
      linewidth=rand(1:20),
      arrowheadlength=rand(1:30),
      arrowheadangle=rand(pi/20:pi/12:pi/3))
    # same length as a non-arrow line?
    line(Point(pos.x, pos.y + 10), Point(pos.x + 200, pos.y+10), :stroke)
    grestore()
  end
  @test finish() == true

end

Random.seed!(42)
fname = "arrow-line-test.pdf"
arrow_test(fname)
println("...finished test: output in $(fname)")


text:
function firefox(p; fname="tmp.html", fdir=normpath(homedir(), ".tmp"))
    # mkpath(dir)
    filepath = normpath(fdir, fname)
    savefig(p, filepath)
    run(`firefox $(filepath)`, wait=false)
    return p
end

function save_results(fdir, ps, data, chain)
    ignores = [k for (k,v) in pairs(ps) if ((v isa String) && isempty(v))]
    mkpath(fdir)

    @info "Saving at: $fdir"
    let
        dic = Dict( zip( keys(ps), values(ps) ) )
        fname = normpath(fdir, savename("PARAMS", ps, "csv"; ignores) )
        safesave( fname, DataFrame( dic ) )
        bson( normpath(fdir, "params.bson") ,  dic )
	end
	let
        fname = normpath( fdir, savename("DATA", ps, "bson"; ignores) )
        bson( fname, Dict("data"=>data) )
	end
    fname = normpath( fdir, savename("CHAIN", ps, "jls"; ignores) )
    safesave( fname, chain )
    return fname
end

# ============================================================================
# post-processing

struct PostProcessing3{F,P,I,C,D,M}
	fdir::F
	ps::P
	ignores::I
	chain::C
	data::D
	model::M
end

function PostProcessing3(fdir, ps, ignores, exclude, fname)
	chain = let
	    chain = read(fname, Chains)
	    not(chain, exclude)
	end
	data = let
		fn_data  = replace(fname, "CHAIN_"=>"DATA_")
		fn_data  = replace(fn_data, ".jls"=>".bson")
		BSON.load(fn_data)["data"]
	end

	model = National.name2model[ps.model]
	PostProcessing3(fdir, ps, ignores, chain, data, model)
end

function not(c::Chains, exclude::AbstractVector=[])
    isempty(exclude) && return c
    n = size(c,3)
    s = filter(x->x ∉ exclude, 1:n)
    return c[:,:,s]
end

function parse_fname(fname; warmup = nothing)
	ignores = []
	fdir_raw, parsed_args, _ = parse_savename(fname, parsetypes = (Int, Float64, Bool))
	!isnothing(warmup) && (parsed_args["warmup"] = warmup)
	!("preds" in keys(parsed_args)) && (parsed_args["preds"] = nothing; push!(ignores, :preds))
	!("prefix" in keys(parsed_args)) && (parsed_args["prefix"] = ""; push!(ignores, :prefix))
	ps = NamedTuple(Symbol(k)=>v for (k,v) in pairs(parsed_args))

	folders = split(fdir_raw, "/")
	i = findfirst(==("CovidSurvey"), folders)
	fdir = normpath(homedir(),"dev/CovidSurvey", folders[i+1:end-1]...)
	@info "save post-processing in $fdir"
	return (; fdir, ps, ignores)
end

function savechain(pp::PostProcessing3)
	@unpack fdir, ps, ignores, chain = pp
    fname = normpath( fdir, savename(ps.prefix*"CHAIN", ps, "jls"; ignores) )
    safesave( fname, chain )
end

function plot_chains(pp::PostProcessing3; plot_results = false)
	@unpack fdir, ps, ignores, chain = pp
	n = filter( x->!occursin(r"latent_Rt", x), String.(names(chain)))
    p = plot(chain[n]);
    fname = normpath( fdir, savename(ps.prefix*"FIG-CHAINSPLOT", ps, "html"; ignores) )
    savefig(p, fname )
    plot_results && (display(p); run(`firefox $(fname)`, wait=false))
	return nothing
end

function plot_means(pp::PostProcessing3; plot_results=false)
	@unpack fdir, ps, ignores, chain = pp
	n = filter( x->!occursin(r"latent_Rt", x), String.(names(chain)))
	p = meanplot(chain[n]);
	fname = normpath( fdir, savename(ps.prefix*"FIG-MEANPLOT", ps, "html"; ignores) )
	savefig(p, fname )
	plot_results && (display(p); run(`firefox $(fname)`, wait=false))
	return nothing
end

function skip_warmup(pp::PostProcessing3)
	@unpack fdir, ps, ignores, chain, data, model = pp
	chain = chain[ps.warmup+1:end,:,:]
	PostProcessing3(fdir, ps, ignores, chain, data, model)
end

function diagnostics(pp::PostProcessing3)
	@unpack fdir, ps, ignores, chain = pp
	if ps.chains > 1
        @info "gelman diagnostics"
        diagnostics = gelmandiag(chain)
        fname = normpath( fdir, savename(ps.prefix*"GELMANDIAG", ps, "csv"; ignores) )
        safesave( fname, diagnostics)
        pretty_table(diagnostics; crop=:none)
    end
end

function generate_posterior(pp::PostProcessing3)
	@unpack fdir, ps, ignores, chain, model, data = pp

	generated_posterior = let
	    m_pred = model(data.turing_data, true)
	    gq = Turing.generated_quantities(m_pred, chain)
	    vectup2tupvec(reshape(gq, length(gq)))
	end

    fname = normpath( fdir, savename(ps.prefix*"GENERATED-QUANTITIES", ps, "bson"; ignores) )
    dic = Dict( zip( keys(generated_posterior), values(generated_posterior) ) )
    bson( fname ,  dic )
	return generated_posterior
end

function plot_prediction(pp, gp; plot_results = false)
	@unpack fdir, ps, ignores, chain, data = pp
	p = National.plot(data, gp)
	fname = normpath( fdir, savename(ps.prefix*"FIG-PREDICTION", ps, "html"; ignores) )
    savefig( p, fname )
    plot_results && (display(p); run(`firefox $(fname)`, wait=false))
end

function save_rt(p, gp)
	@unpack fdir, ps, ignores, chain, data = p
    df = rtstats(
        gp.Rt,
        data.dates;
        startdate=Date("2020-05-15"),
        stopdate=data.observations_end
    )
    fname = normpath( fdir, savename(ps.prefix*"Rt", ps, "csv"; ignores) )
    save(fname, df)
end

function rtstats(Rt, dates; startdate=Date("2020-05-15"), stopdate=nothing)
    ib = findfirst(==(startdate), dates)
    ie = isnothing(stopdate) ? length(dates) : findfirst(==(Date(stopdate)), dates)
    Rt_array = hcat(Rt...)[ib:ie,:]
    qs = [quantile(v, [0.025, 0.25, 0.5, 0.75, 0.975]) for v in eachrow(Rt_array)]
    llq, lq, mq, uq, uuq = (eachrow(hcat(qs...))..., )

    return DataFrame((; dates = dates[ib:ie], llq, lq, mq, uq, uuq))
end

function save_generated_posterior(p, gp)
	@unpack fdir, ps, ignores, chain, data = p
    fname = normpath( fdir, savename(ps.prefix*"GENERATED-QUANTITIES", ps, "bson"; ignores) )
    dic = Dict( zip( keys(gp), values(gp) ) )
    bson( fname ,  dic )
end

function postprocessing(fname; plot_results = false, exclude = [], warmup = nothing)
	## ==========================================================================
    @info "load data"
    fdir, ps, ignores = parse_fname(fname; warmup)
	p = PostProcessing3(fdir, ps, ignores, exclude, fname)
	savechain(p)

    ## ==========================================================================
    @info "plot chain"
    plot_chains(p; plot_results)

    ## ==========================================================================
    @info "meanplot"
    plot_means(p; plot_results)
    ## ==========================================================================
    # perform diagnostics
	p = skip_warmup(p)
	diagnostics(p)
    ## ==========================================================================
    @info "make predictions"
    gp = generate_posterior(p)
	## ==========================================================================
    @info "save generated quantities"
	save_generated_posterior(p, gp)
    ## ==========================================================================
    @info "plot prediction"
    plot_prediction(p, gp; plot_results)
    # -----------------------------------------------------------------------------
    @info "store reproduction number"
    save_rt(p, gp)
	return nothing
end

name2model = Dict(
    "hospit" => National.model_hospit,
	"deaths" => National.model_deaths,
	"cases"  => National.model_cases,
)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
@recipe function f(bg::BondGraph)
    # Default attributes
    nodecolor --> nodecolours(bg.nodes)
    title --> bg.name
    names --> bg.nodes

    # Forced attributes
    curves := false
    nodeshape := :rect

    GraphRecipes.GraphPlot([bg])
end

function nodecolours(nodes)
    colours = []
    for n in nodes
        type_n = type(n)
        # colour i is the i-th default plotting colour
        c = if type_n in [:Se, :Sf, :SS]
            3
        elseif n isa Component{1}
            1
        elseif n isa Component # more than 1 port
            2
        elseif n isa Junction
            :lightgray
        end
        push!(colours, c)
    end
    colours
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
@test has_cuda(true)
@test has_cuda_gpu(true)

# the API shouldn't have been initialized
@test CuCurrentContext() == nothing

context_cb = Union{Nothing, CuContext}[nothing for tid in 1:Threads.nthreads()]
CUDA.atcontextswitch() do tid, ctx
    context_cb[tid] = ctx
end

task_cb = Union{Nothing, Task}[nothing for tid in 1:Threads.nthreads()]
CUDA.attaskswitch() do tid, task
    task_cb[tid] = task
end

# now cause initialization
ctx = context()
@test CuCurrentContext() == ctx
@test device() == CuDevice(0)
@test context_cb[1] == ctx
@test task_cb[1] == current_task()

fill!(context_cb, nothing)
fill!(task_cb, nothing)

# ... on a different task
task = @async begin
    context()
end
@test ctx == fetch(task)
@test context_cb[1] == nothing
@test task_cb[1] == task

device!(CuDevice(0))
device!(CuDevice(0)) do
    nothing
end

context!(ctx)
context!(ctx) do
    nothing
end

@test_throws AssertionError device!(0, CUDA.CU_CTX_SCHED_YIELD)

device_reset!()

device!(0, CUDA.CU_CTX_SCHED_YIELD)

# test the device selection functionality
if length(devices()) > 1
    device!(0)
    device!(1) do
        @test device() == CuDevice(1)
    end
    @test device() == CuDevice(0)

    device!(1)
    @test device() == CuDevice(1)
end

# test that each task can work with devices independently from other tasks
if length(devices()) > 1
    device!(0)
    @test device() == CuDevice(0)

    task = @async begin
        device!(1)
        @test device() == CuDevice(1)
    end
    fetch(task)

    @test device() == CuDevice(0)

    # reset on a task
    task = @async begin
        device!(1)
        device_reset!()
    end
    fetch(task)
    @test device() == CuDevice(0)

    # tasks on multiple threads
    Threads.@threads for d in 0:1
        for x in 1:100  # give threads a chance to trample over each other
            device!(d)
            yield()
            @test device() == CuDevice(d)
            yield()

            sleep(rand(0.001:0.001:0.01))

            device!(1-d)
            yield()
            @test device() == CuDevice(1-d)
            yield()
        end
    end
    @test device() == CuDevice(0)
end


text:
using Random
using InteractBase, Observables, DataStructures, Colors, WebIO, CSSUtil
using Widgets
using Dates
import InteractBase: widgettype
import Widgets: components

using Test

include("test_observables.jl")


text:
const Transformations{T} = NamedTuple{(:rotation, :translation), Tuple{SMatrix{3, 3, T, 9}, CartesianVector{T}}}
# Transformation rule:
#   1. Rotate
#   2. Translate
# Not sure yet if Rotations will always be of type SMatrix{3, 3, T, 9} -> Any

# Transformations{T}() where {T} = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(CartesianVector{T}))

rotate(p::P, r::AbstractMatrix) where {P <: AbstractPrimitive} = P(p, origin = r * p.origin, rotation = r * p.rotation)
(*)(r::AbstractMatrix, p::AbstractPrimitive) = rotate(p, r)

translate(p::P, v::CartesianVector) where {P <: AbstractPrimitive} = P(p, origin = p.origin + v, rotation = p.rotation)
(+)(p::AbstractPrimitive, v::CartesianVector) = translate(p, v)

transform(g::AbstractPrimitive, t::Transformations) =
    translate(rotate(g, t.rotation), t.translation)

transform(csg::G, t::Transformations) where {G <: AbstractConstructiveGeometry} = G(transform(csg.a, t), transform(csg.b, t))

function combine_transformations(current::Transformations, new::Transformations)
    r = new.rotation * current.rotation
    t = new.rotation * current.translation
    t = t + new.translation
    (rotation = r, translation = t)
end

function Dictionary(m::SMatrix{3,3,T,9}) where {T}
    dict = OrderedDict{String, Any}()
    mat = RotXYZ(m)
    if mat.theta1 == 0 && mat.theta2 == 0
        if mat.theta3 != 0
            dict["Z"] = string(mat.theta3)*"rad"
        end
    else
        dict["M"] = m[:]
    end
    dict
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ f8eb6a38-8b98-11eb-1b3f-73fb0325650f
using Images, ImageView, ImageMorphology, TestImages,MosaicViews

# ╔═╡ 36955aa4-8b99-11eb-073c-1d9d54369888
	img = load("hello-world.png")

# ╔═╡ 954f2008-8b98-11eb-25d0-25bf34030eef
begin

img_erode = Gray.(Gray.(img) .< 0.8); # keeps white objects white
img_erosion1 = erode(img_erode)
img_erosion2 = erode(erode(img_erode))
img_erosion3 = erode(erode(erode(img_erode)))
mosaicview(img_erode, img_erosion1, img_erosion2, img_erosion3; nrow=1)

end

# ╔═╡ abc35624-8b98-11eb-29ce-b33cad8c7338
begin
img_dilate = Gray.(1 * Gray.(img) .> 0.9);
img_dilate1 = dilate(img_dilate)
img_dilate2 = dilate(dilate(img_dilate))
img_dilate3 = dilate(dilate(dilate(img_dilate)))
mosaicview(img_dilate, img_dilate1, img_dilate2, img_dilate3; nrow=1)
end

# ╔═╡ af9519e0-8b98-11eb-1a78-d1757da9252f
begin
img_opening = Gray.(1 * Gray.(img) .> 0.5);
img_opening1 = opening(img_opening)
img_opening2 = opening(opening(img_opening))
img_opening3 = opening(opening(opening(img_opening)))
mosaicview(img_opening, img_opening1, img_opening2, img_opening3; nrow=1)
end

# ╔═╡ b3218a1c-8b98-11eb-0e4c-a58efefb1950
begin
img_closing = Gray.(1 * Gray.(img) .> 0.5);
img_closing1 = closing(img_closing)
img_closing2 = closing(closing(img_closing))
img_closing3 = closing(closing(closing(img_closing)))
mosaicview(img_closing1, img_closing1, img_closing2, img_closing3; nrow=1)
end

# ╔═╡ b69f2e06-8b98-11eb-1776-ebd6d564b297
begin
img_tophat = Gray.(1 * Gray.(img) .> 0.2);
img_tophat1 = tophat(img_tophat)
img_tophat2 = tophat(tophat(img_tophat))
# img_tophat3 = tophat(tophat(tophat(img_tophat)))
mosaicview(img_tophat, img_tophat1, img_tophat2; nrow=1)
end

# ╔═╡ b9c3cf40-8b98-11eb-1896-5bef68754153
begin
img_bothat = Gray.(1 * Gray.(img) .> 0.5);
img_bothat1 = bothat(img_tophat)
img_bothat2 = bothat(bothat(img_tophat))
img_bothat3 = bothat(bothat(bothat(img_tophat)))
mosaicview(img_bothat, img_bothat1, img_bothat2; nrow=1)
end

# ╔═╡ bea4a6ee-8b98-11eb-1a7e-17cb2dd9aa9d
begin
img_gray = Gray.(0.8 * Gray.(img) .> 0.7);
img_morphograd = morphogradient(Gray.(0.8 * Gray.(img) .> 0.4))
img_morpholap = morpholaplace(img_gray)
mosaicview(img_gray, img_morphograd, img_morpholap; nrow=1)
end

# ╔═╡ Cell order:
# ╠═f8eb6a38-8b98-11eb-1b3f-73fb0325650f
# ╠═36955aa4-8b99-11eb-073c-1d9d54369888
# ╠═954f2008-8b98-11eb-25d0-25bf34030eef
# ╠═abc35624-8b98-11eb-29ce-b33cad8c7338
# ╠═af9519e0-8b98-11eb-1a78-d1757da9252f
# ╠═b3218a1c-8b98-11eb-0e4c-a58efefb1950
# ╠═b69f2e06-8b98-11eb-1776-ebd6d564b297
# ╠═b9c3cf40-8b98-11eb-1896-5bef68754153
# ╠═bea4a6ee-8b98-11eb-1a7e-17cb2dd9aa9d

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
using MomentBasedEstimators
using ModelsGenerators

srand(1)
y, x, z = ModelsGenerators.sim_iv_d01(CP = 100);

g(theta) = z.*(y-x*theta);

gmm_fad = GMMEstimator(g, [.0]);
estimate!(gmm_fad);

gmm_ana = GMMEstimator(g, [.0], grad = θ -> -z'x);
estimate!(gmm_ana);

kl_base = MDEstimator(g, [.0]);
estimate!(kl_base);

## Example instrumental Variables

## Dsn = ∂gn/ ∂θ' (m x k)
## Dws = ∂∑wᵢgᵢ/∂θ' (m x k)
## Dsl = ∂∑gᵢ*λ (k x 1)

Dsn(θ)    = -z'x;
Dws(θ, p) = -z'*scale(p,x);
Dsl(θ, λ) = -x.*(z*λ);
Hwsl(θ, p, λ) = zeros(1,1);

kl_ana_grad = MDEstimator(g, [.0], grad = (Dsn, Dws, Dsl));
estimate!(kl_ana_grad);
@time estimate!(kl_ana_grad);
kl_ana_full = MDEstimator(g, [.0], grad = (Dsn, Dws, Dsl, Hwsl));
estimate!(kl_ana_full);
@time estimate!(kl_ana_full);

begin
    srand(1)
    sim = 100;
    n = 100
    k = 4
    coeff = zeros(sim);
    statu = Array(Symbol, sim);
    stder = Array(Float64, sim);
    x_add = randn(n, k-1)
    Dsn(θ)    = -z'x;
    Dws(θ, p) = -z'*scale(p,x);
    Dsl(θ, λ) = -x.*(z*λ);
    Hwsl(θ, p, λ) = zeros(1,1);

    @time for j = 1:sim
        y, x, z = ModelsGenerators.sim_iv_d01(CP = 100, n = n);
        x = [x x_add];
        z = [z x_add];
        kl_base = MDEstimator(g, zeros(k));
        estimate!(kl_base);
        coeff[j] = coef(kl_base)[1];
        statu[j] = status(kl_base);
        stder[j] = stderr(kl_base)[1];
    end
end

begin
    srand(1)
    sim = 1000;
    n = 1000
    k = 4
    coeff = zeros(sim);
    x_add = randn(n, k-1)
    Dsn(θ)    = -z'x;
    Dws(θ, p) = -z'*scale(p,x);
    Dsl(θ, λ) = -x.*(z*λ);
    Hwsl(θ, p, λ) = zeros(k,k);

    @time for j = 1:sim
        y, x, z = ModelsGenerators.sim_iv_d01(CP = 100, n = n);
        x = [x x_add];
        z = [z x_add];
        kl_base = MDEstimator(g, zeros(k), grad = (Dsn, Dws, Dsl, Hwsl));
        estimate!(kl_base);
        coeff[j] = coef(kl_base)[1];
    end
end

# W = inv(z'z);

# W = eye(5);
# out = gmm(g, [.1], W)

# W = pinv(GMM.mfvcov(out, HC0())*100)
# out = GMM.gmm(g, [.0], W)

# coef(out)
# sqrt(vcov(out))

# sqrt(vcov(out, ParzenKernel(3))

# W = eye(5);
# out = GMM.gmm(g, [.1], W, solver = NLoptSolver(algorithm = :LD_LBFGS))

# W = pinv(GMM.mfvcov(out, HC1()))
# out = GMM.gmm(g, [.1], W, solver = NLoptSolver(algorithm = :LD_LBFGS))


text:
of_mean_type(p::Colorant) = ImageSegmentation.meantype(typeof(p))(p)

@testset "Seeded Region Growing" begin
    # 2-D image
    img = zeros(Gray{N0f8}, 10, 10)
    img[6:10,4:8] .= 0.5
    img[3:7,2:6] .= 0.8
    seeds = [(CartesianIndex(3,9),1), (CartesianIndex(5,2),2), (CartesianIndex(9,7),3)]

    expected = ones(Int, 10, 10)
    expected[6:10,4:8] .= 3
    expected[3:7,2:6] .= 2
    expected_labels = [1,2,3]
    expected_means = Dict(1 => of_mean_type(img[3,9]), 2 => of_mean_type(img[5,2]), 3 => of_mean_type(img[9,7]))
    expected_count = Dict(1 => 56, 2 => 25, 3 => 19)

    seg = seeded_region_growing(img, seeds)
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # Custom neighbourhood using a function
    seg = seeded_region_growing(img, seeds, c->[CartesianIndex(c[1]-1,c[2]), CartesianIndex(c[1]+1,c[2]), CartesianIndex(c[1],c[2]-1), CartesianIndex(c[1],c[2]+1)])
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # Offset image
    img = centered(img)
    seeds = [(CartesianIndex(-2,4),1), (CartesianIndex(0,-3),2), (CartesianIndex(4,2),3)]
    expected = centered(expected)
    seg = seeded_region_growing(img, seeds)
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # Custom neighbourhood using a [3,3] vs [5,5] kernel
    img = zeros(Gray{N0f8}, 5, 5)
    img[2:4,2:4] .= 1
    img[3,3] = 0
    seeds = [(CartesianIndex(3,3),1), (CartesianIndex(2,3),2)]

    expected = fill(2,(5,5))
    expected[3,3] = 1
    expected_labels = [1,2]
    expected_means = Dict(1=>Gray{Float64}(0.0), 2=>Gray{Float64}(1/3))
    expected_count = Dict(1=>1, 2=>24)

    seg = seeded_region_growing(img, seeds, (3,3))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
    @test seg.image_indexmap == expected

    expected = ones(Int, 5, 5)
    expected[2:4,2:4] .= 2
    expected[3,3] = 1
    expected_labels = [1,2]
    expected_means = Dict(1=>Gray{N0f8}(0.0), 2=>Gray{N0f8}(1.0))
    expected_count = Dict(1=>17, 2=>8)

    seg = seeded_region_growing(img, seeds, (5,5))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # element-type and InexactError
    img = zeros(Int, 4, 4)
    img[2:end, 2:end] .= 10
    img[end,end] = 11
    seeds = [(CartesianIndex(1,1), 1), (CartesianIndex(2,2), 2)]

    expected = ones(Int, size(img))
    expected[2:end,2:end] .= 2
    expected_labels = [1,2]
    expected_means = Dict(1=>0.0, 2=>91/9)
    expected_count = Dict(1=>7, 2=>9)

    seg = seeded_region_growing(img, seeds, (3,3), (v1,v2)->(δ = abs(v1-v2); return δ > 1 ? δ : zero(δ)))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # 3-d image
    img = zeros(RGB{N0f8},(9,9,9))
    img[3:7,3:7,3:7] .= RGB{N0f8}(0.5,0.5,0.5)
    img[2:5,5:9,4:6] .= RGB{N0f8}(0.8,0.8,0.8)
    seeds = [(CartesianIndex(1,1,1),1), (CartesianIndex(6,4,4),2), (CartesianIndex(3,6,5),3)]

    expected = ones(Int, (9,9,9))
    expected[3:7,3:7,3:7] .= 2
    expected[2:5,5:9,4:6] .= 3
    expected_labels = [1,2,3]
    expected_means = Dict([(i, of_mean_type(img[seeds[i][1]])) for i in 1:3])
    expected_count = Dict(1=>571, 2=>98, 3=>60)

    seg = seeded_region_growing(img, seeds)
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # custom diff_fn
    img = zeros(RGB{N0f8},(3,3))
    img[1:3,1] .= RGB{N0f8}(0.4,1,0)
    img[1:3,2] .= RGB{N0f8}(0.2,1,0)
    seeds = [(CartesianIndex(2,1),1), (CartesianIndex(2,3),2)]

    expected = ones(Int, (3,3))
    expected[1:3,3] .= 2
    expected_labels = [1,2]
    expected_means = Dict(1=>RGB{Float64}(0.3,1.0,0.0), 2=>RGB{Float64}(0.0,0.0,0.0))
    expected_count = Dict(1=>6, 2=>3)

    seg = seeded_region_growing(img, seeds)
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    expected = ones(Int, (3,3))
    expected[1:3,2] .= 0
    expected[1:3,3] .= 2
    expected_labels = [0,1,2]
    expected_means = Dict(1=>RGB{Float64}(0.4,1.0,0.0), 2=>RGB{Float64}(0.0,0.0,0.0))
    expected_count = Dict(0=>3, 1=>3, 2=>3)

    seg = seeded_region_growing(img, seeds, [3,3], (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected
end

@testset "Unseeded Region Growing" begin
    # 2-D image
    img = zeros(Gray{N0f8}, 10, 10)
    img[6:10,4:8] .= 0.5
    img[3:7,2:6] .= 0.8

    expected = ones(Int, 10, 10)
    expected[6:10,4:8] .= 2
    expected[3:7,2:6] .= 3
    expected_labels = [1,2,3]
    expected_means = Dict(1 => of_mean_type(img[3,9]), 3 => of_mean_type(img[5,2]), 2 => of_mean_type(img[9,7]))
    expected_count = Dict(1 => 56, 3 => 25, 2 => 19)

    seg = unseeded_region_growing(img, 0.2)
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # Custom neighbourhood using a function
    seg = unseeded_region_growing(img, 0.2, c->[CartesianIndex(c[1]-1,c[2]), CartesianIndex(c[1]+1,c[2]), CartesianIndex(c[1],c[2]-1), CartesianIndex(c[1],c[2]+1)])
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # Offset image
    img = centered(img)
    expected = centered(expected)
    seg = unseeded_region_growing(img, 0.2)
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # Custom neighbourhood using a [5,5] kernel and varying threshold
    img = zeros(Gray{N0f8}, 5, 5)
    img[2:4,2:4] .= 1
    img[3,3] = 0.8

    expected = fill(1,(5,5))
    expected[2:4,2:4] .= 2
    expected_labels = [1,2]
    expected_means = Dict(1=>Gray{Float64}(0.0), 2=>Gray{Float64}(8.8/9))
    expected_count = Dict(2=>9, 1=>16)

    seg = unseeded_region_growing(img, 0.2, (5,5))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
    @test seg.image_indexmap == expected

    expected = ones(Int, 5, 5)
    expected[2:4,2:4] .= 3
    expected[3,3] = 2
    expected_labels = [1,2,3]
    expected_means = Dict(1=>Gray{N0f8}(0.0), 3=>Gray{N0f8}(1.0), 2=>Gray{N0f8}(0.8))
    expected_count = Dict(1=>16, 2=>1, 3=>8)

    seg = unseeded_region_growing(img, 0.1, (5,5))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    expected = ones(Int, 5, 5)
    expected[2:4,2:4] .= 2
    expected[3,3] = 3
    expected_labels = [1,2,3]
    expected_means = Dict(1=>Gray{N0f8}(0.0), 2=>Gray{N0f8}(1.0), 3=>Gray{N0f8}(0.8))
    expected_count = Dict(1=>16, 3=>1, 2=>8)

    seg = unseeded_region_growing(img, 0.1, (3,3))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # 3-d image
    img = zeros(RGB{N0f8},(9,9,9))
    img[3:7,3:7,3:7] .= RGB{N0f8}(0.5,0.5,0.5)
    img[2:5,5:9,4:6] .= RGB{N0f8}(0.8,0.8,0.8)

    expected = ones(Int, (9,9,9))
    expected[3:7,3:7,3:7] .= 2
    expected[2:5,5:9,4:6] .= 3
    expected_labels = [1,2,3]
    expected_means = Dict(1=>of_mean_type(img[1,1,1]), 2=>of_mean_type(img[3,3,3]), 3=>of_mean_type(img[2,5,4]))
    expected_count = Dict(1=>571, 2=>98, 3=>60)

    seg = unseeded_region_growing(img, 0.2)
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

    # custom diff_fn
    img = zeros(RGB{N0f8},(3,3))
    img[1:3,1] .= RGB{N0f8}(0.4,1,0)
    img[1:3,2] .= RGB{N0f8}(0.2,1,0)

    expected = ones(Int, (3,3))
    expected[1:3,2] .= 2
    expected[1:3,3] .= 3
    expected_labels = [1,2,3]
    expected_means = Dict(1=>of_mean_type(img[1,1]), 3=>of_mean_type(img[1,3]), 2=>of_mean_type(img[1,2]))
    expected_count = Dict(1=>3, 2=>3, 3=>3)

    seg = unseeded_region_growing(img, 0.2, [3,3], (c1,c2)->abs(of_mean_type(c1).r - of_mean_type(c2).r))
    @test all(label->(label in expected_labels), seg.segment_labels)
    @test all(label->(label in seg.segment_labels), expected_labels)
    @test expected_count == seg.segment_pixel_count
    @test expected_means == seg.segment_means
    @test seg.image_indexmap == expected

end


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

JLLWrappers.@generate_main_file_header("MSTM")
JLLWrappers.@generate_main_file("MSTM", UUID("2c0bff41-3d2c-5a3a-b138-70c96ef824a5"))
end  # module MSTM_jll


text:
@testset "adaptively sampled distance fields" begin
    include("../examples/adaptive_distance_fields.jl")
    import .AdaptivelySampledDistanceFields: ASDF, evaluate
    import StaticArrays: SVector

    signed_distance = x -> norm(x - SVector(0, 0))
    origin = SVector(-1, -1)
    widths = SVector(2, 2)
    root = ASDF(signed_distance, origin, widths, 1e-2, 1e-2)
    @test evaluate(root, [0, 0]) == 0
    @test evaluate(root, SVector(0, 0)) == 0
    @test findleaf(root, [0.01, -0.01]).boundary.origin == [0.0, -0.0625]
    @test findleaf(root, [0.01, -0.01]).boundary.widths == [0.0625, 0.0625]

    for leaf in allleaves(root)
        @test isleaf(leaf)
    end
    @test_broken length(collect(allleaves(root))) == 112
    @test_broken length(collect(allcells(root))) == 149
end


text:
"""
    boolean_indicator(f)

If `f` is a boolean indicator where 0's correspond to
background and 1s correspond to foreground then mark
background pixels with large number `1e10`
"""
boolean_indicator(f) = @. ifelse(f == 0, 1.0f10, 0.0f0)


text:
using SummationByParts
include("nodecalc3.jl")

function nodecalc(sbp::TriSBP, isDG::Bool)
  vtx = [0.0 0; 1 0; 0 1]  # why is this 0 to 1 but the other method is -1 to 1
  r1 = vtx[1, :]
  r2 = vtx[2, :]
  r3 = vtx[3, :]

  T = zeros(2,2)
  T[:, 1] = r2 - r1
  T[:, 2] = r3 - r1

  if isDG
    coords = SummationByParts.SymCubatures.calcnodes(sbp.cub, vtx)
  else
    coords = calcnodes(sbp, vtx)
  end

  xi = zeros(coords)

  for i=1:size(coords,2)
    xi[:, i] = T\(coords[:, i] - r1)
  end

  return xi, coords
end

function nodecalc(sbp::TetSBP, isDG::Bool)

  vtx = [-1.0 -1.0 -1.0
        1.0 -1.0 -1.0
       -1.0 1.0 -1.0
       -1.0 -1.0 1.0]

  r1 = vtx[1, :]
  r2 = vtx[2, :]
  r3 = vtx[3, :]
  r4 = vtx[4, :]
  T = zeros(3,3)
  T[:, 1] = r2 - r1
  T[:, 2] = r3 - r1
  T[:, 3] = r4 - r1

  if isDG
    coords = SummationByParts.SymCubatures.calcnodes(sbp.cub, vtx)
  else
    coords = calcnodes(sbp, vtx)
  end

  xi = zeros(coords)

  for i=1:size(coords,2)
    xi[:, i] = T\(coords[:, i] - r1)
  end

  return xi, coords
end

function minNodeDist(sbp, isDG::Bool)
# get the minimum distance between nodes on a reference element of degree p

    xi, coords = nodecalc(sbp, isDG)
    min_dist = typemax(Float64)
#    println("coords = ", coords)
    for i=1:size(coords, 2)
      for j=(i+1):size(coords, 2)
	# calculate distance between node i and node j
	dist_j = norm(coords[:, i] - coords[:, j])

	if dist_j < min_dist
	  min_dist = dist_j
	end

      end  # end loop j
    end  # end loop i

#    println("for p=$p elements, min node distance = ", min_dist)

  return min_dist
end

# sbp = TriSBP{Float64}(degree=1, reorder=false, internal=false)
# sbp = getTriSBPWithDiagE(degree=4)
# xi, coords = nodecalc(sbp, true)
# printCaseStatement(xi)
# printBaryCoords(xi, coords)
# writedlm("coordsout.dat", coords, ' ')
#minNodeDist(2)


text:
@testset "argmin_x" begin

#test argmin_x.jl for CSC sparse matrix format
A=sprandn(100,100,0.01)+sparse(I, 100, 100);
A=A'*A;
while rank(Array(A))<100
  A=A+sparse(I, 100, 100)
end
xt=randn(100)
b=A*xt;

maxit=100;
pp=2
p=2
TF=Float64
log_PARSDMM = log_type_PARSDMM(zeros(maxit,pp),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),
zeros(maxit),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),TF(0),TF(0),TF(0),TF(0),TF(0),TF(0),TF(0));

#test for zero initial guess
x_min_solver="CG_normal"
(x,iter,relres)=argmin_x(A,b,zeros(100),x_min_solver,1e-5,1,log_PARSDMM)
@test norm(A*x-b)/norm(b) <= 1e-5
@test isapprox(norm(A*x-b)/norm(b),relres,rtol=1e-5)

#test for good initial guess
x_min_solver="CG_normal"
(x2,iter2,relres2)=argmin_x(A,b,xt+randn(100).*1e-4,x_min_solver,1e-5,1,log_PARSDMM)
@test iter2<iter
@test norm(A*x2-b)/norm(b) <= 1e-5
@test isapprox(norm(A*x2-b)/norm(b),relres2,rtol=1e-5)

#test if 10*eps(Float32/Float64) is achieved
x_min_solver="CG_normal"
(x,iter,relres)=argmin_x(A,b,zeros(100),x_min_solver,10*eps(),1,log_PARSDMM)
@test norm(A*x-b)/norm(b) <= 20*eps()

#test argmin_x.jl for CDS sparse matrix format
A=sprandn(100,100,0.01)+sparse(I, 100, 100);
A=A'*A;
while rank(Array(A))<100
  A=A+sparse(I, 100, 100)
end
xt=randn(100)
b=A*xt;

(R,offset) = mat2CDS(A)

x_min_solver="CG_normal"
y=zeros(Float64,length(b))
(x,iter,relres)=argmin_x(A,b,zeros(100),x_min_solver,1e-5,1,log_PARSDMM,offset,y)

@test norm(A*x-b)/norm(b) <= 1e-5
@test norm(A*x-b)/norm(b) <= 2.0*relres

#test for more accurate solutions with CDS
y=zeros(Float64,length(b))
(x,iter,relres)=argmin_x(A,b,zeros(100),x_min_solver,1e-10,1,log_PARSDMM,offset,y)

@test norm(A*x-b)/norm(b) <= 1e-10
@test norm(A*x-b)/norm(b) <= 2.0*relres

end


text:
q = quote
    x = 1
    y = 2
    c = x +2
    print(c)
end
typeof(q)

max_stars_repo_name: zqsheng/snippet

max_issues_repo_name: zqsheng/snippet

max_forks_repo_name: zqsheng/snippet


text:
"""
	ImplicitTraj{T}
This structure holds the trajectory of evaluations and Jacobians of the implicit dynamics.
These evaluations and Jacobians are computed using a linearizedimation computed around `lin`.
"""
mutable struct ImplicitTraj{T}
	H::Int
	lin::Vector{LinearizedStep{T}}
	d::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}} # dynamics violation
	dq2::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}}
	dγ1::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}}
	db1::Vector{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int}},true}}
	# δq0::Vector{SubArray{T,2,SparseArrays.SparseMatrixCSC{T,Int},Tuple{UnitRange{Int},UnitRange{Int}},false}}  # q0 solution gradient length=H
	# δq1::Vector{SubArray{T,2,SparseArrays.SparseMatrixCSC{T,Int},Tuple{UnitRange{Int},UnitRange{Int}},false}}  # q1 solution gradient length=H
	# δu1::Vector{SubArray{T,2,SparseArrays.SparseMatrixCSC{T,Int},Tuple{UnitRange{Int},UnitRange{Int}},false}}  # u1 solution gradient length=H
	δq0::Vector{SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}} # q0 solution gradient length=H
	δq1::Vector{SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}}  # q1 solution gradient length=H
	δu1::Vector{SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}}  # u1 solution gradient length=H
	ip::Vector{<:AbstractIPSolver}
	mode::Symbol
end

function ImplicitTraj(ref_traj::ContactTraj, s::Simulation;
	κ = ref_traj.κ[1],
	max_time = 1e5,
	mode = :configurationforce,
	ip_type::Symbol = :interior_point,
	opts = eval(interior_point_options(ip_type))(
			undercut = 5.0,
			γ_reg = 0.1,
			κ_tol = κ[1],
			r_tol = 1.0e-8,
			diff_sol = true,
			solver = :empty_solver,
			max_time = max_time))

	model = s.model
	env = s.env

	H = ref_traj.H

	nq = model.dim.q
	nu = model.dim.u
	nc = model.dim.c
	nb = nc * friction_dim(env)
	if mode == :configurationforce
		nd = nq + nc + nb
	elseif mode == :configuration
		nd = nq
	else
		@error "invalid mode"
	end
	nz = num_var(model, env)
	nθ = num_data(model)

	lin = [LinearizedStep(s, ref_traj.z[t], ref_traj.θ[t], κ) for t = 1:H]

	ip =  [eval(ip_type)(
			 deepcopy(ref_traj.z[t]),
			 deepcopy(ref_traj.θ[t]),
			 oss = OptimizationSpace13(model, env),
			 r! = r!,
			 rz! = rz!,
			 rθ! = rθ!,
			 r  = RLin(s, lin[t].z, lin[t].θ, lin[t].r, lin[t].rz, lin[t].rθ),
			 rz = RZLin(s, lin[t].rz),
			 rθ = RθLin(s, lin[t].rθ),
			 opts = opts) for t = 1:H]

	# views
	d = [view(ip[t].z, 1:nd) for t = 1:H]
	dq2 = [view(ip[t].z, 1:nq) for t = 1:H]
	if mode == :configurationforce
		dγ1 = [view(ip[t].z, nq .+ (1:nc)) for t = 1:H]
		db1 = [view(ip[t].z, nq + nc .+ (1:nb)) for t = 1:H]
	elseif mode == :configuration
		dγ1 = [view(ip[t].z, 1:0) for t = 1:H]
		db1 = [view(ip[t].z, 1:0) for t = 1:H]
	else
		@error "invalid mode"
	end

	off = 0
	δq0 = [view(ip[t].δz, 1:nd, off .+ (1:nq)) for t = 1:H]; off += nq
	δq1 = [view(ip[t].δz, 1:nd, off .+ (1:nq)) for t = 1:H]; off += nq
	δu1 = [view(ip[t].δz, 1:nd, off .+ (1:nu)) for t = 1:H]; off += nu

	return ImplicitTraj(H, lin, d, dq2, dγ1, db1, δq0, δq1, δu1, ip, mode)
end

function update!(im_traj::ImplicitTraj, ref_traj::ContactTraj,
	s::Simulation,
	 alt::Vector
	; κ = ref_traj.κ[1])

	H = ref_traj.H
	for t = 1:H
		fill!(im_traj.ip[t].κ, κ)
	end
	for t = 1:H-1
		im_traj.lin[t]   = im_traj.lin[t+1]
		im_traj.ip[t].r  = im_traj.ip[t+1].r
		im_traj.ip[t].rz = im_traj.ip[t+1].rz
		im_traj.ip[t].rθ = im_traj.ip[t+1].rθ

		# altitude
		im_traj.ip[t].r.alt = alt
	end

	update!(im_traj.lin[H], s, ref_traj.z[H], ref_traj.θ[H])

	z0  = im_traj.lin[H].z
	θ0  = im_traj.lin[H].θ
	r0  = im_traj.lin[H].r
	rz0 = im_traj.lin[H].rz
	rθ0 = im_traj.lin[H].rθ

	update!(im_traj.ip[H].r,  z0, θ0, r0, rz0, rθ0)
	update!(im_traj.ip[H].rz, rz0)
	update!(im_traj.ip[H].rθ, rθ0)

	# altitude
	im_traj.ip[H].r.alt = alt
	return nothing
end

function set_altitude!(im_traj::ImplicitTraj, alt::Vector)
	H = im_traj.H
	for t = 1:H
		# altitude
		set_altitude!(im_traj.ip[t], alt)
	end
	return nothing
end

function set_altitude!(ip::InteriorPoint, alt::Vector)
	# altitude
	ip.r.alt = alt
	return nothing
end

"""
	implicit_dynamics!(im_traj::ImplicitTraj, model::ContactModel,
		traj::ContactTraj; κ = traj.κ)
Compute the evaluations and Jacobians of the implicit dynamics on the trajectory 'traj'. The computation is
linearized since it relies on a linearization about a reference trajectory.
"""
function implicit_dynamics!(im_traj::ImplicitTraj, s::Simulation,
	traj; κ = traj.κ) where {T, D}

	model = s.model
	env = s.env

	# Threads.@threads for t = 1:traj.H
	for t = 1:traj.H
		# initialized solver
		z_initialize!(im_traj.ip[t].z, model, env, copy(traj.q[t+2])) #TODO: try alt. schemes
		im_traj.ip[t].θ .= traj.θ[t]

		# solve
		status = interior_point_solve!(im_traj.ip[t])

		# !status && error("implicit dynamics failure (t = $t)")
		!status && (@warn "implicit dynamics failure (t = $t)")

		# compute dynamics violation
		# println("before!: im_traj.dq2[t] = ", scn.(im_traj.dq2[t]))
		im_traj.dq2[t] .-= traj.q[t+2]
		# println("implicit_dynamics!:    traj.q[t+2] = ", scn.(traj.q[t+2]))
		# println(" after!: im_traj.dq2[t] = ", scn.(im_traj.dq2[t]))

		if im_traj.mode == :configurationforce
			im_traj.dγ1[t] .-= traj.γ[t]
			im_traj.db1[t] .-= traj.b[t]
		elseif im_traj.mode == :configuration
			nothing
		end
	end
	return nothing
end


text:
using StatisticalRethinking

ways  = [0, 3, 8, 9, 0];

ways/sum(ways)

d = Binomial(9, 0.5)

pdf(d, 6)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
struct FluidParams
    eta2d :: Float64 # length scale introduced by difference in viscosity of of surface fluid and external fluids
end

abstract type Object end


text:
module RFFT

using FFTW, LinearAlgebra

export RCpair, plan_rfft!, plan_irfft!, rfft!, irfft!, normalization

import Base: real, complex, copy, copy!

mutable struct RCpair{T<:AbstractFloat,N,RType<:AbstractArray{T,N},CType<:AbstractArray{Complex{T},N}}
    R::RType
    C::CType
    region::Vector{Int}
end

function RCpair{T}(::UndefInitializer, realsize::Dims{N}, region=1:length(realsize)) where {T<:AbstractFloat,N}
    sz = [realsize...]
    firstdim = region[1]
    sz[firstdim] = realsize[firstdim]>>1 + 1
    sz2 = copy(sz)
    sz2[firstdim] *= 2
    R = Array{T,N}(undef, (sz2...,)::Dims{N})
    C = unsafe_wrap(Array, convert(Ptr{Complex{T}}, pointer(R)), (sz...,)::Dims{N}) # work around performance problems of reinterpretarray
    RCpair(view(R, map(n->1:n, realsize)...), C, [region...])
end

RCpair(A::Array{T}, region=1:ndims(A)) where {T<:AbstractFloat} = copy!(RCpair{T}(undef, size(A), region), A)

real(RC::RCpair)    = RC.R
complex(RC::RCpair) = RC.C

copy!(RC::RCpair, A::AbstractArray{T}) where {T<:Real} = (copy!(RC.R, A); RC)
function copy(RC::RCpair{T,N}) where {T,N}
    C = copy(RC.C)
    R = reshape(reinterpret(T, C), size(parent(RC.R)))
    RCpair(view(R, RC.R.indices...), C, copy(RC.region))
end

# New API
rplan_fwd(R, C, region, flags, tlim) =
    FFTW.rFFTWPlan{eltype(R),FFTW.FORWARD,true,ndims(R)}(R, C, region, flags, tlim)
rplan_inv(R, C, region, flags, tlim) =
   FFTW.rFFTWPlan{eltype(R),FFTW.BACKWARD,true,ndims(R)}(R, C, region, flags, tlim)
function plan_rfft!(RC::RCpair{T}; flags::Integer = FFTW.ESTIMATE, timelimit::Real = FFTW.NO_TIMELIMIT) where T
    p = rplan_fwd(RC.R, RC.C, RC.region, flags, timelimit)
    return Z::RCpair -> begin
        FFTW.assert_applicable(p, Z.R, Z.C)
        FFTW.unsafe_execute!(p, Z.R, Z.C)
        return Z
    end
end
function plan_irfft!(RC::RCpair{T}; flags::Integer = FFTW.ESTIMATE, timelimit::Real = FFTW.NO_TIMELIMIT) where T
    p = rplan_inv(RC.C, RC.R, RC.region, flags, timelimit)
    return Z::RCpair -> begin
        FFTW.assert_applicable(p, Z.C, Z.R)
        FFTW.unsafe_execute!(p, Z.C, Z.R)
        rmul!(Z.R, 1 / prod(size(Z.R)[Z.region]))
        return Z
    end
end
function rfft!(RC::RCpair{T}) where T
    p = rplan_fwd(RC.R, RC.C, RC.region, FFTW.ESTIMATE, FFTW.NO_TIMELIMIT)
    FFTW.unsafe_execute!(p, RC.R, RC.C)
    return RC
end
function irfft!(RC::RCpair{T}) where T
    p = rplan_inv(RC.C, RC.R, RC.region, FFTW.ESTIMATE, FFTW.NO_TIMELIMIT)
    FFTW.unsafe_execute!(p, RC.C, RC.R)
    rmul!(RC.R, 1 / prod(size(RC.R)[RC.region]))
    return RC
end

@deprecate RCpair(realtype::Type{T}, realsize, region=1:length(realsize)) where T<:AbstractFloat RCpair{T}(undef, realsize, region)

end


text:
using Gadfly

# Issue #678
x = repeat(collect(1:10), inner=[10])
y = repeat(collect(1:10), outer=[10])
plot(x=x,y=y,color=x+y, Geom.rectbin,
            Scale.ContinuousColorScale(Scale.lab_gradient(colorant"green",
                                                          colorant"white",
                                                          colorant"red")))


text:
using MultivariateMedians
using Test

import Documenter

@testset "MultivariateMedians.jl" begin
    @testset "unit tests" begin
        @testset "geometric.jl" begin
        end

        @testset "marginal.jl" begin
        end

        @testset "norm.jl" begin
            @testset "L2 norm and L2 distance" begin
                @test MultivariateMedians.l2_norm([3.0, 4.0]) ≈ 5.0
                @test MultivariateMedians.l2_distance([3.0, 4.0], [0.0, 0.0]) ≈ 5.0
            end
        end
    end

    @testset "doctests" begin
        Documenter.doctest(MultivariateMedians)
    end
end


text:
mutable struct Snake
    body::Vector
    dir::Vector
    @init function Snake()
        this.body = [[w/2, h/2]]
        this.dir = [1, 0]
    end
end

setdir!(this::Snake, dir::Vector) = (this.dir = dir;)

function update!(this::Snake)
    head = [0., 0.]
    copy!(head, last(this.body))
    popfirst!(this.body)
    head += this.dir
    push!(this.body, head)
end

function draw(this::Snake)
    sethue(colors.snake)
    foreach(this.body) do pos
        rect(pos..., 1, 1, :fill)
    end
end

function eat(this::Snake, food::Food)
    x, y = last(this.body)
    return x == food.x && y == food.y
end

function grow(this::Snake)
    head = [0., 0.]
    copy!(head, last(this.body))
    push!(this.body, head)
end

function endgame(this::Snake)
    x, y = last(this.body)

    if x < 0 || x > w - 1 || y < 0 || y > h - 1
        return true
    end

    for pos in this.body[1:end-1]
        if pos == [x, y]
            return true
        end
    end

    return false
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
using RoadGraphs

let
    R0 = ((-36.9323, 174.8616), (-36.8878, 174.9565))
    G = read_graph("graph.json", R0...)

    # 13 configs x 5 minutes = 1 hour
    for λ in [0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1, 2.3, 2.5]
        MP = ModelParams(
            lambda = λ,
            bidir = true,
            junctions = true,
            )
        SP = SolveParams(
            time_limit = 60 * 5,
            )

        run_model(G, MP, SP, "lambda-$(λ)", draw = true)
    end
end


text:
using ControlSystems
using FFTW
using LaTeXStrings
using PyPlot
using Printf

function dsdemo5()
    show_usage = true
    osr = 25
    ntf = synthesizeNTF(6, osr, 1, 4)
    M = 16
    N = 2^14
    sigma_d = 0.01  # 1% mismatch
    window = ds_hann(N) / (M*N/8)
    windowM = repeat(window', outer=(M,1))

    mtf1 = zpk([1], [0], 1, 1)                  # First-order shaping
    mtf2 = zpk([1, 1], [0.3, 0.3], 1, 1)        # Second-order shaping
    mtf2b = synthesizeNTF(2, osr, 1, 2)         # Second-order shaping
    mtf4 = synthesizeNTF(4, osr*0.9, 1, 1.3)    # Fourth-order shaping

    f  = 0.01
    dw = []
    cases = [
        ( undbv(-3), nothing, 0.0, "Thermometer"        , "Thermometer"),
        ( undbv(-3),    mtf1, 0.0, "Rotation"           , "Rotation"),
        ( undbv(-3),    mtf2, 0.0, "2nd-order"          , string(L"2^{nd}", "-order")),
        (undbv(-30),    mtf1, 0.0, "Rotation"           , "Rotation"),
        (undbv(-30),    mtf1, 0.5, "Rot + dither"       , "Rot + dither"),
        ( undbv(-3),   mtf2b, 0.0, "2nd-order with zero", string(L"2^{nd}", "-order with zero")),
        ( undbv(-3),    mtf4, 0.0, "4th-order"          , string(L"4^{th}", "-order")),
    ]
    comparisons = [[1,2], [2,3], [4,5], [3,6], [6,7]]

    sv = Array{Matrix{Float64}}(undef, length(cases))
    leg = Array{String}(undef, length(cases))
    ltx_leg = Array{String}(undef, length(cases))
    Svv = Array{Vector{Float64}}(undef, length(cases))
    Sdd = Array{Vector{Float64}}(undef, length(cases))
    fin = zeros(Int, length(cases))
    for i = 1:length(cases)
        A = cases[i][1]
        mtf = cases[i][2]
        dither = cases[i][3]
        leg_i = cases[i][4]
        ltx_leg_i = cases[i][5]

        fin[i] = round(Int, f*N)
        inband = setdiff(2:(1+ceil(Int, 0.5*N/osr)), 1 .+ [0; 1; fin[i] .+ [-1, 0, 1]])
        w = (2*pi/N) * fin[i]
        u = M * A * sin.(w * (0:N-1))
        v, = simulateDSM(u, ntf, M+1)    # M unit elements requires an (M+1)-level quant.
        Svv[i] = abs.(fft(v .* window)).^2
        if isnothing(mtf)
            sv[i] = ds_therm(v, M)
        else
            sv[i], = simulateMS(v, M, mtf, dither, dw)
        end
        Sdd[i] = sigma_d^2 * sum(abs.(fft(sv[i] .* windowM, 2)).^2, dims=1)[:]
        mnp = sum(Sdd[i][inband]) / 1.5
        leg[i] = @sprintf("%s (MNP = %.0f dBFS)", leg_i, dbp(mnp))
        ltx_leg[i] = @sprintf("%s (MNP = %.0f dBFS)", ltx_leg_i, dbp(mnp))
    end

    # Plot results
    for comp_i = 1:length(comparisons)
        case_nums = comparisons[comp_i]
        println("Mismatch-Shaping Unit-Element DAC")
        @printf("Comparing %s vs. %s.\n", leg[case_nums[1]], leg[case_nums[2]])
        nc = length(case_nums)

        if show_usage
            fig = figure(1)
            fig.clf()
            fig.canvas.set_window_title("Element Usage")
            T = 25
            for i = 1:nc
                subplot(nc, 1, i)
                ci = case_nums[i]
                plotUsage(sv[ci][:,1:T])
            end
            fig.canvas.manager.window.raise_()
        end

        fig = figure(2)
        fig.clf()
        fig.canvas.set_window_title("Error Spectra")
        cols = ["b", "m", "r"]
        cleg = Array{String}(undef, nc)
        for i = 1:nc
            ci = case_nums[i]
            plotSpectrum(sqrt.(Sdd[ci]), fin[ci], cols[i], 8, 4)
            cleg[i] = ltx_leg[ci]
        end
        ci = case_nums[nc]
        plotSpectrum(sqrt.(Svv[ci]), fin[ci], "g", 8, 5)
        axis([1e-3, 0.5, -140, -50])
        plot([1e-3, 0.5/osr], -140*[1, 1], color="k", linewidth=4)
        text(0.5, -140, @sprintf("NBW=%.1e", 1.5/N), horizontalalignment="right", verticalalignment="bottom")
        grid()
        grid(which="minor", linestyle="dotted")
        ylabel("Error PSD")
        xlabel("Normalized Frequency")
        title(@sprintf("A = %.0fdBFS", dbp(Svv[ci][fin[ci]])))
        legend(cleg, loc="upper right")

        fig.canvas.manager.window.raise_()

        println("paused")
        readline()
    end

end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
import Common.test_clone

function run_test_clone()
    test_clone(Atom(index=1), Dict(:residue_id => 2), Atom(index=1, residue_id=2), "Atom test_clone 1")
end


text:
print("Importing packages")
include("CartPendulum.jl")
import .CartPendulum
include("../src/MinimaxAdaptiveControl.jl")
using .MinimaxAdaptiveControl
using Test
using ControlSystems
using LinearAlgebra
using DifferentialEquations
using Plots

function f!(dx::AbstractVector{T},
    x::AbstractVector{T},
    p::Tuple{CartPendulum.Pendulum{T}, MAController}, t::T) where T<:Real
    return CartPendulum.f!(dx, x, p[1], t)
end

function g!(dx::AbstractVector{T},
    x::AbstractVector{T},
    p::Tuple{CartPendulum.Pendulum{T}, MAController}, t::T) where T<:Real
    return CartPendulum.g!(dx, x, p[1], t)
end

function controller(integrator::DiffEqBase.DEIntegrator)
    x=copy(integrator.u[1:4])
    x[3] = mod(x[3],2*pi) - pi
    u=integrator.u[5]
    (_,mac) = integrator.p
    update!(mac,x,[u])
    integrator.u[5] = (-K(mac)*x)[1]
    integrator.u[6] = mac.candidates[1].hist[]
    integrator.u[7] = mac.candidates[2].hist[]
end

function generategif(sol, sys,tspan)
    dt = 0.05
    tu = tspan[1]:dt:tspan[2]
    (M, m, l, b, ulims) = (sys.M, sys.m, sys.l, sys.b, sys.ulims)
    θ = sol(tu)[3,:]

    x1 = sol(tu)[1, :]
    y1 = zeros(length(tu))
    x2 = x1 .+ l*sin.(θ)
    y2 = -l*cos.(θ)

    gr()
    previous_GKSwstype = get(ENV, "GKSwstype", "")
    ENV["GKSwstype"] = "100"
    anim = Animation()
    for i = 1:length(tu)
        str = string("Time = ", round(tu[i]), " sec")
        plot(
            [x1[i], x2[i]],
            [y1[i], y2[i]],
            markersize = 10,
            markershape = :circle,
            label = "",
            title = str,
            title_location = :left,
            aspect_ratio = :equal,
            xlim = (x1[i] - 1, x1[i] + 1),
            ylim = (-1, 1),
            size = (400, 300),
        )
        if i > 9
            plot!(
                [x2[i-3:i]],
                [y2[i-3:i]],
                alpha = 0.15,
                linewidth = 2,
                color = :red,
                label = "",
            )
            plot!(
                [x2[i-5:i-3]],
                [y2[i-5:i-3]],
                alpha = 0.08,
                linewidth = 2,
                color = :red,
                label = "",
            )
            plot!(
                [x2[i-7:i-5]],
                [y2[i-7:i-5]],
                alpha = 0.04,
                linewidth = 2,
                color = :red,
                label = "",
            )
            plot!(
                [x2[i-9:i-7]],
                [y2[i-9:i-7]],
                alpha = 0.01,
                linewidth = 2,
                color = :red,
                label = "",
            )
        end
        frame(anim)
    end
    ENV["GKSwstype"] = previous_GKSwstype
    return anim

end
print("Setting up system dynamics")
sys = CartPendulum.Pendulum()
Δt = 0.01;
x0 = [0; 0; π - 0.05; 0; 0]
tfinal = 20.0

linsys1 = CartPendulum.linearize(sys)
linsys1_d = c2d(linsys1, Δt, :zoh)[1]
(A1,B1,_, _) = ControlSystems.ssdata(linsys1_d)

As = [A1,A1]
Bs = [B1, -B1]
gamma = 1e4
q = [1 1 1e3 1e3];
Q = diagm(q[:])
R = 1e-2*I(1)
mac = MAController(As, Bs, Q, R, gamma, x0[1:4],.9)

tspan = (0,tfinal)
sys.nx = 0.1
sys.nθ = 0.02
condition(u,t,integrator) = t==10
affect!(integrator) = integrator.p[1].sgn = -1.

cb = PeriodicCallback(controller, Δt,
initial_affect=true, save_positions=(false,false))
cb2 = DiscreteCallback(condition,affect!)

prob = SDEProblem(f!, g!, [x0;0;0;0], tspan, (sys, mac), callback = CallbackSet(cb,cb2))

@time sol = solve(prob)

gifname = "gain_switch_fgfac.gif"
@time anim = generategif(sol, sys, tspan)
gif(anim, gifname, fps=20)
plot(1:0.1:tfinal,sol(1:0.1:tfinal)[6:7,:]',
labels=["sum(λ^k||*Ax_k + Bu_k - x_{k+1}||)" "sum(λ^k*||Ax_k + -Bu_k - x_{k+1}||)"],
title="Inverted Pendulum - Gain switch @ 10 s",
legend=:right)
savefig("gain_switch_fgfac.png")


text:
using Printf

include("constants.jl")
include("alias.jl")
include("utils.jl")
include("PGBF.jl")
include("CGBF.jl")
include("overlap.jl")
include("Atoms.jl")
include("sto3g.jl")
include("BasisSet.jl")

function test_old(filename)
    atoms = init_atoms_xyz(filename)
    println(atoms)

    basis = build_basis(atoms)
    Nbasis = length(basis)

    println("\nNbasis = ", Nbasis)
    for i = 1:Nbasis
        @printf("------------\n")
        @printf("basis #%d\n", i)
        @printf("------------\n")
        println(basis[i])
    end
end

function test_new(filename)
    atoms = init_atoms_xyz(filename)
    println(atoms)

    basis = build_basis_new(atoms)
    Nbasis = length(basis)

    println("\nNbasis = ", Nbasis)
    for i = 1:Nbasis
        @printf("------------\n")
        @printf("basis #%d\n", i)
        @printf("------------\n")
        println(basis[i])
    end
end

#test_old("H2O.xyz")
test_new("H2O.xyz")


text:
# Script to automatically insert Markdown footnotes for all [#xxxx] issue
# cross-references in the CHANGELOG file.

NEWS = get(ARGS, 1, "CHANGELOG.md")

s = read(NEWS, String)

m = match(r"\[#[0-9]+\]:", s)
if m !== nothing
    s = s[1:m.offset-1]
end

footnote(n) = "[#$n]: https://github.com/JuliaLang/Pkg.jl/issues/$n"
N = map(m -> parse(Int,m.captures[1]), eachmatch(r"\[#([0-9]+)\]", s))
foots = join(map(footnote, sort!(unique(N))), "\n")

open(NEWS, "w") do f
    println(f, s, "<!--- Generated by NEWS-update.jl --->\n", foots)
end


text:
function to_ast(::Type{T}) where {T}
    io = IOBuffer()
    print(io, T)
    str = String(take!(io))
    ast = Meta.parse(str)
    return ast
end

function to_pascal_case(s::Symbol)
    str = String(s)
    new_str = ""
    next_upper = true
    for letter in str
        if next_upper
            new_str *= uppercase(letter)
            next_upper = false
        elseif letter == '_'
            next_upper = true
        else
            new_str *= letter
        end
    end

    if new_str[end] == 's'
        new_str = new_str[1:end-1]
    end

    return Symbol(new_str)
end

function write_exprs(exprs::Vector, fname::AbstractString)
    open(fname, "w") do io
        for expr in exprs
            str = repr(expr)[3:end-1] # removes :( and )
            str = replace(str, "\n  " => "\n")
            write(io, str)
            write(io, "\n\n")
        end
    end
end

function to_exprs(t, n)
    exprs = []
    to_expr(t, n, exprs)
    return exprs
end

function to_expr(::Type{NamedTuple{N,T}}, root_name::Symbol, exprs) where {N,T<:Tuple}
    sub_exprs = []
    for (n, t) in zip(N, T.types)
        push!(sub_exprs, to_field_expr(t, n, exprs))
    end
    struct_name = to_pascal_case(root_name)
    push!(exprs, Expr(:struct, false, struct_name, Expr(:block, sub_exprs...)))
    return struct_name
end

function to_expr(::Type{Array{T,N}}, root_name::Symbol, exprs) where {T<:NamedTuple,N}
    return to_expr(T, root_name, exprs)
end

function to_expr(t::Type{T}, root_name::Symbol, exprs) where {T}
    if T isa Union
        return Expr(
            :curly,
            :Union,
            to_expr(t.a, root_name, exprs),
            to_expr(t.b, root_name, exprs),
        )
    else
        return to_ast(T)
    end
end

# given the type of a field of a struct, return a node for that field's name/type
function to_field_expr(t::Type{NamedTuple{N,T}}, root_name::Symbol, exprs) where {N,T}
    to_expr(t, root_name, exprs)
    return Expr(:(::), root_name, to_pascal_case(root_name))
end

function to_field_expr(::Type{Array{T,N}}, root_name::Symbol, exprs) where {T,N}
    return Expr(:(::), root_name, Expr(:curly, :Array, to_expr(T, root_name, exprs), 1))
end

function to_field_expr(::Type{T}, root_name, exprs) where {T}
    Expr(:(::), root_name, to_expr(T, root_name, exprs))
end

str_code = """
struct MyStruct
    a::Int
    b::OtherType
    c::Vector{Int,1}
end
"""


text:
function FT_forward(alpha,beta,c0,c1,c2,c3,T,times,rmr0,cf_irad=false)
    ind_oldest = 1
    n_times = length(times)
    dt_equiv = 0.0
    rcb2_base = 0.0
    Tvv = eltype(T)
    N_t_segs = length(T)+1
    rho_r = zeros(Tvv,(N_t_segs,N_t_segs))
    rho_rm = zeros(Tvv,(N_t_segs))
    r_m = zeros(Tvv,(N_t_segs))
    prior = zeros(Tvv,(N_t_segs))
    rho_r_orig = zeros(Tvv,(N_t_segs))
    rcb2 = zeros(Tvv,(N_t_segs))
    taus = zeros(Tvv,(N_t_segs))

    kappa = 1.04 - rmr0
    rdaam2nd_root = 0.0#0.5274 # cutoff bellow which rho_r values can be negative
    rdaam2nd_root_cutoff = (rdaam2nd_root^(1/kappa))*(1-rmr0)+rmr0 # also prevents
    min_red_tracks = 0.13 # after Ketcham, 2000; no tracks shorter than 2.18 um
    for ind in eachindex(times)
        if ind<n_times
            t_dam = times[n_times-ind]-times[n_times-ind+1]+dt_equiv
            rcb2_base = (c0 + c1*(log(t_dam)-c2)/(log((1.0/T[n_times-ind]))-c3))
            rcb2_fill_val = ((rcb2_base^(1.0/alpha))+1)^(-1)
            rcb2[ind] = rcb2_fill_val

            if rcb2[ind]>=rdaam2nd_root_cutoff
                rho_rm[n_times-ind] = ((rcb2_fill_val-rmr0)/(1.0-rmr0))^kappa
            else
                rho_rm[n_times-ind] = 0.0
            end

            rho_r_mean = rho_rm[n_times-ind]#(rho_rm[n_times-ind]+rho_rm[n_times-ind+1])/2

            if cf_irad==true
                r_m[n_times-ind+1] = 1.396*rho_r_mean - 0.4017
            else
                r_m[n_times-ind+1] = -1.499*rho_r_mean^2 + 4.150*rho_r_mean - 1.656
            end
            prior[n_times-ind] = (exp(lambda_38*times[n_times-ind]) -
                exp(lambda_38*times[n_times-ind+1]))/lambda_38
            if (ind>1) & (rho_rm[n_times-ind+1]>0.0)
                if rho_r_mean>=0.765
                    rho_r[n_times-ind,n_times] = 1.6*rho_r_mean - 0.6
                elseif rho_r_mean<0.765
                    rho_r[n_times-ind,n_times] = 9.205*rho_r_mean^2 - 9.157*rho_r_mean + 2.269
                end

                rho_r_orig[n_times-ind+1] = rho_r[n_times-ind,n_times]
                taus[n_times-ind+1] += (1.0./(0.893))*(rho_r[n_times-ind,n_times]*
                    (times[n_times-ind]-times[n_times-ind+1]))
            elseif (ind>1) & iszero(ind_oldest)
                ind_oldest = n_times - ind + 1
            end

        else
            if (ind>1)
                rho_r_mean = rho_rm[n_times-1]

                if rho_r_mean>=0.765
                    rho_r = 1.6*rho_r_mean - 0.6
                elseif rho_r_mean<0.765
                    rho_r = 9.205*rho_r_mean^2 - 9.157*rho_r_mean + 2.269
                end
                rho_r_orig[end] = rho_r

                taus[end] += 0.5*(1.0./(0.893))*(rho_r*
                    (times[n_times-1]-times[n_times]))

            end
        end
              # use concept of equivalent time

        if ind<(n_times-1)
        dt_equiv = exp((((((1.0/rcb2[ind])-1.0)^alpha-c0)*(log(1.0/T[n_times-ind-1])-c3))/c1)+c2)
        end
    end

    return times[ind_oldest],taus, r_m, rho_rm, rcb2, rho_r_orig, prior
end

function FT_forward_age_distribution()
    return nothing
end


text:
using Pkg, DataFrames, Distributions, DelimitedFiles, CSV
using JuMP, Gurobi, KNITRO
Solver = "KNITRO"
setT = 1:24

κ = 0.7
testsystemT = "NYISO_v2_70_rCoal"
testsystemD = "Manhattan"

σ0_pv = 20e-2
σ0_wind = 20e-2/2
υ0 = 0.0
UserGap = 1e-6
D_ratio = 0.2
D_scale = 1
D_new = 1
busT_gep = 10
busD_gep = 4
SD_tol = 0
DCRF = round(0.0002198427046320309 * length(setT) / 24, digits=5)
ScaleObj = 1

semi_cut_MP = 0
semi_cut_SP = 0
round_MP = 1
round_SP = 0

policy_switch = 1
IA_switch = 1
IA_LMP = 0
WMP = 0
CTI = 0
ED =  50
γCT = 1
WHA = 0

η = 0.03

include("StructType_v2.jl")
bound = bounds(0, 200, 0, 200, 0, 200, 0, 200, 0, 0)

include("MPini.jl")
include("MP.jl")
include("SP.jl")

include("NetworkDataType.jl")
include("NetworkLoad_v3.jl")
include("Demand24Load_NYISO.jl")
filename_Load24 = string(pwd(), "/data/load24data-NYISO.csv")
temp_loadcurve = Demand24Load(filename_Load24)
loadcurveT = D_scale * temp_loadcurve'
loadcurveD = D_scale * [0.5915717782270937, 0.5603916143450037, 0.5383633074385442, 0.524103437267213, 0.5188890071299351, 0.5265510269234862, 0.5592210279876556, 0.6044482281579228, 0.6518037671597319, 0.6848994359902096, 0.7139512610407577, 0.7381079067787591, 0.7567308715547515, 0.7688624028945408, 0.7789720123443652, 0.7858891135468766, 0.7866340321379164, 0.7777003298925189, 0.7514100244758966, 0.7289560498031287, 0.7203362775353837, 0.7102266680855592, 0.6857507715228265, 0.6504203469192296]

busesT, linesT, generatorsT, datamatT, setL_T, setN_T, DTp, DTq = NetworkLoad(testsystemT)
busesD, linesD, generatorsD, datamatD, setL_D, setN_D, DDp, DDq = NetworkLoad(testsystemD)
Γpf = map((x) -> if isnan(x) 0 else x end, DDq./DDp)

include("SetGen_v3.jl")
wind_LBsetT, wind_OSsetT, pvsetT, pv_BTMsetT, hydrosetT, nucsetT, coalsetT, oilsetT, gassetT, B_gT, B_gnT, setI_TS, setI_TS_hat, setI_RT, setI_RT_hat, setI_CT, setI_CT_hat = SetGen(generatorsT, busesT, busT_gep)
wind_LBsetD, wind_OSsetD, pvsetD, pv_BTMsetD, hydrosetD, nucsetD, coalsetD, oilsetD, gassetD, B_gD, B_gnD, setI_DS, setI_DS_hat, setI_RD, setI_RD_hat, setI_CD, setI_CD_hat = SetGen(generatorsD, busesD, busD_gep)

include("GenCurveLoad.jl")
include("GenParam_v3.jl")
CinvgT, CopergT, βFC_T, σT, υT = GenParam(generatorsT, setI_TS, setI_TS_hat, wind_LBsetT, wind_OSsetT, pvsetT, pv_BTMsetT, hydrosetT, nucsetT, coalsetT, oilsetT, gassetT, σ0_pv, σ0_wind, υ0)
CinvgD, CopergD, βFC_D, σD, υD = GenParam(generatorsD, setI_DS, setI_DS_hat, wind_LBsetD, wind_OSsetD, pvsetD, pv_BTMsetD, hydrosetD, nucsetD, coalsetD, oilsetD, gassetD, σ0_pv, σ0_wind, υ0)

wind_LBsetD = map((x) -> x + 1000, wind_LBsetD)
wind_OSsetD = map((x) -> x + 1000, wind_OSsetD)
pvsetD = map((x) -> x + 1000, pvsetD)
pv_BTMsetD = map((x) -> x + 1000, pv_BTMsetD)
hydrosetD = map((x) -> x + 1000, hydrosetD)
nucsetD = map((x) -> x + 1000, nucsetD)
coalsetD = map((x) -> x + 1000, coalsetD)
oilsetD = map((x) -> x + 1000, oilsetD)
gassetD = map((x) -> x + 1000, gassetD)
setI_DS = map((x) -> x + 1000, setI_DS)
setI_DS_hat = map((x) -> x + 1000, setI_DS_hat)
setI_RD = map((x) -> x + 1000, setI_RD)
setI_RD_hat = map((x) -> x + 1000, setI_RD_hat)
setI_CD = map((x) -> x + 1000, setI_CD)
setI_CD_hat = map((x) -> x + 1000, setI_CD_hat)

Pmax = merge(Dict(i => generatorsT[i].Pmax for i = 1:length(generatorsT)),
            Dict(i+1000 => generatorsD[i].Pmax for i = 1:length(generatorsD)))
Qmax = merge(Dict(i => generatorsT[i].Qmax for i = 1:length(generatorsT)),
            Dict(i+1000 => generatorsD[i].Qmax for i = 1:length(generatorsD)))

B_gT = Dict(i => B_gT[i] for i = 1:length(B_gT))
B_gD = Dict(i + 1000 => B_gD[i] for i = 1:length(B_gD))
B_g = merge(B_gT, B_gD)
B_gnD = merge(Dict(i => B_gnD[i] + 1000 * ones(Int64, length(B_gnD[i]),1) for i = 1:length(B_gnD) if ~isempty(B_gnD[i])),
   Dict(i => B_gnD[i] for i = 1:length(B_gnD) if isempty(B_gnD[i])))

setI_R = union(setI_RT, setI_RD)
setI_R_hat = union(setI_RT_hat, setI_RD_hat)
setI_C = union(setI_CT, setI_CD)
setI_C_hat = union(setI_CT_hat, setI_CD_hat)
setI_T = union(setI_TS, setI_DS)
setI_T_hat = union(setI_TS_hat, setI_DS_hat)
setI_D = setI_DS
setI_D_hat = setI_DS_hat
setI_DER = union(setI_DS, B_gnT[busT_gep])
setI_DER_hat = union(setI_DS_hat, intersect(setI_T_hat, B_gnT[busT_gep]))

CCT = 50
H = Dict(i =>  if i in union(gassetT, gassetD) 0.419
        elseif i in union(coalsetT, coalsetD) 1.002
        elseif i in union(oilsetT, oilsetD) 0.9606
        else 0 end
    for i in setI_C)

Cinvg = merge(Dict(i => CinvgT[i] for i = 1:length(CinvgT)), Dict(i + 1000 => CinvgD[i] for i = 1:length(CinvgD)))
Coperg = merge(Dict(i => CopergT[i] for i = 1:length(CopergT)), Dict(i + 1000 => CopergD[i] for i = 1:length(CopergD)))
βFC = merge(Dict([i,t] => βFC_T[i,t] for i in 1:size(βFC_T,1), t in setT),
            Dict([i + 1000, t] => βFC_D[i,t] for i in 1:size(βFC_D,1), t in setT))
σ = merge(Dict(i => σT[i] for i = 1:length(σT)), Dict(i + 1000 => σD[i] for i = 1:length(σD)))
υ = merge(Dict([i,t] => υT[i,t] for i = 1:size(υT,1), t in setT),
         Dict([i + 1000,t] => υD[i,t] for i = 1:size(υD,1), t in setT))

spi_ratio = 0.05
π0 = 50
s0 = spi_ratio * π0
d0 = DDp * loadcurveD'
N = 0.25
M = N * d0 .+ π0 .+ s0
M_new = 11(π0+s0)
N_new = 10(π0+s0)ones(size(d0))./d0
ST_res = s0 * ones(length(setN_T), maximum(setT)) # surcharge for transmission nodes

DTp_inter = DTp[busT_gep] * loadcurveT[busT_gep,:]
M_new = 11(π0+s0)
N_new_DTp = 10(π0+s0)ones(size(DTp_inter))./DTp_inter

time_M = @elapsed tempM = MPini(testsystemT, testsystemD, busT_gep, busD_gep, κ, setT, η, σ0_pv, σ0_wind, υ0, Solver, bound, D_ratio, ScaleObj, SD_tol, M, N, ED, γCT, ST_res, WMP, CTI)
if round_SP == 1
   time_S = @elapsed tempS = SP(testsystemT, testsystemD, busT_gep, busD_gep, κ, setT, η, σ0_pv, σ0_wind, υ0,
      round.(value.(tempM[:πDER]), digits=5), round.(value.(tempM[:πC]), digits=5), round.(value.(tempM[:s]), digits=5),
      round.(value.(tempM[:τe]), digits=5), round.(value.(tempM[:τc]), digits=5),
      Solver, D_ratio, ScaleObj, SD_tol, M, N, ED, γCT, ST_res, WMP, CTI)
else
   time_S = @elapsed tempS = SP(testsystemT, testsystemD, busT_gep, busD_gep, κ, setT, η, σ0_pv, σ0_wind, υ0, value.(tempM[:πDER]), value.(tempM[:πC]), value.(tempM[:s]),
      value.(tempM[:τe]), value.(tempM[:τc]), Solver, D_ratio, ScaleObj, SD_tol, M, N, ED, γCT, ST_res, WMP, CTI)
end
Niter = 100
conv = 777*ones(Niter,4)
time_iter = zeros(Niter,2)
MPbound_iter = zeros(Niter)
SPbound_iter = zeros(Niter)

πDER_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempM[:πDER])}}(undef,Niter)
πC_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempM[:πC])}}(undef,Niter)
s_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempM[:s])}}(undef,Niter)
τe_iter = 777*ones(Niter)
τc_iter = 777*ones(Niter)
ObjReg_M_iter = 777*ones(Niter)
ObjDER_M_iter = 777*ones(Niter)
ObjDM_M_iter = 777*ones(Niter)
ObjWM_M_iter = 777*ones(Niter)
ObjReg_S_iter = 777*ones(Niter)
ObjDER_S_iter = 777*ones(Niter)
ObjDM_S_iter = 777*ones(Niter)
ObjWM_S_iter = 777*ones(Niter)

gpmax_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:gpmax])}}(undef,Niter)
gTpo_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:gTpo])}}(undef,Niter)
gDpo_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:gDpo])}}(undef,Niter)
gDqo_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:gDqo])}}(undef,Niter)
gTp_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:gTp])}}(undef,Niter)
gDp_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:gDp])}}(undef,Niter)
gDq_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:gDq])}}(undef,Niter)
λT_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:λT])}}(undef,Niter)
λD_iter = Vector{JuMP.Containers.DenseAxisArray{Float64,ndims(tempS[:λD])}}(undef,Niter)
status_iter = Array{Union{Nothing, String}}(nothing, Niter, 2)

ii = 1
time_iter[ii,1] = time_M
time_iter[ii,2] = time_S
status_iter[ii,1] = string(termination_status(tempM))
status_iter[ii,2] = string(termination_status(tempS))

πDER_iter[ii] = value.(tempM[:πDER])
πC_iter[ii] = value.(tempM[:πC])
s_iter[ii] = value.(tempM[:s])
τe_iter[ii] = value.(tempM[:τe])
τc_iter[ii] = value.(tempM[:τc])
ObjReg_M_iter[ii] = value.(tempM[:Obj_REG])
ObjDER_M_iter[ii] = value.(tempM[:Obj_DER])
ObjDM_M_iter[ii] = sum(value.(tempM[:Obj_D_prim_act]))
ObjWM_M_iter[ii] = sum(value.(tempM[:Obj_W_prim]))
ObjReg_S_iter[ii] = value.(tempS[:Obj_REG])
ObjDER_S_iter[ii] = value.(tempS[:Obj_DER])
ObjDM_S_iter[ii] = sum(value.(tempS[:Obj_D_prim_act]))
ObjWM_S_iter[ii] = sum(value.(tempS[:Obj_W_prim]))

if round_MP == 1
   gpmax_iter[ii] = abs.(round.(value.(tempS[:gpmax]), digits=5))
   gTpo_iter[ii] = abs.(round.(value.(tempS[:gTpo]), digits=5))
   gDpo_iter[ii] = abs.(round.(value.(tempS[:gDpo]), digits=5))
   gDqo_iter[ii] = abs.(round.(value.(tempS[:gDqo]), digits=5))
else
   gpmax_iter[ii] = abs(value.(tempS[:gpmax]))
   gTpo_iter[ii] = abs(value.(tempS[:gTpo]))
   gDpo_iter[ii] = abs(value.(tempS[:gDpo]))
   gDqo_iter[ii] = abs(value.(tempS[:gDqo]))
end

gTp_iter[ii] = value.(tempS[:gTp])
gDp_iter[ii] = value.(tempS[:gDp])
gDq_iter[ii] = value.(tempS[:gDq])
λT_iter[ii] = value.(tempS[:λT])
λD_iter[ii] = value.(tempS[:λD])

conv_temp = (value.(tempM[:Obj_REG]) - value.(tempS[:Obj_REG]))
global conv[ii,:] = [round.(value.(tempM[:Obj_REG]), digits = 3), round.(value.(tempS[:Obj_REG]), digits = 3),
   round.(value.(tempM[:Obj_REG]) - value.(tempS[:Obj_REG]), digits = 3), (value.(tempM[:Obj_REG]) - value.(tempS[:Obj_REG])) / value.(tempM[:Obj_REG])]
global SPbound_iter[ii] = value.(tempS[:SPbound])
global MPbound_iter[ii] = value.(tempM[:MPbound])
df = DataFrame(τc = round.(τc_iter[1:ii], digits = 3), τe = round.(τe_iter[1:ii], digits = 3),
   ObjReg_M = round.(ObjReg_M_iter[1:ii], digits=3),
   ObjDER_M = round.(ObjDER_M_iter[1:ii], digits=3),
   ObjDM_M = round.(ObjDM_M_iter[1:ii], digits=3),
   ObjWM_M = round.(ObjWM_M_iter[1:ii], digits=3),
   ObjReg_S = round.(ObjReg_S_iter[1:ii], digits=3),
   ObjDER_S = round.(ObjDER_S_iter[1:ii], digits=3),
   ObjDM_S = round.(ObjDM_S_iter[1:ii], digits=3),
   ObjWM_S = round.(ObjWM_S_iter[1:ii], digits=3)
   )

for ii in 2:Niter
   time_M = @elapsed global tempM = MP(testsystemT, testsystemD, busT_gep, busD_gep, κ, setT, η, σ0_pv, σ0_wind, υ0, Solver, bound, D_ratio, ScaleObj, SD_tol, M, N, ED, γCT, ST_res, WMP, CTI, gpmax_iter, gTpo_iter, gDpo_iter, gDqo_iter, ii-1)
   global time_iter[ii,1] = time_M
   global status_iter[ii,1] = string(termination_status(tempM))

   global πDER_iter[ii] = value.(tempM[:πDER])
   global πC_iter[ii] = value.(tempM[:πC])
   global s_iter[ii] = value.(tempM[:s])
   global τe_iter[ii] = value.(tempM[:τe])
   global τc_iter[ii] = value.(tempM[:τc])
   global ObjReg_M_iter[ii] = value.(tempM[:Obj_REG])
   global ObjDER_M_iter[ii] = value.(tempM[:Obj_DER])
   global ObjDM_M_iter[ii] = sum(value.(tempM[:Obj_D_prim_act]))
   global ObjWM_M_iter[ii] = sum(value.(tempM[:Obj_W_prim]))

   time_S = @elapsed global tempS = SP(testsystemT, testsystemD, busT_gep, busD_gep, κ, setT, η, σ0_pv, σ0_wind, υ0, value.(tempM[:πDER]), value.(tempM[:πC]), value.(tempM[:s]),value.(tempM[:τe]), value.(tempM[:τc]), Solver, D_ratio, ScaleObj, SD_tol, M, N, ED, γCT, ST_res, WMP, CTI)

   global time_iter[ii,2] = time_S
   global status_iter[ii,2] = string(termination_status(tempS))

   global ObjReg_S_iter[ii] = value.(tempS[:Obj_REG])
   global ObjDER_S_iter[ii] = value.(tempS[:Obj_DER])
   global ObjDM_S_iter[ii] = sum(value.(tempS[:Obj_D_prim_act]))
   global ObjWM_S_iter[ii] = sum(value.(tempS[:Obj_W_prim]))
   global gpmax_iter[ii] = abs(value.(tempS[:gpmax]))
   global gTpo_iter[ii] = abs(value.(tempS[:gTpo]))
   global gDpo_iter[ii] = abs(value.(tempS[:gDpo]))
   global gDqo_iter[ii] = abs(value.(tempS[:gDqo]))

   global gTp_iter[ii] = value.(tempS[:gTp])
   global gDp_iter[ii] = value.(tempS[:gDp])
   global gDq_iter[ii] = value.(tempS[:gDq])
   global λT_iter[ii] = value.(tempS[:λT])
   global λD_iter[ii] = value.(tempS[:λD])

   conv_temp = (value.(tempM[:Obj_REG]) - value.(tempS[:Obj_REG]))
   global conv[ii,:] = [round.(value.(tempM[:Obj_REG]), digits = 3), round.(value.(tempS[:Obj_REG]), digits = 3),
      round.(value.(tempM[:Obj_REG]) - value.(tempS[:Obj_REG]), digits = 3), (value.(tempM[:Obj_REG]) - value.(tempS[:Obj_REG])) / value.(tempM[:Obj_REG])]
   global SPbound_iter[ii] = value.(tempS[:SPbound])
   global MPbound_iter[ii] = value.(tempM[:MPbound])
   df = DataFrame(τc = round.(τc_iter[1:ii], digits = 3), τe = round.(τe_iter[1:ii], digits = 3),
      ObjReg_M = round.(ObjReg_M_iter[1:ii], digits=3),
      ObjDER_M = round.(ObjDER_M_iter[1:ii], digits=3),
      ObjDM_M = round.(ObjDM_M_iter[1:ii], digits=3),
      ObjWM_M = round.(ObjWM_M_iter[1:ii], digits=3),
      ObjReg_S = round.(ObjReg_S_iter[1:ii], digits=3),
      ObjDER_S = round.(ObjDER_S_iter[1:ii], digits=3),
      ObjDM_S = round.(ObjDM_S_iter[1:ii], digits=3),
      ObjWM_S = round.(ObjWM_S_iter[1:ii], digits=3)
   )
end

max_stars_repo_name: JipKim/IA

max_issues_repo_name: JipKim/IA

max_forks_repo_name: JipKim/IA


text:
module SSLTestClip
using SimpleSDMLayers
using Test

M = rand(Bool, (10, 10))
S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)

cl1 = clip(S; left=0.2, right=0.6, bottom=0.5, top=1.0)
@test typeof(cl1) == typeof(S)
@test cl1.top ≈ 1.0
@test cl1.bottom ≈ 0.5
@test cl1.right ≈ 0.6
@test cl1.left ≈ 0.2
@test clip(S; left=0.19).left <= 0.2

cl2 = clip(S; left=0.2, bottom=0.5)
@test typeof(cl2) == typeof(S)
@test cl2.top ≈ 1.0
@test cl2.bottom ≈ 0.5
@test cl2.right ≈ 1.0
@test cl2.left ≈ 0.2

end


text:
#In this file I prototype a simple hospital model
#Only input is number and age of each new hospitalised person on each day
#Only looking at Hosp -> ICU and ICU -> death
include("hospitalisation_data.jl");
function total_incidence_and_cum_incidence_for_each_sim(sims)
    n = length(sims.u)
    T = length(sims.u[1])
    cum_incidence_for_each_sim = zeros(n,T)
    for k = 1:n,t = 1:T
        cum_incidence_for_each_sim[k,t] = sum(sims.u[k][t][:,:,1:3])
    end
    return cum_incidence_for_each_sim,diff(cum_incidence_for_each_sim,dims = 2)
end

function find_peak_time_by_sim(sims)
    n = length(sims.u)
    T = length(sims.u[1])
    cum_incidence_for_each_sim = zeros(n,T)
    for k = 1:n,t = 1:T
        cum_incidence_for_each_sim[k,t] = sum(sims.u[k][t][:,:,1:3])
    end
    incidence = diff(cum_incidence_for_each_sim,dims = 2)
    return [findmax(incidence[k,:]) for k = 1:n]
end

function find_peaks(incidence_array)
    n,T = size(incidence_array)
    return [findmax(incidence_array[k,:])[1] for k = 1:n]
end

function cum_incidence_for_each_sim_by_type(sims)
    n = length(sims.u)
    T = length(sims.u[1])
    cum_incidence_A = zeros(n,T)
    cum_incidence_M = zeros(n,T)
    cum_incidence_V = zeros(n,T)
    cum_incidence_H = zeros(n,T)

    for k = 1:n,t = 1:T
        cum_incidence_A[k,t] = sum(sims.u[k][t][:,:,1])
        cum_incidence_M[k,t] = sum(sims.u[k][t][:,:,2])
        cum_incidence_V[k,t] = sum(sims.u[k][t][:,:,3])
        cum_incidence_H[k,t] = sum(sims.u[k][t][:,:,4])
    end
    return cum_incidence_A,cum_incidence_M,cum_incidence_V,cum_incidence_H
end
"""
generate_hospitalisation_outcomes(H)
"""
death_rate_critical_cases = mean([0.5,0.75])

function split_incoming_hospitalisations_into_outcomes(H,age,area)
    T = length(H)
    #randomly split population into hospital and ICU
    new_hosp = zeros(Int64,T);new_ICU = zeros(Int64,T); new_to_die = zeros(Int64,T);
    for t = 1:T
        outcomes = rand(Multinomial(Int64.(H[t]),[1 - ICU_rate_by_age_cond_hosp[age],
                                            ICU_rate_by_age_cond_hosp[age]*(1-death_rate_critical_cases),
                                            ICU_rate_by_age_cond_hosp[age]*death_rate_critical_cases]) )
        new_hosp[t] = outcomes[1]
        new_ICU[t] = outcomes[2]
        new_to_die[t] = outcomes[3]
    end
    return new_hosp,new_ICU,new_to_die
end

function generate_hospitalisation_outcomes(H,age,area)
    T = length(H)
    hosp_occup = zeros(Int64,T)
    ICU_occup = zeros(Int64,T)
    death_incidence = zeros(Int64,T)
    #randomly split population into hospital and ICU
    new_hosp,new_ICU,new_to_die = split_incoming_hospitalisations_into_outcomes(H,age,area)
    #Assign population outcomes
    for t = 1:T
        for k = 1:new_hosp[t]
            tₑ = rand(6:10)
            timesinhosp = t:(min(T,t+tₑ))#6-10 days in Hospital then discharge
            hosp_occup[timesinhosp] .+= 1
        end
        for k = 1:new_ICU[t]
            tₑ = rand(7:14)
            timesinICU = t:(min(T,t+tₑ))#7-14 days in ICU then discharge
            ICU_occup[timesinICU] .+= 1
        end
        for k = 1:new_to_die[t]
            tₑ = rand(7:14)
            timesinICU = t:(min(T,t+tₑ))#7-14 days in ICU then death
            ICU_occup[timesinICU] .+= 1
            death_incidence[min(T,t+tₑ)] += 1
        end
    end

    return hosp_occup, ICU_occup, new_ICU,death_incidence
end

function total_hospital_outcomes_per_sim(sims)
    n = length(sims.u)
    T = length(sims.u[1])
    hosp_occup_per_sim = zeros(Int64,n,T-1)
    ICU_occup_per_sim = zeros(Int64,n,T-1)
    new_ICU_per_sim = zeros(Int64,n,T-1)
    death_incidence_per_sim = zeros(Int64,n,T-1)
    for k = 1:n
        for a = 1:17
            H = diff([sum(sims.u[k][t][:,a,4]) for t = 1:366])
            hosp_occup, ICU_occup, new_ICU,death_incidence = generate_hospitalisation_outcomes(H,a,1)
            hosp_occup_per_sim[k,:] .+= hosp_occup
            ICU_occup_per_sim[k,:] .+= ICU_occup
            new_ICU_per_sim[k,:] .+= new_ICU
            death_incidence_per_sim[k,:] .+= death_incidence
        end
    end
    return hosp_occup_per_sim,ICU_occup_per_sim,new_ICU_per_sim,death_incidence_per_sim
end