text:
using MarketData

facts("Data quality checks") do

    cl_incomplete = TimeArray(cl.timestamp[1:5], [cl.values[1:4]; NaN], cl.colnames)

    context("Determine uniformity of observations") do
        @fact cl[1:5]                         --> is_equally_spaced
        @fact cl[1:10]                        --> not(is_equally_spaced)
        @fact require_equally_spaced(cl[1:5]) --> nothing
        @fact_throws require_equally_spaced(cl[1:10])
    end

    context("Determine completeness of data") do
        @fact cl[1:5]                         --> is_equally_spaced
        @fact cl[1:10]                        --> not(is_equally_spaced)
        @fact require_equally_spaced(cl[1:5]) --> nothing
        @fact_throws require_equally_spaced(cl[1:10])
    end

    context("Determine cleanliness of data") do
        @fact cl[1:5]                 --> is_clean
        @fact cl[1:10]                --> not(is_clean)
        @fact cl_incomplete           --> not(is_clean)
        @fact require_clean(cl[1:5])  --> nothing
        @fact_throws require_clean(cl[1:10])
        @fact_throws require_clean(cl_incomplete)
    end

end


text:
@testset "Parsing" begin
    # Nice one
    molstringA = """
    C   -2.131551124300    2.286168823700    0.000000000000
    H   -1.061551124300    2.286168823700    0.000000000000
    H   -2.488213906200    1.408104616400    0.496683911300
    H   -2.488218762100    2.295059432700   -1.008766153900
    H   -2.488220057000    3.155340844300    0.512081313000"""

    # A bit messy
    molstringB = """
    C      -2.1315511243  2.2861688237    .0000000000
    H    -1.0615511243      2.2861688237     0.
    H    -2.4882139062   1.4081046164  0.4966839113
    H      -2.4882187621    2.2950594327     -1.0087661539
    H  -2.4882200570  3.1553408443      .5120813130"""

    # Messy with atomic number and masses
    molstringC = """
    6 12.011     -2.1315511243  2.2861688237    .0000000000
    H    -1.0615511243      2.2861688237     0.
    H    -2.4882139062   1.4081046164  0.4966839113
    1 1.008     -2.4882187621    2.2950594327     -1.0087661539
    H  -2.4882200570  3.1553408443      .5120813130"""

    atoms = [
        Atom(6, 12.011, [-2.1315511243,     2.2861688237,     0.0000000000])
        Atom(1,  1.008, [-1.0615511243,     2.2861688237,     0.0000000000])
        Atom(1,  1.008, [-2.4882139062,     1.4081046164,     0.4966839113])
        Atom(1,  1.008, [-2.4882187621,     2.2950594327,    -1.0087661539])
        Atom(1,  1.008, [-2.4882200570,     3.1553408443,     0.5120813130])
    ]

    @test Molecules.parse_string(molstringA, I=Int64) == atoms
    @test Molecules.parse_string(molstringB, I=Int64) == atoms
    @test Molecules.parse_string(molstringC, I=Int64) == atoms

    atoms = [
        Atom(8, 15.999,  [1.2091536548,      1.7664118189,     -0.0171613972])
        Atom(1, 1.008,  [2.1984800075,      1.7977100627,      0.0121161719])
        Atom(1, 1.008,  [0.9197881882,      2.4580185570,      0.6297938832])
    ]

    file_path = joinpath(@__DIR__, "xyz/water.xyz")
    @test Molecules.parse_file(file_path) == atoms

    atoms = Molecules.parse_string(molstringC)
    @test Molecules.get_xyz(atoms) |> chomp == molstringA
end

ext: jl


text:
module RaggedArrays

export RaggedArray , getindex, setindex!, getsubarraysize, getsubarray, length, size, print, copy

import Base.length
import Base.size
import Base.print
import Base.copy

type RaggedArray{T}
    data::Array{T,1};

    #To avoid unnecessary branching on corner cases:
    #offs[0] = 0
    #offs[n_columns+1] = length(data)
    offs::Array{Int64,1};
end

#Create initialized RaggedArray from regular array.
function RaggedArray{T}(A::Array{T,2})
    (m,n) = size(A);
    sizes = Array(Int64,(n,));
    sizes[:] = m;
    out = RaggedArray(T,sizes);
    out.data=copy(A[:]);

    return out;
end

#Create RaggedArray from Array{Array{T}}.
#Note: This will transpose indices, so that A[i][j] = ragged_A[j][i].
function RaggedArray{T}(A::Array{Array{T},1})
    n = length(A);
    sizes = Array(Int64,(n,));
    for i = 1 : n
        sizes[i] = length(A[i]);
    end

    ragged_A = RaggedArray(T,sizes);
    for i = 1 : n
        J = length(A[i]);
        for j = 1 : J
            ragged_A[j,i] = A[i][j];
        end
    end

    return ragged_A;
end

#Create uninitialized RaggedArray of subarrays with sizes given in
#"sizes" array.
function RaggedArray(T::Type,sizes::Array{Int64,1})
    totallen = sum(sizes);
    J = length(sizes);

    #Offsets array.
    offs=Array(Int64,(J+1,));
    #Contiguous data array.
    data=Array(T,(totallen,));

    #Compute offsets.
    offs[1]=0;
    for j = 1 : J
        offs[j+1]=offs[j]+sizes[j];
    end

    return RaggedArray(data,offs);
end

#Return number of columns.
function length{T}(A::RaggedArray{T})
    return length(A.offs)-1;
end

#Return (sizes,ncols) where sizes is array of column lengths.
function size{T}(A::RaggedArray{T})
    #Compute sizes array.
    n = length(A);
    sizes = Array(Int64,(n,));
    for i = 1 : n
        sizes[i] = getsubarraysize(A,i);
    end

    return (sizes,n);
end

#Convenience 2D array getter notation.
function getindex{T}(A::RaggedArray{T},i::Int64,j::Int64)
    check2dbounds(A,i,j);
    return A.data[A.offs[j]+i];
end

#Convenience 2D array setter notation.
function setindex!{T}(A::RaggedArray{T},val::T,i::Int64,j::Int64)
    check2dbounds(A,i,j);
    A.data[A.offs[j]+i]=val;
end

function print{T}(A::RaggedArray{T})
    (ms,n) = size(A);
    for i = 1 : n
        print("[");
        for j = 1 : ms[i]
            print(A[j,i]);print(" ");
        end
        print("\n");
    end
    print("]");
end

#Get j-th subarray
function getsubarray{T}(A::RaggedArray{T},j::Int64)
    check1dbounds(A,j);
    return (sub(A.data,(A.offs[j]+1):(A.offs[j+1])));
end

#Get the length of subarray j of RaggedArray A.
function getsubarraysize{T}(A::RaggedArray{T},j::Int64)
    check1dbounds(A,j);
    return (A.offs[j+1]-A.offs[j]);
end

function copy{T}(A::RaggedArray{T})
    data=copy(A.data);
    offs=copy(A.offs);
    return RaggedArray(data,offs);
end

function check1dbounds{T}(A::RaggedArray{T},j::Int64)
    if(j==length(A.offs))
        throw(BoundsError());
    end
end

function check2dbounds{T}(A::RaggedArray{T},i::Int64,j::Int64)
    jcolsz = getsubarraysize(A,j);
    if i > jcolsz
        throw(BoundsError());
    end
end

end


text:
## utilities

"""
    dimension(sys::BlochHamiltonian{D})  where {D}

Returns the spacial dimensional of a tight-binding Hamiltonian.
"""
function dimension(sys::BlochHamiltonian{D})  where {D}
    return D
end

"""
    norbs(sys::BlochHamiltonian)

Returns the number of orbitals within a unit cell of a BlochHamiltonian.
"""
function norbs(sys::BlochHamiltonian)
    return length(sys.orbitals)
end

## Initialize Hamiltonian
"""
    BlochHamiltonian(A, orbitals::Vector{Orbital}, hamk!::Function)

Given a vector of vector `A` = [a1, a2, ...], and  vector of `orbitals and a function `hamk!` , initializes a  BlochHamiltonian of dimension `D`, with no orbitals or hamiltonian function store
"""
function BlochHamiltonian(A, orbitals::Vector{Orbital}, hamk!::Function)

    return BlochHamiltonian(latticebasis(A), orbitals, hamk!)

end


text:
using LightGraphs
using LightGraphsMatching
using Test
using Cbc
using JuMP
using LinearAlgebra: I

@testset "LightGraphsMatching" begin

@testset "maximum_weight_matching" begin
    g = complete_graph(3)
    w = [
         1 2 1
         1 1 1
         3 1 1
        ]
    match = maximum_weight_matching(g, with_optimizer(Cbc.Optimizer, logLevel=0), w)
    @test match.mate[1] == 3
    @test match.weight ≈ 3

    g = complete_graph(3)
    w = zeros(3,3)
    w[1,2] = 1
    w[3,2] = 1
    w[1,3] = 1
    match = maximum_weight_matching(g,with_optimizer(Cbc.Optimizer, logLevel=0),w)
    @test match.weight ≈ 1

    g = Graph(4)
    add_edge!(g, 1,3)
    add_edge!(g, 1,4)
    add_edge!(g, 2,4)

    w =zeros(4,4)
    w[1,3] = 1
    w[1,4] = 3
    w[2,4] = 1

    match = maximum_weight_matching(g,with_optimizer(Cbc.Optimizer, logLevel=0),w)
    @test match.weight ≈ 3
    @test match.mate[1] == 4
    @test match.mate[2] == -1
    @test match.mate[3] == -1
    @test match.mate[4] == 1

    g = Graph(4)
    add_edge!(g, 1,2)
    add_edge!(g, 2,3)
    add_edge!(g, 3,1)
    add_edge!(g, 3,4)
    match = maximum_weight_matching(g,with_optimizer(Cbc.Optimizer, logLevel=0))
    @test match.weight ≈ 2
    @test match.mate[1] == 2
    @test match.mate[2] == 1
    @test match.mate[3] == 4
    @test match.mate[4] == 3

    w = zeros(4,4)
    w[1,2] = 1
    w[2,3] = 1
    w[1,3] = 1
    w[3,4] = 1

    match = maximum_weight_matching(g,with_optimizer(Cbc.Optimizer, logLevel=0), w)
    @test match.weight ≈ 2
    @test match.mate[1] == 2
    @test match.mate[2] == 1
    @test match.mate[3] == 4
    @test match.mate[4] == 3

    w = zeros(4,4)
    w[1,2] = 1
    w[2,3] = 1
    w[1,3] = 5
    w[3,4] = 1

    match = maximum_weight_matching(g,with_optimizer(Cbc.Optimizer, logLevel=0),w)
    @test match.weight ≈ 5
    @test match.mate[1] == 3
    @test match.mate[2] == -1
    @test match.mate[3] == 1
    @test match.mate[4] == -1
end

@testset "maximum_weight_maximal_matching" begin

    g = complete_bipartite_graph(2,2)
    w = zeros(4,4)
    w[1,3] = 10.
    w[1,4] = 1.
    w[2,3] = 2.
    w[2,4] = 11.
    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=with_optimizer(Cbc.Optimizer, logLevel=0))
    @test match.weight ≈ 21
    @test match.mate[1] == 3
    @test match.mate[3] == 1
    @test match.mate[2] == 4
    @test match.mate[4] == 2

    g =complete_bipartite_graph(2,4)
    w =zeros(6,6)
    w[1,3] = 10
    w[1,4] = 0.5
    w[2,3] = 11
    w[2,4] = 1
    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=with_optimizer(Cbc.Optimizer, logLevel=0))
    @test match.weight ≈ 11.5
    @test match.mate[1] == 4
    @test match.mate[4] == 1
    @test match.mate[2] == 3
    @test match.mate[3] == 2

    g =complete_bipartite_graph(2,6)
    w =zeros(8,8)
    w[1,3] = 10
    w[1,4] = 0.5
    w[2,3] = 11
    w[2,4] = 1
    w[2,5] = -1
    w[2,6] = -1
    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=with_optimizer(Cbc.Optimizer, logLevel=0), cutoff=0)
    @test match.weight ≈ 11.5
    @test match.mate[1] == 4
    @test match.mate[4] == 1
    @test match.mate[2] == 3
    @test match.mate[3] == 2

    g =complete_bipartite_graph(4,2)
    w = zeros(6,6)
    w[3,5] = 10
    w[3,6] = 0.5
    w[2,5] = 11
    w[1,6] = 1
    w[1,5] = -1

    match = maximum_weight_maximal_matching(g, w, algorithm=LPAlgorithm(), solver=with_optimizer(Cbc.Optimizer, logLevel=0), cutoff=0)
    @test match.weight ≈ 12
    @test match.mate[1] == 6
    @test match.mate[2] == 5
    @test match.mate[3] == -1
    @test match.mate[4] == -1
    @test match.mate[5] == 2
    @test match.mate[6] == 1

    g = complete_bipartite_graph(2, 2)
    w = zeros(4, 4)
    w[1, 3] = 10.
    w[1, 4] = 1.
    w[2, 3] = 2.
    w[2, 4] = 11.
    match = maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())
    @test match.weight ≈ 21
    @test match.mate[1] == 3
    @test match.mate[3] == 1
    @test match.mate[2] == 4
    @test match.mate[4] == 2

    g = complete_graph(3)
    w = zeros(3, 3)
    @test ! is_bipartite(g)
    @test_throws ErrorException maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())

    g = complete_bipartite_graph(2, 4)
    w = zeros(6, 6)
    w[1, 3] = 10
    w[1, 4] = 0.5
    w[2, 3] = 11
    w[2, 4] = 1
    match = maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())
    @test match.weight ≈ 11.5

    g = Graph(4)
    add_edge!(g, 1, 3)
    add_edge!(g, 1, 4)
    add_edge!(g, 2, 4)
    w = zeros(4, 4)
    w[1, 3] = 1
    w[1, 4] = 3
    w[2, 4] = 1
    match = maximum_weight_maximal_matching(g, w, algorithm=HungarianAlgorithm())
    @test match.weight ≈ 2

end

@testset "minimum_weight_perfect_matching" begin

    w = Dict(Edge(1,2)=> 500)
    g =Graph(2)
    add_edge!(g,1,2)
    match = minimum_weight_perfect_matching(g, w)
    @test match.mate[1] == 2

    w=Dict( Edge(1,2)=>500,
           Edge(1,3)=>600,
           Edge(2,3)=>700,
           Edge(3,4)=>100,
           Edge(2,4)=>1000)

    g = complete_graph(4)
    match = minimum_weight_perfect_matching(g, w)
    @test match.mate[1] == 2
    @test match.mate[2] == 1
    @test match.mate[3] == 4
    @test match.mate[4] == 3
    @test match.weight ≈ 600

    w = Dict(
             Edge(1, 2) => 500,
             Edge(1, 3) => 400,
             Edge(2, 3) => 300,
             Edge(3, 4) => 1000,
             Edge(2, 4) => 1000
            )
    g = complete_graph(4)
    match = minimum_weight_perfect_matching(g, w)
    @test match.mate[1] == 3
    @test match.mate[2] == 4
    @test match.mate[3] == 1
    @test match.mate[4] == 2
    @test match.weight ≈ 1400

    g =complete_bipartite_graph(2,2)
    w =Dict{Edge,Float64}()
    w[Edge(1,3)] = -10
    w[Edge(1,4)] = -0.5
    w[Edge(2,3)] = -11
    w[Edge(2,4)] = -1

    match = minimum_weight_perfect_matching(g, w)
    @test match.mate[1] == 4
    @test match.mate[4] == 1
    @test match.mate[2] == 3
    @test match.mate[3] == 2
    @test match.weight ≈ -11.5

    g = complete_graph(4)
    w = Dict{Edge,Float64}()
    w[Edge(1,3)] = 10
    w[Edge(1,4)] = 0.5
    w[Edge(2,3)] = 11
    w[Edge(2,4)] = 2
    w[Edge(1,2)] = 100

    match = minimum_weight_perfect_matching(g, w, 50)
    @test match.mate[1] == 4
    @test match.mate[4] == 1
    @test match.mate[2] == 3
    @test match.mate[3] == 2
    @test match.weight ≈ 11.5
end

end


text:
using Documenter, CanonicalTraits

makedocs(;
    modules=[CanonicalTraits],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/thautwarm/CanonicalTraits.jl/blob/{commit}{path}#L{line}",
    sitename="CanonicalTraits.jl",
    authors="thautwarm",
    assets=String[],
)

deploydocs(;
    repo="github.com/thautwarm/CanonicalTraits.jl",
)


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, DataFrames, CSV
	using Distributions, Plots, LaTeXStrings
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.16"

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	alpha, beta = 8, 2
	data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

	newAlpha, newBeta = alpha + sum(data), beta + length(data)
	post = Gamma(newAlpha, 1/newBeta)

	xGrid = quantile(post,0.01):0.001:quantile(post,0.99)
	significance = 0.9; halfAlpha = (1-significance)/2
end;

# ╔═╡ cd726796-3b30-11eb-15c8-85e4fb45fb10
coverage(l,u) = cdf(post,u) - cdf(post,l)

# ╔═╡ cd7297e8-3b30-11eb-1536-955db37fc8fc
function classicalCI(dist)
    l, u = mode(dist),mode(dist)
    bestl, bestu = l, u
    while  coverage(l,u) < significance
        l -= 0.00001; u += 0.00001
    end
    (l,u)
end

# ╔═╡ cd732cee-3b30-11eb-3079-f59edcdcac91
equalTailCI(dist) = (quantile(dist,halfAlpha), quantile(dist,1-halfAlpha))

# ╔═╡ cd83989a-3b30-11eb-21ab-8d2e63f6cee0
function highestDensityCI(dist)
    height = 0.999 * maximum(pdf.(dist,xGrid))
    l,u = mode(dist),mode(dist)
    while coverage(l,u) <= significance
        range = filter(theta -> pdf(dist,theta) > height, xGrid)
        l,u = minimum(range), maximum(range)
        height -= 0.00001
    end
    (l,u)
end

# ╔═╡ cd844a56-3b30-11eb-189e-430dff8d2db3
begin
	l1, u1 = classicalCI(post)
	l2, u2 = equalTailCI(post)
	l3, u3 = highestDensityCI(post)
	Text("Classical: $((l1,u1)), Width: $(u1-l1), Coverage: $(coverage(l1,u1))")
end

# ╔═╡ cd8fc7de-3b30-11eb-3ead-afc3afc38c5a
Text("Equal tails: $((l2,u2)), Width: $(u2-l2), Coverage: $(coverage(l2,u2))")

# ╔═╡ cd98da2a-3b30-11eb-336a-3fa4180d88fe
Text("Highest density: $((l3,u3)), Width: $(u3-l3), Coverage: $(coverage(l3,u3))")

# ╔═╡ cd9a50bc-3b30-11eb-2f56-159d2cb40db5
begin
	plot(xGrid,pdf.(post,xGrid),  yticks=(0:0.25:1.25),
		c=:black, label="Gamma Posterior Distribution",
			xlims=(1.4, 2.9), ylims=(-0.4,1.25))
	plot!([l1,u1],[-0.1,-0.1], label="Classic CI",
			c=:blue, shape=:vline, ms=16)
	plot!([l2,u2],[-0.2,-0.2], label="Equal Tail CI",
			c=:red, shape=:vline, ms=16)
	plot!([l3,u3],[-0.3,-0.3], label="Highest Density CI",
		c=:green, shape=:vline, ms=16, xlabel=L"\lambda", ylabel="Density")
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.16"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═cd726796-3b30-11eb-15c8-85e4fb45fb10
# ╠═cd7297e8-3b30-11eb-1536-955db37fc8fc
# ╠═cd732cee-3b30-11eb-3079-f59edcdcac91
# ╠═cd83989a-3b30-11eb-21ab-8d2e63f6cee0
# ╠═cd844a56-3b30-11eb-189e-430dff8d2db3
# ╠═cd8fc7de-3b30-11eb-3ead-afc3afc38c5a
# ╠═cd98da2a-3b30-11eb-336a-3fa4180d88fe
# ╠═cd9a50bc-3b30-11eb-2f56-159d2cb40db5
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using DifferentialEquations
using LinearAlgebra
using Plots
using ForwardDiff
using ProgressMeter
using DSP
using Serialization

N = 4
A = zeros(N,N)

xs = 1.0:30.0
ys = 1.0:30.0
x = zeros(length(xs),length(ys))
positions = [10.0,10.0,13.0,13.0] # Two components, one at (20.0,30.0) and one at (30.0,40.0)

stencil = [0.0  1.0 0.0;
           1.0 -4.0 1.0;
           0.0  1.0 0.0]

# generate_u_prime(x, dx)

# function generate1DStencil!(A, N, Δx)
#     # Generate the centered difference matrix manually
#     for i in 1:N, j in 1:N

#     if abs(i-j)<=1
#         (A[i,j]+=1)
#     end
#     if i==j
#         (A[i,j]-=3)
#     end

#     end

#     A = A/(Δx^2)
# end

function double_partial(x)
    dx = zeros(size(x))
    for i in 2:length(dx)-1
      dx[i] = 1*x[i-1] - 2*x[i] + 1*x[i+1]
    end
    dx[1] = -(2*x[1] + -1*x[2])
    dx[end] = -(-1*x[end-1] + 2*x[end])
end

function double_partial_2D(x)
    padded_du = conv2(x, stencil)
    du = padded_du[2:end-1,2:end-1]
    du
end

x = reshape(collect(1.0:16.0),(4,4))

function double_partial_2D_manual(x)
    nx = zeros(size(x))
    nx = convert.(eltype(x[1,1]),nx)

    for i=1:size(x,1)
        for j=1:size(x,2)
            nx[i,j] += conv_helper(i,j,x,stencil)
        end
    end
    nx
end

function conv_helper(i,j,x,stencil)
    val = 0
    for dx=-1:1
        for dy=-1:1
            ni = i+dx
            nj = j+dy
            if 1 <= ni <= size(x,1) && 1 <= nj <= size(x,2)
                val += x[ni,nj] * stencil[dx+2,dy+2]
            end
        end
    end
    val
end

f(u, p, t) = gen_forcing_2D(p)

function gen_forcing_2D(positions)
    out_vec = create_2D_gaussian_matrix(xs,ys,positions)
    out_vec
end

function heat_eq(u, p, t)
    du = double_partial_2D_manual(u) + f(u,p,t)
    du
end

function gaussian(x, mu, sig)
    return exp(-(x - mu)^2.0 / (2 * sig^2.0))
end

function gaussian_2D(loc, mu, sig)
    x, y = loc
    x₀, y₀ = mu
    σ₁, σ₂ = sig
    return gaussian(x, x₀, σ₁) * gaussian(y, y₀, σ₂)
end

function gaussian_2D_multipos(loc, positions, sig)
    x, y = loc
    σ₁, σ₂ = sig
    val = 0
    for i=1:2:length(positions)
        x₀, y₀ = positions[i], positions[i+1]
        val += gaussian(x, x₀, σ₁) * gaussian(y, y₀, σ₂)
    end
    val
end

# Loss: Given a list of positions, what is the cost / loss of the
function loss(positions)
    simulation = predict(positions)
    # display(simulation)
    final_state = simulation[end]
    # display(size(final_state)[1])
    # display(final_state)
    return norm(final_state, 2)
end

function predict(positions)
    # Simply roll this forward in time and we get our solution
    tspan = (0.0, 1.0)
    u₀ = zeros(size(xs,1),size(ys,1))
    u_dual = convert.(eltype(positions[1]),u₀)
    prob = ODEProblem(heat_eq, u_dual, tspan, positions)
    sol = solve(prob)
    # display("Finished ODE Solve")
    sol
end

# function create_2D_gaussian_matrix_orig(xs,ys,positions)
#     out = zeros(size(xs)[1],size(ys)[1])
#     out = convert.(eltype(positions[1]),out)
#     x_idx = 1
#     y_idx = 1

#     for pos in positions
#         x_p, y_p = pos
#         for x in xs
#             for y in ys
#                 out[x_idx, y_idx] += gaussian_2D([x,y],[x_p,y_p],[5.0,5.0])
#                 y_idx += 1
#             end
#             y_idx = 1
#             x_idx += 1
#         end
#         x_idx = 1

#     end
#     out
# end

function create_2D_gaussian_matrix(xs,ys,positions)
    xsz = length(xs)
    ysz = length(ys)
    out = zeros(xsz,ysz)
    out = convert.(eltype(positions[1]),out)

    for i=1:xsz
        for j=1:ysz
            x, y = xs[i], ys[j]
            out[i, j] += gaussian_2D_multipos([x,y],positions,[5.0,5.0])
        end
    end

    out
end

function plot_gaussian()
    positions = [[20.0,30.0],[40.0,50.0]]
    xs = 1.0:100.0
    ys = 1.0:100.0
    gm = create_2D_gaussian_matrix(xs,ys,positions)
    gm_f(x, y) = gm[x, y]
    display(plot(1:length(xs), 1:length(ys), gm_f, st=:heatmap, color=:viridis))
end

function plot_heatmap(data)
    display(plot(1:size(data,1), 1:size(data,2), data, st=:heatmap, color=:viridis))
end

function plot_gaussian_2()
    # positions = [20.0,30.0,40.0,50.0,66.0,3.3]
    positions = [-3.46019, 1.13722, 12.5033, 20.4116, 35.1203, 5.80054]
    xs = 1.0:0.5:30.0
    ys = 1.0:0.5:30.0
    gf = (x,y) -> gaussian_2D_multipos([x,y],positions,[5.0,5.0])
    display(plot(xs, ys, gf, st=:heatmap, color=:viridis))
end

function optimize_and_save_thermal(positions, η=0.1)
    solns = []
    append!(solns, [deepcopy(positions)])
    @showprogress for idx in 1:100
        grads = ForwardDiff.gradient(loss, positions) # Magic
        positions .-= η*grads
        # display(positions)
        # display(loss(positions))

        if loss(positions) < 0.1
            break
        end

        append!(solns, [deepcopy(positions)])

    end

    solns
end

function test_solns()
    init_pos = [10.0,10.0,13.0,13.0,20.0,10.0,20.0,20.0]
    @time ps = optimize_and_save_thermal(init_pos, 0.5)
    serialize("positions.bin", ps)

    # display("Creating gif")
    # @gif for pos in ps
    #     # gm = create_2D_gaussian_matrix(xs,ys,pos)
    #     display(pos)
    #     sol = predict(pos)[end]
    #     plot(1:size(sol,1), 1:size(sol,2), sol, st=:heatmap, color=:viridis)
    # end

end

function plot_intermediates()
    ps = deserialize("positions.bin")
    ints = []
    # Final positions: [-3.46019, 1.13722, 12.5033, 20.4116, 35.1203, 5.80054]
    for j=1:10:101
        append!(ints, [ps[j]])
    end

    i = 0
    for positions in ints
        # gm = create_2D_gaussian_matrix(xs,ys,positions)
        sol = predict(positions)[end]'
        # mf(x,y) = sol[x,y]

        display(plot(1:size(sol,1), 1:size(sol,2), sol,
                st=:heatmap,
                color=:viridis,
                title="Thermal Optimizer Iteration "*string(i*100),
                xlabel="x position (meters)",
                ylabel="y position (meters)"))
        # plot_heatmap(sol')
        savefig("final_project/plots/thermal_2D_plots/thermal_2D_plot_4C_"*string(i))
        i += 1
    end
end


text:
function trust_region(x0, deltahat, delta0, eta, f, g, h, max_iter=100, method="cauchy", tol=10e-7)

    """
    x0: vector of the initial starting point of the algorithm
    deltahat: bound on the step length
    delta0: starting trust region radius
    eta: acceptance threshold
    f: objective function
    g: gradient function
    h: Hessian function
	method: method to use for solving subproblem
        "cauchy" - Cauchy point algorithm for Newton
        "dogleg" - dogleg algorithm for Newton
        "cg_steihaug" - Steihaug algorithm, CG method
        "SR1" - symmetric rank 1 updating for quasi-Newton
    max_iter: maximum number of iterations
    tol: when to break out
    """

    @assert deltahat > 0
    @assert eta >= 0 && eta < 1/4
    @assert delta0 > 0 && delta0 < deltahat

    delta = delta0
    n = length(x0)
    xvals = Vector{Vector{Float64}}(); push!(xvals, x0)
    ncalls = 0

    for i in 1:max_iter
        i % 100 == 0 && println("Iteration: ", i)

        xcurr = xvals[end]

        fk = f(xcurr)
        gk = g(xcurr)
        Bk = h(xcurr)

        if !check_posdef(Bk)
            Bk = cholesky_mod(1e-3, Bk)
            ncalls += 1
        end

        p = solve_subproblem(delta, gk, Bk, method)
        rho = compute_rho(xcurr, f, gk, Bk, p)

        # decrease trust region radius
        if rho < 1/4
            delta *= 1/4

        else
            # increase trust region radius
            if rho > 3/4 && norm(p,2) == delta
                delta = min(2*delta, deltahat)
            end
        end

        # accept the step
        if rho > eta
            xnew = xcurr + p
        # reject the step
        else
            xnew = xcurr
        end

        push!(xvals, xnew)

        if abs(mean(xnew - xcurr)) <= tol
            println("Number of indefinite fixes ", ncalls)
            println("Number of iterations: ", i)
            return xvals
        end

    end
    println("Finished algorithm without converging.")
    return xvals

end

function model(p, xcurr, fk, gk, Bk)
    """
    p: vector in which to step in
    xcurr: vector of the current iterate
    fk: scalar, f(xk)
    gk: vector, gradient(xk)
    Bk: symmetric matrix, e.g., Hessian(xk)
    """
    return fk + gk' * p + 1/2 * p' * Bk * p
end

function compute_rho(xcurr, f, gk, Bk, p)
    """
    xcurr: vector of the current iterate
    f: the objective function
    gk: gradient(xk)
    Bk: symmetric matrix
    p: step
    """
    fk = f(xcurr)
    n = length(p)
    return ((fk - f(xcurr + p)) ./ (model(zeros(n), xcurr, fk, gk, Bk) - model(p, xcurr, fk, gk, Bk)))[1]
end

function cholesky_mod(beta, H, max_iter=1000)
    @assert beta > 0
    if minimum(diag(H)) > 0
        tau = 0
    else
        tau = -minimum(diag(H)) + beta
    end
    #while true
    for i in 1:max_iter
        candidate = H + tau * eye(H)
        check_posdef(candidate) && return candidate
        tau = max(2*tau, beta)
    end

    error("Infinite loop encountered")
end

function check_posdef(A)
    try
        L = chol(A)
        return true
    catch
        return false
    end
end

function solve_subproblem(delta, gk, Bk, method="cauchy")
	"""
    Solves the trust-region subproblem, returning a search direction pk.

    delta: trust-region radius
    gk: gradient at current iterate
    Bk: Hessian (or symmetric PD matrix) at current iterate
    method: method to use for solving subproblem
        "cauchy" - Cauchy point algorithm for Newton
        "dogleg" - dogleg algorithm for Newton
        "cg_steihaug" - Steihaug algorithm, CG method
        "SR1" - symmetric rank 1 updating for quasi-Newton
    """

    if method == "cauchy"
        return cauchy_point(delta, gk, Bk)
    elseif method == "dogleg"
        p_U = -gk' * gk ./ ((gk' * Bk * gk)' .* gk)
        p_B = -Bk \ gk
        pnorm = norm(p_B, 2)
        if pnorm > delta
            return cauchy_point(delta, gk, Bk)
        elseif pnorm < delta
            return p_B
        else
            return dogleg(delta, gk, Bk, p_U)
        end
	elseif method == "cg_steihaug"
        return cg_steihaug(delta, gk, Bk)
    elseif method == "SR1"
        println("")
    else
        println("defaulting to Cauchy point")
        return cauchy_point(delta, gk, Bk)
    end
end

function cauchy_point(delta, gk, Bk)
    """
    delta: positive scalar, trust-region radius at current iteration
    gk: gradient at current iterate
    Bk: symmetric matrix, e.g., Hessian

    Returns the cauchy point, a vector.
    """
    tau = 1
    if (gk' * Bk * gk)[1] > 0
        tau = min(norm(gk,2)^3 ./ (delta * gk'*Bk*gk), 1)[1]
    end

    # Cauchy point equation
    return -tau * delta / norm(gk, 2) * gk
end

function dogleg(delta, gk, Bk, p_U)
    p_B = -Bk \ gk
    #p_U = -gk' * gk ./ ((gk' * Bk * gk)' .* gk)

    # compute tau
    diff = p_B - p_U
    a = norm(diff, 2)^2
    b = (2 * p_U' * diff)[1]
    c = norm(p_U, 2)^2 - delta^2
    tau = (1 + (-b + sqrt(b^2 - 4 * a *c)) / (2*a))

    try
        @assert tau >=0 && tau <= 2
    catch
        println(tau)
    end

    if tau <= 1
        return tau * p_U
    else
        return p_U + (tau-1) * diff
    end

end

function cg_steihaug(delta, gk, Bk)
    gradnorm = norm(gk, 2)[1]
    tol = min(0.5, sqrt(gradnorm)) * gradnorm
    z = zeros(gk); r = gk; d = -gk

    if norm(r, 2) < tol
        return z
    end

    while true
        # negative curvature
        if (d' * Bk * d)[1] <= 0
            return cauchy_point(delta, gk, Bk)
        end
        a = (r'*r ./ (d' * Bk * d))[1]
        z += a * d

        # violates trust region bound
        if norm(z, 2) >= delta
            return cauchy_point(delta, gk, Bk)
        end
        rold = r
        r += a * Bk * d

        # convergence tolerance satisfied
        if norm(r, 2) <= tol
            return z
        end
        beta = ((r' * r) ./ (rold' * rold))[1]
        d = -r + beta * d
    end
end


text:
using SymbolicUtils
using SymbolicUtils: Symbolic, term

function symify(s::Symbolic)
    p(s) = s isa UnitRange
    r = @rule ~x::p => term(UnitRange, getproperty(~x, :start), getproperty(~x, :stop); type=UnitRange{Int})
    RW.Prewalk(RW.PassThrough(r))(s)
end

function symify(expr::Expr)
    function branch(head, newargs)
        expr = Expr(head, newargs...)
        @match expr begin
            Expr(:call, :(:), a, b) => :(Soss.term(UnitRange, $a, $b; type=UnitRange))
            # TODO: Make slices work
            :($x[$(i...)]) => :(Soss.term(getindex, $x, $(i...); type=eltype($x)))
            _ => expr
        end
    end

    return foldall(identity, branch)(expr)
end

function symify(m :: Model)
    args = m.args :: Vector{Symbol}
    vals  = map(symify, m.vals)
    dists = map(symify, m.dists)
    retn = m.retn
    Model(getmodule(m), args, vals, dists, retn)
end


text:
# Largely copied from Ti's example in ApproximateGPs.
using Pkg
Pkg.activate(@__DIR__)
using ApproximateGPs
using LinearAlgebra
using Distributions
using FillArrays
using Optim
using ParameterHandling
using Zygote
using ChainRulesCore
using PDMats
using Plots

using Random
Random.seed!(1);
ChainRulesCore.@opt_out ChainRulesCore.rrule(::Type{<:Matrix}, ::PDMats.PDMat)

# Specify parameters.
N_tr = 200;
x_tr = collect(range(-10.0, 10.0; length=N_tr));
N_z = 50;
z = collect(range(-10.0, 10.0; length=N_z))
θ_init = (
    gp = (
        σ² = positive(1.0),
        l = positive(3.0),
    ),
    z = z,
    jitter = fixed(1e-6),
);

function init_var_q(z, ::BernoulliLikelihood, x_tr)
    qω = Vector{PolyaGamma{Int,Float64}}(undef, length(x_tr))
    μ = zeros(length(z))
    Σ = Matrix{Float64}(I(length(z)))
    return μ, Σ, qω
end

μ, Σ, qω = init_var_q(z, BernoulliLikelihood(), x_tr)

# Specify functions to build LatentGP.
build_gp(θ) = GP(θ.σ² * with_lengthscale(SEKernel(), θ.l))

build_latent_gp(θ) = LatentGP(build_gp(θ.gp), BernoulliLikelihood(), θ.jitter)

# Generate synthetic data.

(f_tr, y_tr) = rand(build_latent_gp(ParameterHandling.value(θ_init))(x_tr));

# SVGP : Approximate posterior at f(z)

function build_fz_lf(θ)
    lf = build_latent_gp(θ)
    Zygote.@show lf.f.kernel
    fz = lf.f(θ.z)
    return fz, lf
end

function marginals_to_aug_posterior!(qω::AbstractVector, ps::AbstractVector{<:Normal})
    map!(qω, ps) do p
        PolyaGamma(1, sqrt(abs2(mean(p)) + var(p)))
    end
end

function aug_optimize(fz::AbstractGPs.FiniteGP, x_tr, y_tr, μ, Σ, qω; niter=3)
    K = ApproximateGPs._chol_cov(fz)
    @show κ = K \ cov(fz.f, fz.x, x_tr)
    y_sign = sign.(y_tr .- 0.5)
    @show first(μ)
    for _ in 1:niter
        pf = marginals(posterior(SparseVariationalApproximation(Centered(), fz, MvNormal(μ, Σ)))(x_tr))
        marginals_to_aug_posterior!(qω, pf)
        Σ .= Symmetric(inv(inv(K) + project(Diagonal(mean.(qω)), κ)))
        μ .= Σ * (project(y_sign / 2, κ) - K \ mean(fz))
    end
    return μ, Σ, qω
end

function loss!(θ, μ, Σ, qω)
    fz, f = build_gpzf(θ)
    fx = f(x_tr)
    μ, Σ, qω = Zygote.@ignore aug_optimize(fz, x_tr, y_tr, μ, Σ, qω)
    augsvgp = SparseVariationalApproximation(Centered(), fz, MvNormal(μ, Σ))
    return -aug_elbo(augsvgp, fx, y_tr, qω)
end

project(x::AbstractVector, κ) = κ * x
project(x::AbstractMatrix, κ) = κ * x * κ'

##
scatter(x_tr, y_tr, label="data", alpha=0.5)
plot!(x_tr, f_tr, label="")
plot!(posterior(u)(x_tr))
##
function aug_elbo(
    sva::SparseVariationalApproximation,
    lfx::AbstractGPs.LatentFiniteGP,
    y::AbstractVector,
    aug_variables;
    num_data=length(y),
)
    @assert sva.fz.f === lfx.fx.f
    return _aug_elbo(sva, lfx.fx, y, lfx.lik, aug_variables, num_data)
end

function _aug_elbo(sva, fx, y, lik, aug_variables, num_data)
    @assert sva.fz.f === fx.f
    post = posterior(sva)
    q_f = marginals(post(fx.x))
    variational_exp = expected_aug_loglik(y, q_f, lik, aug_variables)

    n_batch = length(y)
    scale = num_data / n_batch
    return sum(variational_exp) * scale - Zygote.@ignore(kl_term(lik, aug_variables)) - ApproximateGPs.kl_term(sva, post)
    # return -ApproximateGPs.kl_term(sva, post)
end

function kl_term(::BernoulliLikelihood{<:LogisticLink}, aug_variables::AbstractVector{<:PolyaGamma})
    sum(aug_variables) do qω
        c = qω.c
        - abs2(c) * mean(qω) + 2log(cosh(c / 2))
    end
end

function expected_aug_loglik(y, q_f, lik, aug_variables)
    return map(y, q_f, aug_variables) do y, q, aug
        _expected_aug_loglik(lik, q, aug, y)
    end
end

function _expected_aug_loglik(::BernoulliLikelihood{<:LogisticLink}, qf, qω, y)
    m = mean(qf)
    return  m / 2 - (abs2(m) + var(qf)) * mean(qω) / 2 - log(2)
end

# Optimise everything.
θ_flat_init, unflatten = ParameterHandling.value_flatten(θ_init);
unflatten(θ_flat_init)

function hp_loss(μ, Σ, qω)
    return θ -> loss!(unflatten(θ), μ, Σ, qω)
end

hp_loss(μ, Σ, qω)(θ_flat_init)

Zygote.gradient(hp_loss(μ, Σ, qω), θ_flat_init)

# L-BFGS parameters chosen because they seems to work well empirically.
# You could also try with the defaults.
μ, Σ, qω = init_var_q(z, BernoulliLikelihood(), x_tr)

optimisation_result = optimize(
    hp_loss(μ, Σ, qω),
    θ -> only(Zygote.gradient(hp_loss(μ, Σ, qω), θ)),
    θ_flat_init,
    LBFGS(;
        alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
        linesearch=Optim.LineSearches.BackTracking(),
    ),
    Optim.Options(; iterations=1_000);
    inplace=false,
);

θ_opt = unflatten(optimisation_result.minimizer);

#
x_pr = range(-15.0, 15.0; length=250);
fz, lf = build_gpzf(θ_opt);
svgp_opt = SparseVariationalApproximation(Centered(), fz, MvNormal(μ, Σ))
# Abstract these lines? Would be nice to write
# lf_post_opt = posterior(lf_opt, svgp_opt)
# or something.
post_opt = posterior(svgp_opt);
lf_post_opt = LatentGP(post_opt, lf.lik, lf.Σy);

f = [rand(lf_post_opt(x_pr)).f for _ in 1:20];
p = map(_f -> lf.lik.invlink.(_f), f);

let
    plt = plot()
    plot!(plt, x_pr, p; color=:blue, linealpha=0.2, label="")
    plot!(plt, x_pr, mean(p); color=:blue, linalpha=0.5, linewidth=2, label="")
    scatter!(plt, x_tr, y_tr; color=:red, label="")
end

##
u, lf = build_svgp(ParameterHandling.value(θ_opt))

m, S, qω = aug_optimize(u, x_tr, y_tr)
unew = SparseVariationalApproximation(Centered(), u.fz, MvNormal(m, S))
##
scatter(x_tr, y_tr, label="data", alpha=0.5)
plot!(x_tr, f_tr, label="")
plot!(posterior(unew)(x_tr), label="Prediction on f")
scatter!(posterior(unew)(z), label="Inducing Points")

## Comparison With AGP.jl

using AugmentedGaussianProcesses

kernel_init = build_gp(ParameterHandling.value(θ_init.gp)).kernel
agp = AGP.SVGP(kernel_init, LogisticLikelihood(), AnalyticVI(), z; verbose=3)

train!(agp, x_tr, y_tr, 1000)

μ, sig = AGP.predict_f(agp, x_tr; cov=true)

kernel_init
agp[1].prior.kernel
kernel_opt = build_gp(ParameterHandling.value(θ_opt.gp)).kernel

plot!(x_tr, μ)


text:
"Tableau for variational Gauss-Lobatto IIIA-IIIB method with two stages"
function getTableauVPLobIIIA2()
    d = [+1.0, -1.0]
    R∞ = -1
    TableauVPRK(:LobIIIA2, 2, getCoefficientsLobIIIA2(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with three stages"
function getTableauVPLobIIIA3()
    d = [+1.0, -2.0, +1.0]
    R∞ = +1
    TableauVPRK(:LobIIIA3, 4, getCoefficientsLobIIIA3(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with four stages"
function getTableauVPLobIIIA4()
    d = [+1.0, -√5, +√5, -1.0]
    R∞ = -1
    TableauVPRK(:LobIIIA4, 6, getCoefficientsLobIIIA4(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with two stages"
function getTableauVPLobIIIB2()
    d = [+1.0, -1.0]
    R∞ = -1
    TableauVPRK(:LobIIIB2, 2, getCoefficientsLobIIIB2(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with three stages"
function getTableauVPLobIIIB3()
    d = [+1.0, -2.0, +1.0]
    R∞ = +1
    TableauVPRK(:LobIIIB3, 4, getCoefficientsLobIIIB3(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIA-IIIB method with four stages"
function getTableauVPLobIIIB4()
    d = [+1.0, -√5, +√5, -1.0]
    R∞ = -1
    TableauVPRK(:LobIIIB4, 6, getCoefficientsLobIIIB4(), R∞, d)
end

"Tableau for variational Gauss-Lobatto IIIC-III method with two stages"
function getTableauVPLobIIIC2()
    R∞ = -1
    TableauVPRK(:LobIIIC2, 2, getCoefficientsLobIIIC2(), R∞)
end

"Tableau for variational Gauss-Lobatto IIIC-III method with three stages"
function getTableauVPLobIIIC3()
    R∞ = +1
    TableauVPRK(:LobIIIC3, 4, getCoefficientsLobIIIC3(), R∞)
end

"Tableau for variational Gauss-Lobatto IIIC-III method with four stages"
function getTableauVPLobIIIC4()
    R∞ = -1
    TableauVPRK(:LobIIIC4, 6, getCoefficientsLobIIIC4(), R∞)
end

"Tableau for variational Gauss-Lobatto IIID method with two stages"
function getTableauVPLobIIID2()
    lobD = getCoefficientsLobIIID2()
    R∞ = +1
    TableauVPRK(:LobIIID2, 2, lobD, R∞)
end

"Tableau for variational Gauss-Lobatto IIID method with three stages"
function getTableauVPLobIIID3()
    lobD = getCoefficientsLobIIID3()
    R∞ = -1
    TableauVPRK(:LobIIID3, 4, lobD, R∞)
end

"Tableau for variational Gauss-Lobatto IIID method with four stages"
function getTableauVPLobIIID4()
    lobD = getCoefficientsLobIIID4()
    R∞ = +1
    TableauVPRK(:LobIIID4, 6, lobD, R∞)
end

"Tableau for variational Gauss-Lobatto IIIE method with two stages"
function getTableauVPLobIIIE2()
    lobE = getCoefficientsLobIIIE2()
    R∞ = +1
    TableauVPRK(:LobIIIE2, 2, lobE, R∞)
end

"Tableau for variational Gauss-Lobatto IIIE method with three stages"
function getTableauVPLobIIIE3()
    lobE = getCoefficientsLobIIIE3()
    R∞ = -1
    TableauVPRK(:LobIIIE3, 4, lobE, R∞)
end

"Tableau for variational Gauss-Lobatto IIIE method with four stages"
function getTableauVPLobIIIE4()
    lobE = getCoefficientsLobIIIE4()
    R∞ = +1
    TableauVPRK(:LobIIIE4, 6, lobE, R∞)
end

"Tableau for variational Gauss-Lobatto IIIF method with two stages"
function getTableauVPLobIIIF2()
    lobF = getCoefficientsLobIIIF2()
    R∞ = +1
    TableauVPRK(:LobIIIF2, lobF.o, lobF, R∞)
end

"Tableau for variational Gauss-Lobatto IIIF method with three stages"
function getTableauVPLobIIIF3()
    lobF = getCoefficientsLobIIIF3()
    R∞ = -1
    TableauVPRK(:LobIIIF3, lobF.o, lobF, R∞)
end

"Tableau for variational Gauss-Lobatto IIIF method with four stages"
function getTableauVPLobIIIF4()
    lobF = getCoefficientsLobIIIF4()
    R∞ = +1
    TableauVPRK(:LobIIIF4, lobF.o, lobF, R∞)
end

"Tableau for variational Gauss-Lobatto IIIG method with two stages"
function getTableauVPLobIIIG2()
    lobG = getCoefficientsLobIIIG2()
    R∞ = +1
    TableauVPRK(:LobIIIG2, 4, lobG, lobG, R∞)
end

"Tableau for variational Gauss-Lobatto IIIG method with three stages"
function getTableauVPLobIIIG3()
    lobG = getCoefficientsLobIIIG3()
    R∞ = -1
    TableauVPRK(:LobIIIG3, 6, lobG, lobG, R∞)
end

"Tableau for variational Gauss-Lobatto IIIG method with four stages"
function getTableauVPLobIIIG4()
    lobG = getCoefficientsLobIIIG4()
    R∞ = +1
    TableauVPRK(:LobIIIG4, 8, lobG, lobG, R∞)
end

"Tableau for variational Gauss-Legendre method with s stages"
function getTableauVPGLRK(s; T=Float64)
    glrk = getCoefficientsGLRK(s, T=T)
    R∞ = (-1)^s
    TableauVPRK(Symbol("GLRK", s), glrk.o, glrk, R∞)
end

"Tableau for variational symmetric Runge-Kutta method with 3 stages"
function getTableauVPSRK3()
    srk = getCoefficientsSRK3()
    R∞ = (-1)^3
    TableauVPRK(:SRK3, srk.o, srk, R∞)
end

"Tableau for Gauss-Lobatto IIIA-IIIA method with two stages"
function getTableauVPLobIIIAIIIA2()
    R∞ = -1
    TableauVPRK(:LobIIIAIIIA2, 2, getCoefficientsLobIIIA2(), getCoefficientsLobIIIA2(), R∞)
end

"Tableau for Gauss-Lobatto IIIA-IIIA method with three stages"
function getTableauVPLobIIIAIIIA3()
    R∞ = +1
    TableauVPRK(:LobIIIAIIIA3, 4, getCoefficientsLobIIIA3(), getCoefficientsLobIIIA3(), R∞)
end

"Tableau for Gauss-Lobatto IIIA-IIIA method with four stages"
function getTableauVPLobIIIAIIIA4()
    R∞ = -1
    TableauVPRK(:LobIIIAIIIA4, 6, getCoefficientsLobIIIA4(), getCoefficientsLobIIIA4(), R∞)
end

"Tableau for Gauss-Radau IIA-IIA method with two stages"
function getTableauVPRadIIAIIA2()
    R∞ = -1
    TableauVPRK(:RadIIAIIA2, 2, getCoefficientsRadIIA2(), getCoefficientsRadIIA2(), R∞)
end

"Tableau for Gauss-Radau IIA-IIA method with three stages"
function getTableauVPRadIIAIIA3()
    R∞ = +1
    TableauVPRK(:RadIIAIIA3, 4, getCoefficientsRadIIA3(), getCoefficientsRadIIA3(), R∞)
end


text:
module Alt

using KLIEPInference
using LinearAlgebra
import JuMP, SCS, MathOptInterface, MosekTools

export
  SCS_KLIEP,
  Mosek_KLIEP

struct SCS_KLIEP <: KLIEPSolver end
struct Mosek_KLIEP <: KLIEPSolver end

KLIEP(Ψx, Ψy, ::SCS_KLIEP) =
  _KLIEP(Ψx, Ψy, JuMP.optimizer_with_attributes(SCS.Optimizer, "verbose" => 1))
KLIEP(Ψx, Ψy, ::Mosek_KLIEP) =
  _KLIEP(Ψx, Ψy, JuMP.optimizer_with_attributes(MosekTools.Mosek.Optimizer, "QUIET" => false))

spKLIEP(Ψx, Ψy, λ, ::SCS_KLIEP) =
  _spKLIEP!(Vector{Float64}(undef, size(Ψx, 1)), Ψx, Ψy, λ, JuMP.optimizer_with_attributes(SCS.Optimizer, "verbose" => 1))
spKLIEP(Ψx, Ψy, λ, ::Mosek_KLIEP) =
  _spKLIEP!(Vector{Float64}(undef, size(Ψx, 1)), Ψx, Ψy, λ, JuMP.optimizer_with_attributes(MosekTools.Mosek.Optimizer, "QUIET" => false))

### solvers

function _KLIEP(Ψx, Ψy, solver)
  m, nx = size(Ψx)
  ny = size(Ψy, 2)
  @assert m == size(Ψy, 1)

  lny = log(ny)

  problem = JuMP.Model(solver)
  JuMP.@variable(problem, Θ[1:m])
  JuMP.@variable(problem, t)
  JuMP.@variable(problem, u[1:ny])

  JuMP.@constraint(problem, sum(u) <= 1.)
  for j=1:ny
      JuMP.@constraint(problem, [dot(Θ, Ψy[:, j]) - t - lny, 1., u[j]] in MathOptInterface.ExponentialCone())
  end

  JuMP.@objective(problem, Min, -sum(Ψx'*Θ) / nx + t)
  JuMP.optimize!(problem)

  JuMP.value.(Θ)
end

function _spKLIEP!(θhat, Ψx, Ψy, λ, solver)

  m, nx = size(Ψx)
  ny = size(Ψy, 2)
  @assert m == size(Ψy, 1) == length(θhat)

  lny = log(ny)

  problem = JuMP.Model(solver)
  JuMP.@variable(problem, x[1:m])
  JuMP.@variable(problem, t)
  JuMP.@variable(problem, u[1:ny])
  JuMP.@variable(problem, b[1:m])

  # log-sum-exp constraints
  JuMP.@constraint(problem, sum(u) <= 1.)
  for j=1:ny
      JuMP.@constraint(problem, [dot(x, Ψy[:, j]) - t - lny, 1., u[j]] in MathOptInterface.ExponentialCone())
  end

  # l1 constraints
  JuMP.@constraint(problem,   x .<= b )
  JuMP.@constraint(problem,  -x .<= b )

  JuMP.@objective(problem, Min, -sum(Ψx'*x) / nx + t + λ*sum(b))
  JuMP.optimize!(problem)

  θhat .= JuMP.value.(x)
end

end


text:
@testset "baytapfile write" begin
	timeout = collect(DateTime(2010,1,1,0):Dates.Hour(1):DateTime(2010,1,1,13));
	gravout = ones(Float64,length(timeout));
	gravout[[3,end]] .= NaN;
	gravout[[4,7]] .= 10.123;
	dataout = DataFrame(datetime=timeout, grav=gravout);
	output_file = joinpath(dirname(@__DIR__),"test","output","baytap_dataseries.txt");
	writebaytap(dataout,:grav,(14.123,45.888,100.0,982.024), # position+mean gravity
				output_file,header="writebaytap unit test");
	# check manually
	open(output_file,"r") do fid
		@test readline(fid) == "writebaytap unit test"
		@test readline(fid) == "  14.12300  45.88800    100.00  982.0240"
		@test readline(fid) == " 2010    1    1   0.00000"
		@test readline(fid) == "   14    1   9000.D0"
		@test readline(fid) == "(6F10.2)"
		@test readline(fid) == "      1.00      1.00   9999.99     10.12      1.00      1.00"
		@test readline(fid) == "     10.12      1.00      1.00      1.00      1.00      1.00"
		@test readline(fid) == "      1.00   9999.99"
		@test eof(fid)
	end
end

@testset "baytapfile2tsoft conversion" begin
	file_results = joinpath(dirname(@__DIR__),"test","input","baytap08.out");
	file_output = joinpath(dirname(@__DIR__),"test","output","baytap2tsoft.txt");
	isfile(file_output) ? rm(file_output) : nothing
	baytap2tsoft(file_results,file_output,site="Cantlay",name="test",
				file_groups=joinpath(dirname(@__DIR__),"test","input","baytap_tamura_waves.txt"));
	@test isfile(file_output)
end


text:
# Autogenerated wrapper script for MSTM_jll for x86_64-w64-mingw32-libgfortran4
export mstm

using MicrosoftMPI_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("MSTM")
JLLWrappers.@declare_executable_product(mstm)
function __init__()
    JLLWrappers.@generate_init_header(MicrosoftMPI_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        mstm,
        "bin\\mstm.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MPIReco

@testset "sparse single- and multi-channel in-memory reconstruction" begin
  bSF = MPIFile("./data/mdf/systemMatrix.mdf")
  b = MPIFile("./data/mdf/measurement.mdf")
  redFactor = 0.01
  names = (:color, :x, :y, :z, :time)
  values = (1:1,
	    -21.5u"mm":1.0u"mm":21.5u"mm",
	    -21.5u"mm":1.0u"mm":21.5u"mm",
	    0.0u"mm":1.0u"mm":0.0u"mm",
	    0.0u"ms":65.28u"ms":0.0u"ms")
  valuesDF = (1:1,
	      -21.51304347826087u"mm":0.9739130434782609u"mm":-0.08695652173913043u"mm",
	      -21.51304347826087u"mm":0.9739130434782609u"mm":-0.08695652173913043u"mm",
	      -0.5u"mm":1.0u"mm":-0.5u"mm",
	      0.0u"ms":65.28u"ms":0.0u"ms")

  c1 = reconstruction(bSF, b; lambd=0.1,
		      SNRThresh = 2, frames=1:100, minFreq=80e3, nAverages=100,
		      recChannels =1:2, iterations = 3,
		      spectralLeakageCorrection=false)
  @test axisnames(c1) == names
  @test axisvalues(c1) == values
  exportImage("./img/Sparse1.png", arraydata(c1[1,:,:,1,1]))

  c2 = reconstruction(bSF, b; lambd=0.1,SNRThresh = 2, frames=1:100,
		      minFreq=80e3, nAverages=100, sparseTrafo="FFT",
		      recChannels =1:2, iterations = 3,
		      spectralLeakageCorrection=false, redFactor = redFactor,
		      useDFFoV = false)
  @test axisnames(c2) == names
  @test axisvalues(c2) == values
  exportImage("./img/Sparse2.png", arraydata(c2[1,:,:,1,1]))

  c3 = reconstruction(bSF, b; lambd=0.1, SNRThresh = 2, frames=1:100,
		      minFreq=80e3, nAverages=100, sparseTrafo="FFT",
		      recChannels =1:2, iterations = 3,
		      spectralLeakageCorrection=false, redFactor = redFactor,
		      useDFFoV = true)
  @test axisnames(c3) == names
  @test axisvalues(c3) == valuesDF
  exportImage("./img/Sparse3.png", arraydata(c3[1,:,:,1,1]))

  c4 = reconstruction(bSF, b; lambd=0.01, SNRThresh = 3, frames=1:100,
		      minFreq=80e3, nAverages=100, sparseTrafo="DCT-IV",
		      recChannels =1:2, iterations = 1,
		      spectralLeakageCorrection=false, redFactor = redFactor,
		      useDFFoV = false)
  @test axisnames(c4) == names
  @test axisvalues(c4) == values
  exportImage("./img/Sparse4.png", arraydata(c4[1,:,:,1,1]))

  c5 = reconstruction(bSF, b; lambd=0.01, SNRThresh = 3, frames=1:100,
		      minFreq=80e3, nAverages=100, sparseTrafo="DCT-IV",
		      recChannels =1:2, iterations = 1,
		      spectralLeakageCorrection=true, redFactor = redFactor,
		      useDFFoV = true)
  @test axisnames(c5) == names
  @test axisvalues(c5) == valuesDF
  exportImage("./img/Sparse5.png", arraydata(c5[1,:,:,1,1]))

  c6 = reconstruction(bSF, b; lambd=0.01, SNRThresh = 3, frames=1:100,
		      minFreq=80e3, nAverages=100, sparseTrafo="DST",
		      recChannels =1:2, iterations = 1,
		      spectralLeakageCorrection=false, redFactor = redFactor,
		      useDFFoV = false)
  @test axisnames(c6) == names
  @test axisvalues(c6) == values
  exportImage("./img/Sparse6.png", arraydata(c6[1,:,:,1,1]))

  c7 = reconstruction(bSF, b; lambd=0.01, SNRThresh = 3, frames=1:100,
		      minFreq=80e3, nAverages=100, sparseTrafo="DST",
		      recChannels =1:2, iterations = 1,
		      spectralLeakageCorrection=true, redFactor = redFactor,
		      useDFFoV = true)
  @test axisnames(c7) == names
  @test axisvalues(c7) == valuesDF
  exportImage("./img/Sparse7.png", arraydata(c7[1,:,:,1,1]))
end


text:
using TermWin
using Base.Test

@test TermWin.substr_by_width( "abc", 0, 3 ) == "abc"
@test TermWin.substr_by_width( "abc", 0, 1 ) == "a"
@test TermWin.substr_by_width( "abc", 0, -1 ) == "abc"
@test TermWin.substr_by_width( "abc", 1, -1 ) == "bc"
@test TermWin.substr_by_width( "abc", 2, -1 ) == "c"
@test TermWin.substr_by_width( "abc", 3, -1 ) == ""
@test TermWin.substr_by_width( "abc", 3, 1 ) == ""
@test TermWin.substr_by_width( "abc", 2, 1 ) == "c"

# test decorators
ohat = "o" * string( '\U302' )
ahat = "a" * string( '\U302' )
ehat = "e" * string( '\U302' )
str = ohat * ahat * ehat
@test length(str)==6
@test TermWin.substr_by_width( str, 0, 3 ) == str
@test TermWin.substr_by_width( str, 0, 1 ) == ohat
@test TermWin.substr_by_width( str, 0, -1 ) == str
@test TermWin.substr_by_width( str, 1, -1 ) == ahat * ehat
@test TermWin.substr_by_width( str, 2, -1 ) == ehat
@test TermWin.substr_by_width( str, 3, 1 ) == ""
@test TermWin.substr_by_width( str, 3, 1 ) == ""
@test TermWin.substr_by_width( str, 2, 1 ) == ehat

@test TermWin.insertstring( str, "a", 1, false ) == "a" * str
@test TermWin.insertstring( str, "a", 1, true ) == "a" * ahat * ehat
@test TermWin.insertstring( str, "a", 2, false ) == ohat * "a" * ahat * ehat
@test TermWin.insertstring( str, "a", 2, true ) == ohat * "a" * ehat
@test TermWin.insertstring( str, "a", 3, false ) == ohat * ahat * "a" * ehat
@test TermWin.insertstring( str, "a", 3, true ) == ohat * ahat * "a"

@test TermWin.delete_char_at( str, 1 ) == ahat * ehat
@test TermWin.delete_char_at( str, 2 ) == ohat * ehat
@test TermWin.delete_char_at( str, 3 ) == ohat * ahat

@test TermWin.delete_char_before( "abc", 1 ) == ("abc",1)
@test TermWin.delete_char_before( "abc", 2 ) == ("bc",1)
@test TermWin.delete_char_before( "abc", 3 ) == ("ac",2)
@test TermWin.delete_char_before( "abc", 4 ) == ("ab",3)

@test TermWin.delete_char_before( str, 1 ) == (str,1)
@test TermWin.delete_char_before( str, 2 ) == ("o" * ahat * ehat,2)
@test TermWin.delete_char_before( str, 3 ) == (ohat * "a" * ehat,3)
@test TermWin.delete_char_before( str, 4 ) == (ohat * ahat * "e",4)


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("datasets", "iris"), x=:SepalLength, y=:SepalWidth, color=:Species,
     Stat.x_jitter, Stat.y_jitter, Geom.point)


text:
"""
calculates how much `child` needs to move to not touch `parent`
"""
function move_from_touch(
        parent::GeometryPrimitive{N, T}, child::GeometryPrimitive{N},
        pad::Vec{N}
    ) where {N, T}
    pmini, cmini = minimum(parent), minimum(child) .- pad
    pmaxi, cmaxi = maximum(parent), maximum(child) .+ pad

    move = ntuple(Val{N}) do i
        posdir = ifelse(cmini[i] < pmini[i], (pmini[i] - cmini[i]), zero(T)) #always positive
        negdir = ifelse(cmaxi[i] > pmaxi[i], (pmaxi[i] - cmaxi[i]), zero(T)) #always minus
        ifelse(posdir > abs(negdir), posdir, negdir) # move in the bigger direction
    end
    Vec{N, T}(move)
end

"""
Moves `child` so that it doesn't touch parent. Leaves a gap to parent defined by `pad`.
"""
function dont_touch(
        parent::GeometryPrimitive{N}, child::GeometryPrimitive{N},
        pad::Vec{N}
    ) where N
    child + move_from_touch(parent, child, pad)
end

"""
    fit_factor_stretch(rect, lims::NTuple{N}) where N

Calculates the stretch factor to fill `rect` in all dimension.
Returns a stretch `N` dimensional fit factor.
"""
function fit_factor_stretch(rect, lims::NTuple{N, Any}) where N
    w = widths(rect)
    stretches = ntuple(Val{N}) do i
        from, to = lims[i]
        w[i] / abs(to - from)
    end
    stretches
end

"""
    fit_factor(rect, lims::NTuple{N}) where N

Calculates the scaling one needs to apply to lims to fit `rect` without changing aspect ratio.
Returns float scaling and the full strech as given by [`fit_factor_stretch`](@ref)
"""
function fit_factor(rect, lims::NTuple{N, Any}) where N
    stretches = fit_factor_stretch(rect, lims)
    minimum(stretches), stretches
end

"""
    fit_ratio(rect, lims)

Calculates the ratio one needs to stretch `lims` in order to get the same aspect ratio
"""
function fit_ratio(rect, lims)
    s, stretches = fit_factor(rect, lims)
    stretches ./ s
end

# using GeometryTypes
#
# using Base.Test
#
# x = HyperRectangle(Vec3f0(-2), Vec3f0(2))
# c = HyperRectangle(Vec3f0(-2), Vec3f0(2))
# @test dont_touch(x, c, Vec3f0(0)) == x
# c = HyperRectangle(Vec3f0(-2), Vec3f0(1.5))
# @test dont_touch(x, c, Vec3f0(0.25)) == HyperRectangle(Vec3f0(-1.75), Vec3f0(1.5))
# c = HyperRectangle(Vec3f0(0), Vec3f0(1, 1.75, 1))
# @test dont_touch(x, c, Vec3f0(0.25)) == HyperRectangle(Vec3f0(-1.25, -2.0, -1.25), Vec3f0(1.0, 1.75, 1.0))
# x = SimpleRectangle(0, 0, 1, 1)
# SimpleRectangle(HyperRectangle(x))


text:
"""

$(SIGNATURES)

Helper infrastructure to compile and sample models using `cmdstan`.
"""
module StanVariational

using StanBase
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF

import StanBase: stan_sample, get_cmdstan_home
import StanBase: cmdline, read_summary, stan_summary
import StanBase: RandomSeed, Init, Output, StanModelError

include("stanmodel/variational_types.jl")
include("stanmodel/VariationalModel.jl")
include("stanrun/cmdline.jl")
include("stansamples/read_variational.jl")

stan_variational = stan_sample

export
  VariationalModel,
  StanModelError,
  stan_variational,
  read_variational,
  read_summary,
  stan_summary

end # module


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

using Random

dir = joinpath(Sys.BINDIR, Base.DOCDIR, "examples")

include(joinpath(dir, "bubblesort.jl"))
a = rand(1:100,100)
@test issorted(sort!(a;alg=BubbleSort))

include(joinpath(dir, "lru.jl"))
include(joinpath(dir, "lru_test.jl"))

include(joinpath(dir, "ModInts.jl"))
b = ModInts.ModInt{10}(2)
c = ModInts.ModInt{10}(4)
@test b + c == ModInts.ModInt{10}(6)
@test c - b == ModInts.ModInt{10}(2)
x = ModInts.ModInt{256}(13)
y = inv(x)
@test y == ModInts.ModInt{256}(197)
@test x*y == ModInts.ModInt{256}(1)
@test_throws DomainError inv(ModInts.ModInt{8}(4))

include(joinpath(dir, "ndgrid.jl"))
r = repmat(1:10,1,10)
r1, r2 = ndgrid(1:10, 1:10)
@test r1 == r
@test r2 == r'
r3, r4 = meshgrid(1:10,1:10)
@test r3 == r'
@test r4 == r

include(joinpath(dir, "queens.jl"))
@test trysolve(8, 8, 1) == (Queen(1,1),)
@test trysolve(8, 8, 7) ==
    (Queen(1,1), Queen(2,3), Queen(3,5), Queen(4,2), Queen(5,8), Queen(7,4), Queen(8,7))

# Different cluster managers do not play well together. Since
# the test infrastructure already uses LocalManager, we will test the simple
# cluster manager example through a new Julia session.
if Sys.isunix()
    script = joinpath(dir, "clustermanager/simple/test_simple.jl")
    cmd = `$(Base.julia_cmd()) --startup-file=no $script`
    if !success(pipeline(cmd; stdout=STDOUT, stderr=STDERR)) && ccall(:jl_running_on_valgrind,Cint,()) == 0
        error("UnixDomainCM failed test, cmd : $cmd")
    end
end

using Distributed
dc_path = joinpath(dir, "dictchannel.jl")
# Run the remote on pid 1, since runtests may terminate workers
# at any time depending on memory usage
main_ex = quote
    myid() == 1 || include($dc_path)
    remotecall_fetch(1, $dc_path) do f
        include(f)
        nothing
    end
    RemoteChannel(()->DictChannel(), 1)
end
dc = eval(Main, main_ex)
@test typeof(dc) == RemoteChannel{Main.DictChannel}

@test isready(dc) == false
put!(dc, 1, 2)
put!(dc, "Hello", "World")
@test isready(dc) == true
@test isready(dc, 1) == true
@test isready(dc, "Hello") == true
@test isready(dc, 2) == false
@test fetch(dc, 1) == 2
@test fetch(dc, "Hello") == "World"
@test take!(dc, 1) == 2
@test isready(dc, 1) == false

# At least make sure code loads
include(joinpath(dir, "wordcount.jl"))

include(joinpath(dir, "clustermanager/0mq/ZMQCM.jl"))

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
using Base.Threads

"""
    nextid(g, dep::Pair)

Find the next available id such that a dead end (a node with no outgoing paths)
along the dependency chain (`dep`) is continued. If there is no such case, it
gives the maximum id (see [`walkdep`](@ref)).
"""
function nextid(g, dep::Pair)
    dep_end, cpaths = walkdep(g, dep)
    @debug "Paths compatible with the dependency chain" dep_end, cpaths
    !haskey(g.index, dep_end) && return get_prop(g)
    @debug "Last node on dependency chain found in index"
    v = g[dep_end]
    length(cpaths) == 0 && return get_prop(g)
    @debug "The number of compatible paths is not 0"
    if outdegree(g, v) > 0
        @debug "Not a dead end. Asigning path $(get_prop(g))"
        return get_prop(g)
    else
        neighbors = inneighbors(g, v)
        @debug neighbors
        # there is only one possible edge
        previ = findfirst(n->on_path(g, n, cpaths, dir=:out), neighbors)
        # check if the node is isolated and there are no ingoing edges
        previ === nothing && return get_prop(g)
        e = Edge(neighbors[previ], v)
        id = g.paths[e] ∩ cpaths
        @debug "Continuing path $id"
        # There cannot be more than one path since ids are unique and a different
        # path id would be neended only if there were a difference "further down"
        # the graph, but this is not the case since this node has no outgoing paths.
        @assert length(id) == 1
        return first(id)
    end
end

"""
    on_path(g, v, path)

Check if the vertex is on the given path.
"""
function on_path(g, v, path; dir=:in)
    !isempty(intersect!(paths_through(g, v, dir=dir), path))
end

"""
    function walkdep(g, dep::Pair; stopcond=(g,v)->false)

Walk along the dependency chain, but only on already existing paths, and return
the last node and the compatible paths.
"""
function walkdep(g, dep::Pair; stopcond=(g,v)->false)
    current_node = dep[1]
    remaining = dep[2]
    pset = Set{eltype(g)}()
    compatible_paths = paths_through(g, current_node, dir=:both)
    # @debug compatible_paths
    while !stopcond(g, current_node)
        if remaining isa Pair
            node = remaining[1]
            possible_paths = paths_through!(empty!(pset), g, node, dir=:in)
            intersect!(possible_paths, compatible_paths)
            if !isempty(possible_paths)
                current_node = node
                intersect!(compatible_paths, possible_paths)
            else
                return current_node, compatible_paths
            end
            remaining = remaining[2]
        else
            # we have reached the end of the dependency chain
            possible_paths = paths_through!(empty!(pset), g, remaining, dir=:in)
            intersect!(possible_paths, compatible_paths)
            if !isempty(possible_paths)
                return remaining, possible_paths
            else
                return current_node, compatible_paths
            end
        end
    end
    return current_node, compatible_paths
end

"""
    walkpath(g, paths, start; dir=:out, stopcond=(g,v)->false)

Walk on the given `paths` starting from `start` and return the last nodes.
If `dir` is specified, use the corresponding edge direction
(`:in` and `:out` are acceptable values).
"""
function walkpath(g, paths, start::Integer; dir=:out, kwargs...)
    if dir == :out
        walkpath(g, paths, start, outneighbors; kwargs...)
    else
        walkpath(g, paths, start, inneighbors; kwargs...)
    end
end

function walkpath(g, paths, start::Integer, neighborfn; stopcond=(g,v)->false,
        parallel_type=:threads)
    length(paths) == 0 && return Set{eltype(g)}()
    result = Vector{eltype(g)}(undef, length(paths))
    p = [paths...]
    if parallel_type == :threads
        @threads for i in eachindex(p)
            result[i] = walkpath(g, p[i], start, neighborfn, stopcond=stopcond)
        end
    else
        for i in eachindex(p)
            result[i] = walkpath(g, p[i], start, neighborfn, stopcond=stopcond)
        end
    end

    return Set{eltype(g)}(result)
end

function walkpath(g, path::Integer, start::Integer, neighborfn; stopcond=(g,v)->false)
    walkpath!(g, path, start, neighborfn, (g,v,n)->nothing, stopcond=stopcond)
end

"""
    walkpath!(g, path, start, neighborfn, action!; stopcond=(g,v)->false)

Walk on the given `path` and take an action at each node. The action is specified
by a function `action!(g, v, neighbors)` and it can modify the graph.
"""
function walkpath!(g, path, start, neighborfn, action!; stopcond=(g,v)->false)
    while !stopcond(g, start)
        neighbors = neighborfn(g, start)
        action!(g, start, neighbors)
        nexti = findfirst(n->on_path(g, n, path), neighbors)
        if nexti isa Nothing
            return start
        end
        start = neighbors[nexti]
    end
    return start
end


text:
using Base: Test

module  Dependency
  using Compat
  import Base: convert
  export abst, pabst, convert, convertsave

  @compat abstract type abst end
  @compat abstract type pabst end

  convert{P <: pabst, T <: abst}(::Type{P}, ::T) =
          getfield(T.name.module, Symbol("Packed$(T.name.name)"))
  convertsave(t) = convert(pabst, t)
end

module Extend
  using Dependency
  import Dependency: convert

  export T1, PackedT1, convertsave

  mutable struct T1 <: abst  end
  mutable struct PackedT1 <: pabst  end
end

using Extend

@test convertsave(T1()) == Extend.PackedT1

#


text:
@testset "Bump Function Tests" begin

  @test isdefined(Julz, :bump) == true

end


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# It implements the Actor-model
#

"""
    LinkParams(pid=myid(), size=32; taskref=nothing, spawn=false)

Parameters for setting up an [`Actor`](@ref).

- `pid::Int`: process identification,
- `size::Int`: channel buffer size, must be `size ≥ 10`,
- `taskref::Union{Nothing, Ref{Task}}`: If you need a reference to the created task,
    pass a `Ref{Task}` object via the keyword argument `taskref`.
- `spawn::Bool`: If spawn = true, the Task created may be scheduled on another
    thread in parallel, equivalent to creating a task via `Threads.@spawn`.
"""
struct LinkParams
    pid::Int
    size::Int
    taskref::Union{Nothing, Ref{Task}}
    spawn::Bool

    function LinkParams(pid=myid(), size=32; taskref=nothing, spawn=false)
        @assert size ≥ 10 "Link buffer size < 10 not allowed"
        new(pid, size, taskref, spawn)
    end
end

"""
    parallel(size=32; taskref=nothing)

Return [`LinkParams`](@ref) with `spawn=true`.

# Example

```julia
julia> using YAActL, .Threads

julia> myactor = Actor(parallel(), threadid);

julia> call!(myactor)
2
```
"""
parallel(size=32; taskref=nothing) = LinkParams(myid(), size, taskref=taskref, spawn=true)

"""
    Link(size=32)

Create a local Link with a buffered `Channel` `size ≥ 1`.
"""
Link(size=32) = Link(Channel{Message}(max(1, size)), myid(), :local)

# create a local Link
function _link(  func::Function, type::Symbol=:local;
                size=32, taskref=nothing, spawn=false)
    Link(
        Channel{Message}(func, size, taskref=taskref, spawn=spawn),
        myid(),
        type
        )
end

# create a remote Link
function _link(  pid::Int, func::Function;
    size=32, taskref=nothing, spawn=false)
    Link(
        RemoteChannel(()->Channel{Message}(func, size, taskref=taskref, spawn=spawn), pid),
        pid,
        :remote
    )
end

# make a remote link from a local one
_rlink(lk::Link) = lk.chn isa Channel ?
        Link(RemoteChannel(()->lk.chn),myid(),:remote) : lk
_rlink(x) = x

"User remote channel for interacting with actors."
USR = nothing

# Get the local channel to yourself from inside an actor.
_self() = current_task().code.chnl

"""
```
islocal(lk::Link)
islocal(name::Symbol)
```
Returns `true` if the actor `lk` or `name` has the same
pid as the caller, else false.
"""
islocal(lk::Link) = lk.pid == myid()
function islocal(name::Symbol)
    lk = whereis(name)
    return !ismissing(lk) ? islocal(lk) : lk
end


text:
mixer = quote
function MULT_3M5A(DO, DI_X, DI_Y; WL_X = 8, WL_Y = 8)
    WL_O = WL_X + WL_Y

    DO   = Output()[2*WL_O-1:0]
    DI_X =  Input()[2*WL_X-1:0]
    DI_Y =  Input()[2*WL_Y-1:0]

    A  := DI_X[2*WL_X-1:WL_X]
    B  := DI_X[  WL_X-1:   0]
    C  := DI_Y[2*WL_Y-1:WL_Y]
    D  := DI_Y[  WL_Y-1:   0]

    # (a+jb)(c+jd) = [c(a-b) + b(c-d)] + j[d(a+b) + b(c-d)]
    AmB_t[WL_X:0] := [A[WL_X-1] ; A] - [B[WL_X-1] ; B] # A minus B
    ApB_t[WL_X:0] := [A[WL_X-1] ; A] + [B[WL_X-1] ; B] # A plus B
    CmD_t[WL_Y:0] := [C[WL_X-1] ; C] - [D[WL_X-1] ; D] # C minus D

    CtAmB_t[WL_O-1:0] := [C[WL_Y-1](WL_O-WL_Y) ; C] * [AmB_t[WL_X](WL_O-WL_X-1) ; AmB_t]
    # C times A minus B
    DtApB_t[WL_O-1:0] := [D[WL_Y-1](WL_O-WL_Y) ; D] * [ApB_t[WL_X](WL_O-WL_X-1) ; ApB_t]
    # D times A plus B
    BtCmD_t[WL_O-1:0] := [B[WL_X-1](WL_O-WL_X) ; B] * [CmD_t[WL_Y](WL_O-WL_Y-1) ; CmD_t]
    # B times C minus D

    DO_RE[WL_O-1:0] := CtAmB_t + BtCmD_t
    DO_IM[WL_O-1:0] := DtApB_t + BtCmD_t

    DO := [DO_RE ; DO_IM]
end

function SINCOS(COS, SIN, X)
    SIN = Output()[11:0]
    COS = Output()[11:0]
    X   =  Input()[13:0]

    X_w := (X[11] == 1) ? ~X[11:0] : X[11:0] #X_w[11] = 0.1250 *2pi ; X_w[10] = 0.0625 *2pi

    X_4x_w     := X_w
    X_9d2_w    := X_w + [0;0;0;X_w[11:3]]
    X_9d4_w    := [0; X_9d2_w[11:1]]
    X_9d8_w    := [0; X_9d4_w[11:1]]

    Tmp_Sin_w  := X_4x_w + ( (X_w[10] == 1) ? X_9d8_w : X_9d4_w )

    Sin_w = Wire()[10:0]
    Sin_w[ 3:0] := Tmp_Sin_w[4:1]
    Sin_w[10:4] := Tmp_Sin_w[11:5] + [X_w[10] ;0;0; X_w[10]]

    # Cosine Wave Synthesis
    Sin_3d4_w  := [0; Sin_w[10:1]] + [0;0; Sin_w[10:2]]
    Sin_3d8_w  := [0; Sin_3d4_w[10:1]]
    Sin_3d16_w := [0; Sin_3d8_w[10:1]]
    Sin_1d16_w := [0;0;0;0; Sin_w[10:4]]

    Mux_w = Wire()[1:0]
    Mux_w[1] := Sin_w[10] | Sin_w[9]
    Mux_w[0] := Sin_w[10] | (~Sin_w[9] & Sin_w[8])

    Tmp_Cos_r = Wire()[10:0]
    @async begin
        if     Mux_w == [0;0]
            Tmp_Cos_r = Sin_1d16_w
        elseif Mux_w == [0;1]
            Tmp_Cos_r = Sin_3d16_w
        elseif Mux_w == [1;0]
            Tmp_Cos_r = Sin_3d8_w
        elseif Mux_w == [1;1]
            Tmp_Cos_r = Sin_3d4_w
        end
    end

    Tmp_Cos_w := ~Tmp_Cos_r

    C_11_w := ( Mux_w[1]) & ( Mux_w[0])
    C_10_w := ( Mux_w[1]) & (~Mux_w[0])
    C_01_w := (~Mux_w[1]) & ( Mux_w[0])

    Cos_w       = Wire()[10:0]
    Cos_w[ 4:0] := Tmp_Cos_w[4:0]
    Cos_w[10:5] := Tmp_Cos_w + [0; C_11_w; 0; C_10_w; 0; C_01_w]

    # Quadratal Conversion
    Swap_w := X[12]^X[11] # if 0.125 <= X/2pi < 0.375, sin & cos are swaped
    NegS_w := X[13]       # if 0.500 <= X/2pi        , sin is negated
    NegC_w := X[13]^X[12] # if 0.250 <= X/2pi < 0.750, cos is negated

    Sin_Swp_w := Swap_w ? Cos_w : Sin_w # MSB: 2^-1 ; LSB: 2^-11 (No sign bit yet)
    Cos_Swp_w := Swap_w ? Sin_w : Cos_w # MSB: 2^-1 ; LSB: 2^-11 (No sign bit yet)

    SIN := NegS_w ? [1 ; ~Sin_Swp_w] : [0 ; Sin_Swp_w] # Add the sign bit
    COS := NegC_w ? [1 ; ~Cos_Swp_w] : [0 ; Cos_Swp_w] # Add the sign bit
end

function MIXER(DO, DI, THETA, ENAB, CLK, RST ; WL=10, PL=2)
    # "W"ord-"L"ength
    # "P"arallelism "L"evel (Can only be power of 2, and <= 16)

    CL=12    # "C"oefficient Word"L"ength

    DO    = Output()[2*WL-1:0]
    DI    = Input()[PL*WL*2-1:0]

    THETA = Input()[13:0]

    ENAB  = Input()[0]
    CLK   = Input()[0]
    RST   = Input()[0]

    DO_RE = Wire()[WL-1:0]
    DO_IM = Wire()[WL-1:0]

    DO    = [DO_RE ; DO_IM]

    DI_r  = Wire()[WL*2-1:0][PL-1:0]

    for i = 0:PL-1
        @posedge CLK begin
            if ENAB == 1
                b_s = WL*2*i
                b_e  = WL*2*(i+1) - 1
                DI_r[i] = DI[b_e:b_s]
            end
        end
    end

    Theta_r = Wire()[13:0]
    @reg begin
        if RST == 1
            Theta_r = 0
        elseif CLK == 1
            if ENAB == 1
                Theta_r = THETA
            end
        end
    end

    Acc_r  = Wire()[13:0]
    @reg begin
        if RST == 1
            Acc_r = 0
        elseif CLK == 1
            if ENAB == 1
                Acc_r = Acc_r + (Theta_r << Int(log2(PL)))
            end
        end
    end

    cos_w = Wire()[CL-1:0][PL-1:0]
    sin_w = Wire()[CL-1:0][PL-1:0]
    cos_r = Wire()[CL-1:0][PL-1:0]
    sin_r = Wire()[CL-1:0][PL-1:0]
    re_w  = Wire()[WL+CL-1:0][PL-1:0]
    im_w  = Wire()[WL+CL-1:0][PL-1:0]

    if PL >= 4
        X3  = Wire()[13:0]
        X3  = (Theta_r<<1) + Theta_r
    end
    if PL >= 8
        X5  = Wire()[13:0]
        X7  = Wire()[13:0]
        X5  = (Theta_r<<2) + Theta_r
        X7  = (Theta_r<<3) - Theta_r
    end
    if PL == 16
        X9  = Wire()[13:0]
        X15 = Wire()[13:0]
        X9  = (Theta_r<<3) + Theta_r
        X15 = (Theta_r<<4) - Theta_r
    end

    # Figure out the order on this one
    Tmp_w = Wire()[13:0][PL-1:0]
    for i= 0:PL-1
        if i == 0
            Tmp_w[i] = Acc_r
        elseif i == 1
            Tmp_w[i] = Acc_r +  Theta_r
        elseif i == 2
            Tmp_w[i] = Acc_r +  Theta_r << 1
        elseif i == 3
            Tmp_w[i] = Acc_r +       X3
        elseif i == 4
            Tmp_w[i] = Acc_r +  Theta_r << 2
        elseif i == 5
            Tmp_w[i] = Acc_r +       X5
        elseif i == 6
            Tmp_w[i] = Acc_r +       X3 << 1
        elseif i == 7
            Tmp_w[i] = Acc_r +       X7
        elseif i == 8
            Tmp_w[i] = Acc_r +  Theta_r << 3
        elseif i == 9
            Tmp_w[i] = Acc_r +       X9
        elseif i ==10
            Tmp_w[i] = Acc_r +       X5 << 1
        elseif i ==11
            Tmp_w[i] = Acc_r + (Theta_r << 1) + X9
        elseif i ==12
            Tmp_w[i] = Acc_r +       X3 << 2
        elseif i ==13
            Tmp_w[i] = Acc_r + (Theta_r << 2) + X9
        elseif i ==14
            Tmp_w[i] = Acc_r +       X7 << 1
        elseif i ==15
            Tmp_w[i] = Acc_r +      X15
        end

        @block SINCOS "SC_$(i)" (cos_w[i], sin_w[i], Tmp_w[i])

        @posedge CLK begin
            if ENAB == 1
                cos_r[i] = cos_w[i]
            end
        end

        @posedge CLK begin
            if ENAB == 1
                sin_r[i] = sin_w[i]
            end
        end

        @block MULT_3M5A "Mult_$(i)" ([re_w[i] ; im_w[i]], DI_r[i], [Cos_r[i] ; Sin_r[i]] )

        @posedge CLK begin
            if ENAB == 1
                DO_RE[WL*(i+1)-1:WL*i] = re_w[i][WL+CL-1:CL-1]
            end
        end

        @posedge CLK begin
            if ENAB == 1
                DO_IM[WL*(i+1)-1:WL*i] = im_w[i][WL+CL-1:CL-1]
            end
        end
    end
end
end

max_stars_repo_name: trevortheblack/Juliog

max_issues_repo_name: trevortheblack/Juliog

max_forks_repo_name: trevortheblack/Juliog


text:
using JuMP
using SumOfSquares

export getlyap, soslyap, soslyapb, soslyapbs

# Storing the Lyapunov
function setlyap!(s, lyap::Lyapunov)
    d = lyap.d
    lyaps = getlyaps(s)
    if length(lyaps) < d
        sizehint!(lyaps, d)
        while length(lyaps) < d
            push!(lyaps, nothing)
        end
    end
    lyaps[d] = lyap
end
function getlyap(s::AbstractSwitchedSystem, d::Int; kws...)
    lyaps = getlyaps(s)
    if d > length(lyaps) || lyaps[d] === nothing
        soslyapb(s, d; cached=true, kws...)
    end
    return lyaps[d]
end

function getsoslyapinitub(s::AbstractDiscreteSwitchedSystem, d::Integer)
    #_, sosub = pradiusb(s, 2*d)
    #sosub
    return Inf
end
#function getsoslyapinitub(s::AbstractContinuousSwitchedSystem, d::Integer)
#    return Inf
#end

function getsoslyapinit(s, d)
    lyaps = getlyaps(s)
    if d <= length(lyaps) && lyaps[d] !== nothing
        lyap = lyaps[d]
        return lyap.soslb, lyap.dual, lyap.sosub, lyap.primal
    else
        # The SOS ub is greater than the JSR hence also greater than any of its lower bound.
        # Hence getlb(s) can be used as an initial lowerbound
        return getlb(s), nothing, getsoslyapinitub(s, d), nothing
    end
end

function superset(s::HybridSystems.AbstractHybridSystem, mode, d)
    q = GramMatrix(SOSDecomposition(variables(s, mode).^d))
    return SetProg.Sets.PolySet(2d, q)
end

# For values of γ far from 1.0, it is better to divide A_i's by γ,
# it results in a problem that is better conditioned.
# This is clearly visible in [Example 5.4, PJ08] for which the JSR is ≈ 8.9
# So we ask the user to give a scaled system, i.e. do `p(A/γx) ≥ p(x)` rather
# than do `p(Ax) ≥ γ p(x)`
"""
    soslyap(s::AbstractSwitchedSystem, d; optimizer_constructor=nothing)

Find Sum-of-Squares Lyapunov functions; i.e. solves [(5), PJ08]
or gives moment matrices certifying the infeasibility of the problem.
Use [`ScaledHybridSystem`](@ref) to use a different growth rate than 1.

[PJ08] P. Parrilo and A. Jadbabaie.
*Approximation of the joint spectral radius using sum of squares*.
Linear Algebra and its Applications, Elsevier, **2008**, 428, 2385-2402
"""
function soslyap(s::HybridSystems.AbstractHybridSystem, d; optimizer_constructor=nothing)
    sets = HybridSystems.state_property(s, PolynomialLyapunov{Float64})
    infeasibility_certificates = HybridSystems.transition_property(s, MeasureLyapunov{Float64})
    set_variables = PolySet[PolySet(symmetric=true, degree=2d, superset=superset(s, mode, d))
                            for mode in states(s)]
    status = invariant_sets!(sets, 1:nstates(s), s, optimizer_constructor, set_variables,
                             volume_heuristic = nothing,
                             infeasibility_certificates = infeasibility_certificates,
                             verbose=0)
    if isinfeasible(status, true)
        @assert !isfeasible(status, true)
        if status[3] == MOI.NO_SOLUTION
            return status, nothing, nothing
        else
            return status, nothing, infeasibility_certificates
        end
    elseif isfeasible(status, true)
        status, sets, nothing
    else
        @assert !isdecided(status, true)
        status, nothing, nothing
    end
end

function increaselb(s::AbstractDiscreteSwitchedSystem, lb, step)
    lb *= step
end

soschecktol(soslb, sosub) = sosub - soslb
soschecktol(s::AbstractDiscreteSwitchedSystem, soslb, sosub) = soschecktol(log(soslb), log(sosub))
tol_diff_str(s::AbstractDiscreteSwitchedSystem) = "Log-diff   "
#soschecktol(s::AbstractContinuousSwitchedSystem, soslb, sosub) = soschecktol(soslb, sosub)

sosshift(s::AbstractDiscreteSwitchedSystem, b, shift) = exp(log(b) + shift)
function sosshift(s::AbstractDiscreteSwitchedSystem, b, shift, scaling)
    return sosshift(s, b / scaling, shift) * scaling
end
#sosshift(s::AbstractContinuousSwitchedSystem, b, shift) = b + shift

function sosmid(soslb, sosub, step)
    if isfinite(soslb) && isfinite(sosub)
        mid = (soslb + sosub) / 2
    elseif isfinite(soslb)
        mid = soslb + step
    elseif isfinite(sosub)
        mid = sosub - step
    else
        mid = 0
    end
end
usestep(soslb, sosub) = isfinite(soslb) ⊻ isfinite(sosub)
sosmid(s::AbstractDiscreteSwitchedSystem, soslb, sosub, step) = exp(sosmid(log(soslb), log(sosub), step))
usestep(s::AbstractDiscreteSwitchedSystem, soslb, sosub) = usestep(log(soslb), log(sosub))
#sosmid(s::AbstractContinuousSwitchedSystem, soslb, sosub, step) = sosmid(soslb, sosub, step)
function sosmid(s::AbstractDiscreteSwitchedSystem, soslb, sosub, step, scaling)
    sosmid(s, soslb / scaling, sosub / scaling, step) * scaling
end

function soslb2lb(s::AbstractDiscreteSwitchedSystem, soslb, d)
    n = maximum(statedim.(s, states(s)))
    η = min(ρA(s), binomial(n+d-1, d))
    soslb / η^(1/(2*d))
end
#soslb2lb(s::AbstractContinuousSwitchedSystem, soslb, d) = -Inf

function showbs(s, soslb, sosub, tol, verbose, ok::Bool)
    if verbose >= 2 || (ok && verbose >= 1)
        println("Lower bound: $soslb")
        println("Upper bound: $sosub")
        println("$(tol_diff_str(s)): $(soschecktol(s, soslb, sosub)) $(ok ? '≤' : '>') $tol")
    end
end

function showmid(γ, status, verbose)
    if verbose >= 3
        println("  Trial value of γ: $γ")
        println("Termination status: $(status[1])")
        println("     Primal status: $(status[2])")
        println("       Dual status: $(status[3])")
        println("        Raw status: $(status[4])")
        if !isdecided(status, true)
            problem_status = "Unknown"
        elseif isfeasible(status, true)
            problem_status = "Feasible"
        else
            @assert isinfeasible(status, true)
            problem_status = "Infeasible"
        end
        println("    Problem status: $problem_status")
    end
end

# Binary Search

"""
    soslyapbs(s::AbstractSwitchedSystem, d::Integer,
              soslb, dual,
              sosub, primal;
              verbose=0, tol=1e-5, step=0.5, scaling=quickub(s),
              ranktols=tol, disttols=tol, kws...)

Find the smallest `γ` such that [`soslyap`](@ref) is feasible.
"""
function soslyapbs(s::AbstractSwitchedSystem, d::Integer,
                   soslb, dual,
                   sosub, primal;
                   verbose=0, tol=1e-5, step=0.5, scaling=quickub(s),
                   ranktols=tol, disttols=tol, kws...)
    _lyap(γ) = soslyap(ScaledHybridSystem(s, γ), d; kws...)
    while soschecktol(s, soslb, sosub) > tol
        showbs(s, soslb, sosub, tol, verbose, false)
        mid = sosmid(s, soslb, sosub, step, scaling)
        status, curprimal, curdual = _lyap(mid)
        showmid(mid, status, verbose)
        if !isdecided(status, true)
            if usestep(s, soslb, sosub)
                step *= 2
                continue
            end
            # If mid-tol/2 and mid+tol/2 also Stall, there would be an interval of length tol of Stall -> impossible to satisfy requirements
            # the distance between soslb and mid is at least tol/2.
            # Sometimes, mid is far from soslb and is at a point where the solver Stall even if it is far from the optimum point.
            # In that case, it is better to take (mid + soslb)/2
            midlb = min(sosmid(s, soslb, mid, step, scaling),
                        sosshift(s, mid, -tol/2, scaling))
            # If mid-tol/2 is too close to soslb, we would not make progress!
            # So we ensure we make a progress of at least tol/8. If dual is nothing, then that would still be progress to find a dual
            if dual !== nothing
                midlb = max(midlb, sosshift(s, soslb, tol/8))
            end
            statuslb, curprimallb, curduallb = _lyap(midlb)
            showmid(midlb, statuslb, verbose)
            if isdecided(statuslb, true)
                mid = midlb
                status = statuslb
                curprimal = curprimallb
                curdual = curduallb
            else
                midub = max(sosmid(s, mid, sosub, step, scaling),
                            sosshift(s, mid, tol/2, scaling))
                if primal !== nothing
                    midub = min(midub, sosshift(s, sosub, -tol/8))
                end
                statusub, curprimalub, curdualub = _lyap(midub)
                showmid(midub, statusub, verbose)
                if isdecided(statusub, true)
                    mid = midub
                    status = statusub
                    curprimal = curprimalub
                    curdual = curdualub
                end
            end
        end
        if isinfeasible(status, true)
            dual = curdual
            if dual !== nothing
                sosextractcycle(s, dual, d, ranktols=ranktols, disttols=disttols)
            end
            soslb = mid
        elseif isfeasible(status, true)
            if !(curprimal === nothing) # FIXME remove
                primal = curprimal
            end
            sosub = mid
        else
            @warn("Solver returned with status : $statuslb for γ=$midlb, $status for γ=$mid and $statusub for γ=$midub. Stopping bisection with $(soschecktol(s, soslb, sosub)) > $tol (= tol)")
            break
        end
    end
    if soschecktol(s, soslb, sosub) ≤ tol # it is not guaranteed because of the break
        showbs(s, soslb, sosub, tol, verbose, true)
    end
    soslb, dual, sosub, primal
end

"""
    soslyapbs(s::AbstractSwitchedSystem, d::Integer,
              soslb, dual,
              sosub, primal;
              verbose=0, tol=1e-5, step=0.5, scaling=quickub(s),
              ranktols=tol, disttols=tol, kws...)

Find upper bounds to the (constrained) Joint Spectral Radius [(5), PJ08].
Lower bounds a obtained using guarantees in [LPJ19].

[LPJ19] B. Legat, P. A. Parrilo and R. M. Jungers
*An entropy-based bound for the computational complexity of a switched system*.
IEEE Transactions on Automatic Control, **2019**.

[PJ08] P. Parrilo and A. Jadbabaie.
*Approximation of the joint spectral radius using sum of squares*.
Linear Algebra and its Applications, Elsevier, **2008**, 428, 2385-2402
"""
function soslyapb(s::AbstractSwitchedSystem, d::Integer; optimizer_constructor=nothing, tol=1e-5, cached=true, kws...)
    soslb, dual, sosub, primal = soslyapbs(s::AbstractSwitchedSystem, d::Integer, getsoslyapinit(s, d)...; optimizer_constructor=optimizer_constructor, tol=tol, kws...)
    _lyap(γ) = soslyap(ScaledHybridSystem(s, γ), d, optimizer_constructor=optimizer_constructor)
    if cached
        if primal === nothing
            if isfinite(sosub)
                status, primal, _ = _lyap(sosub)
                @assert isfeasible(status, true)
                @assert primal !== nothing
            else
                error("Bisection ended with infinite sosub=$sosub")
            end
        end
        if dual === nothing
            if isfinite(soslb)
                status, _, dual = _lyap(soslb)
                if !isinfeasible(status, true)
                    soslb = sosshift(s, soslb, -tol)
                    status, _, dual = _lyap(soslb)
                    if !isinfeasible(status, true)
                        @warn("We ignore getlb and start from scratch. tol was probably set too small and soslb is too close to the JSR so soslb-tol is too close to the JSR")
                        soslb = 0. # FIXME fix for continuous
                    end
                    soslb, dual, sosub, primal = soslyapbs(s::AbstractSwitchedSystem, d::Integer, soslb, dual, sosub, primal; optimizer_constructor=optimizer_constructor, tol=tol, kws...)
                    @assert dual !== nothing
                end
            else
                error("Bisection ended with infinite soslb=$soslb")
            end
        end
        setlyap!(s, Lyapunov(d, soslb, dual, sosub, primal))
    end
    ub = sosub
    lb = soslb2lb(s, soslb, d)
    updateb!(s, lb, ub)
end


text:
using Statistics
using LinearAlgebra
using Polynomials: Poly
import PyPlot; const plt = PyPlot

function leastsquares(x, y, Ω; order::Integer = 1)
    X = ones((length(x), order+1))
    for n in 1:order
        X[:, n+1] = x.^n
    end

    if ndims(Ω) == 1
        Ω = diagm(0 => Ω)
    end

    w = inv(Ω)  # Weights defined as 1/σ²
    Xᵀ = transpose(X)
    kω = inv(Xᵀ*w*X)*Xᵀ*w
    coeffs = kω*y
    Σ = kω*Ω*kω'
    coeffs, Σ
end

function leastsquares(x, y; order::Integer = 1)
    X = ones((length(x), order+1))
    for n in 1:order
        X[:, n+1] = x.^n
    end

    Xᵀ = transpose(X)
    inv(Xᵀ*X)*Xᵀ*y
end

function plotresiduals(x, y, coefficients; ax=nothing)
    if ax ≡ nothing
        fig, ax = plt.subplots(ncols=2)
    end
    p = Poly(coefficients)
    ŷ = p.(y)
    xrange = range(minimum(x), maximum(x), length=100)
    ŷrange = p.(xrange)

    ê = ŷ - y
    ax[1].scatter(x, y)
    ax[1].scatter(x, ŷ)
    ax[1].plot(xrange, ŷrange, "--")
    ax[1].set_title("Fit")
    ax[2].scatter(x, ê)
    ax[2].set_title("Residuals")
    return ax
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
_loss_xguide(loss::MarginLoss) = "y * h(x)"
_loss_xguide(loss::DistanceLoss) = "h(x) - y"
_loss_yguide(loss::SupervisedLoss) = "L("*_loss_xguide(loss)*")"

@recipe function plot(loss::SupervisedLoss, range=-2:0.05:2; fun=value)
    xguide --> _loss_xguide(loss)
    yguide --> _loss_yguide(loss)
    label  --> string(loss)
    l(a) = fun(loss, a)
    l, range
end


text:
"""
```
Module      : FinanceLib.FixedIncomes
Description : Implement Fixed Incomes modules for the FinanceLib library
Copyright   : (c) 2022 Kishaloy Neogi
License     : MIT
Maintainer  : Kishaloy Neogi
Email       : nkishaloy@yahoo.com
```
The module describes the base modules of FixedIncomes.

You may see the github repository at <https://github.com/n-kishaloy/FinanceLib.jl>
"""
module FixedIncomes

include("Bonds/mod.jl")
include("MoneyMarkets.jl")

end


text:
# TODO: These rules are incomplete - need to figure out how
# to mimic Julia's builtin promotion rules.
Base.promote_rule(::Type{K_real}, ::Type{Int64}) = Float64
Base.promote_rule(::Type{K_real}, ::Type{Int32}) = Float32

Base.promote_rule(::Type{K_float}, ::Type{Int64}) = Float64
Base.promote_rule(::Type{K_float}, ::Type{Int32}) = Float64

Base.promote_rule(::Type{K_long}, ::Type{Float64}) = Float64
Base.promote_rule(::Type{K_long}, ::Type{Int32}) = Int64

Base.promote_rule(::Type{K_int}, ::Type{Float64}) = Float64
Base.promote_rule(::Type{K_int}, ::Type{Int64}) = Int64

Base.promote_rule(::Type{K_short}, ::Type{Float64}) = Float64
Base.promote_rule(::Type{K_short}, ::Type{Int64}) = Int64
Base.promote_rule(::Type{K_short}, ::Type{Int32}) = Int32

# K - K rules
Base.promote_rule(::Type{K_float}, ::Type{K_float}) = Float64
Base.promote_rule(::Type{K_float}, ::Type{K_real}) = Float64
Base.promote_rule(::Type{K_float}, ::Type{K_int}) = Int64
Base.promote_rule(::Type{K_float}, ::Type{K_short}) = Int64
Base.promote_rule(::Type{K_float}, ::Type{K_boolean}) = Int64

Base.promote_rule(::Type{K_long}, ::Type{K_float}) = Float64
Base.promote_rule(::Type{K_long}, ::Type{K_real}) = Float64
Base.promote_rule(::Type{K_long}, ::Type{K_int}) = Int64
Base.promote_rule(::Type{K_long}, ::Type{K_short}) = Int64
Base.promote_rule(::Type{K_long}, ::Type{K_boolean}) = Int64


text:
module Trees

using StatsFuns
import Base.show

export Edge, Tree
export addedge!, path2root

function example(logp1, logp2)
    res = logaddexp(logp1, logp2)
    @info """
    This is simply to show how to declare dependencies:
    `logaddexp` is defined in the StatsFun package.
    If log(p1)= $logp1 and log(p2) = $logp2 then
       log(p1+p2) = $res,
    quite precisely (no underflow)
    """
    return(res)
end

include("types.jl")
include("treemanipulate.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 96d2b82a-37d0-11eb-06f7-77d3f8dbd9c9
using DataStructures

# ╔═╡ 3a2812d8-3770-11eb-34b0-7fd565a38a5f

# ╔═╡ 1a5162ec-376c-11eb-2960-5d2de63d1513
file = open(f->read(f, String), "day_4.txt")

# ╔═╡ 47348628-37cf-11eb-18cf-a3e7b88eceb2
begin
	md"separing every person passaport and its characteristics"
	person_pas = split(file, "\n\n",keepempty=false)
end

# ╔═╡ 3c3986fa-37d0-11eb-0c3b-fbcff15461da
person_pas[1]

# ╔═╡ af8f8f44-37d2-11eb-01ed-8772934a512f
md"validating function"

# ╔═╡ 43215040-37d2-11eb-35a6-4391952f9ef9
begin
	charac = ["ecl","pid","eyr","hcl","byr","iyr","hgt"]
	valid_charac(pass) = ( sum([occursin(c,pass) for c in charac]) == 7)
end

# ╔═╡ ece1fb94-37d1-11eb-2b9f-8df65e2cda06
sum([valid_charac(p)  for p in person_pas])

# ╔═╡ fcd1f052-376f-11eb-2e54-611deedd7c74
md"# Question 2"

# ╔═╡ ee13f122-37d2-11eb-2e70-c3b9fc16d1df
begin
	val_byr(yr::SubString) = 1920 ≤ parse(Int,yr) ≤ 2002
	val_iyr(yr::SubString) = 2010 ≤ parse(Int,yr) ≤ 2020
	val_eyr(yr::SubString) = 2020 ≤ parse(Int,yr) ≤ 2030
	function val_hgt(hgt::SubString)
		value = match(r"([0-9]+)",hgt).captures[1]
		if(occursin(r"cm",hgt))
			return 150 ≤ parse(Int,value) ≤ 193
		elseif (occursin(r"in",hgt))
			return 59 ≤ parse(Int,value) ≤ 76
		else return false
		end

	end
	val_hcl(hcl::SubString) = occursin(r"#[0-9a-f]{6}(?![0-9a-f])",hcl)
	val_ecl(ecl::SubString) = ecl in ("amb","blu","brn","gry","grn","hzl","oth")
	val_pid(pid::SubString) = occursin(r"^[0-9]{9}(?![0-9])",pid)
	function valid_charac2(pass)
		pass_c = [match(Regex("$(c):([^ \n]+)") ,pass).captures[1] for c in charac]
		pass_f = [val_ecl,val_pid,val_eyr,val_hcl,val_byr,val_iyr,val_hgt]
		return all([f(c) for (c,f) in zip(pass_c,pass_f)])
	end
end

# ╔═╡ 083fc8c4-37d4-11eb-0c3e-87633b5edeb6
sum([valid_charac(p) && valid_charac2(p) for p in person_pas])

# ╔═╡ 7a8b10c8-37d4-11eb-01c6-47332288862f
person_pas[1]

# ╔═╡ 04515e12-3770-11eb-2212-c55609da65f4

# ╔═╡ 58485228-3770-11eb-188a-e3f1060f6f89

# ╔═╡ 832f9fd6-3771-11eb-10cb-dd40d51e2e4c

# ╔═╡ bc646fec-3772-11eb-07a7-a963d0f51465

# ╔═╡ e59c62c0-3772-11eb-3887-c52aa4dcf981

# ╔═╡ Cell order:
# ╟─3a2812d8-3770-11eb-34b0-7fd565a38a5f
# ╠═96d2b82a-37d0-11eb-06f7-77d3f8dbd9c9
# ╠═1a5162ec-376c-11eb-2960-5d2de63d1513
# ╠═47348628-37cf-11eb-18cf-a3e7b88eceb2
# ╠═3c3986fa-37d0-11eb-0c3b-fbcff15461da
# ╠═af8f8f44-37d2-11eb-01ed-8772934a512f
# ╠═43215040-37d2-11eb-35a6-4391952f9ef9
# ╠═ece1fb94-37d1-11eb-2b9f-8df65e2cda06
# ╠═fcd1f052-376f-11eb-2e54-611deedd7c74
# ╠═ee13f122-37d2-11eb-2e70-c3b9fc16d1df
# ╠═083fc8c4-37d4-11eb-0c3e-87633b5edeb6
# ╠═7a8b10c8-37d4-11eb-01c6-47332288862f
# ╟─04515e12-3770-11eb-2212-c55609da65f4
# ╟─58485228-3770-11eb-188a-e3f1060f6f89
# ╟─832f9fd6-3771-11eb-10cb-dd40d51e2e4c
# ╟─bc646fec-3772-11eb-07a7-a963d0f51465
# ╟─e59c62c0-3772-11eb-3887-c52aa4dcf981


text:
export compute_fft

using FFTW: fft

function compute_fft(A::AbstractArray)
    in = complex(A)
    n = size(in, 1)
    out = Array{Complex{Float64}}(undef, n)
    out = fft(A)
    return out
end


text:
ambient = Scenario(
    missing,   # mass emission rate, kg/s
    missing,  # release duration, s
    missing,  # jet diameter, m
    missing, # jet velocity, m/s
    missing,   # jet density, kg/m^3
    101325.0,# release_pressure, Pa
    298.15,   # release temperature, K
    missing,   # release height, m
    missing,   # windspeed, m/s
    1.225, # ambient density, kg/m^3
    101325.0,# ambient pressure, Pa
    298.15,# ambient temperature, K
    missing,   # pasquill stability class
)

test_scenario = Scenario(
    1.0,   # mass emission rate, kg/s
    10.0,  # release duration, s
    0.25,  # jet diameter, m
    15.67, # jet velocity, m/s
    1.3,   # jet density, kg/m^3
    101325.0,# release_pressure, Pa
    450,   # release temperature, K
    1.0,   # release height, m
    1.5,   # windspeed, m/s
    1.225, # ambient density, kg/m^3
    101325.0,# ambient pressure, Pa
    298.15,# ambient temperature, K
    "F",   # pasquill stability class
)

bad_class = Scenario(
    test_scenario.mass_emission_rate,
    test_scenario.release_duration,
    test_scenario.jet_diameter,
    test_scenario.jet_velocity,
    test_scenario.jet_density,
    test_scenario.release_pressure,
    test_scenario.release_temperature,
    test_scenario.release_height,
    test_scenario.windspeed,
    test_scenario.ambient_density,
    test_scenario.ambient_pressure,
    test_scenario.ambient_temperature,
    "error",
)

# I have shamelessly stolen this from the tests for show()
replstr(x, kv::Pair...) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80), kv...), MIME("text/plain"), x), x)
test_scenario_str = "Release scenario:\n    mass_emission_rate: 1.0 kg/s \n    release_duration: 10.0 s \n    jet_diameter: 0.25 m \n    jet_velocity: 15.67 m/s \n    jet_density: 1.3 kg/m^3 \n    release_pressure: 101325.0 Pa \n    release_temperature: 450 K \n    release_height: 1.0 m \n    windspeed: 1.5 m/s \n    ambient_density: 1.225 kg/m^3 \n    ambient_pressure: 101325.0 Pa \n    ambient_temperature: 298.15 K \n    pasquill_gifford: F  \n    "

@testset "Scenario constructor tests" begin
    d = Dict{Symbol,Union{Missing, Number, String}}([
        :release_pressure => 101325.0,   # Pa
        :release_temperature => 298.15,# K
        :ambient_density => 1.225,     # kg/m^3
        :ambient_pressure => 101325.0,   # Pa
        :ambient_temperature => 298.15,# K
    ])

    @test Scenario(d) == ambient

    @test Scenario(test_scenario) == test_scenario
    @test Scenario(test_scenario; pasquill_gifford="error") == bad_class

    @test replstr(test_scenario) == test_scenario_str

end

@testset "Scenario Builder tests" begin
    # Liquid jet example, *Guidelines for Consequence Analysis of Chemical
    # Releases* CCPS, 1999, pg 40
    jet = Dict([
        :mass_emission_rate => 0.21691154763598,
        :jet_diameter => 0.01,
        :jet_velocity => 5.636333880812954,
        :jet_density => 490.0,
        :release_pressure => 120935.0368,
        :release_temperature => 298.15,
        :release_height => 0.0,
        :windspeed => 1.5,
        :ambient_density => 1.225,
        :ambient_pressure => 101325.0,
        :ambient_temperature => 298.15,
        :pasquill_gifford => "F"
    ])

    # using default ambient properties
    @test Scenario(jet) == scenario_builder(jet[:release_pressure], jet[:release_temperature];
                           model=:jet, phase=:liquid,
                           stability=jet[:pasquill_gifford], windspeed=jet[:windspeed],
                           liquid_density=jet[:jet_density], hole_diameter=jet[:jet_diameter],
                           discharge_coeff=0.63)

   # getting ambient properties from a dict
   @test Scenario(jet) == scenario_builder(jet[:release_pressure], jet[:release_temperature], jet;
                          model=:jet, phase=:liquid,
                          stability=jet[:pasquill_gifford], windspeed=jet[:windspeed],
                          liquid_density=jet[:jet_density], hole_diameter=jet[:jet_diameter],
                          discharge_coeff=0.63)

    # getting ambient properties from an existing scenario
    js = Scenario(jet)
    @test Scenario(jet) == scenario_builder(jet[:release_pressure], jet[:release_temperature], js;
                         model=:jet, phase=:liquid,
                         stability=jet[:pasquill_gifford], windspeed=jet[:windspeed],
                         liquid_density=jet[:jet_density], hole_diameter=jet[:jet_diameter],
                         discharge_coeff=0.63)

    # missing windspeed, using default
    s = scenario_builder(jet[:release_pressure], jet[:release_temperature];
                           model=:jet, phase=:liquid,
                           stability=jet[:pasquill_gifford], #windspeed=missing,
                           liquid_density=jet[:jet_density], hole_diameter=jet[:jet_diameter],
                           discharge_coeff=0.63)
    @test s.windspeed == 3.0

    # missing model params, throw error
    @test_throws ErrorException scenario_builder(jet[:release_pressure], jet[:release_temperature];
                           model=:jet, phase=:liquid,
                           stability=jet[:pasquill_gifford], windspeed=missing,
                           liquid_density=jet[:jet_density], #hole_diameter=jet[:jet_diameter],
                           discharge_coeff=0.63)

    # missing model, throw error
    @test_throws ErrorException scenario_builder(jet[:release_pressure], jet[:release_temperature];
                           model=:somethingelse, phase=:liquid,
                           stability=jet[:pasquill_gifford], windspeed=jet[:windspeed])

end


text:
connectivity(network::RegularNetwork) = network.k

function isregular(network::AbstractNetwork)
	(typeof(network) <: RegularNetwork) && (return true)

	mat = adjMat(network)
	k = sum(mat[:,1])
	return mat == adjMat(RegularNetwork(network.N, k))
end

# ---------- Calculators ----------
calcNumConnections!(network::RegularNetwork) =
	network._props.numConnections = Int(network.k * network.N / (isdirected(network) ? 1 : 2))

function calcConnectivity(network::RegularNetwork; dir_behaviour::Symbol)
	isdirected || (return network.k)

	(dir_behaviour == :total) && (return 2 * network.k)
	(dir_behaviour == :both) && (return (network.k, network.k))
	return network.k
end

function calcAdjMat!(network::RegularNetwork)
	if network.k >= network.N-1
		network._adjMat = adjMat(GlobalNetwork(network.N))
		return
	end

	network._adjMat = adjMat(EmptyNetwork(network.N), sparse=true)
	(network.k <= 0) && (return)

	for i in 1:network.N
		c = 1
		while c <= floor(Int, network.k / 2)
			network._adjMat[i, i_plus_x(i,  c, network.N)] = 1
			network._adjMat[i, i_plus_x(i, -c, network.N)] = 1
			c += 1
		end
	end
end


text:
function measure(data, pattern)
    start = time()
    count = length(collect(eachmatch(pattern, data)))
    elapsed = time() - start
    elapsed_ms = 1000 * elapsed

    println(string(elapsed_ms) * " - " * string(count))
end

if length(ARGS) < 1
    println("Usage: julia benchmark.jl <filename>")
    exit(1)
end

data = open(f->read(f, String), ARGS[1])

measure(data, r"[\w.+-]+@[\w.-]+\.[\w.-]+")
measure(data, r"[\w]+://[^/\s?#]+[^\s?#]+(?:\?[^\s#]*)?(?:#[^\s]*)?")
measure(data, r"(?:(?:25[0-5]|2[0-4][0-9]|[01]?[0-9][0-9])\.){3}(?:25[0-5]|2[0-4][0-9]|[01]?[0-9][0-9])")

max_stars_repo_name: mariomka/regex-benchmark