text:
function semiconttest(model::MOI.ModelLike, config::MOIT.TestConfig{T}) where T
    atol = config.atol
    rtol = config.rtol

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.MathOptInterface.Semicontinuous{T})

    # 2 variables
    # min  x
    # st   x >= y
    #      x ∈ {0.0} U [2.0,3.0]
    #      y = 0.0

    MOI.empty!(model)
    @test MOI.is_empty(model)

    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    vc1 = MOI.add_constraint(model, MOI.SingleVariable(v[1]), MOI.Semicontinuous(T(2), T(3)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Semicontinuous{T}}()) == 1

    vc2 = MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.EqualTo(zero(T)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.EqualTo{T}}()) == 1

    cf = MOI.ScalarAffineFunction{T}(MOI.ScalarAffineTerm{T}.([one(T), -one(T)], v), zero(T))
    c = MOI.add_constraint(model, cf, MOI.GreaterThan(zero(T)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}}()) == 1

    objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.0], v), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) in [ MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT ]

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0,0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 0.0
    end

    # Change y fixed value

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(one(T)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.0,1.0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(2)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.0,2.0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(5//2)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.5 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.5,2.5] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 2.5
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(3)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [3.0,3.0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 3.0
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(5)))

     if config.solve

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE ||
            MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE_OR_UNBOUNDED
    end
end

function semiinttest(model::MOI.ModelLike, config::MOIT.TestConfig{T}) where T
    atol = config.atol
    rtol = config.rtol

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.MathOptInterface.Semiinteger{T})

    # 2 variables
    # min  x
    # st   x >= y
    #      x ∈ {0.0} U {2.0} U {3.0}
    #      y = 0.0

    MOI.empty!(model)
    @test MOI.is_empty(model)

    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    vc1 = MOI.add_constraint(model, MOI.SingleVariable(v[1]), MOI.Semiinteger(T(2), T(3)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Semiinteger{T}}()) == 1

    vc2 = MOI.add_constraint(model, MOI.SingleVariable(v[2]), MOI.EqualTo(zero(T)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable,MOI.EqualTo{T}}()) == 1

    cf = MOI.ScalarAffineFunction{T}(MOI.ScalarAffineTerm{T}.([one(T), -one(T)], v), zero(T))
    c = MOI.add_constraint(model, cf, MOI.GreaterThan(zero(T)))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},MOI.GreaterThan{T}}()) == 1

    objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.0], v), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) in [ MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT ]

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0,0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 0.0
    end

    # Change y fixed value

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(one(T)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.0,1.0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 1.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(2)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [2.0,2.0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(5//2)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [3.0,2.5] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.5 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 3.0
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(3)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.ResultCount()) >= 1

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [3.0,3.0] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 0.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ObjectiveBound()) <= 3.0
    end

    MOI.set(model, MOI.ConstraintSet(), vc2, MOI.EqualTo(T(4)))

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE ||
            MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE_OR_UNBOUNDED
    end
end


text:
"""
Library which implements a FixedSize variant of Dictionaries.
These can be stack allocated and have `O(1)` indexing performance without boundcheck
It implements most parts of the `Base.Dict` interface.
This package is usefull, when you basically want anonymous composite types.
You should be a bit carefull with generating a lot of FixedSizeDict's, since
it will compile a unique set of functions for every field of a Dict.
"""
module FixedSizeDictionaries

    import Base: start
    import Base: next
    import Base: done
    import Base: getindex
    import Base: setindex!
    import Base: get
    import Base: haskey
    import Base: values
    import Base: keys
    import Base: ==
    import Base: length

    include("core.jl")

    export AbstractFixedSizeDict
    export FixedKeyDict
    export FixedKeyValueDict
    export @get

end # module

ext: jl


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Deep learning with Flux")

#

# ## Going deep: Deep neural networks
#
# So far, we've learned that if we want to classify more than two fruits, we'll need to go beyond using a single neuron and use *multiple* neurons to get multiple outputs. We can think of stacking these multiple neurons together in a single neural layer.
#
# Even so, we found that using a single neural layer was not enough to fully distinguish between bananas, grapes, **and** apples. To do this properly, we'll need to add more complexity to our model. We need not just a neural network, but a *deep neural network*.
#
# There is one step remaining to build a deep neural network. We have been saying that a neural network takes in data and then spits out `0` or `1` predictions that together declare what kind of fruit the picture is. However, what if we instead put the output of one neural network layer into another neural network layer?
#
# This gets pictured like this below:
#
#
#
# On the left we have 3 data points in blue. Those 3 data points each get fed into 4 neurons in purple. Each of those 4 neurons produces a single output, but those output are each fed into three neurons (the second layer of purple). Each of those 3 neurons spits out a single value, and those values are fed as inputs into the last layer of 6 neurons. The 6 values that those final neurons produce are the output of the neural network. This is a deep neural network.

#-

# ### Why would a deep neural network be better?
#
# This is a little perplexing when you first see it. We used neurons to train the model before: why would sticking the output from neurons into other neurons help us fit the data better? The answer can be understood by drawing pictures. Geometrically, the matrix multiplication inside of a layer of neurons is streching and rotating the axis that we can vary:

## using Plots, Random; Random.seed!(0); xs = range(0, stop=10, length=200); ys = randn.() .+ xs .+ 0.5.*xs.^1.5 .+ 6.0.*sin.(clamp.(xs .- 5, 0, Inf)) .+ 4.5.*clamp.(.-abs.(xs .- 2), -2, 0); scatter(xs, ys, label="", ticks = false, title="Original data", ylim=extrema(ys).+[-.5,.5]); savefig("data/17-raw.png"); fit1 = [xs ones(size(xs))] \ ys; linear1(x) = fit1[1]*x + fit1[2]; plot!(linear1, label="Linear transform", linewidth=4, legend=:topleft); savefig("data/17-raw-withfit.png");
HTML("""""")

# A nonlinear transformation—such as the sigmoid function or a piecewise linear function—then adds a bump to this linearly-transformed data:

## scatter(xs, ys .- linear1.(xs), ticks = false, label = "", title="data with linear transform", ylim=extrema(ys).+[-.5,.5]); savefig("data/17-linear1.png"); nonlinearity(x) = clamp.(4.0.*x .- 4.5, -4, 4); plot!(nonlinearity, label = "nonlinearity", linewidth=4); savefig("data/17-linear1-withfit.png");
HTML("""""")

# Resulting in a bit more of the data accounted for:

## ys2 = ys .- linear1.(xs) .- nonlinearity.(xs); scatter(xs, ys2, ticks = false, label = "", title="data after first 'layer'", ylim=extrema(ys).+[-.5,.5]); savefig("data/17-nonlinear1.png"); fit2 = [xs ones(size(xs))] \ ys2; linear2(x) = fit2[1]*x + fit2[2]; plot!(linear2, label="second linear tranform", linewidth=4); savefig("data/17-nonlinear1-withfit.png");
HTML("""""")

# Now let's repeat this process. When we send the data through another layer of neurons, we get another rotation and another "bump":

## ys3 = ys2 .- linear2.(xs); scatter(xs .- 4.5, ys3, ticks = false, label = "", title="data after first layer + second linear transform", ylim=extrema(ys).+[-.5,.5]); savefig("data/17-linear2.png"); plot!(nonlinearity, label = "nonlinearity", linewidth=4); savefig("data/17-linear2-withfit.png");
HTML("""""")

#-

## ys4 = ys3 .- nonlinearity.(xs .- 4.5); scatter(xs, ys4, ticks = false, label = "", title="data after second 'layer'", ylim=extrema(ys).+[-.5,.5]); savefig("data/17-nonlinear2.png");
HTML("""""")

#
# Visually, we see that if we keep doing this process we can make the axis line up with any data. What this means is that **if we have enough layers, then our neural network can approximate any model**.
#
# The trade-off is that with more layers we have more parameters, so it may be harder (i.e. computationally intensive) to train the neural network. But we have the guarantee that the model has enough freedom such that there are parameters that will give the correct output.
#
# Because this model is so flexible, the problem is reduced to that of learning: do the same gradient descent method on this much larger model (but more efficiently!) and we can make it classify our data correctly. This is the power of deep learning.

#-

# **Important caveat**: The above is a _cartoon_ that loosely describes the concepts behind a deep neural network. In reality, remember that we're not successively fitting "the best" linear transformation and then subsequently fitting a nonlinearity; the two are happening simultaneously through gradient descent. This allows the linear transform to — in effect — rotate and stretch the data such that the nonlinearity "falls" in the appropriate place to best minimize the loss function.

#-

# # Deep learning with Flux

#-

# Let's load the same datasets from the previous lecture and pre-process them in the same way:

using CSV, DataFrames, Flux, Plots
apples1 = DataFrame(CSV.File("data/Apple_Golden_1.dat", delim='\t', allowmissing=:none, normalizenames=true))
apples2 = DataFrame(CSV.File("data/Apple_Golden_2.dat", delim='\t', allowmissing=:none, normalizenames=true))
apples3 = DataFrame(CSV.File("data/Apple_Golden_3.dat", delim='\t', allowmissing=:none, normalizenames=true))
apples = vcat(apples1, apples2, apples3)
bananas = DataFrame(CSV.File("data/Banana.dat", delim='\t', allowmissing=:none, normalizenames=true))
grapes1 = DataFrame(CSV.File("data/Grape_White.dat", delim='\t', allowmissing=:none, normalizenames=true))
grapes2 = DataFrame(CSV.File("data/Grape_White_2.dat", delim='\t', allowmissing=:none, normalizenames=true))
grapes = vcat(grapes1, grapes2)
## Extract out the features and construct the corresponding labels
x_apples  = [ [apples[i, :red], apples[i, :blue]] for i in 1:size(apples, 1) ]
x_bananas  = [ [bananas[i, :red], bananas[i, :blue]] for i in 1:size(bananas, 1) ]
x_grapes = [ [grapes[i, :red], grapes[i, :blue]] for i in 1:size(grapes, 1) ]
xs = vcat(x_apples, x_bananas, x_grapes)
ys = vcat(fill(Flux.onehot(1, 1:3), size(x_apples)),
          fill(Flux.onehot(2, 1:3), size(x_bananas)),
          fill(Flux.onehot(3, 1:3), size(x_grapes)));

# In the previous lecture, we used a `Dense(2, 3, σ)` as our model. Now we want to construct multiple layers and chain them together:

layer1 = Dense(2, 4, σ)
layer2 = Dense(4, 3, σ)

#-

layer2(layer1(xs[1]))

#-

#nb ?Chain
#jl @doc Chain

#-

m = Chain(layer1, layer2)
m(xs[1])

#-

xs[1] |> layer1 |> layer2

# ### The core algorithm from the last lecture

model = Chain(Dense(2, 3, σ)) # Update this!
L(x,y) = Flux.mse(model(x), y)
opt = SGD(params(model))
Flux.train!(L, zip(xs, ys), opt)

#-

## Recall his is how we repeatedly walked down our gradient previously...
for _ in 1:1000
    Flux.train!(L, zip(xs, ys), opt)
end
## But our model is now more complicated and this will take more time!

#-

data = zip(xs, ys)
@time Flux.train!(L, data, opt)
@time Flux.train!(L, data, opt)

# ### Improving efficiency by batching

length(data)

#-

first(data)

# Recall our matrix-vector multiplication from the previous lecture:

W = [10 1;
     20 2;
     30 3]
x = [3;
     2]
W*x

#-

Flux.batch(xs)

#-

model(Flux.batch(xs))

#-

databatch = (Flux.batch(xs), Flux.batch(ys))
@time Flux.train!(L, (databatch,), opt)
@time Flux.train!(L, (databatch,), opt)

#-

Flux.train!(L, Iterators.repeated(databatch, 10000), opt)

#-

L(databatch[1], databatch[2])

# ### Visualization

using Plots
function plot_decision_boundaries(model, x_apples, x_bananas, x_grapes)
    plot()

    contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[1], levels=[0.5, 0.501], color = cgrad([:blue, :blue]), colorbar=:none)
    contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[2], levels=[0.5,0.501], color = cgrad([:green, :green]), colorbar=:none)
    contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[3], levels=[0.5,0.501], color = cgrad([:red, :red]), colorbar=:none)

    scatter!(first.(x_apples), last.(x_apples), m=:cross, label="apples", color = :blue)
    scatter!(first.(x_bananas), last.(x_bananas), m=:circle, label="bananas", color = :green)
    scatter!(first.(x_grapes), last.(x_grapes), m=:square, label="grapes", color = :red)
end
plot_decision_boundaries(model, x_apples, x_bananas, x_grapes)

# ### Further improvements with a better loss function and normalization of outputs

scatter([0],[0], label="correct answer", xlabel="model output: [1-x,x]", ylabel="loss against [1, 0]", legend=:topleft, title="Loss function behavior")
plot!(x->Flux.mse([1-x, x/2], [1,0]), -1.5, 1.5, label="mse")
## plot!(x->Flux.crossentropy([1-x, x/2], [1,0]), 0, 1, label="crossentropy")

#-

sum(model(xs[1]))

#-

Flux.mse([0.01,0.98,0.01], [1.0,0,0])

#-

softmax([1.0,-3,0])

# ### The new algorithm
#
# Use `softmax` as a final normalization and change the loss function to `crossentropy`:

model = Chain(Dense(2, 4, σ), Dense(4, 3, identity), softmax)
L(x,y) = Flux.crossentropy(model(x), y)
opt = SGD(params(model))

#-

Flux.train!(L, Iterators.repeated(databatch,5000), opt)

#-

plot_decision_boundaries(model, x_apples, x_bananas, x_grapes)

max_stars_repo_path: Courses/Deep learning with Flux/1700.Introduction-to-deep-learning.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
x0 = randn(10)
it, x, sol = ProximalAlgorithms.Template(x0)
@test it == 10
@test norm(x - x0) == 0.0
println(sol)


text:
@testset "Raw solver: internals" begin
    models = [
        sudoku(2; modeler=:raw),
    ]

    for m in models
        @info describe(m)
        s = solver(m; options=Options(print_level=:verbose, time_limit = Inf, iteration=Inf, info_path="info.json"))
        for x in keys(get_variables(s))
            @test get_name(s, x) == "x$x"
            for c in get_cons_from_var(s, x)
                @test x ∈ get_vars_from_cons(s, c)
            end
            @test constriction(s, x) == 3
            @test draw(s, x) ∈ get_domain(s, x)
        end

        for c in keys(get_constraints(s))
            @test length_cons(s, c) == 4
        end

        for x in keys(get_variables(s))
            add_var_to_cons!(s, 3, x)
            delete_var_from_cons!(s, 3, x)
            @test length_var(s, x) == 4
        end

        for c in keys(get_constraints(s))
            add_var_to_cons!(s, c, 17)
            @test length_cons(s, c) == 5
            @test 17 ∈ get_constraint(s, c)
            delete_var_from_cons!(s, c, 17)
            @test length_cons(s, c) == 4
        end
        solve!(s)
        solution(s)

        # TODO: temp patch for coverage, make it nice
        for x in keys(LocalSearchSolvers.tabu_list(s))
            LocalSearchSolvers.tabu_value(s, x)
        end
        # LocalSearchSolvers._values!(s, Dictionary{Int,Int}())

        # display(solution(s))
        @info time_info(s)
        rm("info.json")
    end
end

@testset "Raw solver: sudoku" begin
    sudoku_instance = collect(Iterators.flatten([
        9  3  0  0  0  0  0  4  0
        0  0  0  0  4  2  0  9  0
        8  0  0  1  9  6  7  0  0
        0  0  0  4  7  0  0  0  0
        0  2  0  0  0  0  0  6  0
        0  0  0  0  2  3  0  0  0
        0  0  8  5  3  1  0  0  2
        0  9  0  2  8  0  0  0  0
        0  7  0  0  0  0  0  5  3
    ]))

    s = solver(sudoku(3; start = sudoku_instance, modeler = :raw); options = Options(print_level = :minimal, iteration = Inf, time_limit = 10))
    display(Dictionary(1:length(sudoku_instance), sudoku_instance))
    solve!(s)
    display(solution(s))
end

@testset "Raw solver: golomb" begin
    s = solver(golomb(5, modeler = :raw); options = Options(print_level = :minimal, iteration = 1000))
    solve!(s)

    @info "Results golomb!"
    @info "Values: $(get_values(s))"
    @info "Sol (val): $(best_value(s))"
    @info "Sol (vals): $(!isnothing(best_value(s)) ? best_values(s) : nothing)"
end

@testset "Raw solver: mincut" begin
    graph = zeros(5, 5)
    graph[1,2] = 1.0
    graph[1,3] = 2.0
    graph[1,4] = 3.0
    graph[2,5] = 1.0
    graph[3,5] = 2.0
    graph[4,5] = 3.0
    s = solver(mincut(graph, source=1, sink=5), options =  Options(print_level = :minimal))
    solve!(s)
    @info "Results mincut!"
    @info "Values: $(get_values(s))"
    @info "Sol (val): $(best_value(s))"
    @info "Sol (vals): $(!isnothing(best_value(s)) ? best_values(s) : nothing)"

    s = solver(mincut(graph, source=1, sink=5, interdiction=1), options =  Options(print_level = :minimal))
    solve!(s)
    @info "Results 1-mincut!"
    @info "Values: $(get_values(s))"
    @info "Sol (val): $(best_value(s))"
    @info "Sol (vals): $(!isnothing(best_value(s)) ? best_values(s) : nothing)"

    s = solver(mincut(graph, source=1, sink=5, interdiction=2); options = Options(print_level=:minimal, time_limit = 15, iteration=Inf))
    @info describe(s)
    solve!(s)
    @info "Results 2-mincut!"
    @info "Values: $(get_values(s))"
    @info "Sol (val): $(best_value(s))"
    @info "Sol (vals): $(!isnothing(best_value(s)) ? best_values(s) : nothing)"
    @info time_info(s)
end


text:
module Kalman

using Distributions, GaussianDistributions
using Trajectories, DynamicIterators
using Random, LinearAlgebra

import DynamicIterators: evolve, dyniterate, @returnnothing
import Base: iterate, IteratorSize, IteratorEltype, eltype, length

import Random.rand
import GaussianDistributions: ⊕

meancov(G) = mean(G), cov(G)
meancov(G::Tuple) = G

struct Condition{T} <: Message2
    u::T
    ll::Float64
end
struct Filter{T} <: Message2
    u::T
    ll::Float64
end

macro NT(args...)
    :(NamedTuple{($(args)...,)})
end

include("dyniterate.jl")

export LinearEvolution, Evolution
include("evolution.jl")

export LinearObservation, LinearObservationModel, GenericLinearObservationModel, Observation, Observe
include("observation.jl")

export LinearStateSpaceModel, StateObs
include("statespacemodel.jl")

export kalmanfilter
include("filter.jl")

export rts_smoother
include("smoother.jl")
include("backwardsampler.jl")

include("combinator.jl")

end # module


text:
"""
Abstract supertype of properties that can be checked.

Every concrete subtype should provide the following functions:
  - `dim(𝑃::Property)::Int`
  - `check(𝑃::Property, X::LazySet; witness::Bool=false)`
"""
abstract type Property end


text:
#!/usr/bin/env julia

"""
# Problem 4: Largest palindrome product

The prime factors of 13195 are 5, 7, 13 and 29.

A palindromic number reads the same both ways. The largest palindrome made from
the product of two 2-digit numbers is 9009 = 91 × 99.

Find the largest palindrome made from the product of two `d`-digit numbers.
"""
module Problem004

include("Digits.jl")

"""
    solve(; d::Integer = 3)

Returns the solution for [`Problem004`](@ref) with the given parameters.

## Preconditions
- `d > 0`
"""
function solve(; d::Integer = 3)
    # calculate max and min d-digit numbers
    base = convert(typeof(d), 10)
    max_factor = base^d - oneunit(d)
    min_factor = base^(d - oneunit(d))

    # search for palindromic products of d-digit numbers
    max_product = -oneunit(d)
    for n in max_factor:-oneunit(d):min_factor
        for m in n:-oneunit(d):min_factor
            product = n * m

            # stop checking products with n when we drop below max_product
            if max_product ≥ product
                break
            end

            # if product is a palindrome, it's the max possible for this n
            if Digits.is_palindrome(product)
                max_product = product
                break
            end
        end
    end

    return max_product
end

end

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
@testset "encoder" begin
    @testset "Real image encode test" begin
        @info "Check if all real images are visually good"
        img_gray = testimage("mandril_gray")
        img_color = testimage("mandril_color")
        img3d = testimage("mri-stack")
        test_sixel_display() do
            sixel_encode(img_gray)
            sixel_encode(img_color)
            sixel_encode(img3d)
            # higher dimensional array is stacked as if it is 3d
            sixel_encode(reshape(img3d, size(img3d)[1:2]..., 3, :))
        end
    end

    @testset "1d vector" begin
        for img in (
            repeat(Gray.(0:0.1:0.9), inner=10),
            repeat(distinguishable_colors(10), inner=(10, ))
        )
            @info "1d vector: $(eltype(img))"
            @test size(img) == (100, )

            enc = Sixel.default_encoder(img)
            @test enc isa Sixel.LibSixel.LibSixelEncoder
            test_sixel_display() do
                sixel_encode(img, enc)
                println()
            end

            io = IOBuffer()
            sixel_encode(io, img, enc)
            bufferdata = String(take!(io))
            h, w = size(img, 1), 1
            @test startswith(bufferdata, "\ePq\"1;1;$w;$h") && endswith(bufferdata, "\e\\")
            test_sixel_display() do
                println(bufferdata)
            end

            tmp_file = tempname()
            open(tmp_file, "w") do io
                sixel_encode(io, img, enc)
            end
            filedata = read(tmp_file, String)
            @test startswith(filedata, "\ePq\"1;1;$w;$h") && endswith(filedata, "\e\\")
            test_sixel_display() do
                println(filedata)
            end

            @info "transpose check: the first two should look the same, while the third is transposed."
            test_sixel_display() do
                sixel_encode(img, enc)
                sixel_encode(img, enc; transpose=false)
                sixel_encode(img, enc; transpose=true)
            end
        end
    end

    @testset "2d matrix" begin
        for img in (
            repeat(Gray.(0:0.1:0.9), inner=(10, 50)),
            repeat(distinguishable_colors(10), inner=(10, 50))
        )
            @info "2d matrix: $(eltype(img))"
            h, w = size(img)
            @test size(img) == (100, 50)

            enc = Sixel.default_encoder(img)
            @test enc isa Sixel.LibSixel.LibSixelEncoder
            test_sixel_display() do
                sixel_encode(img, enc)
                println()
            end
            io = IOBuffer()
            sixel_encode(io, img, enc)
            bufferdata = String(take!(io))
            @test startswith(bufferdata, "\ePq\"1;1;$w;$h") && endswith(bufferdata, "\e\\")
            test_sixel_display() do
                println(bufferdata)
            end

            tmp_file = tempname()
            open(tmp_file, "w") do io
                sixel_encode(io, img, enc)
            end
            filedata = read(tmp_file, String)
            @test startswith(filedata, "\ePq\"1;1;$w;$h") && endswith(filedata, "\e\\")
            test_sixel_display() do
                println(filedata)
            end

            @info "transpose check: the first two should look the same, while the third is transposed."
            test_sixel_display() do
                sixel_encode(img, enc)
                sixel_encode(img, enc; transpose=false)
                sixel_encode(img, enc; transpose=true)
            end
        end
    end

    @testset "3d array" begin
        for img in (
            repeat(Gray.(0:0.1:0.9), inner=(10, 50, 3)),
            repeat(distinguishable_colors(5), inner=(20, 50, 3))
        )
            @info "3d array: $(eltype(img))"
            @test size(img) == (100, 50, 3)

            enc = Sixel.default_encoder(img)
            @test enc isa Sixel.LibSixel.LibSixelEncoder
            test_sixel_display() do
                sixel_encode(img, enc)
                println()
            end
            io = IOBuffer()
            sixel_encode(io, img, enc)
            bufferdata = String(take!(io))
            h, w, c = size(img)
            w = w * c
            @test startswith(bufferdata, "\ePq\"1;1;$w;$h") && endswith(bufferdata, "\e\\")
            test_sixel_display() do
                println(bufferdata)
            end

            tmp_file = tempname()
            open(tmp_file, "w") do io
                sixel_encode(io, img, enc)
            end
            filedata = read(tmp_file, String)
            h, w, c = size(img)
            w = w * c
            @test startswith(filedata, "\ePq\"1;1;$w;$h") && endswith(filedata, "\e\\")
            test_sixel_display() do
                println(filedata)
            end

            @info "transpose check: the first two should look the same, while the third is transposed."
            test_sixel_display() do
                sixel_encode(img, enc)
                sixel_encode(img, enc; transpose=false)
                sixel_encode(img, enc; transpose=true)
            end

            io = IOBuffer()
            sixel_encode(io, img, enc; transpose=true)
            bufferdata = String(take!(io))
            h, w, c = size(img)
            h = h * c
            h, w = w, h # transpose
            @test startswith(bufferdata, "\ePq\"1;1;$w;$h") && endswith(bufferdata, "\e\\")
        end
    end

    @testset "various eltype" begin
        enc = Sixel.LibSixel.LibSixelEncoder()
        img = repeat(distinguishable_colors(5), inner=(20, 50))

        # test different color spaces
        @info "The images below should visually look the same"
        for T in (
            RGB{N0f8}, RGB{Float32}, RGB24,
            HSV, Lab, ARGB, RGBA, BGR, BGRA
        )
            test_sixel_display() do
                sixel_encode(T.(img), enc)
            end
        end

        # test different storage format
        @info "The images below should visually look the same"
        img = gray.(Gray.(img))
        for T in (Gray, N0f8, Float32, Float64, x->round(Int, 255x))
            test_sixel_display() do
                sixel_encode(T.(img), enc)
            end
        end
    end

    @testset "untypical array types" begin
        # lazy array that does not occupy full memory
        img = Diagonal(repeat(distinguishable_colors(5), inner=20))
        w, h = size(img)
        io = IOBuffer()
        sixel_encode(io, img)
        bufferdata = String(take!(io))
        @test startswith(bufferdata, "\ePq\"1;1;$w;$h") && endswith(bufferdata, "\e\\")
        @info "Test Diagonal array"
        test_sixel_display() do
            println(bufferdata)
        end

        # array that has un-contiguous memory layout
        ori_img = repeat(distinguishable_colors(5), inner=(20, 50))
        img = PermutedDimsArray(permutedims(ori_img, (2, 1)), (2, 1))
        @test size(img) == (100, 50) == size(ori_img)
        @test stride(img, 1) == 50
        @test stride(ori_img, 1) == 1

        io = IOBuffer()
        sixel_encode(io, ori_img)
        bufferdata = String(take!(io))
        h, w = size(ori_img)
        @test startswith(bufferdata, "\ePq\"1;1;$w;$h") && endswith(bufferdata, "\e\\")
        io = IOBuffer()
        sixel_encode(io, img)
        @test bufferdata == String(take!(io))

        io = IOBuffer()
        sixel_encode(io, ori_img, transpose=false)
        bufferdata = String(take!(io))
        h, w = size(ori_img)
        @test startswith(bufferdata, "\ePq\"1;1;$w;$h") && endswith(bufferdata, "\e\\")
        io = IOBuffer()
        sixel_encode(io, img, transpose=false)
        @test bufferdata == String(take!(io))
    end

    @testset "OffsetArray" begin
        enc = Sixel.LibSixelEncoder()
        for img in [
            repeat(distinguishable_colors(10), inner=(10, )),
            repeat(distinguishable_colors(10), inner=(10, 50)),
            repeat(distinguishable_colors(5), inner=(20, 50, 3))
        ]
            io = IOBuffer()
            sixel_encode(io, img, enc)
            ref = String(take!(io))

            imgo = OffsetArray(img, OffsetArrays.Origin(0))
            io = IOBuffer()
            sixel_encode(io, imgo, enc)
            actual = String(take!(io))
            @test ref == actual
        end
    end
end

@testset "decoder" begin
    @testset "API" begin
        dec = Sixel.default_decoder()
        @test dec isa Sixel.LibSixelDecoder

        tmp_file = tempname()
        img = repeat(Gray.(0:0.1:0.9), inner=(10, 50))
        sz = size(img)
        open(tmp_file, "w") do io
            sixel_encode(io, img)
        end

        # IO
        img_readback = open(tmp_file, "r") do io
            sixel_decode(io)
        end
        @test sz == size(img_readback)
        @test img_readback isa IndirectArray
        @test eltype(img_readback) == RGB{N0f8}

        img_readback = open(tmp_file, "r") do io
            sixel_decode(io, dec)
        end
        @test sz == size(img_readback)
        @test img_readback isa IndirectArray
        @test eltype(img_readback) == RGB{N0f8}

        img_readback = open(tmp_file, "r") do io
            sixel_decode(Gray{N0f8}, io)
        end
        @test sz == size(img_readback)
        @test img_readback isa IndirectArray
        @test eltype(img_readback) == Gray{N0f8}

        # String and Vector{UInt8}
        bytes_data = open(tmp_file, "r") do io
            read(io)
        end
        string_data = read(tmp_file, String)
        for src in (bytes_data, string_data)
            img_readback = sixel_decode(src, dec)
            @test sz == size(img_readback)
            @test img_readback isa IndirectArray
            @test eltype(img_readback) == RGB{N0f8}

            img_readback = sixel_decode(Gray{N0f8}, src, dec)
            @test sz == size(img_readback)
            @test img_readback isa IndirectArray
            @test eltype(img_readback) == Gray{N0f8}
        end
    end

    @testset "Quality test" begin
        enc = Sixel.LibSixelEncoder()
        dec = Sixel.LibSixelDecoder()
        tmp_file = tempname()

        for img in (
            repeat(Gray.(0:0.1:0.9), inner=10),
            repeat(distinguishable_colors(10), inner=(10, )),
            repeat(Gray.(0:0.1:0.9), inner=(10, 50)),
            repeat(distinguishable_colors(10), inner=(10, 50)),
            repeat(Gray.(0:0.1:0.9), inner=(10, 50, 3)),
            repeat(distinguishable_colors(5), inner=(20, 50, 3))
        )
            open(tmp_file, "w") do io
                sixel_encode(io, img, enc)
            end

            img_readback = open(tmp_file, "r") do io
                sixel_decode(eltype(img), io, dec)
            end

            img_readback = reshape(img_readback, size(img))
            # 30 is actually pretty good given that sixel encode always do quantization
            @test assess_psnr(img, img_readback) > 30
        end
    end
end


text:
using DocOpt

doc="""Voice conversion based on WORLD-based speech analysis and
synthesis framework

Usage:
    vc.jl [options] <input_wav> <model_jld> <dst_wav>
    vc.jl --version
    vc.jl -h | --help

Options:
    -h --help        show this message
    --period=PERIOD  frame period in msec [default: 5.0]
    --order=ORDER    order of mel cepsrum [default: 25]
    --alpha=ALPHA    all-pass constant [default: 0.0]
    --gv=MODEL       global variance [default: ]
    --trajectory     trajectory-based parameter conversion
    --T=t            maximum length for one-step trajectory conversion [default: 100]
"""

using VoiceConversion
using MelGeneralizedCepstrums
using WAV
using HDF5, JLD
using WORLD

using Logging
@Logging.configure(level=DEBUG, output=STDOUT)

function main()
    args = docopt(doc, version=v"0.0.1")

    x, fs = wavread(args["<input_wav>"])
    @assert size(x, 2) == 1 "The input data must be monoral."
    x = map(Float64, vec(x))
    fs = Int(fs)
    @info("length of input signal is $(length(x)/fs) sec.")

    period = parse(Float64, args["--period"])
    order = parse(Int, args["--order"])
    α = parse(Float64, args["--alpha"])
    if α == 0.0
        α = mcepalpha(fs)
    end
    trajectory = args["--trajectory"]
    gvmodel = string(args["--gv"])

    # Load mapping model
    gmm = load(args["<model_jld>"])
    if gmm["diff"]
        @warn("The model seem to be trained on differencial features")
    end

    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])
    if trajectory
        mapper = TrajectoryGMMMap(mapper, parse(Int, args["--T"]))
    end

    if trajectory && !isempty(gvmodel)
        gv = load(gvmodel)
        gv["n_components"] == 1 || error("only single Gaussian for GV is supported")
        μᵛ = gv["means"][:,1]
        Σᵛᵛ = gv["covars"][:,:,1]
        mapper = TrajectoryGVGMMMap(mapper, μᵛ, Σᵛᵛ)
    end

    elapsed_fe = @elapsed begin
        f0, timeaxis = dio(x, fs, DioOption(period=period))
        f0 = stonemask(x, fs, timeaxis, f0)
        spectrogram = cheaptrick(x, fs, timeaxis, f0)
        src = sp2mc(spectrogram, order, α)
        ap = d4c(x, fs, timeaxis, f0)
    end

    @info("elapsed time in feature extraction is $(elapsed_fe) sec.")
    if trajectory
        # add delta feature
        src = [reshape(src[1,:], 1, size(src,2)); push_delta(src[2:end,:])]
    end

    # Perform conversion
    elapsed_vc = @elapsed converted = vc(mapper, src)
    @info("elapsed time in conversion process is $(elapsed_vc) sec.")

    # Mel-Cepstrum to spectrum
    fftlen = size(spectrogram,1)*2-1
    converted_spectrogram = mc2sp(converted, α, fftlen)

    # Waveform synthesis using WORLD
    elapsed_syn = @elapsed begin
        y = synthesis(f0, converted_spectrogram, ap, period, fs, length(x))
    end
    @info("elapsed time in waveform synthesis is $(elapsed_syn) sec.")

    wavwrite(map(Float64, y), args["<dst_wav>"], Fs=fs)
    @info("Dumped to ", args["<dst_wav>"])
end

main()


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# miscellaneous initialisiations for graph
function initGraph!(graph::Graph)
	# should already have nodes, arcs
	numNodes = length(graph.nodes)
	@assert(numNodes > 0 && length(graph.arcs) > 0)

	graph.arcDists = [arc.distance for arc in graph.arcs]

	graph.light = LightGraphs.DiGraph(numNodes) # create LightGraph version of graph
	for arc in graph.arcs
		LightGraphs.add_edge!(graph.light, arc.fromNodeIndex, arc.toNodeIndex) # add edges to light graph
	end
	graph.fadjList = graph.light.fadjlist
	graph.badjList = graph.light.badjlist

	if !graph.isReduced
		graph.nodePairArcIndex = spzeros(Int, numNodes, numNodes) # sparse, to save on memory
		for arc in graph.arcs
			@assert(graph.nodePairArcIndex[arc.fromNodeIndex, arc.toNodeIndex] == 0) # check that there is only one arc for this ordered node pair
			graph.nodePairArcIndex[arc.fromNodeIndex, arc.toNodeIndex] = arc.index
		end
	end
end

# find nearest node by evaluating distances to all nodes from given location
# should not be used except for asserting results of faster algorithms
# can set offRoadAccessRequired in order to search nodes with off-road access, or all nodes
function findNearestNode(map::Map, nodes::Vector{Node}, location::Location;
	offRoadAccessRequired::Bool = true)

	numNodes = length(nodes)

	# compare square distances
	bestSqrDist = Inf
	chosenNode = nullIndex
	sqrDist = 0.0 # init
	for i = 1:numNodes
		node = nodes[i]
		if offRoadAccessRequired && !node.offRoadAccess
			continue
		end
		sqrDist = squareDist(map, node.location, location)
		if sqrDist < bestSqrDist
			bestSqrDist = sqrDist
			chosenNode = i
		end
	end
	dist = sqrt(bestSqrDist)

	return chosenNode, dist
end

# Set arc distance with normDist function,
# for arcs that pass arcFilter.
function setArcDistances!(graph::Graph, map::Map;
	arcFilter::Function = (arc->isnan(arc.distance)))
	@assert(!graph.isReduced)
	for arc in filter(arcFilter, graph.arcs)
		arc.distance = normDist(map, graph.nodes[arc.fromNodeIndex].location, graph.nodes[arc.toNodeIndex].location)
	end
end

# for graph (set of nodes, arcs), check that:
# nodes are inside map borders,
# nodes are attached to at least one arc each,
# arcs are connected to two nodes,
# arcs have non-negative travel time,
# graph is strongly connected
function checkGraph(graph::Graph, map::Map)

	# shorthand
	@unpack nodes, arcs = graph
	numNodes = length(nodes)
	numArcs = length(arcs)

	# check node and arc indices
	for i = 1:numNodes
		@assert(nodes[i].index == i)
	end
	for i = 1:numArcs
		@assert(arcs[i].index == i)
	end

	# check arc distances
	for i = 1:numArcs
		@assert(!isnan(arcs[i].distance), "arc $i has distance NaN")
		@assert(arcs[i].distance >= 0, "arc $i has negative distance")
	end

	# check that nodes are inside map borders
	for i = 1:numNodes
		@assert(map.xMin <= nodes[i].location.x, "node $i outside border")
		@assert(map.xMax >= nodes[i].location.x, "node $i outside border")
		@assert(map.yMin <= nodes[i].location.y, "node $i outside border")
		@assert(map.yMax >= nodes[i].location.y, "node $i outside border")
	end

	# check that all nodes are attached to at least one arc
	# and check that arcs are all attached to two nodes
	nodeUsed = Vector{Bool}(undef, numNodes)
	nodeUsed[1:end] .= false
	for i = 1:numArcs
		@assert(1 <= arcs[i].fromNodeIndex && arcs[i].fromNodeIndex <= numNodes, "arc $i fromNodeIndex is outside range")
		@assert(1 <= arcs[i].toNodeIndex && arcs[i].toNodeIndex <= numNodes, "arc $i toNodeIndex is outside range")
		@assert(arcs[i].fromNodeIndex != arcs[i].toNodeIndex, "arc $i fromNodeIndex and toNodeIndex are same")
		nodeUsed[arcs[i].fromNodeIndex] = true
		nodeUsed[arcs[i].toNodeIndex] = true
	end
	for i = 1:numNodes
		@assert(nodeUsed[i], "node $i not used")
	end

	@assert(LightGraphs.is_weakly_connected(graph.light), "graph is not weakly connected")
	@assert(LightGraphs.is_strongly_connected(graph.light), "graph is not strongly connected")
end


text:
function CR3BPSystem(arg0::CelestialBody, arg1::CelestialBody, arg2::jdouble)
    return CR3BPSystem((CelestialBody, CelestialBody, jdouble), arg0, arg1, arg2)
end

function CR3BPSystem(arg0::CelestialBody, arg1::CelestialBody, arg2::jdouble, arg3::jdouble)
    return CR3BPSystem((CelestialBody, CelestialBody, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function get_ddim(obj::CR3BPSystem)
    return jcall(obj, "getDdim", jdouble, ())
end

function get_gamma(obj::CR3BPSystem, arg0::LagrangianPoints)
    return jcall(obj, "getGamma", jdouble, (LagrangianPoints,), arg0)
end

function get_l_position(obj::CR3BPSystem, arg0::LagrangianPoints)
    return jcall(obj, "getLPosition", Vector3D, (LagrangianPoints,), arg0)
end

function get_mass_ratio(obj::CR3BPSystem)
    return jcall(obj, "getMassRatio", jdouble, ())
end

function get_name(obj::CR3BPSystem)
    return jcall(obj, "getName", JString, ())
end

function get_primary(obj::CR3BPSystem)
    return jcall(obj, "getPrimary", CelestialBody, ())
end

function get_real_apv(obj::CR3BPSystem, arg0::AbsolutePVCoordinates, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getRealAPV", AbsolutePVCoordinates, (AbsolutePVCoordinates, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_rotating_frame(obj::CR3BPSystem)
    return jcall(obj, "getRotatingFrame", Frame, ())
end

function get_secondary(obj::CR3BPSystem)
    return jcall(obj, "getSecondary", CelestialBody, ())
end

function get_tdim(obj::CR3BPSystem)
    return jcall(obj, "getTdim", jdouble, ())
end

function get_vdim(obj::CR3BPSystem)
    return jcall(obj, "getVdim", jdouble, ())
end


text:
function generate(n_samples)
    mixing = [0.2, 0.5, 0.3]
    model = MixtureModel(map(θ -> Normal(θ, 1.0), [-4.5, 0.0, 4.]), mixing)
    return rand(model, n_samples)
end

mutable struct SS
    μ::Float64
    count::Int
end

struct HyperParameters
    α::Float64
    λ1::Float64
    λ2::Float64
end

function sufficient_statistics(data, assignment)
    @assert length(data) == length(assignment)
    cluster_ids = unique(assignment)
    results = Dict{Int,SS}()
    for cluster_id in cluster_ids
        xs = data[assignment .== cluster_id]
        μ = mean(xs)
        count = length(xs)
        results[cluster_id] = SS(μ, count)
    end
    return results
end

mutable struct GibbsState
    data::Vector{Float64}
    assignment::Vector{Int}
    cluster_ids::Vector{Int}
    π::Vector{Float64}
    θ::Vector{Float64}
    v::Float64
    hp::HyperParameters
    suffstats::Dict{Int,SS}

    function GibbsState(data, n_clusters; α=1)
        cluster_ids = Vector(1:n_clusters)
        assignment = [rand(cluster_ids) for _ in 1:length(data)]
        π = ones(n_clusters) / n_clusters
        θ = randn(n_clusters)
        v = 1.
        hp = HyperParameters(α, 0.0, 1.0)
        suffstats = sufficient_statistics(data, assignment)
        state = new(data, assignment, cluster_ids, π, θ, v, hp, suffstats)
        return state
    end
end


text:
# An example of the Burgers code for the
# propagation of Burgers equation in time
# Creates a colormap showing u(x, t) of Burgers equation
# initialized with the sine function used in lecture

using Burgers
using Printf
using PyPlot
include("Utils.jl")
using .Utils

Ttot   = 2.0   # total duration
Nt     = 1024  # number of time divisions
Ncells = 256   # number of position divisions
Nghost = 2     # padding on position divisions for boundary conditions

res   = Resolution{Float64, Int64}(Ttot, Nt, Ncells, Nghost)
sc    = minmod::SlopeCalc
u2fun = get_evolved(get_sinewave, res, sc)

# u(x, t) can be plotted as a colour map
xlabel("t")
ylabel("x")
pcolormesh(u2fun.tdat, u2fun.xdat, u2fun.u2dat)
colorbar()
savefig("images/colormap2.png", dpi = 300)


text:
function fulltest(m, num_stages, objval, solval, ws, wsσ, testniter, solver)
    isclp(solver) && return
    for K in [-1, 40]
        for maxncuts in [-1, 7]
            for newcut in [:InvalidateSolver]#[:AddImmediately, :InvalidateSolver]
                for cutmode in [StructProg.MultiCutGenerator(), StructProg.AvgCutGenerator()]
                    for detectlb in [false, true]
                        for pruningalgo in [AvgCutPruningAlgo(maxncuts), DecayCutPruningAlgo(maxncuts), DeMatosPruningAlgo(maxncuts)]
                            isclp(solver) && K == -1 && maxncuts == 7 && cutmode == :MultiCut && !detectlb && isa(pruningalgo, DeMatosPruningAlgo) && continue
                            sp = SOI.stochasticprogram(m, num_stages, solver, pruningalgo, cutmode, detectlb, newcut)

                            μ, σ = SOI.optimize!(sp, WaitAndSee.Algorithm(solver, K))
                            @test abs(μ - ws) / ws < (K == -1 ? 1e-6 : .03)
                            @test abs(σ - wsσ) / wsσ <= (K == -1 ? 1e-6 : 1.)

                            if K == -1
                                stopcrit = SOI.CutLimit(0)
                            else
                                stopcrit = SOI.Pereira()
                            end
                            stopcrit |= SOI.IterLimit(42)
                            algo = SDDP.Algorithm(K = K, forwardcuts = true, backwardcuts = false)
                            sol = SOI.optimize!(sp, algo, stopcrit, 0)
                            testniter(SOI.niterations(sol), K, maxncuts, cutmode, detectlb)
                            @test SOI.last_result(sol).status == :Optimal
                            lb = SOI.last_result(sol).lowerbound
                            @test abs(lb - objval) / objval < (K == -1 ? 1e-6 : .03)
                            # FIXME the info is missing the SOI result
                            #v11value = sol.sol[1:4]
                            #@test norm(v11value - solval) / norm(solval) < (K == -1 ? 1e-6 : .3)

                            μ, σ = SOI.optimize!(sp, WaitAndSee.Algorithm(solver, K))
                            @test abs(μ - ws) / ws < (K == -1 ? 1e-6 : .03)
                            @test abs(σ - wsσ) / wsσ <= (K == -1 ? 1e-6 : 1.)

                            StructProg.clear(m)
                        end
                    end
                end
            end
        end
    end
end

@testset "Problem 5.2" begin
    include("prob5.2_2stages.jl")
    include("prob5.2_3stages.jl")
    include("prob5.2_3stages_serial.jl")
end


text:
# # Stabilization of nonlinear systems

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Systems and Control/stabilization_of_nonlinear_systems.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Systems and Control/stabilization_of_nonlinear_systems.ipynb)
# **Adapted from**: Examples 1, 2 and 3 of [PPA04]
#
# [PPA04] Prajna, Stephen, Pablo A. Parrilo, and Anders Rantzer.
# *Nonlinear control synthesis by convex optimization*.
# IEEE Transactions on Automatic Control 49.2 (2004): 310-314.

using Test #src
using DynamicPolynomials
@polyvar x[1:2]

using SumOfSquares
using CSDP
using LinearAlgebra # for ⋅
using MultivariatePolynomials
divergence(f) = sum(differentiate.(f, x))
function controller(f, g, b, α, degs)
    solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
    model = SOSModel(solver)
    a = 1
    monos = monomials(x, degs)
    N = length(monos) + 1
    @variable(model, c[1:N] in MOI.NormOneCone(N))
    c_poly = polynomial(c[2:end], monos)
    fagc = f * a + g * c_poly
    @constraint(model, b * divergence(fagc) - α * differentiate(b, x) ⋅ fagc in SOSCone())
    @objective(model, Min, c[1])
    optimize!(model)
    if termination_status(model) != MOI.OPTIMAL
        @warn("Termination status $(termination_status(model)): $(raw_status(model))")
    end
    u = value(c_poly) / value(a)
    return MultivariatePolynomials.mapcoefficientsnz(coef -> abs(coef) < 1e-6 ? 0.0 : coef, u)
end

import DifferentialEquations, Plots
function phase_plot(f, quiver_scaling, Δt, X0, solver = DifferentialEquations.Tsit5())
    ∇(vx, vy) = [fi(x[1] => vx, x[2] => vy) for fi in f]
    ∇pt(v, p, t) = ∇(v[1], v[2])
    function traj(v0)
        tspan = (0.0, Δt)
        prob = DifferentialEquations.ODEProblem(∇pt, v0, tspan)
        return DifferentialEquations.solve(prob, solver, reltol=1e-8, abstol=1e-8)
    end
    ticks = -5:0.5:5
    X = repeat(ticks, 1, length(ticks))
    Y = X'
    Plots.quiver(X, Y, quiver = (x, y) -> ∇(x, y) / quiver_scaling, linewidth=0.5)
    for x0 in X0
        Plots.plot!(traj(x0), vars=(1, 2), label = nothing)
    end
    Plots.plot!(xlims = (-5, 5), ylims = (-5, 5))
end

# ## Example 1

g = [0, 1]
f = [x[2] - x[1]^3 + x[1]^2, 0]
b = 3x[1]^2 + 2x[1]*x[2] + 2x[2]^2
u = controller(f, g, b, 4, 0:3)
@test monomials(u) == [x[2]^3, x[1], x[2]] #src

# We find the controller above which gives the following phase plot.

phase_plot(f + g * u, 200, 10.0, [[x1, x2] for x1 in -5:5:5, x2 in -5:5:5 if x1 != 0 || x2 != 0])

# ## Example 2

g = [0, 1]
f = [2x[1]^3 + x[1]^2*x[2] - 6x[1]*x[2]^2 + 5x[2]^3, 0]
b = x[1]^2 + x[2]^2
u = controller(f, g, b, 2.5, 0:3)
@test monomials(u) == [x[1]^3, x[1]^2*x[2], x[1]*x[2]^2, x[2]^3] #src

# We find the controller above which gives the following phase plot.

phase_plot(f + g * u, 2000, 5.0, [[-1.0, -5.0], [1.0, 5.0]])

# ## Example 3

g = [0, x[2]]
f = [-6x[1]*x[2]^2 - x[1]^2*x[2] + 2x[2]^3, 0]
b = x[1]^2 + x[2]^2
u = controller(f, g, b, 3, 0:2)
@test monomials(u) == [x[1]^2, x[2]^2] #src

# We find the controller above which gives the following phase plot.

X0 = [Float64[x1, x2] for x1 in -5:5:5, x2 in -5:5:5 if x2 != 0]
ε = 1e-4 # We separate the starting point slightly from the hyperplane `x2 = 0` which is invariant.
push!(X0, [-4,  ε])
push!(X0, [-3, -ε])
push!(X0, [ 3,  ε])
push!(X0, [ 4, -ε])
phase_plot(f + g * u, 2000, 10.0, X0)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
const tab = " "^4

struct ZSHCompletionCtx end

function codegen(ctx::ZSHCompletionCtx, cmd::EntryCommand)
    name = cmd_name(cmd)
    "#compdef _$name $name \n" * codegen(ctx, "_", cmd.root, true)
end

function codegen(ctx::ZSHCompletionCtx, prefix::String, cmd::NodeCommand, entry::Bool)
    lines = [
        "# These are set by _arguments",
        "local context state state_descr line",
        "typeset -A opt_args",
        "",
    ]

    args = basic_arguments(entry)

    push!(args, "\"1: :$(actions(cmd.subcmds))\"")
    push!(args, "\"*:: :->args\"")

    push!(lines, "_arguments -C \\")
    append!(lines, map(x -> tab * x * " \\", args))

    push!(lines, "")
    push!(lines, raw"case $state in")
    push!(lines, tab * "(args)")
    push!(lines, tab * tab * raw"case ${words[1]} in")

    commands = []
    for each in cmd.subcmds
        name = cmd_name(each)
        push!(commands, name * ")")
        push!(commands, tab * prefix * cmd_name(cmd) * "_" * name)
        push!(commands, ";;")
    end
    append!(lines, map(x -> tab^3 * x, commands))
    push!(lines, tab^2 * "esac")
    push!(lines, "esac")

    body = join(map(x -> tab * x, lines), "\n")

    script = []
    push!(
        script,
        """
function $prefix$(cmd_name(cmd))() {
$body
}
""",
    )

    for each in cmd.subcmds
        push!(script, codegen(ctx, prefix * cmd_name(cmd) * "_", each, false))
    end

    return join(script, "\n\n")
end

function codegen(ctx::ZSHCompletionCtx, prefix::String, cmd::LeafCommand, entry::Bool)
    lines = ["_arguments \\"]
    args = basic_arguments(entry)

    for (i, each) in enumerate(cmd.options)
        name = cmd_name(each)
        doc = cmd_doc(each).first
        token = "--$name"
        if each.short
            token = "{-$(short_name(each)),--$name}"
        end
        push!(args, "$token'[$doc]'")
    end

    append!(lines, map(x -> tab * x * " \\", args))
    body = join(map(x -> tab * x, lines), "\n")

    return """
    function $prefix$(cmd_name(cmd))() {
    $body
    }
    """
end

function basic_arguments(entry)
    args = ["'(- 1 *)'{-h,--help}'[show help information]'"]

    if entry
        push!(args, "'(- 1 *)'{-V,--version}'[show version information]'")
    end
    return args
end

function actions(args)
    hints = map(args) do x
        str = cmd_name(x) * "\\:"
        str *= "'" * cmd_doc(x).first * "'"
    end

    return "((" * join(hints, " ") * "))"
end


text:
# broadcasting behavior
Broadcast.broadcastable(loss::SupervisedLoss) = Ref(loss)

# fallback to unary evaluation
value(loss::DistanceLoss, target::Number, output::Number)  = value(loss, output - target)
deriv(loss::DistanceLoss, target::Number, output::Number)  = deriv(loss, output - target)
deriv2(loss::DistanceLoss, target::Number, output::Number) = deriv2(loss, output - target)

value(loss::MarginLoss, target::Number, output::Number)  = value(loss, target * output)
deriv(loss::MarginLoss, target::Number, output::Number)  = target * deriv(loss, target * output)
deriv2(loss::MarginLoss, target::Number, output::Number) = deriv2(loss, target * output)

# result type when applying the loss to a single pair of objects
result_type(loss::SupervisedLoss, t::Type, o::Type) = typeof(value(loss, zero(t), zero(o)))

# ------------------
# AVAILABLE LOSSES
# ------------------
include("supervised/distance.jl")
include("supervised/margin.jl")
include("supervised/other.jl")

# meta-losses
include("supervised/ordinal.jl")
include("supervised/scaled.jl")
include("supervised/weighted.jl")

# helper macro (for devs)
macro dimcheck(condition)
    :(($(esc(condition))) || throw(DimensionMismatch("Dimensions of the parameters don't match: $($(string(condition)))")))
end

# ------------------------------
# DEFAULT AGGREGATION BEHAVIOR
# ------------------------------
for FUN in (:value, :deriv, :deriv2)
    @eval begin
        # by default compute the element-wise result
        @inline function ($FUN)(
                loss::SupervisedLoss,
                targets::AbstractArray,
                outputs::AbstractArray)
            ($FUN)(loss, targets, outputs, AggMode.None())
        end

        # translate ObsDim.Last to the correct ObsDim.Constant (for code reduction)
        @inline function ($FUN)(
                loss::SupervisedLoss,
                targets::AbstractArray,
                outputs::AbstractArray{T,N},
                agg::AggregateMode,
                ::ObsDim.Last = ObsDim.Last()) where {T,N}
            ($FUN)(loss, targets, outputs, agg, ObsDim.Constant{N}())
        end

        # -------------------
        # AGGREGATION: NONE
        # -------------------
        @generated function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,M},
                output::AbstractArray{T,N},
                ::AggMode.None) where {Q,M,T,N}
            quote
                $(Expr(:meta, :inline))
                ($($FUN)).(loss, target, output)
            end
        end

        # ------------------
        # AGGREGATION: SUM
        # ------------------
        @generated function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,M},
                output::AbstractArray{T,N},
                ::AggMode.Sum) where {Q,M,T,N}
            bigger = M > N ? :target : :output
            S, B = min(M,N), max(M,N)
            quote
                @nexprs $S (n)->@dimcheck(size(target, n) == size(output, n))
                out = zero(result_type(loss, Q, T))
                @inbounds @simd for I in CartesianIndices(size($bigger))
                    @nexprs $B n->(i_n = I[n])
                    out += ($($FUN))(loss, @nref($M,target,i), @nref($N,output,i))
                end
                out
            end
        end

        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                agg::AggMode.Sum,
                obsdim::ObsDim.Constant{O}) where {Q,T,N,O}
            N == 1 && throw(ArgumentError("Sum per observation non sensible for two Vectors. Try omitting the obsdim"))
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            S = result_type(loss, Q, T)
            out = zeros(S, size(output, O))
            @inbounds @simd for I in CartesianIndices(size(output))
                out[I[O]] += ($FUN)(loss, target[I], output[I])
            end
            out
        end

        # -------------------
        # AGGREGATION: MEAN
        # -------------------
        @generated function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,M},
                output::AbstractArray{T,N},
                ::AggMode.Mean) where {Q,M,T,N}
            bigger = M > N ? :target : :output
            S, B = min(M,N), max(M,N)
            quote
                @nexprs $S (n)->@dimcheck(size(target, n) == size(output, n))
                nrm = 1 / length($bigger)
                out = zero(result_type(loss, Q, T)) * nrm
                @inbounds @simd for I in CartesianIndices(size($bigger))
                    @nexprs $B n->(i_n = I[n])
                    out += ($($FUN))(loss, @nref($M,target,i), @nref($N,output,i)) * nrm
                end
                out
            end
        end

        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                agg::AggMode.Mean,
                obsdim::ObsDim.Constant{O}) where {Q,T,N,O}
            N == 1 && throw(ArgumentError("Mean per observation non sensible for two Vectors. Try omitting the obsdim"))
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            S = result_type(loss, Q, T)
            out = zeros(S, size(output, O))
            nrm = 1 / S(prod(size(output,n) for n in 1:N if n != O))
            @inbounds @simd for I in CartesianIndices(size(output))
                out[I[O]] += ($FUN)(loss, target[I], output[I]) * nrm
            end
            out
        end

        # ---------------------------
        # AGGREGATION: WEIGHTED SUM
        # ---------------------------
        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                agg::AggMode.WeightedSum,
                ::ObsDim.Constant{O}) where {Q,T,N,O}
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            @dimcheck size(output, O) == length(agg.weights)
            nrm = agg.normalize ? inv(sum(agg.weights)) : inv(one(sum(agg.weights)))
            out = zero(result_type(loss, Q, T)) * (agg.weights[1] * nrm)
            @inbounds @simd for I in CartesianIndices(size(output))
                out += ($FUN)(loss, target[I], output[I]) * (agg.weights[I[O]] * nrm)
            end
            out
        end

        # ----------------------------
        # AGGREGATION: WEIGHTED MEAN
        # ----------------------------
        function ($FUN)(
                loss::SupervisedLoss,
                target::AbstractArray{Q,N},
                output::AbstractArray{T,N},
                agg::AggMode.WeightedMean,
                ::ObsDim.Constant{O}) where {Q,T,N,O}
            O > N && throw(ArgumentError("The specified obsdim is larger as the available dimensions."))
            @dimcheck size(target) == size(output)
            @dimcheck size(output, O) == length(agg.weights)
            k = prod(n != O ? size(output,n) : 1 for n in 1:N)::Int
            nrm = agg.normalize ? inv(k * sum(agg.weights)) : inv(k * one(sum(agg.weights)))
            out = zero(result_type(loss, Q, T)) * (agg.weights[1] * nrm)
            @inbounds @simd for I in CartesianIndices(size(output))
                out += ($FUN)(loss, target[I], output[I]) * (agg.weights[I[O]] * nrm)
            end
            out
        end
    end
end

# convenient functor interface
if VERSION ≥ v"1.3.0"
    (loss::SupervisedLoss)(target::AbstractArray, output::AbstractArray) = value(loss, target, output)
else
    # add method manually to all subtypes
    for L in Iterators.flatten([subtypes(DistanceLoss), subtypes(MarginLoss)])
        (loss::L)(target::AbstractArray, output::AbstractArray) = value(loss, target, output)
    end
end


text:
"""
    get_SOA_table(id)

Given the id of a `mort.SOA.org` table, grab it and return it as a `MortalityTable`.

!!! Remember that not all tables have been tested to work.
"""
function get_SOA_table(id::Int)
    path = "https://mort.soa.org/Export.aspx?Type=xml&TableIdentity=$id"
    r = HTTP.request("GET", path, require_ssl_verification=false)

    # Why skip the first three bytes of the response?

    # From https://docs.python.org/3/library/codecs.html
    # To increase the reliability with which a UTF-8 encoding can be detected,
    # Microsoft invented a variant of UTF-8 (that Python 2.5 calls "utf-8-sig")
    # for its Notepad program: Before any of the Unicode characters is written
    # to the file, a UTF-8 encoded BOM (which looks like this as a byte sequence:
    # 0xef, 0xbb, 0xbf) is written.
    if r.body[1:3] == [0xef, 0xbb, 0xbf]
        xml = getXML(String(r.body[4:end]))
    else
        xml = getXML(String(r.body))
    end
    return XTbML_Table_To_MortalityTable(parseXTbMLTable(xml, path))

end

"""
    get_SOA_table!(dict,id)

Will lookup the given mortality table and add it to the given dict, with the name of
the table acting as the added key in the dictionary.

This modifies the given `dict` (as is indicated by the conventional `!` at the end
of the function name).

!! Remember that not all tables have been tested to work.
"""
function get_SOA_table!(dict, id::Int)
    tbl = get_SOA_table(id)
    merge!(dict, Dict(tbl.metadata.name => tbl))

end


text:
@testset "Data compression" begin
    N = 100
    Ws = []
    for k in 1:2
        W = randn(N,N)
        for i in 1:N
            W[i,i] = 0
            for j in 1:i
                W[i,j] = - W[j,i]
            end
        end
        push!(Ws, W)
    end
    D = Diagonal([2.0^(-j) for j in 1:N])
    Y(t) = exp(t*Ws[1])*exp(t)*D*exp(t*Ws[2])
    Y0 = Y(0.0)
    Yf = Y(1.0)
    X0 = truncated_svd(Y0, tol=1e-4)
    prob = MatrixDataProblem(Y, X0, (0.0, 1.0))
    data = [Y(t) for t in range(0, 1.0, step = 0.01)]
    discrete_prob = MatrixDataProblem(data, X0)
    solvers = [UnconventionalAlgorithm(),
               ProjectorSplitting(PrimalLieTrotter()),
               ProjectorSplitting(DualLieTrotter()),
               RankAdaptiveUnconventionalAlgorithm(1e-8, rmax = 20),
               GreedyIntegrator()]
    for solver in solvers
        sol = LowRankIntegrators.solve(prob, solver, 1e-2)
        discrete_sol = LowRankIntegrators.solve(discrete_prob, solver)
        @test all(Matrix(sol.Y[end]) .≈ Matrix(discrete_sol.Y[end]))
    end
end

#sol = LowRankIntegrators.solve(prob, solver, 1e-2)


text:
@test Benitor.fibonacciTail(0,1,0) == 0
@test Benitor.fibonacciTail(1,1,0) == 1
@test Benitor.fibonacciTail(2,1,0) == 1
@test Benitor.fibonacciTail(4,1,0) == 3
@test Benitor.fibonacciTail(10,1,0) == 55
@test Benitor.fibonacciTail(50,1,0) == 12586269025
@test Benitor.fibonacciTail(100,1,0) == 3736710778780434371


text:
"""
    dominate(p, q)

Returns `1` if `p` dominates `q`, `-1` if otherwise, and `0` if dominance cannot be determined.
"""
function dominate(p::T, q::T) where {T<:AbstractArray}
    ret = 0
    for (i, j) in zip(p, q)
        if i < j
            ret == -1 && return 0
            ret = 1
        elseif j < i
            ret == 1 && return 0
            ret = -1
        end
    end
    return ret
end
dominate(p::T, q::T) where {T<:Individual} = dominate(objectives(p), objectives(q))

"""
dominations(P::AbstractVector)

Returns a domination matrix of all elements in the input collection `P`.
"""
function dominations(P::AbstractVector{T}) where {T<:AbstractArray}
    l = length(P)
    D = zeros(Int8, l, l)
    for i = 1:l
        for j = (i+1):l
            D[i, j] = dominate(P[i], P[j])
            D[j, i] = -D[i, j]
        end
    end
    D
end
dominations(P::Population) = dominations(objectives(P))

# todo: 把 P 转置一下
"""
    nondominatedsort!(R, F)

Calculate fronts for fitness values `F`, and store ranks of the individuals into `R`.
"""
function nondominatedsort!(R, P)
    n = size(P, 2)
    @assert length(R) == n "Ranks must be defined for the whole population"

    Sₚ = Dict(i => Set() for i = 1:n)
    C = zeros(Int, n)

    # construct first front
    F = [Int[]]
    for i = 1:n
        for j = i+1:n
            r = dominate(view(P, :, i), view(P, :, j)) #M[i,j]
            if r == 1
                push!(Sₚ[i], j)
                C[j] += 1
            elseif r == -1
                push!(Sₚ[j], i)
                C[i] += 1
            end
        end
        if C[i] == 0
            R[i] = 1
            push!(F[1], i)
        end
    end

    # construct rest of the fronts
    while !isempty(last(F))
        Q = Int[]
        for i in last(F)
            for j in Sₚ[i]
                C[j] -= 1
                if C[j] == 0
                    push!(Q, j)
                    R[j] = length(F) + 1
                end
            end
        end
        push!(F, Q)
    end
    isempty(last(F)) && pop!(F)

    F #, R #, Sₚ
end
function nondominatedsort!(P::Population)
    R = fill(0, length(P))
    nondominatedsort!(R, objectives(P)')
    for i = 1:length(P)
        P[i].rank = R[i]
    end
    P
end

"""
    crowding_distance!((D, P)

Calculate crowding distance for individuals and save the results into `D`
given the fitness values `P` and collection of `F`.
"""
function crowding_distance!(D::AbstractVector, I::AbstractMatrix) where {T}
    N, _ = size(I)

    for m = eachcol(I)
        s = sortperm(m)
        D[s[1]] = Inf
        D[s[end]] = Inf

        max, _ = findmax(m)
        min, _ = findmin(m)

        for i = 2:N-1
            D[i] += (m[s[i+1]] - m[s[i-1]]) / (max - min)
        end
    end
end
function crowding_distance!(P::Population)
    maxF, _ = findmax(map(rank, P))
    F = [Vector{Int64}() for _ in 1:maxF]
    for i = 1:length(P)
        push!(F[rank(P[i])], i)
    end
    for f = F
        C = fill(0.0, length(f))
        crowding_distance!(C, objectives(P[f]))
        for i = 1:length(f)
            P[f][i].distance = C[i]
        end
    end
    P
end

"""
    compare(a, b)
    compares whether two vectors are dominated or not.
    Output:
    `1` if argument 1 (a) dominates argument 2 (b).
    `2` if argument 2 (b) dominates argument 1 (a).
    `3` if both arguments 1 (a) and 2 (b) are incomparable.
    `0` if both arguments 1 (a) and 2 (b) are equal.
"""
function compare(a::Vector, b::Vector)
    k = length(a)
    @assert k == length(b)

    i = 1
    while i <= k && a[i] == b[i]
        i += 1;
    end

    if i > k
        return 0 # equals
    end

    if a[i] < b[i]

        for j = i+1:k# (j = i+1; j <= k; ++j)
            if b[j] < a[j]
                return 3 #a and b are incomparable
            end
        end

        return 1; #  a dominates b
    end

    for j = i+1:k #(j = i+1; j < k; ++j)
        if (a[j] < b[j])
            return 3 #; // a and b are incomparable
        end
    end

    return 2 # b dominates a

end

function compare(a::AbstractIndividual, b::AbstractIndividual)
    compare(objectives(a), objectives(b))
end

"""
    get_non_dominated_solutions_perm(population)
Return a vector of integers `v` such that `population[v]` are the non dominated
solutions contained in `population`.
"""
function get_non_dominated_solutions_perm(population)
    ids = Int[1]
    n = length(population)

    for i in 2:n
        j = 1
        while j <= length(ids)
            jj = ids[j]
            relation = compare(population[i], population[jj])

            if relation == 2 # j dominates i
                break
            elseif relation == 1 # i dominates j
                deleteat!(ids, j)
                continue
            end

            j += 1
        end

        if j > length(ids)
            push!(ids, i)
        end

    end

    return ids
end

"""
    get_non_dominated_solutions(population)
Return the non dominated solutions contained in `population`.
"""
function get_non_dominated_solutions(population)
    return population[get_non_dominated_solutions_perm(population)]
end


text:
using Logging
using ArgParse
using Distributed
using Random

@everywhere using CGP

@everywhere CGP.Config.init("cfg/test.yaml")

@everywhere function xorfit(c::Chromosome, nbits::Int64, nsamples::Int64)
    accuracy = 0
    for d in 1:nsamples
        ins = bitrand(nbits)
        outputs = process(c, Float64.(ins))
        if outputs[1] == xor(ins...)
            accuracy += 1
        end
    end
    accuracy /= nsamples
    accuracy
end

function get_args()
    argtable = ArgParseSettings()
    @add_arg_table(
        argtable,
        "--seed", arg_type = Int, default = 0,
        "--log", arg_type = String, default = "xor.log",
        "--nbits", arg_type = Int, default = 2,
        "--nsamples", arg_type = Int, default = 10
    )
    parse_args(argtable)
end

function run_xor(args::Dict)
    Random.seed!(args["seed"])
    global_logger(SimpleLogger(open(args["log"], "a+")))
    nin = args["nbits"]
    nout = 1
    fitness = x->xorfit(x, args["nbits"], args["nsamples"])
    maxfit, best = GA(nin, nout, fitness; id="xor")
end

if ~isinteractive()
    run_xor(get_args())
end


text:
module SplittablesTesting

import SplittablesBase

for n in SplittablesBase.Testing._PUBLIC_API
    @eval const $n = SplittablesBase.Testing.$n
    d = @eval @doc(SplittablesBase.Testing.$n)
    @eval @doc($d, $n)
end

end


text:
using Convex
using Convex: DotMultiplyAtom
using Test
using SDPA_GMP
using MathOptInterface
using Random
const MOI = MathOptInterface
using GenericLinearAlgebra
using LinearAlgebra
import Random.shuffle
import Statistics.mean

TOL = 1e-3
eye(n) = Matrix(1.0I, n, n)

# Seed random number stream to improve test reliability
Random.seed!(2)

solvers = Any[]

push!(solvers, SDPA_GMP.Optimizer{BigFloat}())

# If Gurobi is installed, uncomment to test with it:
#using Gurobi
#push!(solvers, GurobiSolver(OutputFlag=0))

# If Mosek is installed, uncomment to test with it:
#using Mosek
#push!(solvers, MosekSolver(LOG=0))

@testset "Convex" begin
    # include("test_utilities.jl")
    # include(joinpath("test","test_const.jl"))
    include(joinpath("test_affine.jl"))
    include(joinpath("test_lp.jl"))
    # include("test_socp.jl")
    # include(joinpath("test_sdp.jl"))
    # include("test_exp.jl")
    # include("test_sdp_and_exp.jl")
    # include("test_mip.jl")
end

max_stars_repo_path: test/runtests (2).jl

max_issues_repo_path: test/runtests (2).jl

max_forks_repo_path: test/runtests (2).jl


text:
# https://adventofcode.com/2020/day/10
using DataStructures
using BenchmarkTools

input = read("2020/data/day_10.txt", String)

test_input = """
28
33
18
42
31
14
46
20
48
47
24
23
49
45
19
38
39
11
1
32
25
35
8
17
7
9
4
2
34
10
3
"""

function parse_input(input)
    numbers = [0]
    inp = sort(parse.(Int, split(strip(input), "\n")))
    append!(numbers, inp)
    append!(numbers, numbers[end] + 3)
    numbers
end

function part_1(input)
    numbers = parse_input(input)
    counts = counter(numbers[2:end] - numbers[1:end - 1])
    counts[3] * counts[1]
end
# @info part_1(test_input)
@info part_1(input)

function part_2(input)
    numbers = parse_input(input)
    routes = zeros(Int, length(numbers))
    for (idx, n) in enumerate(numbers)
        if idx == 1
            routes[idx] = 1
        else
            is_connector = numbers[idx] - 3 .<= numbers .< numbers[idx]
            routes[idx] = sum(routes[is_connector])
        end
    end
    routes[end]
end
# @info part_2(test_input)
@info part_2(input)
# @btime part_2(input)

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
include(joinpath(pwd(), "src/models/acrobot.jl"))

# Horizon
T = 101

# Time step
tf = 5.0
h = tf / (T - 1)

# Bounds

# ul <= u <= uu
ul, uu = control_bounds(model, T, -10.0, 10.0)

# Initial and final states
x1 = [0.0; 0.0; 0.0; 0.0]
xT = [π; 0.0; 0.0; 0.0]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(ones(model.n)) : Diagonal(ones(model.n)) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T], [zeros(model.m) for t = 1:T])

# Problem
prob = trajectory_optimization_problem(model,
           obj,
           T,
           h = h,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu,
           )

# Trajectory initialization
X0 = linear_interp(x1, xT, T) # linear interpolation on state
U0 = random_controls(model, T, 0.001) # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

# Solve
@time Z̄ = solve(prob, copy(Z0))

# Visualize
using Plots
X̄, Ū = unpack(Z̄, prob)
plot(hcat(X̄...)', width = 2.0)
plot(hcat(Ū...)', width = 2.0, linetype = :steppost)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
import QuantumOpticsBase:
    tracenorm, tracenorm_h, tracenorm_nh, tracedistance, tracedistance_h, tracedistance_nh
import LinearAlgebra: ishermitian

"""
    tracenorm(rho)
Trace norm of `rho`.
It is defined as
```math
T(ρ) = \\operatorname{tr}\\{\\sqrt{ρ^† ρ}\\}.
```
Depending if `rho` is hermitian either [`tracenorm_h`](@ref) or
[`tracenorm_nh`](@ref) is called.
"""
function tracenorm(rho::DenseSuperOpType)
    ishermitian(rho) ? tracenorm_h(rho) : tracenorm_nh(rho)
end

#=
function tracenorm!(rho::DenseSuperOpType)
    if ishermitian(rho)
        rho.data ./= tracenorm_h(rho)
    else
        rho.data ./= tracenorm_nh(rho)
    end
end
=#

"""
    tracenorm_h(rho)
Trace norm of `rho`.
It uses the identity
```math
T(ρ) = \\operatorname{tr}\\{\\sqrt{ρ^† ρ}\\} = \\sum_i |λ_i|
```
where ``λ_i`` are the eigenvalues of `rho`.
"""
function tracenorm_h(rho::DenseSuperOpType{B,B,T}) where {B<:Tuple{Basis,Basis},T}
    s = eigvals(Hermitian(rho.data))
    sum(abs.(s))
end

#=
function tracenorm_h(rho::T) where {T<:AbstractOperator}
    throw(ArgumentError("tracenorm_h not implemented for $(typeof(rho)). Use dense operators instead."))
end
=#

"""
    tracenorm_nh(rho)
Trace norm of `rho`.
Note that in this case `rho` doesn't have to be represented by a square
matrix (i.e. it can have different left-hand and right-hand bases).
It uses the identity
```math
    T(ρ) = \\operatorname{tr}\\{\\sqrt{ρ^† ρ}\\} = \\sum_i σ_i
```
where ``σ_i`` are the singular values of `rho`.
"""
tracenorm_nh(rho::SuperOperator{B,B,T}) where {B<:Tuple{Basis,Basis},T} =
    sum(svdvals(rho.data))

"""
    tracedistance(rho, sigma)
Trace distance between `rho` and `sigma`.
It is defined as
```math
T(ρ,σ) = \\frac{1}{2} \\operatorname{tr}\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\}.
```
It calls [`tracenorm`](@ref) which in turn either uses [`tracenorm_h`](@ref)
or [`tracenorm_nh`](@ref) depending if ``ρ-σ`` is hermitian or not.
"""
tracedistance(
    rho::SuperOperator{B,B,T},
    sigma::SuperOperator{B,B,T},
) where {B<:Tuple{Basis,Basis},T} = 0.5 * tracenorm(rho - sigma)

"""
    tracedistance_h(rho, sigma)
Trace distance between `rho` and `sigma`.
It uses the identity
```math
T(ρ,σ) = \\frac{1}{2} \\operatorname{tr}\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\} = \\frac{1}{2} \\sum_i |λ_i|
```
where ``λ_i`` are the eigenvalues of `rho` - `sigma`.
"""
tracedistance_h(
    rho::SuperOperator{B,B,T},
    sigma::SuperOperator{B,B,T},
) where {B<:Tuple{Basis,Basis},T} = 0.5 * tracenorm_h(rho - sigma)

"""
    tracedistance_nh(rho, sigma)
Trace distance between `rho` and `sigma`.
Note that in this case `rho` and `sigma` don't have to be represented by square
matrices (i.e. they can have different left-hand and right-hand bases).
It uses the identity
```math
    T(ρ,σ) = \\frac{1}{2} \\operatorname{tr}\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\}
         = \\frac{1}{2} \\sum_i σ_i
```
where ``σ_i`` are the singular values of `rho` - `sigma`.
"""
tracedistance_nh(
    rho::SuperOperator{B1,B2,T},
    sigma::SuperOperator{B1,B2,T},
) where {B1<:Tuple{Basis,Basis},B2<:Tuple{Basis,Basis},T} = 0.5 * tracenorm_nh(rho - sigma)

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
# Local compute version

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.precompile()

using ArgParse
using Gadfly
using DelimitedFiles
using Dates, Statistics
using YAML
using JLD2
using CoordinateTransformations, Rotations, StaticArrays
using ImageCore
using Images, ImageDraw
# using ImageView
using AprilTags
using DataInterpolations

# using CuArrays
using Flux

using DistributedFactorGraphs
using IncrementalInference
using RoME
using RoMEPlotting
using Caesar

## Load all required packages
using Distributed

if nprocs() == 1
  addprocs(7) # make sure there are 4 processes waiting before loading packages
  # machines = [(ENV["JLCLST02"],3);]
end

@everywhere begin
  using Pkg
  Pkg.activate(@__DIR__)
  Pkg.instantiate()
end

using Flux
using RoME
using RoMEPlotting

WP = nprocs == 1 ? WorkerPool([1]) : WorkerPool(2:nprocs() |> collect )

# using DistributedFactorGraphs
for i in procs()
  fetch( Distributed.@spawnat i @eval using DistributedFactorGraphs )
  fetch( Distributed.@spawnat i @eval using IncrementalInference )
  fetch( Distributed.@spawnat i @eval using Flux )
  fetch( Distributed.@spawnat i @eval using RoME )
  fetch( Distributed.@spawnat i @eval using Caesar )
  # fetch( Distributed.@spawnat i @eval using CuArrays )
end

for i in procs()
  fetch( Distributed.@spawnat i @eval using RoMEPlotting )
end

@everywhere using JLD2
@everywhere using DistributedFactorGraphs
@everywhere using IncrementalInference
@everywhere using RoME
@everywhere using Caesar
@everywhere using RoMEPlotting
@everywhere using Flux

# @show ARGS
include("parsecommands.jl")

@everywhere begin
# using Fontconfig
# using Compose
using Gadfly
using Cairo
# using FileIO
# using GeometryTypes # using MeshCat
end

# setup configuration
include(joinpath(@__DIR__,"configParameters.jl") )

# Figure export folder
global currdirtime = now()
if parsed_args["previous"] != ""
  # currdirtime = "2018-11-07T01:36:52.274"
  currdirtime = parsed_args["previous"]
end

resultsparentdir = joinpath(datadir, "results")
resultsdir = joinpath(resultsparentdir, "$(currdirtime)")

if parsed_args["previous"] == ""
  # When running fresh from new data
  include(joinpath(@__DIR__,"createResultsDirs.jl"))
end

# Utils required for this processing script
include(joinpath(@__DIR__,"racecarUtils.jl") )
include(joinpath(@__DIR__,"cameraUtils.jl") )
# include(joinpath(dirname(@__FILE__),"visualizationUtils.jl") )

@everywhere function plotRacecarInterm(fgl::G, resultsdirl, psyml::Symbol)::Nothing where G <: AbstractDFG
  @show ls(fgl)
  pl = drawPosesLandms(fgl, spscale=0.1, drawhist=false, meanmax=:mean) #,xmin=-3,xmax=6,ymin=-5,ymax=2);
  Gadfly.draw(PNG(joinpath(resultsdirl, "images", "$(psyml).png"),15cm, 10cm),pl)
  pl = drawPosesLandms(fgl, spscale=0.1, meanmax=:mean) # ,xmin=-3,xmax=3,ymin=-2,ymax=2);
  Gadfly.draw(PNG(joinpath(resultsdirl, "images", "hist_$(psyml).png"),15cm, 10cm),pl)
  # pl = plotPose2Vels(fgl, Symbol("$(psyml)"), coord=Coord.Cartesian(xmin=-1.0, xmax=1.0))
  # Gadfly.draw(PNG(joinpath(resultsdirl, "images", "vels_$(psyml).png"),15cm, 10cm),pl)
  # save combined image with tags
  nothing
end

global tag_bag = Dict()
if parsed_args["previous"] == ""
  tag_bag = detectTagsInImgs(datafolder, imgfolder, resultsdir, camidxs, iterposes=parsed_args["iterposes"])
  # save the tag bag file for future use
  @save resultsdir*"/tag_det_per_pose.jld2" tag_bag
else
  @load resultsdir*"/tag_det_per_pose.jld2" tag_bag
end

## Load joystick time series data

runnr = parse(Int, parsed_args["folder_name"][end])
joyTimeseries = readdlm(joinpath(datadir,parsed_args["folder_name"],"labRun$(runnr)_joy.csv"), ',')
joyTimeseries = joyTimeseries[2:end,[1,6,8]]
joyTimeseries[:,1] .= joyTimeseries[:,1]*1e-9 # .|> unix2datetime;

# load the detections file to get timestamps
detcData = readdlm(joinpath(datadir,parsed_args["folder_name"],"labRun$(runnr)Detections.csv"), ',')
detcData = detcData[2:end,:]
detcPoseTs = detcData[:,4]*1e-9 .+ detcData[:,3] .|> unix2datetime

## find timeseries segments that go with each pose
joyTsDict = Dict{Symbol, Array{Float64,2}}()
poseTimes = Dict{Symbol, DateTime}(:x0 => detcPoseTs[1])
mask = joyTimeseries[:,1] .< datetime2unix(detcPoseTs[1])
joyTsDict[:x0] = joyTimeseries[mask,:]

for i in 2:length(detcPoseTs)
  mask = datetime2unix(detcPoseTs[i-1]) .<= joyTimeseries[:,1] .< datetime2unix(detcPoseTs[i])
  joyTsDict[Symbol("x$(i-1)")] = joyTimeseries[mask,:]
  poseTimes[Symbol("x$(i-1)")] = detcPoseTs[i-1] # later helper
end

## interpolate up to PyQuest values...

intJoyDict = Dict{Symbol,Vector{Vector{Float64}}}()
# for sym, lclJD = :x1, joyTsDict[:x1]
for (sym, lclJD) in joyTsDict
  intJoyDict[sym] = interpTo25x4(lclJD)
  # if 1 < size(lclJD,1)
  #   tsLcl = range(lclJD[1,1],lclJD[end,1],length=25)
  #   intrTrTemp = DataInterpolations.LinearInterpolation(lclJD[:,2],lclJD[:,1])
  #   intrStTemp = DataInterpolations.LinearInterpolation(lclJD[:,3],lclJD[:,1])
  #   newVec = Vector{Vector{Float64}}()
  #   for tsL in tsLcl
  #     newVal = zeros(4)
  #     newVal[1] = intrTrTemp(tsL)
  #     newVal[2] = intrStTemp(tsL)
  #     push!(newVec, newVal)
  #   end
  #   # currently have no velocity values
  #   intJoyDict[sym] = newVec
  # else
  #   intJoyDict[sym] = [zeros(4) for i in 1:25]
  # end
end

## More code required

inlcude(joinpath(dirname(pathof(Caesar)), "examples", "learning", "hybrid", "FluxModelsPose2Pose2", "FluxModelsPose2Pose2.jl"))

include(joinpath(@__DIR__, "LoadPyNNTxt.jl"))

## load the required models into common predictor

allModels = []
for i in 0:99
# /home/dehann/data/racecar/results/conductor/models/retrained_network_weights0
  push!(allModels, loadPose2OdoNNModelIntoFlux(ENV["HOME"]*"/data/racecar/results/conductor/models/retrained_network_weights$i") ) # loadTfModelIntoFlux
end

## run the solution

fg = main(WP, resultsdir, camidxs, tag_bag, jldfile=parsed_args["jldfile"], failsafe=parsed_args["failsafe"], show=parsed_args["show"], odopredfnc=(x)->JlOdoPredictorPoint2(x, allModels), joyvel=intJoyDict, poseTimes=poseTimes, multiproc=true  )

## development

0

# fg = initfg()
#
# addVariable!(fg, :x0, Pose2, timestamp=poseTimes[:x0])
# addVariable!(fg, :x1, Pose2, timestamp=poseTimes[:x1])
#
# pp2 = PriorPose2(MvNormal(zeros(3),LinearAlgebra.diagm([0.01;0.01;0.005].^2)))
# addFactor!(fg, [:x0], pp2)
#
# # the neural factor
# DXmvn = MvNormal(zeros(3),LinearAlgebra.diagm([0.01;0.01;0.005].^2))
# odopredfnc=(x)->JlOdoPredictorPoint2(x, allModels)
# joyvel=intJoyDict
#
# joyVals = zeros(25,4)
# for i in 1:25
#   joyVals[i,:] .= joyvel[:x0][i]
# end
#
# pp = FluxModelsPose2Pose2(odopredfnc,joyVals, DXmvn,0.4)  # joyvel[:x0]
# addFactor!(fg, [:x0;:x1], pp)
#
# pts = approxConv(fg, :x0x1f1, :x1)

# fails
# key 1 not found
# julia101 -p 4 apriltag_and_zed_slam.jl --previous "2018-11-09T01:42:33.279" --jldfile "racecar_fg_x299.jld2" --folder_name "labrun7" --failsafe

#


text:
#---------------------------------------------------------------------#
#This function plots the solution
#Written by F.X. Giraldo on April 19 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function plot_solution(q0,qb,coord,space_method,time,case)

    #Exact Solution
    Npoin=size(q0,1)
    (qe,qb,gravity) = exact_solution(coord,Npoin,time,case,DFloat)

    #Plot Solution
    p1 = plot(coord,qe[:,1],linecolor = :blue,label = "Analytic Solution",legend = :bottomright,legendfont=font(5))
    plot!(p1,coord,q0[:,1],linestyle = :dashdot,linecolor = :red,label = "$space_method Solution")
    title!(p1,"Height (hs)")
    xlabel!(p1,"x")
    ylabel!(p1,"hs(x,t)")
    plot!(p1,xlims = (first(coord),last(coord)+.1),ylims = (-1,0.5))

    p2 = plot(coord,qe[:,2],linecolor = :blue,label = "Analytic Solution",legend = :bottomright,legendfont=font(5))
    plot!(p2,coord,q0[:,2],linestyle = :dashdot,linecolor = :red,label = "$space_method Solution")
    title!(p2,"Momentum (U)")
    xlabel!(p2,"x")
    ylabel!(p2,"U(x,t)")
    plot!(p2,xlims = (first(coord),last(coord)+.1),ylims = (-1,0.5))

    figure_1 = plot(p1,p2,layout = (1,2))
    display(figure_1)
end