text:
using HorizonSideRobots

function main(r::Robot)
    num_vert = moves!(r, Sud)
    num_hor = moves!(r, West)

    paint()

    vygol()
    home(num_vert,num_hor)
end

function home(num_vert,num_hor)
    moves!(r, Nord, num_vert)
    moves!(r, Ost, num_hor)
end

function paint()
    i = 0
    while isborder(r,Nord)==false

        moves!(r,Ost)

        for i in 1:i
            move!(r,West)
        end
        putmarker!(r)
        putmarkers!(r,West)

        move!(r,Nord)
        putmarker!(r)
        i = i + 1
    end
end

function vygol()
    moves!(r, Sud)
    moves!(r, West)
end

function putmarkers!(r::Robot, side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end

function moves!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps # символ "_" заменяет фактически не используемую переменную
        move!(r,side)
    end
end

r=Robot(5,5,animate=true)

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
module Examples
export UnlabeledModule

using ..Theories

UnlabeledModule = UnlabeledPresentation(
  [(),()],
  [Operation([1,1],1), Operation([1,2],2)],
  [UnlabeledEquation(
    UnlabeledContext([1,1,2]),
    Appl(2, [Appl(1, [Var(1), Var(2)]), Var(3)]),
    Appl(2, [Var(1), Appl(2, [Var(2), Var(3)])])
  )]
)

end

ext: jl


text:
"""
Use a circular buffer to contain used RuleNodes to be recycled.
First check the recycle bin to see if there are available RuleNodes, if not allocated a new one.  This helps minimize memory allocations and improves performance
"""
const NodeRecycler = CircularBuffer{RuleNode}

RuleNode(bin::Nothing, ind::Int) = RuleNode(ind)
RuleNode(bin::Nothing, ind::Int, children::Vector{RuleNode}) = RuleNode(ind, children)
RuleNode(bin::Nothing, ind::Int, _val::Any) = RuleNode(ind, _val)
RuleNode(bin::Nothing, ind::Int, _val::Any, children::Vector{RuleNode}) = RuleNode(ind, _val, children)
function RuleNode(bin::NodeRecycler, ind::Int)
    isempty(bin) && return RuleNode(ind)
    node = pop!(bin)
    node.ind, node._val = ind, nothing
    empty!(node.children)
    return node
end
function RuleNode(bin::NodeRecycler, ind::Int, children::Vector{RuleNode})
    isempty(bin) && return RuleNode(ind, children)
    node = pop!(bin)
    node.ind, node._val, node.children = ind, nothing, children
    return node
end
function RuleNode(bin::NodeRecycler, ind::Int, _val::Any)
    isempty(bin) && return RuleNode(ind, _val)
    node = pop!(bin)
    node.ind, node._val = ind, _val
    empty!(node.children)
    return node
end
function RuleNode(bin::NodeRecycler, ind::Int, _val::Any, children::Vector{RuleNode})
    isempty(bin) && return RuleNode(ind, _val, children)
    node = pop!(bin)
    node.ind, node._val, node.children = ind, _val, children
    return node
end

Base.deepcopy(bin::Nothing, root::RuleNode) = deepcopy(root)
function Base.deepcopy(bin::NodeRecycler, root::RuleNode)
    root2 = RuleNode(bin, root.ind, root._val)
    for child in root.children
        push!(root2.children, deepcopy(bin, child))
    end
    root2
end

recycle!(bin::Nothing, root::RuleNode) = root
function recycle!(bin::NodeRecycler, root::RuleNode)
    if isfull(bin)
        empty!(root.children)
    else
        push!(bin, root)
        while !isempty(root.children)
            child = pop!(root.children)
            recycle!(bin, child)
        end
    end
end


text:
for T in [History, QHistory]
    @testset "$(T): Basic functions" begin
        _history = T(Float64)

        @test push!(_history, 1, 10.) == Float64(10.)
        @test push!(_history, 2, Float64(21.)) == Float64(21.)

        @test_throws ArgumentError push!(_history, 1, Float64(11.))
        @test_throws ArgumentError push!(_history, 2, 10)
        @test_throws ArgumentError push!(_history, 3, Float32(10))

        _history = T(Float64)

        numbers = collect(1:100)
        push!(_history, 10, Float64(5))
        for i = numbers
            @test push!(_history, Float64(i + 1)) == Float64(i + 1)
        end

        @test first(_history) == (10, Float64(5.))
        @test last(_history) == (110, Float64(101.))

        _history = T(Float64)

        numbers = collect(1:2:200)
        for i = numbers
            @test push!(_history, i, Float64(i + 1)) == Float64(i + 1)
        end

        println(_history)
        show(_history); println()

        @test first(_history) == (1, Float64(2.))
        @test last(_history) == (199, Float64(200.))

        for (i, v) in enumerate(_history)
            @test in(i, numbers)
            @test Float64(i + 1) == v
        end

        a1, a2 = get(_history)
        @test typeof(a1) <: Vector{Int} && typeof(a2) <: Vector{Float64}
        @test length(a1) == length(a2) == length(numbers) == length(_history)
        @test convert(Vector{Float64}, a1 .+ 1) == a2
    end

    @testset "$(T): No explicit iteration" begin
        _history = T(Float64)
    end

    @testset "$(T): Storing arbitrary types" begin
        _history = T(String, UInt8)

        for i = 1:100
            @test push!(_history, i % UInt8, string("i=", i + 1)) == string("i=", i+1)
        end

        show(_history); println()

        a1, a2 = get(_history)
        @test typeof(a1) <: Vector{UInt8}
        @test typeof(a2) <: Vector{String}
    end
end

@testset "History: increment!" begin
    _history = History(Float64)
    val = 1.
    @test increment!(_history, 0, val) == val
    @test increment!(_history, 0, val) == 2val
    @test increment!(_history, 2, 4val) == 4val
    @test increment!(_history, 10, 5val) == 5val
    _history2 = QHistory(Float64)
    @test_throws MethodError increment!(_history2, 1, val) == val
end


text:
using TargetFlow
using Test

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


text:
export RandomKooshballTrajectory, randomKooshballNodes, randomKooshballDensity

"""
    KooshballTrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , kargs...)

returns a 3d kooshball trajectory with randomized order of the profiles

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
"""
function RandomKooshballTrajectory(::Type{T}, numProfiles, numSamplingPerProfile
                  ; TE=0.0
                  , AQ=1.e-3
                  , kargs...) where T
  nodes = randomKooshballNodes(T, numProfiles, numSamplingPerProfile; numSlices=1)
  times = kooshballTimes(T, numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return  Trajectory("RandomKooshball", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, false, true)
end

function randomKooshballNodes(::Type{T}, numProfiles, numSamplingPerProfile; kargs...) where T
  nodes = zeros(T, 3, numSamplingPerProfile, numProfiles)
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  profileList = shuffle_vector([i for i=1:numProfiles])

  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      nodes[1,k,l] = (-1)^l * pos[k]*xax(profileList[l],2*numProfiles)
      nodes[2,k,l] = (-1)^l * pos[k]*yax(profileList[l],2*numProfiles)
      nodes[3,k,l] = (-1)^l * pos[k]*zax(profileList[l],2*numProfiles)
    end
  end
  return reshape(nodes, 3, numSamplingPerProfile*numProfiles)
end

function randomKooshballDensity(::Type{T}, numProfiles::Int64, numSamplingPerProfile::Int64) where T
  density = zeros(T, numSamplingPerProfile, numProfiles)
  pos = collect((0:numSamplingPerProfile-1)/numSamplingPerProfile .- 0.5)
  for l = 1:numProfiles
    for k = 1:numSamplingPerProfile
      density[k,l] = abs(pos[k])
    end
  end
  return vec(density)
end


text:
module TrafficAssignment

include("TAFW/TAFW.jl")
include("TAPAS/TAPAS.jl")

using .TAFW
using .TAPAS
using DataFrames
using Plots

"""
    assigntraffic(network; method=:FW, assignment=:UE, tol=1e-5, maxiters=20, maxruntime=300, log=:off)

Fukushima Frank-Wolfe method for traffic assignment.

# Returns
a named tuple with keys `:metadata`, `:report`, and `:output`
- `metadata::String`  : Text defining the traffic assignment run
- `report::DataFrame` : A log of total network flow, total network cost, and run time for every iteration
- `output::DataFrame` : Flow and cost for every arc from the final iteration

# Arguments
- `network::String`         : Network
- `method::Symbol=:FW`      : One of `:FW`, `:fukushimaFW`, `:conjugateFW`
- `assignment::Symbol=:UE`  : Assignment type; one of `:UE`, `:SO`
- `tol::Float64=1e-5`       : Tolerance level for relative gap
- `maxiters::Int64=20`      : Maximum number of iterations
- `maxruntime::Int64=300`   : Maximum algorithm run time (seconds)
"""
function assigntraffic(network; method=:FW, assignment=:UE, tol=1e-5, maxiters=20, maxruntime=300, log=:off)
    if method == :FW
        G = TAFW.build(network, assignment)
        return FW(G, tol, maxiters, maxruntime, log)
    elseif method == :fukushimaFW
        G = TAFW.build(network, assignment)
        return fukushimaFW(G, tol, maxiters, maxruntime, log)
    elseif method == :conjugateFW
        G = TAFW.build(network, assignment)
        return conjugateFW(G, tol, maxiters, maxruntime, log)
    elseif method == :TAPAS
        G = TAPAS.build(network, assignment)
        return itapas(G, tol, maxiters, maxruntime, log)
    else return (metadata="", report=DataFrame(), solution=DataFrame())
    end
end

"""
    compare(network; methods, assignment=:UE, tol=1e-5, maxiters=20, maxruntime=300, log=:on)

Compare assignment methods.
"""
function compare(network; methods, assignment=:UE, tol=1e-5, maxiters=20, maxruntime=300, log=:on)
    fig = plot()

    for method in methods
        _, report, = assigntraffic(network; method=method, assignment=assignment, tol=tol, maxiters=maxiters, maxruntime=maxruntime, log=log)
        y = report[!,:LOG₁₀RG]
        x = 0:length(y)-1
        plot!(x,y, label=String(method))
    end
    display(fig)
end

"""
    getrg(network, solution::DataFrame)

Returns relative gap for `network` given traffic assignment `solution`.
"""
function getrg(network, solution::DataFrame)
    build, djk, path = TAFW.build, TAFW.djk, TAFW.path  # NOTE: Can chose any assignment method tools to calculate relative gap
    G = build(network, :none)
    N, A, O, K = G.N, G.A, G.O, G.K

    for row in 1:nrow(solution)
        i = solution[row, :FROM]::Int64
        j = solution[row, :TO]::Int64
        x = solution[row, :FLOW]
        c = solution[row, :COST]
        k = K[i,j]
        a = A[k]
        a.x = x
        a.c = c
    end

    num = 0.0
    for (i,o) in enumerate(O)
        r = o.n
        Lᵖ = djk(G, o)
        for (j,s) in enumerate(o.S)
            qᵣₛ = o.Q[j]
            pᵣₛ = path(G, Lᵖ, r, s)
            for a in pᵣₛ num += qᵣₛ * a.c end
        end
    end

    den = 0.0
    for (k,a) in enumerate(A) den += a.x * a.c end

    rg = 1 - num/den
    return rg
end

export assigntraffic

end
# TODO: Save output from traffic assignment solution


text:
module MCBN

    import Base: show, length
    using DAI
    import DAI: clearBackups!
    using Distributions
    using Graphs
    import Graphs: add_edge!

    using StatsBase
    import StatsBase: sample
    #using Graphs
    #import Graphs: vertex_index, make_vertex, add_edge!

    using SAMC
    import SAMC: propose!, energy, reject!, save!, record, Sampler

    include("bayesnet.jl")

end # module


text:
function try_move!(robot, side)::Bool
    if isborder(robot, side)
        false
    else
        move!(robot, side)
        true
    end
end

function kross_n_na_n(r::Robot, n::Int)
    while !(isborder(r, Nord))

        first(r)

        second(r)

    end
end

function first(r)
    for i in 1:n
        while !(isborder(r, Ost))
            for i in 1:n
                putmarker!(r)
                try_move!(r, Ost)
            end

            for i in 1:n
                try_move!(r, Ost)
            end
        end

        while !(isborder(r, West))
            move!(r, West)
        end
        try_move!(r, Nord)
    end
end

function second(r)
    for i in 1:n
        while !(isborder(r, Ost))
            for i in 1:n
                try_move!(r, Ost)
            end

            for i in 1:n
                putmarker!(r)
                try_move!(r, Ost)
            end

        end

        while !(isborder(r, West))
            move!(r, West)
        end
        try_move!(r, Nord)
    end
end

kross_n_na_n(r, 3)


text:
using TopOptMakie
using Test

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


text:
using DomainSets: convert_numtype, convert_prectype,
    promote_numtype, promote_prectype

function test_dimension()
    @test DomainSets.euclideandimension(Int) == 1
    @test DomainSets.euclideandimension(Float64) == 1
    @test DomainSets.euclideandimension(ComplexF64) == 1
    @test DomainSets.euclideandimension(SVector{2,Float64}) == 2
    @test DomainSets.euclideandimension(MVector{2,Float64}) == 2
    @test DomainSets.euclideandimension(Tuple{Int,Int}) == 2
    @test DomainSets.euclideandimension(Tuple{Int,Float64}) == 2
    @test_throws MethodError DomainSets.euclideandimension(Vector{Float64})
end

function test_components()
    @test components([1,2,3]) == ()
    d = UnionDomain(Point(1),Point(2))
    @test iscomposite(d)
    @test ncomponents(d) == length(components(d))
end

function test_prectype()
    @test prectype(1.0) == Float64
    @test prectype(big(1.0)) == BigFloat
    @test prectype(1) == typeof(float(1))
    @test prectype(SVector(1,2)) == typeof(float(1))
    @test prectype(SVector(1,big(2))) == typeof(float(big(2)))
    @test prectype(1.0+2.0im) == Float64
    @test prectype([1.0+2.0im, 3.0]) == Float64
    @test prectype(NTuple{2,Int}) == Float64
    @test prectype((1.0,)) == Float64
    @test prectype((1.0, 2.0)) == Float64
    @test prectype((1.0, 2.0, 3.0)) == Float64
    @test prectype((1.0, big(2.0), 3.0+im)) == BigFloat
    @test prectype(NTuple{4,Int}) == Float64
    @test @inferred(prectype(1, 2.0)) == Float64
    @test @inferred(prectype((1, 2.0, 3, 40+im))) == Float64

    @test convert_prectype(2, Float64) == 2
    @test convert_prectype(2, Float64) isa Float64
    @test convert_prectype(1.0+im, BigFloat) == 1+im
    @test convert_prectype(1.0+im, BigFloat) isa Complex{BigFloat}
    @test convert_prectype(SA[1,2], Float64) == SA[1.0,2.0]
    @test convert_prectype(SA[1,2], Float64) isa SVector{2,Float64}
    @test convert_prectype(SA[1,2+im], BigFloat) isa SVector{2,Complex{BigFloat}}
    @test_throws ErrorException convert_prectype("a", BigFloat)

    @test promote_prectype(2) == 2
    @test promote_prectype(2, 3.0) isa Tuple{Float64,Float64}
    @test promote_prectype(2, 3.0+im, big(4)) isa Tuple{BigFloat,Complex{BigFloat},BigFloat}
end

function test_numtype()
    @test numtype(1.0) == Float64
    @test numtype(big(1.0)) == BigFloat
    @test numtype(1) == Int
    @test numtype([1,2,3], [4,5,6]) == Int
    @test numtype(SVector(1,2)) == Int
    @test numtype(SVector(1,big(2))) == BigInt
    @test numtype(Array{Float64,2}) == Float64
    @test numtype(1.0+2.0im) == Complex{Float64}
    @test numtype([1.0+2.0im, 3.0]) == Complex{Float64}
    @test numtype(NTuple{2,Complex{Int}}) == Complex{Int}
    @test numtype((1.0,)) == Float64
    @test numtype((1.0, 2.0)) == Float64
    @test numtype((1.0, 2.0, 3.0)) == Float64
    @test numtype((1.0, 2.0, 3.0, 4.0)) == Float64
    @test numtype(1.0, big(2.0), 3.0+im) == Complex{BigFloat}
    @test numtype((1.0, big(2.0), 3.0+im)) == Complex{BigFloat}
    @test @inferred(numtype(1, 2.0)) == Float64
    @test @inferred(numtype((1, 2.0, 3, 40+im))) == Complex{Float64}

    @test convert_numtype(2, Float64) == 2
    @test convert_numtype(2, Float64) isa Float64
    @test convert_numtype(SA[1,2], Float64) == SA[1,2]
    @test convert_numtype(SA[1,2], Float64) isa SVector{2,Float64}
    @test_throws ErrorException convert_numtype("a", BigFloat)

    @test promote_numtype(2) == 2
    @test promote_numtype(2, 3.0) isa Tuple{Float64,Float64}
    @test promote_numtype(2, 3.0+im, big(4)) isa Tuple{Complex{BigFloat},Complex{BigFloat},Complex{BigFloat}}
end

@testset "common functionality" begin
    @testset "dimension" begin
        test_dimension()
    end
    @testset "elements" begin
        test_components()
    end
    @testset "prectype" begin
        test_prectype()
    end
    @testset "numtype" begin
        test_numtype()
    end
end


text:
export Poisson

"""
Description:

    Poisson factor node

    Real scalars
    l > 0 (rate)

    f(out, l) = Poisson(out|l) = 1/(x!) * l^x * exp(-l)

Interfaces:

    1. out
    2. l

Construction:

    Poisson(id=:some_id)
"""
mutable struct Poisson <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Poisson(out, l; id=generateId(Poisson))
        @ensureVariables(out, l)
        self = new(id, Array{Interface}(undef, 2), Dict{Symbol, Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:l] = self.interfaces[2] = associate!(Interface(self), l)
        return self
    end
end

slug(::Type{Poisson}) = "Poisson"

format(dist::ProbabilityDistribution{Univariate, Poisson}) = "$(slug(Poisson))(l=$(format(dist.params[:l])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{Poisson}; l=1.0) = ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>l))
ProbabilityDistribution(::Type{Poisson}; l=1.0) = ProbabilityDistribution{Univariate, Poisson}(Dict(:l=>l))

dims(dist::ProbabilityDistribution{Univariate, Poisson}) = 1

vague(::Type{Poisson}) = ProbabilityDistribution(Univariate, Poisson, l=huge)

isProper(dist::ProbabilityDistribution{Univariate, Poisson}) = (0 < dist.params[:l] < huge)

unsafeMean(dist::ProbabilityDistribution{Univariate, Poisson}) = Float64(dist.params[:l])

unsafeVar(dist::ProbabilityDistribution{Univariate, Poisson}) = Float64(dist.params[:l])

logPdf(dist::ProbabilityDistribution{Univariate, Poisson}, x) = x*log(dist.params[:l]) - dist.params[:l] - logfactorial(x)

sample(dist::ProbabilityDistribution{Univariate, Poisson}) = poisinvcdf(dist.params[:l], rand())

# ∑ [λ^k*log(k!)]/k! from k=0 to inf
# Approximates the above sum for calculation of averageEnergy and differentialEntropy
# @ref https://arxiv.org/pdf/1708.06394.pdf
function apprSum(l, j=100)
    sum([(l)^(k)*logfactorial(k)/exp(logfactorial(k)) for k in collect(0:j)])
end

# Entropy functional
# @ref https://en.wikipedia.org/wiki/Poisson_distribution
function differentialEntropy(dist::ProbabilityDistribution{Univariate, Poisson})
    l = clamp(dist.params[:l], tiny, huge)
    l*(1-log(l)) + exp(-l)*apprSum(l)
end

# Average energy functional
function averageEnergy(::Type{Poisson}, marg_out::ProbabilityDistribution{Univariate}, marg_l::ProbabilityDistribution{Univariate})
    unsafeMean(marg_l) -
    unsafeMean(marg_out)*unsafeLogMean(marg_l) +
    exp(-unsafeMean(marg_out))*apprSum(unsafeMean(marg_out))
end


text:
using LowRankModels, FullRankModels, Compat

srand(1)
# parameters
m, n, k = 20, 20, 3
lambda = 10

# generate data
X0 = randn(m,k)
Y0 = randn(k,n)
A = X0*Y0 + .1*randn(m,n)
P = sparse(float(randn(m,n) .>= .5)) # observed entries

I,J = findn(P) # observed indices (vectors)
@compat obs = Tuple{Int,Int}[(I[a],J[a]) for a = 1:length(I)]
observed_features, observed_examples = sort_observations(obs,size(P)...)

losses = fill(QuadLoss(1), n)
reg = MaxNormReg(lambda)
gfrm = GFRM(A, losses, reg, k, observed_features, observed_examples,
	zeros(m,n), zeros(m+n,m+n))

fit!(gfrm, PrismaParams(maxiter=100))

# compare with SDP solver
using Convex

W = Variable(m+n, m+n)
obj = sumsquares(P.*(W[1:m, m+1:end]-A)) + lambda*maximum(diag(W))
p = minimize(obj,
	W in :SDP)
solve!(p)

println("mean square error is ", vecnorm(W.value - gfrm.W) / vecnorm(W.value))
W.value = gfrm.W
println("objective of Convex problem evaluated at prisma solution is $(Convex.evaluate(obj))")

# switch out regularizer
gfrm.r = MaxNormReg(1)


text:
"""Helper code for running tests."""
module Testing

using LinearAlgebra: diagind
using TriMatrices
using TriMatrices: TriLayout, nelems, hasdiag
using TriMatrices.Indexing: check_tri_index, car2lin

export LAYOUT_TYPES, LAYOUT_TYPES_P, make_test_matrix, make_test_matrix_pair

const LAYOUT_TYPES = [TriLower, TriUpper, TriSymmetric]
const LAYOUT_TYPES_P = [L{D} for D in [true, false] for L in LAYOUT_TYPES]

default_datafunc(i) = 10i

"""
Make an `n` by `n` `Array` to test TriLayout and TriMatrix methods.

First creates a 1D data vector of the appropriate size for the given layout and
value of `n` and files with values `datafunc(i)` for `i` in `1:nelems(layout, n)`.
The mapping function is to make sure we're not confusing the value with its index.

Then creates an `n` by `n` matrix with the folling values:
* Elements which are in the stored data region for the given `layout` have the
  value of the corresponding element of the data array.
* Elements along the diagonal have value `diag` if the diagonal is not stored
  for the given layout.
* All other elements are zero.

Returns `(data, matrix)`, where `matrix` is the `n` by `n` `Array` and `data`
is the data array for the corresponding `TriMatrix`.
"""
function make_test_matrix(layout::TriLayout, n::Int; diag=-1, datafunc=default_datafunc, T::Type=Int)
	data = collect(T, datafunc(i) for i in 1:nelems(layout, n))

	matrix = zeros(T, n, n)

	for i in 1:n, j in 1:n
		if check_tri_index(Bool, layout, i, j)
			matrix[i, j] = data[car2lin(layout, i, j)]
		elseif i == j
			matrix[i, j] = diag
		end
	end

	return data, matrix
end

"""
Create an `Array` using `make_test_matrix` along with a `TriMatrix` that should
contain the same values.

Returns a `(::Matrix, ::TriMatrix)` tuple.
"""
function make_test_matrix_pair(layout::TriLayout, n::Int; diag=-1, kw...)
	data, arraymat = make_test_matrix(layout, n; diag=diag, kw...)
	trimat = TriMatrix(layout, n, data, diag=diag)
	return arraymat, trimat
end

end  # module


text:
"""
Package: Forecast

    function p(dx, x0)

Return reverse lagged differences of a given order for Vector, Array.

# Arguments
- `dx`: Array or DataFrame of data.
- `x0`: Initial constants the reverse difference. The default value represents an integration of order one and lag one with initial values at zero. The format for the initial values is Array{Real,3}(order, variable, lag)"

# Returns
Lagged differences Vector or Array of a given order.

# Examples
```julia-repl

# Order two with Lag two
julia> x = repeat(1:2,30);
julia> dx = d(x,2,2);
julia> x0 = zeros(2,1,2); # lag 2, 1 variable, order 1
julia> x0[1,:,:] = collect(1:2);
julia> p(dx,x0) ≈ x
true

# Calculation of π
julia> x = 0:0.001:1;
julia> y = sqrt.(1 .- x.^2);
julia> isapprox(4*p(y)[end]/1000 , π, atol = 0.01)
true
```
"""
function p(dx::AbstractArray{<:T},
           x0::AbstractArray{<:T} = reshape(repeat([0],size(dx,2)),1,:,1)) where T<:Real

    format_ol = "x0 format is Array{Real,3}(order, variable, lag)"
    @assert ndims(x0) <= 3 format_ol
    @assert isnothing(findfirst(isnan, dx)) "NaN values not allowed"

    ndims(x0) == 0 && (x0 = reshape([x0],1,1,1))
    ndims(x0) == 1 && (x0 = reshape(x0,:,1,1))
    ndims(x0) == 2 && (x0 = reshape(x0,:,size(x0,2),1))

    or = size(x0,1)
    la = size(x0,3)
    n = size(dx,1)
    nv = size(dx,2)

    @assert nv == size(x0,2) format_ol
    @assert 0 <= la & la <= (n-1)

    px = similar(dx, nv == 1 ? (n+or*la,) : (n+or*la,nv))
    px[1:size(dx,1),:] = dx
    for i in or:-1:1
        for j in 1:la
            pxj  = px[j:la:end,:]
            npxj = size(pxj,1)
            px[j:la:end,:] = cumsum(vcat(x0[i:i,:,j], pxj), dims=1)[1:npxj,:]
        end
    end

    return px

end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ._ProtoBuf_Top_.opentelemetry

mutable struct ExportMetricsServiceRequest <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ExportMetricsServiceRequest(; kwargs...)
        obj = new(meta(ExportMetricsServiceRequest), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ExportMetricsServiceRequest
const __meta_ExportMetricsServiceRequest = Ref{ProtoMeta}()
function meta(::Type{ExportMetricsServiceRequest})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ExportMetricsServiceRequest)
            __meta_ExportMetricsServiceRequest[] =
                target = ProtoMeta(ExportMetricsServiceRequest)
            allflds = Pair{Symbol,Union{Type,String}}[:resource_metrics=>Base.Vector{
                opentelemetry.proto.metrics.v1.ResourceMetrics,
            }]
            meta(
                target,
                ExportMetricsServiceRequest,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_ExportMetricsServiceRequest[]
    end
end
function Base.getproperty(obj::ExportMetricsServiceRequest, name::Symbol)
    if name === :resource_metrics
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{
            opentelemetry.proto.metrics.v1.ResourceMetrics,
        }
    else
        getfield(obj, name)
    end
end

mutable struct ExportMetricsServiceResponse <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ExportMetricsServiceResponse(; kwargs...)
        obj = new(meta(ExportMetricsServiceResponse), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) ||
                error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct ExportMetricsServiceResponse
const __meta_ExportMetricsServiceResponse = Ref{ProtoMeta}()
function meta(::Type{ExportMetricsServiceResponse})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ExportMetricsServiceResponse)
            __meta_ExportMetricsServiceResponse[] =
                target = ProtoMeta(ExportMetricsServiceResponse)
            allflds = Pair{Symbol,Union{Type,String}}[]
            meta(
                target,
                ExportMetricsServiceResponse,
                allflds,
                ProtoBuf.DEF_REQ,
                ProtoBuf.DEF_FNUM,
                ProtoBuf.DEF_VAL,
                ProtoBuf.DEF_PACK,
                ProtoBuf.DEF_WTYPES,
                ProtoBuf.DEF_ONEOFS,
                ProtoBuf.DEF_ONEOF_NAMES,
            )
        end
        __meta_ExportMetricsServiceResponse[]
    end
end

# service methods for MetricsService
const _MetricsService_methods = MethodDescriptor[MethodDescriptor(
    "Export",
    1,
    ExportMetricsServiceRequest,
    ExportMetricsServiceResponse,
)] # const _MetricsService_methods
const _MetricsService_desc = ServiceDescriptor(
    "opentelemetry.proto.collector.metrics.v1.MetricsService",
    1,
    _MetricsService_methods,
)

MetricsService(impl::Module) = ProtoService(_MetricsService_desc, impl)

mutable struct MetricsServiceStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    MetricsServiceStub(channel::ProtoRpcChannel) =
        new(ProtoServiceStub(_MetricsService_desc, channel))
end # mutable struct MetricsServiceStub

mutable struct MetricsServiceBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    MetricsServiceBlockingStub(channel::ProtoRpcChannel) =
        new(ProtoServiceBlockingStub(_MetricsService_desc, channel))
end # mutable struct MetricsServiceBlockingStub

Export(
    stub::MetricsServiceStub,
    controller::ProtoRpcController,
    inp::ExportMetricsServiceRequest,
    done::Function,
) = call_method(stub.impl, _MetricsService_methods[1], controller, inp, done)
Export(
    stub::MetricsServiceBlockingStub,
    controller::ProtoRpcController,
    inp::ExportMetricsServiceRequest,
) = call_method(stub.impl, _MetricsService_methods[1], controller, inp)

export ExportMetricsServiceRequest,
    ExportMetricsServiceResponse,
    MetricsService,
    MetricsServiceStub,
    MetricsServiceBlockingStub,
    Export


text:
#Note: "mu" in the paper is changed to "nu" in this implementation

const numResources = 20
const numCustomers = 30

const costVector = [1.0089874188151076 1.0879646957939293 0.8741122771451293 0.9359882588306315 0.43864501807295264 0.9542393501622234 1.1989859234269766 1.2621107238860385 0.609998805032588 1.1838454086005032 1.2227967100332289 0.6771786786349856 0.9030070554656082 1.008077118707221 0.6223169009297342 1.0141930286652359 0.7858768828913462 0.7378242112490088 1.1780607622142538 1.01807326104081 ]

const lambda_mean = [76.53002063308048 86.14362820943657 122.88025299278662 146.28185560066802 118.36204557525407 116.89212997847898 121.73133969750795 128.81309347095186 109.97936259748218 100.35448635212869 127.50190424998809 144.94138717832348 138.25502814996332 166.98486647696402 124.56609562174265 131.86045631845732 70.14641190233873 129.60121252244954 110.36540167343803 102.23855037080791 117.59720163537895 140.01971460663194 162.9030886752114 154.9291002784065 136.86488610276282 95.8805622206068 108.32708374354893 98.50876348531479 118.2154550673827 110.06932970773386 ]

const lambda_covariance = zeros(Float64,30,30)
lambda_covariance[1,:] = [125.3930185302341 77.92826152166512 63.271552516838874 33.768151878886634 46.30913827673934 72.46552697530642 101.5709508272422 108.81899339018814 74.12717311718666 97.55278941760795 65.29768241991836 60.37030107862086 98.0614050491232 94.41218286271872 101.54561514825632 75.58072261482403 88.63610384039548 49.17764159493406 37.31226701496719 49.54565921368147 69.10077802061265 97.0500535647756 82.45750348696772 86.5209387841706 87.361176411108 73.27376303738711 72.68911190643732 86.23715819167262 109.91317968027055 47.06007262000254 ]
lambda_covariance[2,:] = [77.92826152166512 241.34466425622176 123.26655967215504 135.85002543834102 77.69581014774252 137.9257785083876 154.45397745383215 157.4276920872834 66.17318922028987 65.31718588815153 113.75632863644134 88.74858099382928 156.89320451348104 78.9540303171795 164.55062977418578 139.0990608973585 83.0360390704805 60.56909136421002 109.77173322221971 114.39180865661662 105.10413853892037 106.09606555606173 158.3781620476783 92.33221155079374 105.41721793414318 63.31532179652784 76.59433741247877 90.85223309537729 133.67987394225645 107.55622563843521 ]
lambda_covariance[3,:] = [63.271552516838874 123.26655967215504 194.25900816527476 117.22446490792592 60.1011054327771 108.97503191731747 102.46462681689752 94.81243704481719 57.490712911608824 68.88684943611749 76.92299459365339 69.65781271877995 135.47651349137632 91.09599055638469 113.1514186154446 114.00781512870165 109.83987474584998 97.21693698115662 99.07274938398089 61.91574790884258 48.81055594741975 122.70995140185765 132.45322516874643 101.28991370163848 65.45687802210718 48.37629287396258 84.18744456152609 102.59507772963212 87.44936835719288 109.09527613185146 ]
lambda_covariance[4,:] = [33.768151878886634 135.85002543834102 117.22446490792592 252.86200843965864 84.98560484223748 142.7114661271209 146.8570991536446 140.56768176231733 135.85412306318239 90.19515099461506 112.25623838328013 96.69731091550987 137.58720085380233 103.03646518209119 173.2730472599489 135.05919682078394 78.91279561199606 134.27403593476927 109.3293404093198 105.72570077333346 77.08838723329784 129.9224103074779 174.1724018089613 84.22224814012912 57.19019966874822 105.99995730320177 89.71688985457838 107.6003708015894 75.6194872726754 154.81603581461314 ]
lambda_covariance[5,:] = [46.30913827673934 77.69581014774252 60.1011054327771 84.98560484223748 107.13270854253622 87.73035788501699 62.03150712544537 53.76851339573702 79.98813949708854 45.94289925982661 22.269594231772608 23.73111834490463 88.01798445606941 89.3286003327234 120.10069675893348 69.20440336117244 54.55766332776034 42.88288860679177 70.1129057002168 42.01092858068692 59.219600928930824 42.773303122391596 79.43728980210871 24.671199562469866 74.15155681068951 22.198915997472916 59.51626269015065 77.01111734117143 60.10382898897866 87.78972689961464 ]
lambda_covariance[6,:] = [72.46552697530642 137.9257785083876 108.97503191731747 142.7114661271209 87.73035788501699 196.81156731351402 129.73199192288692 118.51527742609233 99.89185285964962 64.014141543296 65.22265167026764 41.57814323591422 112.30256261364312 107.24449245970433 169.159644032315 136.10043774466817 112.67651491675228 81.34756593989223 83.16130249300898 87.76451912933427 86.93210109231543 116.30359078968738 152.94977668133797 84.05254047271808 113.90544421320413 103.44502430536937 124.10097509282218 141.66643661902717 90.40825977618357 108.94151833990824 ]
lambda_covariance[7,:] = [101.5709508272422 154.45397745383215 102.46462681689752 146.8570991536446 62.03150712544537 129.73199192288692 246.78587212093564 159.816333206852 113.33945936728237 119.33819580151494 150.54483111934005 83.715348894916 135.95042142851236 92.19765014511928 131.39091955202895 155.25184872714766 101.04981875685186 107.21528862617915 102.03579153555103 103.55905941192042 108.67149805022271 154.3494865577637 171.93735904351246 123.44380609818286 140.57366745592802 108.97721149681982 122.32718858409935 122.62272188709112 131.81933536773212 124.05711559790461 ]
lambda_covariance[8,:] = [108.81899339018814 157.4276920872834 94.81243704481719 140.56768176231733 53.76851339573702 118.51527742609233 159.816333206852 229.4000450899269 133.26273697500721 123.18299941763523 125.42277048837666 140.63887890946356 143.14757430008171 116.61926047606192 163.80844617445038 147.79914811484616 120.55387824692136 91.32202290237825 88.63344171509814 90.13164610517698 101.18720431784811 163.04528177603024 163.14962123231697 122.03492070329558 101.19986987123252 119.24451379862802 119.40031359424987 101.84230987849037 143.75332514297773 98.2834085001364 ]
lambda_covariance[9,:] = [74.12717311718666 66.17318922028987 57.490712911608824 135.85412306318239 79.98813949708854 99.89185285964962 113.33945936728237 133.26273697500721 195.7201628818312 108.20812199797342 113.54025482324464 79.58918094245374 130.43701926813216 143.3430654992352 122.74996401214017 96.47363448295651 54.13930709343893 105.98898820013238 75.0752642886231 83.4365814549703 69.475940356554 115.45552301472146 129.86689797715573 86.9572513779304 69.0232589342559 114.9377268730485 109.80174246774094 88.99626244262136 107.23059296721462 118.14790803553956 ]
lambda_covariance[10,:] = [97.55278941760795 65.31718588815153 68.88684943611749 90.19515099461506 45.94289925982661 64.014141543296 119.33819580151494 123.18299941763523 108.20812199797342 166.80086176296027 76.07693572307545 80.26811311703734 113.65728474580696 111.22492107113902 111.74959460831212 111.68415443574753 69.3568004342304 74.84297375124505 73.60853894936437 74.90560279368077 76.82970352508781 93.21202934881335 97.77826918787798 96.56697770478569 64.89053452508826 89.3239441518822 82.2741820294904 82.09140240017425 85.50457130900543 69.98534638939023 ]
lambda_covariance[11,:] = [65.29768241991836 113.75632863644134 76.92299459365339 112.25623838328013 22.269594231772608 65.22265167026764 150.54483111934005 125.42277048837666 113.54025482324464 76.07693572307545 214.94591812150125 68.02605462380691 125.99838732601181 41.10790805593369 61.136046541101024 88.83125425704084 35.46875966095065 93.6317042497511 69.95183665082978 109.04651074800726 51.08594178194177 110.18111699842564 138.5216597929755 112.19098719130103 64.2046627983725 89.67280697294089 66.85389433748807 47.053779930066014 114.76092500933024 79.74429519505637 ]
lambda_covariance[12,:] = [60.37030107862086 88.74858099382928 69.65781271877995 96.69731091550987 23.73111834490463 41.57814323591422 83.715348894916 140.63887890946356 79.58918094245374 80.26811311703734 68.02605462380691 162.4664658162561 96.63820684228517 90.88134548128345 97.29082329382234 96.95704629238494 89.32792035275538 69.75736694717743 47.288804820190386 73.62626632177022 69.0179383069737 109.89371269359617 87.36049937979853 83.71193456408744 35.72524353388308 67.98938858156919 42.66038402822041 57.36348683116531 95.5158096207522 68.59048459183069 ]
lambda_covariance[13,:] = [98.0614050491232 156.89320451348104 135.47651349137632 137.58720085380233 88.01798445606941 112.30256261364312 135.95042142851236 143.14757430008171 130.43701926813216 113.65728474580696 125.99838732601181 96.63820684228517 236.11684095190793 127.79426090582622 156.21861495807775 129.03158590956014 94.84122633028504 112.69942627897169 113.25924452490982 109.08890327663268 109.13719894377545 132.49753995973956 127.83072619794666 120.14618848504415 83.46265897913469 97.50764555651743 95.35720979514976 101.03539327220483 151.6624666390998 132.10950240941185 ]
lambda_covariance[14,:] = [94.41218286271872 78.9540303171795 91.09599055638469 103.03646518209119 89.3286003327234 107.24449245970433 92.19765014511928 116.61926047606192 143.3430654992352 111.22492107113902 41.10790805593369 90.88134548128345 127.79426090582622 202.33422202153932 147.3414353557927 102.70255329790933 88.71129652615635 105.38759642668168 79.36849957827465 57.94410750888447 107.7329766891065 135.03498206185427 110.92546205091234 96.22633814509082 98.32290051595548 93.13178074590148 112.33821625297912 123.05591112672042 111.80288696962616 104.43161650433635 ]
lambda_covariance[15,:] = [101.54561514825632 164.55062977418578 113.1514186154446 173.2730472599489 120.10069675893348 169.159644032315 131.39091955202895 163.80844617445038 122.74996401214017 111.74959460831212 61.136046541101024 97.29082329382234 156.21861495807775 147.3414353557927 267.5910314676565 156.02217618065848 120.98645049520043 115.81883477712316 80.5278815700784 100.1982813430703 122.78777706117441 143.18887356334815 140.22765944384537 98.22713408348761 113.36488089683063 119.38685931643155 118.15631866799498 168.37031827147538 129.77141187726525 139.56132456049755 ]
lambda_covariance[16,:] = [75.58072261482403 139.0990608973585 114.00781512870165 135.05919682078394 69.20440336117244 136.10043774466817 155.25184872714766 147.79914811484616 96.47363448295651 111.68415443574753 88.83125425704084 96.95704629238494 129.03158590956014 102.70255329790933 156.02217618065848 244.54156405262924 121.4072672269005 93.27252235897106 97.35737382731983 107.88767756227664 101.64933282827589 145.70570163690397 129.76925560029318 95.14276671824551 94.16414055759252 82.09213816184129 88.95865104970083 117.48656631903424 102.15171527419699 96.54540883944415 ]
lambda_covariance[17,:] = [88.63610384039548 83.0360390704805 109.83987474584998 78.91279561199606 54.55766332776034 112.67651491675228 101.04981875685186 120.55387824692136 54.13930709343893 69.3568004342304 35.46875966095065 89.32792035275538 94.84122633028504 88.71129652615635 120.98645049520043 121.4072672269005 169.74739112249176 88.81290629215607 52.37927918806878 44.36600753422735 82.4898941460954 119.8487775142126 121.78563172544561 75.73069765437182 93.9626131342698 70.1670695381429 106.05481891231253 135.18338298688656 98.54973357658835 83.78520106456992 ]
lambda_covariance[18,:] = [49.17764159493406 60.56909136421002 97.21693698115662 134.27403593476927 42.88288860679177 81.34756593989223 107.21528862617915 91.32202290237825 105.98898820013238 74.84297375124505 93.6317042497511 69.75736694717743 112.69942627897169 105.38759642668168 115.81883477712316 93.27252235897106 88.81290629215607 222.67679008451339 47.81949283030747 79.35397794698883 73.39741920002591 141.70084771144016 108.32822326555257 119.52565577759258 47.56151220134656 112.86556623756398 99.18660257453666 126.39524291101235 78.63085704611365 130.27159288235117 ]
lambda_covariance[19,:] = [37.31226701496719 109.77173322221971 99.07274938398089 109.3293404093198 70.1129057002168 83.16130249300898 102.03579153555103 88.63344171509814 75.0752642886231 73.60853894936437 69.95183665082978 47.288804820190386 113.25924452490982 79.36849957827465 80.5278815700784 97.35737382731983 52.37927918806878 47.81949283030747 122.83302648396234 53.549191980219554 68.08911348855743 63.843484052661985 117.51560122807862 52.52800575974621 66.01863846445356 38.35565661652843 60.73257359279769 64.58877489565346 47.22394282110494 95.83304669087741 ]
lambda_covariance[20,:] = [49.54565921368147 114.39180865661662 61.91574790884258 105.72570077333346 42.01092858068692 87.76451912933427 103.55905941192042 90.13164610517698 83.4365814549703 74.90560279368077 109.04651074800726 73.62626632177022 109.08890327663268 57.94410750888447 100.1982813430703 107.88767756227664 44.36600753422735 79.35397794698883 53.549191980219554 148.28221585251177 77.71921291238297 79.85845848127805 97.03417546439846 95.71634637300191 61.919517449918146 80.3791262297805 42.78857328077803 70.65187066263188 109.3509006212648 90.68638908926752 ]
lambda_covariance[21,:] = [69.10077802061265 105.10413853892037 48.81055594741975 77.08838723329784 59.219600928930824 86.93210109231543 108.67149805022271 101.18720431784811 69.475940356554 76.82970352508781 51.08594178194177 69.0179383069737 109.13719894377545 107.7329766891065 122.78777706117441 101.64933282827589 82.4898941460954 73.39741920002591 68.08911348855743 77.71921291238297 140.5501594093011 85.21895575223107 93.5807423440966 79.44612964851731 114.7122109593612 92.38305386276072 76.07435275857495 115.90943112307608 99.32586831807845 94.092878895825 ]
lambda_covariance[22,:] = [97.0500535647756 106.09606555606173 122.70995140185765 129.9224103074779 42.773303122391596 116.30359078968738 154.3494865577637 163.04528177603024 115.45552301472146 93.21202934881335 110.18111699842564 109.89371269359617 132.49753995973956 135.03498206185427 143.18887356334815 145.70570163690397 119.8487775142126 141.70084771144016 63.843484052661985 79.85845848127805 85.21895575223107 249.0471988631022 120.93694157050453 128.7258881956602 84.48163878420513 106.57265249262821 91.2908243419166 117.59194920874815 139.09668889813443 101.95786170683083 ]
lambda_covariance[23,:] = [82.45750348696772 158.3781620476783 132.45322516874643 174.1724018089613 79.43728980210871 152.94977668133797 171.93735904351246 163.14962123231697 129.86689797715573 97.77826918787798 138.5216597929755 87.36049937979853 127.83072619794666 110.92546205091234 140.22765944384537 129.76925560029318 121.78563172544561 108.32822326555257 117.51560122807862 97.03417546439846 93.5807423440966 120.93694157050453 265.93204520620844 109.60843646375064 100.21288812171693 125.68362956971448 132.86724735828585 150.08201798668418 130.61963475147442 129.70744377071938 ]
lambda_covariance[24,:] = [86.5209387841706 92.33221155079374 101.28991370163848 84.22224814012912 24.671199562469866 84.05254047271808 123.44380609818286 122.03492070329558 86.9572513779304 96.56697770478569 112.19098719130103 83.71193456408744 120.14618848504415 96.22633814509082 98.22713408348761 95.14276671824551 75.73069765437182 119.52565577759258 52.52800575974621 95.71634637300191 79.44612964851731 128.7258881956602 109.60843646375064 169.9668415944231 64.74183952156866 103.90287614274993 74.5303411870715 113.86834442800112 105.82478427044636 92.74218523997548 ]
lambda_covariance[25,:] = [87.361176411108 105.41721793414318 65.45687802210718 57.19019966874822 74.15155681068951 113.90544421320413 140.57366745592802 101.19986987123252 69.0232589342559 64.89053452508826 64.2046627983725 35.72524353388308 83.46265897913469 98.32290051595548 113.36488089683063 94.16414055759252 93.9626131342698 47.56151220134656 66.01863846445356 61.919517449918146 114.7122109593612 84.48163878420513 100.21288812171693 64.74183952156866 182.58476369114072 76.92462038269538 123.56912670204416 105.99310347093034 113.82448383257186 68.65470911665571 ]
lambda_covariance[26,:] = [73.27376303738711 63.31532179652784 48.37629287396258 105.99995730320177 22.198915997472916 103.44502430536937 108.97721149681982 119.24451379862802 114.9377268730485 89.3239441518822 89.67280697294089 67.98938858156919 97.50764555651743 93.13178074590148 119.38685931643155 82.09213816184129 70.1670695381429 112.86556623756398 38.35565661652843 80.3791262297805 92.38305386276072 106.57265249262821 125.68362956971448 103.90287614274993 76.92462038269538 181.3711003309361 99.43616214915832 118.12905617014775 92.30396108237744 86.55555290999305 ]
lambda_covariance[27,:] = [72.68911190643732 76.59433741247877 84.18744456152609 89.71688985457838 59.51626269015065 124.10097509282218 122.32718858409935 119.40031359424987 109.80174246774094 82.2741820294904 66.85389433748807 42.66038402822041 95.35720979514976 112.33821625297912 118.15631866799498 88.95865104970083 106.05481891231253 99.18660257453666 60.73257359279769 42.78857328077803 76.07435275857495 91.2908243419166 132.86724735828585 74.5303411870715 123.56912670204416 99.43616214915832 181.79951480895073 114.51909136716336 98.47150429569807 76.64320790333844 ]
lambda_covariance[28,:] = [86.23715819167262 90.85223309537729 102.59507772963212 107.6003708015894 77.01111734117143 141.66643661902717 122.62272188709112 101.84230987849037 88.99626244262136 82.09140240017425 47.053779930066014 57.36348683116531 101.03539327220483 123.05591112672042 168.37031827147538 117.48656631903424 135.18338298688656 126.39524291101235 64.58877489565346 70.65187066263188 115.90943112307608 117.59194920874815 150.08201798668418 113.86834442800112 105.99310347093034 118.12905617014775 114.51909136716336 222.6221420096601 88.7147359276616 139.6354775418283 ]
lambda_covariance[29,:] = [109.91317968027055 133.67987394225645 87.44936835719288 75.6194872726754 60.10382898897866 90.40825977618357 131.81933536773212 143.75332514297773 107.23059296721462 85.50457130900543 114.76092500933024 95.5158096207522 151.6624666390998 111.80288696962616 129.77141187726525 102.15171527419699 98.54973357658835 78.63085704611365 47.22394282110494 109.3509006212648 99.32586831807845 139.09668889813443 130.61963475147442 105.82478427044636 113.82448383257186 92.30396108237744 98.47150429569807 88.7147359276616 232.93452871103383 84.27441198232314 ]
lambda_covariance[30,:] = [47.06007262000254 107.55622563843521 109.09527613185146 154.81603581461314 87.78972689961464 108.94151833990824 124.05711559790461 98.2834085001364 118.14790803553956 69.98534638939023 79.74429519505637 68.59048459183069 132.10950240941185 104.43161650433635 139.56132456049755 96.54540883944415 83.78520106456992 130.27159288235117 95.83304669087741 90.68638908926752 94.092878895825 101.95786170683083 129.70744377071938 92.74218523997548 68.65470911665571 86.55555290999305 76.64320790333844 139.6354775418283 84.27441198232314 194.7830600905878 ]

const rho_mean = [0.9472892124877721 0.9915658908360583 0.9765528953112144 0.9460146562824032 0.9249960304903848 0.9354324125501662 0.9695132760046203 0.9139711330580986 0.9929738943042129 0.9644243337786643 0.9845852479466539 0.9962986084718584 0.9325704935414124 0.9506730882496088 0.9747173790210193 0.9069298888437348 0.9987146209474168 0.9770695459053615 0.967590781303612 0.956876465969918 ]

const rho_stdev = 0.05
const rho_max = 1.0

const nu = zeros(Float64,20,30)
nu[1,:] = [1.6429381092448916 0.0 1.7306094385198372 1.8450454397960225 1.8522939124210736 0.0 0.0 0.0 2.0332291084784018 1.9879196513663446 2.1293957854191943 2.1364920206265143 2.154271980078755 0.0 0.0 2.310194075025306 0.0 0.0 0.0 2.2096409083226005 0.0 0.0 0.0 0.0 1.942368934095366 1.87413025575146 0.0 1.9649073255513967 0.0 2.259340581099563 ]
nu[2,:] = [0.0 0.0 0.0 0.0 1.9312711893998953 1.9807446547340612 0.0 2.1750155587621975 2.1122063854572235 0.0 2.208373062398016 2.2154692976053365 2.233249257057577 2.02110849736243 1.8739781815097214 0.0 2.1424308727868526 0.0 2.3377309171264935 0.0 2.3280401118562466 0.0 2.5803455770337553 0.0 0.0 1.9531075327302818 2.1543977350273122 2.0438846025302184 1.9181074400942295 2.3383178580783843 ]
nu[3,:] = [1.5080629675749133 0.0 0.0 1.7101702981260445 1.7174187707510953 0.0 0.0 1.9611631401133975 1.8983539668084235 0.0 1.994520643749216 2.0016168789565363 0.0 0.0 1.6601257628609214 2.175318933355328 0.0 1.6845712829304351 2.1238784984776937 2.0747657666526225 0.0 2.3883119266376043 2.366493158384955 0.0 0.0 1.7392551140814818 0.0 0.0 1.7042550214454295 2.1244654394295845 ]
nu[4,:] = [0.0 0.0 1.657610278535361 1.7720462798115466 1.7792947524365976 1.8287682177707636 0.0 2.0230391217988997 0.0 0.0 2.056396625434718 0.0 2.081272820094279 0.0 1.7220017445464237 0.0 1.9904544358235547 0.0 2.185754480163196 0.0 2.1760636748929487 2.4501879083231066 2.4283691400704575 1.9278297458773332 1.86936977411089 1.801131095766984 2.002421298064015 0.0 0.0 2.186341421115087 ]
nu[5,:] = [1.0725957085027367 1.3461592148582655 0.0 1.2747030390538678 1.2819515116789186 0.0 1.0795550156867026 0.0 0.0 1.4175772506241897 1.5590533846770394 0.0 1.5839295793366002 1.3717888196414532 0.0 0.0 1.4931111950658758 1.2491040238582585 1.688411239405517 1.6392985075804456 1.67872043413527 1.9528446675654276 0.0 1.4304865051196543 1.3720265333532111 1.3037878550093052 0.0 0.0 0.0 1.688998180357408 ]
nu[6,:] = [1.5881900405920075 1.8617535469475361 1.6758613698669529 1.7902973711431385 1.7975458437681895 0.0 1.5951493477759735 0.0 0.0 1.9331715827134603 0.0 2.0817439519736305 2.099523911425871 0.0 1.7402528358780156 2.2554460063724218 2.0087055271551466 0.0 2.204005571494788 2.1548928396697162 2.1943147662245406 0.0 2.4466202314020493 0.0 0.0 1.8193821870985758 2.0206723893956067 1.9101592568985126 0.0 0.0 ]
nu[7,:] = [0.0 0.0 1.9206079431317062 2.035043944407892 2.0422924170329426 0.0 1.8398959210407266 2.286036786395245 0.0 2.1779181559782135 2.319394290031063 0.0 2.3442704846906244 0.0 0.0 2.500192579637175 2.2534521004198997 0.0 2.448752144759541 0.0 2.439061339489294 2.7131855729194516 0.0 0.0 0.0 0.0 0.0 0.0 2.029128667727277 0.0 ]
nu[8,:] = [0.0 2.1696249206713514 1.9837327435907681 2.0981687448669537 2.1054172174920045 0.0 0.0 2.349161586854307 2.2863524135493325 2.2410429564372754 2.382519090490125 2.3896153256974455 2.4073952851496863 2.195254525454539 2.0481242096018306 2.5633173800962368 0.0 0.0 0.0 2.4627642133935312 0.0 0.0 2.7544916051258643 0.0 0.0 2.127253560822391 0.0 2.218030630622328 2.092253468186339 2.5124638861704938 ]
nu[9,:] = [1.243949495462372 0.0 0.0 1.446056826013503 1.453305298638554 1.50277876397272 1.250908802646338 1.6970496680008562 0.0 0.0 1.7304071716366747 0.0 1.7552833662962355 1.5431426066010885 1.39601229074838 0.0 1.664464982025511 1.420457810817894 1.8597650263651524 1.810652294540081 1.8500742210949053 0.0 0.0 0.0 1.5433803203128464 0.0 1.6764318442659711 0.0 1.4401415493328882 1.8603519673170432 ]
nu[10,:] = [0.0 0.0 1.9054674283052329 2.019903429581418 0.0 0.0 0.0 0.0 0.0 2.1627776411517403 2.30425377520459 2.3113500104119105 0.0 0.0 1.9698588943162954 2.4850520648107017 2.2383115855934266 0.0 2.4336116299330675 0.0 2.4239208246628206 2.698045058092978 0.0 0.0 0.0 2.0489882455368558 0.0 2.1397653153367924 0.0 2.4341985708849583 ]
nu[11,:] = [1.856747400463013 0.0 1.9444187297379585 0.0 2.066103203639195 0.0 0.0 2.3098475730014973 2.247038399696523 0.0 0.0 2.350301311844636 2.3680812712968766 0.0 2.008810195749021 0.0 2.277262887026152 0.0 2.4725629313657933 0.0 0.0 2.736996359525704 2.7151775912730547 2.2146381970799305 0.0 0.0 0.0 2.178716616769518 0.0 0.0 ]
nu[12,:] = [1.3111293690647696 1.5846928754202985 1.3988006983397152 0.0 0.0 1.5699586375751178 0.0 1.7642295416032538 1.7014203682982798 1.6561109111862227 1.7975870452390723 0.0 0.0 1.6103224802034861 0.0 1.978385334845184 0.0 1.4876376844202914 0.0 0.0 0.0 0.0 2.1695595598748114 0.0 0.0 1.5423215155713381 0.0 1.633098585371275 1.5073214229352858 0.0 ]
nu[13,:] = [1.5369577458953922 1.810521252250921 0.0 1.7390650764465234 0.0 1.7957870144057404 1.5439170530793582 0.0 1.9272487451289024 1.8819392880168453 0.0 0.0 2.048291616729256 1.8361508570341087 0.0 2.2042137116758065 1.9574732324585313 0.0 2.1527732767981727 0.0 0.0 0.0 2.395387936705434 0.0 0.0 0.0 1.9694400946989914 1.8589269622018976 0.0 2.1533602177500635 ]
nu[14,:] = [0.0 0.0 0.0 1.8441351396881362 1.851383612313187 1.9008570776473532 1.648987116320971 0.0 0.0 0.0 0.0 2.135581720518628 2.1533616799708684 1.9412209202757216 0.0 0.0 2.062543295700144 0.0 2.2578433400397855 2.2087306082147142 0.0 0.0 0.0 1.9999186057539227 0.0 0.0 2.0745101579406042 1.9639970254435104 0.0 2.2584302809916763 ]
nu[15,:] = [1.2562675913595183 0.0 0.0 1.4583749219106492 1.4656233945357002 1.5150968598698662 0.0 0.0 0.0 1.6012491334809713 1.742725267533821 1.7498215027411412 0.0 0.0 1.4083303866455263 1.9235235571399327 1.6767830779226574 1.4327759067150403 1.8720831222622987 0.0 0.0 2.1365165504222094 2.11469778216956 1.6141583879764358 0.0 1.4874597378660868 0.0 0.0 1.4524596452300345 1.8726700632141895 ]
nu[16,:] = [0.0 1.9217072254505485 1.7358150483699655 1.8502510496461508 1.8574995222712019 0.0 0.0 0.0 0.0 1.993125261216473 2.1346013952693226 2.141697630476643 0.0 0.0 0.0 0.0 0.0 0.0 2.2639592499978 0.0 2.254268444727553 0.0 2.506573909905062 2.0060345157119377 0.0 1.8793358656015884 0.0 1.970112935401525 1.844335772965536 0.0 ]
nu[17,:] = [1.4198275733211303 0.0 0.0 1.6219349038722612 0.0 0.0 0.0 0.0 0.0 1.7648091154425831 0.0 1.9133814847027533 1.9311614441549936 1.7190206844598466 1.5718903686071384 2.0870835391015445 1.8403430598842694 0.0 0.0 1.986530372398839 0.0 0.0 2.278257764131172 0.0 0.0 1.6510197198276986 1.8523099221247294 1.7417967896276354 0.0 2.0362300451758015 ]
nu[18,:] = [1.3717749016787928 0.0 0.0 1.573882232229924 0.0 0.0 0.0 0.0 0.0 1.7167564438002458 1.8582325778530955 1.8653288130604158 0.0 1.6709680128175093 0.0 0.0 0.0 0.0 1.9875904325815732 0.0 1.977899627311326 2.2520238607414838 2.2302050924888346 0.0 1.6712057265292672 1.6029670481853613 0.0 1.6937441179852981 1.567966955549309 1.988177373533464 ]
nu[19,:] = [1.8120114526440378 0.0 0.0 2.014118783195169 2.02136725582022 2.070840721154386 1.8189707598280038 2.265111625182522 0.0 0.0 2.2984691288183408 0.0 2.323345323477901 2.1112045637827546 0.0 0.0 2.232526939207177 0.0 2.4278269835468183 2.378714251721747 2.418136178276571 0.0 2.6704416434540796 0.0 2.111442277494512 2.0432035991506066 0.0 2.133980668950543 2.008203506514554 0.0 ]
nu[20,:] = [1.652023951470594 1.925587457826123 0.0 1.8541312820217253 1.861379754646776 0.0 1.65898325865456 0.0 2.042314950704104 1.9970054935920472 0.0 2.145577862852217 0.0 1.9512170626093106 1.8040867467566022 2.3192799172510083 2.072539438033733 1.828532266826116 0.0 2.218726750548303 0.0 0.0 2.510454142280636 0.0 1.9514547763210686 1.8832160979771626 2.0845063002741933 1.9739931677770994 0.0 2.2684264233252653 ]

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
module XGBoost

using XGBoost_jll

using Printf
using Random: randperm, seed!
using SparseArrays: SparseMatrixCSC, nnz
using Statistics: mean, std

export DMatrix, Booster
export xgboost, predict, save, nfold_cv, slice, get_info, set_info, dump_model, importance
export rabit_init, rabit_finalize, rabit_is_distributed, rabit_get_rank, rabit_get_world_size, rabit_get_version_number

global const build_version = "0.82"

include("xgboost_wrapper_h.jl")
include("rabit_wrapper.jl")
include("xgboost_lib.jl")

end # module XGBoost


text:
using PyPlot, transD_GP, Random, LinearAlgebra, Statistics

## model
Random.seed!(3)
nlayers, dz = 50, 5.
ρbg        = 1000
zfixed     = [-1e5,   (0:dz:(nlayers-1)*dz )...]
rho        = [1e12,   10 .^(log10(ρbg) .+ randn(nlayers))...]
ianom      = 5
rho[ianom] = 10.
δ          = 1e-4  # fd step size in log10

##  geometry and frequencies for VMD in frequency domain
nFreqsPerDecade     = 5
freqLowLimit        = 1e-3
freqHighLimit       = 1e5
freqs = 10 .^(log10(freqLowLimit):1/nFreqsPerDecade:log10(freqHighLimit))
rRx = 100.
zRx = -37.
zTx = -35.
nkᵣeval = 50
modelprimary = true

Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      calcjacobian = true,
                      nmax = length(zfixed),
                      zTx    = zTx,
                      rRx    = rRx,
                      freqs  = freqs,
                      zRx    = zRx,
                      nkᵣeval = nkᵣeval,
                      modelprimary = modelprimary)
transD_GP.AEM_VMD_HMD.getfieldFD!(Fvmd, zfixed, rho)

## plot
f, ax = plt.subplots(1,3, sharey=true, figsize=(8,6))
delz, delf = diff(zfixed[2:end]), diff(log10.(freqs))
flabel = [(log10.(freqs) .- delf[1]/2)..., log10.(freqs)[end]+delf[end]/2]
zlab = [zfixed[2:end]..., zfixed[end]+delz[end]]
ax[1].step(log10.(rho[2:end]), zfixed[2:end])
im2 = ax[2].pcolormesh(flabel, zlab, abs.(Fvmd.HFD_z_J[2:end,:]), shading="auto")
cb2 = colorbar(im2, ax=ax[2])
im3 = ax[3].pcolormesh(flabel, zlab, log10.(abs.(Fvmd.HFD_z_J[2:end,:])))
cb2 = colorbar(im3, ax=ax[3])
ax[1].set_xlabel("log10 Ω-m")
ax[1].invert_xaxis()
ax[2].invert_yaxis()
ax[2].set_xlabel("log10 Hz")
ax[3].set_xlabel("log10 Hz")
ax[1].set_ylabel("Depth m")
ax[2].set_title("∂f/∂(log10σ)")
ax[3].set_title("∂f/∂(log10σ)")

## approximate partials with central differences and compare
J = copy(Fvmd.HFD_z_J)

function getJ!(Fin, zin, ρ, Δ)
    # sensitivities computed in log10σ
    # same as in AEM_VMD_HMD code
    Jout = similar(Fin.HFD_z_J)
    modelfd = copy(ρ)
    jflag = Fin.calcjacobian
    Fin.calcjacobian = false
    for i = 2:length(zin)
        modelfd[:] = ρ
        modelfd[i] = 10 .^-(-log10(ρ[i]) .+ Δ/2)
        transD_GP.AEM_VMD_HMD.getfieldFD!(Fin, zin, modelfd)
        HFD_z1 = copy(Fin.HFD_z)
        modelfd[i] = 10 .^-(-log10(ρ[i]) .- Δ/2)
        transD_GP.AEM_VMD_HMD.getfieldFD!(Fin, zin, modelfd)
        Jout[i,:] = (HFD_z1 - Fin.HFD_z) / Δ
    end
    Fin.calcjacobian = jflag
    Jout
end
J_fd = getJ!(Fvmd, zfixed, rho, δ)

## comparison plots
f, ax = plt.subplots(1,3, sharey=true, figsize=(10,6))
ax[1].step(log10.(rho[2:end]), zfixed[2:end])
im2 = ax[2].pcolormesh(flabel, zlab, abs.(J[2:end,:]))
cb2 = colorbar(im2, ax=ax[2])
im3 = ax[3].pcolormesh(flabel, zlab, abs.(J_fd[2:end,:]))
cb2 = colorbar(im3, ax=ax[3])
ax[1].set_xlabel("log10 Ω-m")
ax[1].invert_xaxis()
ax[2].invert_yaxis()
ax[2].set_xlabel("log10 Hz")
ax[3].set_xlabel("log10 Hz")
ax[1].set_ylabel("Depth m")
ax[2].set_title("∂f/∂(log10σ) analytic")
ax[3].set_title("∂f/∂(log10σ) FD")

## more detailed
figure()
nlayer = 3
subplot(121)
semilogx(freqs, real(J_fd[nlayer,:]), "--k")
semilogx(freqs, real(J[nlayer,:]), "-")
subplot(122)
semilogx(freqs, imag(J_fd[nlayer,:]), "--k")
semilogx(freqs, imag(J[nlayer,:]), "-")


text:
@enum PARSE_MODE READ_COMMENT READ_SCRIPT READ_STYLE READ_START READ_TAG READ_ATTR READ_ATTR_CONTENT READ_CLOSE_TAG READ_STRING READ_NUM READ_BOOLEAN READ_TAG_CONTENT

whitespaces = r"[\n\t\r]"
alphanum = r"[0-9a-zA-Z]"
num = r"[1-9]"
singleton = r"area|base|br|col|command|embed|hr|img|input|keygen|link|meta|param|source|track|wbr|!doctype"
alphabets = r"[a-zA-Z]"

mutable struct TokenStack
    array::Array{String,1}
    size::Int
end

function push(stack::TokenStack, value::String)
    push!(stack.array, value)
    stack.size += 1
end

function pop(stack::TokenStack)::String
    if stack.size > 0
        stack.size -= 1
        return pop!(stack.array)
    else
        return ""
    end
end

function tokenizer_keep_contents(html::String)::TokenStack
    stack = TokenStack(String[], 0)
    mode = READ_TAG_CONTENT
    store = ""
    store2 = ""
    for c in html
        if mode == READ_TAG_CONTENT
            if c == '<'
                if store != ""
                    push(stack, store)
                end
                store = "<"
                store2 = ""
                mode = READ_TAG
            elseif occursin(whitespaces, string(c)) || (c == ' ' && store == "")
                # Ignore
            elseif c == ' '
                store2 = string(store2, " ")
            else
                if store2 != ""
                    store = string(store, store2)
                    store2 = ""
                end
                store = string(store, c)
            end
        elseif mode == READ_TAG
            if c == '>'
                store = string(store, c)
                push(stack, store)
                store = ""
                mode = READ_TAG_CONTENT
            else
                store = string(store, c)
            end
        end
    end
    if store != ""
        push(stack, store)
    end
    return stack
end

function tokenizer(html::String)::TokenStack
    stack = TokenStack(String[], 0)
    mode = READ_TAG_CONTENT
    store = ""
    store2 = ""
    counter = 0
    has_content = false
    has_error = false
    has_temp_error = false
    has_exclamation = false
    is_singleton = false
    is_closed = false
    is_script = false
    is_style = false
    for c in html
        if mode == READ_START
            if c == '<'
                if has_error
                    push(stack, "BAD_TOKEN")
                end
                store = "TOKEN_"
                has_content = false
                mode = READ_TAG
            else
                has_error = true
            end
        elseif mode == READ_TAG_CONTENT
            if c == '<'
                if has_content
                    push(stack, "CONTENT")
                end
                store = "TOKEN_"
                has_content = false
                mode = READ_TAG
            else
                has_content = !occursin(whitespaces, string(c)) || has_content
            end
        else
            if mode == READ_TAG
                if c == '>'
                    if !is_closed && store2 == "script"
                        mode = READ_SCRIPT
                    elseif !is_closed && store2 == "style"
                        mode = READ_STYLE
                    else
                        mode = READ_TAG_CONTENT
                    end
                    if occursin(singleton, store2) || is_closed || has_exclamation
                        store = string("SINGLE_", store)
                    end
                    if has_content && !has_error && counter == 0
                        push(stack, store)
                    else
                        push(stack, "BAD_TOKEN")
                    end
                    store = ""
                    store2 = ""
                    counter = 0
                    has_error = false
                    has_content = false
                    has_exclamation = false
                    is_closed = false
                elseif is_closed
                    has_error = true
                elseif c == '/'
                    if store == "TOKEN_"
                        store = string(store, "END_")
                        store2 = string(store2, c)
                    elseif !occursin(r"TOKEN_END", store)
                        is_closed = true
                    else
                        has_error = true
                    end
                elseif occursin(alphanum, string(c))
                    store = string(store, uppercase(c))
                    store2 = string(store2, c)
                    has_content = true
                elseif c == '!' && !has_content
                    has_exclamation = true
                elseif c == '-' && !has_content && has_exclamation
                    counter += 1
                    if counter == 2
                        has_exclamation = false
                        counter = 0
                        mode = READ_COMMENT
                    end
                elseif c == ' ' && has_content && store2[1] != '/' && !is_script && !is_style
                    is_script = store2 == "script"
                    is_style = store2 == "style"
                    has_content = false
                    is_singleton = occursin(singleton, store2)
                    store2 = ""
                    mode = READ_ATTR
                elseif c == ' ' && (is_script || is_style)
                    store = string(store, uppercase(c))
                else
                    has_error = true
                end
            elseif mode == READ_COMMENT
                if c == '-'
                    counter += 1
                else
                    if c == '>' && counter >= 2
                        push(stack, "COMMENT")
                        mode = READ_TAG_CONTENT
                    end
                    counter = 0
                end
            elseif mode == READ_SCRIPT
                if c == '<'
                    has_content = store2 != "" || has_content
                    store2 = ""
                elseif c == '>'
                    if store2 == "/script"
                        is_script = false
                        if has_content
                            push(stack, "CONTENT")
                        end
                        push(stack, "TOKEN_END_SCRIPT")
                        mode = READ_TAG_CONTENT
                        has_content = false
                    else
                        has_content = true
                    end
                    store2 = ""
                else
                    store2 = string(store2, c)
                    has_content = true
                end
            elseif mode == READ_STYLE
                if c == '<'
                    has_content = store2 != "" || has_content
                    store2 = ""
                elseif c == '>'
                    if store2 == "/style"
                        is_style = false
                        if has_content
                            push(stack, "CONTENT")
                        end
                        push(stack, "TOKEN_END_STYLE")
                        mode = READ_TAG_CONTENT
                        has_content = false
                    else
                        has_content = true
                    end
                    store2 = ""
                else
                    store2 = string(store2, c)
                    has_content = true
                end
            elseif c == '>'
                if is_singleton || has_exclamation
                    store = string("SINGLE_", store)
                end
                if (store != "TOKEN_" && !has_error && !has_temp_error && mode != READ_ATTR_CONTENT
                    && (mode != READ_STRING || has_content) && counter == 0)
                    push(stack, store)
                else
                    push(stack, "BAD_TOKEN")
                end
                if !is_closed && is_script
                    mode = READ_SCRIPT
                elseif !is_closed && is_style
                    mode = READ_STYLE
                else
                    mode = READ_TAG_CONTENT
                end
                store = ""
                store2 = ""
                counter = 0
                has_error = false
                has_temp_error = false
                has_content = false
                has_exclamation = false
                is_singleton = false
                is_closed = false
                is_script = false
                is_style = false
            elseif c == '/' && (mode != READ_STRING || has_content)
                is_closed = true
                is_singleton = true
            elseif is_closed
                has_error = true
            elseif mode == READ_ATTR
                if occursin(alphanum, string(c))
                    has_content = true
                elseif c == '=' && has_content
                    mode = READ_ATTR_CONTENT
                    has_content = false
                elseif !(c == ' ')
                    has_error = true
                end
            elseif mode == READ_ATTR_CONTENT
                if occursin(alphabets, string(c))
                    has_temp_error = true
                    store2 = string(c)
                    mode = READ_BOOLEAN
                elseif c == '"' || c == '\''
                    mode = READ_STRING
                    store2 = string(c)
                elseif occursin(num, string(c))
                    mode = READ_NUM
                else
                    has_error = true
                end
            elseif mode == READ_BOOLEAN
                if c != ' '
                    store2 = string(store2, c)
                    has_temp_error = store2 != "false" && store2 != "true"
                else
                    store2 = ""
                    has_error = has_error || has_temp_error
                    has_temp_error = false
                    mode = READ_ATTR
                end
            elseif mode == READ_NUM
                if c == '.' && !has_content
                    has_content = true
                elseif c == ' ' && !has_error
                    has_content = false
                    mode = READ_ATTR
                elseif !occursin(num, string(c))
                    has_error = true
                end
            elseif mode == READ_STRING
                if c == store2[begin]
                    if !has_temp_error
                        has_content = true
                    end
                elseif c == ' '
                    has_error = has_error || has_temp_error
                    has_temp_error = false
                    has_content = false
                    mode = READ_ATTR
                elseif c == '\\'
                    has_temp_error = !has_temp_error
                elseif occursin(r"[abfnrtv]", string(c)) && has_temp_error
                    has_temp_error = false
                elseif has_content
                    has_error = true
                end
            end
        end
    end
    if mode == READ_TAG_CONTENT
        if has_content
            push(stack, "CONTENT")
        end
    else
        push(stack, "BAD_TOKEN")
    end
    return stack
end

function check(stack::TokenStack)::Bool
    store_stack = TokenStack(String[], 0)
    has_error = false
    for token in stack.array
        if token == "BAD_TOKEN"
            has_error = true
            break
        elseif occursin(r"TOKEN_END_", token)
            if pop(store_stack) != string("TOKEN_", token[11:end])
                has_error = true
                break
            end
            has_error = false
        elseif token != "CONTENT" && token != "COMMENT" && !occursin(r"SINGLE_", token)
            push(store_stack, token)
            has_error = true
        end
    end
    return !has_error && store_stack.size == 0
end

struct HtmlJSON
    tag::String
    attributes::Dict{String,Any}
    children::Array{Any,1}
end

function init_HtmlJSON(tag::String)::HtmlJSON
    return HtmlJSON(tag, Dict(), Any[])
end

function convert_tojson(html::String)::HtmlJSON
    temp, _ = convert_tojson_rec("WRAPPER", html, 0, READ_START)
    return temp
end

function convert_tojson_rec(tag::String, html::String, position::Int, mode::PARSE_MODE)::Tuple{HtmlJSON, Int}
    result = init_HtmlJSON(tag)
    store = ""
    store2 = ""
    counter = 0
    has_content = false
    has_temp_error = false
    has_exclamation = false
    is_singleton = occursin(singleton, tag)
    is_closed = false
    while position < length(html)
        position += 1
        c = html[position]
        if mode == READ_START
            if c == '<'
                mode = READ_TAG
            else
                throw(ParseError(position, mode))
            end
        elseif mode == READ_TAG_CONTENT
            if c == '<'
                if store != ""
                    push!(result.children, store)
                end
                store = ""
                store2 = ""
                mode = READ_TAG
            elseif occursin(whitespaces, string(c)) || (c == ' ' && store == "")
                # Ignore
            elseif c == ' '
                store2 = string(store2, " ")
            else
                if store2 != ""
                    store = string(store, store2)
                    store2 = ""
                end
                store = string(store, c)
            end
        else
            if mode == READ_TAG
                if c == '>'
                    if !is_closed && store == "script"
                        mode = READ_SCRIPT
                    elseif !is_closed && store == "style"
                        mode = READ_STYLE
                    else
                        if !has_content
                            if occursin(singleton, store) || is_closed
                                push!(result.children, init_HtmlJSON(store))
                            else
                                temp, position = convert_tojson_rec(store, html, position, READ_TAG_CONTENT)
                                push!(result.children, temp)
                            end
                        else
                            if store != result.tag
                                throw(CloseTokenError(position))
                            end
                            return (result, position)
                        end
                        mode = READ_TAG_CONTENT
                    end
                    store = ""
                    has_temp_error = false
                    has_content = false
                    has_exclamation = false
                    is_closed = false
                elseif c == '/' && !has_content && !is_closed
                    if store == ""
                        has_content = true
                    elseif store != ""
                        is_closed = true
                    end
                elseif c == '!' && !has_content
                    has_exclamation = true
                elseif c == '-' && !has_content && has_exclamation
                    counter += 1
                    if counter == 2
                        has_exclamation = false
                        counter = 0
                        mode = READ_COMMENT
                    end
                elseif occursin(alphanum, string(c)) || (c == '!' && store == "")
                    store = string(store, lowercase(c))
                elseif c == ' ' && store != "" && !has_content
                    temp, position = convert_tojson_rec(store, html, position, READ_ATTR)
                    push!(result.children, temp)
                    has_content = false
                    store = ""
                    mode = READ_TAG_CONTENT
                else
                    throw(ParseError(position, mode))
                end
            elseif mode == READ_COMMENT
                if c == '-'
                    counter += 1
                    store2 = "--"
                elseif !occursin(whitespaces, string(c)) || store != ""
                    if c == '>' && counter >= 2
                        for _ in 1:counter-2
                            store = string(store, "-")
                        end
                        temp = init_HtmlJSON("comment")
                        push!(temp.children, store)
                        push!(result.children, temp)
                        store = ""
                        mode = READ_TAG_CONTENT
                    else
                        store = string(store, store2, c)
                    end
                    store2 = ""
                    counter = 0
                end
            elseif mode == READ_SCRIPT
                if c == '<'
                    store = string(store, store2)
                    store2 = "<"
                elseif c == '>'
                    if store2 == "</script"
                        if store != ""
                            temp = init_HtmlJSON("script")
                            push!(temp.children, store)
                            push!(result.children, temp)
                        end
                        mode = READ_TAG_CONTENT
                        store = ""
                    else
                        store = string(store, store2, c)
                    end
                    store2 = ""
                elseif !occursin(whitespaces, string(c)) || store2 != ""
                    store2 = string(store2, c)
                end
            elseif mode == READ_STYLE
                if c == '<'
                    store = string(store, store2)
                    store2 = "<"
                elseif c == '>'
                    if store2 == "</style"
                        if store != ""
                            temp = init_HtmlJSON("style")
                            push!(temp.children, store)
                            push!(result.children, temp)
                        end
                        mode = READ_TAG_CONTENT
                        store = ""
                    else
                        store = string(store, store2, c)
                    end
                    store2 = ""
                elseif !occursin(whitespaces, string(c)) || store2 != ""
                    store2 = string(store2, c)
                end
            elseif c == '>'
                if has_temp_error || mode == READ_ATTR_CONTENT
                        throw(ParseError(position, mode))
                end

                if mode == READ_STRING
                    if !has_content
                        throw(ParseError(position, mode))
                    end
                    result.attributes[store] = store2[begin+1:end-1]
                elseif mode == READ_NUM
                    if has_content
                        result.attributes[store] = parse(Float64, store2)
                    else
                        result.attributes[store] = parse(Int, store2)
                    end
                elseif mode == READ_BOOLEAN
                    result.attributes[store] = store2 == "true"
                elseif mode == READ_ATTR && store != ""
                    result.attributes[store] = ""
                end
                if is_singleton
                    return (result, position)
                end
                mode = READ_TAG_CONTENT
                store = ""
                store2 = ""
                has_temp_error = false
                has_content = false
                has_exclamation = false
            elseif c == '/' && (mode != READ_STRING || has_content)
                is_closed = true
                is_singleton = true
            elseif mode == READ_ATTR
                if occursin(alphanum, string(c))
                    store = string(store, c)
                elseif c == '=' && store != ""
                    mode = READ_ATTR_CONTENT
                elseif c == ' '
                    if store != ""
                        result.attributes[store] = ""
                        store = ""
                    end
                else
                    throw(ParseError(position, mode))
                end
            elseif mode == READ_ATTR_CONTENT
                if occursin(alphabets, string(c))
                    store2 = string(c)
                    has_temp_error = true
                    mode = READ_BOOLEAN
                elseif c == '"' || c == '\''
                    store2 = string(c)
                    mode = READ_STRING
                elseif occursin(num, string(c))
                    store2 = string(c)
                    mode = READ_NUM
                else
                    throw(ParseError(position, mode))
                end
            elseif mode == READ_BOOLEAN
                if c != ' '
                    store2 = string(store2, c)
                    has_temp_error = store2 != "false" && store2 != "true"
                else
                    if has_temp_error
                        throw(ParseError(position, mode))
                    end
                    result.attributes[store] = store2 == "true"
                    store = ""
                    store2 = ""
                    mode = READ_ATTR
                end
            elseif mode == READ_NUM
                if c == '.' && !has_content
                    store2 = string(store2, c)
                    has_content = true
                elseif c == ' ' && !has_temp_error
                    if has_temp_error
                        throw(ParseError(position, mode))
                    end
                    if has_content
                        result.attributes[store] = parse(Float64, store2)
                    else
                        result.attributes[store] = parse(Int, store2)
                    end
                    store = ""
                    store2 = ""
                    has_content = false
                    mode = READ_ATTR
                elseif occursin(num, string(c))
                    store2 = string(store2, c)
                else
                    has_temp_error = true
                end
            elseif mode == READ_STRING
                if c == store2[begin]
                    if !has_temp_error
                        has_content = true
                    else
                        has_temp_error = false
                    end
                    store2 = string(store2, c)
                elseif c == ' '
                    if has_temp_error || !has_content
                        throw(ParseError(position, mode))
                    end
                    result.attributes[store] = store2[begin+1:end-1]
                    store = ""
                    store2 = ""
                    has_content = false
                    mode = READ_ATTR
                elseif c == '\\'
                    store2 = string(store2, c)
                    has_temp_error = !has_temp_error
                elseif occursin(r"[abfnrtv]", string(c)) && has_temp_error
                    store2 = string(store2, c)
                    has_temp_error = false
                elseif has_content
                    throw(ParseError(position, mode))
                else
                    store2 = string(store2, c)
                end
            end
        end
    end
    if tag != "WRAPPER" || mode != READ_TAG_CONTENT || store != ""
        throw(ParseError(position, mode))
    end
    return (result, position)
end

max_stars_repo_name: haverzard/Husk

max_issues_repo_name: haverzard/Husk

max_forks_repo_name: haverzard/Husk


text:
using Documenter, MLJOpenML

makedocs(
    modules = [MLJOpenML],
    sitename = "MLJOpenML.jl",
)

deploydocs(
    repo = "github.com/alan-turing-institute/MLJOpenML.jl.git",
)


text:
export HitRecord

abstract type Hittable end

struct HitRecord{T, M}
    p::Point3{T}
    normal::Vec3{T}
    material::M
    t::T
    front_face::Bool
end
function HitRecord(r::Ray, p, outward_normal, material, t)
    front_face = r.dir ⋅ outward_normal < 0
    normal = front_face ? outward_normal : -outward_normal
    return HitRecord(p, normal, material, t, front_face)
end

function hit(hittables::AbstractArray, r::Ray, t_min, t_max)
    rec = nothing
    closest_so_far = t_max

    for object in hittables
        temp_rec = hit(object, r, t_min, closest_so_far)
        if !isnothing(temp_rec)
            closest_so_far = temp_rec.t
            rec = temp_rec
        end
    end

    return rec
end


text:
include("../benchmark.jl")
using CUTEst, OnePhase
#nlp_raw2 = CUTEstModel("RK23")
#nlp_raw2 = CUTEstModel("FLOSP2TH")

my_pars = OnePhase.Class_parameters()
nlp_raw2 = CUTEstModel("CAMSHAPE")

#nlp_raw2 = CUTEstModel("NET4")

if false
  my_pars.term.max_it = 100
  new_pars = OnePhase.autotune(nlp_raw2, my_pars)
end
#my_pars.term.tol_opt = 1e-6

iter, status, hist, t, err = OnePhase.one_phase_solve(nlp_raw2,my_pars);
#get_fval(iter)
#eval_f(nlp,iter.point.x)
#obj(nlp_raw2,iter.point.x) - get_fval(iter)
if false
using Ipopt
sol = IpoptSolve(nlp_raw2);
end

finalize(nlp_raw2)


text:
# Base.convert(::Type{<:SisoTf}, b::Real) = Base.convert(SisoRational, b)
# Base.convert{T<:Real}(::Type{<:SisoZpk}, b::T) = SisoZpk(T[], T[], b)
# Base.convert{T<:Real}(::Type{<:SisoRational}, b::T) = SisoRational([b], [one(T)])
# Base.convert{T1}(::Type{SisoRational{Vector{T1}}}, t::SisoRational) =  SisoRational(Poly(T1.(t.num.a)),Poly(T1.(t.den.a)))
# Base.convert(::Type{<:StateSpace}, t::Real) = ss(t)
#

#
# function Base.convert{T<:AbstractMatrix{<:Number}}(::Type{StateSpace{T}}, s::StateSpace)
#     AT = promote_type(T, arraytype(s))
#     StateSpace{AT}(AT(s.A),AT(s.B),AT(s.C),AT(s.D), s.Ts, s.statenames, s.inputnames, s.outputnames)
# end

# TODO Fix these to use proper constructors
# NOTE: no real need to convert numbers to transfer functions, have addition methods..
# How to convert a number to either Continuous or Discrete transfer function
# In this case it would be motivated with a "Static" Type
# Base.convert(::Type{<:TransferFunction}, b::Number) = tf([b])
# Base.convert(::Type{<:TransferFunction{<:SisoRational}}, b::Number) = tf(b)
# Base.convert(::Type{<:TransferFunction{<:SisoZpk}}, b::Number) = zpk(b)

Base.convert(::Type{<:TransferFunction{<:SisoZpk{T1, TR1}}}, b::AbstractMatrix{T2}) where {T1, TR1, T2<:Number} = zpk(T1.(b))
Base.convert(::Type{<:TransferFunction{<:SisoRational{T1}}}, b::AbstractMatrix{T2}) where {T1, T2<:Number} = tf(T1.(b))
function convert(::Type{StateSpace{T,MT}}, D::AbstractMatrix{<:Number}) where {T, MT}
    (ny, nu) = size(D)
    A = MT(fill(zero(T), (0,0)))
    B = MT(fill(zero(T), (0,nu)))
    C = MT(fill(zero(T), (ny,0)))
    D = convert(MT, D)
    Ts = 0.0
    return StateSpace{T,MT}(A,B,C,D,Ts)
end

Base.convert(::Type{<:TransferFunction{<:SisoRational{T}}}, b::Number) where {T} = tf(T(b))
Base.convert(::Type{<:TransferFunction{<:SisoZpk{T, TR}}}, b::Number) where {T, TR} = zpk(T(b))
Base.convert(::Type{StateSpace{T, MT}}, b::Number) where {T, MT} = convert(StateSpace{T, MT}, fill(b, (1,1)))

#Base.convert(::Type{<:TransferFunction{<:SisoZpk}}, s::TransferFunction) = zpk(s)
#Base.convert(::Type{<:TransferFunction{<:SisoRational}}, s::TransferFunction) = tf(s)

#
# function Base.convert{T<:Real,S<:TransferFunction}(::Type{S}, b::VecOrMat{T})
#     r = Matrix{S}(size(b,2),1)
#     for j=1:size(b,2)
#         r[j] = vcat(map(k->convert(S,k),b[:,j])...)
#     end
#     hcat(r...)
# end
#

function convert(::Type{TransferFunction{S}}, G::TransferFunction) where S
    Gnew_matrix = Matrix{S}(undef, size(G))
    for i in eachindex(G.matrix)
        Gnew_matrix[i] = convert(S, G.matrix[i])
    end
    return TransferFunction{S}(Gnew_matrix, G.Ts)
end

function convert(::Type{StateSpace{T,MT}}, sys::StateSpace) where {T, MT}
    return StateSpace{T,MT}(convert(MT, sys.A), convert(MT, sys.B), convert(MT, sys.C), convert(MT, sys.D), sys.Ts)
end

function Base.convert(::Type{StateSpace}, G::TransferFunction{<:SisoTf{T0}}) where {T0<:Number}
    T = Base.promote_op(/,T0,T0)
    convert(StateSpace{T,Matrix{T}}, G)
end

function Base.convert(::Type{StateSpace{T,MT}}, G::TransferFunction) where {T<:Number, MT<:AbstractArray{T}}
    if !isproper(G)
        error("System is improper, a state-space representation is impossible")
    end

    # TODO : These are added due to scoped for blocks, but is a hack. This
    # could be much cleaner.
    #T = Base.promote_op(/, T0, T0)

    Ac = Bc = Cc = Dc = A = B = C = D = Array{T}(undef, 0, 0)
    for i=1:ninputs(G)
        for j=1:noutputs(G)
            a, b, c, d = siso_tf_to_ss(T, G.matrix[j, i])
            if j > 1
                # vcat
                Ac = blockdiag(Ac, a)
                Bc = vcat(Bc, b)
                Cc = blockdiag(Cc, c)
                Dc = vcat(Dc, d)
            else
                Ac, Bc, Cc, Dc = a, b, c, d
            end
        end
        if i > 1
            # hcat
            A = blockdiag(A, Ac)
            B = blockdiag(B, Bc)
            C = hcat(C, Cc)
            D = hcat(D, Dc)
        else
            A, B, C, D = Ac, Bc, Cc, Dc
        end
    end
    # A, B, C = balance_statespace(A, B, C)[1:3] NOTE: Use balance?
    return StateSpace{T,MT}(A, B, C, D, G.Ts)
end

siso_tf_to_ss(T::Type, f::SisoTf) = siso_tf_to_ss(T, convert(SisoRational, f))

# Conversion to statespace on controllable canonical form
function siso_tf_to_ss(T::Type, f::SisoRational)

    num0, den0 = numvec(f), denvec(f)
    # Normalize the numerator and denominator,
    # To allow realization of transfer functions that are proper, but now strictly proper
    num = num0 / den0[1]
    den = den0 / den0[1]

    N = length(den) - 1 # The order of the rational function f

    # Get numerator coefficient of the same order as the denominator
    bN = length(num) == N+1 ? num[1] : 0

    if N==0 #|| num == zero(Poly{T})
        A = fill(zero(T), 0, 0)
        B = fill(zero(T), 0, 1)
        C = fill(zero(T), 1, 0)
    else
        A = diagm(1 => fill(one(T), N-1))
        A[end, :] .= -reverse(den)[1:end-1]

        B = fill(zero(T), N, 1)
        B[end] = one(T)

        C = fill(zero(T), 1, N)
        C[1:min(N, length(num))] = reverse(num)[1:min(N, length(num))]
        C[:] -= bN * reverse(den)[1:end-1] # Can index into polynomials at greater inddices than their length
    end
    D = fill(bN, 1, 1)

    return A, B, C, D
end

"""
`A, B, C, T = balance_statespace{S}(A::Matrix{S}, B::Matrix{S}, C::Matrix{S}, perm::Bool=false)`

`sys, T = balance_statespace(sys::StateSpace, perm::Bool=false)`

Computes a balancing transformation `T` that attempts to scale the system so
that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal.
If `perm=true`, the states in `A` are allowed to be reordered.

This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal`
"""
function balance_statespace(A::AbstractMatrix{P}, B::AbstractMatrix{P}, C::AbstractMatrix{P}, perm::Bool=false) where P <: BlasFloat
    nx = size(A, 1)
    nu = size(B, 2)
    ny = size(C, 1)

    # Compute the transformation matrix
    mag_A = abs.(A)
    mag_B = max.(abs.(B), zero(P))
    mag_C = max.(abs.(C), zero(P))
    T = balance_transform(mag_A, mag_B, mag_C, perm)

    # Perform the transformation
    A = T*A/T
    B = T*B
    C = C/T

    return A, B, C, T
end

#TODO Throw a warning here at least?
# Fallback mehod for systems with exotic matrices (i.e. TrackedArrays)
balance_statespace(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, perm::Bool=false) = A,B,C,I

# First try to promote and see of we get BlasFloat, otherwise, fall beack on function above
function balance_statespace(A::Matrix{<:Number}, B::Matrix{<:Number}, C::Matrix{<:Number}, D::Matrix{<:Number}, perm::Bool=false)
    T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
    A2, B2, C2, D2 = promote(A,B,C,D, fill(zero(T)/one(T),0,0)) # If Int, we get Float64
    balance_statespace(A2, B2, C2, perm)
end

function balance_statespace(sys::StateSpace, perm::Bool=false)
    A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm)
    return ss(A,B,C,sys.D), T
end

"""
`T = balance_transform{R}(A::Matrix{R}, B::Matrix{R}, C::Matrix{R}, perm::Bool=false)`

`T = balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(A,B,C,perm)`

Computes a balancing transformation `T` that attempts to scale the system so
that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal.
If `perm=true`, the states in `A` are allowed to be reordered.

No balancing will be done if `A, B C` are not BLAS compatible

This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal`
See also `balance_statespace`, `balance`
"""
function balance_transform(A::AbstractArray, B::AbstractArray, C::AbstractArray, perm::Bool=false)
    nx = size(A, 1)
    # Compute a scaling of the system matrix M
    R = promote_type(eltype(A), eltype(B), eltype(C), Float32) # Make sure we get at least BlasFloat
    T = R[A B; C zeros(R, size(C*B))]

    size(T,1) < size(T,2) && (T = [T; zeros(R, size(T,2)-size(T,1),size(T,2))])
    size(T,1) > size(T,2) && (T = [T zeros(R, size(T,1),size(T,1)-size(T,2))])
    S = diag(balance(T, false)[1])
    Sx = S[1:nx]
    Sio = S[nx+1]
    # Compute permutation of x (if requested)
    pvec = perm ? balance(A, true)[2] * [1:nx;] : [1:nx;]
    # Compute the transformation matrix
    T = zeros(R, nx, nx)
    T[pvec, :] = Sio * diagm(0 => R(1)./Sx)
    return T
end

balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(sys.A,sys.B,sys.C,perm)

convert(::Type{TransferFunction}, sys::StateSpace) = convert(TransferFunction{SisoRational}, sys)

function convert(::Type{TransferFunction{SisoRational{T}}}, sys::StateSpace) where {T<:Number}
    matrix = Matrix{SisoRational{T}}(undef, size(sys))

    A, B, C, D = ssdata(sys)

    # The following follows from the matrix inversion lemma:
    # det(X + uᵀv) = det(X)(1 + vᵀX⁻¹u), or
    # det((sI-A)+BC) = det(sI-A)(1 + C(si-A)⁻¹B)
    # C(si-A)⁻¹B) + D = 1/det(sI-A) * (det((sI-A)+BC) - I + D*det(sI-A))
    charpolyA = charpoly(A)
    for i=1:ninputs(sys), j=1:noutputs(sys)
        num = charpoly(A-B[:,i:i]*C[j:j,:]) - charpolyA + D[j, i]*charpolyA
        matrix[j, i] = SisoRational{T}(num, charpolyA)
    end
    TransferFunction{SisoRational{T}}(matrix, get_Ts(sys))
end
function convert(::Type{TransferFunction{SisoRational}}, sys::StateSpace{T0}) where {T0<:Number}
    T = typeof(one(T0)/one(T0))
    convert(TransferFunction{SisoRational{T}}, sys)
end

function convert(::Type{TransferFunction{SisoZpk{T,TR}}}, sys::StateSpace) where {T<:Number, TR <: Number}
    matrix = Matrix{SisoZpk{T,TR}}(undef, size(sys))

    A, B, C, D = ssdata(sys)

    for i=1:noutputs(sys), j=1:ninputs(sys)
        z, p, k = siso_ss_to_zpk(sys, i, j)
        matrix[i, j] = SisoZpk{T,TR}(z, p, k)
    end
    TransferFunction{SisoZpk{T,TR}}(matrix, get_Ts(sys))
end
function convert(::Type{TransferFunction{SisoZpk}}, sys::StateSpace{T0}) where {T0<:Number}
    T = typeof(one(T0)/one(T0))
    convert(TransferFunction{SisoZpk{T,complex(T)}}, sys)
end

"""
Convert get zpk representation of sys from input j to output i
"""
function siso_ss_to_zpk(sys, i, j)
    A, B, C = struct_ctrb_obsv(sys.A, sys.B[:, j:j], sys.C[i:i, :])
    D = sys.D[i:i, j:j]
    z = tzero(A, B, C, D)
    nx = size(A, 1)
    nz = length(z)
    k = nz == nx ? D[1] : (C*(A^(nx - nz - 1))*B)[1]
    return z, eigvals(A), k
end

# TODO: Could perhaps be made more accurate. See: An accurate and efficient
# algorithm for the computation of the # characteristic polynomial of a general square matrix.
function charpoly(A::AbstractMatrix{<:Number})
    Λ = eigvals(A)

    return prod(roots2poly_factors(Λ)) # Compute the polynomial factors directly?
end
function charpoly(A::AbstractMatrix{<:Real})
    Λ = eigvals(A)
    return prod(roots2real_poly_factors(Λ))
end

# function charpoly(A)
#     λ = eigvals(A);
#     T = promote_type(primitivetype(A), Float64)
#     I = one(T)
#     p = reduce(*,Poly([I]), Poly[Poly([I, -λᵢ]) for λᵢ in λ]);
#     if maximum(imag.(p[:])./(I+abs.(real.(p[:])))) < sqrt(eps(T))
#         for i = 1:length(p)
#             p[i] = real(p[i])
#         end
#     else
#         error("Characteristic polynomial should be real")
#     end
#     p
# end


text:
using Test
using POMDPs
using POMDPTesting
using POMDPModelTools

import POMDPs:
    transition,
    observation,
    initialstate,
    updater,
    states,
    actions,
    observations

struct TestPOMDP <: POMDP{Bool, Bool, Bool} end
updater(problem::TestPOMDP) = DiscreteUpdater(problem)
initialstate(::TestPOMDP) = BoolDistribution(0.0)
transition(p::TestPOMDP, s, a) = BoolDistribution(0.5)
observation(p::TestPOMDP, a, sp) = BoolDistribution(0.5)
states(p::TestPOMDP) = (true, false)
actions(p::TestPOMDP) = (true, false)
observations(p::TestPOMDP) = (true, false)

@testset "model" begin
    m = TestPOMDP()
    @test has_consistent_initial_distribution(m)
    @test has_consistent_transition_distributions(m)
    @test has_consistent_observation_distributions(m)
    @test has_consistent_distributions(m)
end

@testset "old model" begin
    probability_check(TestPOMDP())
end

@testset "support mismatch" begin
    struct SupportMismatchPOMDP <: POMDP{Int, Int, Int} end
    POMDPs.states(::SupportMismatchPOMDP) = 1:2
    POMDPs.actions(::SupportMismatchPOMDP) = 1:2
    POMDPs.observations(::SupportMismatchPOMDP) = 1:2
    POMDPs.initialstate(::SupportMismatchPOMDP) = Deterministic(3)
    POMDPs.transition(::SupportMismatchPOMDP, s, a) = SparseCat([1, 2, 3], [1.0, 0.0, 0.1])
    POMDPs.observation(::SupportMismatchPOMDP, s, a, sp) = SparseCat([1, 2, 3], [1.0, 0.0, 0.1])
    @test !has_consistent_transition_distributions(SupportMismatchPOMDP())
    @test !has_consistent_observation_distributions(SupportMismatchPOMDP())
    @test !has_consistent_distributions(SupportMismatchPOMDP())
end