text:
function func_U(t)
    e = exp(1)
    return 2.0*((e^(3.0*t))-(e^3))/((e^(3.0*t/2.0))*(2.0+(e^3)))
end

function func(X,t)
    F = [0.5*X[1]+ func_U(t), X[1]*X[1]+0.5*(func_U(t)*func_U(t))]
    return F
end

function func_adj(BF, X)
    BX = [0.5*BF[1]+2.0*X[1]*BF[2], 0.0]
    return BX
end

function advance(F_H,t,h)
    k0 = func(F_H,t)
    G = [F_H[1] + h/2.0*k0[1], F_H[2] + h/2.0*k0[2]]
    k1 = func(G,t+h/2.0)
    F = [F_H[1] + h*k1[1], F_H[2] + h*k1[2]]
    return F
end

function adjoint(F_H,L_H,F,L,t,h)
    k0 = func(F_H,t)
    G = similar(F_H)
    Bk1 = similar(L_H)
    G = [F_H[1] + h/2.0*k0[1], F_H[2] + h/2.0*k0[2]]
    k1 = func(G,t+h/2.0)
    L = [L_H[1], L_H[2]]
    Bk1 = [h*L_H[1], h*L_H[2]]
    BG = func_adj(Bk1,G)
    Bk0 = similar(BG)
    L = [L[1] + BG[1], L[2] + BG[2]]
    Bk0 = [h/2.0*BG[1], h/2.0*BG[2]]
    BH = func_adj(Bk0,F_H)
    L = [L[1] + BH[1], L[2] + BH[2]]
    return L
end

function opt_sol(Y,t)
    e = exp(1)
    Y[1] = (2.0*e^(3.0*t)+e^3)/(e^(3.0*t/2.0)*(2.0+e^3))
    Y[2] = (2.0*e^(3.0*t)-e^(6.0-3.0*t)-2.0+e^6)/((2.0+e^3)^2)
    return
end

function opt_lambda(L,t)
    e = exp(1)
    L[1] = (2.0*e^(3-t)-2.0*e^(2.0*t))/(e^(t/2.0)*(2+e^3))
    L[2] = 1.0
    return
end


text:
# Create example data
x = LinRange(-2, 2, 10)
y = x.^3 .+ x.^2 .+ x .+ 6
err = LinRange(0.5, 3, 10)
# Draw symmetric, horizontal error bars
errorbar(x, y, err, horizontal=true)
# Draw asymmetric error bars with markers
errorbar(x, y, err, err./2, "-o")

ext: jl


text:
######################################################################
# examples.jl
#
# This used to be in constructors of ergodic managers
# It added a lot of clutter and I didn't like it there
######################################################################
function ErgodicManagerR2(example_name::String; K::Int=5, bins::Int=100)
	L = 1.0
	d = Domain([0.,0.], [L, L], [bins,bins])

	if example_name == "single gaussian"
		mu = [L/2.0, L/2.0]
		Sigma = 0.03 * [1 0; 0 1]
		phi = gaussian(d, mu, Sigma)
		return ErgodicManagerR2(d, phi, K)
	elseif example_name == "double gaussian"
		mu1 = [0.3, 0.7]
		Sigma1 = 0.025 * [1 0; 0 1]
		mu2 = [0.7, 0.3]
		Sigma2 = 0.025 * [1 0; 0 1]
		#phi = gaussian(d, [mu1,mu2], [Sigma1, Sigma2], [.5,.5])
		phi = gaussian(d, [mu1,mu2], [Sigma1, Sigma2])
		return ErgodicManagerR2(d, phi, K)
	else
		error("example name not recognized")
	end
end

function ErgodicManagerSE2(example_name::String; K::Int=5,bins::Int=50)
	L = 1.0
	d = Domain([0.,0.,-pi], [L,L,pi], [bins,bins,bins])
	println("h2")
	println("h3")

	if example_name == "single gaussian"
		mu = [L/2.0, L/2.0, 0]
		Sigma = 0.03 * [1 0 0; 0 1 0; 0 0 1]
		phi = gaussian(d, mu, Sigma)
		return ErgodicManagerSE2(d, phi, K)

	elseif example_name == "double gaussian"
		# How I should do it...
		mu1 = [0.3, 0.7]
		Sigma1 = 0.025* [1 0; 0 1]
		mu2 = [0.7, 0.3]
		Sigma2 = 0.025* [1 0; 0 1]
		phi = gaussian(d, mu1, Sigma1, mu2, Sigma2)
		return ErgodicManagerR2(d, phi, K)

		# Create Gaussian distribution and its coefficients
		#mu1 = [0.3, 0.7]
		#Sigma1 = 0.025* [1 0; 0 1]
		#mu2 = [0.7, 0.3]
		#Sigma2 = 0.025* [1 0; 0 1]
		#phi!(em, mu1, Sigma1, mu2, Sigma2)
		#decompose!(em)
	else
		error("example name not recognized")
	end
	return em
end


text:
using CUTEst
using JSOSolvers
using Projeto2Solvers
using LinearAlgebra
using NLPModels
using Plots
using SolverBenchmark
using JuMP
using Ipopt
# using ForwardDiff

#%%
function runcutest()

  #pnames = CUTEst.select(max_var=2, max_con=0, only_free_var=true)
  # pnames = CUTEst.select(max_var=5, max_con=0, only_free_var=true)

  pnames = CUTEst.select(min_var=1, max_var = 1000,
                          max_con=0, only_free_var=true)

  # pnames = CUTEst.select(min_var=501, max_var = 1000,
  #                         max_con=0, only_free_var=true)

  sort!(pnames)
  problems = (CUTEstModel(p) for p in pnames) # Generator of problems

  # If you need to define different arguments, you can wrap
  trunk_wrapper(nlp; kwargs...) = trunk(nlp, max_time=10.0; kwargs...)
  my_newton_wrapper(nlp; kwargs...) = newton(nlp, max_time=10.0; kwargs...)
  my_gradiente_wrapper(nlp; kwargs...) = gradiente(nlp, max_time=10.0; kwargs...)

  solvers = Dict(
    :lbfgs => lbfgs,
    # :trunk => trunk_wrapper,
    # :bfsg_bl => bfgs_bl,
    # :newton => my_newton_wrapper,
    :newcombusca => newtoncombusca,
    # :LBFGS_StrongWolfe => LBFGS_StrongWolfe,
    # :Newton_rc_bissec => Newton_rc_bissec,
    #:bfgs_rc => bfgs_rc,
    :l_bfgs_rcst => l_bfgs_rcst
    # :gradiente => my_gradiente_wrapper

    )

  stats = bmark_solvers(solvers, problems)
end

#%%
function table_and_plots(stats)
  cols = [:name, :nvar, :status, :objective, :dual_feas, :elapsed_time, :neval_obj, :neval_grad, :neval_hess]
  for (k,v) in stats
    open("tabelas/tabela-$k.md", "w") do io
      pretty_stats(io, v[:,cols])
    end
    open("tabelas/tabela-$k.tex", "w") do io
      pretty_latex_stats(io, v[:,cols])
    end
  end

  # Performance profile

  # First, define metrics
  not_first_order(df) = df.status .!= :first_order
  Fx = hcat([df.objective for df in values(stats)]...)
  Fx[isnan.(Fx)] .= Inf
  fmin = minimum(Fx, dims=2)
  not_best_fx(df) = df.dual_feas .> fmin * 1.01 .+ 1e-6 # 1% tolerance, with safeguard for 0.
  costs = [
    df -> df.elapsed_time,
    df -> not_first_order(df) * Inf + df.elapsed_time,
    df -> not_best_fx(df) * Inf + df.elapsed_time,
    df -> not_first_order(df) * Inf + df.neval_obj + df.neval_grad + df.neval_hess + df.neval_hprod
  ]
  costnames = [
    "Pure Δt",
    "Δt checking 1st order",
    "Δt checking fx ≈ fmin",
    "Eval checking 1st order"
  ]
  p = profile_solvers(stats, costs, costnames)
  png(p, "plots/aaprofile1000")
end

#%%
# You can comment out this after the tests are run
stats = runcutest()
p = table_and_plots(stats)

#%%
# Salvar resultados em data frame - MUDAR OS NOMES
CSV.write("todos-1000.csv", stats)
CSV.write("ncb-1000.csv", stats[:newcombusca])
CSV.write("lst-1000.csv", stats[:l_bfgs_rcst])
CSV.write("lbfgs-1000.csv", stats[:lbfgs])


text:
# This file reads the career paths from the excel file. Career paths developed by Guy

#XlsxFile = "C:/Users/Administrator/Dropbox/JuliaManpowerPlanning/GuyCareerPathsAllPop.xlsx"

requiredTypes = [ "AcademicLevel", "Affiliation", "Assignment", "CareerStatus", "Job", "Rank", "SubJob", "CareerPath" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( dirname( Base.source_path() ), "..", "functions", reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end

if !isdefined(:Taro)
    using Taro, JavaCall
    Taro.init()
end

if !isdefined(:CPworkbook)
    CPworkbook = Workbook(XlsxFile)
    GIsheet = getSheet(CPworkbook, "GeneralInfo")
    CPsheet = getSheet(CPworkbook, "CareerPaths")
    CPDurationsheet = getSheet(CPworkbook, "CPDurations")
    CPAttritionsheet = getSheet(CPworkbook, "AttritionRate")
    InitMPsheet = getSheet(CPworkbook, "InitMP")
    Objsheet = getSheet(CPworkbook, "Objective")
    RecruitmentSheet = getSheet(CPworkbook, "Recruitment")
    SubpopulationSheet = getSheet(CPworkbook, "SubPopulations")
    DependenciesSheet = getSheet(CPworkbook, "CPDependecies")
end

InitManpower = Vector{InitMPcluster}()

MPObjectives = Vector{ManpowerObjective}()

    # Reading General Information __________________________________________________
SimulationName = getCellValue(getCell(getRow(GIsheet, 0), 1))
SaveInputs = (getCellValue(getCell(getRow(GIsheet, 1), 1)) == "Yes")
IntegerSolution = (getCellValue(getCell(getRow(GIsheet, 2), 1)) == "Yes")
Tolerances_MIPGap = getCellValue(getCell(getRow(GIsheet, 3), 1))
SimulationYear = Int(getCellValue(getCell(getRow(GIsheet, 4), 1)))
PlottingResults = (getCellValue(getCell(getRow(GIsheet, 5), 1)) == "Yes")
AgeDistPlot = (getCellValue(getCell(getRow(GIsheet, 6), 1)) == "Yes")
println("Simulation : $SimulationName, saving input data : $SaveInputs.")

    # Reading Sub populations __________________________________________________

SubPopulations = Vector{Subpopulation}()

subpopnumber = jcall(SubpopulationSheet, "getLastRowNum", jint, ())

for i in 0:subpopnumber
    Row = getRow(SubpopulationSheet, i)
    push!(SubPopulations, Subpopulation( nm = getCellValue(getCell(Row, 0))))
end

    # Reading Career paths _________________________________________________________

    CPnumber = jcall(CPsheet, "getLastRowNum", jint, ()) + 1

    GuyCareerPaths = Vector{CareerPath}()

    for i in 0:CPnumber-1
        WorkingCP = getRow(CPsheet, i)

        for aff in 1:length(SubPopulations)
            if getCellValue(getCell(WorkingCP, 0)) == 1
                CPDuration = getRow(CPDurationsheet, i)
                CPAttrition = getRow(CPAttritionsheet, i)
                j = 3
                push!(GuyCareerPaths, CareerPath())
                GuyCareerPaths[end].RecruitmentAge = Int(getCellValue(getCell(CPDuration, 12)))
                while true
                    StLevel = getCellValue(getCell(WorkingCP, j))
                    j+=1
                    NdLevel = getCellValue(getCell(WorkingCP, j))
                    j+=1
                    StatusDuration = Int(getCellValue(getCell(CPDuration, length(GuyCareerPaths[end].Path))))
                    attr = Float16(getCellValue(getCell(CPAttrition, length(GuyCareerPaths[end].Path))))
                    ActivateLevel(GuyCareerPaths[end], AcademicLevel)
                    ActivateLevel(GuyCareerPaths[end], Assignment)
                    addStateToCareerPath(GuyCareerPaths[end], CareerStatus(academicLevel=AcademicLevel(StLevel), affiliation=Affiliation(SubPopulations[aff].Name), assignment=Assignment(NdLevel), min=StatusDuration, max=StatusDuration, attrition=attr))
                    Cellj = getCell(WorkingCP, j)
                    trans = getCellValue(Cellj)
                    j+=1
                    # if the end of the career path reached
                    if trans == "PE"
                        # if it finnishes with a pension we set the career path length to MaxCareerPathLength
                        break
                    elseif trans == "B-"
                        # if it finnishes with a B- we set the career path length to 7
                        break
                    elseif trans == nothing
                        # if it has an undefined end we set the career path length to 0
                        break
                    end
                end
            end
        end
    end
    BasicCPNb = Int(length(GuyCareerPaths)/length(SubPopulations))
println("End CP")
    # Reading Initial manpower _____________________________________________________

    InitMPnumber = jcall(InitMPsheet, "getLastRowNum", jint, ())

    if InitMPnumber < 1
        warn("Init were not defined.")
    else
        for i in 1:InitMPnumber
            MPRow = getRow(InitMPsheet, i)
            push!(InitManpower, InitMPcluster(getCellValue(getCell(MPRow, 1)), getCellValue(getCell(MPRow, 2)), Int(SimulationYear - getCellValue(getCell(MPRow, 0))), Int(getCellValue(getCell(MPRow, 3)))))
        end
    end

    # Reading Recruitment Max/min __________________________________________________

    NBYears = 0

    MaxRecruitment = Vector{Int}()
    MinRecruitment = Vector{Int}()
    recruitMaxRow = getRow(RecruitmentSheet, 1)
    recruitMinRow = getRow(RecruitmentSheet, 2)

    while true
        try
            CellMax = getCellValue(getCell(recruitMaxRow, NBYears+1))
            CellMin = getCellValue(getCell(recruitMinRow, NBYears+1))
            NBYears += 1
            push!(MaxRecruitment, CellMax)
            push!(MinRecruitment, CellMin)
        catch
            break
        end
    end

    AllowableDeviation = getCellValue(getCell(getRow(RecruitmentSheet, 4), 1))

    # Reading Objective manpower ___________________________________________________

    Objnumber = jcall(Objsheet, "getLastRowNum", jint, ())
    #Cellj = 0
    for i in 1:Objnumber
        ObjRow = getRow(Objsheet, i)
        PriorI = getCellValue(getCell(ObjRow, 0))
        NbDemnded = Int(getCellValue(getCell(ObjRow, 1)))
        tolInit = getCellValue(getCell(ObjRow, 2))
        TolEnd =  getCellValue(getCell(ObjRow, 3))
        plotnb = getCellValue(getCell(ObjRow, 4))
        Alfa = (tolInit-TolEnd)/(NBYears-1)
        push!(MPObjectives, ManpowerObjective(Int(PriorI), nb = NbDemnded, initT = tolInit, endT = TolEnd, alfa = Alfa, plot = plotnb))
        indexObj = length(MPObjectives)
        j = 5
        while true
            try
                Cellj = getCellValue(getCell(ObjRow, j))
                j += 1
            catch
                if length(MPObjectives[indexObj].Objectives) == 0
                    push!(MPObjectives[indexObj].Objectives, "")
                    push!(MPObjectives[indexObj].targetLevel, AcademicLevel)
                end
                break
            end
            if Cellj == nothing
                #push!(MPObjectives[indexObj].Objectives, Requirements)
                if length(MPObjectives[indexObj].Objectives) == 0
                    push!(MPObjectives[indexObj].Objectives, "")
                    push!(MPObjectives[indexObj].targetLevel, AcademicLevel)
                end
                break
            end

            if !(TypesDict[Cellj] in MPObjectives[indexObj].targetLevel)
                push!(MPObjectives[indexObj].targetLevel, TypesDict[Cellj])
            end
            Cellj = getCellValue(getCell(ObjRow, j))
            j += 1
            push!(MPObjectives[indexObj].Objectives, Cellj)
        end
    end
#end

# Reading Dependencies  ___________________________________________________
DepNB = Int(getCellValue(getCell(getRow(DependenciesSheet, 0), 1)))
NBConstDepend = 0
Dependencies = Vector{Dependency}()
#println("Dependencies", Dependencies)
RowIndex = 2
tempoSet = Vector{Int}()

for i in 1:DepNB
    nbSets = Int(getCellValue(getCell(getRow(DependenciesSheet, RowIndex), 1)))
    NBConstDepend += nbSets
    TempoDep = Dependency()
    RowIndex += 1
    for j in 1:nbSets
        SetRow = getRow(DependenciesSheet, RowIndex)
        RowIndex += 1
        cellindex = 1
        while true
            try
                CPinx = Int(getCellValue(getCell(SetRow, cellindex))); cellindex += 1
                SupPopName = getCellValue(getCell(SetRow, cellindex)); cellindex += 1
                SubPopIndex = findfirst(x -> x.Name == SupPopName, SubPopulations)
                push!(tempoSet, CPinx + ((SubPopIndex-1)*BasicCPNb))
            catch
                push!(TempoDep.Sets, tempoSet)
                tempoSet = Vector{Int}()
                break
            end
        end

    end

    SetRow = getRow(DependenciesSheet, RowIndex)
    RowIndex += 1
    for i in 1:nbSets
        push!(TempoDep.Percentages, getCellValue(getCell(SetRow, i)))
    end

    TempoDep.Priority = Int(getCellValue(getCell(getRow(DependenciesSheet, RowIndex), 1)))
    RowIndex += 1

    push!(Dependencies, TempoDep)

    RowIndex += 1
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
using Distributed
#addprocs(3; exeflags="--project")
@everywhere using FermiCG
@everywhere using PyCall
@everywhere using Printf

@everywhere pydir = joinpath(dirname(pathof(FermiCG)), "python")
@everywhere pushfirst!(PyVector(pyimport("sys")."path"), pydir)
@everywhere ENV["PYTHON"] = Sys.which("python")

#for i in 1:10
@sync @distributed for i in 1:10
    pytest = pyimport("test_pytest")
    a = pytest.test_print(i)
    @printf("%4i %4i\n",i,a)
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
"""
    probplot_data(fm::fittedModel)

Return the probability plot data in a DataFrame.
"""
function probplot_data end

"""
    qqplot_data(fm::fittedModel)

Return the quantile-quantile plot data in a DataFrame.
"""
function qqplot_data end

"""
    returnlevelplot_data(fm::fittedModel)

Return the return level plot data in a DataFrame.
"""
function returnlevelplot_data end

"""
    histplot_data(fm::fittedModel)

Return the histogram plot data in a Dictionary.
"""
function hisplot_data end

function probplot_data(fm::BayesianEVA)::DataFrame

    checkstationarity(fm.model)

    y, p̂ = ecdf(fm.model.data.value)

    dist = getdistribution(fm)

    p_sim = Array{Float64}(undef, length(dist), length(y))

    i = 1

    for d in eachslice(dist, dims=1)
        p_sim[i,:] = cdf.(d, y)
        i +=1
    end

    p = vec(mean(p_sim, dims=1))

    return DataFrame(Model = p, Empirical = p̂)

end

function probplot_data(fm::MaximumLikelihoodEVA)::DataFrame

    checkstationarity(fm.model)

    y, p̂ = ecdf(fm.model.data.value)

    dist = getdistribution(fm.model, fm.θ̂)[]

    p = cdf.(dist, y)

    return DataFrame(Model = p, Empirical = p̂)

end

function probplot_data(fm::pwmEVA)::DataFrame

    checkstationarity(fm.model)

    y, p̂ = ecdf(fm.model.data.value)

    dist = getdistribution(fm.model, fm.θ̂)[]

    p = cdf.(dist, y)

    return DataFrame(Model = p, Empirical = p̂)

end

function qqplot_data(fm::BayesianEVA)::DataFrame

    checkstationarity(fm.model)

    y, p = ecdf(fm.model.data.value)

    dist = getdistribution(fm)

    q_sim = Array{Float64}(undef, length(dist), length(y))

    i = 1

    for d in eachslice(dist, dims=1)
        q_sim[i,:] = quantile.(d, p)
        i +=1
    end

    q = vec(mean(q_sim, dims=1))

    return DataFrame(Model = q, Empirical = y)

end

function qqplot_data(fm::MaximumLikelihoodEVA)::DataFrame

    checkstationarity(fm.model)

    y, p = ecdf(fm.model.data.value)

    n = length(y)

    q = Vector{Float64}(undef, n)

    for i in eachindex(p)
       q[i] = quantile(fm, p[i])[1]
    end

    return DataFrame(Model = q, Empirical = y)

end

function qqplot_data(fm::pwmEVA)::DataFrame

    checkstationarity(fm.model)

    y, p = ecdf(fm.model.data.value)

    n = length(y)

    q = Vector{Float64}(undef, n)

    for i in eachindex(p)
       q[i] = quantile(fm, p[i])[1]
    end

    return DataFrame(Model = q, Empirical = y)

end

function returnlevelplot_data(fm::BayesianEVA)::DataFrame

    checkstationarity(fm.model)

    y, p = ecdf(fm.model.data.value)

    T = 1 ./ (1 .- p)

    dist = getdistribution(fm)

    q_sim = Array{Float64}(undef, length(dist), length(y))

    i = 1
    for d in eachslice(dist, dims=1)
        q_sim[i,:] = quantile.(d, p)
        i +=1
    end

    q = vec(mean(q_sim, dims=1))

    return DataFrame(Data = y, Period = T, Level = q)

end

function returnlevelplot_data(fm::MaximumLikelihoodEVA)::DataFrame

    checkstationarity(fm.model)

    df = qqplot_data(fm)

    y, p = ecdf(fm.model.data.value)

    T = 1 ./ (1 .- p)

    n = length(y)

    q = Vector{Float64}(undef, n)

    for i in eachindex(p)
       q[i] = quantile(fm, p[i])[1]
    end

    return DataFrame(Data = y, Period = T, Level = q)

end

function returnlevelplot_data(fm::pwmEVA)::DataFrame

    checkstationarity(fm.model)

    y, p = ecdf(fm.model.data.value)

    T = 1 ./ (1 .- p)

    n = length(y)

    q = Vector{Float64}(undef, n)

    for i in eachindex(p)
       q[i] = quantile(fm, p[i])[1]
    end

    return DataFrame(Data = y, Period = T, Level = q)

end

function histplot_data(fm::BayesianEVA)::Dict

    checkstationarity(fm.model)

    x = fm.model.data.value
    n = length(x)
    nbin = Int64(ceil(sqrt(n)))

    dist = getdistribution(fm)

    xmin = quantile(dist[1], 1/1000)
    xmax = quantile(dist[1], 1 - 1/1000)
    xp = range(xmin, xmax, length=1000)

    h_sim = Array{Float64}(undef, length(dist), length(xp))

    i = 1
    for d in eachslice(dist, dims=1)
        h_sim[i,:] = pdf.(d, xp)
        i +=1
    end

    h = vec(mean(h_sim, dims=1))

    return Dict(:h => DataFrame(Data = x), :d => DataFrame(DataRange = xp, Density = h),
        :nbin => nbin, :xmin => xmin, :xmax => xmax)

end

function histplot_data(fm::MaximumLikelihoodEVA)::Dict

    checkstationarity(fm.model)

    dist = getdistribution(fm)[1]

    x = fm.model.data.value
    n = length(x)
    nbin = Int64(ceil(sqrt(n)))

    xmin = quantile(dist, 1/1000)
    xmax = quantile(dist, 1 - 1/1000)
    xp = range(xmin, xmax, length=1000)

    return Dict(:h => DataFrame(Data = x), :d => DataFrame(DataRange = xp, Density = pdf.(dist, xp)),
        :nbin => nbin, :xmin => xmin, :xmax => xmax)

end

function histplot_data(fm::pwmEVA)::Dict

    checkstationarity(fm.model)

    dist = getdistribution(fm)[1]

    x = fm.model.data.value
    n = length(x)
    nbin = Int64(ceil(sqrt(n)))

    xmin = quantile(dist, 1/1000)
    xmax = quantile(dist, 1 - 1/1000)
    xp = range(xmin, xmax, length=1000)

    return Dict(:h => DataFrame(Data = x), :d => DataFrame(DataRange = xp, Density = pdf.(dist, xp)),
        :nbin => nbin, :xmin => xmin, :xmax => xmax)

end

"""
    mrlplot_data(y::Vector{<:Real}, steps::Int = 100)::DataFrame

Compute the mean residual life from vector `y`.

The set of thresholds ranges from `minimum(y)` to the second-to-last larger
value in `steps` number of steps.
"""
function mrlplot_data(y::Vector{<:Real}, steps::Int = 100)::DataFrame

    @assert steps > 0 "the number of steps should be positive"

    umin = minimum(y)
    umax = sort(y)[end-2]

    threshold = range(umin, stop = umax, length = steps)

    df = DataFrame(
        Threshold = Float64[],
        mrl = Float64[],
        lbound = Float64[],
        ubound = Float64[],
    )

    for u in threshold
        z = y[y.>u] .- u

        n = length(z)
        m = mean(z)
        s = std(z)

        push!(df, [u; m; m - 1.96 * s / sqrt(n); m + 1.96 * s / sqrt(n)])
    end

    return df

end


text:
# BLAST+ Wrapper
# ==============
#
# Wrapper for BLAST+ command line functions.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

struct BLASTResult
    bitscore::Float64
    expect::Float64
    queryname::String
    hitname::String
    hit::BioSequence
    alignment::AlignedSequence
end

"""
    readblastXML(blastrun::AbstractString)
Parse XML output of a blast run. Input is an XML string eg:
```julia
results = readstring(open("blast_results.xml"))
readblastXML(results)
```
Returns Vector{BLASTResult} with the sequence of the hit, the Alignment with query sequence, bitscore and expect value
"""
function readblastXML(blastrun::AbstractString; seqtype="nucl")
    dc = EzXML.parsexml(blastrun)
    rt = EzXML.root(dc)
    results = BLASTResult[]
    for iteration in findall("/BlastOutput/BlastOutput_iterations/Iteration", rt)
        queryname = EzXML.nodecontent(findfirst("Iteration_query-def", iteration))
        for hit in findall("Iteration_hits", iteration)
            if EzXML.countelements(hit) > 0
                hitname = EzXML.nodecontent(findfirst("./Hit/Hit_def", hit))
                hsps = findfirst("./Hit/Hit_hsps", hit)
                if seqtype == "nucl"
                    qseq = DNASequence(EzXML.nodecontent(findfirst("./Hsp/Hsp_qseq", hsps)))
                    hseq = DNASequence(EzXML.nodecontent(findfirst("./Hsp/Hsp_hseq", hsps)))
                elseif seqtype == "prot"
                    qseq = AminoAcidSequence(EzXML.nodecontent(findfirst("./Hsp/Hsp_qseq", hsps)))
                    hseq = AminoAcidSequence(EzXML.nodecontent(findfirst("./Hsp/Hsp_hseq", hsps)))
                else
                    throw(error("Please use \"nucl\" or \"prot\" for seqtype"))
                end

                aln = AlignedSequence(qseq, hseq)
                bitscore = parse(Float64, EzXML.nodecontent(findfirst("./Hsp/Hsp_bit-score", hsps)))
                expect = parse(Float64, EzXML.nodecontent(findfirst("./Hsp/Hsp_evalue", hsps)))
                push!(results, BLASTResult(bitscore, expect, queryname, hitname, hseq, aln))
            end
        end
    end
    return results
end

"""
`readblastXML(blastrun::Cmd)`
Parse command line blast query with XML output. Input is the blast command line command, eg:
```julia
blastresults = `blastn -query seq1.fasta -db some_database -outfmt 5`
readblastXML(blastresults)
```
Returns Vector{BLASTResult} with the sequence of the hit, the Alignment with query sequence, bitscore and expect value
"""
function readblastXML(blastrun::Cmd; seqtype="nucl")
    return readblastXML(read(blastrun, String), seqtype=seqtype)
end

"""
`blastn(query, subject, flags...)``
Runs blastn on `query` against `subject`.
    Subjects and queries may be file names (as strings), DNASequence type or
    Array of DNASequence.
    May include optional `flag`s such as `["-perc_identity", 95,]`. Do not use `-outfmt`.
"""
function blastn(query::AbstractString, subject::AbstractString, flags=[]; db::Bool=false)
    if db
        results = readblastXML(`blastn -query $query -db $subject $flags -outfmt 5`)
    else
        results = readblastXML(`blastn -query $query -subject $subject $flags -outfmt 5`)
    end
    return results
end

function blastn(query::DNASequence, subject::DNASequence, flags=[])
    querypath, subjectpath = makefasta(query), makefasta(subject)
    return blastn(querypath, subjectpath, flags)
end

function blastn(query::DNASequence, subject::Vector{DNASequence}, flags=[])
    querypath, subjectpath = makefasta(query), makefasta(subject)
    blastn(querypath, subjectpath, flags)
end

function blastn(query::DNASequence, subject::AbstractString, flags=[]; db::Bool=false)
    querypath = makefasta(query)
    if db
        return blastn(querypath, subject, flags, db=true)
    else
        return blastn(querypath, subject, flags)
    end
end

function blastn(query::Vector{DNASequence}, subject::Vector{DNASequence}, flags=[])
    querypath, subjectpath = makefasta(query), makefasta(subject)
    return blastn(querypath, subjectpath, flags)
end

"""
`blastp(query, subject, flags...)``
Runs blastn on `query` against `subject`.
    Subjects and queries may be file names (as strings), `BioSequence{AminoAcidSequence}` type or
    Array of `BioSequence{AminoAcidSequence}`.
    May include optional `flag`s such as `["-perc_identity", 95,]`. Do not use `-outfmt`.
"""
function blastp(query::AbstractString, subject::AbstractString, flags=[]; db::Bool=false)
    if db
        results = readblastXML(`blastp -query $query -db $subject $flags -outfmt 5`, seqtype = "prot")
    else
        results = readblastXML(`blastp -query $query -subject $subject $flags -outfmt 5`, seqtype = "prot")
    end
    return results
end

function blastp(query::AminoAcidSequence, subject::AminoAcidSequence, flags=[])
    querypath, subjectpath = makefasta(query), makefasta(subject)
    return blastp(querypath, subjectpath, flags)
end

function blastp(query::AminoAcidSequence, subject::Vector{AminoAcidSequence}, flags=[])
    querypath, subjectpath = makefasta(query), makefasta(subject)
    return blastp(querypath, subjectpath, flags)
end

function blastp(query::AminoAcidSequence, subject::AbstractString, flags=[]; db::Bool=false)
    querypath = makefasta(query)
    if db
        return blastp(querypath, subject, flags, db=true)
    else
        return blastp(querypath, subject, flags)
    end
end

function blastp(query::Vector{AminoAcidSequence}, subject::Vector{AminoAcidSequence}, flags=[])
    querypath, subjectpath = makefasta(query), makefasta(subject)
    return blastp(querypath, subjectpath, flags)
end

# Create temporary fasta-formated file for blasting.
function makefasta(sequence::BioSequence)
    path, io = mktemp()
    write(io, ">$path\n$(convert(String, sequence))\n")
    close(io)
    return path
end

# Create temporary multi fasta-formated file for blasting.
function makefasta(sequences::Vector{T}) where T <: BioSequence
    path, io = mktemp()
    counter = 1
    for sequence in sequences
        write(io, ">$path$counter\n$(convert(String, sequence))\n")
        counter += 1
    end
    close(io)
    return path
end


text:
function ball_on_support(obj::SupportBallGenerator, arg0::List)
    return jcall(obj, "ballOnSupport", EnclosingBall, (List,), arg0)
end


text:
#=
# 6. Variable free stream
In this notebook we will simulate the flow with a time-varying free stream past a
stationary body. To demonstrate this, we will solve for oscillatory flow past a
rectangular object, in which the $x$ component of the free stream is

$$U_\infty(t) = A \sin(\Omega t)$$
=#
using ViscousFlow
#-
using Plots

#=
### Problem specification
We will set the Reynolds number to 200
=#
Re = 200 # Reynolds number

#=
In order to set the time-varying free stream, we will use the `OscillationX` function,
which generates a set of oscillatory kinematics. Type `?OscillationX` to
learn more.

Since we are using this function to specify the velocity (the first derivative of
the position), then we set
=#
Ux = 0.0
Ω  = 1.0
Ax = 1.0
fs = OscillationX(Ux,Ω,Ax,π/2)

#=
As with body motion, it is useful to verify that this provides the expected behavior
by plotting it:
=#
plot(fs)

# Now let us carry on with the other usual steps:

# ### Discretize
xlim = (-2.0,2.0)
ylim = (-1.5,1.5)

Δx, Δt = setstepsizes(Re,gridRe=3.0)

#=
### Set up bodies
Here, we will set up a rectangle in the center of the domain
=#
body = Rectangle(0.25,0.5,1.5Δx)
T = RigidTransform((0.0,0.0),0.0)
T(body)
#-
plot(body,xlim=xlim,ylim=ylim)

#=
### Construct the system structure
This step is like the previous notebook, but now we also provide the body and the freestream:
=#
sys = NavierStokes(Re,Δx,xlim,ylim,Δt,body,freestream = fs)

#=
### Initialize
Now, we initialize with zero vorticity
=#
u0 = newstate(sys)
# and create the integrator
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)

#=
### Solve
Now we are ready to solve the problem. Let's advance the solution to $t = 1$.
=#
@time step!(integrator,1.0)
sol = integrator.sol;

#=
### Examine
Let's look at the flow field at the end of this interval
=#
plot(
plot(vorticity(integrator),sys,title="Vorticity",clim=(-10,10),levels=range(-10,10,length=30), color = :RdBu,ylim=ylim),
plot(streamfunction(integrator),sys,title="Streamlines",ylim=ylim,color = :Black),
    size=(700,300)
    )

# Now let's make a movie
@gif for t in sol.t
    #plot(streamfunction(u,sys,t),sys, color = :Black)
    plot(vorticity(sol,sys,t),sys,clim=(-10,10),levels=range(-10,10,length=30), color = :RdBu)
end every 10

#=
#### Compute the force history
Just as we did for the stationary body in a constant free stream
=#
fx, fy = force(sol,sys,1);
# Plot them
plot(
plot(sol.t,2*fx,xlim=(0,Inf),ylim=(-6,6),xlabel="Convective time",ylabel="\$C_D\$",legend=:false),
plot(sol.t,2*fy,xlim=(0,Inf),ylim=(-6,6),xlabel="Convective time",ylabel="\$C_L\$",legend=:false),
    size=(800,350)
)

# The mean drag and lift coefficients are
meanCD = GridUtilities.mean(2*fx)
#-
meanCL = GridUtilities.mean(2*fy)


text:
#---------------------------------------------------------------------#
#This function computes the global Mass and Differentiation Matrices.
#Written by F.X. Giraldo on April 19 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function global_matrices(Me,De,Dwe,Fe,intma,coord,Npoin,Ne,Np,periodicity,DFloat)

    #Initialize
    M=zeros(DFloat,Npoin,Npoin)
    D=zeros(DFloat,Npoin,Npoin)
    Dw=zeros(DFloat,Npoin,Npoin)
    F=zeros(DFloat,Npoin,Npoin)
    inode=zeros(Int64,Np)
    x=zeros(DFloat,Np)

    #Loop through elements
    for e=1:Ne
        #Store Coordinates
        for i=1:Np
            I=intma[i,e]
            inode[i]=periodicity[I]
            x[i]=coord[I]
        end

        #Element Length & Jacobian
        dx=x[Np]-x[1]
        jac=dx/2

        #Perform DSS
        for j=1:Np
            J=inode[j]
            for i=1:Np
                I=inode[i]
                M[I,J]+=jac*Me[i,j]
                D[I,J]+=De[i,j]
                Dw[I,J]+=Dwe[i,j]
                F[I,J]+=Fe[i,j]
            end #j
        end #i
    end #ie

    #Periodicity
    if periodicity[Npoin] == periodicity[1]
        M[Npoin,Npoin]=1
    end

    #check Fmatrix
    for i=1:Npoin
        for j=1:Npoin
#            println("i = ",i," j = ",j," F = ",F[i,j])
        end
    end

    #Return arrays
    return(M,D,Dw,F)
end


text:
#=

Testing of syntheticimages.jl
Hard to test these functions other than visually.  This script simply
runs then all to make sure that they at least run

Set 'disp' = true to display the images
=#

disp = false

using ImagePhaseCongruency, Images, TestImages

if disp
    using PyPlot
end

println("Testing syntheticimages...")

function mypause()
    println("Hit return to continue ")
    a = readline()
    return nothing
end

disp ? PyPlot.set_cmap(PyPlot.ColorMap("gray")) : nothing

#lena = Float64.(testimage("lena_gray"))
# We need the test image to be square for the phase amplitude tests.
testimg = Float64.(Gray.(testimage("lighthouse")))[1:512, 1:512]

sze = 512
img = step2line(sze; nscales=50, ampexponent=-1, ncycles=1.5, phasecycles=0.25)
if disp
    imshow(img); title("step2line() test image")
    mypause()
end

img = circsine(sze; wavelength = 40, nscales = 50, ampexponent = -1,
               offset = 0, p = 2, trim = false)
if disp
    imshow(img); title("circsign() test image")
    mypause()
end

img = starsine(sze; ncycles=10, nscales=50, ampexponent=-1, offset=0)
if disp
    imshow(img); title("starsign() test image")
    mypause()
end

img = noiseonf(sze, 1.5)
if disp
    imshow(img); title("Noise with amplitude spectrum f^{-1.5}")
    mypause()
end

newimg = nophase(testimg)
if disp
    imshow(newimg); title("Testimg with randomized phase")
    mypause()
end

newimg = quantizephase(testimg, 4)
if disp
    imshow(newimg); title("Testimg with phase quantized to 4 levels")
    mypause()
end

(newimg1, newimg2) = swapphase(testimg, testimg')
if disp
    imshow(newimg1); title("Testimg with phase of Testimg transposed")
    mypause()
    imshow(newimg2); title("Testimg transposed with phase of Testimg")
    mypause()
end

nothing


text:
# Empirical estimation of CDF and PDF

## Empirical CDF

struct ECDF{T <: AbstractVector{<:Real}}
    sorted_values::T
end

function (ecdf::ECDF)(x::Real)
    searchsortedlast(ecdf.sorted_values, x) / length(ecdf.sorted_values)
end

function (ecdf::ECDF)(v::RealVector)
    ord = sortperm(v)
    m = length(v)
    r = similar(ecdf.sorted_values, m)
    r0 = 0
    i = 1
    n = length(ecdf.sorted_values)
    for x in ecdf.sorted_values
        while i <= m && x > v[ord[i]]
            r[ord[i]] = r0
            i += 1
        end
        r0 += 1
        if i > m
            break
        end
    end
    while i <= m
        r[ord[i]] = n
        i += 1
    end
    return r / n
end

"""
    ecdf(X)

Return an empirical cumulative distribution function (ECDF) based on a vector of samples
given in `X`.

Note: this function that returns a callable composite type, which can then be applied to
evaluate CDF values on other samples.

`extrema`, `minimum`, and `maximum` are supported to for obtaining the range over which
function is inside the interval ``(0,1)``; the function is defined for the whole real line.
"""
ecdf(X::RealVector{T}) where T<:Real = ECDF(sort(X))

minimum(ecdf::ECDF) = first(ecdf.sorted_values)

maximum(ecdf::ECDF) = last(ecdf.sorted_values)

extrema(ecdf::ECDF) = (minimum(ecdf), maximum(ecdf))


text:
# function create_tables(db)
#     create_snippets = """
#     CREATE TABLE IF NOT EXISTS snippets
#     (
#         code TEXT,
#         user_id INTEGER,
#         snippet TEXT,
#         tags TEXT,
#         created TEXT
#     )
#     """
#     DBInterface.execute(db, create_snippets)
#     SQLite.createindex!(db, "snippets", "code_index", "code"; unique = true, ifnotexists = true)
# end

function up(conn)
    create_table(conn, TimedMessage)
    create_table(conn, Sender)
end


text:
module TensorOperations

using Compat  # for sizehint!

export LabelError
export tensorcopy, tensoradd, tensortrace, tensorcontract, tensorproduct, scalar, reset_tcbuffer

# LabelError
#------------
# Tensor operations, either using methods or via index notation,
# are specified by assigning labels to the different indices of
# an array. All errors related to invalid label configurations
type LabelError <: Exception
    msg::String
end

# Constants that define base case in recursive algorithms
#---------------------------------------------------------
# for tensorcopy, tensoradd and tensortrace
const TBASELENGTH=512
# note: total number elements involved = 2*512 = 1024

# for tensorcontract
const OBASELENGTH=16 # total size of all open dimensions in one of the two contraction partners
const CBASELENGTH=24 # total size of all contraction dimensions
# note: total number elements involved = 16*24*2+16*16 = 1024

# Tensor Operations
#-------------------
using Cartesian

include("kernels.jl")
include("tensorcopy.jl")
include("tensoradd.jl")
include("tensortrace.jl")
include("tensorcontract.jl")
include("tensorproduct.jl")

# Scalar
#--------
scalar(C::StridedArray)=length(C)==1 ? C[1] : throw(DimensionMismatch())

end # module

# Index Notation
#----------------
include("indexnotation.jl")


text:
using OrdinaryDiffEq,DiffEqProblemLibrary, DiffEqDevTools, Test

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

prob = prob_ode_2Dlinear
sol =solve(prob,Rosenbrock32(),dt=1/2^4)

sol =solve(prob,ExplicitRK(tableau = constructBogakiShampine3()),dt=1/2^4)
val1 = maximum(abs.(sol.u[end] - sol.u_analytic[end]))

sol2 =solve(prob,ExplicitRK(tableau = constructDormandPrince()),dt=1/2^4)
val2 = maximum(abs.(sol2.u[end] - sol2.u_analytic[end]))

sol3 =solve(prob,ExplicitRK(tableau = constructRKF8(Float64)),dt=1/2^4)
val3 = maximum(abs.(sol3.u[end] - sol3.u_analytic[end]))

@test length(sol.t)>length(sol2.t)>=length(sol3.t)
@test max(val1,val2,val3)<2e-3

function lorenz(u,p,t)
  [10.0(u[2]-u[1])
  u[1]*(28.0-u[3]) - u[2]
  u[1]*u[2] - (8/3)*u[3]]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem{false}(lorenz,u0,tspan)
sol = solve(prob,QNDF())
@test length(sol.t) < 5000

function lorenz(du,u,p,t)
  du[1] = 10.0(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem{true}(lorenz,u0,tspan)
sol = solve(prob,QNDF())
@test length(sol.t) < 5000


text:
"""
    MOBIL
See Treiber & Kesting, 'Modeling Lane-Changing Decisions with MOBIL'

# Constructor

`MOBIL(;kwargs...)`
The keyword arguments are the fields highlighted below.

# Fields
- `dir::Int = DIR_MIDDLE`
- `mlon::LaneFollowingDriver=IntelligentDriverModel()`
- `safe_decel::Float64=2.0`
- `politeness::Float64=0.35`
- `advantage_threshold::Float64=0.1`
"""
@with_kw mutable struct MOBIL <: LaneChangeModel{LaneChangeChoice}
    dir::Int64 = DIR_MIDDLE
    mlon::LaneFollowingDriver = IntelligentDriverModel()
    safe_decel::Float64 = 2.0 # safe deceleration (positive value)
    politeness::Float64 = 0.35 # politeness factor (suggested p ∈ [0.2,0.5])
    advantage_threshold::Float64 = 0.1 # Δaₜₕ
end

"""
Set the desired speed of the longitudinal model within MOBIL
"""
function set_desired_speed!(model::MOBIL, v_des::Float64)
    set_desired_speed!(model.mlon, v_des)
    model
end
function observe!(model::MOBIL, scene::Scene{Entity{S, D, I}}, roadway::Roadway, egoid::I) where {S, D, I}

    vehicle_index = findfirst(egoid, scene)
    veh_ego = scene[vehicle_index]
    v = vel(veh_ego.state)
    egostate_M = veh_ego.state

    ego_lane = get_lane(roadway, veh_ego)

    fore_M = find_neighbor(scene, roadway, veh_ego, targetpoint_ego=VehicleTargetPointFront(), targetpoint_neighbor=VehicleTargetPointRear())
    rear_M = find_neighbor(scene, roadway, veh_ego, rear=true, targetpoint_ego=VehicleTargetPointRear(), targetpoint_neighbor=VehicleTargetPointFront())

    # accel if we do not make a lane change
    accel_M_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid)).a
    model.dir = DIR_MIDDLE

    advantage_threshold = model.advantage_threshold

    if n_lanes_left(roadway, ego_lane) > 0

        rear_L = find_neighbor(scene, roadway, veh_ego,
                              lane=leftlane(roadway, veh_ego),
                              rear=true,
                              targetpoint_ego=VehicleTargetPointRear(),
                              targetpoint_neighbor=VehicleTargetPointFront())

        # candidate position after lane change is over
        footpoint = get_footpoint(veh_ego)
        lane = get_lane(roadway, veh_ego)
        lane_L = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        roadproj = proj(footpoint, lane_L, roadway)
        frenet_L = Frenet(RoadIndex(roadproj), roadway)
        egostate_L = VehicleState(frenet_L, roadway, vel(veh_ego.state))

        Δaccel_n = 0.0
        passes_safety_criterion = true
        if rear_L.ind != nothing
            id = scene[rear_L.ind].id
            accel_n_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

            # update ego state in scene
            scene[vehicle_index] = Entity(veh_ego, egostate_L)
            accel_n_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

            body = inertial2body(get_rear(scene[vehicle_index]), get_front(scene[rear_L.ind])) # project ego to be relative to target
            s_gap = body.x

            # restore ego state
            scene[vehicle_index] = veh_ego
            passes_safety_criterion = accel_n_test ≥ -model.safe_decel && s_gap ≥ 0
            Δaccel_n = accel_n_test - accel_n_orig
        end

        if passes_safety_criterion

            Δaccel_o = 0.0
            if rear_M.ind != nothing
                id = scene[rear_M.ind].id
                accel_o_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

                # update ego state in scene
                scene[vehicle_index] = Entity(veh_ego, egostate_L)
                accel_o_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

                # restore ego state
                scene[vehicle_index] = veh_ego
                Δaccel_o = accel_o_test - accel_o_orig
            end

            # update ego state in scene
            scene[vehicle_index] = Entity(veh_ego, egostate_L)
            accel_M_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid)).a
            # restore ego state
            scene[vehicle_index] = veh_ego

            Δaccel_M = accel_M_test - accel_M_orig

            Δaₜₕ = Δaccel_M + model.politeness*(Δaccel_n + Δaccel_o)

            if Δaₜₕ > advantage_threshold
                model.dir = DIR_LEFT
                advantage_threshold = Δaₜₕ
            end
        end
    end

    if n_lanes_right(roadway, ego_lane) > 0

        rear_R = find_neighbor(scene, roadway, veh_ego, lane=rightlane(roadway, veh_ego), targetpoint_ego=VehicleTargetPointRear(), targetpoint_neighbor=VehicleTargetPointFront())
        # candidate position after lane change is over
        footpoint = get_footpoint(veh_ego)
        lane = roadway[veh_ego.state.posF.roadind.tag]
        lane_R = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
        roadproj = proj(footpoint, lane_R, roadway)
        frenet_R = Frenet(RoadIndex(roadproj), roadway)
        egostate_R = VehicleState(frenet_R, roadway, vel(veh_ego.state))

        Δaccel_n = 0.0
        passes_safety_criterion = true
        if rear_R.ind != nothing
            id = scene[rear_R.ind].id
            accel_n_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

            # update ego vehicle in scene
            scene[vehicle_index] = Entity(veh_ego, egostate_R)
            accel_n_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

            body = inertial2body(get_rear(scene[vehicle_index]), get_front(scene[rear_R.ind])) # project ego to be relative to target
            s_gap = body.x

            # restore ego vehicle state
            scene[vehicle_index] = veh_ego

            passes_safety_criterion = accel_n_test ≥ -model.safe_decel && s_gap ≥ 0
            Δaccel_n = accel_n_test - accel_n_orig
        end

        if passes_safety_criterion

            Δaccel_o = 0.0
            if rear_M.ind != nothing
                id = scene[rear_M.ind].id
                accel_o_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a

                # update ego vehicle in scene
                scene[vehicle_index] = Entity(veh_ego, egostate_R)
                accel_o_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id)).a
                # restore ego vehicle state
                scene[vehicle_index] = veh_ego

                Δaccel_o = accel_o_test - accel_o_orig
            end

            # update ego vehicle in scene
            scene[vehicle_index] = Entity(veh_ego, egostate_R)
            accel_M_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid)).a
            # restor ego vehicle state
            scene[vehicle_index] = veh_ego

            Δaccel_M = accel_M_test - accel_M_orig

            Δaₜₕ = Δaccel_M + model.politeness*(Δaccel_n + Δaccel_o)
            if Δaₜₕ > advantage_threshold
                model.dir = DIR_RIGHT
                advantage_threshold = Δaₜₕ
            elseif Δaₜₕ == advantage_threshold
                # in case of tie, if we are accelerating we go left, else, right
                if Δaccel_M > model.politeness*(Δaccel_n + Δaccel_o)
                    model.dir = DIR_LEFT
                else
                    model.dir = DIR_RIGHT
                end
            end
        end
    end

    model
end
Base.rand(rng::AbstractRNG, model::MOBIL) = LaneChangeChoice(model.dir)


text:
function fqPredictionPar(x::Array{Float64}, T::Int64, gibbsIter::Int64, burnIn::Int64, gEps1::Float64, gEps2::Float64, thLow::Float64,
  thUp::Float64, zLow::Float64, zUp::Float64, thin::Int64, samplerSeed::Int64)

	srand(samplerSeed)

  filename = "/Results/Par"

  savelocation = string(pwd(), filename, "/seed$samplerSeed")
  mkpath(savelocation)

	# initialization
	n = length(x)
  nn = Int((gibbsIter - burnIn) / thin)

  if T .> 0 # forward-prediction
    x = append!(x, zeros(T))
    for i in 1:T
      x[n+i] = rand(Uniform(-1.,1.))
    end
    ss = n + T
    xpred = zeros(nn, T)
  else # simple GSBR
    ss = n
  end

  theta = zeros(6) # polynomial coefficients vector
  sampledTheta = zeros(nn, 6) # matrix to store sample thetas

  x0s = zeros(nn, 2)
  sx01, sx02 = 0.5, 0.5

  lambdas = zeros(nn)
  slam = 0.5

  toler = 1e-06

	display("Starting MCMC ...")

	##########################################################################
	# MCMC
	##########################################################################

	for iter in 1:gibbsIter

      ## Sample precision
      ##########################################################################

      temp = (x[1] - henon(theta, sx01, sx02)) ^ 2 + (x[2] - henon(theta, x[1], sx01)) ^ 2
      for i in 3:ss
        temp += (x[i] - henon(theta, x[i-1], x[i-2])) ^ 2
      end
      shapel = gEps1 + 0.5 * ss
      ratel = gEps2 + 0.5 * temp
      # println("shape: $shapel")
      # println("rate: $ratel")
      slam = rand(Gamma(shapel, 1. / ratel))

      ## Sample θ₀ = theta[1]
      ##########################################################################

      muth = (x[1] - (theta[5] * sx01^2 + theta[4] * sx01 * sx02 + theta[6] * sx02^2 + theta[2] * sx01 + theta[3] * sx02)) +
                      (x[2] - (theta[5] * x[1]^2 + theta[4] * x[1] * sx01 + theta[6] * sx01^2 + theta[2] * x[1] + theta[3] * sx01))

      tauth = 1. + 1.

      for j = 3:ss
        muth +=  x[j] - (theta[5] * x[j-1]^2 + theta[4] * x[j-1] * x[j-2] + theta[6] * x[j-2]^2 + theta[2] * x[j-1] + theta[3] * x[j-2])
        tauth += 1.
      end
      muth = muth / tauth
      tauth = tauth * slam

      temp = -2. / tauth * log.(rand()) + (theta[1] - muth) ^ 2

      theta[1] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₁ = theta[2]
      # ##########################################################################
      #
      muth = (x[1] * sx01 - (theta[5] * sx01^3 + theta[4] * sx01^2 * sx02 + theta[6] * sx01 * sx02^2 +
                              theta[3] * sx01 * sx02 + theta[1] * sx01)) +
             (x[2] * x[1] - (theta[5] * x[1]^3 + theta[4] * x[1]^2 * sx01 + theta[6] * x[1] * sx01^2 +
                              theta[3] * x[1] * sx01 + theta[1] * x[1]))

      tauth = sx01 ^ 2 + x[1] ^ 2

      for j = 3:ss
        muth += (x[j] * x[j-1] - (theta[5] * x[j-1]^3 + theta[4] * x[j-1]^2 * x[j-2] + theta[6] * x[j-1] * x[j-2]^2 +
                                theta[3] * x[j-1] * x[j-2] + theta[1] * x[j-1]))
        tauth += x[j - 1] ^ 2
      end
      muth = muth / tauth
      tauth = tauth * slam

      temp = -2. / tauth * log.(rand()) + (theta[2] - muth) ^ 2

      theta[2] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₂ = theta[3]
      # ##########################################################################
      #
      muth =  (x[1] * sx02 - (theta[5] * sx01^2 * sx02 + theta[4] * sx01 * sx02^2 + theta[6] * sx02^3 +
                              theta[2] * sx01 * sx02 + theta[1] * sx02)) +
              (x[2] * sx01 - (theta[5] * x[1]^2 * sx01 + theta[4] * x[1] * sx01^2 + theta[6] * sx01^3 +
                              theta[2] * x[1] * sx01 + theta[1] * sx01))

      tauth = sx02 ^ 2 + sx01 ^ 2

      for j = 3:ss
        muth += (x[j] * x[j-2] - (theta[5] * x[j-1]^2 * x[j-2] + theta[4] * x[j-1] * x[j-2]^2 + theta[6] * x[j-2]^3 +
                                theta[2] * x[j-1] * x[j-2] + theta[1] * x[j-2]))
        tauth += x[j - 2] ^ 2
      end

      muth = muth / tauth
      tauth = tauth * slam

      temp = -2. / tauth * log.(rand()) + (theta[3] - muth) ^ 2

      theta[3] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₃ = theta[4]
      # ##########################################################################
      #
      muth = (x[1] * sx01 * sx02 - (theta[5] * sx01^3 * sx02 + theta[6] * sx01 * sx02^3 + theta[2] * sx01^2 * sx02 +
                              theta[3] * sx01 * sx02^2 + theta[1] * sx01 * sx02)) +
             (x[2] * x[1] * sx01 - (theta[5] * x[1]^3 * sx01 + theta[6] * x[1] * sx01^3 + theta[2] * x[1]^2 * sx01 +
                              theta[3] * x[1] * sx01^2 + theta[1] * x[1] * sx01))

      tauth =  sx01 ^ 2 * sx02 ^ 2 + x[1] ^ 2 *sx01 ^ 2

      for j = 3:ss
        muth += (x[j] * x[j-1] * x[j-2] - (theta[5] * x[j-1]^3 * x[j-2] + theta[6] * x[j-1] * x[j-2]^3 + theta[2] * x[j-1]^2 * x[j-2] +
                                theta[3] * x[j-1] * x[j-2]^2 + theta[1] * x[j-1] * x[j-2]))
        tauth += x[j - 1] ^ 2 * x[j - 2] ^ 2
      end

      muth = muth / tauth
      tauth = tauth * slam

      temp = -2. / tauth * log.(rand()) + (theta[4] - muth) ^ 2

      theta[4] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₄ = theta[5]
      # ##########################################################################
      #
      muth = (x[1] * sx01^2 - (theta[4] * sx01^3 * sx02 + theta[6] * sx01^2 * sx02^2 + theta[2] * sx01^3 +
                              theta[3] * sx01^2 * sx02 + theta[1] * sx01^2)) +
             (x[2] * x[1]^2 - (theta[4] * x[1]^3 * sx01 + theta[6] * x[1]^2 * sx01^2 + theta[2] * x[1]^3 +
                              theta[3] * x[1]^2 * sx01 + theta[1] * x[1]^2))

      tauth = sx01 ^ 4 + x[1] ^ 4

      for j = 3:ss
        muth += (x[j] * x[j-1]^2 - (theta[4] * x[j-1]^3 * x[j-2] + theta[6] * x[j-1]^2 * x[j-2]^2 + theta[2] * x[j-1]^3 +
                                theta[3] * x[j-1]^2 * x[j-2] + theta[1] * x[j-1]^2))
        tauth += x[j - 1] ^ 4
      end

      muth = muth / tauth
      tauth = tauth * slam

      temp = -2. / tauth * log.(rand()) + (theta[5] - muth) ^ 2

      theta[5] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₅ = theta[6]
      # ##########################################################################
      #
      muth = (x[1] * sx02^2 - (theta[5] * sx01^2 * sx02^2 + theta[4] * sx01 * sx02^3 + theta[2] * sx01 * sx02^2 +
                              theta[3] * sx02^3 + theta[1] * sx02^2)) +
             (x[2] * sx01^2 - (theta[5] * x[1]^2 * sx01^2 + theta[4] * x[1] * sx01^3 + theta[2] * x[1] * sx01^2 +
                              theta[3] * sx01^3 + theta[1] * sx01^2))

      tauth = sx02 ^ 4 + sx01 ^ 4

      for j = 3:ss
        muth += (x[j] * x[j-2]^2 - (theta[5] * x[j-1]^2 * x[j-2]^2 + theta[4] * x[j-1] * x[j-2]^3 + theta[2] * x[j-1] * x[j-2]^2 +
                                theta[3] * x[j-2]^3 + theta[1] * x[j-2]^2))
        tauth += x[j - 2] ^ 4
      end

      muth = muth / tauth
      tauth = tauth * slam

      temp = -2. / tauth * log.(rand()) + (theta[6] - muth) ^ 2

      theta[6] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      ## Sample past unobserved observations
      ##########################################################################

      if T .> 0

        if T .> 1
          for i in 1:1:(T-2)

            a4 = theta[5] ^ 2 + theta[6] ^ 2

            a3 = 2. * (theta[4] * theta[5] * x[n+i-1] + theta[2] * theta[5]) + 2. * (theta[4] * theta[6] * x[n+i+1] + theta[3] * theta[6])

            a2 = (theta[4]^2 * x[n+i-1]^2 + 2. * theta[5] * theta[6] * x[n+i-1]^2 + 2. * theta[2] * theta[4] * x[n+i-1] + 2. * theta[3] * theta[5] * x[n+i-1] -
                                2. * theta[5] * x[n+i+1] + 2. * theta[1] * theta[5] + theta[2]^2) +
                  (theta[4]^2 * x[n+i+1]^2 + 2. * theta[5] * theta[6] * x[n+i+1]^2 + 2 * theta[2] * theta[6] * x[n+i+1] + 2 * theta[3] * theta[4] * x[n+i+1] -
                                2. * theta[6] * x[n+i+2] + 2. * theta[1] * theta[6] + theta[3]^2) +
                   1.

            a1 = 2. * ( (theta[4] * theta[6] * x[n+i-1]^3 + (theta[2] * theta[6] + theta[3] * theta[4]) * x[n+i-1]^2 + (-theta[4] * x[n+i+1] + theta[1] * theta[4] + theta[2] * theta[3]) * x[n+i-1] -
                                      theta[2] * x[n+i+1] + theta[1] * theta[2]) +
                       (theta[4] * theta[5] * x[n+i+1]^3 + (theta[2] * theta[4] + theta[3] * theta[5]) * x[n+i+1]^2 -
                                    theta[4] * x[n+i+1] * x[n+i+2] + (theta[1] * theta[4] + theta[2] * theta[3]) * x[n+i+1] - theta[3] * x[n+i+2] + theta[1] * theta[3]) -
                      (theta[4] * x[n+i-2] * x[n+i-1] + theta[5] * x[n+i-1]^2 + theta[6] * x[n+i-2]^2 + theta[2] * x[n+i-1] + theta[3] * x[n+i-2] + theta[1])  )

            aux = -2. / slam * log.(rand()) + a1 * x[n+i] + a2 * x[n+i] ^ 2 + a3 * x[n+i] ^ 3 + a4 * x[n+i] ^ 4

            poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
            allRoots = Polynomials.roots(Poly(poly))
            sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
            realRoots = sort(real(allRoots[sel]))
            intervals = rangeIntersection(realRoots, [zLow; zUp])

            x[n+i] = unifmixrnd(intervals)

          end

          # 2nd last future unobserved observation
            a4 = theta[5] ^ 2

            a3 = 2. * (theta[4] * theta[5] * x[n+T-2] + theta[2] * theta[5])

            a2 = (theta[4]^2 * x[n+T-2]^2 + 2. * theta[5] * theta[6] * x[n+T-2]^2 + 2. * theta[2] * theta[4] * x[n+T-2] + 2. * theta[3] * theta[5] * x[n+T-2] -
                                2. * theta[5] * x[n+T] + 2. * theta[1] * theta[5] + theta[2]^2) +
                 1.

            a1 = 2. * ( (theta[4] * theta[6] * x[n+T-2]^3 + (theta[2] * theta[6] + theta[3] * theta[4]) * x[n+T-2]^2 + (-theta[4] * x[n+T] + theta[1] * theta[4] + theta[2] * theta[3]) * x[n+T-2] -
                                      theta[2] * x[n+T] + theta[1] * theta[2]) -
                       (theta[4] * x[n+T-3] * x[n+T-2] + theta[5] * x[n+T-2]^2 + theta[6] * x[n+T-3]^2 + theta[2] * x[n+T-2] + theta[3] * x[n+T-3] + theta[1])  )

            aux = -2. / slam * log.(rand()) + a1 * x[n+T-1] + a2 * x[n+T-1] ^ 2 + a3 * x[n+T-1] ^ 3 + a4 * x[n+T-1] ^ 4

            poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
            allRoots = Polynomials.roots(Poly(poly))
            sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
            realRoots = sort(real(allRoots[sel]))
            intervals = rangeIntersection(realRoots, [zLow; zUp])

            x[n+T-1] = unifmixrnd(intervals)

        end

        # last future unobserved observation

        a2 = slam

        a1 = - 2. * slam * (theta[4] * x[n+T-2] * x[n+T-1] + theta[5] * x[n+T-1]^2 + theta[6] * x[n+T-2]^2 + theta[2] * x[n+T-1] + theta[3] * x[n+T-2] + theta[1])

        aux = -2. * log.(rand()) + a1 * x[n+T] + a2 * x[n+T] ^ 2

        poly = [-aux / a2; a1 / a2; 1.0]
        allRoots = Polynomials.roots(Poly(poly))
        sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
        realRoots = sort(real(allRoots[sel]))
        intervals = rangeIntersection(realRoots, [zLow; zUp])

        x[n+T] = unifmixrnd(intervals)

        # a2 = 1.

        # a1 = - (theta[4] * x[n+T-2] * x[n+T-1] + theta[5] * x[n+T-1]^2 + theta[6] * x[n+T-2]^2 + theta[2] * x[n+T-1] + theta[3] * x[n+T-2] + theta[1])

        # aux = -2. / slam * log.(rand()) + a1 * x[n+T] + a2 * x[n+T] ^ 2

        # poly = [-aux / a2; a1 / a2; 1.0]
        # allRoots = Polynomials.roots(Poly(poly))
        # sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
        # realRoots = sort(real(allRoots[sel]))
        # intervals = rangeIntersection(realRoots, [zLow; zUp])

        # x[n+T] = unifmixrnd(intervals)

        # new_mu = henon(theta, x[n+T-1], x[n+T-1])
        # new_var = 1./slam
        # x[n+T] = rand(Normal(new_mu, sqrt(new_var)))

      end

      ## Sample initial conditions x₀,y₀ = x0,x00
      ##########################################################################

      # x₀

      a4 = theta[5] ^ 2 + theta[6] ^ 2

      a3 = 2. * (theta[4] * theta[5] * sx02 + theta[2] * theta[5]) + 2. * (theta[4] * theta[6] * x[1] + theta[3] * theta[6])

      a2 = (theta[4]^2 * sx02^2 + 2. * theta[5] * theta[6] * sx02^2 + 2. * theta[2] * theta[4] * sx02 + 2. * theta[3] * theta[5] * sx02 -
                            2. * theta[5] * x[1] + 2. * theta[1] * theta[5] + theta[2]^2) +
            (theta[4]^2 * x[1]^2 + 2. * theta[5] * theta[6] * x[1]^2 + 2 * theta[2] * theta[6] * x[1] + 2 * theta[3] * theta[4] * x[1] -
                            2. * theta[6] * x[2] + 2. * theta[1] * theta[6] + theta[3]^2)

      a1 = 2. * ( (theta[4] * theta[6] * sx02^3 + (theta[2] * theta[6] + theta[3] * theta[4]) * sx02^2 + (-theta[4] * x[1] + theta[1] * theta[4] + theta[2] * theta[3]) * sx02 -
                                  theta[2] * x[1] + theta[1] * theta[2]) +
                  (theta[4] * theta[5] * x[1]^3 + (theta[2] * theta[4] + theta[3] * theta[5]) * x[1]^2 -
                                theta[4] * x[1] * x[2] + (theta[1] * theta[4] + theta[2] * theta[3]) * x[1] - theta[3] * x[2] + theta[1] * theta[3]))

      aux = -2. / slam * log.(rand()) + a1 * sx01 + a2 * sx01 ^ 2 + a3 * sx01 ^ 3 + a4 * sx01 ^ 4

      poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
      allRoots = Polynomials.roots(Poly(poly))
      sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
      realRoots = sort(real(allRoots[sel]))
      intervals = rangeIntersection(realRoots, [zLow; zUp])
      sx01 = unifmixrnd(intervals)

      #
      # x00 = y₀

      a4  = theta[6] ^ 2

      a3 = 2. * theta[6] * (theta[4] * sx01 + theta[3])

      a2 = ( -2. * theta[6] * (x[1] - theta[1] - theta[2] * sx01 - theta[5] * sx01^2) + (theta[4] * sx01 + theta[3])^2 )

      a1 = -2. * (x[1] - theta[1] - theta[5] * sx01^2 - theta[2] * sx01) * (theta[4] * sx01 + theta[3])

      aux = -2. / slam * log.(rand()) + a1 * sx02 + a2 * sx02 ^ 2 + a3 * sx02 ^ 3 + a4 * sx02 ^ 4

      poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
      allRoots = Polynomials.roots(Poly(poly))
      sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
      # sel = imag(allRoots) .== 0.0 # treat as real the roots with imagine part < ϵ
      realRoots = sort(real(allRoots[sel]))
      intervals = rangeIntersection(realRoots, [zLow; zUp])
      sx02 = unifmixrnd(intervals)

      ## After Burn-In period
      ###############################wwwwwwwww###########################################

      if (iter .> burnIn) & ((iter-burnIn) % thin .== 0)

        ii = Int((iter - burnIn)/thin)

        # Store values
        sampledTheta[ii, :] = theta
        x0s[ii, :] = [sx01 sx02]
        lambdas[ii] = slam

        if T .> 0
          xpred[ii,:] = x[(n+1):(n+T)]
        end

      end

      if iter % 50000 .== 0
         println("MCMC Iterations: $iter")
      end

	end

	display("... MCMC finished !")

	## Write values in .txt files - specific path
	##########################################################################

  if T .> 0
    writedlm(string(savelocation,"/xpred.txt"), xpred)
  end
	writedlm(string(savelocation, "/thetas.txt"), sampledTheta)
	writedlm(string(savelocation,"/x0s.txt"), x0s)
  writedlm(string(savelocation,"/noise.txt"), lambdas, '\n')

 return 0

end

max_stars_repo_path: Complete quadratic/fqPredictionPar.jl

max_stars_repo_name: KosKal/GSBR-DPR

max_issues_repo_path: Complete quadratic/fqPredictionPar.jl

max_issues_repo_name: KosKal/GSBR-DPR

max_forks_repo_path: Complete quadratic/fqPredictionPar.jl

max_forks_repo_name: KosKal/GSBR-DPR


text:
using CVChannel
using Test

println("This function verifies the 3-sym CV SDP works.")
println("It does this by making sure it gets the right answer in cases that we don't need 3sym to be exact.")
println("It is in its own script to not slow down any other script, since it takes a while.")
@testset "3SymCV Works" begin
    println("Tests starting.")
    test1 = 3*wernerState(3,0)
    cv1 = pptCVPrimal(test1,3,3)
    cv1_3_sym = threeSymCVPrimal(test1,3,3)
    println("Test 1 done.")
    @test isapprox(cv1[1],cv1_3_sym[1], atol=1e-6)
    test2 = vec([1 0 0 ; 0 1 0 ; 0 0 1])*vec([1 0 0 ; 0 1 0 ; 0 0 1])'
    cv2 = pptCVPrimal(test2,3,3)
    cv2_3_sym = threeSymCVPrimal(test2,3,3)
    @test isapprox(cv2[1],cv2_3_sym[1], atol=2e-6)
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
module TerminalCalendar

using Dates
using PrettyTables

### helper functions ###

# turn month of `dt` into nx7 array of dates with missing values
function calendarize(dt::Date)
    first_day = firstdayofmonth(dt)
    last_day = lastdayofmonth(dt)

    first_missings = fill(missing, dayofweek(first_day) - 1)
    last_missings = fill(missing, 7 - dayofweek(last_day))

    dates = first_day:Day(1):last_day
    padded_dates = vcat(first_missings, dates, last_missings)

    return permutedims(reshape(padded_dates, 7, :))
end

# formatter function to replace missings, thanks to @ronisbr
ft_missing = (v, i, j) -> ismissing(v) ? "" : v

# center month-year header for each monthly calendar
function create_month_year_header(dt::Date, table_str::AbstractString)
    mo = monthname(dt)
    yr = year(dt)
    byte_index = findfirst('\n', table_str)
    table_width = length(table_str, 1, byte_index)
    month_year_width = length(mo) + length(string(yr)) + 1
    front_pad_width = Int(floor((table_width - month_year_width) / 2)) - 1
    front_pad = repeat(' ', front_pad_width)
    return string('\n', front_pad, "$mo $yr", '\n')
end

### exported functions ###

"""
    calendar(dt::Date)

Pretty prints the monthly calendar of `dt`.

# Example

```jldoctest
julia> calendar(Date(2021, 12))

              December 2021
┌─────┬─────┬─────┬─────┬─────┬─────┬─────┐
│ Mon │ Tue │ Wed │ Thu │ Fri │ Sat │ Sun │
├─────┼─────┼─────┼─────┼─────┼─────┼─────┤
│     │     │   1 │   2 │   3 │   4 │   5 │
│   6 │   7 │   8 │   9 │  10 │  11 │  12 │
│  13 │  14 │  15 │  16 │  17 │  18 │  19 │
│  20 │  21 │  22 │  23 │  24 │  25 │  26 │
│  27 │  28 │  29 │  30 │  31 │     │     │
└─────┴─────┴─────┴─────┴─────┴─────┴─────┘
```

    calendar(yr::Int, mo::Int)

Pretty prints the monthly calendar for `yr` and `mo`.

# Example

```jldoctest
julia> calendar(2021, 12)

              December 2021
┌─────┬─────┬─────┬─────┬─────┬─────┬─────┐
│ Mon │ Tue │ Wed │ Thu │ Fri │ Sat │ Sun │
├─────┼─────┼─────┼─────┼─────┼─────┼─────┤
│     │     │   1 │   2 │   3 │   4 │   5 │
│   6 │   7 │   8 │   9 │  10 │  11 │  12 │
│  13 │  14 │  15 │  16 │  17 │  18 │  19 │
│  20 │  21 │  22 │  23 │  24 │  25 │  26 │
│  27 │  28 │  29 │  30 │  31 │     │     │
└─────┴─────┴─────┴─────┴─────┴─────┴─────┘
```

And there is more! Input different things and see what works!
"""
function calendar end

# date type arguments
function calendar(dt::Date = today(); tf=tf_unicode)
    calenderized_dt = calendarize(dt)

    table_str = pretty_table(
        String,
        day.(calenderized_dt),
        header=first.(dayname.(1:7), 3),
        formatters=ft_missing,
        tf=tf,
    )

    month_year_header = create_month_year_header(dt, table_str)

    printstyled(month_year_header, bold=true)
    print(table_str)
end

function calendar(dts::Union{Vector{Date}, StepRange{Date, Month}}; tf=tf_unicode)
    for dt in dts
        calendar(dt; tf=tf)
    end
end

# integer type arguments
function calendar(year::Integer, month::Integer; tf=tf_unicode)
    dt = Date(year, month)
    calendar(dt; tf=tf)
end

function calendar(year::Integer; tf=tf_unicode)
    for month in 1:12
        calendar(Date(year, month); tf=tf)
    end
end

function calendar(
    years::Union{Vector{Int64}, UnitRange{Int64}, StepRange{Int64, Int64}};
    tf=tf_unicode
    )
    for year in years
        calendar(year; tf=tf)
    end
end

function calendar(
    years::Union{Vector{Int64}, UnitRange{Int64}, StepRange{Int64, Int64}},
    month::Integer;
    tf=tf_unicode
    )
    for year in years
        calendar(year, month; tf=tf)
    end
end

function calendar(
    year::Integer,
    months::Union{Vector{Int64}, UnitRange{Int64}, StepRange{Int64, Int64}};
    tf=tf_unicode
    )
    for month in months
        calendar(year, month; tf=tf)
    end
end

function calendar(
    years::Union{Vector{Int64}, UnitRange{Int64}, StepRange{Int64, Int64}},
    months::Union{Vector{Int64}, UnitRange{Int64}, StepRange{Int64, Int64}};
    tf=tf_unicode
    )
    for year in years, month in months
        calendar(year, month; tf=tf)
    end
end

include("precompile_statements.jl")

export calendar

# re-export table formatting from PrettyTables.jl
export tf_ascii_dots,
    tf_ascii_rounded,
    tf_borderless,
    tf_compact,
    tf_markdown,
    tf_matrix,
    tf_mysql,
    tf_simple,
    tf_unicode,
    tf_unicode_rounded

end


text:
@testset "capture" begin
    @testset "basic" begin
        @system SCapture(Controller) begin
            a => 1 ~ track
            b(a) => a + 1 ~ capture
            c(a) => a + 1 ~ accumulate
        end
        s = instance(SCapture)
        @test s.b' == 0 && s.c' == 0
        update!(s)
        @test s.b' == 2 && s.c' == 2
        update!(s)
        @test s.b' == 2 && s.c' == 4
    end

    @testset "time" begin
        @system SCaptureTime(Controller) begin
            t(x=context.clock.time) => 2x ~ track(u"hr")
            a => 1 ~ track
            b(a) => a + 1 ~ capture(time=t)
            c(a) => a + 1 ~ accumulate(time=t)
        end
        s = instance(SCaptureTime)
        @test s.b' == 0 && s.c' == 0
        update!(s)
        @test s.b' == 4 && s.c' == 4
        update!(s)
        @test s.b' == 4 && s.c' == 8
    end

    @testset "unit hour" begin
        @system SCaptureUnitHour(Controller) begin
            a => 1 ~ capture(u"hr")
        end
        s = instance(SCaptureUnitHour)
        @test iszero(s.a')
        update!(s)
        @test s.a' == 1u"hr"
        update!(s)
        @test s.a' == 1u"hr"
    end

    @testset "unit day" begin
        @system SCaptureUnitDay(Controller) begin
            a => 1 ~ capture(u"d")
        end
        s = instance(SCaptureUnitDay)
        @test iszero(s.a')
        update!(s)
        @test s.a' == 1u"hr"
        update!(s)
        @test s.a' == 1u"hr"
    end

    @testset "when" begin
        @system SCaptureWhen(Controller) begin
            t(context.clock.tick) ~ track::int
            f ~ preserve(parameter)
            w(t, f) => t <= f ~ flag
            a => 1 ~ capture
            b => 1 ~ capture(when=w)
            c => 1 ~ capture(when=!w)
        end
        n = 5
        s = instance(SCaptureWhen; config=:0 => :f => n)
        simulate!(s, stop=n)
        @test s.a' == 1
        @test s.b' == 1
        @test s.c' == 0
        simulate!(s, stop=n)
        @test s.a' == 1
        @test s.b' == 0
        @test s.c' == 1
    end
end


text:
using PooksoftOptionsKit
using PyPlot
using Dates
using Distributions

# simulation logic -
function simulate(pathToSimulationFile::String)

    # initialize -
    lattice_model = nothing
    asset_set = nothing
    asset_price_array = nothing
    data_array_at_exp = nothing

    # load the lattice model -
    result = build_simulation_lattice_data_structure(pathToSimulationFile)
    if (typeof(result.value) == PSError)
        return result.value
    else
        lattice_model = result.value
    end

    # load the asset set -
    result = build_simulation_asset_set(pathToSimulationFile)
    if (typeof(result.value) == PSError)
        return result.value
    else
        asset_set = result.value
    end

    # load asset price array -
    result = build_simulation_price_array(pathToSimulationFile)
    if (typeof(result.value) == PSError)
        return result.value
    else
        asset_price_array = result.value
    end

    # main loop -
    number_of_time_steps = length(asset_price_array)
    data_array = zeros(number_of_time_steps,2)
    for (index,asset_price) in enumerate(asset_price_array)

        # build the pricing tree -
        tree = build_ternary_price_tree(asset_set, lattice_model, asset_price)

        # compute -
        local result = option_contract_price(tree, lattice_model)

        # grab -
        data_array[index,1] = asset_price
        data_array[index,2] = result.value
    end

    # run simulation at expiration -
    # compute -
    result = compute_option_profit_and_loss_at_expiration(asset_set, asset_price_array)
    if (typeof(result.value) == PSError)
        return result.value
    else
        data_array_at_exp = result.value
    end

    # return -
    return (data_array, data_array_at_exp)
end

# setup -
path_to_design_file = "$(pwd())/tests/config/IronCondor.json"
(data_array, data_array_at_exp) = simulate(path_to_design_file)

# make a plot from the simulation -


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# Camera.jl: interface to Camera objects
export serial, model, vendor, isrunning, start!, stop!, getimage, getimage!, saveimage,
       triggermode, triggermode!,
       triggersource, triggersource!,
       trigger!,
       exposure, exposure!, exposure_limits,
       autoexposure_limits, autoexposure_limits!,
       framerate, framerate!, framerate_limits,
       gain, gain!, gain_limits,
       adcbits, adcbits!,
       gammaenable!,
       pixelformat, pixelformat!,
       acquisitionmode, acquisitionmode!,
       sensordims, imagedims, imagedims!, imagedims_limits, offsetdims, offsetdims!, offsetdims_limits,
       buffercount, buffercount!, buffermode, buffermode!, bufferunderrun, bufferfailed,
       reset!

"""
 Spinnaker SDK Camera object

 Create a camera object by referencing a CameraList, e.g.,
 `cl = CameraList()
  cam = cl[0]`

  The camera is initialised when created and deinitialised when garbage collected.
"""
mutable struct Camera
  handle::spinCamera
  names::Dict{String, String}

  function Camera(handle)
    @assert spinsys.handle != C_NULL
    @assert handle != C_NULL
    spinCameraDeInit(handle)
    spinCameraInit(handle)
    names = Dict{String, String}()
    cam = new(handle, names)
    finalizer(_release!, cam)

    # Activate chunk mode
    set!(SpinBooleanNode(cam, "ChunkModeActive"), true)
    _chunkselect(cam, ["FrameID", "FrameCounter"], "frame indentification")
    _chunkselect(cam, ["ExposureTime"], "exposure time")
    _chunkselect(cam, ["Timestamp"], "timestamps")

    # Discover ambiguous names
    cam.names["AutoExposureTimeLowerLimit"] = "AutoExposureTimeLowerLimit"
    cam.names["AutoExposureTimeUpperLimit"] = "AutoExposureTimeUpperLimit"
    cam.names["AcquisitionFrameRateEnabled"] = "AcquisitionFrameRateEnabled"

    try
      Spinnaker.get(Spinnaker.SpinFloatNode(cam, "AutoExposureTimeLowerLimit"))
    catch
      cam.names["AutoExposureTimeLowerLimit"] = "AutoExposureExposureTimeLowerLimit"
      cam.names["AutoExposureTimeUpperLimit"] = "AutoExposureExposureTimeUpperLimit"
    end
    try
      Spinnaker.get(Spinnaker.SpinBooleanNode(cam, "AcquisitionFrameRateEnabled"))
    catch
      cam.names["AcquisitionFrameRateEnabled"] = "AcquisitionFrameRateEnable"
    end

    return cam
  end
end

# Attempt to activate chunk data for each entry in chunknames
# - this allows chunk names to differ on cameras
function _chunkselect(cam::Camera, chunknames::Vector{String}, desc::String)

  fail = true
  i = 1
  while fail == true
    try
      fail = false
      set!(SpinEnumNode(cam, "ChunkSelector"), chunknames[i])
      set!(SpinBooleanNode(cam, "ChunkEnable"), true)
    catch e
      fail = true
    end
    i += 1
  end

  if fail
    @warn "Unable to enable chunk data for $(desc), tried $(chunknames), metadata may be incorrect"
  end

end

unsafe_convert(::Type{spinCamera}, cam::Camera) = cam.handle
unsafe_convert(::Type{Ptr{spinCamera}}, cam::Camera) = pointer_from_objref(cam)

function _reinit(cam::Camera)
  spinCameraDeInit(cam)
  spinCameraInit(cam)
  return cam
end

# Release handle to system
function _release!(cam::Camera)
  if cam.handle != C_NULL
    try
      stop!(cam)
    catch e
    end
    spinCameraDeInit(cam)
    spinCameraRelease(cam)
    cam.handle = C_NULL
  end
  return nothing
end

"""
  reset!(cam::Camera)

Immediately reset and reboot the camera.
"""
function reset!(cam::Camera)
  hNodeMap = Ref(spinNodeMapHandle(C_NULL))
  spinCameraGetNodeMap(cam, hNodeMap)

  hDeviceReset = Ref(spinNodeHandle(C_NULL))
  spinNodeMapGetNode(hNodeMap[], "DeviceReset", hDeviceReset);
  spinCommandExecute(hDeviceReset[])
end

# Include subfiles
include(joinpath("camera", "acquisition.jl"))
include(joinpath("camera", "analog.jl"))
include(joinpath("camera", "format.jl"))
include(joinpath("camera", "stream.jl"))
include(joinpath("camera", "digitalio.jl"))

#
# Device metadata
#

"""
  serial(::Camera) -> String

  Return camera serial number (string)
"""
serial(cam::Camera) = get(SpinStringNode(cam, "DeviceSerialNumber", CameraTLDeviceNodeMap()))

"""
  vendor(::Camera) -> String

  Return vendor name of specified camera.
"""
vendor(cam::Camera) = get(SpinStringNode(cam, "DeviceVendorName", CameraTLDeviceNodeMap()))

"""
  model(::Camera) -> String

  Return model name of specified camera.
"""
model(cam::Camera) = get(SpinStringNode(cam, "DeviceModelName", CameraTLDeviceNodeMap()))

"""
  show(::IO, ::Camera)

  Write details of camera to supplied IO.
"""
function show(io::IO, cam::Camera)
  vendorname = vendor(cam)
  modelname = model(cam)
  serialno = serial(cam)
  write(io, "$vendorname $modelname ($serialno)")
end

#
# Device status
#

"""
  devicetemperature(cam::Camera, location::String) -> Float

  Return the temperature of the specified device location.
"""
function devicetemperature(cam::Camera, location::String)
  set!(SpinEnumNode(cam, "DeviceTemperatureSelector"), location)
  return Spinnaker.get(Spinnaker.SpinFloatNode(cam, "DeviceTemperature"))
end

#
# Image acquistion
#

"""
  acquisitionmode(::Camera) -> String

  Return camera acquistion mode.
"""
acquisitionmode(cam::Camera) = get(SpinEnumNode(cam, "AcquisitionMode"))

"""
  acquisitionmode!(::Camera, ::AbstractString) -> String

  Set camera acquistion mode, returns set mode.
"""
acquisitionmode!(cam::Camera, mode) = set!(SpinEnumNode(cam, "AcquisitionMode"), mode)

function _isimagecomplete(himage_ref)
  isIncomplete = Ref(bool8_t(false))
  spinImageIsIncomplete(himage_ref[], isIncomplete);
  if isIncomplete == true
    imageStatus = Ref(spinImageStatus(IMAGE_NO_ERROR))
    spinImageGetStatus(himage_ref[], imageStatus)
    spinImageRelease(himage_ref[])
    @warn "Image incomplete with error $(imageStatus)"
    return false
  else
    return true
  end
end

#
# Image retrieval -> SpinImage
#

"""
  getimage(::Camer; release=true, timeout=-1) -> Image

  Copy the next image from the specified camera, blocking until available unless a timeout of >= 0 (ms) is specified. If release
  is false, the image buffer is not released.
"""
getimage(cam::Camera; release=true, timeout=-1) = getimage!(cam, SpinImage(), release=release, timeout=timeout)

"""
  getimage!(::Camera, ::SpinImage; release=true, timeout=-1) -> Image

  Copy the next image from the specified camera, blocking until available unless a timeout of >= 0 (ms) is specified, overwriting existing.
  If releaseis false, the image buffer is not released.
"""
function getimage!(cam::Camera, image::SpinImage; release=true, timeout=-1)

  # Get image handle and check it's complete
  himage_ref = Ref(spinImage(C_NULL))
  if timeout == -1
    spinCameraGetNextImage(cam, himage_ref);
  else
    spinCameraGetNextImageEx(cam, timeout, himage_ref);
  end
  @assert _isimagecomplete(himage_ref)

  # Create output image, copy and release buffer
  spinImageDeepCopy(himage_ref[], image)
  if release
    spinImageRelease(himage_ref[])
  end

  return image

end

#
# Image retrieval -> CameraImage
#

"""
  getimage(::Camera, ::Type{T}; normalize=true; release=true, timeout=-1) -> CameraImage

  Copy the next image from the specified camera, converting the image data to the specified array
  format, blocking until available unless a timeout of >= 0 (ms) is specified. If release is false, the image buffer is not released.

  If `normalize == false`, the input data from the camera is interpreted as a number in the range of
  the underlying type, e.g., for a camera operating in Mono8 pixel format, a call
  `getimage!(cam, Float64, normalize=false)` will return an array of dobule precision numbers in
  the range [0, 255]. `If normalize == true` the input data is interpreted as an associated fixed point
  format, and thus the array will be in the range [0,1].

  To return images compatible with Images.jl, one can request a Gray value, e.g.,
  `getimage!(cam, Gray{N0f8}, normalize=true)`.

  Function also returns image ID and timestamp metadata.
"""
function getimage(cam::Camera, ::Type{T}; normalize=true, release=true, timeout=-1) where T

  himage_ref, width, height, id, timestamp, exposure = _pullim(cam, timeout=timeout)
  imdat = Array{T,2}(undef, (width,height))
  camim = CameraImage(imdat, id, timestamp, exposure)
  _copyimage!(himage_ref[], width, height, camim, normalize)
  if release
    spinImageRelease(himage_ref[])
  end
  return camim

end

"""
  getimage!(::Camera, ::CameraImage{T,2}; normalize=false; release=true, timeout=-1)

  Copy the next iamge from the specified camera, converting to the format of, and overwriting the
  provided CameraImage, blocking until available unless a timeout of >= 0 (ms) is specified. If release is
  false, the image buffer is not released.

  If `normalize == false`, the input data from the camera is interpreted as a number in the range of
  the underlying type, e.g., for a camera operating in Mono8 pixel format, a call
  `getimage!(cam, Float64, normalize=false)` will return an array of dobule precision numbers in
  the range [0, 255]. `If normalize == true` the input data is interpreted as an associated fixed point
  format, and thus the array will be in the range [0,1].

  To return images compatible with Images.jl, one can request a Gray value, e.g.,
  `getimage!(cam, Gray{N0f8}, normalize=true)`.
"""
function getimage!(cam::Camera, image::CameraImage{T,2}; normalize=true, release=true, timeout=-1) where T

  himage_ref, width, height, id, timestamp, exposure = _pullim(cam, timeout=timeout)
  camim = CameraImage(image.data, id, timestamp, exposure)
  _copyimage!(himage_ref[], width, height, camim, normalize)
  if release
    spinImageRelease(himage_ref[])
  end
  return camim

end

function _pullim(cam::Camera;timeout=-1)

  # Get image handle and check it's complete
  himage_ref = Ref(spinImage(C_NULL))
  if timeout == -1
    spinCameraGetNextImage(cam, himage_ref);
  else
    spinCameraGetNextImageEx(cam, timeout, himage_ref);
  end
  if !_isimagecomplete(himage_ref)
    spinImageRelease(himage_ref[])
    throw(ErrorException("Image not complete"))
  end

  # Get image dimensions, ID and timestamp
  width = Ref(Csize_t(0))
  height = Ref(Csize_t(0))
  id = Ref(Int64(0))
  timestamp = Ref(Int64(0))
  exposure = Ref(Float64(0))
  spinImageGetWidth(himage_ref[], width)
  spinImageGetHeight(himage_ref[], height)
  spinImageChunkDataGetIntValue(himage_ref[], "ChunkFrameID", id);
  spinImageChunkDataGetFloatValue(himage_ref[], "ChunkExposureTime", exposure);
  spinImageChunkDataGetIntValue(himage_ref[], "ChunkTimestamp", timestamp)
  return himage_ref, Int(width[]), Int(height[]), id[], timestamp[], exposure[]

end

#
# Image retrieval -> Array
#

"""
  getimage!(::Camera, ::AbstractArray{T,2}; normalize=false, relase=true, timeout=-1)

  Copy the next iamge from the specified camera, converting to the format of, and overwriting the
  provided abstract array, blocking until available unless a timeout of >= 0 (ms) is specified. If release
  is false, the image buffer is not released.

  If `normalize == false`, the input data from the camera is interpreted as a number in the range of
  the underlying type, e.g., for a camera operating in Mono8 pixel format, a call
  `getimage!(cam, Array{Float64}(undef, dims...), normalize=false)` will return an array of dobule
  precision numbers in the range [0, 255]. `If normalize == true` the input data is interpreted as
  an associated fixed point format, and thus the array will be in the range [0,1].
"""
function getimage!(cam::Camera, image::AbstractArray{T,2}; normalize=true, release=true, timeout=-1) where T

  himage_ref, width, height, id, timestamp, exposure = _pullim(cam, timeout=timeout)
  _copyimage!(himage_ref[], width, height, image, normalize)
  if release
    spinImageRelease(himage_ref[])
  end
  return id, timestamp, exposure

end

#
# Image retrieval -> File
#

"""
    saveimage()::Camera, fn::AbstractString, ::spinImageFileFormat; release=true, timeout=-1)

    Save the next image from the specified camera to file `fn`, blocking until available unless
    a timeout of >= 0 (ms) is specified. If release is false, the image buffer is not released.
"""
function saveimage(cam::Camera, fn::AbstractString, fmt::spinImageFileFormat; relase=true, timeout=-1)

    # Get image handle and check it's complete
    himage_ref = Ref(spinImage(C_NULL))
    if timeout == -1
      spinCameraGetNextImage(cam, himage_ref);
    else
      spinCameraGetNextImageEx(cam, timeout, himage_ref);
    end
    @assert _isimagecomplete(himage_ref)
    spinImageSave(himage_ref[], fn, fmt)
    if release
      spinImageRelease(himage_ref[])
    end

end


text:
"""
Convert edge table into set of synapse-key dictionaries.

Output:
    dict with presynaptic segments
    dict with postsynaptic segments
    dict with synaptic centroid coordinates
"""
function edges_to_syn_dicts(tbl)
    pre = Dict()
    post = Dict()
    for i in 1:size(tbl,1)
        synID = tbl[i,1]
        pre[synID] = tbl[i,2][1]
        post[synID] = tbl[i,2][2]
    end
    return pre, post
end

"""
Find matching synapse in list of synapse locations given location of one synapse

Input:
    coordinate of particular synapse
    list of synase coordinates
Output:
    index of most similar synapse coordinate

For now, assume that synapse coordinates exactly line up
"""
function find_matching_synapse(syn_A, list_of_syn_B)
    return findfirst(list_of_syn_B, syn_A)
end

"""
For two sets of synapses labels on a given volume, map one set to the other

Input:
	synapse set A IDs with centroid location as 3-element coordinate (Ax2 array)
	synapse set B IDs with centroid location as 3-element coordinate (Bx2 array)

Output:
	Nx2 array:
		col 1: synapse set A IDs
		col 2: associated synapse set B IDs
		If there is no associated ID, that column entry will be 0
"""
function map_synapses(tbl_A, tbl_B)
    synA_to_synB = []
    for i in 1:size(tbl_A,1)
        ind_B = find_matching_synapse(tbl_A[i,3], tbl_B[:,3])
        if ind_B != 0
            push!(synA_to_synB, [tbl_A[i,1], tbl_B[ind_B,1]])
        else
            push!(synA_to_synB, [tbl_A[i,1], 0])
        end
    end

    for i in 1:size(tbl_B,1)
        ind_A = find_matching_synapse(tbl_B[i,3], tbl_A[:,3])
        if ind_A == 0
            push!(synA_to_synB, [0, tbl_B[i,1]])
        end
    end
    return hcat(synA_to_synB...)'
end

"""
Construct table for NRI from two edge tables (see Synaptor by N Turner)

Input:
    edge tables for two reconstructions

Output:
    (A+1)x(B+1) sparse count table built according to:
                '160618 - Computing NRI, Documentation.pdf'
        1st column & 1st row are FP
    lookup of row index to segment ID (considers row 2 to be index 1)
    lookup of col index to segment ID (considers col 2 to be index 1)

"""
function build_count_table(tbl_A, tbl_B)
    synA_to_synB = map_synapses(tbl_A, tbl_B)
    A_to_inds, inds_to_A = get_indexed_seg_IDs(tbl_A, 1)
    preA, postA = edges_to_syn_dicts(tbl_A)
    B_to_inds, inds_to_B = get_indexed_seg_IDs(tbl_B, 1)
    preB, postB = edges_to_syn_dicts(tbl_B)
    count_table = zeros(Int64, length(A_to_inds)+1, length(B_to_inds)+1)
    for k in 1:size(synA_to_synB,1)
        synA, synB = synA_to_synB[k,:]
        for (segA,segB) in [(preA,preB), (postA,postB)]
            if synA > 0
                i = A_to_inds[segA[synA]]
            else
                i = 1
            end
            if synB > 0
                j = B_to_inds[segB[synB]]
            else
                j = 1
            end
            count_table[i,j] += 1
        end
    end
    return count_table, A_to_inds, B_to_inds
end

"""
Simulate a merge error in the count table
"""
function merge_columns(count_table, i, j)
    @assert 1 < i <= size(count_table,1)
    @assert 1 < j <= size(count_table,2)
    # @assert i != j
    count_table[:,i] = count_table[:,i] + count_table[:,j]
    return count_table[:,[collect(1:j-1)...,collect(j+1:end)...]]
end

"""
Simulate a split error in the count table
"""
function split_column(count_table, j)
    @assert 1 < j <= size(count_table,2)
    # no need to unnecessarily split a segment that only has one synapse
    if sum(count_table[:,j]) > 1
        orig = count_table[:,j]
        split = count_table[:,j]
        syn_split = [rand(1:syn_count-1) for syn_count in orig[orig .> 1]]
        split[orig .> 1] = syn_split
        orig[orig .> 1] -= syn_split
        return hcat(count_table[:,1:j-1], orig, split, count_table[:,j+1:end])
    end
    return count_table
end

"""
Simulate a missed synapse in the count table
"""
function remove_synapse(count_table, pre_pair, post_pair)
    for (i,j) in [pre_pair, post_pair]
        @assert 1 < i <= size(count_table,1)
        @assert 1 < j <= size(count_table,2)
    end
    for (i,j) in [pre_pair, post_pair]
        if count_table[i,j] > 0
            count_table[i,j] -= 1
            count_table[i,1] += 1
        end
    end
    return count_table
end

"""
Simulate an added synapse in the count table
"""
function add_synapse(count_table, seg_pair)
    for j in seg_pair
        @assert 1 < j <= size(count_table,2)
    end
    for j in seg_pair
        count_table[1,j] += 1
    end
    return count_table
end


text:
macro methods(func)
  :(methods($func))
end

is_wsl() = Sys.islinux() && isfile("/mnt/c/Windows/System32/cmd.exe")
is_windows() = Sys.iswindows()
is_linux() = Sys.islinux()

"""
    path_mnt(path = ".")

Relative path will kept the original format.
"""
function path_mnt(path=".")
  # path = realpath(path)
  n = length(path)
  if is_wsl() && n >= 2 && path[2] == ':'
    pan = "/mnt/$(lowercase(path[1]))"
    path = n >= 3 ? "$pan$(path[3:end])" : pan
  elseif is_windows() && n >= 6 && path[1:5] == "/mnt/"
    pan = "$(uppercase(path[6])):"
    path = n >= 7 ? "$pan$(path[7:end])" : pan
  end
  path
end

"""
    open pdf file in SumatraPDF
"""
function show_pdf(file)
  app = "C:/Program Files/RStudio/bin/sumatra/SumatraPDF.exe"
  if is_wsl()
    app = path_mnt(app)
    run(`$app $file`; wait=false)
  elseif is_windows()
    run(`$app $file`; wait=false)
  end
  nothing
end

function show_file(file, verbose=false)
  file = abspath(file)
  cmd = `cmd /c "$file"`
  verbose && @show cmd
  run(cmd; wait=false)
  if !verbose
    return nothing
  end
end

function write_fig(file="Rplot.pdf", width=10, height=5; show=true)
  plot!(size=(width * 72, height * 72))
  savefig(file)
  if show
    show_file(file)
  end
end

"""
    merge_pdf("*.pdf", output="Plot.pdf")

Please install [pdftk](https://www.pdflabs.com/tools/pdftk-the-pdf-toolkit/) first.
On Linux, `sudo apt install pdftk-java`.

merge multiple pdf files by `pdftk`
"""
function merge_pdf(input, outfile="Plot.pdf"; is_del=false, show=true)
  # input = abspath(input)
  files = glob(input)
  id = str_extract(basename.(files), "\\d{1,}")
  id = parse.(Int32, id) |> sortperm
  files = files[id]

  run(`pdftk $files cat output $outfile`)
  show && show_file(outfile)
  is_del && run(`rm $files`)
  nothing
end

export @methods, is_wsl, is_windows, is_linux,
  path_mnt, show_pdf, show_file, write_fig


text: Dict{String,Any}("answer" => Dict{String,Any}("value" => "42","type" => "integer"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
struct MPCScenario <: AbstractScenario
    settings::Settings
    pvs::Array{MPCPV, 1}
    storage::Storage
    electric_tariff::MPCElectricTariff
    electric_load::MPCElectricLoad
    electric_utility::ElectricUtility
    financial::MPCFinancial
    generator::MPCGenerator
    cooling_load::MPCCoolingLoad
    limits::MPCLimits
end

"""
    MPCScenario(d::Dict)

Method for creating the MPCScenario struct:
```julia
    struct MPCScenario <: AbstractScenario
        settings::Settings
        pvs::Array{MPCPV, 1}
        storage::Storage
        electric_tariff::MPCElectricTariff
        electric_load::MPCElectricLoad
        electric_utility::ElectricUtility
        financial::MPCFinancial
        generator::MPCGenerator
        limits::MPCLimits
    end
```

The Dict `d` must have at a minimum the keys:
    - "ElectricLoad"
    - "ElectricTariff"

Other options include:
    - "PV", which can contain a Dict or Dict[]
    - "ElectricStorage"
    - "Generator"
    - "ElectricUtility"
    - "Settings"
    - "Financial"
    - "Limits"
"""
function MPCScenario(d::Dict)
    if haskey(d, "Settings")
        settings = Settings(;dictkeys_tosymbols(d["Settings"])...)
    else
        settings = Settings()
    end

    pvs = MPCPV[]
    if haskey(d, "PV")
        if typeof(d["PV"]) <: AbstractArray
            for (i, pv) in enumerate(d["PV"])
                if !(haskey(pv, "name"))
                    pv["name"] = string("PV", i)
                end
                push!(pvs, MPCPV(;dictkeys_tosymbols(pv)...))
            end
        elseif typeof(d["PV"]) <: AbstractDict
            push!(pvs, MPCPV(;dictkeys_tosymbols(d["PV"])...))
        else
            error("PV input must be Dict or Dict[].")
        end
    end

    if haskey(d, "Financial")
        financial = MPCFinancial(; dictkeys_tosymbols(d["Financial"])...)
    else
        financial = MPCFinancial()
    end

    if haskey(d, "ElectricUtility")
        electric_utility = ElectricUtility(; dictkeys_tosymbols(d["ElectricUtility"])...)
    else
        electric_utility = ElectricUtility()
    end

    if haskey(d, "ElectricStorage")
        # only modeling electrochemical storage so far
        storage_dict = Dict(dictkeys_tosymbols(d["ElectricStorage"]))
    else
        storage_dict = Dict(:size_kw => 0.0, :size_kwh => 0.0)
    end
    storage = Storage(Dict{String, AbstractStorage}("ElectricStorage" => MPCElectricStorage(; storage_dict...)))

    electric_load = MPCElectricLoad(; dictkeys_tosymbols(d["ElectricLoad"])...)

    electric_tariff = MPCElectricTariff(d["ElectricTariff"])

    if haskey(d, "Generator")
        generator = MPCGenerator(; dictkeys_tosymbols(d["Generator"])...)
    else
        generator = MPCGenerator(; size_kw=0)
    end

    # Placeholder/dummy cooling load set to zeros
    cooling_load = MPCCoolingLoad(; loads_kw_thermal = zeros(length(electric_load.loads_kw)), cop=1.0)
    if haskey(d, "Limits")
        limits = MPCLimits(; dictkeys_tosymbols(d["Limits"])...)
    else
        limits = MPCLimits()
    end

    return MPCScenario(
        settings,
        pvs,
        storage,
        electric_tariff,
        electric_load,
        electric_utility,
        financial,
        generator,
        cooling_load,
        limits
    )
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
using GraphQL, Test

@testset "Client" begin
  @info "Running client tests..."
  @test true # TODO: Implement more than placeholder test
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FirewallRuleProperties <: SwaggerModel
    startIpAddress::Any # spec type: Union{ Nothing, String } # spec name: startIpAddress
    endIpAddress::Any # spec type: Union{ Nothing, String } # spec name: endIpAddress

    function FirewallRuleProperties(;startIpAddress=nothing, endIpAddress=nothing)
        o = new()
        validate_property(FirewallRuleProperties, Symbol("startIpAddress"), startIpAddress)
        setfield!(o, Symbol("startIpAddress"), startIpAddress)
        validate_property(FirewallRuleProperties, Symbol("endIpAddress"), endIpAddress)
        setfield!(o, Symbol("endIpAddress"), endIpAddress)
        o
    end
end # type FirewallRuleProperties

const _property_map_FirewallRuleProperties = Dict{Symbol,Symbol}(Symbol("startIpAddress")=>Symbol("startIpAddress"), Symbol("endIpAddress")=>Symbol("endIpAddress"))
const _property_types_FirewallRuleProperties = Dict{Symbol,String}(Symbol("startIpAddress")=>"String", Symbol("endIpAddress")=>"String")
Base.propertynames(::Type{ FirewallRuleProperties }) = collect(keys(_property_map_FirewallRuleProperties))
Swagger.property_type(::Type{ FirewallRuleProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_FirewallRuleProperties[name]))}
Swagger.field_name(::Type{ FirewallRuleProperties }, property_name::Symbol) =  _property_map_FirewallRuleProperties[property_name]

function check_required(o::FirewallRuleProperties)
    true
end

function validate_property(::Type{ FirewallRuleProperties }, name::Symbol, val)
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export connect!
"""
  `connect!(fst::WFST, [i = get_initial(fst,single=true)])`

Same as `connect` but modifies `fst` in-place.
"""
function connect!(fst::WFST, i = get_initial(fst,single=true))
  if !isempty(fst)
    scc, c, v = get_scc(fst,i)
    idxs = findall( (|).((!).(v),(!).(c)) )
    return rm_state!(fst,idxs)
  else
    return fst
  end
end

export connect
"""
  `connect(fst::WFST, [i = get_initial(fst,single=true)])`

Remove states that do not belong to paths that end in final states when starting from `i`.
"""
function connect(fst::WFST, args...)
  fst_copy = deepcopy(fst)
  connect!(fst_copy, args...)
  return fst_copy
end


text:
directions = []

function reverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function goToAngle()
    while (!isborder(robot, Nord) || !isborder(robot, West))
        while !isborder(robot, Nord)
            move!(robot, Nord)
            push!(directions, Nord)
        end
        while !isborder(robot, West)
            move!(robot, West)
            push!(directions, West)
        end
    end
end

function putMarkersToAngle(side)
    while !isborder(robot, side)
        move!(robot, side)
    end
    putmarker!(robot)
end

function f()
    goToAngle()
    for i in [Sud, Ost, Nord, West]
        putMarkersToAngle(i)
    end

    while length(directions) > 0
        move!(robot, reverse(pop!(directions)))
    end
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
module DAI

#using BinDeps
#@BinDeps.load_dependencies
const libdai = Pkg.dir("DAI", "deps", "libdaiwrap.so")
dlopen(libdai)

include("dai_wrap.jl")

end