text:
# White noise
"
White noise process

$FIELDS
"
struct WhiteNoise{T} <: TimeSeriesModel{T}
  "Standard deviation"
  σ::T
end

WhiteNoise() = WhiteNoise(1.0)

function Base.show(io::IO,wn::WhiteNoise)
  σstr = @sprintf "%0.3g" wn.σ
  print(io,"y = ϵ | σ = $σstr")
end

function Base.show(io::IO,::MIME"text/plain",wn::WhiteNoise)
  print(io,"Scalar white noise\n",wn)
end

order(m::WhiteNoise) = 0

Distributions.std(x::WhiteNoise) = x.σ

"
$SIGNATURES

log pdf of white noise model.
If `has_z` is `true`, the first observations is ignored, as is done for SDE(p)
processes in `logpdf` when `z` is provided. In this way, closer agreement
with the log probability of an SDE(p) process is achieved.
"
function Distributions.logpdf(proc::WhiteNoise,data,has_z=false)
  σ2 = proc.σ^2
  istart = has_z ? 2 : 1
  y = [x.y for x in data]
  Nl = length(y)
  lp = -1/2*(sum(y.^2)/σ2+Nl*log(σ2)+Nl*log(2π))
  return lp
end

function Base.rand(rng::AbstractRNG,proc::WhiteNoise,t::AbstractVector)
  data = map(x->Obs(x,proc.σ*randn(rng)),t)
  return data
end


text:
function Base.:+(A::T,B::T) where T <: CoordinateVector
    A.coord == B.coord || throw(ArgumentError("Vectors must have the same coordinates"))
    return T(A.data .+ B.data, A.coord, A.J, A.g, A.h)
end

function Base.:+(A::T,B::S) where {T<:CoordinateVector, S<:CoordinateVector}
    return A + convert(B)
end

function Base.:+(a::T,x::S) where {T<:Real, S<:CoordinateVector}
    return S(x.data+a,x.coord,x.J,x.g,x.h)
end

function Base.:+(x::T,a::S) where {T<:CoordinateVector, S<:Real}
    return T(x.data + a,x.coord,x.J,x.g,x.h)
end

function Base.:-(A::T,B::T) where {T<:CoordinateVector}
    A.coord == B.coord || throw(ArgumentError("Vectors must have the same coordinates"))
    return T(A.data .- B.data, A.coord, A.J, A.g, A.h)
end

function Base.:-(A::T,B::S) where {T<:CoordinateVector,S<:CoordinateVector}
    return A - convert(B)
end

function Base.:-(a::T,x::S) where {T<:Real,S<:CoordinateVector}
    return S(a - x.data,x.coord,x.J,x.g,x.h)
end

function Base.:-(x::T,a::S) where {T<:CoordinateVector,S<:Real}
    return T(x.data - a,x.coord,x.J,x.g,x.h)
end

function Base.:*(a::T,x::S) where {T<:Real,S<:CoordinateVector}
    return S(a*x.data,x.coord,x.J,x.g,x.h)
end

function Base.:*(x::T,a::S) where {T<:CoordinateVector,S<:Real}
    return T(a*x.data,x.coord,x.J,x.g,x.h)
end

function Base.:/(a::T,x::S) where {T<:Real,S<:CoordinateVector}
    return S(a./ x.data,x.coord,x.J,x.g,x.h)
end

function Base.:/(x::T,a::S) where {T<:CoordinateVector,S<:Real}
    return T(x.data/a,x.coord,x.J,x.g,x.h)
end

function dot(A::T,B::T) where {T<:CoordinateVector}
    A.coord === B.coord || throw(ArgumentError("Vectors must have the same coordinates"))
    dp = sum(A.data .* (A.g*B.data))
    return dp
end

function dot(A::T,B::S) where {T<:CoordinateVector,S<:CoordinateVector}
    A.coord === B.coord || throw(ArgumentError("Vectors must have the same coordinates"))
    dp = sum(A.data .* B.data)
    return dp
end

norm(x::T) where {T<:CoordinateVector} = sqrt(dot(x,x))

function cross(A::T,B::T) where {T<:CoordinateVector}
    A.coord === B.coord || throw(ArgumentError("Vectors must have the same coordinates"))
    length(A) == 3 || throw(DomainError("Cross product only defined for 3-vectors"))
    AxB = zeros(eltype(A),3)
    detJ = det(A.J)

    for i=1:3
        for j=1:3
            v = detJ*A[i]*B[j]
            for k=1:3
                AxB[k] = AxB[k] + levicivita([i,j,k])*v
            end
        end
    end

    AxB = inv(A.g)*AxB
    return T(AxB,A.coord,A.J,A.g,A.h)
end

function cross(A::T,B::S) where {T<:CoordinateVector,S<:CoordinateVector}
    return cross(A,convert(B))
end

ext: jl


text:
using Test
using ScHoLP
using Combinatorics
using LinearAlgebra
using SparseArrays

function test_projected_graph()
    data = example_dataset("example1")
    simplices, nverts = data.simplices, data.nverts
    B = basic_matrices(simplices, nverts)[3]

    @test size(B) == (9, 9)
    @test nnz(B - B') == 0
    @test sum(Diagonal(B)) == 0
    @test B[1, 2] == 2
    @test B[1, 3] == 2
    @test B[1, 4] == 1
    @test B[1, 5] == 1
    @test B[1, 6] == 2
    @test B[1, 7] == 1
    @test B[1, 8] == 1
    @test B[2, 3] == 1
    @test B[2, 4] == 1
    @test B[2, 6] == 2
    @test B[3, 4] == 1
    @test B[3, 5] == 1
    @test B[3, 9] == 1
    @test B[5, 8] == 1
    @test B[7, 8] == 1
    @test nnz(B) == 30
end

function test_basic_matrices()
    for dataset in ["example1", "example2", "email-Enron", "contact-high-school"]
        println("$dataset")
        data = example_dataset(dataset)
        simplices, nverts = data.simplices, data.nverts
        A, At, B = basic_matrices(simplices, nverts)
        @test nnz(A) == length(simplices)
        @test A' == At
        ms = maximum(simplices)
        @test size(B) == (ms, ms)
        @test size(A) == (ms, length(nverts))

        all_pairs = Set{NTuple{2, Int64}}()
        curr_ind = 1
        for nvert in nverts
            simp = simplices[curr_ind:(curr_ind + nvert - 1)]
            curr_ind += nvert
            for (i, j) in combinations(simp, 2)
                if i != j; push!(all_pairs, (i, j), (j, i)); end
            end
        end
        @test nnz(B) == length(all_pairs)
    end
end

function test_open_closed_triangles1()
    for dataset in ["example1", "example2", "email-Enron", "contact-high-school"]
        data = example_dataset(dataset)
        simplices, nverts = data.simplices, data.nverts

        # Using internal implementation
        num_open1, num_closed1 =
            num_open_closed_triangles(basic_matrices(simplices, nverts)...)

        # Naive implementation by looping over simplices
        closed_tris = Set{NTuple{3, Int64}}()
        I, J = Int64[], Int64[]
        curr_ind = 1
        for nvert in nverts
            simp = simplices[curr_ind:(curr_ind + nvert - 1)]
            curr_ind += nvert
            for (i, j, k) in combinations(simp, 3)
                push!(closed_tris, sorted_tuple(i, j, k))
            end
            for (i, j) in combinations(simp, 2)
                push!(I, i, j)
                push!(J, j, i)
            end
        end
        num_closed2 = length(closed_tris)
        @test num_closed1 == num_closed2

        n = maximum(simplices)
        A = sparse(I, J, ones(length(I)), n, n)
        total_triangles = 0
        for (i, j, k) in combinations(1:n, 3)
            if A[i, j] > 0 && A[i, k] > 0 && A[j, k] > 0
                total_triangles += 1
            end
        end
        num_open2 = total_triangles - num_closed2
        @test num_open1 == num_open2
    end
end

function test_open_closed_triangles2()
    function check(simplices::Vector{Int64}, nverts::Vector{Int64},
                   num_open::Int64, num_closed::Int64)
        num_open_check, num_closed_check =
            num_open_closed_triangles(basic_matrices(simplices, nverts)...)
        @test num_open == num_open_check
        @test num_closed == num_closed_check
    end
    simplices = Int64[]
    nverts    = Int64[]
    check(simplices, nverts, 0, 0)

    append!(simplices, [1, 2, 2, 5])
    append!(nverts,    [2,    2   ])
    check(simplices, nverts, 0, 0)

    append!(simplices, [1, 2, 2, 5])
    append!(nverts,    [2,    2   ])
    check(simplices, nverts, 0, 0)

    append!(simplices, [1, 5])
    append!(nverts,    [2])
    check(simplices, nverts, 1, 0)

    append!(simplices, [1, 2, 5])
    append!(nverts,    [3])
    check(simplices, nverts, 0, 1)

    append!(simplices, [1, 2, 5, 1, 2, 5])
    append!(nverts,    [3,       3])
    check(simplices, nverts, 0, 1)

    append!(simplices, [1, 8, 8, 2, 8, 5, 8])
    append!(nverts,    [2,    1, 2,    2   ])
    check(simplices, nverts, 3, 1)

    append!(simplices, [1, 2, 5, 8])
    append!(nverts,    [4         ])
    check(simplices, nverts, 0, 4)

    data = example_dataset("example1")
    simplices, nverts = data.simplices, data.nverts
    check(simplices, nverts, 1, 7)
end

function test_closure_counts1()
    simplices = [1, 2, 3, 1, 2, 3, 1, 2, 4, 2, 3, 4, 1, 3, 4, 8]
    nverts    = [3,       3,       3,       3,       3,       1]
    open_counts = open_types4(simplices, nverts)
    @test open_counts[(WEAK, WEAK, WEAK, STRONG)] == 1
    @test open_counts[(WEAK, WEAK, WEAK, WEAK)] == 0
    @test open_counts[(-3, -3, -3, STRONG)] == 1
    @test open_counts[(-3, -3, -3, WEAK)] == 3
    @test open_counts[(-3, -3, -3, OPEN)] == 0

    # none of the following close any open triangles
    new_simplices = [9, 11, 12, 11, 12, 14, 9, 11, 12]
    new_nverts    = [3,         3,          3]
    closed_counts = newly_closed_types4(simplices, nverts,
                                        new_simplices, new_nverts)
    for (_, val) in closed_counts; @test val == 0; end

    # Close the open tetrahedron
    new_simplices = [1, 2, 3, 4]
    new_nverts    = [4]
    closed_counts = newly_closed_types4(simplices, nverts,
                                        new_simplices, new_nverts)
    @test closed_counts[(WEAK, WEAK, WEAK, STRONG)] == 1
end

function test_equivalent_counters(counter1, counter2)
    sorted_keys(counts) = sort(collect(keys(counts)))
    skeys1 = sorted_keys(counter1)
    skeys2 = sorted_keys(counter2)
    @test skeys1 == skeys2
    for key in skeys1
        @test counter1[key] == counter2[key]
    end
end

function test_closure_counts2()
    for dataset in ["example1", "example2", "email-Enron"]
        println("$dataset...")
        data = example_dataset("example1")
        simplices, nverts = data.simplices, data.nverts
        open_counts3 = open_types3(simplices, nverts)
        open_counts4 = open_types4(simplices, nverts)

        # Close all open simplices. Add an extra simplex which should not be
        # counted since it did not appear in the original data.
        maxsize = maximum(simplices)
        new_simplices = [collect(1:maxsize); collect((maxsize + 1):(maxsize + 5))]
        new_nverts = [maxsize, 5]

        closed_counts3 = newly_closed_types3(simplices, nverts,
                                             new_simplices, new_nverts)
        test_equivalent_counters(open_counts3, closed_counts3)
        closed_counts4 = newly_closed_types4(simplices, nverts,
                                             new_simplices, new_nverts)
        test_equivalent_counters(open_counts4, closed_counts4)
    end
end

function test_closure_counts3()
    # Naive 3 nested loops computation of open type counts.
    function open_types3_naive(simplices::Vector{Int64}, nverts::Vector{Int64})
        # Get set of closed triangles
        closed_tris = Set{NTuple{3, Int64}}()
        curr_ind = 1
        for nvert in nverts
            simp = simplices[curr_ind:(curr_ind + nvert - 1)]
            curr_ind += nvert
            for (i, j, k) in combinations(simp, 3)
                push!(closed_tris, sorted_tuple(i, j, k))
            end
        end

        # naive counting
        type_counts = initialize_type_counter3()
        A, At, B = basic_matrices(simplices, nverts)
        inds = findall(vec(sum(At, dims=1)) .> 0)
        n = length(inds)
        for (i, j, k) in combinations(inds, 3)
            if !((i, j, k) in closed_tris)
                push!(type_counts, simplex_key3(B[i, j], B[i, k], B[j, k]), 1)
            end
        end
        return type_counts
    end

    # Naive implementation of open counting open 4-node simplex types.
    function open_types4_naive(simplices::Vector{Int64}, nverts::Vector{Int64})
        # Get closed triangles and tetrahedra
        closed_tris_counts = triangle_counts(simplices, nverts)
        closed_tetras = Set{NTuple{4, Int64}}()
        curr_ind = 1
        for nv in nverts
            simplex = simplices[curr_ind:(curr_ind + nv - 1)]
            curr_ind += nv
            for (i, j, k, l) in combinations(simplex, 4)
                push!(closed_tetras, sorted_tuple(i, j, k, l))
            end
        end

        # naive counting
        type_counts = initialize_type_counter4()
        A, At, B = basic_matrices(simplices, nverts)
        is_triangle(a::Int64, b::Int64, c::Int64) =
            all(v -> v > 0, [B[a, b], B[a, c], B[b, c]])
        inds = findall(vec(sum(At, dims=1)) .> 0)
        n = length(inds)
        for (i, j, k, l) in combinations(inds, 4)
            if (is_triangle(i, j, k) || is_triangle(i, j, l) ||
                is_triangle(i, k, l) || is_triangle(j, k, l)) &&
                !((i, j, k, l) in closed_tetras)
                num_edges = sum((B[i, j], B[i, k], B[i, l],
                                 B[j, k], B[j, l], B[k, l]) .> 0)
                wijk = closed_tris_counts[(i, j, k)]
                wijl = closed_tris_counts[(i, j, l)]
                wikl = closed_tris_counts[(i, k, l)]
                wjkl = closed_tris_counts[(j, k, l)]
                push!(type_counts, simplex_key4(wijk, wijl, wikl, wjkl, num_edges), 1)
            end
        end
        return type_counts
    end

    for dataset in ["example1", "example2", "email-Enron", "contact-primary-school"]
        println("$dataset...")
        data = example_dataset("example1")
        simplices, nverts = data.simplices, data.nverts
        test_equivalent_counters(open_types3_naive(simplices, nverts),
                                 open_types3(simplices, nverts))
        test_equivalent_counters(open_types4_naive(simplices, nverts),
                                 open_types4(simplices, nverts))
    end
end

function test_simplicial_closure()
    println("test_projected_graph")
    test_projected_graph()

    println("basic_matrices")
    test_basic_matrices()

    println("open_closed_triangles (1/2)")
    test_open_closed_triangles1()
    println("open_closed_triangles (2/2)")
    test_open_closed_triangles2()

    println("closure_counts (1/3)")
    test_closure_counts1()
    println("closure_counts (2/3)")
    test_closure_counts2()
    println("closure_counts (3/3)")
    test_closure_counts3()
end

function test_new_closures()
    data = example_dataset("example1")
    simplices, nverts = data.simplices, data.nverts
    old_simplices, old_nverts = Int64[], Int64[]
    new_tris = new_closures(old_simplices, old_nverts, simplices, nverts)
    @test length(new_tris) == 0
    old_simplices = collect(1:9)
    old_nverts = ones(Int64, 9)
    new_tris = new_closures(old_simplices, old_nverts, simplices, nverts)
    @test length(new_tris) == 7
end

function test_enum_open_triangles()
    data = example_dataset("example1")
    simplices, nverts = data.simplices, data.nverts
    tris = enum_open_triangles(simplices, nverts)
    @test length(tris) == 1
    @test Set{Int64}(tris[1]) == Set{Int64}([1, 5, 8])
end

function test_common_neighbor_map()
    data = example_dataset("example1")
    simplices, nverts = data.simplices, data.nverts
    tris = enum_open_triangles(simplices, nverts)
    A, At, B = basic_matrices(simplices, nverts)
    common_nbrs = common_neighbors_map(B, tris)
    @test sort(collect(common_nbrs[(1, 5)])) == [3, 8]
    @test sort(collect(common_nbrs[(1, 8)])) == [5, 7]
    @test collect(common_nbrs[(5, 8)]) == [1]
end

function test_grad_and_curl()
    data = example_dataset("example1")
    simplices, nverts = data.simplices, data.nverts
    A, At, B = basic_matrices(simplices, nverts)
    grad, curl, edge_map = grad_and_curl(A, At, B)
    @test length(edge_map) == 15
    @test nnz(grad) == 15 * 2
    @test nnz(curl) == 7 * 3

    # test gradient map for edge (1, 3)
    edge_id = edge_map[(1, 3)]
    @test grad[edge_id, 3] == 1
    @test grad[edge_id, 1] == -1
    @test nnz(grad[edge_id, :]) == 2
    degrees = vec(sum(make_sparse_ones(B), dims=1))
    @test nnz(grad[:, 1]) == degrees[1]
    @test nnz(grad[:, 3]) == degrees[3]

    # test curl map for triangle (1, 2, 4)
    edge_id1 = edge_map[(1, 2)]
    edge_id2 = edge_map[(2, 4)]
    edge_id3 = edge_map[(1, 4)]
    inds1 = findnz(curl[:, edge_id1])[1]
    inds2 = findnz(curl[:, edge_id2])[1]
    inds3 = findnz(curl[:, edge_id3])[1]
    @test length(inds1) == 3
    @test length(inds2) == 2
    @test length(inds3) == 2
    common = collect(intersect(inds1, inds2, inds3))
    @test length(common) == 1
    tri_ind = common[1]
    @test curl[tri_ind, edge_id1] == 1
    @test curl[tri_ind, edge_id2] == 1
    @test curl[tri_ind, edge_id3] == -1
    @test nnz(curl[tri_ind, :]) == 3
end

function test_score_functions()
    data = example_dataset("example1")
    simplices, nverts = data.simplices, data.nverts
    tris = enum_open_triangles(simplices, nverts)
    @test length(tris) == 1
    A, At, B = basic_matrices(simplices, nverts)
    @test arithmetic_mean(tris, B)[1] ≈ 1
    @test geometric_mean(tris, B)[1] ≈ 1

    degrees = vec(sum(make_sparse_ones(B), dims=1))
    @test pref_attach3(tris, degrees)[1] ≈ 63
    simp_degrees = vec(sum(At, dims=1))
    @test pref_attach3(tris, simp_degrees)[1] ≈ 20

    common_nbrs = common_neighbors_map(B, tris)
    @test common3(tris, common_nbrs)[1] ≈ 0
    @test jaccard3(tris, common_nbrs, degrees)[1] ≈ 0
    @test adamic_adar3(tris, common_nbrs, degrees)[1] ≈ 0

    function walk_test(func, weighted::Bool)
        scores1 = func(tris, B, weighted, true)[1]
        scores2 = func(tris, B, weighted, false)[1]
        @test scores1[1] ≈ scores2[1] atol=2e-3
    end
    walk_test(PKatz3, true)
    walk_test(PKatz3, false)
    walk_test(PPR3,   true)
    walk_test(PPR3,   false)
end

function all_tests()
    println("test_simplicial_closure")
    test_simplicial_closure()

    println("test_new_closures")
    test_new_closures()

    println("test_enum_open_triangles")
    test_enum_open_triangles()

    println("test_common_neighbor_map")
    test_common_neighbor_map()

    println("test_grad_and_curl")
    test_grad_and_curl()

    println("test_score_functions")
    test_score_functions()
end

all_tests()


text:
using Pkg

ENV["PYTHON"] = ""
ENV["CONDA_JL_HOME"] = "/home/$(ENV["USERNAME"])/miniconda3/envs/julia"

Pkg.add("PyCall")
Pkg.build("PyCall")

using PyCall

np=pyimport("numpy")

println("Numpy version: $(np.__version__)")

max_stars_repo_name: aurelio-amerio/docker-julia-pycall

max_issues_repo_name: aurelio-amerio/docker-julia-pycall

max_forks_repo_name: aurelio-amerio/docker-julia-pycall


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

@testset "LogDet" begin
    bridged_mock = MOIB.Constraint.LogDet{Float64}(mock)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 0, 1, 1, 0, 1, 0, 0, 1])
    MOIT.logdett1vtest(bridged_mock, config)
    MOIT.logdett1ftest(bridged_mock, config)
    # Dual is not yet implemented for LogDet bridge
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.LogDetConeTriangle}()))
    test_delete_bridge(bridged_mock, ci, 5, (
        (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone, 0),
        (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0),
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0)
    ))
end

@testset "RootDet" begin
    bridged_mock = MOIB.Constraint.RootDet{Float64}(mock)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 1, 1, 0, 1])
    MOIT.rootdett1vtest(bridged_mock, config)
    MOIT.rootdett1ftest(bridged_mock, config)
    # Dual is not yet implemented for RootDet bridge
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RootDetConeTriangle}()))
    test_delete_bridge(bridged_mock, ci, 4, (
        (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
        (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone, 0),
        (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0)
    ))
end


text:
module SubpixelRegistrationAF

using ComputationalResources, FFTW, SubpixelRegistration, ArrayFire

include("dftRegAF.jl")

end


text:
#
# An anonymous composite type with the given fields.
#

abstract AbstractRecord

function recordtype(fields::Vector{Symbol})
    namebuf = IOBuffer()
    print(namebuf, "Record")
    for field in fields
        print(namebuf, "#")
        print(namebuf, field)
    end
    name = Symbol(String(namebuf))
    if isdefined(RBT, name)
        return getfield(RBT, name)
    end
    typeline = []
    valline = []
    first = true
    for field in fields
        if !first
            push!(typeline, ", ")
            push!(valline, ", ")
        end
        first = false
        push!(typeline, "$field")
        push!(typeline, :( let T = fieldtype(R, $(QuoteNode(field))); !isa(T, TypeVar) ? "::$T" : ""; end ))
        push!(valline, "$field")
        push!(valline, :( let T = fieldtype(R, $(QuoteNode(field))); typeof(r.$field) != T ? "::$T" : ""; end ))
        push!(valline, " = ")
        push!(valline, :( repr(r.$field) ))
    end
    ex = Expr(
            :block,
            Expr(
                :type,
                false,
                Expr(
                    :(<:),
                    Expr(
                        :curly,
                        name,
                        (Symbol("T", k) for k in eachindex(fields))...),
                    AbstractRecord),
                Expr(
                     :block,
                     (Expr(:(::), field, Symbol("T", k)) for (k, field) in enumerate(fields))...)),
            Expr(
                :(=),
                :( Base.show{R<:$name}(io::IO, ::Type{R}) ),
                :( print(io, "@Record(", $(typeline...), ")") )),
            Expr(
                :(=),
                :( Base.show{R<:$name}(io::IO, r::R) ),
                :( print(io, "@Record(", $(valline...), ")") )))
    eval(RBT, ex)
    return getfield(RBT, name)
end

macro Record(args...)
    nargs = length(args)
    hastypes = false
    hasexprs = false
    fieldnames = Vector{Symbol}(nargs)
    fieldtypes = Vector{Any}(nargs)
    fieldexprs = Vector{Any}(nargs)
    seen = Set{Symbol}()
    for k in 1:nargs
        arg = args[k]
        if isa(arg, Expr) && arg.head == :(kw) && length(arg.args) == 2
            hasexprs = true
            fieldexprs[k] = arg.args[2]
            arg = arg.args[1]
        else
            fieldexprs[k] = nothing
        end
        if isa(arg, Expr) && arg.head == :(::) && length(arg.args) == 2
            hastypes = true
            fieldtypes[k] = arg.args[2]
            arg = arg.args[1]
        else
            fieldtypes[k] = Any
        end
        @assert isa(arg, Symbol) "ill-formed record field ($arg)"
        @assert !(arg in seen) "duplicate field name ($arg)"
        fieldnames[k] = arg
        push!(seen, arg)
    end
    ex = recordtype(fieldnames)
    if hastypes
        ex = Expr(:curly, ex, fieldtypes...)
    end
    if hasexprs
        ex = Expr(:call, ex, fieldexprs...)
    end
    return ex
end

@generated function fieldtypes{R<:AbstractRecord}(::Type{R})
    L = nfields(R)
    Ts = Vector{Type}(L)
    for k = 1:L
        T = fieldtype(R, k)
        Ts[k] = !isa(T, TypeVar) ? T : Any
    end
    return quote
        $(Expr(:meta, :inline))
        Tuple{$(Ts...)}
    end
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
import HorizonSideRobots
using HorizonSideRobots
sides = [Nord, Ost, Sud, West]
nx = [0, 1, 0, -1]
ny = [1, 0, -1, 0]
used = [Int[0 for i in 1:50] for i in 1:50]

mutable struct MyRobot
    r::Robot
    x::Int
    y::Int
    max_x::Int
    min_x::Int
    spawn::Int
    max_y::Int
    min_y::Int
    MyRobot() = new(r, 25, 25, 25, 25, 25, 25, 25)
end

function clear()
    for i in 1:50
        for j in 1:50
            used[i][j] = 0
        end
    end
end

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

function mark(R)
    if (abs(R.x + R.y - R.spawn - R.spawn) % 2 == 0)
        putmarker!(R.r)
    end
end

function up(R)
    R.min_x = min(R.x, R.min_x)
    R.min_y = min(R.y, R.min_y)
    R.max_x = max(R.x, R.max_x)
    R.max_y = max(R.y, R.max_y)
end

function dfs_find_bord(R)
    for i in 1:4
        x1 = R.x + nx[i]
        y1 = R.y + ny[i]
        if (used[x1][y1] == 0 && !isborder(R.r, sides[i]))
            used[x1][y1] = 1
            move!(R.r, sides[i])
            R.x += nx[i]
            R.y += ny[i]
            up(R)
            dfs_find_bord(R)
            R.x -= nx[i]
            R.y -= ny[i]
            move!(r, inverse(sides[i]))
        end
    end
end

function dfs(R)
    for i in 1:4
        x1 = R.x + nx[i]
        y1 = R.y + ny[i]
        if (used[x1][y1] == 0 && !isborder(r, sides[i]))
            used[x1][y1] = 1
            move!(r, sides[i])
            R.x += nx[i]
            R.y += ny[i]
            dfs(R)
            if (ismarker(R.r))
                return
            end
            move!(r, inverse(sides[i]))
            R.x -= nx[i]
            R.y -= ny[i]
        end
    end
end

function main(r)
    R = MyRobot()
    clear()
    dfs_find_bord(R)
    clear()
    dfs(R)
end

r = Robot(animate = false, "mirea-progs/Task9/temp.sit")
main(r)
show(r)

#=
ДАНО: Где-то на неограниченном со всех сторон поле и без внутренних перегородок имеется единственный маркер. Робот - в произвольной клетке поля.
РЕЗУЛЬТАТ: Робот - в клетке с тем маркером.
=#

max_stars_repo_name: Sergej-Kiryakov/mirea-progs

max_issues_repo_name: Sergej-Kiryakov/mirea-progs

max_forks_repo_name: Sergej-Kiryakov/mirea-progs


text:
immutable HTML
   s::String

  function HTML(s::String)
    return new(s)
  end
end

function HTML(tbl::Table)
  buf = IOBuffer()
  render(buf, HTMLCell, HTMLStyleElement, tbl)
  return HTML(takebuf_string(buf))
end

writemime(io::IO, ::MIME"text/html", html::HTML) = print(io, html.s)

const _htmlescape_chars = ['<'=>"<", '>'=>">",
                               '"'=>""", '\''=>"&#39",
                               '`'=>"``", # ' '=>" ",
                               '!'=>"!", '@'=>"@",
                               '$'=>"$", '%'=>"%",
                               '('=>"(", ')'=>")",
                               '='=>"=", '+'=>"+",
                               '{'=>"{", '}'=>"|",
    '['=>"[", ']'=>"["]
function htmlescape(s::String)
    s1 = replace(s, r"&(?!(\w+|\#\d+);)", "&")
    join([get(_htmlescape_chars, ch, ch) for ch in s1])
end


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

mutable struct IoK8sApiExtensionsV1beta1PodSecurityPolicy <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiExtensionsV1beta1PodSecurityPolicySpec } # spec name: spec

    function IoK8sApiExtensionsV1beta1PodSecurityPolicy(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing)
        o = new()
        validate_property(IoK8sApiExtensionsV1beta1PodSecurityPolicy, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiExtensionsV1beta1PodSecurityPolicy, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiExtensionsV1beta1PodSecurityPolicy, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiExtensionsV1beta1PodSecurityPolicy, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        o
    end
end # type IoK8sApiExtensionsV1beta1PodSecurityPolicy

const _property_map_IoK8sApiExtensionsV1beta1PodSecurityPolicy = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"))
const _property_types_IoK8sApiExtensionsV1beta1PodSecurityPolicy = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiExtensionsV1beta1PodSecurityPolicySpec")
Base.propertynames(::Type{ IoK8sApiExtensionsV1beta1PodSecurityPolicy }) = collect(keys(_property_map_IoK8sApiExtensionsV1beta1PodSecurityPolicy))
Swagger.property_type(::Type{ IoK8sApiExtensionsV1beta1PodSecurityPolicy }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiExtensionsV1beta1PodSecurityPolicy[name]))}
Swagger.field_name(::Type{ IoK8sApiExtensionsV1beta1PodSecurityPolicy }, property_name::Symbol) =  _property_map_IoK8sApiExtensionsV1beta1PodSecurityPolicy[property_name]

function check_required(o::IoK8sApiExtensionsV1beta1PodSecurityPolicy)
    true
end

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


text:
module Inno
export iter, Atan, OU, generate, runexamples
using Winston
#using Debug
#using Distributions
include("SDE.jl")
using .SDE
include("misc.jl")
import .SDE.sigma, .SDE.b, .SDE.a #, SDE.UvPro, SDE.CTPro
import .SDE.dotVs, .SDE.Vs

####### dimension

const d = 1

######## data?
generated = false

####### Atan

type Atan <: UvPro
    thC
    thnonC
    d::Int
    Atan(B, beta, sigma) = new([B, beta], [sigma], 1)
end
nphi(_::Atan) = 2
nth(_::Atan) = 3

b(t,x, P::Atan) =  P.thC[1]*phi1(t,x, P) + P.thC[2]*phi2(t, x, P)
sigma(t,x, P::Atan) = P.thnonC[1]

Bof(P::Atan) = P.thC[1]*(cos(-P.thC[2]/P.thC[1]))^2
betaof(P::Atan) = 1/(2P.thC[1])*sin(2*P.thC[2]/P.thC[1])

#Bof(P::Atan) = P.thC[1]
#betaof(P::Atan) = P.thC[2]

muof(P::Atan) = 0.

phi1(t, x, _::Atan) = atan(x)
phi2(t, x, _::Atan) = 1.
phi(t, x, _::Atan) = [atan(x), 1.]

####### Prior
import Base.clamp;
clamp(x::Vector,l::Vector, u::Vector) = map(i -> clamp(x[i], l[i], u[i]), 1:length(x))

function generate(Model_; thtrue = [-2., 0.0, 0.75], u = 0.1, Nobs = 101, N = 400000 + 1, Delta=0.3, SV=true)
    global Model = Model_
    global Ptrue = Model(thtrue...)
    global model = string(Model)

    global M = Nobs - 1 # number of bridges
    # N - full observations i*M*n+1 for natural number i
    global Ttotal = Delta*M      #total time span

    tag = "$model$(M)T$Ttotal"
    println(tag)

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

    ############
    print("Generate X")
    ############

    srand(3)
    tt = collect(linspace(0., Ttotal, N))
    W = sample(tt, Wiener())
    xx = euler(u, W, Ptrue).yy
    println(".")

    global ttf = tt[1:(N-1)//(M*100):end]
    global xxf = xx[1:(N-1)//(M*100):end]
    global ttd = tt[1:(N-1)//M:end]
    global xxd = xx[1:(N-1)//M:end]

    global quxx = (sum(diff(xx[:]).^2)/Ttotal)
    global quxxd = (sum(diff(xxd[:]).^2)/Ttotal)
    println("si2 $( (thtrue[end])^2) qufull $quxx quobs $quxxd")
    global generated = true
    println("beta ", (xxd[end] - xxd[1])/Ttotal)

    if SV
        writecsv("xtrue$(tag).csv",[linspace(0, Ttotal, length(xx)) xx])
        writecsv("xobs$(tag).csv",[linspace(0, Ttotal, length(xxd)) xxd])
    end
end

#randq(thold, dl) = clamp( thold .* [ones(Nphi), exp(dl*(2rand(Nth - Nphi)-1))], [-Inf, -Inf, 0.0], [-0.001, Inf, Inf])
function randq(thold, delta)
        th = copy(thold)
        th[end] *= exp(delta*(2rand()-1))
#        th[1] = min(th[1], -0.001)
#        th[3] += delta*(2rand()-1)
        return th
end

function iter(K, n; TC = false, SV = true, DRAW = false, VERBOSE = 1, proptype = :lin, th1 = [-0.1, -0.1, 2.0], delta=0.1)
# n - timepoints each bridge (excluding excluding endpoint)
# K - number of iterations
    if !generated error("Use Inno.generate(Inno.Atan) or Inno.generate(Inno.OU) to generate observation process.") end

    tag = "$model$K$proptype$(M)x$(n)T$Ttotal"*(TC ? "TC":"")

    TC2 = TC

    Nth = nth(Ptrue)
    Nphi = nphi(Ptrue)

## PRIOR
    #uniform on 1/sigma
    pth(th) = 1. #pdf(InverseGamma(0.001, 0.001), (th[end])^2)
    qth(th, thold) = 1.

    xi = 0.2*ones(Nphi) # si^2 = 5

    ############
    VERBOSE > 0 && println("Generate bridges (using proposals $proptype): $tag.")
    ############
    srand(3)
    global thetas = zeros(Nth, K)
    global bb = zeros(Bool, M, K)
    global bbp = zeros(M, K)
    global aacc = zeros(Bool, K)
    global aaccp = zeros(K)

    llold = zeros(M)
    ll = zeros(M)
    ll2 = 0.

    global llmax = zeros(M)
    global ll2max = 0

    TT = cell(M)
    TTU = cell(M)
    W = Array(CTPath{1},M)
    W° = Array(CTPath{1},M)
    Y = Array(CTPath{1},M)
    Y° = Array(CTPath{1},M)
    Z = Array(CTPath{1},M)
    Z° = Array(CTPath{1},M)
    U = Array(CTPath{1},M)
    U° = Array(CTPath{1},M)
    Bth = zeros(Nth, M)
    Bth° = zeros(Nth, M)

    # allocate arrays for bridges
    dt = Ttotal/M/n
    for m = 1:M
        Tmin = ttd[m]
        Tmax = ttd[m+1]
        TTU[m] = collect(linspace(0.0, Tmax-Tmin, n))
        if TC
                TT[m] = map(s -> tofs(s, Tmin, Tmax-Tmin), TTU[m])
        else
                TT[m] = collect(linspace(Tmin, Tmax, n))
        end
        W[m] = UvPath(TT[m])
        W°[m] = UvPath(TT[m])
        Y[m] = UvPath(TT[m])
        Y°[m] = UvPath(TT[m])
        U[m] = UvPath(TTU[m])
        U°[m] = UvPath(TTU[m])
        Z[m] = UvPath(TTU[m])
        Z°[m] = UvPath(TTU[m])

    end

    th = copy(th1)

    srand(3)
    bbsum = 0
    aaccsum = 0
    m1th = 0*th
    m2th = eps2 .+ 0*th

    # conjugate posterior distribution N(mu, Si)
    mu = zeros(Nphi)
    Si = zeros(Nphi, Nphi)

    for k = 1:K

        m1th[:] = m1th .+ th
        m2th[:] = m2th .+ th.^2
        VERBOSE > 0 && k > 1 && print("k $k $(mc2str(k, m1th, m2th)) $(repr(round(th,3)))")

        #### generate bridges
        P = Model(th...)

        for m = 1:M

            u = xxd[m]
            v = xxd[m+1]

            Tmin = ttd[m]
            Tmax = ttd[m+1]
            T = Tmax-Tmin

            if proptype == :brown || proptype == :brownbrown
                Pt = UvAffPro(muof(P), sigma(Tmax,v, P))
                #Pt  = UvLinPro(-1E-7, 0.0, sigma(Tmax,v, P))
            elseif proptype == :lin || proptype == :linlin
                Pt = UvLinPro(Bof(P),betaof(P), sigma(Tmax,v,P))
            end

            Bth°[:,m] = th
            if TC
                resample!(Z°[m], Wiener())
#                println(Z°[1])
                u0 = uofx(0., u, T, v, Pt)
                eulerU!(U°[m], u0, Z°[m], Tmin, T, v, Pt, P)
                XofU!(Y°[m], U°[m], Tmin, T, v, Pt)
                ll[m] = llikeliU(U°[m], Tmin, T, v, Pt, P)
                llold[m] = llikeliU(U[m], Tmin, T, v, Pt, P)
            else
                resample!(W°[m], Wiener())
                guidedeuler!(Y°[m], u, W°[m], Tmax, v, Pt, P)
                ll[m] = llikeliXcirc(Y°[m], Pt, P)
                llold[m] = llikeliXcirc(Y[m], Pt, P)
            end

            llmax = map(max, ll, llmax)
            if(k == 1)
                    Bth[m] = Bth°[m]
                    Z[m].yy[:] = Z°[m].yy
                    U[m].yy[:] = U°[m].yy
                    Y[m].yy[:] = Y°[m].yy
                    Y[m].tt[:] = Y°[m].tt
                    W[m].yy[:] = W°[m].yy
                    llold[m] = ll[m]
            end

            acc = min(1.0, exp(ll[m]-llold[m]))
            VERBOSE > 1 && println("\t acc ", round(acc, 3), " ", round(llold[m],4), " ", round(ll[m],4))
            bbp[m, k] = acc
            bbsum += acc
            bb[m, k] = false
            if rand() <= acc
                bb[m, k] = true
                Bth[m], Bth°[m] = Bth°[m], Bth[m]
	            Y[m], Y°[m] = Y°[m], Y[m]
                W[m], W°[m] = W°[m], W[m]
	            Z[m], Z°[m] = Z°[m], Z[m]
                U[m], U°[m] = U°[m], U[m]
            end
        end #for m

        if DRAW
             I = 5
             p = FramedPlot()
             setattr(p, "xrange", (Y[1].tt[1], Y[I].tt[end]))
             setattr(p, "yrange", (-5.,5.))
             add(p, Points(ttd[1:I+1],xxd[1:I+1], "color","red"))
             add(p, Curve(ttf[1:100*I+1],xxf[1:100I+1], "color","red","linewidth", 0.5))
        end

        if DRAW
           for i in 1:I
                    if bb[m,k]
                        add(p, Curve(Y[i].tt[:],Y[i].yy[:], "color","black", "linewidth", 0.5))
                        add(p, Curve(Y[i].tt[:],Y°[i].yy[:], "color","grey", "linewidth", 0.5))
                    else
                        add(p, Curve(Y[i].tt[:],Y[i].yy[:], "color","grey", "linewidth", 0.5))
                        add(p, Curve(Y[i].tt[:],Y°[i].yy[:], "color","orange", "linewidth", 0.5))
                    end
             end
             Winston.display(p)
        end

        P = Model(th...)
        th° = randq(th, delta)
		P° = Model(th°...)

		VERBOSE > 0 && k > 1 && print(" ", round(100*bbsum/k/M))

        #### update sigma and theta
        laccth = 0.

        for m in 1:M
            u = xxd[m]
            v = xxd[m+1]

            Tmin = ttd[m]
            Tmax = ttd[m+1]
            T = Tmax-Tmin

            if proptype == :brown || proptype == :brownbrown
                Pt = UvAffPro(muof(P), sigma(Tmax,v, P))
                Pt° = UvAffPro(muof(P°),  sigma(Tmax,v, P°))
                #Pt = UvLinPro(-1E-6, 0.0, sigma(Tmax,v, P))
                #Pt° = UvLinPro(-1E-6, 0.0,  sigma(Tmax,v, P°))
            elseif proptype == :lin || proptype == :linlin
                Pt = UvLinPro(Bof(P),betaof(P), sigma(Tmax,v,P))
                Pt° = UvLinPro(Bof(P°),betaof(P°), sigma(Tmax,v, P°))
            end

            # compute acceptance probability of change in theta with equivalent innovations

            ### bridge Y° (U) using original W (Z)
            Bth°[:,m] = th°
            if TC2
                u0 = uofx(0., u, T, v, Pt)
                eulerU!(U[m], u0, Z[m], Tmin, T, v, Pt, P)
                XofU!(Y[m], U[m], Tmin, T, v, Pt)

                u0 = uofx(0., u, T, v, Pt°)
                eulerU!(U°[m], u0, Z[m], Tmin, T, v, Pt°, P°)
                XofU!(Y°[m], U°[m], Tmin, T, v, Pt°)

                m3 = +llikeliU(U°[m], Tmin, T, v, Pt°, P°)
                m4 = -llikeliU(U[m], Tmin, T, v, Pt, P) # assert that U is a Pt-P-proposal
#                m3 = +llikeliXcirc(Y°[m], Pt°, P°) # for testing: this should yield the same indirectly
#                m4 = -llikeliXcirc(Y[m], Pt, P)
            else
                guidedeuler!(Y[m], u, W[m], Tmax, v, Pt, P)
                guidedeuler!(Y°[m], u, W[m], Tmax, v, Pt°, P°)

                m3 = +llikeliXcirc(Y°[m], Pt°, P°)
                m4 = -llikeliXcirc(Y[m], Pt, P)
            end
            m1 = +lp(Tmin, u, Tmax, v, Pt°)
            m2 = -lp(Tmin, u, Tmax, v, Pt)
            laccth += m1 + m2 + m3 + m4
        end #for m
        ll2max = max(laccth, ll2max)
        accth = min(1., exp(laccth) * pth(th°)/pth(th)*qth(th, th°)/qth(th°, th))
        aacc[k] = false
        aaccp[k] = accth
        VERBOSE > 0 && print(" th° $(repr(round(th°,3))) Pacc $(repr(round(accth,3))) avg ", round(100aaccsum/k,2))

        if rand() <= accth
            aacc[k] = true
            aaccsum += 1

            th[:] = th°
            for m = 1:M
                # W and Z do not change
                Bth[m] = Bth°[m]
                Y[m], Y°[m] = Y°[m], Y[m]
                U[m], U°[m] = U°[m], U[m]
            end

        end
        if false
             for i in 1:I
                    if aacc[k]
                        add(p, Curve(Y[i].tt[:],Y[i].yy[:], "color","black", "linewidth", 0.5))
                        add(p, Curve(Y[i].tt[:],Y°[i].yy[:], "color","grey", "linewidth", 0.5))
                    else
                        add(p, Curve(Y[i].tt[:],Y[i].yy[:], "color","grey", "linewidth", 0.5))
                        add(p, Curve(Y[i].tt[:],Y°[i].yy[:], "color","orange", "linewidth", 0.5))
                    end
             end
             Winston.display(p)
        end

        VERBOSE > 0 && print(aacc[k])

        #### update conjugate parameters
        P = Model(th...)
        mu[:] = 0.
        Si[:] = 0.

        for m in 1:M
            t = Y[m].tt
            yy = Y[m].yy
            if Nphi == 2
                for i in 1:length(t)-1
                    ainv = 1./sigma(t[i], yy[i], P)^2
                    gdy = (yy[i+1]-yy[i])*ainv
                    gdt = (t[i+1]-t[i])*ainv
                    p1 = phi1(t[i], yy[i], P)
                    p2 = phi2(t[i], yy[i], P)
                    mu[1] += p1*gdy
                    mu[2] += p2*gdy
                    Si[1,1] +=  p1*p1 * gdt
                    Si[1,2] +=  p1*p2 * gdt
                    Si[2,2] +=  p2*p2 * gdt
                end
                Si[2,1] = Si[1,2]
            else
                #stop()
                for i in 1:length(t)-1
                    ainv = 1/sigma(t[i], yy[i], P)^2
                    phii = phi(t[i], yy[i], P)
                    dy = yy[i+1]-yy[i]
                    mu[:] += phii*dy*ainv
                    Si[:] = Si +  phii*phii' * (t[i+1]-t[i])*ainv
                end
            end
        end #for m

        # sampling conditional posterior
        WW = Si + diagm(xi)
        WL = chol(WW)
        th[1:Nphi] = WL\(WL'\mu + randn(Nphi))
        thetas[:, k] = th
        # th changed, and so the innovations process

        P = Model(th...)
        for m in 1:M
            u = xxd[m]
            v = xxd[m+1]

            Tmin = ttd[m]
            Tmax = ttd[m+1]
            T = Tmax-Tmin

            if proptype == :brown || proptype == :brownbrown
                Pt = UvAffPro(muof(P), sigma(Tmax,v, P))
                #Pt  = UvLinPro(-1E-7, 0.0, sigma(Tmax,v, P))
            elseif proptype == :lin || proptype == :linlin
                Pt = UvLinPro(Bof(P), betaof(P), sigma(Tmax,v,P))
            end

            Bth°[:,m] = th
            if TC
                Z°[m].yy[1] = 0.

                for i in 1:length(U[m].tt)-1
                    s = U[m].tt[i]
                    u1 = U[m].yy[i]
                    t, y = txofsu(s, u1, Tmin, T, v, Pt)
                    ds = U[m].tt[i+1] - s
                    T1 = sqrt(T*(T-s)/2)
                    T2 = T1*2/T
                    T3 = sqrt((T/2)/(T-s))
                    Z°[m].yy[i+1] = Z°[m].yy[i] -inv(sigma(t, y, P))*(T1*(U[m].yy[i+1] - u1) -
                                T2*(dotVs(s, T, v, Pt) - b(t,y, P))*ds +
                                T3*(u1 - 2*a(t, y, P)*SDE.Ju(s, T, Pt, u1))*ds)
                end

            else
                W°[m].yy[1] = 0.
                for i in 1:length(Y[m].tt)-1
                    t = Y[m].tt[i]
                    dt = Y[m].tt[i+1] - t
                    y = Y[m].yy[i]
                    W°[m].yy[i+1] = W°[m].yy[i] + inv(sigma(t, y, P))*((Y[m].yy[i+1] - y) - b(t,y, P)*dt - a(t, y, P)*r(t, y, Tmax, v, Pt)*dt)
                end
            end
        end
        if false #DRAW && TC
            for i in 1:I
                add(p, Curve(Y[i].tt[:],Z[i].yy[:], "color","green", "linewidth", 0.5))
                add(p, Curve(Y[i].tt[:],Z°[i].yy[:], "color","blue", "linewidth", 0.5))
            end
            Winston.display(p)
        end
        for m in 1:M
            Bth[m] = Bth°[m]
            # Y and U do not change
            W[m], W°[m] = W°[m], W[m]
            Z[m], Z°[m] = Z°[m], Z[m]
        end
        println()

    end # for K

    if SV
        if !isdir(tag)
             mkdir(tag)
        end
        writecsv("$tag/thetas$tag.csv", thetas')
        accs = [aacc aaccp bb' bbp']
        writecsv("$tag/accs$tag.csv", accs)
        open("$tag/acc$tag.txt", "w") do f
            println(f,"qufull $quxx quobs $quxxd")
            println(f,"estimate: $(repr(mean(thetas,2)[:]))")
            println(f,"est: $(mc2str(K, m1th, m2th))")
            println(f,"bridge acc ",  repr(round(100*mean(1.0bb),2)), "%, avg prob  ", repr(round(100*mean(bbp),2)) )
            println(f, "llmax ", maximum(llmax))
            #println(f,"bridge acc ",  repr(round(100*mean(1.0bb,2),2)), "%, avg prob  ", repr(round(mean(bbp, 2),2)) )
            println(f,"sigma acc ",  round(100*sum(aacc)/K,2), "% avg prob ", round(100*sum(aaccp)/K, 2), "% " )
        end
    end

    println("estimate: $(repr(mean(thetas,2)[:]))")
    println("bridge acc ",  round(100*mean(1.0bb,2),2), "%, avg prob  ", repr(round(100*mean(bbp, 2),2)) )
    println("sigma acc ",  round(100*sum(aacc)/K,2), "% avg prob ", round(100*sum(aaccp)/K, 2), "% " )
    return thetas, accs
end

function runexamples()
    generate(Atan)
    for n in [10, 100,1000], t in [:lin], tc in [true, false]
        println("$n $t $tc")
        iter(10000,n; proptype=t, DRAW=true, TC=tc, th1=[-.1, -.1, 2.])
    end
end

end #module


text:
@testset "Persistence" begin
gpath = joinpath(testdir, "testdata")
pdict = loadgraphs(joinpath(gpath, "swgs.jgz"), SWGFormat())
p1 = pdict["g"]
p2 = pdict["h"]

(f, fio) = mktemp()
# test :lg
@test savegraph(f, p1, SWGFormat()) == 1
@test savegraph(f, p1, SWGFormat(); compress=true) == 1
@test savegraph(f, p1, SWGFormat(); compress=true) == 1
@test savegraph(f, p2, SWGFormat(); compress=true) == 1
@test (ne(p2), nv(p2)) == (3, 3)

g2 = loadgraph(f, SWGFormat())
j2 = loadgraph(f, "graph", SWGFormat())
@test g2 == j2
@test (ne(g2), nv(g2)) == (3, 3)

(f, fio) = mktemp()
@test length(sprint(savegraph, p1, SWGFormat())) == 104
@test length(sprint(savegraph, p2, SWGFormat())) == 104
gs = loadgraph(joinpath(gpath, "swgs.jgz"), "h", SWGFormat())
@test gs == p2
@test_throws ErrorException loadgraph(joinpath(gpath, "swgs.jgz"), "badname", SWGFormat())

@test savegraph(f, p1) == 1
d = Dict{String,AbstractGraph}("p1" => p1, "p2" => p2)
@test savegraph(f, d, SWGFormat()) == 2
@test savegraph(f, d) == 2

close(fio)
rm(f)
end


text:
macro myevalpoly(z,a...)
    isempty(a) && error("You forgot to pass coefficients!")
    ex = :($(a[length(a)]))
    for i in 1:length(a)-1
       ex = :($ex * $(z) + $(a[length(a)-i]) )
    end
    println(ex)
    ex
end

@myevalpoly 7 2 3 4 5
@evalpoly 7 2 3 4 5

##
function root2coeff(z::AbstractVector{T}) where T
    N = length(z)
    co = zeros(T, N+1)
    # The last coefficient is always one
    co[end] = 1
    # The outer loop adds one root at a time
    for j in 1:N, i in j:-1:1
        co[end-i] -= z[j]*co[end-i+1]
    end
    co
end
@show typemax(Int), typemax(Int128)
root2coeff(1:20)
root2coeff(Int128(1):20)

##
using LinearAlgebra
function poly_roots(z)
    len = length(z)
    # construct the ones part
    mat = diagm(-1 => ones(len-2))
    # insert coefficients
    mat[:, end] = -z[1:end-1]
    eigvals(mat)
end

## Calculate all the roots and plot it
using Random
Random.seed!(2020)
function wilkinson_poly_roots(n=100)
    # original coefficients
    coeff = root2coeff(Int128(1):20)
    rts = Vector{Complex{Float64}}[]
    # add perturbation
    for i in 1:n
        pert_coeff = coeff.*(1 .+ rand(21)*1e-6)
        push!(rts, poly_roots(pert_coeff))
    end
    rts
end
using Plots; gr()
function plt_wilkinson_roots(rts)
    # plot roots without perturbation
    plt = scatter(1:20, zeros(20), color = :green, markersize = 5, legend = false, grid = false, aspect_ratio = 1)
    for i in eachindex(rts)
        # plot roots with perturbation
        scatter!(plt, real.(rts[i]), imag.(rts[i]), color = :red, mscolor = :red, markersize = 2, shape = :star)
    end
    plt
end
wilkinson_poly_roots(100) |> plt_wilkinson_roots


text:
using OpenEphysLoader
using Main.TestUtilities, Main.TestContinuous

using Test

@testset "Continuous" begin
    ### Tests ###
    samps_per_block = OpenEphysLoader.CONT_REC_N_SAMP
    cont_rec_block_size = OpenEphysLoader.CONT_REC_BLOCK_SIZE
    # sampno_to_block
    @test OpenEphysLoader.sampno_to_block(1) == 1
    @test OpenEphysLoader.sampno_to_block(samps_per_block) == 1
    @test OpenEphysLoader.sampno_to_block(samps_per_block + 1) == 2
    @test OpenEphysLoader.sampno_to_block(0) == 0

    # sampno_to_offset
    @test OpenEphysLoader.sampno_to_offset(1) == 1
    @test OpenEphysLoader.sampno_to_offset(2) == 2
    @test OpenEphysLoader.sampno_to_offset(samps_per_block) ==
        samps_per_block
    @test OpenEphysLoader.sampno_to_offset(samps_per_block + 1) == 1

    # block_start_pos
    @test OpenEphysLoader.block_start_pos(1) == 1024
    @test OpenEphysLoader.block_start_pos(2) == 3094

    # pos_to_blockno
    @test_throws ArgumentError OpenEphysLoader.pos_to_blockno(-1)
    @test OpenEphysLoader.pos_to_blockno(0) == 0
    @test OpenEphysLoader.pos_to_blockno(1023) == 0
    @test OpenEphysLoader.pos_to_blockno(1024) == 1
    @test OpenEphysLoader.pos_to_blockno(1025) == 1
    @test OpenEphysLoader.pos_to_blockno(1024 + cont_rec_block_size) == 2

    # blockno_to_start_sampno
    @test OpenEphysLoader.blockno_to_start_sampno(1) == 1
    @test OpenEphysLoader.blockno_to_start_sampno(2) == 1025

    # verify_tail!
    badtail = b"razzmatazz"
    tailarr = Vector{UInt8}(undef, 10)
    filecontext(writeio -> write(writeio, badtail)) do io
        @test OpenEphysLoader.verify_tail!(io, tailarr) == false
    end
    write_block_end = io -> write(io,
                                  OpenEphysLoader.CONT_REC_END_MARKER)
    filecontext(write_block_end) do io
        @test OpenEphysLoader.verify_tail!(io, tailarr) == true
    end

    testblockhead = OpenEphysLoader.BlockHeader()
    testblock = OpenEphysLoader.DataBlock()
    # Data conversion
    blockdata = rand_block_data()
    copyto!(testblock.body, to_OE_bytes(blockdata))
    OpenEphysLoader.convert_block!(testblock)
    @test testblock.data == blockdata

    # Test continuous interfaces
    # Test sample iterator
    nblock = 2
    testdata = vcat([rand_block_data() for i = 1:nblock]...)
    recno = OpenEphysLoader.CONT_REC_REC_NO_BITTYPE(0)
    startsamp = OpenEphysLoader.CONT_REC_TIME_BITTYPE(1)
    filecontext(write_continuous, testdata, recno, startsamp) do io
        # counting utilities
        @test OpenEphysLoader.count_data(io) == length(testdata)

        @test OpenEphysLoader.count_blocks(io) == nblock

        @test OpenEphysLoader.check_filesize(io)

        # File seeking
        OpenEphysLoader.seek_to_block(io, 2)
        @test position(io) == 3094
        OpenEphysLoader.seek_to_block(io, 1)
        @test position(io) == 1024

        # Block header reading
        @test OpenEphysLoader.read_into!(io, testblockhead)
        verify_BlockBuffer(testblockhead, startsamp, recno)

        # Block buffering
        OpenEphysLoader.seek_to_block(io, 1)
        @test OpenEphysLoader.read_into!(io, testblock, true)
        block_data, block_body, block_startsamp =
            to_block_contents(testdata, 1, startsamp)
        verify_BlockBuffer(
            testblock,
            block_startsamp,
            recno,
            block_body,
            block_data
        )

        # Test array types
        sampletypes = [Int,
                       OpenEphysLoader.CONT_REC_SAMP_BITTYPE,
                       Float64]
        test_OEContArray(io,
                         SampleArray,
                         sampletypes,
                         testdata,
                         nblock,
                         recno,
                         startsamp)
        test_OEContArray(io,
                         TimeArray,
                         sampletypes,
                         testdata,
                         nblock,
                         recno,
                         startsamp)
        recnotypes = [Int]
        test_OEContArray(io,
                         RecNoArray,
                         recnotypes,
                         testdata,
                         nblock,
                         recno,
                         startsamp)
        jointtypes = [Tuple{Float64, Float64, Int}, Tuple{Int, Int, Int}]
        test_OEContArray(io,
                         JointArray,
                         jointtypes,
                         testdata,
                         nblock,
                         recno,
                         startsamp)
    end

    # Totally destroyed file
    OEContArrayTypes = [SampleArray, TimeArray, RecNoArray, JointArray]
    filecontext(bad_file) do io
        # check_filesize
        @test ! OpenEphysLoader.check_filesize(io)

        # read_into!
        @test ! OpenEphysLoader.read_into!(io, testblockhead)
        seekstart(io)
        @test ! OpenEphysLoader.read_into!(io, testblock, true)

        # ContinuousFile constructor
        seekstart(io)
        @test_throws CorruptedException ContinuousFile(io)

        # Array constructors
        for t in OEContArrayTypes
            seekstart(io)
            @test_throws CorruptedException t(io)
        end
    end

    # Test reading from blocks
    filecontext(bad_blockhead) do io
        @test ! OpenEphysLoader.read_into!(io, testblockhead)
        seekstart(io)
        @test ! OpenEphysLoader.read_into!(io, testblock, true)
    end

    filecontext(bad_blocktail) do io
        @test OpenEphysLoader.read_into!(io, testblockhead)
        seekstart(io)
        @test ! OpenEphysLoader.read_into!(io, testblock, true)
        seekstart(io)
        @test OpenEphysLoader.read_into!(io, testblock, false)
    end

    filecontext(damaged_file, testdata, recno, startsamp) do io
        for t in OEContArrayTypes
            seekstart(io)
            A = t(io)
            @test_throws(CorruptedException,
                         OpenEphysLoader.prepare_block!(A, 3073))
            @test_throws CorruptedException A[end]
        end
    end
end


text:
# ----------------------------------------- #
# SyncPlot -- sync Plot object with display #
# ----------------------------------------- #
mutable struct SyncPlot
    plot::PlotlyBase.Plot
    scope::Scope
    window::Union{Nothing,Blink.Window}
end

Base.getindex(p::SyncPlot, key) = p.scope[key] # look up Observables

@WebIO.register_renderable(SyncPlot) do plot
    return WebIO.render(plot.scope)
end

function Base.show(io::IO, mm::MIME"text/html", p::SyncPlot)
    # if we are rendering docs -- short circuit and display html
    if get_renderer() == DOCS
        return show(io, mm, p.plot, full_html=false, include_plotlyjs="require-loaded", include_mathjax=missing)
    end
    show(io, mm, p.scope)
end
Base.show(io::IO, mm::MIME"application/prs.juno.plotpane+html", p::SyncPlot) = show(io, mm, p.scope)

function SyncPlot(
        p::Plot;
        kwargs...
    )
    lowered = JSON.lower(p)
    id = string("#plot-", p.divid)

    # setup scope
    deps = [
        "Plotly" => joinpath(artifact"plotly-artifacts", "plotly.min.js"),
        joinpath(@__DIR__, "..", "assets", "plotly_webio.bundle.js")
    ]
    scope = Scope(imports=deps)
    scope.dom = dom"div"(id=string("plot-", p.divid))

    # INPUT: Observables for plot events
    scope["hover"] = Observable(Dict())
    scope["selected"] = Observable(Dict())
    scope["click"] = Observable(Dict())
    scope["relayout"] = Observable(Dict())
    scope["image"] = Observable("")
    scope["__gd_contents"] = Observable{Any}(Dict())  # for testing

    # OUTPUT: setup an observable which sends modify commands
    scope["_commands"] = Observable{Any}([])
    scope["_toImage"] = Observable(Dict())
    scope["_downloadImage"] = Observable(Dict())
    scope["__get_gd_contents"] = Observable("")

    onjs(scope["_toImage"], @js function (options)
        this.Plotly.toImage(this.plotElem, options).then(function (data)
            $(scope["image"])[] = data
        end)
    end)

    onjs(scope["__get_gd_contents"], @js function (prop)
        if prop == "data"
            $(scope["__gd_contents"])[] = this.plotElem.data
        end

        if prop == "layout"
            $(scope["__gd_contents"])[] = this.plotElem.layout
        end
    end)

    onjs(scope["_downloadImage"], @js function (options)
        this.Plotly.downloadImage(this.plotElem, options)
    end)

    # Do the respective action when _commands is triggered
    onjs(scope["_commands"], @js function (args)
        @var fn = args.shift()
        @var elem = this.plotElem
        @var Plotly = this.Plotly
        args.unshift(elem) # use div as first argument
        Plotly[fn].apply(this, args)
    end)

    onimport(scope, JSExpr.@js function (Plotly, PlotlyWebIO)
        # Add the PlotlyCommands object to the WebIO JS instance
        PlotlyWebIO.init(WebIO)

        # set up container element
        @var gd = this.dom.querySelector($id);

        # save some vars for later
        this.plotElem = gd
        this.Plotly = Plotly
        if (window.Blink !== undefined)
            # set css style for auto-resize
            gd.style.width = "100%";
            gd.style.height = "100vh";
            gd.style.marginLeft = "0%";
            gd.style.marginTop = "0vh";
        end

        window.onresize = () -> Plotly.Plots.resize(gd)

        # Draw plot in container
        Plotly.newPlot(
            gd,
            $(lowered[:data]),
            $(lowered[:layout]),
            $(lowered[:config]),
        )

        # hook into plotly events
        gd.on("plotly_hover", function (data)
            @var filtered_data = WebIO.PlotlyCommands.filterEventData(gd, data, "hover");
            if !(filtered_data.isnil)
                $(scope["hover"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_unhover", () -> $(scope["hover"])[] = Dict())

        gd.on("plotly_selected", function (data)
            @var filtered_data = WebIO.PlotlyCommands.filterEventData(gd, data, "selected");
            if !(filtered_data.isnil)
                $(scope["selected"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_deselect", () -> $(scope["selected"])[] = Dict())

        gd.on("plotly_relayout", function (data)
            @var filtered_data = WebIO.PlotlyCommands.filterEventData(gd, data, "relayout");
            if !(filtered_data.isnil)
                $(scope["relayout"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_click", function (data)
            @var filtered_data = WebIO.PlotlyCommands.filterEventData(gd, data, "click");
            if !(filtered_data.isnil)
                $(scope["click"])[] = filtered_data.out
            end
        end)
    end)

    # create no-op `on` callback for image so it is _always_ sent
    # to us
    on(scope["image"]) do x end

    SyncPlot(p, scope, nothing)
end

function plot(args...; kwargs...)
    SyncPlot(Plot(args...; kwargs...))
end

# Add some basic Julia API methods on SyncPlot that just forward onto the Plot
Base.size(sp::SyncPlot) = size(sp.plot)
Base.copy(sp::SyncPlot) = SyncPlot(copy(sp.plot))

Base.display(::PlotlyJSDisplay, p::SyncPlot) = display_blink(p::SyncPlot)

function display_blink(p::SyncPlot)
    sizeBuffer = 1.15
    plotSize = size(p.plot)
    windowOptions = Dict(
        "width" => floor(Int, plotSize[1] * sizeBuffer),
        "height" => floor(Int, plotSize[2] * sizeBuffer)
    )
    p.window = Blink.Window(windowOptions)
    Blink.body!(p.window, p.scope)
end

function Base.close(p::SyncPlot)
    if p.window !== nothing && Blink.active(p.window)
        close(p.window)
    end
end

function send_command(scope, cmd, args...)
    # The handler for _commands is set up when plot is constructed
    scope["_commands"][] = [cmd, args...]
    nothing
end

function add_trace!(p::SyncPlot, trace::GenericTrace; kw...)
    add_trace!(p.plot, trace; kw...)
    send_command(p.scope, :react, p.plot.data, p.plot.layout)
end

# ----------------------- #
# Plotly.js api functions #
# ----------------------- #

# for each of these we first update the Julia object, then update the display

function restyle!(
        plt::SyncPlot, ind::Union{Int,AbstractVector{Int}},
        update::AbstractDict=Dict();
        kwargs...
    )
    restyle!(plt.plot, ind, update; kwargs...)
    send_command(plt.scope, :restyle, merge(update, prep_kwargs(kwargs)), ind .- 1)
end

function restyle!(plt::SyncPlot, update::AbstractDict=Dict(); kwargs...)
    restyle!(plt.plot, update; kwargs...)
    send_command(plt.scope, :restyle, merge(update, prep_kwargs(kwargs)))
end

function relayout!(plt::SyncPlot, update::AbstractDict=Dict(); kwargs...)
    relayout!(plt.plot, update; kwargs...)
    update!(plt, layout=plt.plot.layout)
end

function react!(plt::SyncPlot, data::AbstractVector{<:AbstractTrace}, layout::Layout)
    react!(plt.plot, data, layout)
    send_command(plt.scope, :react, data, layout)
end

function update!(
        plt::SyncPlot, ind::Union{Int,AbstractVector{Int}},
        update::AbstractDict=Dict();
        layout::Layout=Layout(),
        kwargs...)
    update!(plt.plot, ind, update; layout=layout, kwargs...)
    send_command(plt.scope, :update, merge(update, prep_kwargs(kwargs)), layout, ind .- 1)
end

function update!(
        plt::SyncPlot, update::AbstractDict=Dict(); layout::Layout=Layout(),
        kwargs...
    )
    update!(plt.plot, update; layout=layout, kwargs...)
    send_command(plt.scope, :update, merge(update, prep_kwargs(kwargs)), layout)
end

function addtraces!(plt::SyncPlot, traces::AbstractTrace...)
    addtraces!(plt.plot, traces...)
    send_command(plt.scope, :addTraces, traces)
end

function addtraces!(plt::SyncPlot, i::Int, traces::AbstractTrace...)
    addtraces!(plt.plot, i, traces...)
    send_command(plt.scope, :addTraces, traces, i - 1)
end

function deletetraces!(plt::SyncPlot, inds::Int...)
    deletetraces!(plt.plot, inds...)
    send_command(plt.scope, :deleteTraces, collect(inds) .- 1)
end

function movetraces!(plt::SyncPlot, to_end::Int...)
    movetraces!(plt.plot, to_end...)
    send_command(plt.scope, :moveTraces, traces, collect(to_end) .- 1)
end

function movetraces!(
        plt::SyncPlot, src::AbstractVector{Int}, dest::AbstractVector{Int}
    )
    movetraces!(plt.plot, src, dest)
    send_command(plt.scope, :moveTraces, src .- 1, dest .- 1)
end

function redraw!(plt::SyncPlot)
    redraw!(plt.plot)
    send_command(plt.scope, :redraw)
end

function purge!(plt::SyncPlot)
    purge!(plt.plot)
    send_command(plt.scope, :purge)
end

function to_image(plt::SyncPlot; kwargs...)
    to_image(plt.plot)
    plt.scope["image"][] = ""  # reset
    plt.scope["_toImage"][] = Dict(kwargs)

    tries = 0
    while length(plt.scope["image"][]) == 0
        tries == 10 && error("Could not get image")
        sleep(0.25)
        tries += 1
    end
    return plt["image"].val
end

function download_image(plt::SyncPlot; kwargs...)
    download_image(plt.plot)
    plt.scope["_downloadImage"][] = Dict(kwargs)
    nothing
end

# unexported (by plotly.js) api methods
function extendtraces!(plt::SyncPlot, update::AbstractDict,
              indices::AbstractVector{Int}=[1], maxpoints=-1;)
    extendtraces!(plt.plot, update, indices, maxpoints)
    send_command(
        plt.scope, :extendTraces, prep_kwargs(update), indices .- 1, maxpoints
    )
end

function prependtraces!(plt::SyncPlot, update::AbstractDict,
               indices::AbstractVector{Int}=[1], maxpoints=-1;)
    prependtraces!(plt.plot, update, indices, maxpoints)
    send_command(
        plt.scope, :prependTraces, prep_kwargs(update), indices .- 1, maxpoints
    )
end

for f in [:restyle, :relayout, :update, :addtraces, :deletetraces,
          :movetraces, :redraw, :extendtraces, :prependtraces, :purge, :react]
    f! = Symbol(f, "!")
    @eval function $(f)(plt::SyncPlot, args...; kwargs...)
        out = SyncPlot(deepcopy(plt.plot))
        $(f!)(out, args...; kwargs...)
        out
    end
end

# add extra same extra methods we have on ::Plot for these functions
for f in (:extendtraces!, :prependtraces!)
    @eval begin
        function $(f)(p::SyncPlot, inds::Vector{Int}=[0],
                      maxpoints=-1; update...)
            d = Dict()
            for (k, v) in update
                d[k] = _tovec(v)
            end
            ($f)(p, d, inds, maxpoints)
        end

        function $(f)(p::SyncPlot, ind::Int, maxpoints=-1;
                      update...)
            ($f)(p, [ind], maxpoints; update...)
        end

        function $(f)(p::SyncPlot, update::AbstractDict, ind::Int,
                      maxpoints=-1)
            ($f)(p, update, [ind], maxpoints)
        end
    end
end

for mime in ["text/plain", "application/vnd.plotly.v1+json", "application/prs.juno.plotpane+html"]
    function Base.show(io::IO, m::MIME{Symbol(mime)}, p::SyncPlot, args...)
        show(io, m, p.plot, args...)
    end
end

PlotlyBase.savejson(sp::SyncPlot, fn::String) = PlotlyBase.savejson(sp.plot, fn)


text:
include("shared_types.jl")
include("shared_funcs.jl")
include("conv.jl")
include("arg.jl")

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule eudev_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("eudev")
JLLWrappers.@generate_main_file("eudev", UUID("35ca27e7-8b34-5b7f-bca9-bdc33f59eb06"))
end  # module eudev_jll


text:
"""
Write a function that returns the earned points in a single toss of a Darts game.

Darts is a game where players throw darts to a target.

In our particular instance of the game, the target rewards with 4 different
amounts of points, depending on where the dart lands:

If the dart lands outside the target, player earns no points (0 points).
If the dart lands in the outer circle of the target, player earns 1 point.
If the dart lands in the middle circle of the target, player earns 5 points.
If the dart lands in the inner circle of the target, player earns 10 points.

The outer circle has a radius of 10 units (This is equivalent to the total
radius for the entire target), the middle circle a radius of 5 units, and the
inner circle a radius of 1. Of course, they are all centered to the same point
(That is, the circles are concentric) defined by the coordinates (0, 0).

Write a function that given a point in the target (defined by its real
cartesian coordinates x and y), returns the correct amount earned by a dart
landing in that point.
"""

function distance(x,y)
  return sqrt((x*x) + (y*y))
end

function point(p)
  if 0 <= p <= 1
    return 10
  elseif 1 < p <= 5
    return 5
  elseif 5 < p <= 10
    return 1
  else
    return 0
  end
end

function score(x, y)
  return point(distance(x,y))
end


text:
##########################
# FunctionCell Interface #
##########################
"""
	ica(k [;kwargs])

Constructs an untrained cell that when piped data inside, calculates the ICA
projection matrix of the input data considering `k` signal components to recover.

# Keyword arguments (same as in `MultivariateStats`)
  * `alg` must be `:fastica` (as of June 2017 at least)
  * `fun` is the approx neg-entripy functor. `icafun(:tanh)`, `icafun(:tanh,a)`, `icafun(:gaus)` (default `icafun(:tanh)`)
  * `do_whiten` whether to perform pre-whitening (default `true)`
  * `maxiter` maximum number of iterations (default `100`)
  * `tol` tolerable change of W at convergence (default `1e-6`)
  * `mean` can be 0 e.g. data already centralized, `nothing` e.g. function will compute mean, or a pre-computed vector (default `nothing`)
  * `winit` initial guess of `W` (default `zeros(0,0)`)
  * `verbose` whether to display iteration information (default `false`)

Read the `MultivariateStats.jl` documentation for more information.
"""
ica(k::Int; kwargs...) = FunctionCell(ica, (k,), ModelProperties(), kwtitle("ICA", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	ica(x, k, [;kwargs])

Trains a function cell that when receiving data into will compute the idependent components of
the observations based on the projection matrix W extracted from `x`.
"""
# Training
ica(x::T where T<:CellData, k::Int; kwargs...) = ica(getx!(x), k;kwargs...)
ica(x::T where T<:AbstractVector, k::Int; kwargs...) = ica(mat(x, LearnBase.ObsDim.Constant{2}()), k; kwargs...)
ica(x::T where T<:AbstractMatrix, k::Int; kwargs...) = begin

	icadata = fit(MultivariateStats.ICA, getobs(x), k; kwargs...)

	# Build model properties
	modelprops = ModelProperties(size(icadata.W,1), 					# Wᵀx = y so size(x,1) == size(Wᵀ,2) == size(W,1)
		 		     size(icadata.W,2)						# also, size(y e.g. out,1) == size(Wᵀ,1) == size(W,w)
	)

	# Returned trained cell
	FunctionCell(ica, Model(icadata, modelprops), kwtitle("ICA",kwargs))
end

# Execution
ica(x::T where T<:CellData, model::Model{<:MultivariateStats.ICA}) =
	datacell(ica(getx!(x), model), gety(x))
ica(x::T where T<:AbstractVector, model::Model{<:MultivariateStats.ICA}) =
	ica(mat(x, LearnBase.ObsDim.Constant{2}()), model)
ica(x::T where T<:AbstractMatrix, model::Model{<:MultivariateStats.ICA}) =
	MultivariateStats.transform(model.data, getobs(x))


text:
module NearestNeighbors
	using Base.Collections
	using Distance

	export AbstractNearestNeighborTree, NearestNeighborTree

	export NaiveNeighborTree
	export KDTree #, BallTree, CoverTree
	export nearest, inball

	# export AbstractBall, Ball, EmptyBall
	# export AbstractBallTreeNode, BallTreeNode, EmptyBallTreeNode

	# export intersects

	include("generic.jl")
	include("utils.jl")
	include("naive.jl")
	include("kdtree.jl")
	# include("ball.jl")
	# include("balltree.jl")
	# TODO: Implement cover-trees
	# include("covertree.jl")
end


text:
using Bridge, Distributions, StaticArrays
using Test

n = 10000
T = 2.0
TT = 3.0
tt = 0.0:T/n:T
m = 1000
S = SVector{2,Float64}
M = SArray{Tuple{2,2},Float64,2,4}
B = SMatrix{2,2,Float64}([-1 0.1; -0.2 -1])
mu = 0.0*S([0.2, 0.3])
sigma = SMatrix{2,2,Float64}(2*[-0.212887  0.0687025
  0.193157  0.388997 ])
a = sigma*sigma'

P = LinPro(B, mu, sigma)

u = S([1.0, -0.2])
v = S([.5, 0.0])

@test norm(solve(Bridge.R3(), Bridge.b, tt, u, P) - Bridge.mu(0.0, u, T, P)) < 1e-8
@test norm(solve(Bridge.R3(), Bridge._dK, tt, 0*a, P) - Bridge.K(0.0, T, P)) < 1e-8
@test norm(inv(Bridge.solvebackward(Bridge.R3(), Bridge._dHinv, tt, 0*a, P)) - (Bridge.H(0.0, T, P))) < 1e-8
@test norm(solve(Bridge.R3(), Bridge._dPhi, tt, one(a), P) - Bridge.Phi(0.0, T, P)) < 1e-8

Ps = Bridge.LinProBridge(TT, v, P)
#@test norm(solve(Bridge.R3(), Bridge.b, tt, u, Ps) - Bridge.mu(0.0, u, T, Ps)) < 1e-8
#@test norm(solve(Bridge.R3(), Bridge._dPhi, tt, one(a), Ps) - Bridge.Phi(0.0, T, Ps))< 1e-8
@test norm(solve(Bridge.R3(), Bridge._dPhi, tt, one(a), Ps) - Bridge.Phi(0.0, T, Ps)) < 1e-8

#@show solve(Bridge.R3(), (t,x,P)->Bridge.B(t,P)*x, tt, one(a), Ps)*u

P2 = LinPro(B[1,1], mu[1], sqrt(a[1,1]))
Ps2 = Bridge.LinProBridge(TT, v[1], P2)
@test norm(solve(Bridge.R3(), Bridge.b, tt, u[1], Ps2) - Bridge.mu(0.0, u[1], T, Ps2)) < 1e-8
@test norm(solve(Bridge.R3(), Bridge._dK, tt, 0.0, Ps2) - Bridge.K(0.0, T, Ps2)) < 1e-8


text:
using DashBootstrapComponents, DashHtmlComponents

alerts = html_div([
    dbc_alert("This is a primary alert", color = "primary"),
    dbc_alert("This is a secondary alert", color = "secondary"),
    dbc_alert("This is a success alert! Well done!", color = "success"),
    dbc_alert("This is a warning alert... be careful...", color = "warning"),
    dbc_alert("This is a danger alert. Scary!", color = "danger"),
    dbc_alert("This is an info alert. Good to know!", color = "info"),
    dbc_alert("This is a light alert", color = "light"),
    dbc_alert("This is a dark alert", color = "dark"),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
macro cb_gen(name, num)
    fn = Symbol(name * "_" * string(num))
    quote
        function $(esc(fn))(ext::Dict{Symbol,Any},
                            phd::PH.PHData,
                            winfo::PH.WorkerInf,
                            niter::Int
                            )::Bool
            return true
        end
    end
end

function fake_init(ext::Dict{Symbol,Any}, phd::PHData)
    ext[:fake] = rand() + 1.0
    return
end

@testset "Callback Utilities" begin

    @cb_gen("cb", 1)
    @cb_gen("cb", 2)
    @cb_gen("cb", 3)
    @cb_gen("cb", 4)
    @cb_gen("cb", 5)
    @cb_gen("cb", 6)
    my_cbs = [cb(cb_1),
              cb(cb_2, Dict{Symbol,Any}(:b=>"TESTING!")),
              cb(cb_3, fake_init),
              cb(cb_4, fake_init, Dict{Symbol,Any}(:d => "Still a test")),
              cb("Fifth CB", cb_5, Dict{Symbol,Any}(:e => "More of a test")),
              cb("Sixth CB", cb_6, fake_init, Dict{Symbol,Any}(:f => "This is only a test")),
              ]

    phd = fake_phdata(3,4)
    for (k,cb) in enumerate(my_cbs)
        PH.add_callback(phd, cb)
        if k == 1
            @test isempty(cb.ext)
            @test cb.h === cb_1
        elseif k == 2
            @test cb.ext[:b] == "TESTING!"
            @test cb.h === cb_2
        elseif k == 3
            @test cb.ext[:fake] > 1.0
            @test cb.h === cb_3
        elseif k == 4
            @test cb.ext[:d] == "Still a test"
            @test cb.ext[:fake] > 1.0
            @test cb.h === cb_4
        elseif k == 5
            @test cb.name == "Fifth CB"
            @test cb.ext[:e] == "More of a test"
            @test cb.h === cb_5
        else
            @test cb.name == "Sixth CB"
            @test cb.ext[:f] == "This is only a test"
            @test cb.ext[:fake] > 1.0
            @test cb.h === cb_6
        end
    end

    for cb in my_cbs
        @test cb === get_callback(phd, cb.name)
        @test cb.ext === get_callback_ext(phd, cb.name)
    end

    @test_throws ErrorException get_callback(phd, "'Tis but a flesh wound")
end

@testset "Callback Termination" begin
    function my_term_cb(ext::Dict{Symbol,Any},
                        phd::PH.PHData,
                        winf::PH.WorkerInf,
                        niter::Int
                        )::Bool
        return niter < 5
    end

    the_cb = cb(my_term_cb)
    io = IOBuffer()
    Base.show(io, the_cb)
    @test String(take!(io)) == "my_term_cb"

    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              two_stage_model,
                                              PH.ScalarPenaltyParameter(2.0),
                                              atol=1e-8,
                                              rtol=1e-12,
                                              max_iter=500,
                                              timing=false,
                                              warm_start=false,
                                              callbacks=[the_cb]
                                              )

    @test err > 1e-8
    @test rerr > 1e-12
    @test n == 5
end

@testset "Callback Subproblem" begin

    struct TestSubproblem <: PH.AbstractSubproblem
        scenario::PH.ScenarioID
        names::Dict{PH.VariableID,PH.VariableInfo}
        values::Dict{PH.VariableID,Float64}
    end

    function create_test(scen::PH.ScenarioID)::TestSubproblem
        return TestSubproblem(scen,
                              Dict{PH.VariableID,PH.VariableInfo}(),
                              Dict{PH.VariableID,Float64}())
    end

    function PH.add_ph_objective_terms(ts::TestSubproblem,
                                       vids::Vector{PH.VariableID},
                                       r::Union{Float64,
                                                Dict{VariableID,Float64}},
                                       )::Nothing
        return
    end

    function PH.objective_value(ts::TestSubproblem)::Float64
        return 0.0
    end

    function PH.report_values(ts::TestSubproblem,
                              vars::Vector{PH.VariableID},
                              )::Dict{PH.VariableID,Float64}
        vals = Dict{PH.VariableID,Float64}()
        for vid in vars
            vals[vid] = ts.values[vid]
        end
        return vals
    end

    function PH.report_variable_info(ts::TestSubproblem,
                                     st::ScenarioTree
                                     )::Dict{VariableID,PH.VariableInfo}

        for node in PH.scenario_nodes(st, ts.scenario)
            stid = PH.stage(node)
            stage = PH.value(stid)
            vid = PH.VariableID(ts.scenario, stid, PH.Index(0))
            ts.names[vid] = PH.VariableInfo(string('a' + stage - 1), false)
            ts.values[vid] = rand()
        end

        return ts.names
    end

    function PH.solve_subproblem(ts::TestSubproblem)::MOI.TerminationStatusCode
        return MOI.OPTIMAL
    end

    function PH.update_ph_terms(ts::TestSubproblem,
                                w_vals::Dict{PH.VariableID,Float64},
                                xhat_vals::Dict{VariableID,Float64}
                                )::Nothing
        return
    end

    function to_apply(ts::TestSubproblem, set_values::Dict{VariableID,Float64})
        for (vid, value) in pairs(set_values)
            ts.values[vid] = value
        end
        return
    end

    function my_test_cb(ext::Dict{Symbol,Any},
                        phd::PH.PHData,
                        winf::PH.WorkerInf,
                        niter::Int
                        )::Bool

        if niter == 1

            vid0 = PH.VariableID(PH.ScenarioID(0), PH.StageID(1), PH.Index(0))
            val1 = PH.branch_value(phd, PH.ScenarioID(1), PH.StageID(1), PH.Index(0))
            values = Dict(vid0 => val1)
            ext[:value] = val1
            ext[:iter] = niter

            PH.apply_to_subproblem(to_apply,
                                   phd,
                                   winf,
                                   PH.ScenarioID(0),
                                   (values,)
                                   )
        end

        return true
    end

    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              create_test,
                                              PH.ScalarPenaltyParameter(2.0),
                                              atol=1e-8,
                                              rtol=1e-12,
                                              max_iter=500,
                                              timing=false,
                                              warm_start=false,
                                              callbacks=[PH.Callback(my_test_cb)]
                                              )

    @test err == 0.0
    @test rerr == 0.0
    @test obj == 0.0
    @test n == 3

    cb_ext = phd.callbacks[1].ext
    value = cb_ext[:value]
    @test PH.branch_value(phd, PH.ScenarioID(0), PH.StageID(1), PH.Index(0)) == value
    @test PH.branch_value(phd, PH.ScenarioID(1), PH.StageID(1), PH.Index(0)) == value
    @test cb_ext[:iter] == 1

end

@testset "Variable Fixing Callback" begin

    phd = fake_phdata(3, 4; add_leaves=false)
    winf = fake_worker_info()
    nscen = length(PH.scenarios(phd))
    lag = 2

    PH.add_callback(phd, PH.variable_fixing(lag=lag))

    phd.scenario_map[PH.scid(0)].branch_vars[PH.VariableID(PH.scid(0), PH.stid(1), PH.index(1))] = -1.0
    PH.compute_and_save_xhat(phd)

    for niter in 1:(2*lag*nscen)
        PH._execute_callbacks(phd, winf, niter)
    end
    cb_ext = phd.callbacks[1].ext

    for (xhid, xhat) in pairs(PH.consensus_variables(phd))
        if PH.name(phd, xhid) == "a1"
            @test !(xhid in cb_ext[:fixed])
        else
            @test xhid in cb_ext[:fixed]
        end
    end

    msg_count = 0
    while isready(winf.inputs[1])

        msg = take!(winf.inputs[1])
        msg_count += 1

        @test typeof(msg) == PH.SubproblemAction
        @test string(msg.action) == "fix_variables"

        s = msg.scen
        value_dict = msg.args[1]

        for (vid, value) in pairs(value_dict)
            xhid = PH.convert_to_xhat_id(phd, vid)
            xhat = phd.xhat[xhid]
            @test value == PH.value(xhat)
            @test xhid in cb_ext[:fixed]
        end
    end

    @test msg_count == nscen
end

@testset "Mean Deviation Callback" begin

    phd = fake_phdata(3, 4; add_leaves=false)
    winf = fake_worker_info()
    nscen = length(PH.scenarios(phd))
    tol = 1e-8

    PH.add_callback(phd, PH.mean_deviation(tol=tol, save_deviations=true))
    PH.compute_and_save_xhat(phd)
    user_continue = PH._execute_callbacks(phd, winf, 1)

    @test user_continue == false
    md_ext = PH.get_callback_ext(phd, "mean_deviation")
    @test isapprox(md_ext[:deviations][1], 0.0, atol=1e-8)
    @test md_ext[:tol] == tol

    phd.scenario_map[PH.scid(0)].branch_vars[PH.VariableID(PH.scid(0), PH.stid(1), PH.index(1))] = 4.0
    PH.compute_and_save_xhat(phd)
    user_continue = PH._execute_callbacks(phd, winf, 2)

    @test user_continue == true
    @test isapprox(md_ext[:deviations][2], 0.6)

    phd = fake_phdata(3, 4; add_leaves=false)
    phd.scenario_map[PH.scid(0)].branch_vars[PH.VariableID(PH.scid(0), PH.stid(1), PH.index(1))] = 4.0
    PH.add_callback(phd, PH.mean_deviation(tol=1.0, save_deviations=false))
    PH.compute_and_save_xhat(phd)
    user_continue = PH._execute_callbacks(phd, winf, 1)

    @test user_continue == false

end

macro spcb_gen(name, num)
    fn = Symbol(name * "_" * string(num))
    quote
        function $(esc(fn))(ext::Dict{Symbol,Any},
                            sp::T,
                            niter::Int,
                            scid::PH.ScenarioID
                            )::Bool where T <: AbstractSubproblem
            return
        end
    end
end

@testset "Subproblem Callback Utilities" begin
    @spcb_gen("spcb", 1)
    @spcb_gen("spcb", 2)
    @spcb_gen("spcb", 3)
    my_spcbs = [spcb(spcb_1),
                spcb(spcb_2, Dict{Symbol,Any}(:b=>"TESTING!")),
                spcb("Third CB", spcb_3, Dict{Symbol,Any}(:c => "This is a test.",
                                                          :d => "This is only a test.")
                     ),
              ]

    for (k,spcb) in enumerate(my_spcbs)
        io = IOBuffer()
        Base.show(io, spcb)
        @test String(take!(io)) == spcb.name
        if k == 1
            @test isempty(spcb.ext)
            @test spcb.h === spcb_1
        elseif k == 2
            @test length(spcb.ext) == 1
            @test spcb.ext[:b] == "TESTING!"
            @test spcb.h === spcb_2
        else
            @test length(spcb.ext) == 2
            @test spcb.ext[:c] == "This is a test."
            @test spcb.ext[:d] == "This is only a test."
            @test spcb.name == "Third CB"
            @test spcb.h === spcb_3
        end
    end
end

@testset "Subproblem Callback" begin

    function my_subproblem_callback(ext::Dict{Symbol,Any},
                                    sp::JuMPSubproblem,
                                    niter::Int,
                                    scenario_id::ScenarioID
                                    )

        if niter == 1 && PH.value(scenario_id) == 0
            ext[:check] += 1
            sp.model.ext[:check] = 5
        end

        if niter == 2 && PH.value(scenario_id) == 0
            ext[:sp_check] = sp.model.ext[:check]
        end

        if niter == 3
            for vref in values(sp.vars)
                if JuMP.name(vref) == "u"
                    JuMP.@constraint(sp.model, vref == 1.0)
                elseif JuMP.name(vref) == "x"
                    JuMP.@constraint(sp.model, vref == 3.0)
                end
            end
        end

        return
    end

    ext = Dict{Symbol,Any}(:check => 0)
    my_spcb = SubproblemCallback(my_subproblem_callback, ext)
    (n, err, rerr, obj, soln, phd) = PH.solve(PH.two_stage_tree(2),
                                              two_stage_model,
                                              PH.ScalarPenaltyParameter(1.0),
                                              atol=1e-8,
                                              rtol=1e-8,
                                              max_iter=10,
                                              report=0,
                                              timing=false,
                                              warm_start=false,
                                              subproblem_callbacks=[my_spcb]
                                              )
    @test ext[:check] == 1
    @test ext[:sp_check] == 5
    @test n == 4
    @test err < 1e-8
end


text:
## ExoplanetsSysSim/src/star.jl
## (c) 2015 Eric B. Ford

#using Distributions

@compat abstract type StarAbstract end               # Check does using StarAbstract cause a significant performance hit

immutable Star <: StarAbstract
  radius::Float64
  mass::Float64
  flux::Float64                      # relevant once have multiple stars in one target
  ld::LimbDarkeningParamAbstract         # TODO SCI DETAIL: add limb darkening param?
  id::Int64                          # id for looking up properties in stellar catalog

end
#typealias SingleStar Star
SingleStar = Star

immutable BinaryStar <: StarAbstract
  primary::Star
  secondary::Star
  orbit::Orbit
end

immutable MultipleStar <: StarAbstract                      # Will we want to handle triple, quad systems?
  component::Vector{StarAbstract}
  orbit::Vector{Orbit}
end

flux(s::Star) = s.flux                              # Demo of how to specify function behavior that depends on the derived type
flux(s::BinaryStar) = s.primary.flux + s.secondary.flux
flux(s::MultipleStar) = sum( flux, s.component)

mass(s::Star) = s.mass
mass(s::BinaryStar) = s.primary.mass + s.secondary.mass
mass(s::MultipleStar) = sum( mass, s.component)::Float64

function generate_stars(sim_param::SimParam)
   generate_star = get_function(sim_param,"generate_star")
   num_target_stars = get_int(sim_param,"num_targets_sim_pass_one")
   star_list = Array{StarAbstract}(num_target_stars)
   for i in 1:num_target_stars
     s =  generate_star(sim_param)
     star_list[i] = s
     #star_list[i] = generate_star(sim_param)
   end
  return star_list
end

function generate_star_dumb(sim_param::SimParam)
  r = rand(Uniform(0.8,1.3))::Float64
  m = rand(Normal(r,0.1))::Float64
  while m<0.0
    m = rand(Normal(r,0.1))::Float64
  end
  f = rand(Normal(1.0,0.1))::Float64
  while f<0.0
    f = 1.0+0.1*randn()
  end
  ld = LimbDarkeningParamQuadratic(0.4603,0.2291)   # TODO: Once we implement limb darkening
  return SingleStar(r,m,f,ld,0)
  #return SingleStar(r,m,f,0)
end

function test_star_constructors(sim_param::SimParam)
  star_tmp = generate_star_dumb(sim_param)
  f1 = flux(star_tmp)
  f2 = flux(BinaryStar(star_tmp,star_tmp,Orbit(10.0,0.0,0.0,0.0,0.0,0.0)))
  f4 = flux(MultipleStar([star_tmp for i in 1:4], [Orbit(10.0,0.0,0.0,0.0,0.0,0.0) for i in 1:4]) )
  # println("# Fluxes: ", f1, " ", f2, " ", f4)
  star_list = generate_stars(sim_param)
  return true
end


text:
include("../src/csd.jl")

using DataFrames
# using DelimitedFiles
using CSV

"""
Timing and backward stability testing of the CS Decomposition (special case)

First, we generate an (m+l)-by-(m+l) random matrix A
Then, we compute QR decomposition of A, A = QR
Afterwards, we choose Q1 = Q[1:m, 1:l]
Q2 = Q[m+1:m+l, 1:l], such that:

      Q1' * Q1 + Q2' * Q2 = I

Call csd function, return U, V, Z, cosine, sine.
Construct C and S.

Assume ||Q1' * Q1 + Q2' * Q2 - I|| ≈ eps,
compute:

      res_q1 = ||U * C * Z' - Q1||/(l * ||Q1|| * eps)

      res_q2 = ||V * S * Z' - Q2||/(l * ||Q2|| * eps)

      orthog_cs = ||C' * C + S' * S -I||/(l * eps)

      orthog_u = ||U' * U - I||/(m * eps)

      orthog_v = ||V' * V - I||/(p * eps)

      orthog_z = ||Z' * Z - I||/(l * eps)

argument
m: num of rows of Q1
l: num of cols of Q1, Q2 and num of rows of Q2

no return
"""

function test(m, l)
    if m > l
        throw(ArgumentError("m <= l is required"))
    else
        F = qr(randn(m+l, l))

        Q1 = @views Matrix(F.Q)[1:m, :]
        Q2 = @views Matrix(F.Q)[m+1:m+l, :]

        e = eps(Float64)
        U, V, Z, cosine, sine = csd2by1(Q1, Q2)
        if m == l
            C = Diagonal(cosine)
            S = Diagonal(sine)
        else
            C = [Diagonal(cosine[1:m]) zeros(m, l-m)]
            S = Diagonal([sine[1:m]; ones(l-m)])
        end
        res_q1 = opnorm(U * C * Z' - Q1, 1)/(l*opnorm(Q1, 1)*e)
        res_q2 = opnorm(V * S * Z' - Q2, 1)/(l*opnorm(Q2, 1)*e)
        orthog_cs = opnorm(C' * C + S' * S - I, 1)/(l*e)
        orthog_u = opnorm(U' * U - I, 1)/(m*e)
        orthog_v = opnorm(V' * V - I, 1)/(l*e)
        orthog_z = opnorm(Z' * Z - I, 1)/(l*e)
        return res_q1, res_q2, orthog_cs, orthog_u, orthog_v, orthog_z
    end
end

function test()
    df = DataFrame(m = Int64[], l = Int64[], res_q1 = Float64[],
    res_q2 = Float64[], orthog_cs = Float64[], orthog_u = Float64[],
    orthog_v = Float64[], orthog_z = Float64[])

    # case 1: m < l
    println("===============================================================")
    for i = 1:20
        a, b, c, d, e, f = test(40, 60)
        push!(df, [40 60 a b c d e f])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f = test(200, 300)
        push!(df, [200 300 a b c d e f])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f = test(600, 900)
        push!(df, [600 900 a b c d e f])
        println("---------------------------------------------------------------")
    end
    # May throw LAPACK Exception at line of svd(Q2)
    # for i = 1:20
    #     a, b, c, d, e, f = test(1000, 1500)
    #     push!(df, [1000 1500 a b c d e f])
    #     println("---------------------------------------------------------------")
    # end

    # case 2: m = l
    println("===============================================================")
    for i = 1:20
        a, b, c, d, e, f = test(50, 50)
        push!(df, [50 50 a b c d e f])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f = test(250, 250)
        push!(df, [250 250 a b c d e f])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f = test(750, 750)
        push!(df, [750 750 a b c d e f])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f = test(1250, 1250)
        push!(df, [1250 1250 a b c d e f])
        println("---------------------------------------------------------------")
    end

    CSV.write("rand_matrix_csd.csv", df, append=true)
    return df
end

# Example 2 in Chapter 3 of Jenny's thesis
function test2()
    Q1 = [1/sqrt(7) 0 0 1/sqrt(3);
          1/sqrt(7) -2/sqrt(10) -1/2 -1/(2*sqrt(3));
          1/sqrt(7) -1/sqrt(10) 1/4  3/(4*sqrt(3))]
    Q2 = [1/sqrt(7) 0 3/4 -3/(4*sqrt(3));
          1/sqrt(7) 0 0 0;
          1/sqrt(7) 1/sqrt(10) -1/4 -3/(4*sqrt(3));
          1/sqrt(7) 2/sqrt(10) -1/4 1/(4*sqrt(3))]
    m = 3; l = 4;
    e = eps(Float64)

    return U, V, Z, cosine, sine = csd2by1(Q1, Q2)
end


text:
using FunctionalCollections

include("./bench.jl")

println("\nPersistentHashMap")
println("=================\n")

const rands = rand(1:100000, 100000)
const unos = ones(Int, 100000)

pmap = PersistentHashMap{Int, Int}()
dict = Dict{Int, Int}()

for i=1:100000
    pmap = assoc(pmap, i, i)
    dict[i] = i
end

function getting(map, keys)
    function ()
        for k=keys
            map[k]
        end
    end
end

@bench "Getting (big)" 20 [getting(pmap, rands), getting(dict, rands)]
@bench "Getting (small)" 20 [getting(PersistentHashMap((1, 1)), unos)
                             getting([1 => 1], unos)]

function updating{T}(::Type{PersistentHashMap}, keys::Vector{T})
    pmap = PersistentHashMap{T, T}()
    function ()
        for i=keys
            pmap = assoc(pmap, i, i)
        end
    end
end
function updating{T}(::Type{Dict}, keys::Vector{T})
    dict = Dict{T, T}()
    function ()
        for i=keys
            dict[i] = i
        end
    end
end

@bench "Updating" 20 [updating(PersistentHashMap, rand(1:100000, 10000)),
                      updating(Dict, rand(1:100000, 10000))]