text:
include("../btg.jl")
include("../bayesopt/updateBTG.jl")
using Test

# Test updateBTG.jl in bayesopt folder
# tests if trainingData, testingData, and 3 kinds of buffers are upadated correctly
#

x = [0.5 0.1; 0.15 0.4; 0.4 0.6; 0.2 0.3];
Fx = linear_polynomial_basis(x);
y = [1.0, 3.0, 1.5, 2.6];
x0 = [0.8 0.35];
Fx0 = linear_polynomial_basis(x0);
range_theta = [1.0 5];
range_lambda = [2.0 4]
train_data = extensible_trainingData(x, Fx, y);
test_data = testingData(x0, Fx0);
b1 = btg(train_data, range_theta, range_lambda);
x_new = [.2 .7]
Fx_new = linear_polynomial_basis(x_new)
y_new = [1.6]

lambda_keys = keys(b1.λbuffer_dict)
theta_keys = keys(b1.train_buffer_dict)
arb_lambda_key = first(lambda_keys)
arb_theta_key = first(theta_keys)

randb = rand(5,1)

test_data = testingData(x0, Fx0)

#println(b1.train_buffer_dict[arb_theta_key].choleskyΣθ \ b])

updateBTG!(b1, x_new, Fx_new, y_new) #should update trainingData, and all three buffers

if true #test set
@testset "trainingData update test" begin
    @test b1.trainingData.x[b1.trainingData.n:b1.trainingData.n, :] == x_new
    @test b1.trainingData.Fx[b1.trainingData.n:b1.trainingData.n, :] == Fx_new
    @test [b1.trainingData.y[b1.trainingData.n]] == y_new
end

@testset "λbuffer update test" begin
    g = x -> b1.g(x, arb_lambda_key)
    dg = x -> partialx(b1.g, x, arb_lambda_key)
    #@info b1.λbuffer_dict[arb_lambda_key].dgλz[end]
    #@info dg(y_new)
    #@info b1.λbuffer_dict[arb_lambda_key].gλz[end]
    #@info y_new
    @test [b1.λbuffer_dict[arb_lambda_key].dgλz[end]] == dg(y_new)
    @test [b1.λbuffer_dict[arb_lambda_key].gλz[end]] == g(y_new)
end

@testset "train_buffer update test" begin
    x = getPosition(b1.trainingData)
    Fx = getCovariates(b1.trainingData)
    chol = cholesky( correlation(Gaussian(), arb_theta_key, x) + 1e-12 * Matrix(I, 5, 5)) #assume fixed noise level of 1e-12
    #@info chol\randb
    #@info b1.train_buffer_dict[arb_theta_key].choleskyΣθ\randb
    @test norm(chol\randb - b1.train_buffer_dict[arb_theta_key].choleskyΣθ\randb) <1e-6

    @test norm(b1.train_buffer_dict[arb_theta_key].Σθ_inv_X - chol\Fx) <1e-6
    v = rand(3, 1)
    @test norm(b1.train_buffer_dict[arb_theta_key].choleskyXΣX\v  - (Fx' * (chol\Fx))\v)<1e-6
    @test norm(b1.train_buffer_dict[arb_theta_key].logdetΣθ - logdet(chol))<1e-6
    @test norm(b1.train_buffer_dict[arb_theta_key].logdetXΣX - logdet(Fx' * (chol\Fx)))<1e-6
end

#@testset "θλbuffer update test"

println("test solve b1...")
(pdf, cdf, dpdf, augmented_cdf_deriv, augmented_cdf_hess, quantInfo, tgridpdf, tgridcdf, tgridm, tgridsigma_m) = solve(b1);
pdf(x0, Fx0, 2) #implicitly initializes b1.testingData

#println(b1.testingData.x0)
#println(b1.testingData.Fx0)
#println(b1.test_buffer_dict[arb_theta_key].Bθ)

@info norm(Bθ - Bθ_actual)

@testset "test_buffer update tests" begin
    Bθ = b1.test_buffer_dict[arb_theta_key].Bθ
    Eθ = b1.test_buffer_dict[arb_theta_key].Eθ
    Dθ = b1.test_buffer_dict[arb_theta_key].Dθ
    Hθ = b1.test_buffer_dict[arb_theta_key].Hθ

    Bθ_actual = cross_correlation(Gaussian(), arb_theta_key, getPosition(b1.testingData), getPosition(b1.trainingData))
    Eθ_actual = [1]
    Dθ_actual = Eθ_actual - Bθ_actual * (b1.train_buffer_dict[arb_theta_key].choleskyΣθ \Bθ_actual')
    Hθ_actual = getCovariates(b1.testingData) - Bθ_actual * b1.train_buffer_dict[arb_theta_key].Σθ_inv_X

    @test norm(Bθ - Bθ_actual) <1e-6
    @test norm(Eθ - Eθ_actual) <1e-6
    @test norm(Dθ - Dθ_actual) <1e-6
    @test norm(Hθ - Hθ_actual) <1e-6
end
if false
#@testset "test multiple btg update" begin
    (pdf, cdf, dpdf, augmented_cdf_deriv, augmented_cdf_hess, quantInfo, tgridpdf, tgridcdf, tgridm, tgridsigma_m) = solve(b1);
    x_new = [1.3 .7]
    y_new = [.6]
    @info "cdf(x_new, linear_polynomial_basis(x_new), .4)" cdf(x_new, linear_polynomial_basis(x_new), .4)
    @info "cdf([.5 .9], linear_polynomial_basis([.5 .9]), .3)" cdf([.5 .9], linear_polynomial_basis([.5 .9]), .3)
    updateBTG!(b1, x_new, linear_polynomial_basis(x_new), y_new)
    (pdf, cdf, dpdf, augmented_cdf_deriv, augmented_cdf_hess, quantInfo, tgridpdf, tgridcdf, tgridm, tgridsigma_m) = solve(b1);
    @info "cdf(x_new, linear_polynomial_basis(x_new), .4)" cdf(x_new, linear_polynomial_basis(x_new), .41)
    @info "cdf([.5 .9], linear_polynomial_basis([.5 .9]), .3)" cdf([.5 .9], linear_polynomial_basis([.5 .9]), .3)
#end

end

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
    set_param!(ref::ComponentReference, name::Symbol, value)

Set a component parameter as `set_param!(reference, name, value)`.
This creates a unique name :compname_paramname in the model's external parameter list,
and sets the parameter only in the referenced component to that value.
"""
function set_param!(ref::ComponentReference, name::Symbol, value)
    compdef = find_comp(ref)
    unique_name = Symbol("$(compdef.name)_$name")
    set_param!(parent(ref), compdef, name, unique_name, value)
end

"""
    update_param!(ref::ComponentReference, name::Symbol, value)

Update a component parameter as `update_param!(reference, name, value)`.
This uses the unique name :compname_paramname in the model's external parameter list,
and updates the parameter only in the referenced component to that value.
"""
function update_param!(ref::ComponentReference, name::Symbol, value)
    compdef = find_comp(ref)
    unique_name = Symbol("$(compdef.name)_$name")
    update_param!(parent(ref), unique_name, value)
end

"""
    Base.setindex!(ref::ComponentReference, value, name::Symbol)

Set a component parameter as `reference[name] = value`.
This creates a unique name :compname_paramname in the model's external parameter list,
and sets the parameter only in the referenced component to that value.
"""
function Base.setindex!(ref::ComponentReference, value, name::Symbol)
    compdef = find_comp(ref)
    unique_name = Symbol("$(compdef.name)_$name")
    if has_parameter(parent(ref), unique_name)
        update_param!(ref, name, value)
    else
        set_param!(ref, name, value)
    end
end

"""
    connect_param!(dst::ComponentReference, dst_name::Symbol, src::ComponentReference, src_name::Symbol)

Connect two components as `connect_param!(dst, dst_name, src, src_name)`.
"""
function connect_param!(dst::ComponentReference, dst_name::Symbol, src::ComponentReference, src_name::Symbol)
    _connect_param!(parent(dst), pathof(dst), dst_name, pathof(src), src_name)
end

"""
    connect_param!(dst::ComponentReference, src::ComponentReference, name::Symbol)

Connect two components with the same name as `connect_param!(dst, src, name)`.
"""
function connect_param!(dst::ComponentReference, src::ComponentReference, name::Symbol)
    _connect_param!(parent(dst), pathof(dst), name, pathof(src), name)
end

"""
    Base.getindex(ref::ComponentReference, name::Symbol)

Get a sub-comp, parameter, or variable reference as `ref[name]`.
"""
function Base.getindex(ref::ComponentReference, name::Symbol)
    VariableReference(ref, name)
end

# Methods to convert components, params, and vars to references
# for use with getproperty() chaining.
function _make_reference(obj::AbstractComponentDef, comp_ref::ComponentReference)
    return ComponentReference(parent(obj), nameof(obj))
end

function _make_reference(obj::VariableDef, comp_ref::ComponentReference)
    return VariableReference(comp_ref, obj.name)
end

function _make_reference(obj::ParameterDef, comp_ref::ComponentReference)
    return ParameterReference(comp_ref, obj.name)
end

function Base.getproperty(ref::ComponentReference, name::Symbol)
    comp = find_comp(ref)
    return _make_reference(getfield(comp, name), ref) # might be ref to comp, var, or param
end

function _same_composite(ref1::AbstractComponentReference, ref2::AbstractComponentReference)
    # @info "same_composite($(ref1.comp_path), $(ref2.comp_path))"
    return head(pathof(ref1)) == head(pathof(ref2))
end

"""
    Base.setindex!(comp_ref::ComponentReference, var_ref::VariableReference, var_name::Symbol)

Connect two components as `comp_ref[var_name] = var_ref`.
"""
function Base.setindex!(comp_ref::ComponentReference, var_ref::VariableReference, vname::Symbol)
    _same_composite(comp_ref, var_ref) || error("Cannot connect variables defined in different composite trees")

    _connect_param!(parent(comp_ref), pathof(comp_ref), vname, pathof(var_ref), var_name(var_ref))
end


text:
@testset "test version parser" begin

    @test Pandoc.pandoc_api_version([1]) == v"1.0.0"
    @test Pandoc.pandoc_api_version([1, 2]) == v"1.2"
    @test Pandoc.pandoc_api_version([1, 2, 3]) == v"1.2.3"
    @test Pandoc.pandoc_api_version([1, 2, 3, 4]) == v"1.2.3-4"
    @test Pandoc.pandoc_api_version([1, 2, 3, 4, 5]) == v"1.2.3-4.5"

end

@testset "test conversions" begin

    @test (
           convert(Pandoc.Header, Markdown.Header{1}(Any["Header 1"]))
           ==
           Pandoc.Header(1, Pandoc.Attributes("header-1", [], []), Pandoc.Element[Pandoc.Str("Header"), Pandoc.Space(), Pandoc.Str("1")])
          )

    @test (
           convert(Pandoc.Str, "hello")
           ==
           Pandoc.Str("hello")
          )

    @test (
           convert(Pandoc.Link, Markdown.Link("title", "https://example.com"))
           ==
           Pandoc.Link(
                       Pandoc.Attributes(),
                       Pandoc.Inline[Pandoc.Str("title")],
                       Pandoc.Target("https://example.com", ""),
                      )
          )

    @test (
           convert(Pandoc.Link, Markdown.Link(Any["title"], "https://example.com"))
           ==
           Pandoc.Link(
                       Pandoc.Attributes(),
                       Pandoc.Inline[Pandoc.Str("title")],
                       Pandoc.Target("https://example.com", ""),
                      )
          )

    @test (
           convert(Pandoc.Strong, Markdown.Bold(["bold"]))
           ==
           Pandoc.Strong(Pandoc.Inline[Pandoc.Str("bold")])
          )

    @test (
           convert(Pandoc.Emph, Markdown.Italic(["italic"]))
           ==
           Pandoc.Emph(Pandoc.Inline[Pandoc.Str("italic")])
          )

    @test (
           convert(Pandoc.Para, Markdown.Paragraph(Any["this ", Markdown.Bold(Any["is"]), " ", Markdown.Italic(Any["a"]), " ", Markdown.Link(Any["link"], "https://example.com")]))
           ==
           Pandoc.Para(
                       Pandoc.Inline[
                                     Pandoc.Str("this"),
                                     Pandoc.Space(),
                                     Pandoc.Strong(Pandoc.Inline[Pandoc.Str("is")]),
                                     Pandoc.Space(),
                                     Pandoc.Emph(Pandoc.Inline[Pandoc.Str("a")]),
                                     Pandoc.Space(),
                                     Pandoc.Link(Pandoc.Attributes(), Pandoc.Inline[Pandoc.Str("link")], Pandoc.Target("https://example.com", "")),
                                    ]
                      )
          )

    @test (
           convert(Pandoc.HorizontalRule, Markdown.HorizontalRule())
           ==
           Pandoc.HorizontalRule()
          )

    @test (
           convert(Pandoc.LineBreak, Markdown.LineBreak())
           ==
           Pandoc.LineBreak()
          )

    @test (
           convert(Pandoc.BlockQuote, Markdown.BlockQuote(Any[Markdown.Paragraph(Any["this is a paragraph"])]))
           ==
           Pandoc.BlockQuote(
                             Pandoc.Block[
                                          Pandoc.Para(
                                                      Pandoc.Inline[
                                                                    Pandoc.Str("this"),
                                                                    Pandoc.Space(),
                                                                    Pandoc.Str("is"),
                                                                    Pandoc.Space(),
                                                                    Pandoc.Str("a"),
                                                                    Pandoc.Space(),
                                                                    Pandoc.Str("paragraph")
                                                            ],
                                              ),
                                         ]
                            )
          )

    @test (
           convert(Pandoc.CodeBlock, Markdown.Code("x = π \ny = π / 2"))
           ==
           Pandoc.CodeBlock(
                            Pandoc.Attributes("", [""], []),
                            "x = π \ny = π / 2"
                           )
          )

    @test (
           convert(Pandoc.Code, Markdown.Code("x = π"))
           ==
           Pandoc.Code(
                       Pandoc.Attributes("", [""], []),
                       "x = π"
                      )
          )

    @test (
           convert(Pandoc.BulletList, Markdown.List(Any[Any[Markdown.Paragraph(Any["one"])], Any[Markdown.Paragraph(Any["two"])]], -1, false))
           ==
           Pandoc.BulletList(
                              [
                               [
                                Pandoc.Plain(
                                     Pandoc.Inline[Pandoc.Str("one")],
                                    )
                               ],
                               [
                                Pandoc.Plain(
                                     Pandoc.Inline[Pandoc.Str("two")],
                                    )
                               ]
                              ]
                             )
          )

    @test (
           convert(Pandoc.OrderedList, Markdown.List(Any[Any[Markdown.Paragraph(Any["one"])], Any[Markdown.Paragraph(Any["two"])]], 1, false))
           ==
           Pandoc.OrderedList(
                              Pandoc.ListAttributes(
                                                    1,
                                                    Pandoc.Decimal,
                                                    Pandoc.Period
                                                   ),
                              [
                               [
                                Pandoc.Plain(
                                     Pandoc.Inline[Pandoc.Str("one")],
                                    )
                               ],
                               [
                                Pandoc.Plain(
                                     Pandoc.Inline[Pandoc.Str("two")],
                                    )
                               ]
                              ]
                             )
          )

    @test (
           convert(Pandoc.Math, Markdown.LaTeX("x = π"))
           ==
           Pandoc.Math(Pandoc.InlineMath, "x = π")
          )

    @test (
           convert(Pandoc.Image, Markdown.Image("./image.png", ""))
           ==
           Pandoc.Image(Pandoc.Attributes(), [], Pandoc.Target("./image.png", ""))
          )

    @test (
           convert(Pandoc.Note, Markdown.Footnote("1", nothing))
           ==
           Pandoc.Note(Pandoc.Block[])
          )

    @test (
           convert(Pandoc.BulletList,
                   Markdown.List(
                                 Any[
                                     Any[
                                         Markdown.Paragraph(Any["Item 1"]),
                                         Markdown.List(
                                                       Any[
                                                           Any[
                                                               Markdown.Paragraph(Any["Sub Item 1"]),
                                                               Markdown.List(
                                                                             Any[
                                                                                Any[Markdown.Paragraph(Any["Sub sub item 1"])],
                                                                                Any[Markdown.Paragraph(Any["Lorem ipsum"])]
                                                                                ], -1, false
                                                                            )
                                                              ],
                                                           Any[Markdown.Paragraph(Any["Sub Item 2"])]
                                                          ], -1, false
                                                      )
                                        ],
                                     Any[Markdown.Paragraph(Any["Item 2"])],
                                     Any[Markdown.Paragraph(Any["Item 3"])]
                                    ], -1, false
                                )
                  )
           ==
            Pandoc.BulletList(
                              Array{Pandoc.Block,1}[
                                                    [
                                                     Pandoc.Plain(
                                                                  Pandoc.Inline[
                                                                                Pandoc.Str("Item"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("1")]
                                                                 ),
                                                     Pandoc.BulletList(
                                                                       Array{Pandoc.Block,1}[
                                                                                [
                                                                                Pandoc.Plain(
                                                                                Pandoc.Inline[
                                                                                Pandoc.Str("Sub"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("Item"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("1")
                                                                                ]
                                                                                ),
                                                                                Pandoc.BulletList(
                                                                                Array{Pandoc.Block,1}[
                                                                                [
                                                                                Pandoc.Plain(
                                                                                Pandoc.Inline[
                                                                                Pandoc.Str("Sub"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("sub"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("item"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("1")
                                                                                ]
                                                                                )
                                                                                ],
                                                                                [
                                                                                Pandoc.Plain(
                                                                                Pandoc.Inline[
                                                                                Pandoc.Str("Lorem"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("ipsum"),
                                                                                ]
                                                                                )
                                                                                ]
                                                                                ]
                                                                                )
                                                                                ],
                                                                                [
                                                                                Pandoc.Plain(
                                                                                Pandoc.Inline[
                                                                                Pandoc.Str("Sub"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("Item"),
                                                                                Pandoc.Space(),
                                                                                Pandoc.Str("2")
                                                                                ]
                                                                                )
                                                                                ]
                                                                                ]
                                                                      )
                                                    ],
                                                    [
                                                     Pandoc.Plain(Pandoc.Inline[Pandoc.Str("Item"), Pandoc.Space(), Pandoc.Str("2")])
                                                    ],
                                                    [
                                                     Pandoc.Plain(Pandoc.Inline[Pandoc.Str("Item"), Pandoc.Space(), Pandoc.Str("3")])
                                                    ],
                                                   ]
                             )
          )

end

@testset "test pandoc parser" begin

    @testset "docbook" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "docbook-reader.docbook"))) == Pandoc.Document
    end

    @testset "creole" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "creole-reader.txt"))) == Pandoc.Document
    end

    @testset "markdown" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "markdown-reader-more.txt"))) == Pandoc.Document
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "writer.markdown"))) == Pandoc.Document
    end

    @testset "wiki" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "mediawiki-reader.wiki"))) == Pandoc.Document
    end

    @testset "ipynb" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "simple.ipynb"))) == Pandoc.Document
    end

    @testset "latex" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "latex-reader.latex"))) == Pandoc.Document
    end

    @testset "man" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "man-reader.man"))) == Pandoc.Document
    end

    @testset "tables" begin
        @test typeof(Pandoc.parse_file(joinpath(@__DIR__, "data", "pipe-tables.txt"))) == Pandoc.Document
    end
end


text:
function spdiags(B, d, m, n)
    d = d[:]
    p = length(d)

    len = zeros(p + 1, 1)
    for k in 1:p
        len[k+1] =
            int(len[k] + length(Base.max(1, 1 - d[k]):Base.min(m, n - d[k])))
    end
    a = zeros(int(len[p+1]), 3)
    for k in 1:p
        # Append new d[k]-th diagonal to compact form
        i = Base.max(1, 1 - d[k]):Base.min(m, n - d[k])
        a[(int(len[k])+1):int(len[k+1]), :] = [i i + d[k] B[i+(m>=n)*d[k], k]]
    end

    A = sparse(int(a[:, 1]), int(a[:, 2]), a[:, 3], m, n)

    return A
end

y = readcsv("snp_500.txt")
n = length(y)

e = ones(n, 1)
D = spdiags([e -2 * e e], 0:2, n - 2, n)

lambda = 50
start = time()
D * x
println(time() - start)

x = Variable(n)
p = minimize(lambda * sum(D * x))
println(time() - start)
solve!(p)
p.optval
println(time() - start)


text:
using Test, LinearAlgebra
using QBase

@testset "./src/types/measurements.jl" begin

@testset "is_povm_element()" begin
    @testset "valid cases" begin
        @test is_povm_element([1 0;0 0])
        @test is_povm_element([0.5 -0.5im;0.5im 0.5])
        @test is_povm_element([0.5 0.2;0.3 0.5], atol=0.11)
    end

    @testset "invalid cases" begin
        @test !is_povm_element([0.5 0.2;0.3 0.5])
        @test !is_povm_element([1 0.1;0.1 0])
    end
end

@testset "is_povm()" begin
    @testset "valid cases" begin
        E1 = sqrt(2)/(1 + sqrt(2)) * [0 0;0 1]
        E2 = sqrt(2)/(1 + sqrt(2))*(1/2) * [1 -1; -1 1]
        E3 = [1 0;0 1] - E1 - E2

        @test is_povm([E1,E2,E3])

        @test is_povm([
            [1/3 0; 0 1/3],
            2/3*[-cos(2*π/3) 0.5*sin(2*π/3); 0.5*sin(2*π/3) -cos(2*π/3)],
            2/3*[-cos(2*π/3) -0.5*sin(2*π/3); -0.5*sin(2*π/3) -cos(2*π/3)]
        ])
    end

    @testset "negative test case" begin
        @test !is_povm([[1 0;0 0],[0 0;0 1],[0 0;0 0.1]])
    end

    @testset "atol" begin
        Π = [[0.5 0.6;0.6 0.5],[0.5 -0.5;-0.5 0.5]]
        @test !is_povm(Π)
        @test is_povm(Π,atol=0.2)
    end

    @test is_povm(POVM([[1 0;0 0],[0 0;0 1]]))
end

@testset "POVMel()" begin
    @testset "valid povm elements" begin
        M = POVMel(State([1 0 0;0 0 0;0 0 0]))
        @test M isa POVMel{Int64}
        @test M == [1 0 0;0 0 0;0 0 0]

        M = POVMel([0.5 0.5im;-0.5im 0.5])
        @test M isa POVMel{Complex{Float64}}
        @test M == [0.5 0.5im;-0.5im 0.5]

        M = POVMel([0.5 0.45im;-0.5im 0.5], atol=0.06)
        @test M isa POVMel{Complex{Float64}}
        @test M == [0.5 0.45im;-0.5im 0.5]
    end

    @testset "errors" begin
        @test_throws DomainError POVMel([0.5 0.5im;0.5im 0.5])
    end
end

@testset "kron for POVMel" begin
    M1 = POVMel([1 0;0 0])
    M2 = POVMel([0 0;0 1.])

    M1_M2 =  kron(M1,M2)
    @test M1_M2 isa POVMel{Float64}
    @test M1_M2 == [0 0 0 0;0 1 0 0;0 0 0 0;0 0 0 0]
end

@testset "POVM()" begin
    @testset "valid povms:" begin
        Π = [[1 0;0 0],[0 0;0 1]]
        povm = POVM(Π)
        @test povm isa POVM{Int64}
        @test povm isa Measurement
        @test povm == Π
        @test length(povm) == 2

        Π = [[0.5 (0.1-0.1*im);(0.1+0.1*im) 0.5],[0.5 -(0.1-0.1*im);-(0.1+0.1*im) 0.5]]
        povm = POVM(Π)
        @test povm isa POVM{Complex{Float64}}
        @test povm isa Measurement
        @test povm == Π
        @test length(povm) == 2

        Π = [[2/3 0;0 0],[1/6 1/sqrt(12);1/sqrt(12) 1/2],[1/6 -1/sqrt(12);-1/sqrt(12) 1/2]]
        povm = POVM(Π)
        @test povm isa POVM{Float64}
        @test povm isa Measurement
        @test povm == Π
        @test length(povm) == 3
    end

    @testset "invalid povms" begin
        invalid_povms = [
            [[1 0;0 0],[0 0;0 1],[0.5 0;0 0.5]],
            [[0.5 (0.1+0.1*im);(0.1+0.1*im) 0.5],[0.5 -(0.1+0.1*im);-(0.1+0.1*im) 0.5]],
            [[0.5 0.5;0.5 0.5],[0.5 -0.5*im;-0.5*im 0.5]]
        ]

        for Π in invalid_povms
            @test_throws DomainError POVM(Π)
        end
    end

    @testset "atol" begin
        Π = [[1 0;0 0],[0 0.1;0 1]]
        @test_throws DomainError POVM(Π)
        @test POVM(Π, atol=0.11) isa POVM{Float64}
    end
end

@testset "PVM()" begin
    @testset "valid cases" begin
        Π = [[1,0],[0,1]]
        pvm = PVM(Π)
        @test pvm isa PVM{Int64}
        @test pvm isa Measurement
        @test pvm == Π
        @test length(pvm) == 2

        pvm = PVM(Ket.(Π))
        @test pvm isa PVM{Int64}
        @test pvm isa Measurement
        @test pvm == Π

        Π = [[1,-1,-im,-im],[1,-1,im,im],[1,1,im,-im],[1,1,-im,im]]/2
        pvm = PVM(Π)
        @test pvm isa PVM{Complex{Float64}}
        @test pvm isa Measurement
        @test pvm == Π
        @test length(pvm) == 4
    end

    @testset "invalid cases" begin
        @test_throws MethodError PVM([[1 0;0 0],[0 0;0 1]])
        @test_throws DomainError PVM([[1,0,0],[0,1,0]])
        @test_throws DomainError PVM([[0.5,0.5],[0.5,-0.5]])
    end

    @testset "atol" begin
        @test_throws DomainError PVM([[1,1],[1,1]])
        @test PVM([[1,1],[1,1]], atol=2) isa PVM{Int64}
    end
end

end


text:
function output(transaction_result)
    json = transaction_result.output
    return RelationIterator(TrieWalker(relations(json)))
end

struct RelationIterator  # TODO: Name?
    iterator::TrieWalker
end

function Base.show(io::IO, ::MIME"text/plain", cursor::RelationIterator)
    num_rels = length(cursor.iterator.relations)
    tups = tuples(cursor)
    print(io, "$RelationIterator with $(length(tups)) tuples in $(num_rels) physical relations:")
    _show_trie(io, cursor)
end

Base.keytype(_::RelationIterator) = Any
Base.keys(cursor::RelationIterator) = keys(cursor.iterator)
Base.getindex(cursor::RelationIterator, element) = RelationIterator(cursor.iterator[element])
function Base.getindex(cursor::RelationIterator, elements...)
    for element in elements
        cursor = getindex(cursor, element)
    end
    return cursor
end

function _show_trie(io, cursor::RelationIterator)
    __show_trie(io, cursor)
end
function __show_trie(io, cursor::RelationIterator; indent=0, num_lines = 10)
    for key in keys(cursor)
        num_lines -= 1
        if num_lines < 0
            print("\n" * " "^indent, "⋮")
            return
        end
        print("\n" * " "^indent, repr(key))
        subcursor = cursor[key]
        #if !isempty(tuples(subcursor))
            print(io, " => ")
            __show_trie(io, subcursor;
                indent = indent + length(key)+1,
                num_lines = num_lines)
        #end
    end
end

function tuples(cursor::RelationIterator)
    len = sum((_num_tuples(r) for r in relations(cursor)), init=0)
    tuples = (
        row_getter(i)[(length(cursor.iterator.prefix)+1):end]
        for r in sort(relations(cursor), by=r->r.relpath)
        for row_getter in (RAI._make_getrow(r.relpath, r.columns),)
        for i in 1:_num_tuples(r)
    )
    return TuplesIterator(len, tuples)
end

function relations(cursor::RelationIterator)
    return PhysicalRelation[
        r_iter.relation
        for r_iter in cursor.iterator.relations
    ]
end

function relations(outputs::JSON3.Array)
    return PhysicalRelation[
        PhysicalRelation(_relpath_from_rel_key(relation.rel_key), relation.columns)
        for relation in outputs
    ]
end
function _relpath_str(rel_key::JSON3.Object)
    name_and_keys = reduce(joinpath, rel_key.keys, init=":$(rel_key.name)")
    name_keys_values = reduce(joinpath, rel_key.values, init=name_and_keys)
    return name_keys_values
end
function _relpath_from_rel_key(rel_key::JSON3.Object)
    return String[":$(rel_key.name)", rel_key.keys..., rel_key.values...]
end


text:
# Concatenation

str1 = "Hello";
str2 = "world";

message = string(str1, " ", str2);
println(message);

message = str1 * " " * str2;
println(message);

# Interpolation

message = "$str1 $str2";
println(message);

println("2 + 8 = $(2 + 8)");

vector = [1, 2, 3];
println("Vector: $vector");

# String compare

println("abcde" < "bcde");
println("abcde" > "aacde");
println("2 + 3 == 5" == "2 + 3 == $(2 + 3)");

# Index search of character in string

str = "abcdefghi";

println(findfirst('d', str));
println(findlast('g', str));
println(findnext("gh", str, 3));
println(findprev('b', str, 4));

# Substring

println(occursin("def", str));
println(occursin('d', str));

# Repeat and join

str = "1 million = " * "1" * repeat("0", 6);
println(str);

str = join(["1", "000", "000"], ",");
println(str);

# Length of a string
println(length(str));

# Length of a substring
println(length(str, 2, 5));

# Size of a string (in bytes)
println(sizeof(str));

# reverse String

println(reverse(str));

# split string

str = "a,b,c,d,e";
println(split(str, ","));

# string replace

str = replace("Python is fast and easy", "Python" => "Julia");
println(str);

replace("Quickly Quicker", "Quick" => "Slow", 1) # Specifies max number of replacements to make

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
import Pkg
Pkg.activate(".")

include("./test/runtests.jl")


text:
module Test08

import Automa
import Automa.RegExp: @re_str
using Test

@testset "Test08" begin
    re = Automa.RegExp
    int = re"[-+]?[0-9]+"
    prefloat = re"[-+]?([0-9]+\.[0-9]*|[0-9]*\.[0-9]+)"
    float = prefloat | re.cat(prefloat | re"[-+]?[0-9]+", re"[eE][-+]?[0-9]+")
    number = int | float
    spaces = re.rep(re.space())
    numbers = re.cat(re.opt(spaces * number), re.rep(re.space() * spaces * number), spaces)

    number.actions[:enter] = [:mark]
    int.actions[:exit]     = [:int]
    float.actions[:exit]   = [:float]

    machine = Automa.compile(numbers)

    actions = Dict(
        :mark  => :(mark = p),
        :int   => :(push!(tokens, (:int, String(data[mark:p-1])))),
        :float => :(push!(tokens, (:float, String(data[mark:p-1])))),
    )

    ctx = Automa.CodeGenContext()

    @eval function tokenize(data)
        tokens = Tuple{Symbol,String}[]
        mark = 0
        $(Automa.generate_init_code(ctx, machine))
        p_end = p_eof = lastindex(data)
        $(Automa.generate_exec_code(ctx, machine, actions))
        return tokens, cs == 0 ? :ok : cs < 0 ? :error : :incomplete
    end

    @test tokenize(b"") == ([], :ok)
    @test tokenize(b"  ") == ([], :ok)
    @test tokenize(b"42") == ([(:int, "42")], :ok)
    @test tokenize(b"3.14") == ([(:float, "3.14")], :ok)
    @test tokenize(b"1 -42 55") == ([(:int, "1"), (:int, "-42"), (:int, "55")], :ok)
    @test tokenize(b"12. -22. .1 +10e12") == ([(:float, "12."), (:float, "-22."), (:float, ".1"), (:float, "+10e12")], :ok)
    @test tokenize(b" -3 -1.2e-3  +54 1.E2  ") == ([(:int, "-3"), (:float, "-1.2e-3"), (:int, "+54"), (:float, "1.E2")], :ok)

    @test tokenize(b"e") == ([], :error)
    @test tokenize(b"42,") == ([], :error)
    @test tokenize(b"42 ,") == ([(:int, "42")], :error)

    @test tokenize(b".") == ([], :incomplete)
    @test tokenize(b"1e") == ([], :incomplete)
    @test tokenize(b"1e-") == ([], :incomplete)
end

end


text:
module Ctrl

# Types
import RIOPA: Config, DataSet, DataSetConfig, DataStreamConfig, PayloadGroup
# Functions
import RIOPA: read_config, default_config
# Modules
import RIOPA: DataGen, IO

function configure_stream(streamCfg::Config)
    payloads = map(
        grpCfg -> PayloadGroup(grpCfg[:size_range], grpCfg[:ratio]),
        streamCfg[:proc_payloads],
    )
    return DataStreamConfig(streamCfg[:name], payloads)
end

function configure_dataset(subCfg::Config)
    compute_seconds = get(subCfg, :compute_seconds, 0)
    datagen_tag_str = get(subCfg, :datagen_backend, nothing)
    io_tag_str = get(subCfg, :io_backend, nothing)
    streams = map(configure_stream, subCfg[:data_streams])
    ds = DataSet(
        DataSetConfig(
            subCfg[:name],
            subCfg[:basename],
            DataGen.get_tag(datagen_tag_str),
            IO.get_tag(io_tag_str),
            subCfg[:nsteps],
            compute_seconds,
            streams,
        ),
    )
    DataGen.initialize_streams!(ds.cfg.datagen_backend_tag, ds)
    return ds
end

generate_data!(ds::DataSet) = DataGen.generate!(ds.cfg.datagen_backend_tag, ds)

function perform_io_step!(ds::DataSet)
    IO.perform_step(ds.cfg.io_backend_tag, ds)
    ds.curr_step += 1
    ds.timestamp = time()
end

function steps_remain(ds::DataSet)
    return ds.curr_step <= ds.cfg.nsteps
end

struct Controller
    datasets::Vector{DataSet}
end

function (c::Controller)()
    cfgtimes = map(ds -> ds.cfg.compute_seconds, c.datasets)
    enabled = map(steps_remain, c.datasets)
    get_reset_time = (i, enable) -> enable ? cfgtimes[i] : Inf
    times = map(((i, enable),) -> get_reset_time(i, enable), enumerate(enabled))
    while any(enabled)
        ta = time()
        dur, idx = findmin(times)
        ds = c.datasets[idx]
        generate_data!(ds)
        diff = time() - ta
        dur = max(0, dur - diff)
        sleep(dur)
        times .-= dur
        perform_io_step!(ds)
        enabled[idx] = steps_remain(ds)
        times[idx] = get_reset_time(idx, enabled[idx])
    end
end

function run(config::Config)
    ctrl = Controller(map(configure_dataset, config[:datasets]))
    ctrl()
end

run(filename::AbstractString) = run(read_config(filename))

run(::Nothing) = run(default_config())

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    ModeAsDefined <: AbstractModeEstimator

Constructors:

    ModeAsDefined()

Get the mode as defined by the density, resp. the underlying distribution
(if available), via `StatsBase.mode`.
"""
struct ModeAsDefined <: AbstractModeEstimator end
export ModeAsDefined

function bat_findmode_impl(target::AnySampleable, algorithm::ModeAsDefined)
    (result = StatsBase.mode(target),)
end

function bat_findmode_impl(target::Distribution, algorithm::ModeAsDefined)
    (result = varshape(target)(StatsBase.mode(unshaped(target))),)
end

function bat_findmode_impl(target::DistributionDensity, algorithm::ModeAsDefined)
    bat_findmode_impl(parent(target), algorithm)
end

"""
    MaxDensitySampleSearch <: AbstractModeEstimator

Constructors:

    MaxDensitySampleSearch()

Estimate the mode as the variate with the highest posterior density value
within a given set of samples.
"""
struct MaxDensitySampleSearch <: AbstractModeEstimator end
export MaxDensitySampleSearch

function bat_findmode_impl(target::DensitySampleVector, algorithm::MaxDensitySampleSearch)
    v, i = _get_mode(target)
    (result = v, mode_idx = i)
end

"""
    MaxDensityNelderMead <: AbstractModeEstimator

Constructors:

```julia
MaxDensityNelderMead(init::InitvalAlgorithm = InitFromTarget)
```

Estimate the mode of a probability density using Nelder-Mead optimization
(via [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)).
"""
@with_kw struct MaxDensityNelderMead{
    TR<:AbstractDensityTransformTarget,
    IA<:InitvalAlgorithm
} <: AbstractModeEstimator
    trafo::TR = NoDensityTransform()
    init::IA = InitFromTarget()
end
export MaxDensityNelderMead

function bat_findmode_impl(target::AnySampleable, algorithm::MaxDensityNelderMead)
    density_notrafo = convert(AbstractDensity, target)
    density, trafo = bat_transform(algorithm.trafo, density_notrafo)
    shape = varshape(density)

    x = collect(unshaped(bat_initval(density, algorithm.init).result))
    r_optim = Optim.maximize(p -> eval_logval(density, p), x, Optim.NelderMead())
    mode_trafo = shape(Optim.minimizer(r_optim.res))
    mode_notrafo = inv(trafo)(mode_trafo)
    (result = mode_notrafo, result_trafo = mode_trafo, trafo = trafo, info = r_optim)
end

"""
    MaxDensityLBFGS <: AbstractModeEstimator

Constructors:

```julia
MaxDensityLBFGS(init::InitvalAlgorithm = InitFromTarget)
```

Estimate the mode of a probability density using LBFGS optimization (via
[Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl)). The gradient
of the density is computed using forward-mode auto-differentiation (via
[ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)).
"""
@with_kw struct MaxDensityLBFGS{
    TR<:AbstractDensityTransformTarget,
    IA<:InitvalAlgorithm
} <: AbstractModeEstimator
    trafo::TR = PriorToGaussian()
    init::IA = InitFromTarget()
end
export MaxDensityLBFGS

function bat_findmode_impl(target::AnySampleable, algorithm::MaxDensityLBFGS)
    density_notrafo = convert(AbstractDensity, target)
    density, trafo = bat_transform(algorithm.trafo, density_notrafo)
    shape = varshape(density)

    x = collect(unshaped(bat_initval(density, algorithm.init).result))
    r_optim = Optim.maximize(p -> eval_logval(density, p), x, Optim.LBFGS(); autodiff = :forward)
    mode_trafo = shape(Optim.minimizer(r_optim.res))
    mode_notrafo = inv(trafo)(mode_trafo)
    (result = mode_notrafo, result_trafo = mode_trafo, trafo = trafo, info = r_optim)
end

function bat_marginalmode_impl(samples::DensitySampleVector; nbins::Union{Integer, Symbol} = 200)

    shape = varshape(samples)
    flat_samples = flatview(unshaped.(samples.v))
    n_params = size(flat_samples)[1]
    nt_samples = ntuple(i -> flat_samples[i,:], n_params)
    marginalmode_params = Vector{Float64}()

    for param in Base.OneTo(n_params)
        if typeof(nbins) == Symbol
            number_of_bins = _auto_binning_nbins(nt_samples, param, mode=nbins)
        else
            number_of_bins = nbins
        end

        marginalmode_param = find_marginalmodes(get_marginal_dist(samples, param, bins=number_of_bins).result)

        if length(marginalmode_param[1]) > 1
            @warn "More than one bin with the same weight is found. Returned the first one"
        end
        push!(marginalmode_params, marginalmode_param[1][1])
    end
    (result = shape(marginalmode_params),)
end

# From Plots.jl, original authors Oliver Schulz and Michael K. Borregaard
function _auto_binning_nbins(vs::NTuple{N,AbstractVector}, dim::Integer; mode::Symbol = :auto) where N
    max_bins = 10_000
    _cl(x) = min(ceil(Int, max(x, one(x))), max_bins)
    _iqr(v) = (q = quantile(v, 0.75) - quantile(v, 0.25); q > 0 ? q : oftype(q, 1))
    _span(v) = maximum(v) - minimum(v)

    n_samples = length(LinearIndices(first(vs)))

    # The nd estimator is the key to most automatic binning methods, and is modified for twodimensional histograms to include correlation
    nd = n_samples^(1/(2+N))
    nd = N == 2 ? min(n_samples^(1/(2+N)), nd / (1-cor(first(vs), last(vs))^2)^(3//8)) : nd # the >2-dimensional case does not have a nice solution to correlations

    v = vs[dim]

    if mode == :auto
        mode = :fd
    end

    if mode == :sqrt  # Square-root choice
        _cl(sqrt(n_samples))
    elseif mode == :sturges  # Sturges' formula
        _cl(log2(n_samples) + 1)
    elseif mode == :rice  # Rice Rule
        _cl(2 * nd)
    elseif mode == :scott  # Scott's normal reference rule
        _cl(_span(v) / (3.5 * std(v) / nd))
    elseif mode == :fd  # Freedman–Diaconis rule
        _cl(_span(v) / (2 * _iqr(v) / nd))
    elseif mode == :wand
        _cl(wand_edges(v))  # this makes this function not type stable, but the type instability does not propagate
    else
        error("Unknown auto-binning mode $mode")
    end
end


text:
using BenchmarkTools

import Algopt.LocalDescent

using Algopt.TestFunctions: Rosenbrock

ros = Rosenbrock(a = 1, b = 5)

ls = LocalDescent.LineSearch()
res = @btime LocalDescent.search(ls, ros.f, [-2, -2], [1, 1])
println("LineSearch: ", res)

sbt = LocalDescent.StrongBacktracking()
res = @btime LocalDescent.search(sbt, ros.f, ros.∇f, [-2, -2], [1, 1])
println("StrongBacktracking: ", res)


text:
using Glimmer, Glimmer.FlexUI
using Test

@testset "creating an application" begin
    app = App()
    prop!(app, :title, "Glimmer Example - Controls Gallery")

    labelText = addVariable!(app, Variable(name="labelText", type="string",value="This text comes from a variable"))
    firstName = addVariable!(app, Variable(name="firstName", type="string",value=""))
    height = addVariable!(app, Variable(name="height", type="flota64",value="176"))

    option = addVariable!(app, Variable(name="option", type="string",value="Four"))

    range1 = addVariable!(app, Variable(name="range1", type="flota64",value=0))
    range2 = addVariable!(app, Variable(name="range2", type="flota64",value=0.5))

    ui = VContainer(
        Accordion(
            panels = [
                ExpansionPanel(
                    title="Label",
                    subtitle="Examples of different labels",
                    content = VContainer(
                        HContainerSpace(
                            Label(text="A Normal Label"),
                            Label(text="""Label(text="A Normal Label")"""),
                        ),
                        HContainerSpace(
                            Label(text="Label with style", style="color: red;"),
                            Label(text="""Label(text="Label with style", style="color: red;")"""),
                        ),
                        HContainerSpace(
                            Label(text="Label with style and class H1", class="h1", style="color: red;"),
                            Label(text="""Label(text="Label with style and class H1", class="h1", style="color: red;")"""),
                        ),
                        HContainerSpace(
                            Label(text="Label with style and class H2", class="h2", style="color: blue;"),
                            Label(text="""Label(text="Label with style and class H2", class="h2", style="color: blue;")"""),
                        ),
                        HContainerSpace(
                            Label(text="Label with style and class H1", class="h3", style="background-color:rgb(255, 99, 71);"),
                            Label(text="""Label(text="Label with style and class H1", class="h3", style="background-color:rgb(255, 99, 71);")"""),
                        ),
                        HContainerSpace(
                            Label(text="Label with style and class H1", class="h1", style="background-color:#ff6347;"),
                            Label(text="""Label(text="Label with style and class H1", class="h1", style="background-color:#ff6347;")"""),
                        ),
                        HContainerSpace(
                            Label(text="Allow HTML tags in text.", isHTML=true),
                            Label(text="""Label(text="Allow HTML tags in text.", isHTML=true)"""),
                        ),
                        Divider(),
                        HContainerSpace(
                            Label(variable="labelText", isHTML=true),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Label(variable="labelText", isHTML=true)"""),
                                    Label(text="""You can change the content of the label using the Julia REPL:"""),
                                    Label(text="""REPL: Glimmer.Example.labelText[] = "This text is from Julia" """, style="color: red;" ),
                                ]
                            ),
                        ),
                    ),
                ),

                ExpansionPanel(
                    title="Field",
                    subtitle="Examples of different fields",
                    content = VContainer(
                        HContainerSpace(
                            Field(
                                input="text",
                                label="First Name",
                                hint="No digits please",
                                variable="firstName",
                            ),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Field(
                                        input="text",
                                        label="First Name",
                                        hint="No digits please",
                                        variable="firstName",
                                    )"""),
                                    Label(text="""REPL: Glimmer.Example.firstName[] = "This text is from Julia" """, style="color: red;" ),
                                ]
                            ),
                        ),
                        HContainerSpace(
                            Field(
                                input="number",
                                label="Height",
                                hint="in centimeters",
                                variable="height",
                            ),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Field(
                                        input="number",
                                        label="Height",
                                        hint="in centimeters",
                                        variable="height",
                                    )"""),
                                    Label(text="""REPL: Glimmer.Example.height[] = 183.1 """, style="color: red;" ),
                                ]
                            ),
                        ),

                        HContainerSpace(
                            Field(
                                input="select",
                                label="A Simple ComboBox",
                                hint="Please select an option",
                                variable="option",
                                options = [
                                    Dict(:key=>"One", :value=>"First Option"),
                                    Dict(:key=>"Two", :value=>"Second Option"),
                                    Dict(:key=>"Three", :value=>"Third Option"),
                                    Dict(:key=>"Four", :value=>"Forth Option"),
                                ]
                            ),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Field(
                                        input="select",
                                        label="A Simple ComboBox",
                                        hint="Please select an option",
                                        variable="option",
                                        options = [
                                            Dict(:key=>"One", :value=>"First Option"),
                                            Dict(:key=>"Two", :value=>"Second Option"),
                                            Dict(:key=>"Three", :value=>"Third Option"),
                                            Dict(:key=>"Four", :value=>"Forth Option"),
                                        ]
                                    )"""),
                                    Label(text="""REPL: Glimmer.Example.option[] = "Two" """, style="color: red;" ),
                                ]
                            ),
                        ),
                    ),
                ),

                ExpansionPanel(
                    title="Slider",
                    subtitle="Examples of different sliders",
                    content = VContainer(
                        HContainerSpace(
                            Slider(
                                min=0,
                                max=100,
                                value=10,
                                variable="range1"
                            ),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Slider(
                                        min=0,
                                        max=100,
                                        value=10,
                                        variable="range1"
                                    )"""),
                                ]
                            ),
                        ),
                        HContainerSpace(
                            Slider(
                                text="Some Text",
                                min=0,
                                max=100,
                                value=10,
                                variable="range1"
                            ),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Slider(
                                        text="Some Text",
                                        min=0,
                                        max=100,
                                        value=10,
                                        variable="range1"
                                    )"""),
                                ]
                            ),
                        ),
                        HContainerSpace(
                            Slider(
                                text="Embedded Value [\$()]",
                                trailingText="[\$()cm]",
                                min=0,
                                max=100,
                                value=10,
                                variable="range1"
                            ),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Slider(
                                        text="Embedded Value [\$()]",
                                        trailingText="[\$()cm]",
                                        min=0,
                                        max=100,
                                        value=10,
                                        variable="range1"
                                    )"""),
                                ]
                            ),
                        ),

                        HContainerSpace(
                            Slider(
                                text="Step Size",
                                trailingText="[\$()cm]",
                                min=0,
                                max=1,
                                value=0.5,
                                step=0.01,
                                variable="range2"
                            ),
                            Container(
                                direction = "column",
                                align = "center end",
                                children = [
                                    Label(text="""Slider(
                                        text="Step Size",
                                        trailingText="[\$()cm]",
                                        min=0,
                                        max=1,
                                        value=0.5,
                                        step=0.01,
                                        variable="range2"
                                    )"""),
                                ]
                            ),
                        ),
                    ),
                )
            ],
        ),

    )
    FlexUI.controls!(app, ui)
end


text:
# Represent evenly spaced points on [0,2π)
struct FourierGrid{T} <: AbstractVector{T}
    n::Int
    π_over_n::T
    FourierGrid{T}(n::Int) where T = new{T}(n, convert(T,π)/n)
end

size(g::FourierGrid) = (g.n,)
@inline function getindex(g::FourierGrid, k::Int)
    @boundscheck Base.checkbounds(g, k)
    (2k-2)*g.π_over_n
end
Base.IndexStyle(::Type{<:FourierGrid}) = IndexLinear()

FourierGrid(10)

using BenchmarkTools
g = FourierGrid{Float64}(10_000)

fp(n) = collect(FourierGrid{Float64}(n))
@btime collect(fp(10_000))
@btime fourierpoints(10_000)

n=10
T=Float64
    convert(T,π)*(0:2:2n-2)/n - [convert(T,π)*(2k-2)/n for k=1:n]


text:
# Autogenerated wrapper script for LERC_jll for x86_64-w64-mingw32-cxx11
export libLerc

JLLWrappers.@generate_wrapper_header("LERC")
JLLWrappers.@declare_library_product(libLerc, "libLerc.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libLerc,
        "bin\\libLerc.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
###################################################
#
# Copyright © Akamai Technologies. All rights reserved.
#
# File: Domain.jl
#
# Functions to communicate with the mPulse Repository REST API regarding Domain Objects.
# This file MUST be `include()`d from `mPulseAPI.jl`
#
###################################################

export
    getRepositoryDomain

"""
Fetches a Domain object from the mPulse repository

To fetch a single domain, at least one of `domainID`, `appKey` or `appName` must be passed in to identify the domain.
If none of these are passed in, then all domains that are readable by the specified `token` will be returned as an array.

The domain will be cached in memory for 1 hour, so subsequent calls using a matching `domainID`, `appKey` or `appName` return
quickly without calling out to the API.  This can be a problem if the domain changes in the repository.
You can clear the cache for this domain using [`mPulseAPI.clearDomainCache`](@ref) and passing in one of `domainID`, `appKey` or `appName`.

#### Arguments
`token::AbstractString`
:    The Repository authentication token fetched by calling [`getRepositoryToken`](@ref)

#### Keyword Arguments
`domainID::Int64`
:    The ID of the domain to fetch.  This is the fastest method, but it can be hard to figure out a domain's ID

`appKey::AbstractString`
:    The App Key (formerly known as API key) associated with the domain.  This is available from the mPulse domain configuration dialog.

`appName::AbstractString`
:    The App name in mPulse. This is available from the mPulse domain configuration dialog.

#### Returns
`{Dict|Array{Dict}}` If one of `domainID`, `appKey` or `appName` are passed in, then a single `domain` object is returned as a `Dict`.

If none of these are passed in, then an array of all domains is returned, each is a `Dict`.

The `domain` `Dict` has the following fields:

`name`
:    The app's name

`id::Int64`
:    The app's ID

`body::XMLElement`
:    An XML object representing the app's XML definition

`tenantID::Int64`
:    The ID of the tenant that this app is in

`description::AbstractString`
:    The description of this app entered into mPulse

`created::DateTime`
:    The timestamp when this object was created

`lastModified::DateTime`
:    The timestamp when this object was created

`attributes::Dict`
:    A `Dict` of attributes for this app, including its `AppKey`

`custom_metrics::Dict`
:    A $(mPulseAPI.readdocs("CustomMetricMap-structure", indent=5))

`custom_timers::Dict`
:    A $(mPulseAPI.readdocs("CustomTimerMap-structure", indent=5))

`session_timeout::Int64`
:    The session timeout value in minutes

`resource_timing::Bool`
:    Flag indicating whether resource timing collection is enabled or not

`vertical_market::AbstractString`
:    The vertical market that this domain belongs to

#### Throws
`ArgumentError`
:    if token is empty or domainID, appKey and appName are all empty

`mPulseAPIException`
:    if API access failed for some reason

`Exception`
:    if something unexpected happened while parsing the repository object

"""
function getRepositoryDomain(token::AbstractString; domainID::Int64=0, appKey::AbstractString="", appName::AbstractString="", appID::AbstractString="")
    # Keep appID for backwards compatibility
    if isempty(appKey) && !isempty(appID)
        appKey = appID
    end

    domain_list = getRepositoryObject(
                token,
                "domain",
                Dict{Symbol, Any}(:id => domainID, :apiKey => appKey, :name => appName),
                filterRequired=false
        )

    # Always convert to an array for easier processing
    if !isa(domain_list, Array)
        domain_list = [domain_list]
    end

    for domain in domain_list
        # If the object came out of cache, then these fields have already been populated
        if !haskey(domain, "custom_metrics")
            try
                domain["custom_metrics"]  = getCustomMetricMap(domain)
                domain["custom_timers"]   = getCustomTimerMap(domain)
                domain["session_timeout"] = getNodeContent(domain, "SessionTimeout", 30)
                domain["resource_timing"] = getNodeContent(domain, "CollectResources", false)
                domain["vertical_market"] = getNodeContent(domain, "VerticalMarket", "")
            catch ex
                # If this is an Exception that we are not prepared to deal with
                if !isa(ex, LightXML.XMLParseError)
                    rethrow()
                end
            end
        end

        domain["custom_dimensions"] = "Custom Dimension map for a beacon is stored in the beacon and is more accurate and timely than the repository"

        domain["dswb_table_name"]   = "beacons_$(domain["id"])"

        if haskey(domain, "attributes") && haskey(domain["attributes"], "apiKey")
            domain["attributes"]["appID"] = domain["attributes"]["apiKey"]
            domain["attributes"]["appKey"] = domain["attributes"]["apiKey"]
        end

        delete!(domain, "readOnly")
    end

    # Return the first element only if the caller asked for a unique domain, else
    # return the list even if it only has one element in it
    if domainID != 0 || appKey != "" || appName != ""
        return domain_list[1]
    else
        return domain_list
    end
end


text:
module XtalsPyTools

using Xtals, PyCall, FIGlet, Printf, Graphs, UUIDs

# - lists python dependencies
# - checks for python dependencies at pre-compile, warns on failure
# - provides function for python dependency instantiation
include("pydeps.jl")

include("infer_geometry_based_bonds!.jl")
include("primitive_cell.jl")
include("misc.jl")

function __init__()
    # load Python dependencies
    init_pydeps()
    # check Python dependencies
    for pydep ∈ PYDEPS
        check_pydep(pydep)
    end
end

export infer_geometry_based_bonds!, primitive_cell, rc

end

max_stars_repo_name: SimonEnsemble/XtalsPyTools

max_issues_repo_name: SimonEnsemble/XtalsPyTools


text:
using MDP
using Test

P_01 = Matrix[ [ 0.5  0.5
                 0.0  1.0 ],
               [ 0.0  1.0
                 0.5  0.5 ] ]

r_01 = [ 1 2 ; 3 4 ; 5 6]
@test_throws ErrorException ProbModel(r_01, P_01; objective = :Max)

P_02 = Matrix[ [ 0.5  0.4999
                 0.0  1.0 ],
               [ 0.0  1.0
                 0.5  0.5 ] ]

r_02 = [ 1 2 ; 3 4]

@test_throws ErrorException ProbModel(r_02, P_02; objective = :Max)


text:
function gettril(x::Array{Float64, 2})
  n, m = size(x)
  a = zeros(convert(Int, n.*(n+1)/2))
  k::Int = 1
  for i = 1:n
    for j = 1:i
      a[k] = x[i, j]
      k += 1
    end
  end
  return a
end


text:
using GNSSTools
using Statistics
using ProgressMeter
using JLD
using PyPlot
include("ch4_constellation_hist_structs.jl")
pygui(true)

function get_median_doppler(dopplers, doppler_rates, max_elevations,
                            elev_range; nbins=100)
    elev_min = minimum(elev_range)
    elev_max = maximum(elev_range)
    idxs = (max_elevations .>= elev_min) .== (max_elevations .< elev_max)
    doppler_elev = dopplers[idxs]
    doppler_rate_elev = doppler_rates[idxs]
    x = get_histogram((doppler_rate_elev, doppler_elev); nbins=nbins)
    doppler_elev_mean = x.edges[2][1:end-1] .+ (x.edges[2][2] - x.edges[2][1])/2
    doppler_rate_mean = Array{Float64}(undef, length(doppler_elev_mean))
    for i in 1:size(x.weights)[2]
        # y = x.weights[:,i]
        # edges = x.edges[1][1:end-1]
        # idxs = y .> 0
        # if any(idxs)
        #     doppler_rate_mean[i] = median(edges[y .> 0])
        # else
        #     doppler_rate_mean[i] = NaN
        # end
        doppler_rate_mean[i] = sum(x.weights[:,i].*x.edges[1][1:end-1])/sum(x.weights[:,i])
    end
    return (x, doppler_elev, doppler_rate_elev,
            doppler_elev_mean, doppler_rate_mean)
end

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
# directory = ""
doppler_file_name = string(directory, "ch1_doppler_curves.jld")
doppler_raw_gps, doppler_rate_raw_gps,
elevations_gps, sig_freq, max_elevations_gps = load(
                                                       doppler_file_name,
                                                       "doppler_raw_gps",
                                                       "doppler_rate_raw_gps",
                                                       "elevations_gps",
                                                       "freq",
                                                       "max_elevations_gps")
doppler_raw_iridium, doppler_rate_raw_iridium,
elevations_iridium, max_elevations_iridium = load(
                                                     doppler_file_name,
                                                     "doppler_raw_iridium",
                                                     "doppler_rate_raw_iridium",
                                                     "elevations_iridium",
                                                     "max_elevations_iridium")
doppler_raw_starlink, doppler_rate_raw_starlink,
elevations_starlink, max_elevations_starlink = load(
                                                       doppler_file_name,
                                                       "doppler_raw_starlink",
                                                       "doppler_rate_raw_starlink",
                                                       "elevations_starlink",
                                                       "max_elevations_starlink")
doppler_raw_oneweb, doppler_rate_raw_oneweb,
elevations_oneweb, max_elevations_oneweb = load(doppler_file_name,
                         "doppler_raw_oneweb",
                         "doppler_rate_raw_oneweb",
                         "elevations_oneweb",
                         "max_elevations_oneweb")

gps_hist = initiate_hist(doppler_raw_gps, doppler_rate_raw_gps,
                         elevations_gps, max_elevations_gps;
                         nbins=1000, elev_bins=100)
iridium_hist = initiate_hist(doppler_raw_iridium, doppler_rate_raw_iridium,
                             elevations_iridium, max_elevations_iridium;
                             nbins=1000)
starlink_hist = initiate_hist(doppler_raw_starlink, doppler_rate_raw_starlink,
                              elevations_starlink, max_elevations_starlink;
                              nbins=1000)
oneweb_hist = initiate_hist(doppler_raw_oneweb, doppler_rate_raw_oneweb,
                            elevations_oneweb, max_elevations_oneweb;
                            nbins=1000)

# elevation_bins = Array(range(1, 91, step=2))
# gps_doppler_max, gps_doppler_max = max_doppler_curve_per_elev(
#                                                     gps_hist, elevation_bins)
# iridium_doppler_max, iridium_doppler_rate_max = max_doppler_curve_per_elev(
#                                                     iridium_hist, elevation_bins)
# starlink_doppler_max, starlink_doppler_rate_max = max_doppler_curve_per_elev(
#                                                     starlink_hist, elevation_bins)
# oneweb_doppler_max, oneweb_doppler_rate_max = max_doppler_curve_per_elev(
#                                                     oneweb_hist, elevation_bins)

elev_ranges = [[10, 11],
               [15, 16],
               [30, 31],
               [45, 46],
               [60, 61],
               [90, 89]]
nbins = 100
cmap = "gray_r"
fig = figure(figsize=(7.5, 5))
ax1 = fig.add_subplot(2,2,1)
hist2D(gps_hist.dopplers./1000, gps_hist.doppler_rates, bins=nbins,
          cmap=cmap, density=true, vmax=1)
# for i in 1:length(elev_ranges)
#     elev_min = elev_ranges[i][1]
#     elev_max = elev_ranges[i][2]
#     x, d, dd, d_mean, dd_mean = get_median_doppler(gps_hist.dopplers,
#                                                    gps_hist.doppler_rates,
#                                                    gps_hist.max_elevations,
#                                                    [elev_min, elev_max];
#                                                    nbins=(500, 100))
#     ax1.plot(d_mean./1000, dd_mean, label="$(elev_min)ᵒ",
#              linewidth=1)
# end
xlabel("Doppler (kHz)")
ylabel("Doppler Rate (Hz/s)")
# legend()
title("(a)")

ax2 = fig.add_subplot(2,2,2, aspect="auto")
ax2.hist2d(iridium_hist.dopplers./1000, iridium_hist.doppler_rates, bins=nbins,
            cmap=cmap, density=true, vmax=0.0005)
for i in 1:length(elev_ranges)
    elev_min = elev_ranges[i][1]
    elev_max = elev_ranges[i][2]
    x, d, dd, d_mean, dd_mean = get_median_doppler(iridium_hist.dopplers,
                                                   iridium_hist.doppler_rates,
                                                   iridium_hist.max_elevations,
                                                   [elev_min, elev_max];
                                                   nbins=(500, 100))
    ax2.plot(d_mean./1000, dd_mean, label="$(elev_min)°",
             linewidth=1)
end
xlabel("Doppler (kHz)")
ylabel("Doppler Rate (Hz/s)")
legend(bbox_to_anchor=(1, 1.04), loc="upper left", ncol=1)
title("(b)")

ax3 = fig.add_subplot(2,2,3, aspect="auto")
ax3.hist2d(starlink_hist.dopplers./1000, starlink_hist.doppler_rates, bins=nbins,
           cmap=cmap, density=true, vmax=0.0005)
for i in 1:length(elev_ranges)
    elev_min = elev_ranges[i][1]
    elev_max = elev_ranges[i][2]
    x, d, dd, d_mean, dd_mean = get_median_doppler(starlink_hist.dopplers,
                                                   starlink_hist.doppler_rates,
                                                   starlink_hist.max_elevations,
                                                   [elev_min, elev_max];
                                                   nbins=(500, 100))
    ax3.plot(d_mean./1000, dd_mean,
             linewidth=1)
end
xlabel("Doppler (kHz)")
ylabel("Doppler Rate (Hz/s)")
# legend()
title("(c)")

ax4 = fig.add_subplot(2,2,4, aspect="auto")
ax4.hist2d(oneweb_hist.dopplers./1000, oneweb_hist.doppler_rates, bins=nbins,
           cmap=cmap, density=true, vmax=0.0005)
for i in 1:length(elev_ranges)
    elev_min = elev_ranges[i][1]
    elev_max = elev_ranges[i][2]
    x, d, dd, d_mean, dd_mean = get_median_doppler(oneweb_hist.dopplers,
                                                   oneweb_hist.doppler_rates,
                                                   oneweb_hist.max_elevations,
                                                   [elev_min, elev_max];
                                                   nbins=(500, 100))
    ax4.plot(d_mean./1000, dd_mean,
             linewidth=1)
end
xlabel("Doppler (kHz)")
ylabel("Doppler Rate (Hz/s)")
# legend()
title("(d)")

subplots_adjust(wspace=0.4, hspace=0.5,
                bottom=0.1, left=0.12, right=0.88, top=0.93)

savefig(string(directory, "figures/ch4_joint_doppler-doppler_rate_distributions.pdf"), dpi=300)

fig = figure()
ax = fig.add_subplot(1,1,1)
elev_ranges = [[10, 10.5],
               [15, 15.1],
               [30, 30.1],
               [45, 45.1],
               [60, 60.1],
               [90, 89.9]]
ax.hist2d(doppler_raw_iridium./1000, doppler_rate_raw_iridium, bins=100,
          cmap="gray_r", density=true, vmax=0.0005)
for i in 1:length(elev_ranges)
    elev_min = elev_ranges[i][1]
    elev_max = elev_ranges[i][2]
    x, d, dd, d_mean, dd_mean = get_median_doppler(doppler_raw_iridium,
                                                   doppler_rate_raw_iridium,
                                                   max_elevations_iridium,
                                                   [elev_min, elev_max];
                                                   nbins=(500, 100))
    ax.plot(d_mean./1000, dd_mean, label="$(elev_min)ᵒ",
            linewidth=1)
end
xlabel("Doppler (kHz)")
ylabel("Doppler Rate (Hz/s)")
legend()

hist_data = starlink_hist
fig = figure(figsize=(7.5, 4))
nbins = 100
max_elev_lim = 10
max_elev_idxs = hist_data.max_elevations .>= max_elev_lim
ax1 = fig.add_subplot(2, 3, 1)
idxs = (hist_data.elevations .>= 10) .== (hist_data.elevations .< 15) .== max_elev_idxs
ax1.hist(hist_data.dopplers[idxs]./1000,
         color="k", bins=nbins, density=true);
xlabel("Doppler (kHz)")
ylabel("Probability")
xlim([minimum(hist_data.dopplers), maximum(hist_data.dopplers)]./1000)
ax1.ticklabel_format(axis="y", style="sci", scilimits=[-1,1])
title("(a)")

ax2 = fig.add_subplot(2, 3, 2)
idxs = (hist_data.elevations .>= 60) .== (hist_data.elevations .< 65) .== max_elev_idxs
ax2.hist(hist_data.dopplers[idxs]./1000,
         color="k", bins=nbins, density=true);
xlabel("Doppler (kHz)")
ylabel("Probability")
xlim([minimum(hist_data.dopplers), maximum(hist_data.dopplers)]./1000)
ax2.ticklabel_format(axis="y", style="sci", scilimits=[-1,1])
title("(b)")

ax3 = fig.add_subplot(2, 3, 3)
idxs = (hist_data.elevations .>= 85) .== (hist_data.elevations .< 90) .== max_elev_idxs
ax3.hist(hist_data.dopplers[idxs]./1000,
color="k", bins=nbins, density=true);
xlabel("Doppler (kHz)")
ylabel("Probability")
xlim([minimum(hist_data.dopplers), maximum(hist_data.dopplers)]./1000)
ax3.ticklabel_format(axis="y", style="sci", scilimits=[-1,1])
title("(c)")

ax1 = fig.add_subplot(2, 3, 4)
idxs = (hist_data.elevations .>= 10) .== (hist_data.elevations .< 15) .== max_elev_idxs
ax1.hist(hist_data.doppler_rates[idxs],
         color="k", bins=nbins, density=true);
xlabel("Doppler Rate (Hz/s)")
ylabel("Probability")
xlim([minimum(hist_data.doppler_rates), maximum(hist_data.doppler_rates)])
ax1.ticklabel_format(axis="y", style="sci", scilimits=[-1,1])
title("(d)")

ax2 = fig.add_subplot(2, 3, 5)
idxs = (hist_data.elevations .>= 60) .== (hist_data.elevations .< 65) .== max_elev_idxs
ax2.hist(hist_data.doppler_rates[idxs],
         color="k", bins=nbins, density=true);
xlabel("Doppler Rate (Hz/s)")
ylabel("Probability")
xlim([minimum(hist_data.doppler_rates), maximum(hist_data.doppler_rates)])
ax2.ticklabel_format(axis="y", style="sci", scilimits=[-1,1])
title("(e)")

ax3 = fig.add_subplot(2, 3, 6)
idxs = (hist_data.elevations .>= 85) .== (hist_data.elevations .< 90) .== max_elev_idxs
ax3.hist(hist_data.doppler_rates[idxs],
color="k", bins=nbins, density=true);
xlabel("Doppler Rate (Hz/s)")
ylabel("Probability")
xlim([minimum(hist_data.doppler_rates), maximum(hist_data.doppler_rates)])
ax3.ticklabel_format(axis="y", style="sci", scilimits=[-1,1])
title("(f)")

subplots_adjust(wspace=0.36, hspace=0.6, bottom=0.115,
                left=0.07, right=0.93, top=0.94)
savefig(string(directory, "figures/ch4_starlink_hist_at_different_elevations_above_$(max_elev_lim)_degrees_elevation.pdf"), dpi=300)

######### Skyplots for diff passes with max elevations
# Iridium
δt = 1  # seconds
total_time = 1*24*60*60  # seconds
t_range = Array(0:δt:total_time);
min_elevation = 0  # degrees
user_lla = (40.01, -105.2437, 1655)
a = Rₑ+780*1000;  # meters
plane_num = 6;
sat_per_plane = 11;
incl = 86;
ΔΩ = 360/plane_num/2;
Δf_per_plane = 360/sat_per_plane/2;  # degrees
orbital_period = 2π*sqrt(a^3 / GM)  # seconds
iridium_constellation= define_constellation(a, plane_num,
                                            sat_per_plane, incl,
                                            t_range; obs_lla=user_lla,
                                            show_plot=false,
                                            print_steps=true,
                                            ΔΩ=ΔΩ,
                                            Δf_per_plane=Δf_per_plane);

elev_azimuths = []
elev_elevations = []
elev_max_elevations = []
for satellite in iridium_constellation.satellites
    ts = Array{Float64}(undef, length(satellite.t))
    ts[1] = 0.0
    azimuths = Array{Float64}(undef, length(satellite.t))
    elevations = Array{Float64}(undef, length(satellite.t))
    max_elevations = Array{Float64}(undef, length(satellite.t))
    k_start = 1
    k = 1
    pass = 1
    pass_num = 3
    for j in 1:length(satellite.t)
        sat_range, azimuth, elevation = calcelevation(satellite.r_ecef[j,:],
                                                        user_lla)
        if k > 1
            if elevation > min_elevation
                ts[k] = satellite.t[j]
                azimuths[k] = azimuth
                elevations[k] = elevation
                Δt = abs(ts[k] - ts[k-1])*(60*60*24)
                if Δt > 2*δt
                    # if pass == pass_num
                        max_elevations[k_start:k-1] .= maximum(elevations[k_start:k-1])
                        k_start = k
                    # else
                    #     max_elevations[k_start:j-1] .= -1
                    #     k_start = k
                    # end
                    pass += 1
                end
                k += 1
            end
        else
            if elevation > min_elevation
                ts[k] = satellite.t[j]
                azimuths[k] = azimuth
                elevations[k] = elevation
                k += 1
            end
        end
    end
    if k_start < k
        max_elevations[k_start:k-1] .= -1
    end
    push!(elev_azimuths, azimuths[1:k-1])
    push!(elev_elevations, elevations[1:k-1])
    push!(elev_max_elevations, max_elevations[1:k-1])
end

δ_deg = 1
elev_ranges = [[10, 10+δ_deg],
               [15, 15+δ_deg],
               [30, 30+δ_deg],
               [45, 45+δ_deg],
               [60, 60+δ_deg],
               [90, 90-δ_deg]]
fig = figure(figsize=(7.5, 4))
colors = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
          "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf"]
for i in 1:length(elev_ranges)
    ax = fig.add_subplot(2, 3, i, polar=true)
    color = colors[i]
    elev_min = minimum(elev_ranges[i])
    elev_max = maximum(elev_ranges[i])
    for j in 1:length(elev_max_elevations)
        max_els = elev_max_elevations[j]
        idxs = (max_els .>= elev_min) .== (max_els .<= elev_max)
        if any(idxs)
            ax.plot(elev_azimuths[j][idxs].*(π/180),
                       90 .- elev_elevations[j][idxs],
                       linestyle="-", linewidth=1, color=color)
            break
        end
    end
    ax.set_theta_zero_location("N")
    # Reverse direction of azimuth to clockwise
    ax.set_theta_direction(-1)
    # Reverse y-axis so that 90 degrees elevation is at center
    ax.set_rlim(0, 90, 1)
    ax.set_yticks(range(0, 91, step=30))
    fontsize = 8
    ax.set_yticklabels(reverse(ax.get_yticks()), fontsize=fontsize-2)
    ax.set_xticklabels(["N", "", "E", "", "S", "", "W", ""], fontsize=fontsize)
    # ylabel("Elevation (deg)")
    letter = 'a' + (i-1)
    title("($letter)")
    # title("$(elev_ranges[i][1])°")
end
subplots_adjust(wspace=0.7, hspace=0.8, bottom=0.07,
                left=0.02, right=0.98, top=0.88)
savefig(string(directory, "figures/ch4_diff_max_elevation_passes.pdf"), dpi=300)


text:
#THis module defines the lattice geometry, namely the bonds
#in the lattice with or without periodic boundary condition (Pbc)
#for the Heisenberg model in a L_x x L_y lattice

module Lattice

#Initialization of the lattice (store all the bonds in the Heisenberg Hamiltonian)
function InitLattice(L_x::Int64, L_y::Int64, nspins::Int64)

    #Define arrays that will be filled with nearest neigbourng bonds
    v_bonds= Array{Array{Int64}}(undef,0)
    h_bonds= Array{Array{Int64}}(undef,0)

    #Set boundary conditions
    if L_x > 2 && L_y > 2
        pbc_v = true; pbc_h = true
    elseif L_x <= 2 && L_y > 2
        pbc_v = true; pbc_h = false
    elseif L_x > 2 && L_y == 2
        pbc_v = false; pbc_h = true
    else
        pbc_v = false; pbc_h = false
    end

    #Check if the number of spins is compatible with the given lattice size
    if (L_x*L_y != nspins)
        println("Error , the number of spins is not compabitle with lattice length ")
        error();
    end;

    nn = Array{Array}(undef,nspins)
    for i=1:nspins
        nn[i]= Array{Int64}(undef,4)
    end;

    #Defining the nearest-neighbors
    for i=1:nspins

        if pbc_h
            nn[i][1] = PbcH(L_x,i-1,i)
        else nn[i][1] = PbcH(L_x,i-1,i)
        end

        if pbc_h
            nn[i][2] = PbcH(L_x,i+1,i)
        else nn[i][2] = PbcH(L_x,i+1,i)
       end

       if pbc_v
           nn[i][3] = PbcV(nspins,i-L_x)
       else nn[i][3] = V(nspins,i-L_x)
       end

       if pbc_v
           nn[i][4] = PbcV(nspins,i+L_x)
       else nn[i][4]= V(nspins,i+L_x)
       end
   end

 #Push nearest neighbourong bonds along y into the array v_bonds
   for  i=1:nspins
       for  k=3:4
           j = nn[i][k]
           if i<j
               push!(v_bonds,[i;j])
           end
       end
   end

  #Push nearest neighbourong bonds along x into the array h_bonds
   for  i=1:nspins
       for  k=1:2
           j = nn[i][k];
           if i<j
               push!(h_bonds,[i;j])
           end
       end
   end

   #We only keep inequivalent bonds
    bondsH = unique(h_bonds)
    bondsV = unique(v_bonds)

   return bondsH, bondsV

end

#Small functions to set up the lattice Horizontal (x) Pbc
function PbcH(L_x,nn,s)
    if (s-1)%L_x==0 && nn== (s-1)
      return (s+L_x-1);
  elseif (s)%L_x==0 && nn==(s+1)
      return (s-L_x+1);
    else
      return nn;
  end;
end;

#Small functions to set up the lattice Vertical (y) Pbc
function PbcV(nspins,nn)
    if nn>nspins
      return (nn-nspins);
  elseif nn<=0
      return (nspins+nn);
    else
      return nn;
  end;
end;

#Small functions to set up the lattice Horizontal (x) without Pbc
function H(L_x,nn,s)
    if s%L_x==0 && nn==(s-1)
      return -1
  elseif (s+1)%L_x==0 && nn==(s+1)
      return -1
    else
      return nn
  end;
end;

#Small functions to set up the lattice Vertical (y) without Pbc
function V(nspins,nn)
    if nn>nspins
      return -1;
    elseif nn<0
      return -1;
    else
      return nn;
  end;
end;

end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
#md ---
#md layout: post
#md title: "kinetic models part 3"
#md date: 2019-11-01 00:00:00 +0000
#md categories: blog
#md mathjax: true
#md ---

# # kinetic models part 3

# >In this post I'll introduce **an intuitive approach to kinetic models** and how they can be used for interaction problems.

# >Suppose you want to study how a virus is spreading. When two people come in contact the virus infect the non sick person. We can use a kinetic-based maths model to track: how many sick individuals we have, where are they localised and eventually the stage of the sickness (in other words: for how long they had been infected).

# >While statistics tell you that there will be a Martin Luther King every 100 years, kinetic modelling is trying to tell you who's gonna be.

# I believe this is not the place, nor yet the time, to use a full maths formalism to explain kinetic modelling. Instead, I prefer to walk you through the code so you can build
# the model in your mind via a practical example.
# To simplify the problem, avoind complex scenarios, we assume all people are in a room and their motion is limited therein. We also assume we have an initial stage where some people are sick. We can assume they contracted the virus outside and
# are bringing the infection in the room.

# Here I definie the main variable (with a *mutable struct*):
# * dt             = evolving timestep
# * ρ              = collision-distance, distance below which two individuals interact - in other words the virus is passed over
# * total_infected = a total infected counter
# * people         = the population. People is a tensor with 4 rows: (**x** position on the plane, **y** position on the plane, **flag** that identify if the person is infected on not, **stage** of infection if infected)
using Plots
mutable struct boltzmann
    dt; ρ; total_infected; people;
end

# We initialize the population
function initialize_population(n_people)
    x, y, flag, stage = rand(n_people), rand(n_people), zeros(n_people), ones(n_people)
    people = zeros(n_people,4)
    flag[1]=1
    [people[i,:]=[x[i],y[i],flag[i],stage[i]] for i=1:n_people]
    return people
end

# This is where the people move. Given a position, they are asked to step forward. Here, for simplicity, we assumed that each person moves in accordance with a normal distribution, in other words he/she can pick any angle and just make a step in that direction.
# It is rather obvious to assume that we are dealing with persons, they are probably moving straight until an obstacle stop them. However, that said it includes a concept of *objective*: I am going somewhere to do something. Since we are using, again for simplicity, a square room, such a concept is wider than the present scope.
# We include also some coundary conditions, people cannnot exit the room, they are bounced back in.
function advance(b::boltzmann)
    b.people[:,1]= b.people[:,1] + randn(size(b.people)[1]) * b.dt
    b.people[:,2]= b.people[:,2] + randn(size(b.people)[1]) * b.dt
    b.people[b.people[:,1] .< 0., 1] .*= -1
    b.people[b.people[:,2] .< 0., 2] .*= -1
    b.people[b.people[:,1] .> 1., 1]  .= 1
    b.people[b.people[:,2] .> 1., 2]  .= 1
    # b.people[b.people[:,1] .> 1, 1] .= 1-(b.people[b.people[:,1] .> 1, 1]-1)
end

# When two people are very close apart, the infection spreads. The distance at which the virus spread is ρ.
function infect(b::boltzmann)
    #loop only on infected individuals to speed up calculations
    b.people[b.people[:,3].==1,4] .+= 1
    originally_infected = b.people[b.people[:,3].==1,:]
    not_infected        = b.people[b.people[:,3].!=1,:]
    for idx = 1:size(originally_infected)[1]
        for jdx = 1:size(not_infected)[1]
            d = sqrt( (originally_infected[idx,1]-not_infected[jdx,1])^2 + (originally_infected[idx,2]-not_infected[jdx,2])^2 )
            if d < b.ρ
                not_infected[jdx,3]=1
            end
        end
    end
    b.people[1:size(originally_infected)[1],:] = originally_infected
    b.people[size(originally_infected)[1]+1:end,:] = not_infected
end

# We are also considering healing: after a person has been ill for a given *time*, he/she recovers.
function heal(b::boltzmann, time)
    #after a person has been ill for a while he/she heals
    identify = (b.people[:,4] .> time)
    b.people[identify,3] .= 0
    b.people[identify,4] .= 0
end

function count_infect(b::boltzmann)
    push!(b.total_infected, sum(b.people[:,3]) )
end

function display(b::boltzmann)
    rectangle(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

    p1 = plot(rectangle(0,0,1,1), fill="rgb(22,55,66)", xlim=(-.2,1.2), ylim=(-.2,1.2), background_color=RGB(0.2,0.2,0.2), label="room", xaxis=false, yaxis=false)
    p1 = scatter!( b.people[ findall(x->x==0., b.people[:,3]) ,1], b.people[ findall(x->x==0., b.people[:,3]) ,2], color="rgb(34,31,26)", shape=:circle, ms=4, alpha=1.0, label="ok")
    p1 = scatter!( b.people[ findall(x->x==1., b.people[:,3]) ,1], b.people[ findall(x->x==1., b.people[:,3]) ,2], color="rgb(255,114,217)", markerstrokecolor="rgb(240,229,92)", shape=:square, ms=3 .+ 0.02.*b.people[ findall(x->x==1., b.people[:,3]) ,4], alpha=1.0, label="infected")
    p1 = xlims!(0,1)
    p1 = ylims!(0,1)

    plot(p1)
end

function display2(b1::boltzmann,b2::boltzmann)
    rectangle(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

    p1 = plot(rectangle(0,0,1,1), fill="rgb(22,55,66)", xlim=(-.2,1.2), ylim=(-.2,1.2), background_color=RGB(0.2,0.2,0.2), xaxis=false, yaxis=false, dpi=500)
    p1 = scatter!( b1.people[ findall(x->x==0., b1.people[:,3]) ,1], b1.people[ findall(x->x==0., b1.people[:,3]) ,2], color="rgb(34,31,26)", shape=:circle, ms=4, alpha=1.0)
    p1 = scatter!( b1.people[ findall(x->x==1., b1.people[:,3]) ,1], b1.people[ findall(x->x==1., b1.people[:,3]) ,2], color="rgb(255,114,217)", markerstrokecolor="rgb(240,229,92)", shape=:square, ms=3 .+ 0.02.*b1.people[ findall(x->x==1., b1.people[:,3]) ,4], alpha=1.0)
    p1 = xlims!(0,1)
    p1 = ylims!(0,1)

    p2 = plot(rectangle(0,0,1,1), fill="rgb(22,55,66)", xlim=(-.2,1.2), ylim=(-.2,1.2), background_color=RGB(0.2,0.2,0.2), xaxis=false, yaxis=false, dpi=500)
    p2 = scatter!( b2.people[ findall(x->x==0., b2.people[:,3]) ,1], b2.people[ findall(x->x==0., b2.people[:,3]) ,2], color="rgb(34,31,26)", shape=:circle, ms=4, alpha=1.0)
    p2 = scatter!( b2.people[ findall(x->x==1., b2.people[:,3]) ,1], b2.people[ findall(x->x==1., b2.people[:,3]) ,2], color="rgb(255,114,217)", markerstrokecolor="rgb(240,229,92)", shape=:square, ms=3 .+ 0.02.*b2.people[ findall(x->x==1., b2.people[:,3]) ,4], alpha=1.0)
    p2 = xlims!(0,1)
    p2 = ylims!(0,1)

    p3 = plot(b1.total_infected, dpi=500)
    p3 = xlims!(0,200)
    p3 = ylims!(0,30)
    p4 = plot(b2.total_infected, dpi=500)
    p4 = xlims!(0,200)
    p4 = ylims!(0,300)

    plot(p1,p2,p3,p4, legend=false, layout = grid(2,2,heights=[0.7,0.3,0.7,0.3]) )
end

small = boltzmann((dt = 0.01, ρ = 0.05, total_infected=ones(1), people=initialize_population(30))...)
large = boltzmann((dt = 0.01, ρ = 0.05, total_infected=ones(1), people=initialize_population(300))...)

@mp4 for i=1:200
    advance(small), advance(large)
    infect(small), infect(large)
    heal(small,100), heal(large,2)
    count_infect(small), count_infect(large)
    display2(small,large)
end every 1

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

######-----------------------
using Plots
# define the Lorenz attractor
mutable struct Lorenz
    dt; σ; ρ; β; x; y; z
end

function step!(l::Lorenz)
    dx = l.σ*(l.y - l.x)       ; l.x += l.dt * dx
    dy = l.x*(l.ρ - l.z) - l.y ; l.y += l.dt * dy
    dz = l.x*l.y - l.β*l.z     ; l.z += l.dt * dz
end

attractor = Lorenz((dt = 0.02, σ = 10., ρ = 28., β = 8//3, x = 1., y = 1., z = 1.)...)

# initialize a 3D plot with 1 empty series
plt = plot3d(1, xlim=(-25,25), ylim=(-25,25), zlim=(0,50),
                title = "Lorenz Attractor", marker = 2)

# build an animated gif by pushing new points to the plot, saving every 10th frame
@gif for i=1:1500
    step!(attractor)
    push!(plt, attractor.x, attractor.y, attractor.z)
end every 10
######-----------------------

x = zeros(0)
append!( x, rand(10) )

    for idx = 1:size(people)[1]
        if people[idx,3] > 1.
            originally_infected.append(people[idx,:])
        end
    end

        people[idx,:]= [ people[idx,1]+rand()*delta_movement, people[idx,2]+rand()*delta_movement, people[idx,3] ]
    end
    return people
end

delta_movement = 1e-4
move_(people,delta_movement)

delta_movement = 1e-4
function move_(people)
    for idx = 1:size(people)[1]
        people[idx,:]= [ people[idx,1]+rand()*delta_movement, people[idx,2]+rand()*delta_movement, people[idx,3] ]
        if people[idx,1] > 1.
            people[idx,1] = 1-people[idx,1]
        end
        if people[idx,1] < 0.
            people[idx,1] = -people[idx,1]
        end
        if people[idx,2] > 1.
            people[idx,2] = 1-people[idx,2]
        end
        if people[idx,2] < 0.
            people[idx,2] = -people[idx,2]
        end
    end
    return people
end
# Despite a person trajectory depends on what the person is doing, the person is going to a colleague or to the printer, so shoud have a given trajectory depending on the schope we simplify by assuming that
# a person new position only depends on the present position and his advancment is normal distributed both in *x* and *y*. This assumption, goes, alongside with the assumption (for the moment) that points are volumeless and do not hit one another.
# These simplifications are operated to focus our attention on the *collision operator*, when two points are close enought the two people interact and if one of the persons is aware of the gossip shares it.
function collision_(people)
    new_flag = zeros( size(people)[1] )
    for idx = 1:size(people)[1]
        if people[idx,3] == 1.
            distances = ( (people[:,1].-people[idx,1]).^2 .+ (people[:,2].-people[idx,2]).^2 ).^0.5
            distances = (distances .< 0.05)
            new_flag[distances] .= 1
        end
    end
    people[:,3] = new_flag
    return people
end

# Let's adavance the system and let's check its evolution.
for time = 1:10
    move_(people)
    collision_(people)
end
plot_(people)
savefig(".\\_img\\002_kinetic_models_in_real_life\\002_kinetic_models_in_real_life_02.svg")
# ![time evolution](IMAGEFOLDER/002_kinetic_models_in_real_life_02.svg)

# # A brief disquisition on the underlying maths
# What I personally like about kinetic models is the possibility to explain them in simple words with down to heath example.
# I reckon the gossip example is rather simple, eventually funny, and intuitively simple. The underlying math follow a different formalism known as Bolzmann's statistics.
# A person, that we have represented with a circle or eventually a square, can be modelled mathematically as a Dirac Delta centered where the person is,
# ```math
# p(x_i) = \delta(x-x_i)  \:\:  :: i=1,...,N
# ```
# since each person p(x_i) is also moving is also characterised by a velocity, that using the same notation,
# ```math
# p(x_i) = \delta(x-x_i) \: \delta(v-v_i)  \:\:   :: i=1,...,N
# ```
# and to describe the entire population we sum all over the people,
# ```math
# P = \Sum_{i=0}_N \delta(x-x_i) \delta(v-v_i)
# ```
# The complexity of the problem scales linearly with the number of people *N*, and analytically we quickly encount a calculation obstaacle whether we wish to solve the system evolution by hand.
# To be able to treat analytically such a problem it is convenient to treat the ansemble behaviour of *P*, ansembled over infinitesimal volumes and infinitesimal time-intervals.

# -----------------
cd("/Users/alberto/codes/albz.github.io/_source")                                                                         #src
using Literate                                                                                #src
# preprocess for notebooks                                                                                #src
function setimagefolder(content)                                                                                #src
    content = replace(content, "IMAGEFOLDER" => "$IMAGEFOLDER")                                                           #src
    return content                                                                                #src
end                                                                                #src
# for Jupyter notebook, put images in subfolder                                                                           #src
IMAGEFOLDER = "_img/003_kinetic_models"                                                                                #src
Literate.notebook("_source/003_kinetic_models.jl", "notebooks", preprocess = setimagefolder)                              #src
# # for Markdown/Jekyll notebook, put images in "/images"                                                                 #src
# IMAGEFOLDER = "/_img/002_kinetic_models_in_real_life"                                                                   #src
Literate.markdown("_source/003_kinetic_models.jl", ".", name="_posts/0023_kinetic_models",                                #src
 preprocess = setimagefolder,                                                                                #src
 documenter=false)                                                                                #src


text:
abstract MarkovChain

abstract MarkovChainFitParams

type DTFiniteMarkovChainMleFitParams{I<:Integer,D<:Dict} <: MarkovChainFitParams
  counts_initial::AbstractVector{I}
  counts_transition::AbstractMatrix{I}
  alphabetdict::D
end

type DTFiniteMarkovChain{
      A<:Any,F<:AbstractFloat} <: MarkovChain
  initial::AbstractVector{F}
  transition::AbstractMatrix{F}
  alphabet::AbstractVector{A}
end

type DTFiniteMarkovChainWithFit{
    FMC<:DTFiniteMarkovChain,MCFP<:MarkovChainFitParams} <: MarkovChain
  mc::FMC
  fitparams::MCFP
end