text:
"""
    struct SparseMatrixCSR{Bi,Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}

Matrix type for storing Bi-based sparse matrices
in the Compressed Sparse Row format. The standard
way of constructing SparseMatrixCSR is through the
[`sparsecsr`](@ref) function.
"""
struct SparseMatrixCSR{Bi,Tv,Ti} <: AbstractSparseMatrix{Tv,Ti}
    offset :: Int             # Bi-index offset (Bi-1)
    m      :: Int             # Number of columns
    n      :: Int             # Number of rows
    rowptr :: Vector{Ti}      # Row i is in (rowptr[i]+offset):(rowptr[i+1]-Bi)
    colval :: Vector{Ti}      # Col indices (with base Bi) of stored values
    nzval  :: Vector{Tv}      # Stored values, typically nonzeros

    function SparseMatrixCSR{Bi}(m::Integer, n::Integer, rowptr::Vector{Ti}, colval::Vector{Ti},
                                    nzval::Vector{Tv}) where {Bi,Tv,Ti<:Integer}
        @noinline throwsz(str, lbl, k) =
            throw(ArgumentError("number of $str ($lbl) must be ≥ 0, got $k"))
        m < 0 && throwsz("rows", 'm', m)
        n < 0 && throwsz("columns", 'n', n)
        rowptr[m+1] == length(nzval) && throwsz("nnz", "rowptr[m+1]", length(nzval))
        offset = Bi-1
        rowptr .+= offset
        colval .+= offset
        new{Bi,Tv,Ti}(Int(offset), Int(m), Int(n), rowptr, colval, nzval)
    end
end

SparseMatrixCSR(transpose::SparseMatrixCSC) where {Tv,Ti} =
        SparseMatrixCSR{1}(transpose.n, transpose.m, transpose.colptr, transpose.rowval, transpose.nzval)
SparseMatrixCSR{Bi}(transpose::SparseMatrixCSC{Tv,Ti}) where {Bi,Tv,Ti} =
        SparseMatrixCSR{Bi}(transpose.n, transpose.m, transpose.colptr, transpose.rowval, transpose.nzval)
size(S::SparseMatrixCSR) = (S.m, S.n)

function show(io::IO, ::MIME"text/plain", S::SparseMatrixCSR)
    xnnz = nnz(S)
    print(io, S.m, "×", S.n, " ", typeof(S), " with ", xnnz, " stored ",
              xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        print(io, ":")
        show(IOContext(io, :typeinfo => eltype(S)), S)
    end
end
show(io::IO, S::SparseMatrixCSR) = show(convert(IOContext, io), S::SparseMatrixCSR)
#show(io::IOContext, S::SparseMatrixCSR) = dump(S)

function show(io::IOContext, S::SparseMatrixCSR{Bi}) where{Bi}
    nnz(S) == 0 && return show(io, MIME("text/plain"), S)

    ioc = IOContext(io, :compact => true)
    function _format_line(r, col, padr, padc)
        print(ioc, "\n  [", rpad(col, padr), ", ", lpad(S.colval[r], padc), "]  =  ")
        if isassigned(S.nzval, Int(r))
            show(ioc, S.nzval[r])
        else
            print(ioc, Base.undef_ref_str)
        end
    end

    function _get_cols(from, to)
        idx = eltype(S.rowptr)[]
        c = searchsortedlast(S.rowptr, from)
        for i = from:to
            while i == S.rowptr[c+1]
                c +=1
            end
            push!(idx, c)
        end
        idx
    end

    rows = displaysize(io)[1] - 4 # -4 from [Prompt, header, newline after elements, new prompt]
    if !get(io, :limit, false) || rows >= nnz(S) # Will the whole matrix fit when printed?
        cols = _get_cols(Bi, nnz(S)+S.offset)
        padr, padc = ndigits.((maximum(S.colval[1:nnz(S)]), cols[end]))
        _format_line.(1:nnz(S), cols.+S.offset, padr, padc)
    else
        if rows <= 2
            print(io, "\n  \u22ee")
            return
        end
        s1, e1 = 1, div(rows - 1, 2) # -1 accounts for \vdots
        s2, e2 = nnz(S) - (rows - 1 - e1) + 1, nnz(S)
        cols1, cols2 = _get_cols(s1+S.offset, e1+S.offset), _get_cols(s2+S.offset, e2+S.offset)
        padr = ndigits(max(maximum(S.colval[s1:e1]), maximum(S.colval[s2:e2])))
        padc = ndigits(cols2[end])
        _format_line.(s1:e1, cols1, padr, padc)
        print(io, "\n  \u22ee")
        _format_line.(s2:e2, cols2, padr, padc)
    end
    return
end

function getindex(S::SparseMatrixCSR{Bi,Tv,Ti}, i0::Integer, i1::Integer) where {Bi,Tv,Ti}
    if !(1 <= i0 <= S.m && 1 <= i1 <= S.n); throw(BoundsError()); end
    r1 = Int(S.rowptr[i0]-S.offset)
    r2 = Int(S.rowptr[i0+1]-Bi)
    (r1 > r2) && return zero(Tv)
    i1_= i1+S.offset
    r1 = searchsortedfirst(S.colval, i1_, r1, r2, Base.Order.Forward)
    ((r1 > r2) || (S.colval[r1] != i1_)) ? zero(Tv) : S.nzval[r1]
end

"""
    nnz(S::SparseMatrixCSR{Bi}) where {Bi}

Returns the number of stored (filled) elements in a sparse array.
"""
nnz(S::SparseMatrixCSR{Bi}) where {Bi} = length(nonzeros(S))

"""
    count(pred, S::SparseMatrixCSR) -> Integer

Count the number of elements in nonzeros(S) for which predicate pred returns true.
"""
count(pred, S::SparseMatrixCSR) = count(pred, nonzeros(S))

"""
    nonzeros(S::SparseMatrixCSR)

Return a vector of the structural nonzero values in sparse array S.
This includes zeros that are explicitly stored in the sparse array.
The returned vector points directly to the internal nonzero storage of S,
and any modifications to the returned vector will mutate S as well.
"""
nonzeros(S::SparseMatrixCSR) = S.nzval

"""
    nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi}

Return the range of indices to the structural nonzero values of a
sparse matrix row.
"""
nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi} = S.rowptr[row]-S.offset:S.rowptr[row+1]-Bi

"""
    function findnz(S::SparseMatrixCSR{Bi,Tv,Ti})

Return a tuple (I, J, V) where I and J are the row and column indices
of the stored ("structurally non-zero") values in sparse matrix A,
and V is a vector of the values.
"""
function findnz(S::SparseMatrixCSR{Bi,Tv,Ti}) where {Bi,Tv,Ti}
    numnz = nnz(S)
    I = Vector{Ti}(undef, numnz)
    J = Vector{Ti}(undef, numnz)
    V = Vector{Tv}(undef, numnz)

    count = 1
    @inbounds for row in 1:S.m
        @inbounds for k in nzrange(S,row)
            I[count] = S.colval[k]-S.offset
            J[count] = row
            V[count] = S.nzval[k]
            count += 1
        end
    end

    return (I, J, V)
end

"""
    rowvals(S::SparseMatrixCSR)

Return an error.
CSR sparse matrices does not contain raw row values.
It contains row pointers instead that can be accessed
by using [`nzrange`](@ref).
"""
rowvals(S::SparseMatrixCSR) = error("CSR sparse matrix does not contain raw row values")

"""
    colvals(S::SparseMatrixCSR)

Return a vector of the col indices of S.
Any modifications to the returned vector will mutate S as well.
Providing access to how the co,l indices are stored internally
can be useful in conjunction with iterating over structural
nonzero values. See also [`nonzeros`](@ref) and [`nzrange`](@ref).
"""

colvals(S::SparseMatrixCSR) = S.colval

"""
    sparsecsr(I, J, V, [m, n, combine])

Create a sparse matrix S of dimensions m x n such that S[I[k], J[k]] = V[k].
The combine function is used to combine duplicates.
If m and n are not specified, they are set to
maximum(I) and maximum(J) respectively.
If the combine function is not supplied, combine defaults to +
unless the elements of V are Booleans in which case combine defaults to |.
All elements of I must satisfy 1 <= I[k] <= m,
and all elements of J must satisfy 1 <= J[k] <= n.
Numerical zeros in (I, J, V) are retained as structural nonzeros.
"""
sparsecsr(I,J,args...) =
        SparseMatrixCSR(sparse(J,I,args...))
sparsecsr(::Type{SparseMatrixCSR},I,J,args...) =
        sparsecsr(I,J,args...)
sparsecsr(::Type{SparseMatrixCSR{Bi}},I,J,args...) where {Bi} =
        SparseMatrixCSR{Bi}(sparse(J,I,args...))
sparsecsr(::Type{SparseMatrixCSR{Bi}},I,J,V,m,n,args...) where {Bi} =
        SparseMatrixCSR{Bi}(sparse(J,I,V,n,m,args...))
sparsecsr(::Type{SparseMatrixCSR{Bi,Tv,Ti}},I::Vector{Ti},J::Vector{Ti},args...) where {Bi,Tv,Ti} =
        SparseMatrixCSR{Bi}(sparse(J,I,args...))
sparsecsr(::Type{SparseMatrixCSR{Bi,Tv,Ti}},I::Vector{Ti},J::Vector{Ti},V::Vector{Tv},m,n,args...) where {Bi,Tv,Ti} =
        SparseMatrixCSR{Bi}(sparse(J,I,V,n,m,args...))

"""
    function push_coo!(::Type{SparseMatrixCSR},I,J,V,ik,jk,vk)

Inserts entries in COO vectors for further building a SparseMatrixCSR.
"""
function push_coo!(::Type{SparseMatrixCSR},
    I::Vector,J::Vector,V::Vector,ik::Integer,jk::Integer,vk::Number) where {Bi}
    (push!(I, ik), push!(J, jk), push!(V, vk))
end

push_coo!(::Type{SparseMatrixCSR{Bi}}, I, J, V, ik, jk, vk) where {Bi} =
    push_coo!(SparseMatrixCSR, I, J, V, ik, jk, vk)

push_coo!(::Type{SparseMatrixCSR{Bi,Tv,Ti}}, I, J, V, ik, jk, vk) where {Bi,Tv,Ti} =
    push_coo!(SparseMatrixCSR, I, J, V, ik, jk, vk)

"""
    function finalize_coo!(::Type{SparseMatrixCSR},I,J,V,m,n)

Finalize COO arrays for building a SparseMatrixCSR.
"""
function finalize_coo!(::Type{SparseMatrixCSR},
    I::Vector,J::Vector,V::Vector,m::Integer,n::Integer)
end

finalize_coo!(T::Type{SparseMatrixCSR{Bi}}, I , J, V, m, n)  where {Bi} =
    finalize_coo!(SparseMatrixCSR, I, J, V, m, n)

finalize_coo!(T::Type{SparseMatrixCSR{Bi,Tv,Ti}}, I , J, V, m, n)  where {Bi,Tv,Ti} =
    finalize_coo!(SparseMatrixCSR, I, J, V, m, n)

"""
    function mul!(y::AbstractVector,A::SparseMatrixCSR,v::AbstractVector{T}) where {T}

Calculates the matrix-vector product ``Av`` and stores the result in `y`,
overwriting the existing value of `y`.
"""
function mul!(y::AbstractVector,A::SparseMatrixCSR,v::AbstractVector{T}) where {T}
    A.n == size(v, 1) || throw(DimensionMismatch())
    A.m == size(y, 1) || throw(DimensionMismatch())

    y .= zero(T)
    for row = 1:size(y, 1)
        @inbounds for nz in nzrange(A,row)
            col = A.colval[nz]-A.offset
            y[row] += A.nzval[nz]*v[col]
        end
    end
    return y
end

*(A::SparseMatrixCSR, v::Vector) = (y = similar(v,A.n);mul!(y,A,v))

"""
    function hasrowmajororder(::Type{SparseMatrixCSR})

Check if values are stored in row-major order.
Return true.
"""
hasrowmajororder(::Type{SparseMatrixCSR}) = true
hasrowmajororder(a::SparseMatrixCSR) = hasrowmajororder(SparseMatrixCSR)

"""
    function hascolmajororder(::Type{SparseMatrixCSR})

Check if values are stored in col-major order.
Return false.
"""
hascolmajororder(::Type{SparseMatrixCSR}) = false
hascolmajororder(a::SparseMatrixCSR) = hascolmajororder(SparseMatrixCSR)

"""
    function getptr(S::SparseMatrixCSR)

Return rows pointer.
"""
getptr(S::SparseMatrixCSR) = S.rowptr

"""
    function getindices(S::SparseMatrixCSR)

Return column indices.
"""
getindices(S::SparseMatrixCSR) = colvals(S)

"""
    function convert(::Type{SparseMatrixCSR}, x::AbstractSparseMatrix)

Convert x to a value of type SymSparseMatrixCSR.
"""
convert(::Type{SparseMatrixCSR}, x::AbstractSparseMatrix)  = convert(SparseMatrixCSR{1}, x)

function convert(::Type{SparseMatrixCSR{Bi}}, x::SparseMatrixCSR{xBi}) where {Bi,xBi}
    if Bi == xBi
        return x
    else
        return SparseMatrixCSR{Bi}( x.m,
                                    x.n,
                                    copy(getptr(x)).-x.offset,
                                    copy(getindices(x)).-x.offset,
                                    copy(nonzeros(x)))
    end
end

function convert(::Type{SparseMatrixCSR{Bi,Tv,Ti}}, x::SparseMatrixCSR{xBi,xTv,xTi}) where {Bi,Tv,Ti,xBi,xTv,xTi}
    if (Bi,Tv,Ti) == (xBi,xTv,xTi)
        return x
    else
        return SparseMatrixCSR{Bi}( x.m,
                                    x.n,
                                    convert(Vector{Ti}, copy(getptr(x)).-x.offset),
                                    convert(Vector{Ti}, copy(getindices(x)).-x.offset),
                                    convert(Vector{Tv}, copy(nonzeros(x))))
    end
end

function convert(::Type{SparseMatrixCSR{Bi}}, x::SparseMatrixCSC) where {Bi}
    A = sparse(transpose(x))
    (m, n) = size(A)
    return SparseMatrixCSR{Bi}(m, n, getptr(A), getindices(A), nonzeros(A))
end

function convert(::Type{SparseMatrixCSR{Bi,Tv,Ti}}, x::SparseMatrixCSC) where {Bi,Tv,Ti}
    A = sparse(transpose(x))
    return SparseMatrixCSR{Bi}( A.m,
                                A.n,
                                convert(Vector{Ti}, getptr(A)),
                                convert(Vector{Ti}, getindices(A)),
                                convert(Vector{Tv}, nonzeros(A)))
end

"""
    function convert(::Type{SparseMatrixCSC}, x::SparseMatrixCSR)

Convert x to a value of type SparseMatrixCSC.
"""
function convert(::Type{SparseMatrixCSC{Tv,Ti}}, x::SparseMatrixCSR{xBi,xTv,xTi}) where {Tv,Ti,xBi,xTv,xTi}
    A = sparse(transpose(x))
    return SparseMatrixCSR{Bi}( A.m,
                                A.n,
                                convert(Vector{Ti}, getptr(A)),
                                convert(Vector{Ti}, getindices(A)),
                                convert(Vector{Tv}, nonzeros(A)))
end


text:
function verify_export_results(results, export_path)
    exports = SimulationResultsExport(
        make_export_all(keys(results.problem_results)),
        results.params,
    )
    export_results(results, exports)

    for problem_results in values(results.problem_results)
        problem = problem_results.problem
        rpath = problem_results.results_output_folder
        for timestamp in get_existing_timestamps(problem_results)
            for name in get_existing_duals(problem_results)
                compare_results(rpath, export_path, problem, "duals", name, timestamp)
            end
            for name in get_existing_parameters(problem_results)
                compare_results(rpath, export_path, problem, "parameters", name, timestamp)
            end
            for name in get_existing_variables(problem_results)
                compare_results(rpath, export_path, problem, "variables", name, timestamp)
            end
        end

        # This file is not currently exported during the simulation.
        @test isfile(
            joinpath(
                problem_results.results_output_folder,
                problem_results.problem,
                "optimizer_stats.csv",
            ),
        )
    end
end

function compare_results(rpath, epath, problem, field, name, timestamp)
    filename = string(name) * "_" * PSI.convert_for_path(timestamp) * ".csv"
    df1 = PSI.read_dataframe(joinpath(rpath, problem, field, filename))
    df2 = PSI.read_dataframe(joinpath(epath, problem, field, filename))
    @test df1 == df2
end

function make_export_all(problems)
    return [
        ProblemResultsExport(x, duals = [:all], variables = [:all], parameters = [:all]) for
        x in problems
    ]
end

function test_simulation_results(file_path::String, export_path; in_memory = false)
    @testset "Test simulation results" begin
        template_uc = get_template_hydro_st_uc()
        template_ed = get_template_hydro_st_ed()
        c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_ems_uc")
        c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ems_ed")
        time_series_cache_size = 0  # This is only for test coverage.
        problems = SimulationProblems(
            UC = OperationsProblem(
                template_uc,
                c_sys5_hy_uc;
                optimizer = GLPK_optimizer,
                time_series_cache_size = time_series_cache_size,
            ),
            ED = OperationsProblem(
                template_ed,
                c_sys5_hy_ed;
                optimizer = GLPK_optimizer,
                constraint_duals = [:CopperPlateBalance],
                time_series_cache_size = time_series_cache_size,
            ),
        )

        sequence_cache = SimulationSequence(
            problems = problems,
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                    binary_source_problem = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
                ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                    variable_source_problem = PSI.ACTIVE_POWER,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            cache = Dict(
                ("UC",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON),
                ("UC", "ED") => StoredEnergy(PSY.HydroEnergyReservoir, PSI.ENERGY),
            ),
            ini_cond_chronology = InterProblemChronology(),
        )
        sim = Simulation(
            name = "cache",
            steps = 2,
            problems = problems,
            sequence = sequence_cache,
            simulation_folder = file_path,
        )
        build_out = build!(sim)
        @test build_out == PSI.BuildStatus.BUILT

        exports = Dict(
            "problems" => [
                Dict(
                    "name" => "UC",
                    "variables" => [:all],
                    "parameters" => [:all],
                    "duals" => [:all],
                ),
                Dict(
                    "name" => "ED",
                    "variables" => [:all],
                    "parameters" => [:all],
                    "duals" => [:all],
                ),
            ],
            "path" => export_path,
            "optimizer_stats" => true,
        )
        execute_out = execute!(sim, exports = exports, in_memory = in_memory)
        @test execute_out == PSI.RunStatus.SUCCESSFUL

        results = SimulationResults(sim)
        @test list_problems(results) == ["ED", "UC"]
        results_uc = get_problem_results(results, "UC")
        results_ed = get_problem_results(results, "ED")

        ed_expected_vars = [
            :Sp__HydroEnergyReservoir
            :P__ThermalStandard
            :P__RenewableDispatch
            :E__HydroEnergyReservoir
            :P__InterruptibleLoad
            :P__HydroEnergyReservoir
        ]

        uc_expected_vars = [
            :Sp__HydroEnergyReservoir
            :P__ThermalStandard
            :P__RenewableDispatch
            :start__ThermalStandard
            :stop__ThermalStandard
            :E__HydroEnergyReservoir
            :P__HydroEnergyReservoir
            :On__ThermalStandard
        ]
        if in_memory
            @test IS.get_uuid(get_system(results_uc)) === IS.get_uuid(c_sys5_hy_uc)
            @test IS.get_uuid(get_system(results_ed)) === IS.get_uuid(c_sys5_hy_ed)
        else
            @test get_system(results_uc) === nothing
            @test length(read_realized_variables(results_uc)) == 12 #verifies this works without system
            @test_throws IS.InvalidValue set_system!(results_uc, c_sys5_hy_ed)
            set_system!(results_uc, c_sys5_hy_uc)
            @test IS.get_uuid(get_system!(results_uc)) === IS.get_uuid(c_sys5_hy_uc)
            @test get_system(results_ed) === nothing
            @test IS.get_uuid(get_system!(results_ed)) === IS.get_uuid(c_sys5_hy_ed)

            results_from_file = SimulationResults(joinpath(file_path, "cache"))
            @test list_problems(results) == ["ED", "UC"]
            results_uc_from_file = get_problem_results(results_from_file, "UC")
            results_ed_from_file = get_problem_results(results_from_file, "ED")

            @test isempty(
                setdiff(uc_expected_vars, get_existing_variables(results_uc_from_file)),
            )
            @test isempty(
                setdiff(ed_expected_vars, get_existing_variables(results_ed_from_file)),
            )
        end

        @test isempty(setdiff(uc_expected_vars, get_existing_variables(results_uc)))
        @test isempty(setdiff(ed_expected_vars, get_existing_variables(results_ed)))
        p_thermal_standard_ed = read_variable(results_ed, :P__ThermalStandard)
        @test length(keys(p_thermal_standard_ed)) == 48
        for v in values(p_thermal_standard_ed)
            @test size(v) == (12, 6)
        end

        ren_dispatch_params = read_parameter(
            results_ed,
            :P__max_active_power__RenewableDispatch_max_active_power,
        )
        @test length(keys(ren_dispatch_params)) == 48
        for v in values(ren_dispatch_params)
            @test size(v) == (12, 4)
        end

        network_duals = read_dual(results_ed, :CopperPlateBalance)
        @test length(keys(network_duals)) == 48
        for v in values(network_duals)
            @test size(v) == (12, 2)
        end

        p_variables_uc =
            read_variables(results_uc, names = [:P__RenewableDispatch, :P__ThermalStandard])
        @test length(keys(p_variables_uc)) == 2
        for var_name in values(p_variables_uc)
            for v_ in values(var_name)
                @test size(v_)[1] == 24
            end
        end

        realized_var_uc = read_realized_variables(results_uc)
        @test length(keys(realized_var_uc)) == 12
        for var in values(realized_var_uc)
            @test size(var)[1] == 48
        end

        realized_param_uc = read_realized_parameters(
            results_uc,
            names = [:P__max_active_power__RenewableDispatch_max_active_power],
        )
        @test length(keys(realized_param_uc)) == 1
        for param in values(realized_param_uc)
            @test size(param)[1] == 48
        end

        realized_duals_ed = read_realized_duals(results_ed)
        @test length(keys(realized_duals_ed)) == 1
        for param in values(realized_duals_ed)
            @test size(param)[1] == 576
        end

        # request non sync data
        @test_logs(
            (:error, r"Requested time does not match available results"),
            match_mode = :any,
            @test_throws IS.InvalidValue read_realized_variables(
                results_ed,
                names = [:P__ThermalStandard],
                initial_time = DateTime("2024-01-01T02:12:00"),
                len = 3,
            )
        )

        # request good window
        @test size(
            read_realized_variables(
                results_ed,
                names = [:P__ThermalStandard],
                initial_time = DateTime("2024-01-02T23:10:00"),
                len = 10,
            )[:P__ThermalStandard],
        )[1] == 10

        # request bad window
        @test_logs(
            (:error, r"Requested time does not match available results"),
            (@test_throws IS.InvalidValue read_realized_variables(
                results_ed,
                names = [:P__ThermalStandard],
                initial_time = DateTime("2024-01-02T23:10:00"),
                len = 11,
            ))
        )

        # request bad window
        @test_logs(
            (:error, r"Requested time does not match available results"),
            (@test_throws IS.InvalidValue read_realized_variables(
                results_ed,
                names = [:P__ThermalStandard],
                initial_time = DateTime("2024-01-02T23:10:00"),
                len = 12,
            ))
        )

        load_results!(
            results_ed,
            3,
            initial_time = DateTime("2024-01-01T00:00:00"),
            variables = [:P__ThermalStandard],
        )

        @test !isempty(results_ed.variable_values[:P__ThermalStandard])
        @test length(results_ed.variable_values[:P__ThermalStandard]) == 3
        @test length(results_ed) == 3
        @test length(results) == length(results_ed)

        @test_logs(
            (:error, r"invalid is not stored"),
            @test_throws(IS.InvalidValue, read_parameter(results_ed, :invalid))
        )
        @test_logs(
            (:error, r"invalid is not stored"),
            @test_throws(IS.InvalidValue, read_variable(results_ed, :invalid))
        )
        @test_logs(
            (:error, r"not stored"),
            @test_throws(
                IS.InvalidValue,
                read_variable(results_uc, :P__ThermalStandard; initial_time = now())
            )
        )
        @test_logs(
            (:error, r"not stored"),
            @test_throws(
                IS.InvalidValue,
                read_variable(results_uc, :P__ThermalStandard; count = 25)
            )
        )

        empty!(results_ed)
        @test isempty(results_ed.variable_values[:P__ThermalStandard])

        initial_time = DateTime("2024-01-01T00:00:00")
        load_results!(
            results_ed,
            3,
            initial_time = initial_time,
            variables = [:P__ThermalStandard],
            duals = [:CopperPlateBalance],
            parameters = [:P__max_active_power__RenewableDispatch_max_active_power],
        )

        @test !isempty(results_ed.variable_values[:P__ThermalStandard])
        @test !isempty(results_ed.dual_values[:CopperPlateBalance])
        @test !isempty(
            results_ed.parameter_values[:P__max_active_power__RenewableDispatch_max_active_power],
        )

        @test !isempty(results_ed)
        @test !isempty(results)
        empty!(results)
        @test isempty(results_ed)
        @test isempty(results)

        verify_export_results(results, export_path)

        # Test that you can't read a failed simulation.
        PSI.set_simulation_status!(sim, RunStatus.FAILED)
        PSI.serialize_status(sim)
        @test PSI.deserialize_status(sim) == RunStatus.FAILED
        @test_throws ErrorException SimulationResults(sim)
        @test_logs(
            match_mode = :any,
            (:warn, r"Results may not be valid"),
            SimulationResults(sim, ignore_status = true),
        )

        if in_memory
            container = sim.internal.store.data[:ED].variables[:P__ThermalStandard]
            @test !isempty(container)
            @test !isempty(sim.internal.store.optimizer_stats)
            empty!(sim.internal.store)
            @test isempty(container)
            @test isempty(sim.internal.store.optimizer_stats)
        end
    end

    @testset "Test receding horizon simulation results" begin
        template_uc = get_template_hydro_st_uc()
        c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_ems_uc")
        problems = SimulationProblems(
            UC = OperationsProblem(
                template_uc,
                c_sys5_hy_uc;
                optimizer = GLPK_optimizer,
                constraint_duals = [:CopperPlateBalance],
            ),
        )

        sequence_rh = SimulationSequence(
            problems = problems,
            intervals = Dict("UC" => (Hour(24), RecedingHorizon())),
            ini_cond_chronology = InterProblemChronology(),
        )
        sim = Simulation(
            name = "RH",
            steps = 2,
            problems = problems,
            sequence = sequence_rh,
            simulation_folder = file_path,
        )
        build_out = build!(sim)
        @test build_out == PSI.BuildStatus.BUILT

        execute_out = execute!(sim, in_memory = in_memory)
        @test execute_out == PSI.RunStatus.SUCCESSFUL

        results = SimulationResults(sim)
        @test list_problems(results) == ["UC"]
        results_rh = get_problem_results(results, "UC")

        if in_memory
            @test IS.get_uuid(get_system(results_rh)) === IS.get_uuid(c_sys5_hy_uc)
        else
            @test get_system(results_rh) === nothing
            @test length(read_realized_variables(results_rh)) == 12 #verifies this works without system
            set_system!(results_rh, c_sys5_hy_uc)
            @test IS.get_uuid(get_system!(results_rh)) === IS.get_uuid(c_sys5_hy_uc)
        end

        uc_expected_vars = [
            :Sp__HydroEnergyReservoir
            :P__ThermalStandard
            :P__RenewableDispatch
            :start__ThermalStandard
            :stop__ThermalStandard
            :E__HydroEnergyReservoir
            :P__HydroEnergyReservoir
            :On__ThermalStandard
        ]
        @test isempty(setdiff(uc_expected_vars, get_existing_variables(results_rh)))

        p_thermal_standard_rh = read_variable(results_rh, :P__ThermalStandard)
        @test length(keys(p_thermal_standard_rh)) == 2
        for v in values(p_thermal_standard_rh)
            @test size(v) == (24, 6)
        end

        ren_dispatch_params = read_parameter(
            results_rh,
            :P__max_active_power__RenewableDispatch_max_active_power,
        )
        @test length(keys(ren_dispatch_params)) == 2
        for v in values(ren_dispatch_params)
            @test size(v) == (24, 4)
        end

        if !in_memory
            # this creates a container for duals but doesn't write anything because
            # it's a MIP and the duals are unavailable.
            network_duals = read_dual(results_rh, :CopperPlateBalance)
            @test length(keys(network_duals)) == 2
            for v in values(network_duals)
                @test size(v) == (24, 2)
            end
        end

        realized_var_rh = read_realized_variables(results_rh)
        @test length(keys(realized_var_rh)) == 12
        for var in values(realized_var_rh)
            @test size(var)[1] == 48
            existing_timetsamps = get_existing_timestamps(results_rh)
            for ts in existing_timetsamps
                val_cols = setdiff(propertynames(var), [:DateTime])
                first_row = Matrix(var[var.DateTime .== ts, val_cols])
                all_rows = Matrix(
                    var[
                        (var.DateTime .>= ts) .& (var.DateTime .< ts + existing_timetsamps.step),
                        val_cols,
                    ],
                )
                @test all(first_row .== all_rows)
            end
        end

        realized_param_rh = read_realized_parameters(results_rh)
        @test length(keys(realized_param_rh)) == 4
        for var in values(realized_param_rh)
            @test size(var)[1] == 48
            existing_timetsamps = get_existing_timestamps(results_rh)
            for ts in existing_timetsamps
                val_cols = setdiff(propertynames(var), [:DateTime])
                first_row = Matrix(var[var.DateTime .== ts, val_cols])
                all_rows = Matrix(
                    var[
                        (var.DateTime .>= ts) .& (var.DateTime .< ts + existing_timetsamps.step),
                        val_cols,
                    ],
                )
                @test all(first_row .== all_rows)
            end
        end

        if !in_memory
            realized_duals_rh = read_realized_duals(results_rh)
            @test length(keys(realized_duals_rh)) == 1
            for var in values(realized_duals_rh)
                @test size(var)[1] == 48
                existing_timetsamps = get_existing_timestamps(results_rh)
                for ts in existing_timetsamps
                    val_cols = setdiff(propertynames(var), [:DateTime])
                    first_row = Matrix(var[var.DateTime .== ts, val_cols])
                    all_rows = Matrix(
                        var[
                            (var.DateTime .>= ts) .& (var.DateTime .< ts + existing_timetsamps.step),
                            val_cols,
                        ],
                    )
                    @test all(first_row .== all_rows)
                end
            end
        end

        # request non sync data
        @test_logs(
            (:error, r"Requested time does not match available results"),
            match_mode = :any,
            @test_throws IS.InvalidValue read_realized_variables(
                results_rh,
                names = [:P__ThermalStandard],
                initial_time = DateTime("2024-01-01T02:12:00"),
                len = 3,
            )
        )

        # request good window
        @test size(
            read_realized_variables(
                results_rh,
                names = [:P__ThermalStandard],
                initial_time = DateTime("2024-01-01T23:00:00"),
                len = 10,
            )[:P__ThermalStandard],
        )[1] == 10

        # request bad window
        @test_logs(
            (:error, r"Requested time does not match available results"),
            (@test_throws IS.InvalidValue read_realized_variables(
                results_rh,
                names = [:P__ThermalStandard],
                initial_time = DateTime("2024-01-01T23:00:00"),
                len = 26,
            ))
        )
    end
end

@testset "Test simulation results" begin
    for in_memory in (false, true)
        file_path = mkpath(joinpath(pwd(), "test_simulation_results"))
        export_path = mkpath(joinpath(pwd(), "test_export_path"))
        try
            test_simulation_results(file_path, export_path, in_memory = in_memory)
        finally
            rm(file_path, force = true, recursive = true)
            rm(export_path, force = true, recursive = true)
        end
    end
end

ext: jl


text:
using Neovim
import Neovim: on_notify, on_request, get_buffers, set_line
#test events
struct ChildHandler
end
on_notify(h::ChildHandler, c, name, args) = (println("$name $args"), flush(stdout))
function on_request(h::ChildHandler, c, serial, name, args)
    if name == "run_julia"
        code = args[1]
        reply_result(c, serial, eval(code))
    end
end

nvim = nvim_child(ChildHandler())
buf = get_buffers(nvim)[1]
@assert isa(buf,Buffer)
set_line(buf, 1, "some text")

wait(nvim) #wait for request from nvim


text:
using Test
using OrdinaryDiffEq, Calculus, ForwardDiff

import DiffEqBase: RightRootFind

println("AD Test 1")

function f(du,u,p,t)
  du[1] = -p[1]
  du[2] = p[2]
end

for x in 0:0.001:5
  called = false
  function test_f(p)
    cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(called=true;integrator.p[2]=zero(integrator.p[2])))
    prob = ODEProblem(f,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-14,reltol=1e-14,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end
  p = [2.0, x]
  called = false
  findiff = Calculus.finite_difference_jacobian(test_f,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_f,p)
  @test called
  @test findiff ≈ fordiff
end

println("AD Test 2")

function f2(du,u,p,t)
  du[1] = -u[2]
  du[2] = p[2]
end

for x in 2.1:0.001:5
  called = false
  function test_f2(p)
    cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(called=true;integrator.p[2]=zero(integrator.p[2])))
    prob = ODEProblem(f2,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end
  p = [2.0, x]
  findiff = Calculus.finite_difference_jacobian(test_f2,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_f2,p)
  @test called
  @test findiff ≈ fordiff
end

#=
#x = 2.0 is an interesting case

x = 2.0

function test_f2(p)
  cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(@show(x,integrator.t);called=true;integrator.p[2]=zero(integrator.p[2])))
  prob = ODEProblem(f2,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
  integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
  step!(integrator)
  solve!(integrator).u[end]
end

p = [2.0, x]
findiff = Calculus.finite_difference_jacobian(test_f2,p)
@test called
called = false
fordiff = ForwardDiff.jacobian(test_f2,p)
@test called

# At that value, it shouldn't be called, but a small perturbation will make it called, so finite difference is wrong!
=#

println("AD Test 3")

for x in 1.0:0.001:2.5
  function lotka_volterra(du,u,p,t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α*x - β*x*y
    du[2] = dy = -δ*y + γ*x*y
  end
  u0 = [1.0,1.0]
  tspan = (0.0,10.0)
  p = [x,1.0,3.0,1.0]
  prob = ODEProblem(lotka_volterra,u0,tspan,p)
  sol = solve(prob,Tsit5())

  called=false
  function test_lotka(p)
    cb = ContinuousCallback((u,t,i) -> u[1]-2.5, (integrator)->(called=true;integrator.p[4]=1.5))
    prob = ODEProblem(lotka_volterra,eltype(p).([1.0,1.0]),eltype(p).((0.0,10.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end

  findiff = Calculus.finite_difference_jacobian(test_lotka,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_lotka,p)
  @test called
  @test findiff ≈ fordiff
end

# Gradients and Hessians

println("AD Test 4")

function myobj(θ)
  f(u,p,t) = -θ[1]*u
  u0, _ = promote(10.0, θ[1])
  prob = ODEProblem(f, u0, (0.0, 1.0))
  sol = solve(prob, Tsit5())
  diff = sol.u - 10*exp.(-sol.t)
  return diff'diff
end

ForwardDiff.gradient(myobj, [1.0])
ForwardDiff.hessian(myobj, [1.0])

function myobj2(θ)
  f(du,u,p,t) = (du[1]=-θ[1]*u[1])
  u0, _ = promote(10.0, θ[1])
  prob = ODEProblem(f, [u0], (0.0, 1.0))
  sol = solve(prob, Tsit5())
  diff = sol[:,1] .- 10 .*exp.(-sol.t)
  return diff'diff
end

println("AD Test 5")

ForwardDiff.gradient(myobj2, [1.0])
ForwardDiff.hessian(myobj2, [1.0])

function myobj3(θ)
  f(u,p,t) = -θ[1]*u
  u0, _ = promote(10.0, θ[1])
  tspan_end, _ = promote(1.0, θ[1])
  prob = ODEProblem(f, u0, (0.0, tspan_end))
  sol = solve(prob, Tsit5())
  diff = sol.u - 10*exp.(-sol.t)
  return diff'diff
end

println("AD Test 6")

ForwardDiff.gradient(myobj3, [1.0])
ForwardDiff.hessian(myobj3, [1.0])

function myobj4(θ)
  f(du,u,p,t) = (du[1] = -θ[1]*u[1])
  u0, _ = promote(10.0, θ[1])
  tspan_end, _ = promote(1.0, θ[1])
  prob = ODEProblem(f, [u0], (0.0, tspan_end))
  sol = solve(prob, Tsit5())
  diff = sol[:,1] .- 10 .* exp.(-sol.t)
  return diff'diff
end

println("AD Test 7")

ForwardDiff.gradient(myobj4, [1.0])
ForwardDiff.hessian(myobj4, [1.0])

function terminate_affect!(integrator)
    terminate!(integrator)
end

function terminate_condition(u,t,integrator)
    min(u[1],u[2])
end

terminate_cb = ContinuousCallback(terminate_condition,terminate_affect!)

function damped!(du,u,p,t)
    #Damped harmonic oscillator
    du[1] = u[2]
    du[2] = -p[1]u[2] - p[2]u[1]
end

u0 = [5.,10.] #starting point
tspan = (0.,10.) #time span

function loss(p)
    temp_prob = ODEProblem(damped!, u0, tspan, p)
    temp_solution = solve(temp_prob, Tsit5(), p = p, callback = terminate_cb)
    return temp_solution.t[end]
end

@test ForwardDiff.gradient(loss,[2.,5.]) ≈ [-0.025245268698610393,-0.036292733502239424]

println("Using Right Limit of Rootfind")
println("AD Test 1")

function f(du,u,p,t)
  du[1] = -p[1]
  du[2] = p[2]
end

for x in 0:0.001:5
  called = false
  function test_f(p)
    cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(called=true;integrator.p[2]=zero(integrator.p[2])); rootfind=RightRootFind)
    prob = ODEProblem(f,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-14,reltol=1e-14,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end
  p = [2.0, x]
  called = false
  findiff = Calculus.finite_difference_jacobian(test_f,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_f,p)
  @test called
  @test findiff ≈ fordiff
end

println("AD Test 2")

function f2(du,u,p,t)
  du[1] = -u[2]
  du[2] = p[2]
end

for x in 2.1:0.001:5
  called = false
  function test_f2(p)
    cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(called=true;integrator.p[2]=zero(integrator.p[2])); rootfind=RightRootFind)
    prob = ODEProblem(f2,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end
  p = [2.0, x]
  findiff = Calculus.finite_difference_jacobian(test_f2,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_f2,p)
  @test called
  @test findiff ≈ fordiff
end

#=
#x = 2.0 is an interesting case

x = 2.0

function test_f2(p)
  cb = ContinuousCallback((u,t,i) -> u[1], (integrator)->(@show(x,integrator.t);called=true;integrator.p[2]=zero(integrator.p[2])))
  prob = ODEProblem(f2,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))
  integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
  step!(integrator)
  solve!(integrator).u[end]
end

p = [2.0, x]
findiff = Calculus.finite_difference_jacobian(test_f2,p)
@test called
called = false
fordiff = ForwardDiff.jacobian(test_f2,p)
@test called

# At that value, it shouldn't be called, but a small perturbation will make it called, so finite difference is wrong!
=#

println("AD Test 3")

for x in 1.0:0.001:2.5
  function lotka_volterra(du,u,p,t)
    x, y = u
    α, β, δ, γ = p
    du[1] = dx = α*x - β*x*y
    du[2] = dy = -δ*y + γ*x*y
  end
  u0 = [1.0,1.0]
  tspan = (0.0,10.0)
  p = [x,1.0,3.0,1.0]
  prob = ODEProblem(lotka_volterra,u0,tspan,p)
  sol = solve(prob,Tsit5())

  called=false
  function test_lotka(p)
    cb = ContinuousCallback((u,t,i) -> u[1]-2.5, (integrator)->(called=true;integrator.p[4]=1.5); rootfind=RightRootFind)
    prob = ODEProblem(lotka_volterra,eltype(p).([1.0,1.0]),eltype(p).((0.0,10.0)),copy(p))
    integrator = init(prob,Tsit5(),abstol=1e-12,reltol=1e-12,callback=cb)
    step!(integrator)
    solve!(integrator).u[end]
  end

  findiff = Calculus.finite_difference_jacobian(test_lotka,p)
  @test called
  called = false
  fordiff = ForwardDiff.jacobian(test_lotka,p)
  @test called
  @test findiff ≈ fordiff
end

println("AD Test 4")

function terminate_affect!(integrator)
    terminate!(integrator)
end

function terminate_condition(u,t,integrator)
    min(u[1],u[2])
end

terminate_cb = ContinuousCallback(terminate_condition,terminate_affect!; rootfind=RightRootFind)

function damped!(du,u,p,t)
    #Damped harmonic oscillator
    du[1] = u[2]
    du[2] = -p[1]u[2] - p[2]u[1]
end

u0 = [5.,10.] #starting point
tspan = (0.,10.) #time span

function loss(p)
    temp_prob = ODEProblem(damped!, u0, tspan, p)
    temp_solution = solve(temp_prob, Tsit5(), p = p, callback = terminate_cb)
    return temp_solution.t[end]
end

@test ForwardDiff.gradient(loss,[2.,5.]) ≈ [-0.025245268698610393,-0.036292733502239424]


text:
using Weber
using Base.Test
include("find_timing.jl")

import Weber: poll_events

struct ExtensionA <: Weber.Extension end
struct ExtensionB <: Weber.Extension end
mutable struct ExtensionC <: Weber.Extension
  triggered::Bool
end
ExtensionC() = ExtensionC(false)

import Weber: setup, record, poll_events, run

function setup(fn::Function,e::ExtendedExperiment{ExtensionA})
  setup(next(e)) do
    addcolumn(top(e),:extension_a)
    fn()
  end
end

function record(e::ExtendedExperiment{ExtensionA},code;kwds...)
  record(next(e),code;extension_a = :test,kwds...)
end

function setup(fn::Function,e::ExtendedExperiment{ExtensionB})
  setup(next(e)) do
    addcolumn(top(e),:extension_b)
    fn()
  end
end

function record(e::ExtendedExperiment{ExtensionB},code;kwds...)
  record(next(e),code;extension_b = :test,kwds...)
end

function run(e::ExtendedExperiment{ExtensionB};keys...)
  record(top(e),:extension_b_run)
  run(next(e);keys...)
end

function poll_events(callback,e::ExtendedExperiment{ExtensionC},time::Float64)
  if !extension(e).triggered
    extension(e).triggered = true
    record(top(e),:extension_c_polled)
  end
  poll_events(callback,next(e),time)
end

_,_,extension_events1 = find_timing(extensions=[ExtensionA(),ExtensionB(),ExtensionC()]) do
  addtrial(moment(1ms))
end

_,_,extension_events2 = find_timing(extensions=[ExtensionC(),ExtensionB(),ExtensionA()]) do
  addtrial(moment(1ms))
end

@testset "Extensions" begin
  extension_event1 = filter(r -> r[:code] == :extension_c_polled,extension_events1)
  extension_event2 = filter(r -> r[:code] == :extension_c_polled,extension_events2)
  @test extension_event1[1][:extension_a] == :test
  @test extension_event1[1][:extension_b] == :test
  @test extension_event2[1][:extension_a] == :test
  @test extension_event2[1][:extension_b] == :test

  extension_b_run1 = filter(r -> r[:code] == :extension_b_run,extension_events1)
  extension_b_run2 = filter(r -> r[:code] == :extension_b_run,extension_events2)
  @test length(extension_b_run1) == 1
  @test length(extension_b_run2) == 1
end


text:
# This is where most of the code from PeaceVote will go

# abstract type AbstractPort end

### Could have additional information on regulation
struct BraiderConfig
    N::UInt8
    M::UInt8
    gk::GateKeeper
end

# DHasym, DHsym

struct BraidOfficer <: Officer
    sb::BraiderConfig
    bc::BraidChain
end

### This module actually specifies seal method with fold and unfold. I could use thoose for unfolding data to a Dictionary. It should be arguable whether it would not be better to type each signature

function SynchronicBallot.regulation(bo::BraidOfficer)

    N, M = bo.sb.N, bo.sb.M
    mix = bo.sb.gk.mix.id
    gk = bo.sb.gk.id
    hash = Hash(bo.bc)

    metadata = MetaData(mix, gk, hash)
    mbytes = fold(Dict(metadata))

    return Regulation(N, M, mbytes)
end

Base.in(id, bo::BraidOfficer) = id in pseudonyms(bo.bc)

function SynchronicBallot.audit!(bo::BraidOfficer, ballot::Ballot, signatures)

    braid = SynchronicBraid(ballot)
    contract = Contract(Dict[unfold(i) for i in signatures])
    transaction = Transaction(braid, contract)
    push!(bo.bc, transaction)

    return nothing
end

function validate(sb::SynchronicBraid, bo::BraidOfficer, id)
    bc, tlog = bo.sb, bo.bc

    metadata = sb.metadata

    validate(tlog, metadata.hash) || return false
    bc.gk.id == metadata.gk || return false
    bc.gk.mix.id == metadata.mix || return false

    id in sb.pseudonyms || return false

    return true
end

serve!(bo::BraidOfficer, signer::Signer) = SynchronicBallot.serve!(bo.sb.gk, bo, signer)

function braid(oldsigner::Signer, newsigner::Signer, bo::BraidOfficer)

    newid = id(newsigner)
    idbytes = Vector{UInt8}(newid, base=16, length=bo.sb.M)

    signbraid = (ballot) -> begin
        braid = SynchronicBraid(ballot)
        @assert validate(braid, bo, newid)
        signature = sign(braid, oldsigner)
        fold(Dict(signature))
    end

    vote(idbytes, bo.sb.gk, oldsigner, signbraid)

    return nothing
end


text:
module Bifurcation

using PyPlot

include("model/model.jl")
using .Model
if !isfile("./differential_equation.jl")
    run(`sh create_diffeq.sh`)
end
include("./differential_equation.jl")
include("../../continuation.jl")

function analysis()
    global fp,br

    # BP: name(index) of bifurcation parameter
    global BP = C.S

    global p = param_values()
    new_curve!(p, direction=false)
    fp = readdlm("./Data/fp.dat",'\t',Float64,'\n')
    ev = readdlm("./Data/ev.dat",'\t',Float64,'\n')
    br = get_bistable_regime(ev)
end

function diagram()
    rc("figure",figsize = (8,6))
    rc("font",family = "Arial")
    rc("font",size = 24)
    rc("axes",linewidth = 1)
    rc("xtick.major",width = 1)
    rc("ytick.major",width = 1)
    rc("lines",linewidth = 3)

    plot(fp[1:br[1]-1,VN+1],fp[1:br[1]-1,V.E+1],"k-")
    plot(fp[br,VN+1],fp[br,V.E+1],lw=1.5,"k--")
    plot(fp[br[end]+1:end,VN+1],fp[br[end]+1:end,V.E+1],"k-")

    xlabel("Serum (percentage)")
    ylabel("E2F (μM)")

    xlim(0,2)
    xticks([0,0.5,1,1.5,2])
    yscale("log")
    ylim(1e-4,2)
    yticks([1e-4,1e-2,1])

    show()
end

end  # module

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
###############################################################################
#
# Update xylem pressure and capacitance
#
###############################################################################
"""
    update_PVF!(hs::LeafHydraulics{FT}, Δt::FT) where {FT<:AbstractFloat}
    update_PVF!(hs::StemHydraulics{FT}, Δt::FT) where {FT<:AbstractFloat}
    update_PVF!(hs::RootHydraulics{FT}, Δt::FT) where {FT<:AbstractFloat}
    update_PVF!(hs::RootHydraulics{FT},
                Δt::FT,
                nss::Bool
    ) where {FT<:AbstractFloat}
    update_PVF!(roots::Array{RootHydraulics{FT},1},
                ks::Array{FT,1},
                ps::Array{FT,1},
                qs::Array{FT,1},
                q_sum::FT,
                Δt::FT
    ) where {FT<:AbstractFloat}
    update_PVF!(tree::GrassLikeOrganism{FT}, Δt::FT) where {FT<:AbstractFloat}
    update_PVF!(tree::PalmLikeOrganism{FT}, Δt::FT) where {FT<:AbstractFloat}
    update_PVF!(tree::TreeLikeOrganism{FT}, Δt::FT) where {FT<:AbstractFloat}

Update pressure, capacitance, and flow rates in hydraulic system, given
- `hs` [`AbstractHydraulicOrgan`](@ref) or [`AbstractPlantOrganism`](@ref) type
    struct
- `Δt` Time interval
- `nss` NonSteadyStateMode placeholder (useless)
- `ks` Container for conductance in each root layer
- `ps` Container for end xylem pressure in each layer
- `qs` Container for flow rate out of each layer
- `q_sum` Total flow rate out of the roots
- `tree` [`GrassLikeOrganism`](@ref), [`PalmLikeOrganism`](@ref), or
    [`TreeLikeOrganism`](@ref) type struct

Note that this function only updates the equilibrium pressure in the tissue,
    but not the xylem flow pressure. The difference between the two pressures
    is used to drive water exchange between xylem and capacictance tissues.
"""
function update_PVF!(hs::LeafHydraulics{FT}, Δt::FT) where {FT<:AbstractFloat}
    # unpack values
    @unpack f_vis, p_element, p_leaf, q_out, pv, T_sap, v_maximum,
            v_storage = hs;

    # calculate the flow rate out of the capacitance
    # use buffer_rate here to avoid alloaction issues
    f_cap = (hs.p_storage - p_leaf) * buffer_rate(pv) / f_vis * v_maximum;
    if (f_cap > 0) && (hs.v_storage <= f_cap * Δt)
        f_cap = hs.v_storage / Δt;
    end

    # update flow into the tissue, storage and the tissue pressure (p_storage)
    hs.q_in       = q_out - f_cap;
    hs.v_storage -= f_cap * Δt;
    hs.p_storage  = p_from_volume(pv, hs.v_storage/v_maximum, T_sap);

    return nothing
end

function update_PVF!(hs::StemHydraulics{FT}, Δt::FT) where {FT<:AbstractFloat}
    # unpack values
    @unpack f_vis, N, q_element, p_element, p_storage, pv, q_out, T_sap,
            v_maximum, v_storage = hs;

    # update flow rate in, storage volume and pressure per slice
    f_sum::FT = 0;
    for i in N:-1:1
        f_cap = (p_storage[i] - p_element[i]) * buffer_rate(pv) / f_vis *
                v_maximum[i];
        if (f_cap > 0) && (v_storage[i] <= f_cap * Δt)
            f_cap = v_storage[i] / Δt;
        end
        v_storage[i] -= f_cap * Δt;
        p_storage[i]  = p_from_volume(pv, v_storage[i]/v_maximum[i], T_sap);
        q_element[i]  = q_out - f_sum;
        f_sum        += f_cap;
    end
    hs.q_in = q_out - f_sum;

    return nothing
end

function update_PVF!(hs::RootHydraulics{FT}, Δt::FT) where {FT<:AbstractFloat}
    # unpack values
    @unpack f_vis, N, p_element, p_storage, pv, q_buffer, q_diff, q_element,
            q_out, T_sap, v_maximum, v_storage = hs;

    # update flow rate in, storage volume and pressure per slice
    f_sum::FT = 0;
    for i in N:-1:1
        f_cap = (p_storage[i] - p_element[i]) * buffer_rate(pv) / f_vis *
                v_maximum[i];
        if (f_cap > 0) && (v_storage[i] <= f_cap * Δt)
            f_cap = v_storage[i] / Δt;
        end
        q_buffer[i]   = f_cap;
        q_diff[i]     = f_sum;
        v_storage[i] -= f_cap * Δt;
        p_storage[i]  = p_from_volume(pv, v_storage[i]/v_maximum[i], T_sap);
        q_element[i]  = q_out - f_sum;
        f_sum        += f_cap;
    end
    hs.q_in = q_out - f_sum;

    return nothing
end

function update_PVF!(
            hs::RootHydraulics{FT},
            Δt::FT,
            nss::Bool
) where {FT<:AbstractFloat}
    # unpack values
    @unpack f_vis, N, p_element, p_storage, pv, q_buffer, q_diff, q_element,
            q_out, v_maximum, v_storage = hs;

    # update flow rate in, storage volume and pressure per slice
    f_sum::FT = 0;
    for i in N:-1:1
        f_cap = (p_storage[i] - p_element[i]) * buffer_rate(pv) /f_vis *
                v_maximum[i];
        if (f_cap > 0) && (v_storage[i] <= f_cap * Δt)
            f_cap = v_storage[i] / Δt;
        end
        q_buffer[i] = f_cap;
        q_diff[i]   = f_sum;
        f_sum      += f_cap;
    end

    return nothing
end

function update_PVF!(
            roots::Array{RootHydraulics{FT},1},
            cache_k::Array{FT,1},
            cache_p::Array{FT,1},
            cache_q::Array{FT,1},
            q_sum::FT,
            Δt::FT
) where {FT<:AbstractFloat}
    # update flow rate for roots, partition total flow rates into roots
    # algorithm from function `roots_flow!`
    # 1. update root buffer rate
    for root in roots
        update_PVF!(root, Δt, true);
    end

    # 2. calculate the flow rate in each root layer
    count = 0;
    while true
        count += 1;
        if count > 20
            break
        end

        # 2.1 flush the values to ks, ps, and qs
        for i in eachindex(roots)
            root  = roots[i];
            _q_in = cache_q[i] - sum(root.q_buffer);
            root.q_element .= cache_q[i] .- root.q_diff;
            _p,_k = root_pk(root, _q_in, root.q_element);
            cache_p[i] = _p;
            cache_k[i] = _k;
        end

        # 2.2 use ps and ks to compute the Δq to adjust
        pm = mean(cache_p);
        for i in eachindex(roots)
            cache_q[i] -= (pm - cache_p[i]) * cache_k[i]
        end

        # 2.3 adjust the qs so that sum(qs) = f_sum
        q_diff = sum(cache_q) - q_sum;
        if (abs(q_diff) < max(FT(1e-6), eps(FT))) &&
           (maximum(cache_p) - minimum(cache_p) < 1e-4)
            break
        end
        k_sum  = sum(cache_k);
        for i in eachindex(roots)
            cache_q[i] -= q_diff * cache_k[1] / k_sum;
        end
    end

    # 2.4 update root PVF again
    for i in eachindex(roots)
        root = roots[i];
        root.q_out = cache_q[i];
        update_PVF!(root, Δt);
    end

    return nothing
end

function update_PVF!(
            tree::GrassLikeOrganism{FT},
            Δt::FT
) where {FT<:AbstractFloat}
    @unpack cache_k, cache_p, cache_q, leaves, roots = tree;

    # 0. note that leaf flow rates need to be updated outside this function
    # 1. update leaf PVF for each layer
    f_sum::FT = 0;
    for leaf in leaves
        update_PVF!(leaf, Δt);
        f_sum += leaf.q_in * leaf.area;
    end

    # 2. update flow rate for roots, partition total flow rates into roots
    update_PVF!(roots, cache_k, cache_p, cache_q, f_sum, Δt);

    # 3. update pressure profiles
    pressure_profile!(tree, NonSteadyStateMode(); update=false);

    return nothing
end

function update_PVF!(
            tree::PalmLikeOrganism{FT},
            Δt::FT
) where {FT<:AbstractFloat}
    @unpack cache_k, cache_p, cache_q, leaves, roots, trunk = tree;

    # 0. note that leaf flow rates need to be updated outside this function
    # 1. update leaf PVF for each layer
    f_sum::FT = 0;
    for leaf in leaves
        update_PVF!(leaf, Δt);
        f_sum += leaf.q_in * leaf.area;
    end

    # 2. update PVF for trunk
    trunk.q_out = f_sum;
    update_PVF!(trunk, Δt);

    # 3. update flow rate for roots, partition total flow rates into roots
    update_PVF!(roots, cache_k, cache_p, cache_q, trunk.q_in, Δt);

    # 4. update pressure profiles
    pressure_profile!(tree, NonSteadyStateMode(); update=false);

    return nothing
end

function update_PVF!(
            tree::TreeLikeOrganism{FT},
            Δt::FT
) where {FT<:AbstractFloat}
    @unpack branch, cache_k, cache_p, cache_q, leaves, roots,
            trunk = tree;

    # 0. note that leaf flow rates need to be updated outside this function
    # 1. update leaf and stem PVF for each layer
    f_sum::FT = 0;
    for i in eachindex(leaves)
        leaf = leaves[i];
        stem = branch[i];
        update_PVF!(leaf, Δt);
        stem.q_out = leaf.q_in * leaf.area;
        update_PVF!(stem, Δt);
        f_sum += stem.q_in;
    end

    # 2. update PVF for trunk
    trunk.q_out = f_sum;
    update_PVF!(trunk, Δt);

    # 3. update flow rate for roots, partition total flow rates into roots
    update_PVF!(roots, cache_k, cache_p, cache_q, trunk.q_in, Δt);

    # 4. update pressure profiles
    pressure_profile!(tree, NonSteadyStateMode(); update=false);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#Parameters used to estimate canopy transpiration
mutable struct TranspirationParameter{T<:Float64}
    ρₐ::T #density of dry air (Kg m⁻³)
    cₚₐ::T #Specific heat capacity of dry air (J kg⁻¹ K⁻¹)
    γ::T #Psychrometric constant (Pa K⁻¹)
    λ::T #Latent heat of water vaporization (J kg⁻¹)
end
TranspirationParameter(;ρₐ::T=1.204,cₚₐ::T=1004.0,γ::T=66.1,λ::T=2.454e6) where {T<:Float64} = TranspirationParameter(ρₐ,cₚₐ,γ,λ)

#derivative of saturation Vapor Pressure (Pa) with respect to temperature at T degree celsius
#Derived from Equation 21 in Alduchov and Eskridge (1996)
dSVPₜ(T::AbstractFloat) = 610.94*17.625*243.04*exp(17.625*T/(243.04+T))/(243.04+T)^2

function CalcE_eq(ϕₙₐ::T,Tₐ::T;par::TranspirationParameter=TranspirationParameter()) where {T<:Float64}
    #Equation for calcualting the equilibrium evaporation rate E_eq (Kg s⁻¹ m⁻²)
    #ϕₙₐ net radiation absorbed by the canopy (J s⁻¹ m⁻²)
    #Tₐ air temperature (°C)
    s = dSVPₜ(Tₐ)
    E_eq = s*ϕₙₐ/((s+par.γ)*par.λ)
    return E_eq
end

function CalcE_imp(g_C::T,VPD::T;par::TranspirationParameter=TranspirationParameter()) where {T<:Float64}
    #Equation for calcualting the imposed evaporation rate E_imp (Kg s⁻¹ m⁻²)
    #g_C canopy conductance (water based conductance) (m³ s⁻¹ m⁻²)
    #VPD vapour-pressure deficit (Pa)
    E_imp = par.ρₐ*par.cₚₐ*g_C*VPD/(par.γ*par.λ)
    return E_imp
end

function CalcE(ϕₙₐ::T,Tₐ::T,g_C::T,VPD::T,Ω::T;par::TranspirationParameter=TranspirationParameter()) where {T<:Float64}
    #Equaiton for estimating the canopy tranpiraiton E (Kg s⁻¹ m⁻²)
    #ϕₙₐ net radiation absorbed by the canopy (J s⁻¹ m⁻²)
    #Tₐ air temperature (°C)
    #g_C canopy conductance (water based conductance) (m³ s⁻¹ m⁻²)
    #VPD vapour-pressure deficit (Pa)
    #Ω a dimentionless number in the range 0.0-1.0
    E_eq = CalcE_eq(ϕₙₐ,Tₐ;par=par)
    E_imp = CalcE_imp(g_C,VPD;par=par)
    E = Ω*E_eq+(1- Ω)*E_imp
    return E
end

function Estg_C(E::T,ϕₙₐ::T,Tₐ::T,VPD::T,Ω::T;par::TranspirationParameter=TranspirationParameter()) where {T<:Float64}
    #Equaiton for estimating the canopy conductance g_C (m³ s⁻¹ m⁻²)
    #E canopy tranpiraiton (Kg s⁻¹ m⁻²)
    #ϕₙₐ net radiation absorbed by the canopy (J s⁻¹ m⁻²)
    #Tₐ air temperature (°C)
    #VPD vapour-pressure deficit (Pa)
    #Ω a dimentionless number in the range 0.0-1.0
    E_eq = CalcE_eq(ϕₙₐ,Tₐ;par)
    g_C = (E-Ω*E_eq)/(1-Ω)*(par.γ*par.λ)/(par.ρₐ*par.cₚₐ*VPD)
    return g_C
end

function PAR2Rad(PAR::T;conv_coeff::T=10^6/2.3) where {T<:Float64}
    #Funciton for converting From PAR (mol s⁻¹ m⁻²) to global
    #radiation (Rad, J s⁻¹ m⁻²)
    Rad = PAR*conv_coeff
    return Rad
end

function gₛ2g_C(gₛ::T,PAR::T,LAI::T;cons::CCPH.Constants=CCPH.Constants(),treepar::CCPH.TreePar=CCPH.TreePar(),
    I_gₛ::T=50.0e-6) where {T<:Float64}
    #Estimates canopy conductance g_C (m³ s⁻¹ m⁻²)
    #gₛ stomatal conductanceg (mol CO₂ s⁻¹ m⁻²)
    #PAR at canopy top (mol s⁻¹ m⁻²)
    #LAI Canopy leaf area (m² leaf area m⁻² ground area)
    #I_gₛ Irradiace level at which gₛ is half of canopy top gₛ (mol s⁻¹ m⁻²)

    g_C = cons.r*cons.M_H2O/cons.ρ_H2O*gₛ*log((I_gₛ+PAR)/(I_gₛ+PAR*exp(-treepar.k*LAI)))/treepar.k

    return g_C
end

function Calctranpiration(VPD_z::T,Sₑ::T;E_cm_ref::T=1.812,s_VPD_z::T=3.121,
    s_Sₑ::T=18.342) where {T<:Float64}
    #Tor-Ngern et al. 2017 model for estimating canopy transpiration E_c (mm/d)
    #Standard paramater values (E_cm_ref,s_VPD_z,s_Sₑ) are taken from the
    #original paper for the Rosinedal scots pine stand
    #--Input--
    #VPD_z vapor pressure deficit (Pa)
    #Sₑ effective saturation (-)

    E_cm = E_cm_ref*(1-exp(-s_VPD_z*VPD_z*10^-3))
    E_c = E_cm*(1-exp(-s_Sₑ*Sₑ))

    return E_c
end

function Est_gₛ(LAI::T,daylight::T;cons::CCPH.Constants=CCPH.Constants(),
    env::CCPH.EnvironmentStruct=CCPH.EnvironmentStruct()) where {T<:Float64}
    # Estimate stomatal conductance gₛ (mol C m⁻² leaf area s⁻¹) based on
    # Weather data and leaf area index, using the canopy transpiration from
    # Tor-Ngern et al. 2017
    #--Input--
    # LAI leaf area index (m² leaf area m⁻² ground area)
    # daylight lenght of a single day (s)
    VPD = env.VPD
    θₛ = env.θₛ
    P = env.P
    VPD_z = VPD.*daylight/(24*3600)
    Sₑ = CCPH.CalcSₑ.(θₛ)
    E_c = Calctranpiration.(VPD_z,Sₑ) #mm/day
    E_c *= 10^-3*cons.ρ_H2O/cons.M_H2O/daylight #mol m⁻² s⁻¹
    gₛ = E_c*P/(VPD*cons.r*LAI)
    return gₛ
end


text:
#=

Log-likelihood for ODEs, integrating along the tree

Ignacio Quintero Mächler

t(-_-t)

Created 20 09 2017
Updated 26 03 2019

=#

"""
    make_lhf(llf ::Function,
             lpf ::Function,
             assign_hidfacs!::Function,
             dcp ::Dict{Int64,Int64},
             dcfp::Dict{Int64,Int64})

Make log posterior function with the likelihood, **llf**,
and prior, **lpf**, functions.
"""
function make_lhf(llf            ::Function,
                  lpf            ::Function,
                  assign_hidfacs!::Function,
                  dcp            ::Dict{Int64,Int64},
                  dcfp           ::Dict{Int64,Int64})

  ks  = keys(dcp)
  ksf = keys(dcfp)

  function f(p ::Array{Float64,1},
             fp::Array{Float64,1})

    # constraints
    for wp in ks
      while haskey(dcp, wp)
        tp = dcp[wp]
        p[tp] = p[wp]
        wp = tp
      end
    end

    for wp in ksf
      while haskey(dcfp, wp)
        tp = dcfp[wp]
        fp[tp] = fp[wp]
        wp = tp
      end
    end

   # factors
    assign_hidfacs!(p, fp)

    return llf(p) + lpf(p, fp)
  end

  return f
end

"""
    make_loglik(X        ::Array{Array{Float64,1},1},
                abts1    ::Array{Float64,1},
                abts2    ::Array{Float64,1},
                trios    ::Array{Array{Int64,1},1},
                ode_solve::Function,
                λevent!  ::Function,
                rootll   ::Function,
                k        ::Int64,
                h        ::Int64,
                ns       ::Int64)

Make likelihood function for a tree given an ODE function.
"""
function make_loglik(X        ::Array{Array{Float64,1},1},
                     abts1    ::Array{Float64,1},
                     abts2    ::Array{Float64,1},
                     trios    ::Array{Array{Int64,1},1},
                     int      ::DiffEqBase.DEIntegrator,
                     λevent!  ::Function,
                     rootll   ::Function,
                     k        ::Int64,
                     h        ::Int64,
                     ns       ::Int64,
                     ned      ::Int64)

  # preallocate vectors
  llik = Array{Float64,1}(undef, ns)
  w    = Array{Float64,1}(undef, ns)
  extp = Array{Float64,1}(undef, ns)

  function f(p::Array{Float64,1})

    @inbounds begin

      int.p = p

      llxtra = 0.0

      # loop for integrating over internal branches
      for triad in trios

        pr, d1, d2 = triad::Array{Int64,1}

        ud1 = @views solvef(int, X[d1], abts2[d1], abts1[d1])::Array{Float64,1}

        check_negs(ud1, ns) && return -Inf

        ud2 = @views solvef(int, X[d2], abts2[d2], abts1[d2])::Array{Float64,1}

        check_negs(ud2, ns) && return -Inf

        # update likelihoods with speciation event
        λevent!(abts2[pr], llik, ud1, ud2, p)

        # loglik to sum for integration
        tosum = normbysum!(llik, ns)

        llxtra += log(tosum)

        # assign the remaining likelihoods &
        # assign extinction probabilities and
        # check for extinction of `1.0`
        @views Xpr = X[pr]
        for i in Base.OneTo(ns)
          ud1[i+ns] > 1.0 && return -Inf
          Xpr[i+ns] = ud1[i+ns]
          Xpr[i]    = llik[i]
        end

      end

      # assign root likelihood in non log terms &
      # assign root extinction probabilities
      @views Xned = X[ned]
      for i in Base.OneTo(ns)
        llik[i] = Xned[i]
        extp[i] = Xned[i+ns]
      end

      # estimate likelihood weights
      normbysum!(llik, w, ns)

      # combine root likelihoods
      ll = rootll(abts1[ned], llik, extp, w, p)

      return (log(ll) + llxtra)::Float64
    end
  end

  return f
end

"""
    make_rootll(t   ::Float64,
                llik::Array{Float64,1},
                extp::Array{Float64,1},
                w   ::Array{Float64,1},
                p   ::Array{Float64,1},
                lλs ::Array{Float64,1},
                λts::Array{Float64,1},
                r   ::Array{Float64,1},
                ::Val{k},
                ::Val{h},
                ::Val{ny},
                ::Val{mdS}) where {k,h,ny,mdS}

Generated function for full tree likelihood at the root for normal likelihoods.
"""
@generated function rootll_full(t   ::Float64,
                                llik::Array{Float64,1},
                                extp::Array{Float64,1},
                                w   ::Array{Float64,1},
                                p   ::Array{Float64,1},
                                λts ::Array{Float64,1},
                                r   ::Array{Float64,1},
                                af! ::Function,
                                ::Val{k},
                                ::Val{h},
                                ::Val{ny},
                                ::Val{model})::Float64 where {k, h, ny, model}

  S = create_states(k, h)

  eqs = quote end
  popfirst!(eqs.args)

  if model[1]
   # last non β parameters
    bbase = h*(k+1) + 2*k*h + k*(k-1)*h + h*(h-1)

    # ncov
    ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

    # y per parameter
    yppar = isone(ny) ? 1 : div(ny,ncov)

    # add environmental function
    push!(eqs.args, :(af!(t, r)))

    isone(ny) && push!(eqs.args, :(r1 = r[1]))

    for j = Base.OneTo(h), i = Base.OneTo(k)
      coex = Expr(:call, :+)
      for yi in Base.OneTo(yppar)
        rex = isone(ny) ? :(r1) : :(r[$(yi+yppar*(i-1))])
        push!(coex.args, :(p[$(bbase + yi + yppar*((i-1) + k*(j-1)))] * $rex))
      end
      push!(eqs.args, :(λts[$(i + k*(j-1))] = p[$(i + (k+1)*(j-1))]*exp($coex)))
    end

    eq = Expr(:call, :+)
    for j in Base.OneTo(h)
      # for single areas
      for i in Base.OneTo(k)
        push!(eq.args,
          :(llik[$(i + (2^k-1)*(j-1))]*w[$(i + (2^k-1)*(j-1))] /
            (λts[$(i + k*(j-1))]*(1.0-extp[$(i + (j-1)*(2^k-1))])^2)))
      end
      # for widespread
      for i in k+1:2^k-1
        wl = :(llik[$(i + (2^k-1)*(j-1))]*w[$(i + (2^k-1)*(j-1))] /(le))

        lams = Expr(:call, :+)
        # single area speciation
        for a in S[i + (2^k-1)*(j-1)].g
          push!(lams.args, :(
            λts[$(a + k*(j-1))] *
           (1.0-extp[$(a + (j-1)*(2^k-1))]) * (1.0-extp[$(i + (j-1)*(2^k-1))])))
        end

        # for allopatric speciation
        for (la, ra) = vicsubsets(S[i + (2^k-1)*(j-1)].g)[1:(div(end,2))]
          push!(lams.args,
            :(p[$((k+1)*(1+(j-1)))] *
             (1.0-extp[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*(j-1))]) *
             (1.0-extp[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*(j-1))])))
        end

        # change l in wl
        wl.args[3] = lams

        push!(eq.args, wl)
      end
    end

  else

    eq = Expr(:call, :+)
    for j in Base.OneTo(h)
      # for single areas
      for i in Base.OneTo(k)
        push!(eq.args,
          :(llik[$(i + (2^k-1)*(j-1))]*w[$(i + (2^k-1)*(j-1))] /
            (p[$(i + (k+1)*(j-1))]*(1.0-extp[$(i + (j-1)*(2^k-1))])^2)))
      end
      # for widespread
      for i in k+1:2^k-1
        wl = :(llik[$(i + (2^k-1)*(j-1))]*w[$(i + (2^k-1)*(j-1))] /(le))

        lams = Expr(:call, :+)
        # single area speciation
        for a in S[i + (2^k-1)*(j-1)].g
          push!(lams.args, :(
            p[$(a + (k+1)*(j-1))] *
           (1.0-extp[$(a + (j-1)*(2^k-1))]) * (1.0-extp[$(i + (j-1)*(2^k-1))])))
        end

        # for allopatric speciation
        for (la, ra) = vicsubsets(S[i + (2^k-1)*(j-1)].g)[1:(div(end,2))]
          push!(lams.args,
            :(p[$((k+1)*(1+(j-1)))] *
             (1.0-extp[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*(j-1))]) *
             (1.0-extp[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*(j-1))])))
        end

        # change l in wl
        wl.args[3] = lams

        push!(eq.args, wl)
      end
    end
  end

  push!(eqs.args, :($eq))

  #Core.println(eqs)

  return quote
    @inbounds begin
      $eqs
    end
  end
end

"""
    λevent_full(t   ::Float64,
                llik::Array{Float64,1},
                ud1 ::Array{Float64,1},
                ud2 ::Array{Float64,1},
                lλs ::Array{Float64,1},
                λts::Array{Float64,1},
                p   ::Array{Float64,1},
                r   ::Array{Float64,1},
                ::Val{k},
                ::Val{h},
                ::Val{ny},
                ::Val{S},
                ::Val{mdS})

Generated function for speciation event likelihoods for normal likelihood.
"""
@generated function λevent_full(t   ::Float64,
                                llik::Array{Float64,1},
                                ud1 ::Array{Float64,1},
                                ud2 ::Array{Float64,1},
                                p   ::Array{Float64,1},
                                λts ::Array{Float64,1},
                                r   ::Array{Float64,1},
                                af! ::Function,
                                ::Val{k},
                                ::Val{h},
                                ::Val{ny},
                                ::Val{model}) where {k, h, ny, model}

  eqs = quote end
  popfirst!(eqs.args)

  S = create_states(k, h)

  # if speciation model
  if model[1]

    # last non β parameters
    bbase = h*(k+1) + 2*k*h + k*(k-1)*h + h*(h-1)

    # ncov
    ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

    # y per parameter
    yppar = isone(ny) ? 1 : div(ny,ncov)

    # add environmental function
    push!(eqs.args, :(af!(t, r)))

    isone(ny) && push!(eqs.args, :(r1 = r[1]))

    for j = Base.OneTo(h), i = Base.OneTo(k)
      coex = Expr(:call, :+)
      for yi in Base.OneTo(yppar)
        rex = isone(ny) ? :(r1) : :(r[$(yi+yppar*(i-1))])
        push!(coex.args, :(p[$(bbase + yi + yppar*((i-1) + k*(j-1)))] * $rex))
      end
      push!(eqs.args, :(λts[$(i + k*(j-1))] = p[$(i + (k+1)*(j-1))]*exp($coex)))
    end

    # likelihood for individual areas states
    for j = Base.OneTo(h), i = Base.OneTo(k)
      push!(eqs.args,
          :(llik[$(i+(2^k-1)*(j-1))] =
            ud1[$(i+(2^k-1)*(j-1))] * ud2[$(i+(2^k-1)*(j-1))] * λts[$(i+k*(j-1))]))
    end

    # likelihood for widespread states
    for j = Base.OneTo(h), i = k+1:2^k-1
      s = S[i + (2^k-1)*(j-1)]

      # within region speciation
      ex = Expr(:call, :+)
      for a in s.g
        push!(ex.args, :((ud1[$(a + s.h*(2^k-1))]   * ud2[$(i + (2^k-1)*(j-1))] +
                          ud1[$(i + (2^k-1)*(j-1))] * ud2[$(a + s.h*(2^k-1))]) *
                          0.5 * λts[$(a + s.h*(k))]))
      end

      # between region speciation
      for (la, ra) = vicsubsets(s.g)[1:div(end,2)]
        push!(ex.args,
         :((ud1[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] +
            ud1[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)]) *
            0.5 * p[$((k+1)*(1+s.h))]))
      end

      push!(eqs.args, :(llik[$(i + (2^k-1)*(j-1))] = $ex))
    end

  # *not* speciation model
  else

    # likelihood for individual areas states
    for j = Base.OneTo(h), i = Base.OneTo(k)
      push!(eqs.args,
          :(llik[$(i+(2^k-1)*(j-1))] =
            ud1[$(i+(2^k-1)*(j-1))] * ud2[$(i+(2^k-1)*(j-1))] * p[$(i+(k+1)*(j-1))]))
    end

    # likelihood for widespread states
    for j = Base.OneTo(h), i = k+1:2^k-1
      s = S[i + (2^k-1)*(j-1)]

      # within region speciation
      ex = Expr(:call, :+)
      for a in s.g
        push!(ex.args, :((ud1[$(a + s.h*(k+1))]*ud2[$(i + (2^k-1)*(j-1))] +
                          ud1[$(i + (2^k-1)*(j-1))]*ud2[$(a + s.h*(k+1))]) *
                          0.5 * p[$(a + s.h*(k+1))]))
      end

      # between region speciation
      for (la, ra) = vicsubsets(s.g)[1:div(end,2)]
        push!(ex.args,
         :((ud1[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] +
            ud1[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)]) *
            0.5 * p[$((k+1)*(1+s.h))]))
      end

      push!(eqs.args, :(llik[$(i + (2^k-1)*(j-1))] = $ex))
    end
  end

  #Core.println(eqs)

  return quote
    @inbounds begin
      $eqs
    end
    return nothing
  end
end

"""
    λevent_full(t   ::Float64,
                llik::Array{Array{Float64,1},1},
                ud1 ::Array{Float64,1},
                ud2 ::Array{Float64,1},
                p   ::Array{Float64,1},
                pr  ::Int64,
                λts ::Array{Float64,1},
                r   ::Array{Float64,1},
                af! ::Function,
                ::Val{k},
                ::Val{h},
                ::Val{ny},
                ::Val{model}) where {k, h, ny, model}

Generated function for speciation event likelihoods for flow algorithm.
"""
@generated function λevent_full(t   ::Float64,
                                llik::Array{Array{Float64,1},1},
                                ud1 ::Array{Float64,1},
                                ud2 ::Array{Float64,1},
                                p   ::Array{Float64,1},
                                pr  ::Int64,
                                λts ::Array{Float64,1},
                                r   ::Array{Float64,1},
                                af! ::Function,
                                ::Val{k},
                                ::Val{h},
                                ::Val{ny},
                                ::Val{model}) where {k, h, ny, model}

  eqs = quote end
  popfirst!(eqs.args)

  S = create_states(k, h)

  # if speciation model
  if model[1]

    # last non β parameters
    bbase = h*(k+1) + 2*k*h + k*(k-1)*h + h*(h-1)

    # ncov
    ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

    # y per parameter
    yppar = isone(ny) ? 1 : div(ny,ncov)

    # add environmental function
    push!(eqs.args, :(af!(t, r)))

    isone(ny) && push!(eqs.args, :(r1 = r[1]))

    for j = Base.OneTo(h), i = Base.OneTo(k)
      coex = Expr(:call, :+)
      for yi in Base.OneTo(yppar)
        rex = isone(ny) ? :(r1) : :(r[$(yi+yppar*(i-1))])
        push!(coex.args, :(p[$(bbase + yi + yppar*((i-1) + k*(j-1)))] * $rex))
      end
      push!(eqs.args, :(λts[$(i + k*(j-1))] = p[$(i + (k+1)*(j-1))]*exp($coex)))
    end

    # likelihood for individual areas states
    for j = Base.OneTo(h), i = Base.OneTo(k)
      push!(eqs.args,
          :(llik[pr][$(i+(2^k-1)*(j-1))] =
            ud1[$(i+(2^k-1)*(j-1))] * ud2[$(i+(2^k-1)*(j-1))] * λts[$(i+k*(j-1))]))
    end

    # likelihood for widespread states
    for j = Base.OneTo(h), i = k+1:2^k-1
      s = S[i + (2^k-1)*(j-1)]

      # within region speciation
      ex = Expr(:call, :+)
      for a in s.g
        push!(ex.args, :((ud1[$(a + s.h*(2^k-1))]   * ud2[$(i + (2^k-1)*(j-1))] +
                          ud1[$(i + (2^k-1)*(j-1))] * ud2[$(a + s.h*(2^k-1))]) *
                          0.5 * λts[$(a + s.h*(k))]))
      end

      # between region speciation
      for (la, ra) = vicsubsets(s.g)[1:div(end,2)]
        push!(ex.args,
         :((ud1[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] +
            ud1[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)]) *
            0.5 * p[$((k+1)*(1+s.h))]))
      end

      push!(eqs.args, :(llik[pr][$(i + (2^k-1)*(j-1))] = $ex))
    end

  # *not* speciation model
  else

    # likelihood for individual areas states
    for j = Base.OneTo(h), i = Base.OneTo(k)
      push!(eqs.args,
          :(llik[pr][$(i+(2^k-1)*(j-1))] =
            ud1[$(i+(2^k-1)*(j-1))] * ud2[$(i+(2^k-1)*(j-1))] * p[$(i+(k+1)*(j-1))]))
    end

    # likelihood for widespread states
    for j = Base.OneTo(h), i = k+1:2^k-1
      s = S[i + (2^k-1)*(j-1)]

      # within region speciation
      ex = Expr(:call, :+)
      for a in s.g
        push!(ex.args, :((ud1[$(a + s.h*(k+1))]*ud2[$(i + (2^k-1)*(j-1))] +
                          ud1[$(i + (2^k-1)*(j-1))]*ud2[$(a + s.h*(k+1))]) *
                          0.5 * p[$(a + s.h*(k+1))]))
      end

      # between region speciation
      for (la, ra) = vicsubsets(s.g)[1:div(end,2)]
        push!(ex.args,
         :((ud1[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] +
            ud1[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h)] *
            ud2[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)]) *
            0.5 * p[$((k+1)*(1+s.h))]))
      end

      push!(eqs.args, :(llik[pr][$(i + (2^k-1)*(j-1))] = $ex))
    end
  end

  #Core.println(eqs)

  return quote
    @inbounds begin
      $eqs
    end
    return nothing
  end
end

"""
    normbysum!(v::Array{Float64,1}, r::Array{Float64,1}, ns::Int64)

Return weights by the sum of all elements of the array to sum 1.
"""
function normbysum!(v::Array{Float64,1}, r::Array{Float64,1}, ns::Int64)
  s = sum(v)
  @inbounds begin
    @simd for i in Base.OneTo(ns)
      r[i] = v[i]/s
    end
  end
end

"""
    normbysum!(v::Array{Float64,1}, ns::Int64)

Return weights by the sum of all elements of the array to sum 1 and overwrites.
"""
function normbysum!(v::Array{Float64,1}, ns::Int64)
  s = sum(v)
  @inbounds begin
    @simd for i in Base.OneTo(ns)
      v[i] /= s
    end
  end
  return s
end

"""
    make_lpf(pupd   ::Array{Int64,1},
             phid   ::Array{Int64,1},
             λpriors::Float64,
             μpriors::Float64,
             gpriors::Float64,
             lpriors::Float64,
             qpriors::Float64,
             βpriors::Tuple{Float64,Float64},
             hpriors::Float64,
             k      ::Int64,
             h      ::Int64,
             model  ::Tuple{Bool,Bool,Bool})

Make log-prior function.
"""
function make_lpf(pupd   ::Array{Int64,1},
                  phid   ::Array{Int64,1},
                  λpriors::Float64,
                  μpriors::Float64,
                  gpriors::Float64,
                  lpriors::Float64,
                  qpriors::Float64,
                  βpriors::Tuple{Float64,Float64},
                  hpriors::Float64,
                  k      ::Int64,
                  h      ::Int64,
                  ny     ::Int64,
                  model  ::Tuple{Bool,Bool,Bool})

  λupds = intersect(1:(h*(k+1)), pupd)
  μupds = intersect((h*(k+1)+1):(h*(k+1)+k*h), pupd)
  gupds = intersect((h*(k+1)+k*h+1):(h*(k+1)+k*h+k*(k-1)*h), pupd)
  lupds = intersect((h*(k+1)+k*h+k*(k-1)*h+1):(h*(k+1)+2k*h+k*(k-1)*h), pupd)
  qupds = intersect((h*(k+1)+2k*h+k*(k-1)*h+1):(h*(k+1)+2k*h+k*(k-1)*h+h*(h-1)),
            pupd)
  bbase = (h*(k+1)+2k*h+k*(k-1)*h+h*(h-1))
  yppar = ny == 1 ? 1 : div(ny,
    model[1]*k +
    model[2]*k +
    model[3]*k*(k-1))

  ncov  = model[1]*k*yppar + model[2]*k*yppar + model[3]*k*(k-1)*yppar
  βupds = intersect((bbase+1):(bbase+ncov), pupd)

  βp_m, βp_v = βpriors

  function f(p ::Array{Float64,1},
             fp::Array{Float64,1})
    lq = 0.0

    # speciation priors
    for i in λupds
      lq += logdexp(p[i], λpriors)
    end

    # global extinction priors
    for i in μupds
      lq += logdexp(p[i], μpriors)
    end

    # area colonization priors
    for i in gupds
      lq += logdexp(p[i], gpriors)
    end

    # area loss priors
    for i in lupds
      lq += logdexp(p[i], lpriors)
    end

    # hidden states transition
    for i in qupds
      lq += logdexp(p[i], qpriors)
    end

    # betas
    for i in βupds
      lq += logdnorm(p[i], βp_m, βp_v)
    end

    # hidden states factors
    for i in phid
      lq += logdexp(fp[i], hpriors)
    end

    return lq
  end

  return f
end

"""
    make_assign_hidfacs(::Val{k},
                        ::Val{h},
                        ::Val{ny},
                        ::Val{model}) where {k, h, ny, model})

Generated function to assign factors to parameters given
factor+parameter `fp` vector.
"""
function make_assign_hidfacs(::Val{k},
                             ::Val{h},
                             ::Val{ny},
                             ::Val{model}) where {k, h, ny, model}

  assign_hidfacs! = (p::Array{Float64,1}, fp::Array{Float64,1}) ->
    assign_hidfacs_full(p, fp, Val(k), Val(h), Val(ny), Val(model))

  return assign_hidfacs!
end

"""
    assign_hidfacs(p ::Array{Float64,1},
                   fp::Array{Float64,1},
                   ::Val{k},
                   ::Val{h},
                   ::Val{ny},
                   ::Val{model}) where {k, h, ny, model}

Generated function to assign factors to parameters given
factor+parameter `fp` vector.
"""
@generated function assign_hidfacs_full(p ::Array{Float64,1},
                                        fp::Array{Float64,1},
                                        ::Val{k},
                                        ::Val{h},
                                        ::Val{ny},
                                        ::Val{model}) where {k, h, ny, model}

  # number of covariates
  yppar = ny == 1 ? 1 : div(ny,
    model[1]*k +
    model[2]*k +
    model[3]*k*(k-1))

  # starting indices for models 2 and 3
  m2s = model[1]*k*h*yppar
  m3s = m2s + model[2]*k*h*yppar

  ex = quote end
  pop!(ex.args)

  for j in 1:(h-1)
    for i in 1:k

      # speciation
      push!(ex.args, :(p[$((k+1)*j + i)] =
        p[$((k+1)*(j-1) + i)] + fp[$((k+1)*j + i)]))

      # extinction
      s = (k+1)*h
      push!(ex.args,
        :(p[$(s + k*j + i)] =
            p[$(s + k*(j-1) + i)] + fp[$(s + k*j + i)]))

      # gain
      s = (2k+1)*h
      for a in 1:(k-1)
        push!(ex.args,
          :(p[$(s + k*(k-1)*j + a + (k-1)*(i-1))] =
              p[$(s + k*(k-1)*(j-1) + a + (k-1)*(i-1))] +
              fp[$(s + k*(k-1)*j + a + (k-1)*(i-1))]))
      end

      # loss
      s = (2k+1 + k*(k-1))*h
      push!(ex.args,
        :(p[$(s + k*j + i)] =
            p[$(s + k*(j-1) + i)] + fp[$(s + k*j + i)]))

      # betas
      s = (3k+1+k*(k-1))*h + h*(h-1)
      if model[1]
        for l = Base.OneTo(yppar)
          push!(ex.args,
            :(p[$(s + l + yppar*(i-1) + k*j*yppar)] =
              p[$(s + l + yppar*(i-1) + k*(j-1)*yppar)] +
              fp[$(s + l + yppar*(i-1) + k*j*yppar)]))
        end
      end
      if model[2]
        for l = Base.OneTo(yppar)
          push!(ex.args,
            :(p[$(s + m2s + l + yppar*(i-1) + k*j*yppar)] =
              p[$(s + m2s + l + yppar*(i-1) + k*(j-1)*yppar)] +
              fp[$(s + m2s + l + yppar*(i-1) + k*j*yppar)]))
        end
      end
      if model[3]
        for a = 1:(k-1), l = Base.OneTo(yppar)
          push!(ex.args,
            :(p[$(s + m3s + l + (a-1)*yypar + yppar*(i-1) + k*j*yppar)] =
              p[$(s + m3s + l + (a-1)*yypar + yppar*(i-1) + k*(j-1)*yppar)] +
              fp[$(s + m3s + l + (a-1)*yypar + yppar*(i-1) + k*j*yppar)]))
        end
      end
    end
    # between-region speciation
    push!(ex.args, :(p[$((k+1)*j + (k+1))] =
      (p[$((k+1)*(j-1) + (k+1))] + fp[$((k+1)*j + (k+1))])))
  end

  #Core.println(ex)

  return quote
    @inbounds begin
      $ex
    end
  end
end

"""
    check_negs(x::Array{Float64,1}, k::Int64)

Return `true` if any in x is less or equal than 0.
"""
function check_negs(x::Array{Float64,1}, k::Int64)
  @inbounds begin
    @simd for i in Base.OneTo(k)
      if x[i] < 0.0
        return true
      end
    end
    return false
  end
end

"""
    make_λevent(h    ::Int64,
                k    ::Int64,
                ny   ::Int64,
                flow ::Bool,
                model::NTuple{3,Bool},
                af!  ::Function)

Make function for λevent.
"""
function make_λevent(h    ::Int64,
                     k    ::Int64,
                     ny   ::Int64,
                     flow ::Bool,
                     model::NTuple{3,Bool},
                     af!  ::Function)

  λts = Array{Float64,1}(undef,h*k)
  r   = Array{Float64,1}(undef,ny)

  # make speciation events and closure
  if flow
    λevent! = (t   ::Float64,
               llik::Array{Array{Float64,1},1},
               ud1 ::Array{Float64,1},
               ud2 ::Array{Float64,1},
               p   ::Array{Float64,1},
               pr  ::Int64) ->
    begin
      λevent_full(t, llik, ud1, ud2, p, pr, λts, r, af!,
        Val(k), Val(h), Val(ny), Val(model))
      return nothing
    end

  else
    λevent! = (t   ::Float64,
               llik::Array{Float64,1},
               ud1 ::Array{Float64,1},
               ud2 ::Array{Float64,1},
               p   ::Array{Float64,1}) ->
    begin
      λevent_full(t, llik, ud1, ud2, p, λts, r, af!,
        Val(k), Val(h), Val(ny), Val(model))
      return nothing
    end
  end

  return λevent!
end

"""
    make_rootll(h    ::Int64,
                k    ::Int64,
                ny   ::Int64,
                model::NTuple{3,Bool},
                af!  ::Function)

Make root conditioning likelihood function.
"""
function make_rootll(h    ::Int64,
                     k    ::Int64,
                     ny   ::Int64,
                     model::NTuple{3,Bool},
                     af!  ::Function)

  λts = Array{Float64,1}(undef,h*k)
  r   = Array{Float64,1}(undef,ny)

  rootll = (t   ::Float64,
            llik::Array{Float64,1},
            extp::Array{Float64,1},
            w   ::Array{Float64,1},
            p   ::Array{Float64,1}) ->
    begin
      rootll_full(t, llik, extp, w, p, λts, r, af!,
        Val(k), Val(h), Val(ny), Val(model))::Float64
    end

  return rootll
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
module C99
include("Utls.jl")
using Statistics
using .Utls

"""
    C99.SegmentObject(window_size, similarity_matrix, rank_matrix, sum_matrix, std_coeff, tokenizer)
C99 is a method for determining segment boundaries through divisive clustering.
# Arguments
- `window_size`: window_size is used to create a rank matrix and specifies the range of adjacent sentences to be referenced.
- `similarity_matrix`: Matrix of calculated cosine similarity between sentences.
- `rank_matrix`: Each value in the similarity matrix is replaced by a rank in the local domain. A rank is the number of neighboring elements with lower similarity score.
- `sum_matrix`: Sum of rank matrix in segment regions i to j.
- `std_coeff`: std_coeff is used for the threshold that determines the segment boundary. μ and v are the mean and variance of the gradient δD(n) of the internal density.
- `tokenizer`: Tokenizer for word segmentation.
"""
mutable struct SegmentObject
    window_size::Int
    similarity_matrix::Array{Float64,2}
    rank_matrix::Array{Float64,2}
    sum_matrix::Array{Float64,2}
    std_coeff::Float64
    tokenizer::Any
end

function preprocessing(seg::SegmentObject, document)
    n = length(document)
    @assert n > 0 && length([d for d in document if typeof(d) != String]) == 0
    if n < 3
        return [1] + [0 for _ = 1:n-1]
    end

    seg.window_size = minimum([seg.window_size, n])
    return [Utls.count_elements(seg.tokenizer(document[i])) for i = 1:n]
end

function create_similarity_matrix(seg::SegmentObject, n, preprocessed_document)
    seg.similarity_matrix = zeros(n, n)
    for i = 1:n
        for j = i:n
            seg.similarity_matrix[i, j] = Utls.calculate_cosin_similarity(
                preprocessed_document[i],
                preprocessed_document[j],
            )
            seg.similarity_matrix[j, i] = seg.similarity_matrix[i, j]
        end
    end
end

function create_rank_matrix(seg::SegmentObject, n)
    seg.rank_matrix = zeros(n, n)
    for i = 1:n
        for j = i:n
            r1 = maximum([1, i - seg.window_size + 1])
            r2 = minimum([n, i + seg.window_size - 1])
            c1 = maximum([1, j - seg.window_size + 1])
            c2 = minimum([n, j + seg.window_size - 1])

            s = seg.similarity_matrix[r1:(r2), c1:(c2)]
            k = size(s)
            l = k[1] * k[2]
            sublist = reshape(transpose(s), 1, l)
            lowlist = [x for x in sublist if x < seg.similarity_matrix[i, j]]

            seg.rank_matrix[i, j] = 1.0 * length(lowlist) / ((r2 - r1 + 1) * (c2 - c1 + 1))
            seg.rank_matrix[j, i] = seg.rank_matrix[i, j]
        end
    end
end

function create_sum_matrix(seg::SegmentObject, n)
    seg.sum_matrix = zeros(n, n)
    prefix_sm = zeros(n, n)
    for i = 1:n
        for j = 1:n
            prefix_sm[i, j] = seg.rank_matrix[i, j]
            if i - 1 >= 1
                prefix_sm[i, j] += prefix_sm[i-1, j]
            end
            if j - 1 >= 1
                prefix_sm[i, j] += prefix_sm[i, j-1]
            end
            if i - 1 >= 1 && j - 1 >= 1
                prefix_sm[i, j] -= prefix_sm[i-1, j-1]
            end
        end
    end
    for i = 1:n
        for j = i:n
            if i == 1
                seg.sum_matrix[i, j] = prefix_sm[j, j]
            else
                seg.sum_matrix[i, j] =
                    prefix_sm[j, j] - prefix_sm[i-1, j] - prefix_sm[j, i-1] +
                    prefix_sm[i-1, i-1]
            end
            seg.sum_matrix[j, i] = seg.sum_matrix[i, j]
        end
    end
end

function determine_boundaries(seg::SegmentObject, n)
    D = 1.0 * seg.sum_matrix[1, n] / (n * n)
    g = init_region(1, n, seg)
    darr, region_arr, idx = [D], [g], []
    sum_region, sum_area = float(seg.sum_matrix[1, n]), float(n * n)
    for i = 1:n-1
        mx, pos = -1e9, -1

        for (j, region) in enumerate(region_arr)
            if region.l == region.r
                continue
            end

            split_region(region, seg)
            den = sum_area - region.area + region.lch.area + region.rch.area
            cur = (sum_region - region.tot + region.lch.tot + region.rch.tot) / den
            if cur > mx
                mx, pos = cur, j
            end
        end

        @assert pos >= 1
        tmp = region_arr[pos]
        region_arr[pos] = tmp.rch
        insert!(region_arr, pos, tmp.lch)
        sum_region += tmp.lch.tot + tmp.rch.tot - tmp.tot
        sum_area += tmp.lch.area + tmp.rch.area - tmp.area
        append!(darr, sum_region / sum_area)
        append!(idx, tmp.best_pos)
    end

    dgrad = [(darr[i+1] - darr[i]) for i = 1:length(darr)-1]
    cutoff_threshold = mean(dgrad) + seg.std_coeff * std(dgrad)

    @assert length(idx) == length(dgrad)
    above_cutoff_idx = [i for i in 1:length(dgrad) if dgrad[i] >= cutoff_threshold]

    boundaries = ["0" for _ = 1:n]
    if length(above_cutoff_idx) != 0
        for i in idx[1:maximum(above_cutoff_idx)]
            boundaries[i] = "1"
        end
    end
    return join(boundaries[1:end-1])
end

mutable struct Region
    tot::Float64
    l::Int
    r::Int
    area::Float64
    lch::Any
    rch::Any
    best_pos::Int
end

function init_region(l, r, seg::SegmentObject)
    @assert r >= l
    g = Region(seg.sum_matrix[l, r], l, r, 0, nothing, nothing, -1)
    g.area = (r - l + 1)^2
    return g
end

function split_region(g::Region, seg::SegmentObject)
    if g.best_pos >= 1
        return
    end
    if g.l == g.r
        g.best_pos = g.l
        return
    end

    @assert g.r > g.l
    mx, pos = -1e9, -1
    for i = g.l:g.r-1
        carea = (i - g.l + 1)^2 + (g.r - i)^2
        cur = (seg.sum_matrix[g.l, i] + seg.sum_matrix[i+1, g.r]) / carea
        if cur > mx
            mx, pos = cur, i
        end
    end
    @assert pos >= g.l && pos < g.r
    g.lch = init_region(g.l, pos, seg)
    g.rch = init_region(pos + 1, g.r, seg)
    g.best_pos = pos

end

"""
    C99.segment(seg, document, n) -> String
Performs the splitting of the document entered in the `document` argument.
# Arguments
- `seg`: Segment object.
- `document`: The document to be text segmented.
- `n`: Document Length.

# Examples
```jldoctest
using TextSegmentation

n = length(document)
init_matrix = zeros(n, n)
window_size = 2
std_coeff = 1.2

c99 = C99.SegmentObject(window_size, init_matrix, init_matrix, init_matrix, std_coeff, Utls.tokenize)
result = C99.segment(c99, document, n)
println(result)
000100010000
```
"""
function segment(seg, document, n)
    preprocessed_document = preprocessing(seg, document)
    create_similarity_matrix(seg, n, preprocessed_document)
    create_rank_matrix(seg, n)
    create_sum_matrix(seg, n)
    return determine_boundaries(seg, n)
end
end