text:
"""Initiates a phase"""
function phase_run(j)
  info("initiating run phase")
  r = [@spawnat w runfun(j) for w in workers()]
  map(wait, r)
  info("run phase finished")
end

"""Runs a MapReduce job"""
function runfun(j)
  info("runfun worker $(myid()) initiated")
  length(j.inputs) > 2 && throw(ArgumentError("Cannot operate on more than two inputs."))
  length(j.inputs) == 1 && runfunsinglearg(j)
  info("runfun worker $(myid()) finished")
end

function runfunsinglearg(j)
  info("runfunsinglearg $(myid()) initiated")
  input_filename = j.inputs[1]
  inputfile = joinpath(splitdir, string(myid()), input_filename)
  io_input = open(inputfile)
  io_output = write_sink(j)
  for (linenum, line) in enumerate(eachline(io_input))
    try
      v = j.fun(line)
      v == "" && continue
      write(io_output, json(v) * "\n")
    catch e
      if typeof(e) in [MapperException, UDFException]
        warn("[JOB $(j.jobid)] mapper encountered an exception of type $(typeof(e)) on line $(linenum): $(e.msg)")
      else
        rethrow(e)
      end
    end
  end

  if length(ctx) > 0
    for (k, v) in ctx
      write(io_output, json((k, v)) * "\n")
    end
  end
  close(io_input)
  close(io_output)
  info("runfunsinglearg $(myid()) finished")
end

"""Runs a combiner"""
function runcombiner(filename, fun_combiner)
  r = [@spawnat w readlines(joinpath(splitdir, string(w), filename)) for w in workers()]
  combined = []
  map(wait, r)
  for p in r
    append!(combined, fetch(p))
  end
  fun_combiner(combined)
end

"""Runs a combiner and stores the result in an output file"""
function runcombiner(filename, fun_combiner, outputfilename)
  output = runcombiner(filename, fun_combiner)
  output_file_path = joinpath(datadir, outputfilename)
  isfile(output_file_path) && rm(output_file_path)
  fh = open(output_file_path, "w")
  for l in output
    write(fh, l * "\n")
  end
  close(fh)

  output
end


text:
using BenchmarkTools
nbit=16
Cb = control(nbit, (3,), 5=>X)
Pb = PutBlock{nbit}(CNOT, (3, 5))
@benchmark apply!($(copy(Reg)), $Cb)
@benchmark apply!($(copy(Reg)), $Pb)

ext: jl


text:
module Buttons

Base.Experimental.@optlevel 1

using Redux
using CImGui

# actions
abstract type AbstractButtonAction <: AbstractSyncAction end

get_label(a::AbstractButtonAction) = a.label

"""
    Rename(label, new_label)
Change button's label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractButtonAction
    label::String
    new_label::String
end

"""
    SetTriggeredTo(label, is_triggered)
Update widget's `is_triggered` state.
"""
struct SetTriggeredTo <: AbstractButtonAction
    label::String
    is_triggered::Bool
end

"""
    Resize(label, x, y)
Change button size to (x,y).
"""
struct Resize <: AbstractButtonAction
    label::String
    x::Cfloat
    y::Cfloat
end

"""
    ChangeWidth(label, w)
Change button width to `w`.
"""
struct ChangeWidth <: AbstractButtonAction
    label::String
    width::Cfloat
end

"""
    ChangeHeight(label, h)
Change button height to `h`.
"""
struct ChangeHeight <: AbstractButtonAction
    label::String
    height::Cfloat
end

"""
    AddButton(label::AbstractString, size = (0,0))
Add a [`Button`](@ref).
"""
struct AddButton <: AbstractButtonAction
    label::String
    size::Tuple{Cfloat,Cfloat}
end
AddButton(label::AbstractString) = AddButton(label, (0,0))

"""
    DeleteButton(label::AbstractString)
Delete the [`Button`](@ref).
"""
struct DeleteButton <: AbstractButtonAction
    label::String
end

# state
abstract type AbstractButtonState <: AbstractImmutableState end

get_label(s::AbstractButtonState) = s.label

"""
    Buttons.State(label::AbstractString)
    Buttons.State(label::AbstractString, size)
    Buttons.State(label::AbstractString, size, is_triggered)
A simple button state which contains a label a.k.a the identifier,
a size tuple (x,y) and a flag value `is_triggered` that records the state of
the latest poll events.
"""
struct State <: AbstractButtonState
    label::String
    size::Tuple{Cfloat,Cfloat}
    is_triggered::Bool
end
State(label::AbstractString, size) = State(label, size, false)
State(label::AbstractString) = State(label, (0,0))

get_size(s::State) = s.size
is_triggered(s::State) = s.is_triggered

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::State, a::Rename) = State(a.new_label, s.size, s.is_triggered)
reducer(s::State, a::Resize) = State(s.label, (a.x, a.y), s.is_triggered)
reducer(s::State, a::ChangeWidth) = State(s.label, (a.width, s.size[2]), s.is_triggered)
reducer(s::State, a::ChangeHeight) = State(s.label, (s.size[1], a.height), s.is_triggered)
reducer(s::State, a::SetTriggeredTo) = State(s.label, s.size, a.is_triggered)

reducer(s::Dict{String,<:AbstractButtonState}, a::AbstractButtonAction) =
    Dict(k => (get_label(v) == get_label(a) ? reducer(v, a) : v) for (k, v) in s)

reducer(s::Vector{<:AbstractButtonState}, a::AbstractButtonAction) = map(s) do s
    get_label(s) == get_label(a) ? reducer(s, a) : s
end
reducer(s::Vector{<:AbstractButtonState}, a::AddButton) = [s..., State(a.label, a.size, false)]
reducer(s::Vector{<:AbstractButtonState}, a::DeleteButton) = filter(s -> get_label(s) !== get_label(a), s)

# helper
"""
    Button(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
`get_state` is a router function that tells how to find the target state from `store`. `chain_action` is for chaining upstream actions.
"""
function Button(store::AbstractStore, get_state=Redux.get_state)
    s = get_state(store)
    is_triggered = CImGui.Button(s.label, CImGui.ImVec2(s.size...))
    dispatch!(store, SetTriggeredTo(s.label, is_triggered))
    return is_triggered
end

end # module


text:
"""
	wiggle(cmd0::String="", arg1=nothing; kwargs...)

Reads (length,azimuth) pairs from file and plot a windwiggle diagram.

Full option list at [`pswiggle`]($(GMTdoc)wiggle.html)

Parameters
----------

- $(GMT.opt_R)
- **Z** | **scale** :: [Type => Number | Str]

    Gives anomaly scale in data-units/distance-unit.
    ($(GMTdoc)wiggle.html#a)
- **A** | **azimuth** :: [Type => Str | number]

    Sets the preferred positive azimuth. Positive wiggles will “gravitate” towards that direction.
    ($(GMTdoc)wiggle.html#a)
- $(GMT.opt_B)
- **C** | **center** :: [Type => Number]

    Subtract center from the data set before plotting [0].
    ($(GMTdoc)wiggle.html#c)
- **D** | **scale_bar** :: [Type => Str]

    Defines the reference point on the map for the vertical scale bar using one of four coordinate systems.
    ($(GMTdoc)wiggle.html#d)
- **F** | **box** :: [Type => Str]

    Without further options, draws a rectangular border around the vertical scale bar.
    ($(GMTdoc)wiggle.html#f)
- **G** | **fill** :: [Type => Number | Str]

    Set fill shade, color or pattern for positive and/or negative wiggles [Default is no fill].
    ($(GMTdoc)wiggle.html#g)
- **I** | **fixed_azim** :: [Type => Number]

    Set a fixed azimuth projection for wiggles [Default uses track azimuth, but see -A].
    ($(GMTdoc)wiggle.html#i)
- $(GMT.opt_J)
- $(GMT.opt_P)
- **T** | **track** :: [Type => Number or Str | Tuple | []]

    Draw track [Default is no track]. Append pen attributes to use [Defaults: width = 0.25p, color =
    black, style = solid].
    ($(GMTdoc)wiggle.html#t)
- **W** | **pen** :: [Type => Number | Str | tuple | []]

    Specify outline pen attributes [Default is no outline].
    ($(GMTdoc)wiggle.html#w)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_w)
- $(GMT.opt_swap_xy)
"""
function wiggle(cmd0::String="", arg1=nothing; first=true, kwargs...)

    gmt_proggy = (IamModern[1]) ? "wiggle "  : "pswiggle "
	length(kwargs) == 0 && return monolitic(gmt_proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd, = parse_common_opts(d, cmd, [:c :e :f :g :p :t :w :F :UVXY :params], first)
	cmd  = parse_these_opts(cmd, d, [[:A :azimuth], [:C :center], [:I :fixed_azim], [:S], [:Z :scale]])

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R)

	cmd = parse_type_anchor(d, cmd, [:D :scale_bar],
                            (map=("g", nothing, 1), inside=("j", nothing, 1), anchor=("", arg2str, 2), width=("+w", arg2str), justify="+j", label_left="_+al", labels="+l", label="+l", offset=("+o", arg2str)), 'j')
	cmd *= opt_pen(d, 'T', [:T :track])
	cmd *= opt_pen(d, 'W', [:W :pen])
	cmd = add_opt_fill(cmd, d, [:G :fill], 'G')

	return finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1)
end

# ---------------------------------------------------------------------------------------------------
wiggle!(cmd0::String="", arg1=nothing; kw...) = wiggle(cmd0, arg1; first=false, kw...)
wiggle(arg1,  cmd0::String=""; kw...) = wiggle(cmd0, arg1; first=true, kw...)
wiggle!(arg1, cmd0::String=""; kw...) = wiggle(cmd0, arg1; first=false, kw...)

const pswiggle  = wiggle			# Alias
const pswiggle! = wiggle!			# Alias


text:
# Use
#    @warnpcfail precompile(args...)
# if you want to be warned when a precompile directive fails
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(render2file)),NamedTuple{(:id, :age, :field, :GPA, :courses), Tuple{Int64, Int64, String, Float64, Vector{String}}},typeof(render2file),String,Bool})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Vector{Float64},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(aml_create),Expr,Vector{Union{Expr, Symbol}},Vector{Any},Vector{Union{Expr, Symbol, Type}},Vector{Union{Expr, Symbol}},Vector{Union{Missing, String}},Vector{Union{Missing, Function, Symbol}},Vector{Union{Missing, Type}},String,Type,Array{Union{Missing, Function, Symbol}, 0},Bool,Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Expr})
    Base.precompile(Tuple{typeof(aml_create),Expr,Vector{Union{Expr, Symbol}},Vector{Any},Vector{Union{Expr, Symbol, Type}},Vector{Union{Expr, Symbol}},Vector{Union{Missing, String}},Vector{Union{Missing, Function, Symbol}},Vector{Union{Missing, Type}},String,Type,Array{Union{Missing, Function, Symbol}, 0},Bool,Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Vector{Union{Nothing, Expr}},Symbol})
    Base.precompile(Tuple{typeof(aml_parse),Expr})
    Base.precompile(Tuple{typeof(findcontent),Type{Vector{Float64}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Vector{Int64}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findfirsttext),Int64,Node})
    Base.precompile(Tuple{typeof(multiString),Vector{String}})
    Base.precompile(Tuple{typeof(updatecontent!),Dict{String, String},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Vector{Any},String,Node,Type{AbsNormal}})
end

max_stars_repo_name: aminya/AML


text:
export Reducer, PoolingReducer, CentroidReducer, AbstractReducer, reset!

# using OpenCV
using Statistics
using TiledIteration
using Hungarian
using Images
using ImageMorphology
using LinearAlgebra

abstract type AbstractReducer end

struct Reducer <: AbstractReducer
    reduct::Function
    parameters::Dict
end

"""
    PoolingReducer(f::Function, size::Int64)

Pooling reducer constructor.
"""
function PoolingReducer(f::Function, size::Int64)
    p = Dict("type"=>"pooling", "pooling_function"=>f, "size"=>size)
    Reducer(pooling_reduction, p)
end

"""
    CentroidReducer(
        n_centroids::Int64,
        n_in::Int64,
        img_size::Tuple{Int64,Int64}
    )

Centroid reducer constructor.
"""
function CentroidReducer(n_centroids::Int64, n_in::Int64,
                         img_size::Tuple{Int64,Int64})
    p = Dict(
        "type"=>"centroid",
        "n"=>n_centroids,
        "c_prev"=>Array{Array{Tuple{Float64,Float64},1},1}(undef, n_in),
        "a_prev"=>Array{Array{Int64,1},1}(undef, n_in),
        "img_size"=>img_size
    )
    Reducer(centroid_reduction, p)
end

function reset!(r::Reducer)
    if r.parameters["type"] == "centroid"
        n_in = length(r.parameters["c_prev"])
        r.parameters["c_prev"] = Array{Array{Tuple{Float64,Float64},1},1}(undef, n_in)
        r.parameters["a_prev"] = Array{Array{Int64,1},1}(undef, n_in)
    end
end

"""
    function Reducer(
        reducer_cfg::Dict;
        n_in::Int64=0,
        img_size::Tuple{Int64,Int64}=(0,0)
    )

General constructor for reducer.
"""
function Reducer(
    reducer_cfg::Dict;
    n_in::Int64=0,
    img_size::Tuple{Int64,Int64}=(0,0)
)
    reducer_type = reducer_cfg["type"]
    if reducer_type == "pooling"
        pooling_function = reducer_cfg["pooling_function"]
        if pooling_function == "mean"
            pf = Statistics.mean
        elseif pooling_function == "max"
            pf = maximum
        elseif pooling_function == "min"
            pf = minimum
        else
            throw(ArgumentError("Pooling function $pooling_function not implemented."))
        end
        reducer = PoolingReducer(pf, reducer_cfg["features_size"])
    elseif reducer_type == "centroid"
        reducer = CentroidReducer(reducer_cfg["n_centroids"], n_in, img_size)
    else
        throw(ArgumentError("Reducer type $reducer_type not implemented."))
    end
end

"""
    function remove_nan!(
        c::Array{Tuple{Float64,Float64},1},
        a::Array{Int64,1}
    )

Remove `NaN` values from centroids vector and apply the same removal to the
corresponding areas vector.
"""
function remove_nan!(
    c::Array{Tuple{Float64,Float64},1},
    a::Array{Int64,1}
)
    indexes = Int64[]
    for i in eachindex(c)
        if isnan(c[i][1]) || isnan(c[i][2])
            push!(indexes, i)
        end
    end
    deleteat!(c, indexes)
    deleteat!(a, indexes)
end

"""
    reorder_features(
        c1::Array{Tuple{Float64,Float64},1},
        a1::Array{Int64,1},
        c2::Array{Tuple{Float64,Float64},1},
        a2::Array{Int64,1}
    )

Reorder the features contained in c2, a2 by making a best match with c1, a1.
"""
function reorder_features(
    c1::Array{Tuple{Float64,Float64},1},
    a1::Array{Int64,1},
    c2::Array{Tuple{Float64,Float64},1},
    a2::Array{Int64,1}
)
    weights = zeros(length(c1), length(c1))
    for i in eachindex(c1)
        for j in eachindex(c2)
            weights[i, j] = norm(c1[i] .- c2[j]) + abs(a1[i] - a2[j])
        end
    end
    assignment, cost = hungarian(weights)  # 5.680 μs (40 allocations: 5.44 KiB)
    a2[assignment], c2[assignment]
end

"""
    normalized_centroids(
        c::Array{Tuple{Float64,Float64},1},
        img_size::Tuple{Int64,Int64}
    )

Normalize the input centroid vector to [0.0, 1.0].
Requires the image size as a 2D tuple.
"""
function normalized_centroids(c::Array{Tuple{Float64,Float64},1},
                              img_size::Tuple{Int64,Int64})
    c_norm = Array{Tuple{Float64,Float64},1}(undef, length(c))
    for i in eachindex(c)
        c_norm[i] = c[i] ./ img_size
    end
    c_norm
end

"""
    scaled_centroids(
        c::Array{Tuple{Float64,Float64},1},
        img_size::Tuple{Int64,Int64}
    )

Scale the input centroid vector from [0, 1] to the actual image size.
Requires the image size as a 2D tuple.
"""
function scaled_centroids(
    c::Array{Tuple{Float64,Float64},1},
    img_size::Tuple{Int64,Int64}
)
    c_scaled = Array{Tuple{Float64,Float64},1}(undef, length(c))
    for i in eachindex(c)
        c_scaled[i] = c[i] .* img_size
    end
    c_scaled
end

"""
    centroid_reduction(xs::Array{Array{UInt8,2},1}, parameters::Dict)

Generic centroid reduction function for several images (sequential application).
"""
function centroid_reduction(xs::Array{Array{UInt8,2},1}, parameters::Dict)
    fs = Array{Array{Tuple{Float64,Float64},1},1}(undef, length(xs))
    for i in eachindex(xs)
        if isdefined(parameters["c_prev"], i)
            c_prev = parameters["c_prev"][i]
            a_prev = parameters["a_prev"][i]
        else
            c_prev = nothing
            a_prev = nothing
        end
        c, a = centroid_reduction(xs[i], parameters["n"], c_prev, a_prev)
        fs[i] = normalized_centroids(c, parameters["img_size"])
        parameters["c_prev"][i] = c
        parameters["a_prev"][i] = a
    end
    return fs
end

"""
    centroid_reduction(
        x::Array{UInt8, 2},
        c_prev::Array{Tuple{Float64,Float64},1},
        a_prev::Array{Int64,1}
    )

Given an image, return the centroids and the boxes areas of the `n` largest
connected components, `n` being defined in the parameters dictionary.
Fill with zeros if there are less than `n` components.
"""
function centroid_reduction(
    x::Array{UInt8, 2},
    n::Int64,
    c_prev::Union{Array{Tuple{Float64,Float64},1}, Nothing},
    a_prev::Union{Array{Int64,1}, Nothing}
)
    labels = label_components(x)
    boxes = component_boxes(labels)
    centroids = component_centroids(labels)
    # popfirst!(centroids)  # remove largest
    areas = [abs(b[1][1]-b[2][1]-1) * abs(b[1][2]-b[2][2]-1) for b in boxes]
    p = sortperm(areas, rev=true)
    centroids = centroids[p]
    areas = areas[p]
    remove_nan!(centroids, areas)
    c = fill((0.0, 0.0), n)
    a = fill(0, n)
    for i in eachindex(centroids)
        c[i] = centroids[i]
        a[i] = areas[i]
        if i == n
            break
        end
    end
    if c_prev != nothing && a_prev != nothing
        # a, c = reorder_features(parameters["c_prev"], parameters["a_prev"], c, a)
        a, c = reorder_features(c_prev, a_prev, c, a)
        heatmp = implot(labels)
        display(heatmp)
    end
    c, a
end

"""
    pooling_reduction(xs::Array{Array{UInt8,2},1}, parameters::Dict)

Generic pooling function for several images (sequential application).
"""
function pooling_reduction(xs::Vector{Array{UInt8,2}}, parameters::Dict)
    [pooling_reduction(xs_i, parameters) for xs_i in xs]
end

"""
    pooling_reduction(x::Array{UInt8,2}, parameters::Dict)

Generic pooling function for a single image.
"""
function pooling_reduction(x::Array{UInt8,2}, parameters::Dict)
    outsz = (parameters["size"], parameters["size"])
    tilesz = ceil.(Int, size(x)./outsz)
    R = TileIterator(axes(x), tilesz)
    [parameters["pooling_function"](view(x, tileaxs...)) for tileaxs in R] ./ 255.0
end


text:
@testset "Storage data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type GenericBattery, consider changing the device models"
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5 = PSB.build_system(PSITestSystems, "c_sys5")
    c_sys14 = PSB.build_system(PSITestSystems, "c_sys14")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5)
    @test_logs (:warn, warn_message) match_mode = :any mock_construct_device!(
        op_problem,
        model,
    )
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys14)
    @test_logs (:warn, warn_message) match_mode = :any mock_construct_device!(
        op_problem,
        model,
    )
end

@testset "Storage Basic Storage With DC - PF" begin
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 72, 0, 72, 72, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage Basic Storage With AC - PF" begin
    model = DeviceModel(GenericBattery, BookKeeping)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 96, 0, 96, 96, 24, false)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation DC - PF" begin
    model = DeviceModel(GenericBattery, BookKeepingwReservation)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 96, 0, 72, 72, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation With AC - PF" begin
    model = DeviceModel(GenericBattery, BookKeepingwReservation)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 120, 0, 96, 96, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation DC - PF" begin
    model = DeviceModel(GenericBattery, BatteryAncillaryServices)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 96, 0, 168, 72, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "Storage with Reservation With AC - PF" begin
    model = DeviceModel(GenericBattery, BatteryAncillaryServices)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 120, 0, 192, 96, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "BatteryEMS with EnergyTarget with DC - PF" begin
    model = DeviceModel(BatteryEMS, EnergyTarget)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    op_problem = OperationsProblem(MockOperationProblem, DCPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 144, 0, 72, 72, 48, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "BatteryEMS with EnergyTarget With AC - PF" begin
    model = DeviceModel(BatteryEMS, EnergyTarget)
    c_sys5_bat = PSB.build_system(PSITestSystems, "c_sys5_bat_ems")
    op_problem = OperationsProblem(MockOperationProblem, ACPPowerModel, c_sys5_bat)
    mock_construct_device!(op_problem, model)
    moi_tests(op_problem, false, 168, 0, 96, 96, 48, true)
    psi_checkobjfun_test(op_problem, GAEVF)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1b-2b" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_b_sys")
    op_problem = OperationsProblem(
        EconomicDispatchProblem,
        template,
        c_sys5;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    moi_tests(op_problem, true, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(op_problem, [MOI.OPTIMAL], 5811.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1c-2c" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_c_sys")

    op_problem = OperationsProblem(
        EconomicDispatchProblem,
        template,
        c_sys5;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    moi_tests(op_problem, true, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(op_problem, [MOI.OPTIMAL], -63.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1d-2d" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_d_sys")

    op_problem = OperationsProblem(
        EconomicDispatchProblem,
        template,
        c_sys5;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    moi_tests(op_problem, true, 28, 0, 20, 16, 12, true)
    psi_checksolve_test(op_problem, [MOI.OPTIMAL], -11118.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1e-2e" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_e_sys")
    op_problem = OperationsProblem(
        EconomicDispatchProblem,
        template,
        c_sys5;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    moi_tests(op_problem, true, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(op_problem, [MOI.OPTIMAL], 5547.0, 10.0)
end

@testset "BatteryEMS with EnergyTarget Formulations (energy target - cases 1f-2f" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, BatteryEMS, EnergyTarget)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    c_sys5 = PSB.build_system(PSITestSystems, "batt_test_case_f_sys")

    op_problem = OperationsProblem(
        EconomicDispatchProblem,
        template,
        c_sys5;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    @test build!(op_problem; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT
    moi_tests(op_problem, true, 21, 0, 12, 9, 9, true)
    psi_checksolve_test(op_problem, [MOI.OPTIMAL], -1825.0, 10.0)
end


text:
#!/usr/bin/env julia
__precompile__(true)

module RHEOS

# installed from Julia package repository
using InverseLaplace
using NLopt
using JLD2
using DataStructures
using FunctionWrappers: FunctionWrapper
using Dierckx

import DSP.conv

# Base and stdlib imports
import Base: +, -, *,|
import DelimitedFiles: readdlm, writedlm

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

# This defines the data type for all arrays, parameters and processing
# it is defined as a const to avoid performance penalties.
# See julia docs for more info on this.
const RheoFloat = Float64

# Singularity approximation constant
const singularity_offset = 10.0

# convenience data types used as in many function parameters
const IntOrNone = Union{Integer, Nothing}
const RheovecOrNone = Union{Vector{RheoFloat}, Nothing}

######################################################################
export RheoFloat

# definitions.jl
export RheoLogItem, RheoLog, rheologrun, showlog
export RheoTimeData, RheoFreqData
export RheoTimeDataType, RheoFreqDataType, check_time_data_consistency
export LoadingType, strain_imposed, stress_imposed
export hastime, hasstress, hasstrain
export TimeDataType, time_only, strain_only, stress_only, strain_and_stress
export FreqDataType, invalid_freq_data, freq_only, with_modulus
export RheoModelClass, RheoModel, model_parameters
export rheoconv, invLaplace
export freeze_params
export relaxmod, creepcomp, storagemod, lossmod, dynamicmod

# IO.jl
export importcsv, exportcsv#, savedata, loaddata, savemodel, loadmodel

# models.jl
export null_modulus
export Springpot, Spring, Dashpot
export Fract_Maxwell, FractS_Maxwell, FractD_Maxwell, Maxwell
export Fract_KelvinVoigt, FractS_KelvinVoigt, FractD_KelvinVoigt, KelvinVoigt
export Fract_Zener, FractSLS_Zener, SLS_Zener, FractJeffreys_Zener, Jeffreys_Zener, FractSolid
export Fract_PT, FractSLS_PT, SLS_PT, Jeffreys_PT, FractJeffreys_PT
export BurgersLiquid
export SLS2, PowerLawPlateau

# datagen.jl
export timeline
export strainfunction, stressfunction
export hstep, ramp, stairs, square, sawtooth, triangle
export frequencyspec

# processing.jl
export resample, indexweight, cutting, smooth, extract
export modelfit, modelpredict, modelstepfit, modelsteppredict
export dynamicmodelfit, dynamicmodelpredict

#interface.jl
export Interface
export AFM, Tweezers

######################################################
# bundled dependencies from rheos-cambridge forked repos
MittLeffLiteDir = joinpath(@__DIR__, "..", "deps", "MittLeffLite", "MittLeffLite.jl")
include(MittLeffLiteDir)

include("base.jl")
include("definitions.jl")
include("IO.jl")
include("modeldatabase.jl")
include("datagen.jl")
include("processing.jl")
include("interface.jl")

end


text:
"""
    BoolDistribution(p_true)

Create a distribution over Boolean values (`true` or `false`).

`p_true` is the probability of the `true` outcome; the probability of `false` is 1-`p_true`.
"""
struct BoolDistribution
    p::Float64 # probability of true
end

pdf(d::BoolDistribution, s::Bool) = s ? d.p : 1.0-d.p

rand(rng::AbstractRNG, d::BoolDistribution) = rand(rng) <= d.p

Base.iterate(d::BoolDistribution) = ((d.p, true), true)
function Base.iterate(d::BoolDistribution, state::Bool)
    if state
        return  ((1.0 - d.p, false), false)
    else
        return nothing
    end
end

support(d::BoolDistribution) = [true, false]

==(d1::BoolDistribution, d2::BoolDistribution) = d1.p == d2.p

Base.hash(d::BoolDistribution) = hash(d.p)

Base.length(d::BoolDistribution) = 2


text:
# New version
"""
	CV_solve_time_independent(C_in::Array{T1,2},a_in::Vector{Float64},
													 h::Array{T2,2},
													 d::Array{T3,1},
													 Z::Array{T4,2},
													 t0::AbstractFloat,
													 tf::AbstractFloat;
													 tstep::AbstractFloat=0.,
													 tols::Vector{Float64}=[1e-6,1e-3],
													 alg = Tsit5()) where {T1 <: Number, T2 <: Number,  T3 <: Number, T4 <: Number}

Time independent solver for N-mode Gaussian continuous variable systems
using the covariance matrix approach. Defined in terms of the canonical
X and P quadratures for each, the covariance matrix elements are defined by
C_ab = <cc^⊺> - <c><c^⊺> where c^⊺ = [X_1,P_1,X_2,P_2,...,X_N,P_N] with
X = (a + a^†)/√2, P = -i(a - a^†)/√2 where "a" is the lowering operator
of a harmonic oscilaltor. Note that c is a column vector and c^⊺ is its
transpose, a row vector.

The total Hamiltonian is H = H_quad + H_lin = 1/2*(c^⊺*h*c) + d^⊺*c

## args
* C_in:  Array 2N × 2N, initial condition for the covariance matrix
* a_in:  vector 2N × 1, initial condition for the average value vector

* h:     2N × 2N matrix, describing the quadratic part of the Hamiltonian
         evolution of the system in the X_i, P_i basis, such that H_quad = 1/2*(c^⊺*h*c)
* d:     2N × 1 vector, describing the linear part of the Hamiltonian
         evolution of the system in the X_i, P_i basis, such that H_lin = d^⊺*c
* Z:     Array 2N × 2N, describes the dissipative evolution of the system.
         Given the master equation d/dt ρ = ∑_k γ_k*D[L_k]ρ with D[L_k]ρ the
         standard dissipator for the Lindblad operator L_k, we describe each
         Lindblad operator in the X_i, P_i basis as L_k = R_k^⊺*c Then
         Z = ∑_k γ_k*conj(R_k)*R_k^⊺

* t0:    float, start time
* tf:    float, end time
* tstep: float, time steps at which output data should be saved
* tols:  2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
* alg:   function, algorithm from DifferentialEquations for the solver to use, default is Tsit5

## returns
* sol_C.t, sol_C.u, sol_a.t, sol_a.u
"""
function CV_solve_time_independent(C_in::Array{T1,2},a_in::Vector{Float64},
													 h::Array{T2,2},
													 d::Array{T3,1},
													 Z::Array{T4,2},
													 t0::AbstractFloat,
													 tf::AbstractFloat;
													 tstep::AbstractFloat=0.,
													 tols::Vector{Float64}=[1e-6,1e-3],
													 alg = Tsit5()) where {T1 <: Number, T2 <: Number,  T3 <: Number, T4 <: Number}

	C_in = convert(Array{ComplexF64,2}, C_in)
	h = convert(Array{Float64,2}, h)
	d = convert(Array{Float64,1}, d)
	Z = convert(Array{ComplexF64,2}, Z)

	if tstep < 1e-10
        tstep = (tf-t0)./100
    end
    tspan = (t0,tf)

    σ = zeros(Float64,size(C_in,1),size(C_in,2))
    for nn = 1:2:size(C_in,1)
        σ[nn:1:nn+1,nn:1:nn+1] = [0. 1.;-1. 0.]
    end

    # Covariance matrix
    M = σ*(h + imag.(Z))
    Zs = σ*Z*transpose(σ)

    dif_C(du,u,p,t) = dCV!(du,u,M,Zs,t) # In place

    prob_C = ODEProblem{true}(dif_C,C_in,tspan)

    sol_C = solve(prob_C, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2])

    # Average
    d = σ*d

	dif_a(du,u,p,t) = d_av!(du,u,M,d,t) # In place

    prob_a = ODEProblem{true}(dif_a,a_in,tspan)

    sol_a = solve(prob_a, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2])

    return sol_C.t, sol_C.u, sol_a.t, sol_a.u
end

"""
    CV_solve_time_dependent(C_in::Array{T1,2},a_in::Vector{Float64},
												   h::Array{T2,2},
												   d::Function,
												   Z::Array{T3,2},
												   t0::AbstractFloat,
												   tf::AbstractFloat;
												   tstep::AbstractFloat=0.,
												   tols::Vector{Float64}=[1e-6,1e-3],
												   alg = Tsit5()) where {T1 <: Number, T2 <: Number,  T3 <: Number}

Time dependent (drive only) solver for N-mode Gaussian continuous variable systems using the covariance matrix approach.

args:
REQUIRED
    C_in: Array 2N × 2N, initial condition for the covariance matrix

    a_in: vector 2N × 1, initial condition for the average value vector

    h: 2N × 2N matrix, describing the quadratic part of the Hamiltonian evolution of the system in the X_i, P_i basis, such that
        H_quad = 1/2*(c^⊺*h*c)

    d: Function that returns as output a 2N × 1 vector, describing the linear part of the Hamiltonian evolution of the system in the X_i, P_i basis, such that
        H_lin = d^⊺*c

    Z: Array 2N × 2N, describes the dissipative evolution of the system. Given the master equation
    d/dt ρ = ∑_k γ_k*D[L_k]ρ
    with D[L_k]ρ the standard dissipator for the Lindblad operator L_k, we describe each Lindblad operator in the X_i, P_i basis as
    L_k = R_k^⊺*c
    Then
    Z = ∑_k γ_k*conj(R_k)*R_k^⊺

    t0: float, start time

    tf: float, end time
KEYWORD OPTIONAL
    tstep: float, time steps at which output data should be saved
    tols: 2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
    alg: function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
"""
function CV_solve_time_dependent(C_in::Array{T1,2},a_in::Vector{Float64},
												   h::Array{T2,2},
												   d::Function,
												   Z::Array{T3,2},
												   t0::AbstractFloat,
												   tf::AbstractFloat;
												   tstep::AbstractFloat=0.,
												   tols::Vector{Float64}=[1e-6,1e-3],
												   alg = Tsit5()) where {T1 <: Number, T2 <: Number,  T3 <: Number}

	C_in = convert(Array{ComplexF64,2}, C_in)
	h = convert(Array{Float64,2}, h)
	Z = convert(Array{ComplexF64,2}, Z)

	if tstep < 1e-10
        tstep = (tf-t0)./100
    end
    tspan = (t0,tf)

    σ = zeros(Float64,size(C_in,1),size(C_in,2))
    for nn = 1:2:size(C_in,1)
        σ[nn:1:nn+1,nn:1:nn+1] = [0. 1.;-1. 0.]
    end

    # Covariance matrix
    M = σ*(h + imag.(Z))
    Zs = σ*Z*transpose(σ)

    dif_C(du,u,p,t) = dCV!(du,u,M,Zs,t) # In place

    prob_C = ODEProblem{true}(dif_C,C_in,tspan)

    sol_C = solve(prob_C, alg, saveat = tstep, dense=false,
    										   save_everystep=false,
    										   abstol = tols[1],
    										   reltol = tols[2])

    # Average

    function drv(in_d::AbstractArray,t::AbstractFloat)
        in_d[:] = σ*d(t)
        nothing
    end

	drv_temp = Array{Float64}(undef,size(a_in,1))
	dif_a(du,u,p,t) = d_av!(du,u,M,drv,t,drv_temp) # In place

    prob_a = ODEProblem{true}(dif_a,a_in,tspan)

    sol_a = solve(prob_a, alg, saveat = tstep, dense=false, save_everystep=false, abstol = tols[1], reltol = tols[2])

    return sol_C.t, sol_C.u, sol_a.t, sol_a.u
end

"""
Time dependent solver for N-mode Gaussian continuous variable systems using the covariance matrix approach.

Input Parameters:
REQUIRED
    C_in: Array 2N × 2N, initial condition for the covariance matrix

    a_in: vector 2N × 1, initial condition for the average value vector

    h: Function that returns as output a 2N × 2N matrix, describing the quadratic part of the Hamiltonian evolution of the system in the X_i, P_i basis, such that
        H_quad = 1/2*(c^⊺*h*c)

    d: Function that returns as output a 2N × 1 vector, describing the linear part of the Hamiltonian evolution of the system in the X_i, P_i basis, such that
        H_lin = d^⊺*c

    Z: Array 2N × 2N, describes the dissipative evolution of the system. Given the master equation
    d/dt ρ = ∑_k γ_k*D[L_k]ρ
    with D[L_k]ρ the standard dissipator for the Lindblad operator L_k, we describe each Lindblad operator in the X_i, P_i basis as
    L_k = R_k^⊺*c
    Then
    Z = ∑_k γ_k*conj(R_k)*R_k^⊺

    t0: float, start time

    tf: float, end time
KEYWORD OPTIONAL
    tstep: float, time steps at which output data should be saved
    tols: 2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
    alg: function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
"""
function CV_solve_time_dependent(C_in::Array{T1,2},a_in::Vector{Float64},
												   h::Function,
												   d::Function,
												   Z::Array{T2,2},
												   t0::AbstractFloat,
												   tf::AbstractFloat;
												   tstep::AbstractFloat=0.,
												   tols::Vector{Float64}=[1e-6,1e-3],
												   alg = Tsit5()) where {T1 <: Number, T2 <: Number}

	C_in = convert(Array{ComplexF64,2}, C_in)
	Z = convert(Array{ComplexF64,2}, Z)

	if tstep < 1e-10
        tstep = (tf-t0)./100
    end
    tspan = (t0,tf)

    σ = zeros(Float64,size(C_in,1),size(C_in,2))
    for nn = 1:2:size(C_in,1)
        σ[nn:1:nn+1,nn:1:nn+1] = [0. 1.;-1. 0.]
    end

    # Covariance matrix
    function M(in_M::AbstractArray,t::AbstractFloat)
        in_M[:,:] = σ*(h(t) + imag.(Z))
        nothing
    end

    Zs = σ*Z*transpose(σ)

	M_temp = Array{Float64}(undef,size(C_in,1),size(C_in,2))
    dif_C(du,u,p,t) = dCV!(du,u,M,Zs,t,M_temp) # In place

    prob_C = ODEProblem{true}(dif_C,C_in,tspan)

    sol_C = solve(prob_C, alg, saveat = tstep, dense=false,
    										   save_everystep=false,
    										   abstol = tols[1],
    										   reltol = tols[2])

    # Average

    function drv(in_d::AbstractArray,t::AbstractFloat)
        in_d[:] = σ*d(t)
        nothing
    end

	fill!(M_temp,0.)
	drv_temp = Array{Float64}(undef,size(a_in,1))
	dif_a(du,u,p,t) = d_av!(du,u,M,drv,t,M_temp,drv_temp) # In place

    prob_a = ODEProblem{true}(dif_a,a_in,tspan)

    sol_a = solve(prob_a, alg, saveat = tstep, dense=false,
    										   save_everystep=false,
    										   abstol = tols[1],
    										   reltol = tols[2])
    return sol_C.t, sol_C.u, sol_a.t, sol_a.u
end

"""
Input in the a,a^† basis

NOTE: Due to the overhead of converting from the a,a^† basis to the X, P basis, benchmarking has these implementations at least a factor of 2 slower for the time dependent solvers, and requiring at least twice the memory than the X, P basis implementations.
"""

"""
Input Parameters:
REQUIRED
    C_in: Array 2N × 2N, initial condition for the covariance matrix

    a_in: vector 2N × 1, initial condition for the average value vector

    ω: N × 1 vector of the mode frequencies

    g: N × N array of the flip-flop coupling. Only the lower triangular part is used, and g[k,j] corresponds to the rate of the term a_k^†a_j in the Hamiltonian. The diagonal should be all zeros.

    λ: N × N array of the single and two-mode squeezing rates. Only the lower triangular part (including diagonal) is used, and λ[k,j] corresponds to the rate of the term a_k^†a_j^† in the Hamiltonian.

    α: N × 1 vector of the complex drive amplitudes, where α[k] corresponds to the term a_k in the Hamiltonian

    Z: Array 2N × 2N, describes the dissipative evolution of the system. Given the master equation
    d/dt ρ = ∑_k γ_k*D[L_k]ρ
    with D[L_k]ρ the standard dissipator for the Lindblad operator L_k, we describe each Lindblad operator in the X_i, P_i basis as
    L_k = R_k^⊺*c
    Then
    Z = ∑_k γ_k*conj(R_k)*R_k^⊺

    t0: float, start time

    tf: float, end time
KEYWORD OPTIONAL
    tstep: float, time steps at which output data should be saved
    tols: 2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
    alg: function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
"""
function CV_solve_time_independent(C_in::Array{T1,2},a_in::Vector{Float64},
													 ω::Vector{Float64},
													 g::Array{T2,2},
													 λ::Array{T3,2},
													 α::Array{T4,1},
													 Z::Array{T5,2},
													 t0::AbstractFloat,
													 tf::AbstractFloat;
													 tstep::AbstractFloat=0.,
													 tols::Vector{Float64}=[1e-6,1e-3],
													 alg = Tsit5()) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: Number, T5 <: Number}

	# Conert to X,P basis
	Num = length(ω)
	h = zeros(Float64,size(C_in,1),size(C_in,2))
	d = zeros(Float64,size(C_in,1))

	for jj = 1:1:Num
		for kk = (jj+1):1:Num
			# XX
			h[2*kk-1,2*jj-1] = (real(g[kk,jj]) + real(λ[kk,jj]))

			#XY
			h[2*kk-1,2*jj] = (imag(λ[kk,jj]) - imag(g[kk,jj]))

			#YX
			h[2*kk,2*jj-1] = (imag(λ[kk,jj]) + imag(g[kk,jj]))

			#YY
			h[2*kk,2*jj] = (real(g[kk,jj]) - real(λ[kk,jj]))
		end

		#YX same
		h[2*jj,2*jj-1] = 2*imag(λ[jj,jj])
	end

	h[:,:] = h[:,:] + h[:,:]'

	h[diagind(h)] = vec([(ω+2*real(diag(λ))) (ω-2*real(diag(λ)))]')[:]

	d[:] = sqrt(2)*vec([real.(α) -imag.(α)]')[:]

	return CV_solve_time_independent(C_in,a_in,h,d,Z,t0,tf; tstep=tstep,
															tols=tols,
															alg=alg)
end

"""
Input Parameters:
REQUIRED
    C_in: Array 2N × 2N, initial condition for the covariance matrix

    a_in: vector 2N × 1, initial condition for the average value vector

    ω: N × 1 vector of the mode frequencies

    g: N × N array of the flip-flop coupling. Only the lower triangular part is used, and g[k,j] corresponds to the rate of the term a_k^†a_j in the Hamiltonian. The diagonal should be all zeros.

    λ: N × N array of the single and two-mode squeezing rates. Only the lower triangular part (including diagonal) is used, and λ[k,j] corresponds to the rate of the term a_k^†a_j^† in the Hamiltonian.

    α: Function, returns an N × 1 vector of the complex drive amplitudes, where α[k] corresponds to the term a_k in the Hamiltonian

    Z: Array 2N × 2N, describes the dissipative evolution of the system. Given the master equation
    d/dt ρ = ∑_k γ_k*D[L_k]ρ
    with D[L_k]ρ the standard dissipator for the Lindblad operator L_k, we describe each Lindblad operator in the X_i, P_i basis as
    L_k = R_k^⊺*c
    Then
    Z = ∑_k γ_k*conj(R_k)*R_k^⊺

    t0: float, start time

    tf: float, end time
KEYWORD OPTIONAL
    tstep: float, time steps at which output data should be saved
    tols: 2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
    alg: function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
"""
function CV_solve_time_dependent(C_in::Array{T1,2},a_in::Vector{Float64},
												   ω::Vector{Float64},
												   g::Array{T2,2},
												   λ::Array{T3,2},
												   α::Function,
												   Z::Array{T4,2},
												   t0::AbstractFloat,
												   tf::AbstractFloat;
												   tstep::AbstractFloat=0.,
												   tols::Vector{Float64}=[1e-6,1e-3],
												   alg = Tsit5()) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: Number}

	# Conert to X,P basis
	Num = length(ω)
	h = zeros(Float64,size(C_in,1),size(C_in,2))

	for jj = 1:1:Num
		for kk = (jj+1):1:Num
			# XX
			h[2*kk-1,2*jj-1] = (real(g[kk,jj]) + real(λ[kk,jj]))

			#XY
			h[2*kk-1,2*jj] = (imag(λ[kk,jj]) - imag(g[kk,jj]))

			#YX
			h[2*kk,2*jj-1] = (imag(λ[kk,jj]) + imag(g[kk,jj]))

			#YY
			h[2*kk,2*jj] = (real(g[kk,jj]) - real(λ[kk,jj]))
		end

		#YX same
		h[2*jj,2*jj-1] = 2*imag(λ[jj,jj])
	end

    h[:,:] = h[:,:] + h[:,:]'

    h[diagind(h)] = vec([(ω+2*real(diag(λ))) (ω-2*real(diag(λ)))]')[:]

    function d(t::AbstractFloat)
        return sqrt(2)*vec([real.(α(t)) -imag.(α(t))]')[:]
    end

    return CV_solve_time_dependent(C_in,a_in,h,d,Z,t0,tf; tstep=tstep,tols=tols,alg=alg)
end

"""
The total Hamiltonian is H = ∑_k (ω_k*n_k + α_k*a_k + conj(α_k)*a_k^†) + ∑_kj (g_kj*a_k^†a_j + λ_kj*a_k^†a_j^† + h.c.)
with n = a^†a

Input Parameters:
REQUIRED
    C_in: Array 2N × 2N, initial condition for the covariance matrix

    a_in: vector 2N × 1, initial condition for the average value vector

    ω: Function, returns an N × 1 vector of the mode frequencies

    g: Function, returns an N × N array of the flip-flop coupling. Only the lower triangular part is used, and g[k,j] corresponds to the rate of the term a_k^†a_j in the Hamiltonian. The diagonal should be all zeros.

    λ: Function, returns an N × N array of the single and two-mode squeezing rates. Only the lower triangular part (including diagonal) is used, and λ[k,j] corresponds to the rate of the term a_k^†a_j^† in the Hamiltonian.

    α: Function, returns an N × 1 vector of the complex drive amplitudes, where α[k] corresponds to the term a_k in the Hamiltonian

    Z: Array 2N × 2N, describes the dissipative evolution of the system. Given the master equation
    d/dt ρ = ∑_k γ_k*D[L_k]ρ
    with D[L_k]ρ the standard dissipator for the Lindblad operator L_k, we describe each Lindblad operator in the X_i, P_i basis as
    L_k = R_k^⊺*c
    Then
    Z = ∑_k γ_k*conj(R_k)*R_k^⊺

    t0: float, start time

    tf: float, end time
KEYWORD OPTIONAL
    tstep: float, time steps at which output data should be saved
    tols: 2 x 1 array, vector of solver tolernaces in the order [abstol, reltol]
    alg: function, algorithm from DifferentialEquations for the solver to use, default is Tsit5
"""
function CV_solve_time_dependent(C_in::Array{T1,2},a_in::Vector{Float64},ω::Function,g::Function,λ::Function,α::Function,Z::Array{T2,2},t0::AbstractFloat,tf::AbstractFloat; tstep::AbstractFloat=0.,tols::Vector{Float64}=[1e-6,1e-3],alg = Tsit5()) where {T1 <: Number, T2 <: Number}

    # Conert to X,P basis
    Num = length(ω(0.))
    temp_h1 = zeros(Float64,size(C_in,1),size(C_in,2))
    gs1 = zeros(Float64,Num,Num)
    λs1 = zeros(Float64,Num,Num)

    # function h(t::AbstractFloat;temp_h=temp_h1,gs=gs1,λs=λs1)
    #   # fill!(temp_h,0.)
    #   gs[:,:] = g(t)
    #   λs[:,:] = λ(t)
    #
    #   for jj = 1:1:Num
    #       for kk = (jj+1):1:Num
    #           # XX
    #           temp_h[2*kk-1,2*jj-1] = (real(gs[kk,jj]) + real(λs[kk,jj]))
    #
    #           #XY
    #           temp_h[2*kk-1,2*jj] = (imag(λs[kk,jj]) - imag(gs[kk,jj]))
    #
    #           #YX
    #           temp_h[2*kk,2*jj-1] = (imag(λs[kk,jj]) + imag(gs[kk,jj]))
    #
    #           #YY
    #           temp_h[2*kk,2*jj] = (real(gs[kk,jj]) - real(λs[kk,jj]))
    #       end
    #
    #       #YX same
    #       temp_h[2*jj,2*jj-1] = 2*imag(λs[jj,jj])
    #   end
    #
    #   temp_h[:,:] = temp_h[:,:] + temp_h[:,:]'
    #
    #   # for ii = 1:1:Num
    #   #   temp_h[2*ii-1,2*ii-1] = ω(t)[ii] + 2*real(λ(t)[ii,ii])
    #   #   temp_h[2*ii,2*ii] = ω(t)[ii] - 2*real(λ(t)[ii,ii])
    #   # end
    #   temp_h[diagind(temp_h)] = vec([(ω(t)+2*real(diag(λs))) (ω(t)-2*real(diag(λs)))]')[:]
    #
    #   return temp_h
    # end

    function h(t::AbstractFloat)
        # fill!(temp_h,0.)
        gs = g(t)
        λs = λ(t)
        temp_h = zeros(size(gs,1)*2,size(gs,2)*2)

        for jj = 1:1:Num
            for kk = (jj+1):1:Num
                # XX
                temp_h[2*kk-1,2*jj-1] = (real(gs[kk,jj]) + real(λs[kk,jj]))

                #XY
                temp_h[2*kk-1,2*jj] = (imag(λs[kk,jj]) - imag(gs[kk,jj]))

                #YX
                temp_h[2*kk,2*jj-1] = (imag(λs[kk,jj]) + imag(gs[kk,jj]))

                #YY
                temp_h[2*kk,2*jj] = (real(gs[kk,jj]) - real(λs[kk,jj]))
            end

            #YX same
            temp_h[2*jj,2*jj-1] = 2*imag(λs[jj,jj])
        end

        temp_h[:,:] = temp_h[:,:] + temp_h[:,:]'

        # for ii = 1:1:Num
        #   temp_h[2*ii-1,2*ii-1] = ω(t)[ii] + 2*real(λ(t)[ii,ii])
        #   temp_h[2*ii,2*ii] = ω(t)[ii] - 2*real(λ(t)[ii,ii])
        # end
        temp_h[diagind(temp_h)] = vec([(ω(t)+2*real(diag(λs))) (ω(t)-2*real(diag(λs)))]')[:]

        return temp_h
    end

    function d(t::AbstractFloat)
        return sqrt(2)*vec([real.(α(t)) -imag.(α(t))]')[:]
    end

    return CV_solve_time_dependent(C_in,a_in,h,d,Z,t0,tf; tstep=tstep,tols=tols,alg=alg)
end


text:
using Mocking
Mocking.enable()

using Base.Test
using TimeZones
import TimeZones: PKG_DIR, TZDATA_DIR
import TimeZones.Olson: ZoneDict, RuleDict, tzparse, resolve

const TZFILE_DIR = joinpath(PKG_DIR, "test", "tzfile")

# For testing we'll reparse the tzdata every time to instead of using the serialized data.
# This should make the development/testing cycle simplier since you won't be forced to
# re-build the cache every time you make a change.
#
# Note: resolving only the time zones we want is much faster than running compile which
# recompiles all the time zones.
tzdata = Dict{AbstractString,Tuple{ZoneDict,RuleDict}}()
for name in ("asia", "australasia", "europe", "northamerica")
    tzdata[name] = tzparse(joinpath(TZDATA_DIR, name))
end

include("utils.jl")
include("time.jl")
include("Olson.jl")
include("utcoffset.jl")
include("types.jl")
include("interpret.jl")
include("accessors.jl")
include("arithmetic.jl")
include("io.jl")
include("tzfile.jl")
include("adjusters.jl")
include("conversions.jl")
include("ranges.jl")
include("local.jl")
include("local_mocking.jl")
include("discovery.jl")
VERSION >= v"0.5.0-dev+5244" && include("rounding.jl")
include("TimeZones.jl")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Low level module for mmap (memory mapping of files).
"""
module Mmap

import Base: OS_HANDLE, INVALID_OS_HANDLE

const PAGESIZE = Int(Sys.isunix() ? ccall(:jl_getpagesize, Clong, ()) : ccall(:jl_getallocationgranularity, Clong, ()))

# for mmaps not backed by files
mutable struct Anonymous <: IO
    name::AbstractString
    readonly::Bool
    create::Bool
end

"""
    Mmap.Anonymous(name::AbstractString="", readonly::Bool=false, create::Bool=true)

Create an `IO`-like object for creating zeroed-out mmapped-memory that is not tied to a file
for use in [`Mmap.mmap`](@ref Mmap.mmap). Used by `SharedArray` for creating shared memory arrays.

# Examples
```jldoctest
julia> using Mmap

julia> anon = Mmap.Anonymous();

julia> isreadable(anon)
true

julia> iswritable(anon)
true

julia> isopen(anon)
true
```
"""
Anonymous() = Anonymous("",false,true)

Base.isopen(::Anonymous) = true
Base.isreadable(::Anonymous) = true
Base.iswritable(a::Anonymous) = !a.readonly

# const used for zeroed, anonymous memory
gethandle(io::Anonymous) = INVALID_OS_HANDLE

# platform-specific mmap utilities
if Sys.isunix()

const PROT_READ     = Cint(1)
const PROT_WRITE    = Cint(2)
const MAP_SHARED    = Cint(1)
const MAP_PRIVATE   = Cint(2)
const MAP_ANONYMOUS = Cint(Sys.isbsd() ? 0x1000 : 0x20)
const F_GETFL       = Cint(3)

gethandle(io::IO) = RawFD(fd(io))

# Determine a stream's read/write mode, and return prot & flags appropriate for mmap
function settings(s::RawFD, shared::Bool)
    flags = shared ? MAP_SHARED : MAP_PRIVATE
    if s == INVALID_OS_HANDLE
        flags |= MAP_ANONYMOUS
        prot = PROT_READ | PROT_WRITE
    else
        mode = ccall(:fcntl, Cint, (RawFD, Cint, Cint...), s, F_GETFL)
        systemerror("fcntl F_GETFL", mode == -1)
        mode = mode & 3
        prot = (mode == 0) ? PROT_READ : ((mode == 1) ? PROT_WRITE : (PROT_READ | PROT_WRITE))
        if prot & PROT_READ == 0
            throw(ArgumentError("mmap requires read permissions on the file (open with \"r+\" mode to override)"))
        end
    end
    return prot, flags, (prot & PROT_WRITE) > 0
end

# Before mapping, grow the file to sufficient size
# Note: a few mappable streams do not support lseek. When Julia
# supports structures in ccall, switch to fstat.
grow!(::Anonymous,o::Integer,l::Integer) = return
function grow!(io::IO, offset::Integer, len::Integer)
    pos = position(io)
    filelen = filesize(io)
    if filelen < offset + len
        failure = ccall(:jl_ftruncate, Cint, (Cint, Int64), fd(io), offset+len)
        Base.systemerror(:ftruncate, failure != 0)
    end
    seek(io, pos)
    return
end

elseif Sys.iswindows()

const DWORD = Culong

const PAGE_READONLY          = DWORD(0x02)
const PAGE_READWRITE         = DWORD(0x04)
const PAGE_WRITECOPY         = DWORD(0x08)
const PAGE_EXECUTE_READ      = DWORD(0x20)
const PAGE_EXECUTE_READWRITE = DWORD(0x40)
const PAGE_EXECUTE_WRITECOPY = DWORD(0x80)
const FILE_MAP_COPY          = DWORD(0x01)
const FILE_MAP_WRITE         = DWORD(0x02)
const FILE_MAP_READ          = DWORD(0x04)
const FILE_MAP_EXECUTE       = DWORD(0x20)

function gethandle(io::IO)
    handle = Libc._get_osfhandle(RawFD(fd(io)))
    Base.windowserror(:mmap, handle == INVALID_OS_HANDLE)
    return handle
end

settings(sh::Anonymous) = sh.name, sh.readonly, sh.create
settings(io::IO) = Ptr{Cwchar_t}(0), isreadonly(io), true

else
    error("mmap not defined for this OS")
end # os-test

# core implementation of mmap

"""
    Mmap.mmap(io::Union{IOStream,AbstractString,Mmap.AnonymousMmap}[, type::Type{Array{T,N}}, dims, offset]; grow::Bool=true, shared::Bool=true)
    Mmap.mmap(type::Type{Array{T,N}}, dims)

Create an `Array` whose values are linked to a file, using memory-mapping. This provides a
convenient way of working with data too large to fit in the computer's memory.

The type is an `Array{T,N}` with a bits-type element of `T` and dimension `N` that
determines how the bytes of the array are interpreted. Note that the file must be stored in
binary format, and no format conversions are possible (this is a limitation of operating
systems, not Julia).

`dims` is a tuple or single [`Integer`](@ref) specifying the size or length of the array.

The file is passed via the stream argument, either as an open [`IOStream`](@ref) or filename string.
When you initialize the stream, use `"r"` for a "read-only" array, and `"w+"` to create a
new array used to write values to disk.

If no `type` argument is specified, the default is `Vector{UInt8}`.

Optionally, you can specify an offset (in bytes) if, for example, you want to skip over a
header in the file. The default value for the offset is the current stream position for an
`IOStream`.

The `grow` keyword argument specifies whether the disk file should be grown to accommodate
the requested size of array (if the total file size is < requested array size). Write
privileges are required to grow the file.

The `shared` keyword argument specifies whether the resulting `Array` and changes made to it
will be visible to other processes mapping the same file.

For example, the following code

```julia
# Create a file for mmapping
# (you could alternatively use mmap to do this step, too)
using Mmap
A = rand(1:20, 5, 30)
s = open("/tmp/mmap.bin", "w+")
# We'll write the dimensions of the array as the first two Ints in the file
write(s, size(A,1))
write(s, size(A,2))
# Now write the data
write(s, A)
close(s)

# Test by reading it back in
s = open("/tmp/mmap.bin")   # default is read-only
m = read(s, Int)
n = read(s, Int)
A2 = Mmap.mmap(s, Matrix{Int}, (m,n))
```

creates a `m`-by-`n` `Matrix{Int}`, linked to the file associated with stream `s`.

A more portable file would need to encode the word size -- 32 bit or 64 bit -- and endianness
information in the header. In practice, consider encoding binary data using standard formats
like HDF5 (which can be used with memory-mapping).
"""
function mmap(io::IO,
              ::Type{Array{T,N}}=Vector{UInt8},
              dims::NTuple{N,Integer}=(div(filesize(io)-position(io),sizeof(T)),),
              offset::Integer=position(io); grow::Bool=true, shared::Bool=true) where {T,N}
    # check inputs
    isopen(io) || throw(ArgumentError("$io must be open to mmap"))
    isbitstype(T)  || throw(ArgumentError("unable to mmap $T; must satisfy isbitstype(T) == true"))

    len = prod(dims) * sizeof(T)
    len >= 0 || throw(ArgumentError("requested size must be ≥ 0, got $len"))
    len == 0 && return Array{T}(undef, ntuple(x->0,Val(N)))
    len < typemax(Int) - PAGESIZE || throw(ArgumentError("requested size must be < $(typemax(Int)-PAGESIZE), got $len"))

    offset >= 0 || throw(ArgumentError("requested offset must be ≥ 0, got $offset"))

    # shift `offset` to start of page boundary
    offset_page::Int64 = div(offset, PAGESIZE) * PAGESIZE
    # add (offset - offset_page) to `len` to get total length of memory-mapped region
    mmaplen = (offset - offset_page) + len

    file_desc = gethandle(io)
    # platform-specific mmapping
    @static if Sys.isunix()
        prot, flags, iswrite = settings(file_desc, shared)
        iswrite && grow && grow!(io, offset, len)
        # mmap the file
        ptr = ccall(:jl_mmap, Ptr{Cvoid}, (Ptr{Cvoid}, Csize_t, Cint, Cint, RawFD, Int64),
            C_NULL, mmaplen, prot, flags, file_desc, offset_page)
        systemerror("memory mapping failed", reinterpret(Int, ptr) == -1)
    else
        name, readonly, create = settings(io)
        szfile = convert(Csize_t, len + offset)
        readonly && szfile > filesize(io) && throw(ArgumentError("unable to increase file size to $szfile due to read-only permissions"))
        handle = create ? ccall(:CreateFileMappingW, stdcall, Ptr{Cvoid}, (OS_HANDLE, Ptr{Cvoid}, DWORD, DWORD, DWORD, Cwstring),
                                file_desc, C_NULL, readonly ? PAGE_READONLY : PAGE_READWRITE, szfile >> 32, szfile & typemax(UInt32), name) :
                          ccall(:OpenFileMappingW, stdcall, Ptr{Cvoid}, (DWORD, Cint, Cwstring),
                                readonly ? FILE_MAP_READ : FILE_MAP_WRITE, true, name)
        Base.windowserror(:mmap, handle == C_NULL)
        ptr = ccall(:MapViewOfFile, stdcall, Ptr{Cvoid}, (Ptr{Cvoid}, DWORD, DWORD, DWORD, Csize_t),
                    handle, readonly ? FILE_MAP_READ : FILE_MAP_WRITE, offset_page >> 32, offset_page & typemax(UInt32), mmaplen)
        Base.windowserror(:mmap, ptr == C_NULL)
    end # os-test
    # convert mmapped region to Julia Array at `ptr + (offset - offset_page)` since file was mapped at offset_page
    A = unsafe_wrap(Array, convert(Ptr{T}, UInt(ptr) + UInt(offset - offset_page)), dims)
    finalizer(A) do x
        @static if Sys.isunix()
            systemerror("munmap",  ccall(:munmap, Cint, (Ptr{Cvoid}, Int), ptr, mmaplen) != 0)
        else
            status = ccall(:UnmapViewOfFile, stdcall, Cint, (Ptr{Cvoid},), ptr)!=0
            status |= ccall(:CloseHandle, stdcall, Cint, (Ptr{Cvoid},), handle)!=0
            Base.windowserror(:UnmapViewOfFile, status == 0)
        end
    end
    return A
end

mmap(file::AbstractString,
     ::Type{T}=Vector{UInt8},
     dims::NTuple{N,Integer}=(div(filesize(file),sizeof(eltype(T))),),
     offset::Integer=Int64(0); grow::Bool=true, shared::Bool=true) where {T<:Array,N} =
    open(io->mmap(io, T, dims, offset; grow=grow, shared=shared), file, isfile(file) ? "r" : "w+")::Array{eltype(T),N}

# using a length argument instead of dims
mmap(io::IO, ::Type{T}, len::Integer, offset::Integer=position(io); grow::Bool=true, shared::Bool=true) where {T<:Array} =
    mmap(io, T, (len,), offset; grow=grow, shared=shared)
mmap(file::AbstractString, ::Type{T}, len::Integer, offset::Integer=Int64(0); grow::Bool=true, shared::Bool=true) where {T<:Array} =
    open(io->mmap(io, T, (len,), offset; grow=grow, shared=shared), file, isfile(file) ? "r" : "w+")::Vector{eltype(T)}

# constructors for non-file-backed (anonymous) mmaps
mmap(::Type{T}, dims::NTuple{N,Integer}; shared::Bool=true) where {T<:Array,N} = mmap(Anonymous(), T, dims, Int64(0); shared=shared)
mmap(::Type{T}, i::Integer...; shared::Bool=true) where {T<:Array} = mmap(Anonymous(), T, convert(Tuple{Vararg{Int}},i), Int64(0); shared=shared)

"""
    Mmap.mmap(io, BitArray, [dims, offset])

Create a [`BitArray`](@ref) whose values are linked to a file, using memory-mapping; it has the same
purpose, works in the same way, and has the same arguments, as [`mmap`](@ref Mmap.mmap), but
the byte representation is different.

# Examples
```jldoctest
julia> using Mmap

julia> io = open("mmap.bin", "w+");

julia> B = Mmap.mmap(io, BitArray, (25,30000));

julia> B[3, 4000] = true;

julia> Mmap.sync!(B);

julia> close(io);

julia> io = open("mmap.bin", "r+");

julia> C = Mmap.mmap(io, BitArray, (25,30000));

julia> C[3, 4000]
true

julia> C[2, 4000]
false

julia> close(io)

julia> rm("mmap.bin")
```
This creates a 25-by-30000 `BitArray`, linked to the file associated with stream `io`.
"""
function mmap(io::IOStream, ::Type{<:BitArray}, dims::NTuple{N,Integer},
              offset::Int64=position(io); grow::Bool=true, shared::Bool=true) where N
    n = prod(dims)
    nc = Base.num_bit_chunks(n)
    chunks = mmap(io, Vector{UInt64}, (nc,), offset; grow=grow, shared=shared)
    if !isreadonly(io)
        chunks[end] &= Base._msk_end(n)
    else
        if chunks[end] != chunks[end] & Base._msk_end(n)
            throw(ArgumentError("the given file does not contain a valid BitArray of size $(join(dims, 'x')) (open with \"r+\" mode to override)"))
        end
    end
    B = BitArray{N}(undef, ntuple(i->0,Val(N))...)
    B.chunks = chunks
    B.len = n
    if N != 1
        B.dims = dims
    end
    return B
end

mmap(file::AbstractString, ::Type{T}, dims::NTuple{N,Integer}, offset::Integer=Int64(0);grow::Bool=true, shared::Bool=true) where {T<:BitArray,N} =
    open(io->mmap(io, T, dims, offset; grow=grow, shared=shared), file, isfile(file) ? "r" : "w+")::BitArray{N}

# using a length argument instead of dims
mmap(io::IO, ::Type{T}, len::Integer, offset::Integer=position(io); grow::Bool=true, shared::Bool=true) where {T<:BitArray} =
    mmap(io, T, (len,), offset; grow=grow, shared=shared)
mmap(file::AbstractString, ::Type{T}, len::Integer, offset::Integer=Int64(0); grow::Bool=true, shared::Bool=true) where {T<:BitArray} =
    open(io->mmap(io, T, (len,), offset; grow=grow, shared=shared), file, isfile(file) ? "r" : "w+")::BitVector

# constructors for non-file-backed (anonymous) mmaps
mmap(::Type{T}, dims::NTuple{N,Integer}; shared::Bool=true) where {T<:BitArray,N} = mmap(Anonymous(), T, dims, Int64(0); shared=shared)
mmap(::Type{T}, i::Integer...; shared::Bool=true) where {T<:BitArray} = mmap(Anonymous(), T, convert(Tuple{Vararg{Int}},i), Int64(0); shared=shared)

# msync flags for unix
const MS_ASYNC = 1
const MS_INVALIDATE = 2
const MS_SYNC = 4

"""
    Mmap.sync!(array)

Forces synchronization between the in-memory version of a memory-mapped `Array` or
[`BitArray`](@ref) and the on-disk version.
"""
function sync!(m::Array, flags::Integer=MS_SYNC)
    offset = rem(UInt(pointer(m)), PAGESIZE)
    ptr = pointer(m) - offset
    mmaplen = sizeof(m) + offset
    GC.@preserve m @static if Sys.isunix()
        systemerror("msync",
                    ccall(:msync, Cint, (Ptr{Cvoid}, Csize_t, Cint), ptr, mmaplen, flags) != 0)
    else
        Base.windowserror(:FlushViewOfFile,
            ccall(:FlushViewOfFile, stdcall, Cint, (Ptr{Cvoid}, Csize_t), ptr, mmaplen) == 0)
    end
end
sync!(B::BitArray, flags::Integer=MS_SYNC) = sync!(B.chunks, flags)

@static if Sys.isunix()
const MADV_NORMAL = 0
const MADV_RANDOM = 1
const MADV_SEQUENTIAL = 2
const MADV_WILLNEED = 3
const MADV_DONTNEED = 4
if Sys.islinux()
    const MADV_FREE = 8
    const MADV_REMOVE = 9
    const MADV_DONTFORK = 10
    const MADV_DOFORK = 11
    const MADV_MERGEABLE = 12
    const MADV_UNMERGEABLE = 13
    const MADV_HUGEPAGE = 14
    const MADV_NOHUGEPAGE = 15
    const MADV_DONTDUMP = 16
    const MADV_DODUMP = 17
    const MADV_WIPEONFORK = 18
    const MADV_KEEPONFORK = 19
    const MADV_COLD = 20
    const MADV_PAGEOUT = 21
    const MADV_HWPOISON = 100
    const MADV_SOFT_OFFLINE = 101
elseif Sys.isapple()
    const MADV_FREE = 5
elseif Sys.isfreebsd() || Sys.isdragonfly()
    const MADV_FREE = 5
    const MADV_NOSYNC = 6
    const MADV_AUTOSYNC = 7
    const MADV_NOCORE = 8
    const MADV_CORE = 9
    if Sys.isfreebsd()
        const MADV_PROTECT = 10
    else
        const MADV_INVAL = 10
        const MADV_SETMAP = 11
    end
elseif Sys.isopenbsd() || Sys.isnetbsd()
    const MADV_SPACEAVAIL = 5
    const MADV_FREE = 6
end

"""
    Mmap.madvise!(array, flag::Integer = Mmap.MADV_NORMAL)

Advises the kernel on the intended usage of the memory-mapped `array`, with the intent
`flag` being one of the available `MADV_*` constants.
"""
function madvise!(m::Array, flag::Integer=MADV_NORMAL)
    offset = rem(UInt(pointer(m)), PAGESIZE)
    ptr = pointer(m) - offset
    mmaplen = sizeof(m) + offset
    GC.@preserve m begin
        systemerror("madvise",
                    ccall(:madvise, Cint, (Ptr{Cvoid}, Csize_t, Cint), ptr, mmaplen, flag) != 0)
    end
end
madvise!(B::BitArray, flag::Integer=MADV_NORMAL) = madvise!(B.chunks, flag)
end # Sys.isunix()

end # module

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
#=

  "ESSE: Environmental and State Dependent Diversification" submodule package

=#

module ESSE

using Random: randexp
using DelimitedFiles: readdlm, writedlm
using ProgressMeter: Progress, next!
using DifferentialEquations: ODEProblem, init, reinit!, solve!, Tsit5, DiffEqBase
using LinearAlgebra: BLAS.gemv!, rank, mul!, ldiv!, qr!
using Distributed: @sync, @distributed
using SharedArrays: SharedArray

# other submodules dependencies
using ..Utils

# files
include("esse/musse_eqs.jl")
include("esse/esse_eqs.jl")
include("esse/states_handling.jl")
include("esse/loglik_mu.jl")
include("esse/simulate_esse.jl")
include("esse/approx_fun.jl")
include("esse/model_definitions.jl")
include("esse/simulate_geoesse.jl")
include("esse/sse_wrapper.jl")
include("esse/slice_sampler_utils.jl")
include("esse/slice_sampler.jl")
include("esse/ode_solve.jl")
include("esse/egeohisse_eqs.jl")
include("esse/prepare_data.jl")
include("esse/prepare_ll.jl")
include("esse/loglik_geo.jl")
include("esse/loglik_flow.jl")

end # module ESSE

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
using MAT
using Random
include("../include/FlowSeed.jl")
include("../src/HyperLocal.jl")

## Read in the matrix
s = time()
@time M = matread("../data/stackoverflow_answer_H.mat")
LabelMatrix = M["LabelMatrix"]
LabelNames = M["LabelNames"]
MainLabels = M["MainLabels"]
H = M["H"]
order = round.(Int64,vec(sum(H,dims=2)))
d = vec(sum(H,dims=1))
volA = sum(d)
m,n = size(H)
toc = time()-s
println("Done loading things into memory in $toc seconds.")

# Updated MainLabels be the set of topics between 2000-10000
#   (already stored in current MainLabels), to be only those with conductanc < .2
Tconds = Vector{Float64}()
Tsizes = Vector{Int64}()
labels = Vector{Int64}()
for lab = 1:length(MainLabels)
    label = MainLabels[lab]
    T = findnz(LabelMatrix[:,label])[1]
    nT = length(T)
    condT, volT, cutT = tl_cond(H,T,d,1.0,volA,order)
    if condT < .2
        # println("$nT \t $condT \t"*LabelNames[label])
        push!(Tconds,condT)
        push!(labels,label)
        push!(Tsizes,nT)
    end
end

##

using Plots
using Statistics
ntimes = 1
epsilon = 1.0
delta = 5000.0
outputmat = "Output_Stack/Set45_$(epsilon)_$(delta).mat"

mat = matread(outputmat)
hl_time = round.(mean(mat["hl_time"],dims = 2),digits = 2)
sl_time = round.(mean(mat["sl_time"],dims = 2),digits = 2)
zl_time = round.(mean(mat["zl_time"],dims = 2),digits = 2)
hl_size= round.(mean(mat["hl_size"],dims = 2),digits = 2)
r = round.(mean(mat["r_f1"],dims = 2),digits = 2)
hl = round.(mean(mat["hl_f1"],dims = 2),digits = 2)
b1 = round.(mean(mat["b1_f1"],dims = 2),digits = 2)
b2 = round.(mean(mat["b2_f1"],dims = 2),digits = 2)
nS = round.(mean(mat["newS"],dims = 2),digits = 2)
sl = round.(mean(mat["sl_f1"],dims = 2),digits = 2)
zl = round.(mean(mat["zl_f1"],dims = 2),digits = 2)

ln = LabelNames[labels]
plot()
p = (sortperm(vec(hl)))
y1 = hl[p]
yb1 = b1[p]
yr = r[p]
yb2 = b2[p]
y2 = max.(yb1,yb2)
y4 = max.(sl[p],zl[p])
x = 1:length(y1)
xlabels = ln[p]

## Wide plot
s1 = 750
s2 = 400
ms = 6
stepy = 1
scatter(x,y1,grid = false, markersize = ms, label = "HyperLocal",
markerstrokewidth = 0, markershape = :circle, linewidth = 0,
legend = :topleft, size = (s1, s2), ymirror = false, ylabel = "F1 Scores",
xticks = (1:stepy:length(ln), "\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t".*ln[1:stepy:length(ln)]),xrotation = 40,
xtickfont=font(8),
ytickfont=font(11),
guidefont=font(12),
titlefont=font(10),
legendfont=font(9)
)
scatter!(x,y2, markersize = ms,markerstrokewidth = 0,label = "TN/BN",
 markershape = :circle)
 scatter!(x,y4,markersize = ms,markerstrokewidth = 0, label = "FlowSeed",
  markershape = :circle)

savefig("Plots/StackDots_Wide.pdf")

max_stars_repo_name: nveldt/HypergraphFlowClustering

max_issues_repo_name: nveldt/HypergraphFlowClustering

max_forks_repo_name: nveldt/HypergraphFlowClustering


text:
# const getf_worker = Libdl.dlsym(polo_lib, :getf_worker)
# const getx_worker = Libdl.dlsym(polo_lib, :getx_worker)

struct Worker <: ParameterServer
    fval::Float64
    x::Vector{Float64}
    psopts::ParameterServerOptions

    function (::Type{Worker})(; kw...)
        worker =  new(0., Vector{Float64}(), ParameterServerOptions(; kw...))
    end
end

function initialize!(worker::Worker, x₀::AbstractVector)
    resize!(worker.x,length(x₀))
    worker.x .= x₀
end
getx(worker::Worker) = worker.x

paramserver_handle(::Worker) = POLO.proxgradient_worker
execution_handle(::Worker) = POLO.run_worker
delete_handle(::Worker) = POLO.delete_proxgradient_worker
# getf_handle(::Worker) = getf_worker
# getx_handle(::Worker) = getx_worker


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

# This is the library we care about
products = Product[
    LibraryProduct(prefix, "libpq", :LIBPQ_HANDLE),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/invenia/LibPQBuilder/releases/download/v10.3-1-4"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Windows(:i686)          => ("$bin_prefix/libpq.i686-w64-mingw32.tar.gz", "6acdd37081a3d80313263cac41e161ba6636e0d8a2f61cdae7447717f8f2aadb"),
    Windows(:x86_64)        => ("$bin_prefix/libpq.x86_64-w64-mingw32.tar.gz", "c3af66cf3f3b149dc1fcec4c4e4a6bb42d6514019e90be9b1b01e5e367725354"),
    Linux(:i686, :glibc)    => ("$bin_prefix/libpq.i686-linux-gnu.tar.gz", "89543222abea46936b9fe0426e4dab4d2fea8aff49823999100b1487dee42f6e"),
    Linux(:x86_64, :glibc)  => ("$bin_prefix/libpq.x86_64-linux-gnu.tar.gz", "f31d0f5f43acc458223b0f66711993106e98a75c7666d5839ae561981f22bb5c"),
    MacOS()                 => ("$bin_prefix/libpq.x86_64-apple-darwin14.tar.gz", "d2fa7faa26caa02790eac3743c864c3153998015a3ffdc68a229415cf2840045"),
)
# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=verbose) for p in products)
    if platform_key() in keys(download_info)
        # Download and install binaries
        url, tarball_hash = download_info[platform_key()]
        install(url, tarball_hash; prefix=prefix, force=true, verbose=true)
    else
        error("Your platform $(Sys.MACHINE) is not supported by this package!")
    end

    # Finally, write out a deps.jl file
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
end


text:
include("wrappers/wrappers.jl")
include("examples/examples.jl")
include("non_interactive/non_interactive.jl")
include("3rd_party/structs.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "gmsh"
version = v"4.10.2"

# Collection of sources required to build Gmsh
sources = [
    ArchiveSource("https://gmsh.info/src/gmsh-$(version)-source.tgz",
                  "d921109e6c419f68d5507d9612ee815d40b5239a5db175f28285e0cd0dbf2517"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gmsh-*
if [[ "${target}" == *linux* ]] || [[ "${target}" == *freebsd* ]]; then
    OPENGL_FLAGS="-DOpenGL_GL_PREFERENCE=LEGACY"
fi
mkdir build
cd build
cmake .. \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DENABLE_BUILD_DYNAMIC=1 \
    -DDEFAULT=1 \
    ${OPENGL_FLAGS}
make -j${nproc}
make install
mv ${prefix}/lib/gmsh.jl ${prefix}/lib/gmsh.jl.bak
sed ${prefix}/lib/gmsh.jl.bak \
  -e 's/^\(import Libdl\)/#\1/g' \
  -e 's/^\(const lib.*\)/#\1/g' \
  -e 's/^\(module gmsh\)$/\1\nusing gmsh_jll: libgmsh\nconst lib = libgmsh/g' \
  > ${prefix}/lib/gmsh.jl
install_license ../LICENSE.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct(["libgmsh", "gmsh"], :libgmsh),
    FileProduct("lib/gmsh.jl", :gmsh_api)
]

# Some dependencies are needed or available only on certain platforms
x11_platforms = filter(p ->Sys.islinux(p) || Sys.isfreebsd(p), platforms)
hdf5_platforms = [
    Platform("x86_64", "linux"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows")
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_xorgproto_jll"; platforms=x11_platforms),
    Dependency("Cairo_jll"),
    Dependency("CompilerSupportLibraries_jll"; platforms=filter(!Sys.isbsd, platforms)),
    Dependency("FLTK_jll"),
    Dependency("FreeType2_jll"),
    Dependency("GLU_jll"; platforms=x11_platforms),
    Dependency("GMP_jll"; compat="6.2"),
    Dependency("HDF5_jll"; platforms=hdf5_platforms, compat="1.12.1"),
    Dependency("JpegTurbo_jll"),
    Dependency("Libglvnd_jll"; platforms=x11_platforms),
    Dependency("libpng_jll"),
    Dependency("LLVMOpenMP_jll"; platforms=filter(Sys.isbsd, platforms)),
    Dependency("METIS_jll"),
    Dependency("MMG_jll"),
    Dependency("OCCT_jll"),
    Dependency("Xorg_libX11_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXext_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXfixes_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXft_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXinerama_jll"; platforms=x11_platforms),
    Dependency("Xorg_libXrender_jll"; platforms=x11_platforms),
    Dependency("Zlib_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"9")


text:
# # Model Forecasting

# ## Set the State-Space model
using Plots
using AnalogDataAssimilation
using DifferentialEquations

σ = 10.0
ρ = 28.0
β = 8.0 / 3

dt_integration = 0.01
dt_states = 1
dt_obs = 8
parameters = [σ, ρ, β]
var_obs = [1]
nb_loop_train = 100
nb_loop_test = 10
sigma2_catalog = 0.0
sigma2_obs = 2.0

ssm = StateSpaceModel(
    lorenz63,
    dt_integration,
    dt_states,
    dt_obs,
    parameters,
    var_obs,
    nb_loop_train,
    nb_loop_test,
    sigma2_catalog,
    sigma2_obs,
)

# - compute u0 to be in the attractor space

u0 = [8.0; 0.0; 30.0]
tspan = (0.0, 5.0)
prob = ODEProblem(ssm.model, u0, tspan, parameters)
u0 = last(solve(prob, save_everystep = false))

# ## Generate data

xt, yo, catalog = generate_data(ssm, u0);

plot(xt.t, xt[1])
scatter!(yo.t, yo[1]; markersize = 2)

# ## Data assimilation with model forecasting

np = 100
DA = DataAssimilation(ssm, xt)
@time x̂ = forecast(DA, yo, AnPF(np), progress = false);
println(RMSE(xt, x̂))

# ## Plot the times series

plot(xt.t, x̂[1])
scatter!(xt.t, xt[1]; markersize = 2)
plot!(xt.t, x̂[2])
scatter!(xt.t, xt[2]; markersize = 2)
plot!(xt.t, x̂[3])
scatter!(xt.t, xt[3]; markersize = 2)

# ## Plot the phase-space plot

p = plot3d(
    1,
    xlim = (-25, 25),
    ylim = (-25, 25),
    zlim = (0, 50),
    title = "Lorenz 63",
    marker = 2,
)
for x in eachrow(vcat(x̂.u'...))
    push!(p, x...)
end
p


text:
println("Julia is Alive")
time_julia_alive = time()

using Parsimonious

REPEAT = 20  # Repeat parse several times and take minimum as benchmark

# As a baseline for speed, parse some JSON.
#
# I have no reason to believe that JSON is a particularly representative or
# revealing grammar to test with. Also, this is a naive, unoptimized,
# incorrect grammar, so don't use it as a basis for comparison with other
# parsers. It's just meant to compare across versions of Parsimonious.

father =
"""{
        "id" : 1,
        "married" : true,
        "name" : "Larry Lopez",
        "sons" : null,
        "daughters" : [
          {
            "age" : 26,
            "name" : "Sandra"
            },
          {
            "age" : 25,
            "name" : "Margaret"
            },
          {
            "age" : 6,
            "name" : "Mary"
            }
          ]
        }
"""
more_fathers = join(repeat([father], inner=[60]), ",")

sample_json = """{"fathers" : [""" * more_fathers * "]}"

not_really_json_spec = """
    value = space (string / number / object / array / true_false_null)
            space

    object = "{" members "}"
    members = (pair ("," pair)*)?
    pair = string ":" value
    array = "[" elements "]"
    elements = (value ("," value)*)?
    true_false_null = "true" / "false" / "null"

    dq = '"'
    string = space dq chars dq space
    chars = ~'[^"]*'  # TODO implement the real thing
    number = (int frac exp) / (int exp) / (int frac) / int
    int = "-"? ((digit1to9 digits) / digit)
    frac = "." digits
    exp = e digits
    digits = digit+
    e = "e+" / "e-" / "e" / "E+" / "E-" / "E"

    digit1to9 = ~"[1-9]"
    digit = ~"[0-9]"
    space = ~'\\s*'
"""

grammar = Grammar(not_really_json_spec)

# Compile
parse(grammar, sample_json)

time_compiled = time()
elapsed_compiled = time_compiled - time_julia_alive

@printf "Importing, compiling and initializing grammar done after %.2f s\n\n" elapsed_compiled

# Disabling garbage collection (in Julia) during parse gets more consistent
# results.  The mean time doesn't seem to change much and is sometimes faster
# with gc enabled somehow. Anyway, let's not disable it so we can benchmark
# more real-world.

# gc_disable()
min_time = minimum([@elapsed parse(grammar, sample_json) for i in 1:REPEAT])
# gc_enable()
# gc()

# Show time and memory allocated for one run. At the moment, should be around
# 70 ms and 15 MB.

println("For one pass over the text: @time parse(grammar, sample_json")
@time parse(grammar, sample_json)
println()

kb = length(sample_json)/1024.
ms = min_time*1e3
kbps =(length(sample_json)/min_time/1024.)
@printf "Min time to parse %.1f KB is %.0f ms which is %.0f KB/s\n" kb ms kbps


text:
using ZXCalculus
using YaoPlots

function generate_example()
    zxd = ZXDiagram(4)
    push_gate!(zxd, Val{:Z}(), 1, 3//2)
    push_gate!(zxd, Val{:H}(), 1)
    push_gate!(zxd, Val{:Z}(), 1, 1//2)
    push_gate!(zxd, Val{:H}(), 4)
    push_gate!(zxd, Val{:CZ}(), 4, 1)
    push_gate!(zxd, Val{:CNOT}(), 1, 4)
    push_gate!(zxd, Val{:H}(), 1)
    push_gate!(zxd, Val{:H}(), 4)
    push_gate!(zxd, Val{:Z}(), 1, 1//4)
    push_gate!(zxd, Val{:Z}(), 4, 3//2)
    push_gate!(zxd, Val{:X}(), 4, 1//1)
    push_gate!(zxd, Val{:H}(), 1)
    push_gate!(zxd, Val{:Z}(), 4, 1//2)
    push_gate!(zxd, Val{:X}(), 4, 1//1)

    push_gate!(zxd, Val{:Z}(), 2, 1//2)
    push_gate!(zxd, Val{:CNOT}(), 3, 2)
    push_gate!(zxd, Val{:H}(), 2)
    push_gate!(zxd, Val{:CNOT}(), 3, 2)
    push_gate!(zxd, Val{:Z}(), 2, 1//4)
    push_gate!(zxd, Val{:Z}(), 3, 1//2)
    push_gate!(zxd, Val{:H}(), 2)
    push_gate!(zxd, Val{:H}(), 3)
    push_gate!(zxd, Val{:Z}(), 3, 1//2)
    push_gate!(zxd, Val{:CNOT}(), 3, 2)

    return zxd
end

zxd = generate_example()
ex_zxd = clifford_simplification(zxd)

plot(zxd)
plot(ex_zxd)


text:
using Tries
function Either(x::Vector{<:AbstractString})
    P = Trie{Char,Union{Missing,Nothing}}
    r = P()
    for e in x
        r[e...] = nothing
    end
    Either{NCodeunitsState, String}(r)
end

function Either(x::Dict)
    P = Trie{Char,Union{Missing,valtype(x)}}
    r=P()
    for (e,v) in pairs(x)
        r[e...] = v
    end
    Either{NCodeunitsState, valtype(x)}(r)
end
either_state_type(T::Type{<:Trie}) = NCodeunitsState{T}

"""
    _iterate(p::Trie{Char}, str, till, posi, next_i, ::Nothing)

Match char path in `p` greedily, recording shorter matches in state.
"""
@inline function _iterate(p::Union{Trie{Char},SubTrie{Char}}, str, till, posi, next_i, state::Nothing)
    ni = ni_ = posi
    st = st_ = p
    while st !== nothing && ni <= till
        @inbounds c, ni = iterate(str,ni)
        @inbounds st = get( Tries.nodes(st),c, nothing)
        st === nothing && break
        if get(st) !== missing
            ni_ = ni
            st_ = st
        end
    end
    if get(st_) !== missing
        return ni_, NCodeunitsState(ni_-posi,st_)
    else
        return nothing
    end
end

@inline function _iterate(p::AbstractTrie{Char}, str, till, posi, next_i, state)
    _iterate(p, str, _prevind(str,next_i,2), posi, posi, nothing)
end

function _iterate(p::Either{<:AbstractTrie}, str, till, posi, next_i, state)
    _iterate(p.options, str, till, posi, next_i, state)
end

@inline _rightof(str,i,parser::Either{<:AbstractTrie},x::NCodeunitsState)  = i+x.nc
@inline _leftof(str,i,parser::Either{<:AbstractTrie},x::NCodeunitsState)  = i-x.nc

children(x::Either{<:AbstractTrie}) =
    children(x.options)

function deepmap_parser(f::typeof(lowercase),mem::AbstractDict,x::Either{<:AbstractTrie},a...;kw...)
    get!(mem,x) do
        g = (lowercase.(Tries.path(st))=>get(st)
             for st in PreOrderDFS(x.options)
             if !isempty(Tries.path(st)))
        Either{result_type(x)}(Trie(g))
    end
end

# struct InternedParser{E, I<:Either{<:AbstractTrie}} <: WrappedParser{E,NCodeunitsState,T} where {P}
#     interned::I
#     parser::E
# end


text:
@doc raw"""
    LinearMixingModelKernel(k::Kernel, H::AbstractMatrix)
    LinearMixingModelKernel(Tk::AbstractVector{<:Kernel},Th::AbstractMatrix)

Kernel associated with the linear mixing model, taking a vector of `m` kernels and a `m × p` matrix H for a function with `p` outputs. Also accepts a single kernel `k` for use across all `m` basis vectors.

# Definition

For inputs ``x, x'`` and output dimensions ``p_x, p_{x'}'``, the kernel is defined as[^BPTHST]
```math
k\big((x, p_x), (x, p_{x'})\big) = H_{:,p_{x}}K(x, x')H_{:,p_{x'}}
```
where ``K(x, x') = Diag(k_1(x, x'), ..., k_m(x, x'))`` with zero off-diagonal entries.
``H_{:,p_{x}}`` is the ``p_x``-th column (`p_x`-th output) of ``H \in \mathbb{R}^{m \times p}``
representing ``m`` basis vectors for the ``p`` dimensional output space of ``f``.
``k_1, \ldots, k_m`` are ``m`` kernels, one for each latent process, ``H`` is a
mixing matrix of ``m`` basis vectors spanning the output space.

[^BPTHST]: Wessel P. Bruinsma, Eric Perim, Will Tebbutt, J. Scott Hosking, Arno Solin, Richard E. Turner (2020). [Scalable Exact Inference in Multi-Output Gaussian Processes](https://arxiv.org/pdf/1911.06287.pdf).
"""
struct LinearMixingModelKernel{Tk<:AbstractVector{<:Kernel},Th<:AbstractMatrix} <: MOKernel
    K::Tk
    H::Th
    function LinearMixingModelKernel(Tk::AbstractVector{<:Kernel}, H::AbstractMatrix)
        @assert length(Tk) == size(H, 1) "Number of kernels and number of rows in H must match"
        return new{typeof(Tk),typeof(H)}(Tk, H)
    end
end

function LinearMixingModelKernel(k::Kernel, H::AbstractMatrix)
    return LinearMixingModelKernel(Fill(k, size(H, 1)), H)
end

function (κ::LinearMixingModelKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,Int})
    (px > size(κ.H, 2) || py > size(κ.H, 2) || px < 1 || py < 1) &&
        error("`px` and `py` must be within the range of the number of outputs")
    return sum(κ.H[i, px] * κ.K[i](x, y) * κ.H[i, py] for i in 1:length(κ.K))
end

function Base.show(io::IO, k::LinearMixingModelKernel)
    return print(io, "Linear Mixing Model Multi-Output Kernel")
end

function Base.show(io::IO, mime::MIME"text/plain", k::LinearMixingModelKernel)
    print(io, "Linear Mixing Model Multi-Output Kernel. Kernels:")
    for k in k.K
        print(io, "\n\t")
        show(io, mime, k)
    end
end


text:
using Documenter, SystemSnoop

makedocs(
    modules = [SystemSnoop],
    format = :html,
    sitename = "SystemSnoop.jl",
    html_prettyurls = get(ENV, "CI", nothing) == "true",
    pages = Any[
        "index.md",
        "trace.md",
        "process.md",
        "Measurements" => [
            "Idle Page Tracking" => [
                "measurements/idlepages/idlepages.md",
                "measurements/idlepages/sample.md",
                "measurements/idlepages/vma.md",
                "measurements/idlepages/rangevector.md",
                "measurements/idlepages/utils.md",
            ],
            "measurements/diskio.md",
            "measurements/processio.md",
            "measurements/smaps.md",
            "measurements/statm.md",
            "measurements/uptime.md",
        ],
        "utils.md",
        "hugepages.md",
        "proof-of-concept.md",
        "thoughts.md",
        "docstring_index.md",
    ]
)

deploydocs(
    repo = "github.com/hildebrandmw/SystemSnoop.jl.git",
    target = "build",
    julia = "1.1",
    deps = nothing,
    make = nothing,
)


text:
using PortAudio
using DSP

function create_measure_signal()
    signal = zeros(Float32, 20000)
    for i in 1:3
        signal = vcat(signal, rand(Float32, 100), zeros(Float32, i*10000))
    end
    return signal
end

function measure_latency(in_latency = 0.1, out_latency=0.1; is_warmup = false)

    in_stream = PortAudioStream(1,0; latency=in_latency)
    out_stream = PortAudioStream(0,1; latency=out_latency)

    cond = Base.Event()

    writer_start_time = Int64(0)
    reader_start_time = Int64(0)

    reader = Threads.@spawn begin
        wait(cond)
        writer_start_time = time_ns() |> Int64
        return read(in_stream, 100000)
    end

    signal = create_measure_signal()
    writer = Threads.@spawn begin
            wait(cond)
            reader_start_time = time_ns() |> Int64
            write(out_stream, signal)
    end

    notify(cond)

    wait(reader)
    wait(writer)

    recorded = collect(reader.result)[:,1]

    close(in_stream)
    close(out_stream)

    diff = reader_start_time - writer_start_time |> abs

    diff_in_ms = diff / 10^6 # 1 ms = 10^6 ns

    if !is_warmup && diff_in_ms > 1
        @warn "Threads start time difference $diff_in_ms ms is bigger than 1 ms"
    end

    delay = finddelay(recorded, signal) / 48000

    return trunc(Int, delay * 1000)# result in ms
end

measure_latency(0.1, 0.1; is_warmup = true) # warmup

latencies = [0.1, 0.01, 0.005]
for in_latency in latencies
    for out_latency in latencies
        measure = measure_latency(in_latency, out_latency)
        println("$measure ms latency for in_latency=$in_latency, out_latency=$out_latency")
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 63577d5e-8d41-4bc5-8d4c-c74be034d08c
using Pkg, DrWatson

# ╔═╡ a2d386d3-782f-4042-bae9-9751c73d42b5
begin
    # Use '#@quickactivate ...' to enable Pluto Pkg management.
    # Remove '#' to disable Pluto Pkg management.
    # Careful: It will _no longer_ be a reproducible environment!
    @quickactivate "SR2StanPluto"
end

# ╔═╡ 88faad64-f887-44d6-99e5-d56e9d6df486
begin
	using DataFrames
    using PrettyTables
    using StanSample
    using StatisticalRethinking
    using StatisticalRethinkingPlots
    using PlutoUI
end

# ╔═╡ 865295d7-451e-437b-bc4d-cbd939d3a080
md"#### Below cell enables switching between dev-environments and package environments. In most cases use `#@quickactivate ...`"

# ╔═╡ da61fb27-229d-4885-ba23-580f088db42b
begin
    df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);
    scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
    data = (N=size(df, 1), D=df.Divorce_s, A=df.MedianAgeMarriage_s,
        M=df.Marriage_s)
end;

# ╔═╡ d4cc546f-cb4d-401d-bf7e-5c87e1ac8fcd
df

# ╔═╡ c6ccfa63-fd10-4d2a-b260-8fdba221f6d7
stan5_1 = "
data {
    int < lower = 1 > N; // Sample size
    vector[N] D; // Outcome
    vector[N] A; // Predictor
}
parameters {
    real a; // Intercept
    real bA; // Slope (regression coefficients)
    real < lower = 0 > sigma;    // Error SD
}
transformed parameters {
    vector[N] mu;               // mu is a vector
    for (i in 1:N)
        mu[i] = a + bA * A[i];
}
model {
    a ~ normal(0, 0.2);         //Priors
    bA ~ normal(0, 0.5);
    sigma ~ exponential(1);
    D ~ normal(mu , sigma);     // Likelihood
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(D[i] | mu[i], sigma);
}
";

# ╔═╡ 5041bb2f-ba75-47cc-9a9a-40b6882fe161
begin
    m5_1s = SampleModel("m5.1s", stan5_1)
    rc5_1s = stan_sample(m5_1s; data)
end;

# ╔═╡ c63617a8-aea2-4e01-baca-e000fe1a8268
if success(rc5_1s)
    st5_1s = read_samples(m5_1s, :table)
end

# ╔═╡ 6d3a8e78-0758-463c-aa18-948dee2decef
df5_1s = DataFrame(st5_1s)

# ╔═╡ 4e8dfe63-f39b-4847-8e62-ad184ae222c6
md"#### Select a set of columns."

# ╔═╡ 270a906b-575c-4bcf-adcd-1ddfca3b42a5
df5_1s[:, [:a, :bA]]

# ╔═╡ be3252ef-32a6-411b-8c69-fc0a41b8ca96
md"#### Select a block of variables using Stan '.' notation."

# ╔═╡ dfa87068-887f-4dba-95f7-b4c472de3e37
DataFrame(df5_1s, :log_lik)

# ╔═╡ 984d07a4-e74a-4a56-bae9-8d833aa8f78e
md"#### Do not append multiple chains."

# ╔═╡ 3e26a611-9a74-407a-9ebe-927e9c94d0d7
if success(rc5_1s)
    sts5_1s = read_samples(m5_1s, :tables)
end

# ╔═╡ 149be65d-c4cf-4a7a-a250-a91bd72c6417
md"#### Create a DataFrame for each chain ( Vector{DataFrame} )."

# ╔═╡ 516b91db-7ebf-4699-95c9-1863d5d2b5b6
dfs5_1s = DataFrame.(sts5_1s)

# ╔═╡ 307306b3-e28d-440f-8023-14cb1fd62185
md" #### Extract Stan block of variables for each chain."

# ╔═╡ 0bc25437-a231-47c6-b92e-e4a09c94c962
DataFrame.(dfs5_1s, :log_lik)

# ╔═╡ Cell order:
# ╠═63577d5e-8d41-4bc5-8d4c-c74be034d08c
# ╟─865295d7-451e-437b-bc4d-cbd939d3a080
# ╠═a2d386d3-782f-4042-bae9-9751c73d42b5
# ╠═88faad64-f887-44d6-99e5-d56e9d6df486
# ╠═da61fb27-229d-4885-ba23-580f088db42b
# ╠═d4cc546f-cb4d-401d-bf7e-5c87e1ac8fcd
# ╠═c6ccfa63-fd10-4d2a-b260-8fdba221f6d7
# ╠═5041bb2f-ba75-47cc-9a9a-40b6882fe161
# ╠═c63617a8-aea2-4e01-baca-e000fe1a8268
# ╠═6d3a8e78-0758-463c-aa18-948dee2decef
# ╟─4e8dfe63-f39b-4847-8e62-ad184ae222c6
# ╠═270a906b-575c-4bcf-adcd-1ddfca3b42a5
# ╟─be3252ef-32a6-411b-8c69-fc0a41b8ca96
# ╠═dfa87068-887f-4dba-95f7-b4c472de3e37
# ╟─984d07a4-e74a-4a56-bae9-8d833aa8f78e
# ╠═3e26a611-9a74-407a-9ebe-927e9c94d0d7
# ╟─149be65d-c4cf-4a7a-a250-a91bd72c6417
# ╠═516b91db-7ebf-4699-95c9-1863d5d2b5b6
# ╟─307306b3-e28d-440f-8023-14cb1fd62185
# ╠═0bc25437-a231-47c6-b92e-e4a09c94c962


text:
# common.jl
# Common provides types and helper functions
#
# Copyright (C) 2012   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

## COMMON FUNCTIONS

# Tied ranking, also computing the adjustment for ties
# sum(n^3-n) where n is the number of ties for each data point
function tiedrank_adj!(ord::AbstractVector, v::AbstractArray)
    n     = length(v)
    place = sortperm(v)
    tieadj = 0.0

    i = 1
    while i <= n
        j = i
        while j + 1 <= n && v[place[i]] == v[place[j + 1]]
            j += 1
        end

        if j > i
            t = j - i + 1
            m = sum(i:j) / t
            tieadj += t^3 - t
            for k = i:j
                ord[place[k]] = m
            end
        else
            ord[place[i]] = i
        end

        i = j + 1
    end

    (ord, tieadj)
end

tiedrank_adj(v::AbstractArray) = tiedrank_adj!(Vector{Float64}(undef, length(v)), v)

# Pool covariance matrices, overwriting the first
function poolcov!(Sx::AbstractMatrix, nxm1::Int, Sy::AbstractMatrix, nym1::Int)
    @inbounds for i = eachindex(Sx, Sy)
        Sx[i] = (Sx[i] * nxm1 + Sy[i] * nym1) / (nxm1 + nym1)
    end
    Sx
end


text:
include("C:/Users/stravaglino3/Downloads/julia_code/julia_code/model_fit_global.jl")
include("C:/Users/stravaglino3/Downloads/julia_code/julia_code/models.jl")
include("C:/Users/stravaglino3/Downloads/julia_code/julia_code/helper_functions.jl")
using DataFrames
tcr_mhc_j = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/dissociation data/jurkat_TCR_MHC.csv", DataFrame)
cd4_mhc_j = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/dissociation data/jurkat_CD4_MHC.csv", DataFrame)
tcr_mhc_cd4_j = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/dissociation data/jurkat_TCR_MHC_CD4.csv", DataFrame)

tcr_mhc = clean_dissociation_data(tcr_mhc_j)
cd4_mhc = clean_dissociation_data(cd4_mhc_j)
tcr_mhc_cd4 = clean_dissociation_data(tcr_mhc_cd4_j)

# pyplot()
plot()
plot_dissociations(tcr_mhc)
plot_dissociations(cd4_mhc)
plot_dissociations(tcr_mhc_cd4)

t_tm, lnp_tm = separate_diss_data(tcr_mhc)
t_cm, lnp_cm = separate_diss_data(cd4_mhc)
t_tmc, lnp_tmc = separate_diss_data(tcr_mhc_cd4)
p_tm = lnp2p.(lnp_tm)
p_cm = lnp2p.(lnp_cm)
p_tmc = lnp2p.(lnp_tmc)

plot()
fs_tm = extract_avg_forces(tcr_mhc)
fs_cm = extract_avg_forces(cd4_mhc)
fs_tmc = extract_avg_forces(tcr_mhc_cd4)

scatter(t_tm,lnp_tm)

## TCR-MHC
u₀ = [1.0, 0.]          # Initial condition
tspan = (0.0,10.0)  # Simulation interval
model_opt = (-5.,3.)
model_opt = (0.,100.)
pinit =[11.,0.01,0.14,0.4,3,0.2,20,4]
# pinit = log10.(pinit)
dens_j = fill([],length(t_tm))
model = bimolecular_diss_bell!
# pinit = 10 .^ [0.8934057970002582, 0.6126600385641812, -2.0499405365827617, 0.7136091557046687, -0.6215834038143795, 0.23515760918441408, -0.11706056162578093, -1.5512600092183506]
rates_tm,hes_tm,l_tm,boot_tm = optimization(model,p_tm,
                    dens_j,fs_tm,model_opt,t_tm)

## CD4-MHC

rates_cm,hes_cm,l_cm = optimization(model,p_cm,
                            dens_j,fs_cm,model_opt,t_cm)

## TCR-MHC-CD4

u₀ = [0.8, 0., 0.2]
pinit = [1,1.e-2,1,1.e-2]
bimolecular_cons = 10 .^ rates_tm
model = trimolecular_diss_bell!
rates_tmc,hes_tmc,l_tmc = optimization(model,p_tmc,
                            dens_j,fs_tmc,model_opt,t_tmc)

p = 10 .^ rates
p = rates
p = [1.1385, -1.6272, -2.4367, -0.127488, 0.057953, -3.70439, 0.440775,-0.398936]
# cd4-TCR
p = rates
p = 10 .^ [0.41751687626687767, -0.23866855433293688, 0.06689410353314687, -0.5510289313898404, -0.09606981254493874, -0.1155499030032302, -3.9999999999999964, -1.7908513448187287]
 #tcr Mhc
p = rates_tm
p = 10 .^ rates_tm
p = 10 .^ [1.1854325378992236, -5.237609441955307, -0.2921252076084758, -2.703452812162055, 0.05179350352175943, -3.0727673775338427, -0.9157012289343508, -2.3711064932852954]

## ploting
dens = []

model = bimolecular_diss_bell!

gr()
p = rates_tm
p = 10 .^ rates_tm
p = [11.0,0.1,
        0.14,0.4,
        5.0,0.2,
        20.0,4.0]
show_diss_fit(model,tcr_mhc,fs_tm)
p1 = plot()
show_diss_fit!(p1,M,rates,data;cmap=nothing)
p = rates_cm
p = 10 .^ rates_cm
show_diss_fit(model,cd4_mhc,fs_cm)
p = rates_tmc
p = 10 .^ rates_tmc
p = [5.,1.,1.,0.1]
u₀ = [0.7, 0., 0.3]
model = trimolecular_diss_bell!
show_diss_fit(model,tcr_mhc_cd4,fs_tmc)

plotly()
function show_diss_fit(model, data, forces)
  plot()
  plot_diss_fit(model,forces)
  plot_dissociations(data)
end
function show_diss_fit(model, t, lnp, forces)
  plot()
  plot_diss_fit(model,forces)
  scatter!(t, lnp)
end

Fs = collect(0.:0.5:30.)
plot()
plot_force_lifetime_fit(model,Fs,tcr_mhc)
plot_force_lifetime_fit(model,Fs,cd4_mhc)

plot_F_vs_t(tcr_mhc)
plot_F_vs_t(cd4_mhc)

plot()
# p =[11.,0.,0.14,0.4,3,0.2,20,4]
plot_bell_equations(p)

max_t_tm  = extract_max_lifetimes(tcr_mhc)
scatter(fs_tm,max_t_tm)
max_lifetime_bounds(fs_tm,max_t_tm,0:0.1:30)


text:
using DifferentiableStateSpaceModels,
    SparseArrays,
    LinearAlgebra,
    Parameters,
    Test,
    TimerOutputs,
    BenchmarkTools

@testset "Dense SW First Order" begin
    m = @include_example_module(Examples.SW07)

    p = [
        10,
        0.51,
        10,
        0,
        0.7,
        0.742,
        0,
        0,
        0.24,
        0.2696,
        6.0144,
        0.025,
        1.5,
        0.6361,
        1.5,
        0.3243,
        0.8087,
        0.47,
        0.6,
        1.9423,
        1.5,
        1.488,
        0.2347,
        0.0593,
        0.8762,
        0.9977,
        0.5799,
        0.9957,
        0.7165,
        0,
        0,
        0,
        0.3982,
        0.18,
    ]
    generate_perturbation(m, p)
    reset_timer!()
    sol = generate_perturbation(m, p)
    print_timer()
    @test sol.retcode == :Success
end